// =====================================================================
// CC patch evaluation. JS port of VtMeshDecoder::EvalPatchQuadOrTri
// (VtMeshDecoder.cpp:3150). Given a (patchIter, uv) for a quad or
// triangle source patch, returns the (P, N) on the Catmull-Clark
// limit surface at that parametric location.
//
// Inputs come from Phase 1 (CC topology build), Phase 2 (limit
// positions + normals), and Phase 3 (irregular wedge subdivision):
//
//   meshData.m_subdivData -- m_polyInfo, m_regionData,
//                            m_irregularWedgeInfo, m_weightInflData,
//                            m_vertexInfo, m_numRegularPoints,
//                            m_numIrregularVertexWedges,
//                            m_numIrregularTriangles
//   limitPosData / limitNrmData      -- from Phase 2's
//                                       CalcLimitPointsAndNormals.
//   cpRedirect                       -- weightInflData.m_regularAndIrregularPoints
//   reusableCpData                   -- one Vec3 per runtime CP,
//                                       built by ccBeginFrameCompute
//                                       ReusableCpData below.
//   fullWedgeDstPoints / Num / Start -- Phase 3's per-frame wedge
//                                       subdivision output.
//   irregularLevels                  -- per-wedge requested subdiv
//                                       level (== irregularProcLevel).
//
// Two paths inside the helper:
//   - Regular sub-region: 4x4 bicubic B-spline eval with
//     sharpness-aware basis on whichever border edge is sharp.
//   - Irregular sub-region: routed through the per-wedge fan of
//     subdivision slices (Phase 3 output).
'use strict';

// ---- Vec3 helpers (flat-xyz Float32Array, indexed by Vec3 slot) ----
function cpe_v3GetFromFlat(flat, i) {
  return [flat[i*3+0], flat[i*3+1], flat[i*3+2]];
}
function cpe_v3GetFromArray(arr, i) {
  // arr is amolite::vector<Vec3> JS equivalent -- an Array of [x,y,z].
  return arr[i];
}
function cpe_v3Sub(a, b)  { return [a[0]-b[0], a[1]-b[1], a[2]-b[2]]; }
function cpe_v3Scale(a, s){ return [a[0]*s, a[1]*s, a[2]*s]; }
function cpe_v3Add(a, b)  { return [a[0]+b[0], a[1]+b[1], a[2]+b[2]]; }
function cpe_v3Lerp(a, b, t) {
  return [a[0]+(b[0]-a[0])*t, a[1]+(b[1]-a[1])*t, a[2]+(b[2]-a[2])*t];
}
function cpe_v3Cross(a, b) {
  return [a[1]*b[2] - a[2]*b[1], a[2]*b[0] - a[0]*b[2], a[0]*b[1] - a[1]*b[0]];
}
function cpe_v3Normalize(a) {
  const m2 = a[0]*a[0] + a[1]*a[1] + a[2]*a[2];
  if (m2 === 0) return [0, 0, 0];
  const inv = 1.0 / Math.sqrt(m2);
  return [a[0]*inv, a[1]*inv, a[2]*inv];
}

// ---- pull in the BSpline pieces from cc_top_peology.js ----
// Imported via require() in Node; in the browser, the same names land
// on globalThis so we can use them directly. Each module uses its
// own variable name to avoid collisions when multiple <script> tags
// share global scope.
let cc_top_pe;
if (typeof require !== 'undefined') {
  cc_top_pe = require('./cc_topology.js');
} else {
  // Browser: globalThis carries function declarations but NOT the
  // top-level `const`s (kSubdivPolyType_*, kRegionScale_*, etc.) --
  // those bindings are script-scope only. Build a facade that
  // overlays the consts via their lexical-scope names.
  cc_top_pe = Object.assign({}, globalThis, {
    kSubdivPolyType_QuadRegular, kSubdivPolyType_QuadIrregular1,
    kSubdivPolyType_QuadIrregular2, kSubdivPolyType_Triangle,
    kRegionScale_00, kRegionScale_01, kRegionScale_10, kRegionScale_11,
    kRegionScale_Full, kRegionScale_Num,
    SUBDIV_MAX_SHARPNESS, IRREGULAR_GROUP_SIZE,
    kCpType_Invalid, kCpType_Direct, kCpType_Border,
    kCpType_ConcaveCorner, kCpType_ConvexCorner, kCpType_Num,
  });
}

const SUBDIV_MAX_SHARPNESS_PE = 10;

// ---- ccBeginFrameComputeReusableCpData -----------------------------
// Port of the m_reusableCpData population loop in
// VtMeshDecoder::BeginFrame (VtMeshDecoder.cpp:2312-2318). Builds one
// Vec3 per runtime CP from the source positions, via
// CalcPointFromControlPointsAndWeights (DcciSubdiv.cpp:4806).
//
// posData is the source-vertex positions as a flat Float32Array (3*N).
function ccBeginFrameComputeReusableCpData(weightInflData, posDataFlat) {
  const numTotalCps = weightInflData.m_numBaseCps + weightInflData.m_cpItemStart.length;
  const out = new Array(numTotalCps);

  for (let i = 0; i < numTotalCps; i++) {
    let P;
    if (i < weightInflData.m_numBaseCps) {
      // CalcPointFromControlPointsAndWeights: base CPs map directly
      // to posData[ptIndex].
      P = [posDataFlat[i*3+0], posDataFlat[i*3+1], posDataFlat[i*3+2]];
    } else {
      const cpIndex = i - weightInflData.m_numBaseCps;
      const start   = weightInflData.m_cpItemStart[cpIndex];
      const count   = weightInflData.m_cpItemCount[cpIndex];
      let x = 0, y = 0, z = 0;
      for (let k = 0; k < count; k++) {
        const srcIndex = weightInflData.m_inflData[start + k];
        const w        = weightInflData.m_weightData[start + k];
        x += posDataFlat[srcIndex*3+0] * w;
        y += posDataFlat[srcIndex*3+1] * w;
        z += posDataFlat[srcIndex*3+2] * w;
      }
      P = [x, y, z];
    }
    out[i] = P;
  }
  return out;
}

// ---- SubdivUtil::SubQuadBaryFromTriBary (DcciSubdiv.cpp:6972) ----
function subQuadBaryFromTriBary(triBary) {
  const a = -1.0;
  const b = (3 - 2 * triBary[0] + 2 * triBary[1]);
  const c = -6 * triBary[1];

  const disc     = b * b - 4 * a * c;
  const quadratic= Math.sqrt(disc > 0 ? disc : 0);
  const y0 = (-b + quadratic) / (2.0 * a);
  // amolite keeps y0; y1 is computed but unused. We mirror.
  const denom = 0.5 - (1.0 / 6.0) * y0;
  const x0 = triBary[0] * (denom !== 0 ? (1.0 / denom) : 0);
  return [x0, y0];
}

// ---- RotateTriBaryForCorner (DcciSubdiv.cpp:6992) ----
function rotateTriBaryForCorner(baseCorner, triBary) {
  if (baseCorner === 0) return [triBary[0], triBary[1]];
  if (baseCorner === 1) return [1.0 - (triBary[0] + triBary[1]), triBary[0]];
  if (baseCorner === 2) return [triBary[1], 1.0 - (triBary[0] + triBary[1])];
  throw new Error('rotateTriBaryForCorner: baseCorner ' + baseCorner);
}

// ---- SubdivUtil::SubQuadAndUvFromTriBaryLevel1 (DcciSubdiv.cpp:7019) ----
function subQuadAndUvFromTriBaryLevel1(triBary) {
  // segment0: y = 1 - 2x  (under: y < 1 - 2x)
  // segment1: y = 0.5 - 0.5x
  // segment2: y = x
  const underS0 = triBary[1] < 1.0 - 2.0 * triBary[0];
  const underS1 = triBary[1] < 0.5 - 0.5 * triBary[0];
  const underS2 = triBary[1] < triBary[0];

  let baseCorner = 0;
  if (!underS2 && !underS1) baseCorner = 1;
  if (underS2 && !underS0)  baseCorner = 2;

  const rotateBary = rotateTriBaryForCorner(baseCorner, triBary);
  const subQuadBary = subQuadBaryFromTriBary(rotateBary);
  return { foundCorner: baseCorner, subQuadBary };
}

// ---- RegionInfo::GetSubRegionRange (DcciSubdiv.cpp:3644) ----
// regionScale is the eRegionScale enum: 0=00 1=01 2=10 3=11 4=Full.
function regionInfoGetSubRegionRange(regionScale) {
  if (regionScale < cc_top_pe.kRegionScale_Full) {
    const x = regionScale % 2;
    const y = (regionScale / 2) | 0;
    return {
      uvOffset: [x / 2.0, y / 2.0],
      uvScale:  [0.5, 0.5],
    };
  }
  if (regionScale === cc_top_pe.kRegionScale_Full) {
    return { uvOffset: [0, 0], uvScale: [1, 1] };
  }
  throw new Error('regionInfoGetSubRegionRange: bad regionScale ' + regionScale);
}

// ---- BSplineUtil::EvalBicubicBSplineSharpness (Vec3 variant; DcciSubdiv.cpp:8998) ----
// Returns { foundP, foundN } where foundN = normalize(cross(d/du, d/dv)).
//
// 2-axis sharpness extension: when sharpnessY > 0, the Y axis ALSO
// runs through the sharpness-aware basis (reverseY controls which
// end of Y has the kink). Used by the Fix-3 multi-sharp case for
// 2-adjacent (perpendicular) creases. sharpnessY == 0 falls back to
// the standard smooth-Y B-spline basis -- byte-equal to the original
// 1-axis path so existing single-side sharp regions are unchanged.
function evalBicubicBSplineSharpness(ptOffsets, ptData, sharpness, tx, ty,
                                     sharpnessY, reverseY, reverseX) {
  // reverseX defaults to true to match the original 1-axis path
  // (sharp at low X = side 0 = L after canonical rotation).
  const rx = (reverseX == null) ? true : !!reverseX;
  const sx = cc_top_pe.bsplineUtilCalcBasisAndDerivForSharpnessRef(tx, sharpness, rx);
  const basisX  = sx.basis;
  const basisDx = sx.basisDeriv;

  let basisY, basisDy;
  if (sharpnessY != null && sharpnessY > 0) {
    const sy = cc_top_pe.bsplineUtilCalcBasisAndDerivForSharpnessRef(ty, sharpnessY, !!reverseY);
    basisY  = sy.basis;
    basisDy = sy.basisDeriv;
  } else {
    basisY  = cc_top_pe.bsplineUtilEvalBSplineBasis(ty);
    basisDy = cc_top_pe.bsplineUtilEvalBSplineBasisDeriv(ty);
  }

  const dataX  = [basisX.x,  basisX.y,  basisX.z,  basisX.w];
  const dataY  = [basisY.x,  basisY.y,  basisY.z,  basisY.w];
  const dataDx = [basisDx.x, basisDx.y, basisDx.z, basisDx.w];
  const dataDy = [basisDy.x, basisDy.y, basisDy.z, basisDy.w];

  let sumP  = [0, 0, 0];
  let sumDx = [0, 0, 0];
  let sumDy = [0, 0, 0];

  for (let r = 0; r < 4; r++) {
    for (let c = 0; c < 4; c++) {
      const w_x  = dataX[c];
      const w_dx = dataDx[c];
      const w_y  = dataY[r];
      const w_dy = dataDy[r];

      const offsetVtx = ptOffsets[r * 4 + c];
      const P = ptData[offsetVtx];

      const wxy   = w_x  * w_y;
      const wdxy  = w_dx * w_y;
      const wxdy  = w_x  * w_dy;
      sumP [0] += P[0] * wxy;   sumP [1] += P[1] * wxy;   sumP [2] += P[2] * wxy;
      sumDx[0] += P[0] * wdxy;  sumDx[1] += P[1] * wdxy;  sumDx[2] += P[2] * wdxy;
      sumDy[0] += P[0] * wxdy;  sumDy[1] += P[1] * wxdy;  sumDy[2] += P[2] * wxdy;
    }
  }
  return {
    foundP: sumP,
    foundN: cpe_v3Normalize(cpe_v3Cross(sumDx, sumDy)),
  };
}

// ---- GetPointsOffsetsForIrregularRegion (VtMeshDecoder.cpp:2431) ----
// Computes the 4x4 grid of point offsets into the wedge slice data
// for an irregular sub-region. Region 0 is unreachable (asserts).
// Regions 1/2/3 each have their own layout (region 1 is rotated 90deg
// so the sharp edge is along tx=0).
function getPointsOffsetsForIrregularRegion(targetRegion, currCorner, valence, isBorder,
                                            sliceStartOffset, irregularSharpness, sharpnessStart) {
  const numThreadsForWedge = isBorder ? (valence + 1) : valence;
  const indexPtsInSlice = new Array(16).fill(-1);

  const prevCorner     = (currCorner + (numThreadsForWedge - 1)) % numThreadsForWedge;
  const nextCorner     = (currCorner + 1) % numThreadsForWedge;
  const nextNextCorner = (currCorner + 2) % numThreadsForWedge;

  let foundSharpness = 0.0;
  const isPrevBorder = (isBorder && currCorner === 0);
  const isNextBorder = (isBorder && currCorner === valence - 1);

  if (targetRegion === 0) {
    throw new Error('getPointsOffsetsForIrregularRegion: region 0 unreachable');
  } else if (targetRegion === 3) {
    foundSharpness = 0;
    indexPtsInSlice[0]  = sliceStartOffset + (0);
    indexPtsInSlice[1]  = sliceStartOffset + (isNextBorder ? (1 + nextCorner * 12 + 0) : (1 + nextCorner * 12 + 0));
    indexPtsInSlice[2]  = sliceStartOffset + (isNextBorder ? (1 + nextCorner * 12 + 1) : (1 + nextCorner * 12 + 4));
    indexPtsInSlice[3]  = sliceStartOffset + (isNextBorder ? (1 + nextCorner * 12 + 2) : (1 + nextCorner * 12 + 8));

    indexPtsInSlice[4]  = sliceStartOffset + (1 + currCorner * 12 + 0);
    indexPtsInSlice[5]  = sliceStartOffset + (1 + currCorner * 12 + 1);
    indexPtsInSlice[6]  = sliceStartOffset + (1 + currCorner * 12 + 2);
    indexPtsInSlice[7]  = sliceStartOffset + (1 + currCorner * 12 + 3);

    indexPtsInSlice[8]  = sliceStartOffset + (1 + currCorner * 12 + 4);
    indexPtsInSlice[9]  = sliceStartOffset + (1 + currCorner * 12 + 5);
    indexPtsInSlice[10] = sliceStartOffset + (1 + currCorner * 12 + 6);
    indexPtsInSlice[11] = sliceStartOffset + (1 + currCorner * 12 + 7);

    indexPtsInSlice[12] = sliceStartOffset + (1 + currCorner * 12 + 8);
    indexPtsInSlice[13] = sliceStartOffset + (1 + currCorner * 12 + 9);
    indexPtsInSlice[14] = sliceStartOffset + (1 + currCorner * 12 + 10);
    indexPtsInSlice[15] = sliceStartOffset + (1 + currCorner * 12 + 11);
  } else if (targetRegion === 2) {
    foundSharpness = (sharpnessStart >= 0) ? irregularSharpness[sharpnessStart + currCorner] : 0.0;
    if (isPrevBorder) foundSharpness = SUBDIV_MAX_SHARPNESS_PE;

    indexPtsInSlice[0]  = sliceStartOffset + (isPrevBorder ? (0) : (1 + prevCorner * 12 + 0));
    indexPtsInSlice[1]  = sliceStartOffset + (0);
    indexPtsInSlice[2]  = sliceStartOffset + (isNextBorder ? (1 + nextCorner * 12 + 0) : (1 + nextCorner * 12 + 0));
    indexPtsInSlice[3]  = sliceStartOffset + (isNextBorder ? (1 + nextCorner * 12 + 1) : (1 + nextCorner * 12 + 4));

    indexPtsInSlice[4]  = sliceStartOffset + (isPrevBorder ? (0) : (1 + prevCorner * 12 + 1));
    indexPtsInSlice[5]  = sliceStartOffset + (1 + currCorner * 12 + 0);
    indexPtsInSlice[6]  = sliceStartOffset + (1 + currCorner * 12 + 1);
    indexPtsInSlice[7]  = sliceStartOffset + (1 + currCorner * 12 + 2);

    indexPtsInSlice[8]  = sliceStartOffset + (isPrevBorder ? (0) : (1 + prevCorner * 12 + 2));
    indexPtsInSlice[9]  = sliceStartOffset + (1 + currCorner * 12 + 4);
    indexPtsInSlice[10] = sliceStartOffset + (1 + currCorner * 12 + 5);
    indexPtsInSlice[11] = sliceStartOffset + (1 + currCorner * 12 + 6);

    indexPtsInSlice[12] = sliceStartOffset + (isPrevBorder ? (0) : (1 + prevCorner * 12 + 3));
    indexPtsInSlice[13] = sliceStartOffset + (1 + currCorner * 12 + 8);
    indexPtsInSlice[14] = sliceStartOffset + (1 + currCorner * 12 + 9);
    indexPtsInSlice[15] = sliceStartOffset + (1 + currCorner * 12 + 10);
  } else if (targetRegion === 1) {
    foundSharpness = (sharpnessStart >= 0) ? irregularSharpness[sharpnessStart + nextCorner] : 0.0;
    if (isNextBorder) foundSharpness = SUBDIV_MAX_SHARPNESS_PE;

    // 90 degree rotation so the sharp edge sits along tx=0.
    indexPtsInSlice[0]  = sliceStartOffset + (isNextBorder ? (0) : (1 + nextCorner * 12 + 9));
    indexPtsInSlice[1]  = sliceStartOffset + (isNextBorder ? (1 + nextCorner * 12 + 2) : (1 + nextCorner * 12 + 8));
    indexPtsInSlice[2]  = sliceStartOffset + (1 + currCorner * 12 + 3);
    indexPtsInSlice[3]  = sliceStartOffset + (1 + currCorner * 12 + 7);

    indexPtsInSlice[4]  = sliceStartOffset + (isNextBorder ? (0) : (1 + nextCorner * 12 + 5));
    indexPtsInSlice[5]  = sliceStartOffset + (isNextBorder ? (1 + nextCorner * 12 + 1) : (1 + nextCorner * 12 + 4));
    indexPtsInSlice[6]  = sliceStartOffset + (1 + currCorner * 12 + 2);
    indexPtsInSlice[7]  = sliceStartOffset + (1 + currCorner * 12 + 6);

    indexPtsInSlice[8]  = sliceStartOffset + (isNextBorder ? (0) : (1 + nextCorner * 12 + 1));
    indexPtsInSlice[9]  = sliceStartOffset + (isNextBorder ? (1 + nextCorner * 12 + 0) : (1 + nextCorner * 12 + 0));
    indexPtsInSlice[10] = sliceStartOffset + (1 + currCorner * 12 + 1);
    indexPtsInSlice[11] = sliceStartOffset + (1 + currCorner * 12 + 5);

    indexPtsInSlice[12] = sliceStartOffset + (isNextBorder ? (0) : (1 + nextNextCorner * 12 + 0));
    indexPtsInSlice[13] = sliceStartOffset + (0);
    indexPtsInSlice[14] = sliceStartOffset + (1 + currCorner * 12 + 0);
    indexPtsInSlice[15] = sliceStartOffset + (1 + currCorner * 12 + 4);
  } else {
    throw new Error('getPointsOffsetsForIrregularRegion: bad region ' + targetRegion);
  }

  return { indexPtsInSlice, foundSharpness };
}

// ---- EvaluateIrregularPointFull (VtMeshDecoder.cpp:2574) ----
function evaluateIrregularPointFull(limitP, limitN, allPts, irregularPatchStart,
                                    irregularSharpness, sharpnessStart,
                                    valence, isBorder, patchSliceLevel,
                                    currCorner, uv) {
  if (!(currCorner < valence)) {
    throw new Error('evaluateIrregularPointFull: currCorner ' + currCorner + ' >= valence ' + valence);
  }
  // Fractional-level threshold: the gap boundary uvT slides continuously
  // across the LAST computed slice as patchSliceLevel sweeps from
  // (numSlices-1) up to numSlices. At an integer level N, numSlices=N and
  // minT lands at the inner corner of slice N-1 (uvT=2^-N). As the user
  // raises the level to N+f (0<f<1), numSlices becomes N+1, a new outer-
  // most-meaningful slice gets computed, and minT slides linearly across
  // it from 2^-N (outer corner of slice N) down to 2^-(N+1) (inner corner).
  // Old code used minT = 2^-L which positioned the threshold OUTSIDE the
  // last slice's parametric domain at fractional L, clamping regionUv and
  // making the transition snap to the slice's outer corner instead of
  // sliding through it.
  const numSlices = Math.max(1, Math.ceil(patchSliceLevel));
  const innerFrac = Math.max(0.0, Math.min(1.0, numSlices - patchSliceLevel));
  const minT = (1.0 + innerFrac) * Math.pow(2.0, -numSlices);

  const isUvXLessThanY = uv[0] < uv[1];
  const uvT = Math.max(uv[0], uv[1]);
  const evalT = Math.max(uvT, minT);

  const ratioT  = uvT * (evalT !== 0 ? (1.0 / evalT) : 0);
  const uvScale = ratioT !== 0 ? (1.0 / ratioT) : 0;

  const uvEval = [uvScale * uv[0], uvScale * uv[1]];
  const levelBaseT = Math.max(uvEval[0], uvEval[1]);
  const levelT = -Math.log2(Math.max(1e-5, levelBaseT));

  const sliceLevelForGridCalculation = Math.min(numSlices - 1, Math.floor(levelT));
  const halfGridSize = 1.0 / (1 << (sliceLevelForGridCalculation + 1));

  //  x --- x --- gridSize
  //  |  2  |  3  |
  //  x --- x --- x
  //  |  0  |  1  |
  let regionIndex;
  if (isUvXLessThanY) {
    regionIndex = (uvEval[0] < halfGridSize) ? 2 : 3;
  } else {
    regionIndex = (uvEval[1] < halfGridSize) ? 1 : 3;
  }

  const regionStartUvX = regionIndex % 2;
  const regionStartUvY = (regionIndex / 2) | 0;

  const regionUvX = (uvEval[0] - regionStartUvX * halfGridSize) / halfGridSize;
  const regionUvY = (uvEval[1] - regionStartUvY * halfGridSize) / halfGridSize;

  // Region 1 has its 90-degree-rotated point layout, so the uv flips too.
  let regionUv;
  if (regionIndex === 1) {
    regionUv = [cc_top_pe.saturateF(regionUvY), cc_top_pe.saturateF(1.0 - regionUvX)];
  } else {
    regionUv = [cc_top_pe.saturateF(regionUvX), cc_top_pe.saturateF(regionUvY)];
  }

  const numPtsPerSlice = cc_top_pe.calcNumDstPointsPerSliceForWedge(valence, isBorder);
  const sliceStartOffset = irregularPatchStart + sliceLevelForGridCalculation * numPtsPerSlice;

  const { indexPtsInSlice, foundSharpness: rawSharpness } =
    getPointsOffsetsForIrregularRegion(regionIndex, currCorner, valence, isBorder,
                                       sliceStartOffset, irregularSharpness, sharpnessStart);

  let foundSharpness = rawSharpness;
  if (foundSharpness < SUBDIV_MAX_SHARPNESS_PE) {
    foundSharpness = Math.max(foundSharpness - sliceLevelForGridCalculation, 0.0);
  }
  foundSharpness = Math.max(0.0, foundSharpness - 1.0);

  const { foundP: evalP, foundN: evalN } = evalBicubicBSplineSharpness(
    indexPtsInSlice, allPts, foundSharpness, regionUv[0], regionUv[1]);

  const foundP = cpe_v3Lerp(limitP, evalP, ratioT);
  const foundN = cpe_v3Normalize(cpe_v3Lerp(limitN, evalN, ratioT));
  return { foundP, foundN };
}

// ---- EvalSurfaceUvPointHelper (VtMeshDecoder.cpp:2687) ----
function evalSurfaceUvPointHelper(meshData, regionGrid, regularSharpness, subdivLevel,
                                  limitPosData, limitNrmData, cpRedirect,
                                  reusableCpData, fullWedgeDstPoints,
                                  wedgeNumDstPoints, wedgeStartPoints,
                                  irregularLevels, patchUv) {
  const regionDim = (1 << subdivLevel);

  const baseUvX = patchUv[0] * regionDim;
  const baseUvY = patchUv[1] * regionDim;
  function clampInt(lo, hi, v) { return v < lo ? lo : (v > hi ? hi : v); }
  const regionC = clampInt(0, regionDim - 1, Math.floor(baseUvX));
  const regionR = clampInt(0, regionDim - 1, Math.floor(baseUvY));
  const regionUv = [baseUvX - regionC, baseUvY - regionR];

  const regionBaseIndex = regionR * regionDim + regionC;
  const redirectIndex = cc_top_pe.regionGridGetRegionOffset(regionGrid, regionBaseIndex);

  if (redirectIndex < 0) {
    // ---- regular sub-region: bicubic B-spline over 4x4 patch ----
    const sharpnessSides = [0, 0, 0, 0];
    const isRegionL = (regionC === 0);
    const isRegionR = (regionC === regionDim - 1);
    const isRegionD = (regionR === 0);
    const isRegionU = (regionR === regionDim - 1);

    let rotation = 0;
    let numAdded = 0;
    if (isRegionL && regularSharpness[0] > 0.0) { rotation = 0; numAdded++; sharpnessSides[0] = regularSharpness[0]; }
    if (isRegionU && regularSharpness[1] > 0.0) { rotation = 1; numAdded++; sharpnessSides[1] = regularSharpness[1]; }
    if (isRegionR && regularSharpness[2] > 0.0) { rotation = 2; numAdded++; sharpnessSides[2] = regularSharpness[2]; }
    if (isRegionD && regularSharpness[3] > 0.0) { rotation = 3; numAdded++; sharpnessSides[3] = regularSharpness[3]; }

    let ptStart, ptDirX, ptDirY, localUv;
    if (rotation === 0)      { ptStart = [0, 0]; ptDirX = [ 1,  0]; ptDirY = [ 0,  1]; localUv = [regionUv[0],            regionUv[1]]; }
    else if (rotation === 1) { ptStart = [0, 3]; ptDirX = [ 0, -1]; ptDirY = [ 1,  0]; localUv = [1.0 - regionUv[1],      regionUv[0]]; }
    else if (rotation === 2) { ptStart = [3, 3]; ptDirX = [-1,  0]; ptDirY = [ 0, -1]; localUv = [1.0 - regionUv[0], 1.0 - regionUv[1]]; }
    else                     { ptStart = [3, 0]; ptDirX = [ 0,  1]; ptDirY = [-1,  0]; localUv = [regionUv[1],       1.0 - regionUv[0]]; }

    // AMOLITE_DEFERRED_FIXES Fix 3 -- multi-sharp regions. amolite's
    // single-axis rotation-cascade basis can't handle two-sharp
    // regions directly. We dispatch by configuration:
    //
    //   numAdded == 1: original 1-axis sharpness path (rotation
    //     cascade already done above).
    //   numAdded == 2 adjacent (L+U, U+R, R+D, D+L): apply sharpness
    //     on BOTH X and Y axes of the canonical bicubic; no window
    //     rotation. The cross-product (NDF) of the two sharpnesses
    //     models the perpendicular crease intersection.
    //   numAdded == 2 opposite (L+R or U+D): both ends of the same
    //     axis are sharp. Eval twice (one sharpness at each end)
    //     and average; the averaged basis stays C0-continuous with
    //     adjacent patches and converges to the correct limit on
    //     flat-input meshes (grid3x3_parallel_creases).
    //   numAdded >= 3: smooth fallback (sharpness = 0 on canonical
    //     window). Same as before -- conservative.
    let sharpness;
    let multiSharpMode = 0;   // 0=normal, 1=adjacent, 2=opposite, 3=fallback
    let multiSharpY = 0, multiSharpReverseY = false;
    let multiSharpReverseX = true;
    let oppositeSharpA = 0, oppositeSharpB = 0;   // for the average path
    let oppositeAxis = 0;     // 0=X (L+R), 1=Y (U+D)
    const sL = sharpnessSides[0], sU = sharpnessSides[1];
    const sR = sharpnessSides[2], sD = sharpnessSides[3];
    // Gate: multi-sharp basis math is disabled by default while we
    // investigate why fixtures with multi-sharp cases (single_crease,
    // suzanne, bigguy, monsterfrog) look broken at tess=0 even
    // though CC_PerVertex and CC_Gpu_Patch are bit-exact.
    const useMultiSharp = (typeof window !== 'undefined' &&
                            window.STATE && window.STATE.useMultiSharpBasis);
    if (numAdded === 1) {
      sharpness = sharpnessSides[rotation];
    } else if (numAdded === 2 && useMultiSharp) {
      // Always use canonical orientation for multi-sharp.
      rotation = 0;
      ptStart = [0, 0]; ptDirX = [1, 0]; ptDirY = [0, 1];
      localUv = [regionUv[0], regionUv[1]];
      if (sL > 0 && sR > 0) {
        multiSharpMode = 2; oppositeAxis = 0;
        oppositeSharpA = sL; oppositeSharpB = sR;
        sharpness = 0;
      } else if (sU > 0 && sD > 0) {
        multiSharpMode = 2; oppositeAxis = 1;
        oppositeSharpA = sU; oppositeSharpB = sD;
        sharpness = 0;
      } else {
        // adjacent: pick X / Y axis sharpness + reverse flags.
        multiSharpMode = 1;
        const sx = (sL > 0) ? sL : sR;
        multiSharpReverseX = (sL > 0);   // reverse=true => sharp at low X
        const sy = (sU > 0) ? sU : sD;
        multiSharpReverseY = (sD > 0);   // reverse=true => sharp at low Y
        sharpness = sx;
        multiSharpY = sy;
      }
    } else {
      // 3+ sharp: smooth fallback (pre-existing behavior).
      multiSharpMode = 3;
      rotation = 0;
      ptStart = [0, 0]; ptDirX = [1, 0]; ptDirY = [0, 1];
      localUv = [regionUv[0], regionUv[1]];
      sharpness = 0;
    }

    const start = regionGrid.m_pointsStart;
    const fullCpOffsets = new Array(16);
    for (let r = 0; r < 4; r++) {
      for (let c = 0; c < 4; c++) {
        const offX = ptStart[0] + c * ptDirX[0] + r * ptDirY[0];
        const offY = ptStart[1] + c * ptDirX[1] + r * ptDirY[1];
        const srcC = regionC + offX;
        const srcR = regionR + offY;
        const localOffset = srcR * (regionDim + 3) + srcC;
        const redirect = cpRedirect[start + localOffset];
        fullCpOffsets[r * 4 + c] = redirect;
      }
    }

    let foundP, foundN;
    if (multiSharpMode === 2) {
      // 2-opposite: eval twice with sharpness at each end, average.
      // Axis 0 (L+R): sharp on X at both ends.
      // Axis 1 (U+D): sharp on Y at both ends.
      let p1, p2, n1, n2;
      if (oppositeAxis === 0) {
        const r1 = evalBicubicBSplineSharpness(
          fullCpOffsets, reusableCpData, oppositeSharpA, localUv[0], localUv[1],
          0, false, true);   // sharpX=A, reverseX=true (sharp at low X = L)
        const r2 = evalBicubicBSplineSharpness(
          fullCpOffsets, reusableCpData, oppositeSharpB, localUv[0], localUv[1],
          0, false, false);  // sharpX=B, reverseX=false (sharp at high X = R)
        p1 = r1.foundP; p2 = r2.foundP; n1 = r1.foundN; n2 = r2.foundN;
      } else {
        // Y axis. Pass sharpX=0 (smooth) so X axis is plain B-spline.
        // For the Y sharpness we use the 2-axis path's sharpnessY/reverseY.
        const r1 = evalBicubicBSplineSharpness(
          fullCpOffsets, reusableCpData, 0, localUv[0], localUv[1],
          oppositeSharpA, false, true);  // sharpY=A, reverseY=false (sharp at high Y = U)
        const r2 = evalBicubicBSplineSharpness(
          fullCpOffsets, reusableCpData, 0, localUv[0], localUv[1],
          oppositeSharpB, true, true);   // sharpY=B, reverseY=true (sharp at low Y = D)
        p1 = r1.foundP; p2 = r2.foundP; n1 = r1.foundN; n2 = r2.foundN;
      }
      foundP = [(p1[0]+p2[0])*0.5, (p1[1]+p2[1])*0.5, (p1[2]+p2[2])*0.5];
      // The second eval uses reverseX=false (or reverseY=true on the
      // Y axis), which flips the corresponding derivative sign vs the
      // first eval -- the cross-product normal then points the
      // opposite way. Averaging would cancel. Use the first eval's N
      // (canonically oriented to match the surrounding 1-side-sharp
      // verts) and trust the avg-position to be the right limit point.
      foundN = n1;
    } else if (multiSharpMode === 1) {
      // 2-adjacent: sharpness on both axes simultaneously.
      const r = evalBicubicBSplineSharpness(
        fullCpOffsets, reusableCpData, sharpness, localUv[0], localUv[1],
        multiSharpY, multiSharpReverseY, multiSharpReverseX);
      foundP = r.foundP;
      // When (reverseX XOR reverseY) is true, the cross product of the
      // two-axis basis derivatives is flipped once -> N inverted vs
      // the canonical outward direction. Compensate.
      if (multiSharpReverseX !== multiSharpReverseY) {
        foundN = [-r.foundN[0], -r.foundN[1], -r.foundN[2]];
      } else {
        foundN = r.foundN;
      }
    } else {
      // numAdded == 0 or 1 or 3+ : single-axis sharpness path.
      const r = evalBicubicBSplineSharpness(
        fullCpOffsets, reusableCpData, sharpness, localUv[0], localUv[1]);
      foundP = r.foundP; foundN = r.foundN;
    }
    return { calcP: foundP, calcN: foundN };
  }

  // ---- irregular sub-region: route through Phase 3 wedge data ----
  let halfGridShift;
  if      (subdivLevel === 1) halfGridShift = 0;
  else if (subdivLevel === 2) halfGridShift = 1;
  else throw new Error('evalSurfaceUvPointHelper: irregular path subdivLevel ' + subdivLevel);

  const halfGridIndex = (regionR >> halfGridShift) * 2 + (regionC >> halfGridShift);
  let rotatedUv;
  if      (halfGridIndex === 0) rotatedUv = [regionUv[0],            regionUv[1]];
  else if (halfGridIndex === 1) rotatedUv = [regionUv[1],            1.0 - regionUv[0]];
  else if (halfGridIndex === 2) rotatedUv = [1.0 - regionUv[1],      regionUv[0]];
  else                          rotatedUv = [1.0 - regionUv[0], 1.0 - regionUv[1]];

  const actualRegion = meshData.m_subdivData.m_regionData[regionGrid.m_regionDataStart + redirectIndex];
  const irregularIndex = actualRegion.m_irregularIndex;
  if (!(irregularIndex >= 0)) throw new Error('evalSurfaceUvPointHelper: irregularIndex < 0');

  const irregularInfo = meshData.m_subdivData.m_irregularWedgeInfo[irregularIndex];
  const runtimeInfo   = meshData.m_subdivData.m_weightInflData.m_runtimeIrregularInfo[irregularIndex];

  // pick the limit pos/nrm source: triangle-center wedges have their
  // own slot at numRegular + numIrregularVertexWedges + triCornerLimit
  // Offset; vertex wedges point at the source vertex.
  let srcLimitIndex = -1;
  if (irregularInfo.m_centerTrianglePoly >= 0) {
    const baseOffset = meshData.m_subdivData.m_numRegularPoints +
                       meshData.m_subdivData.m_numIrregularVertexWedges;
    srcLimitIndex = baseOffset + irregularInfo.m_centerTriangleLimitOffset;
  } else {
    srcLimitIndex = irregularInfo.m_centerVertex;
  }

  const limitN = limitNrmData[srcLimitIndex];
  const limitP = limitPosData[srcLimitIndex];

  if (!(actualRegion.m_irregularPolyIndex < runtimeInfo.m_valence)) {
    throw new Error('evalSurfaceUvPointHelper: irregularPolyIndex >= valence');
  }
  const regionCorner = actualRegion.m_irregularPolyIndex;
  const valence  = runtimeInfo.m_valence;
  const isBorder = !!runtimeInfo.m_isBorder;

  // Sub-region uv offset/scale (the per-region partition of the
  // wedge's parametric domain).
  const { uvOffset, uvScale } = regionInfoGetSubRegionRange(actualRegion.m_irregularRegionScale);

  const wedgeStart = wedgeStartPoints[irregularIndex];
  const wedgeNum   = wedgeNumDstPoints[irregularIndex];

  // amolite copies the wedge's portion of fullWedgeDstPoints into a
  // local buffer + uses irregularPatchStart=0. Match that here so
  // the offset indices stay small.
  const wedgeDstPoints = new Array(wedgeNum);
  for (let i = 0; i < wedgeNum; i++) {
    wedgeDstPoints[i] = fullWedgeDstPoints[wedgeStart + i];
  }

  const baseUv = [rotatedUv[0], rotatedUv[1]];
  const uv = [uvOffset[0] + uvScale[0] * baseUv[0], uvOffset[1] + uvScale[1] * baseUv[1]];

  const currIrregularLevel = irregularLevels[irregularIndex];

  const { foundP, foundN } = evaluateIrregularPointFull(
    limitP, limitN,
    wedgeDstPoints,
    /*irregularPatchStart=*/ 0,
    meshData.m_subdivData.m_weightInflData.m_irregularEdgeSharpness,
    runtimeInfo.m_edgeSharpnessStart,
    valence, isBorder, currIrregularLevel,
    regionCorner, uv);

  return { calcP: foundP, calcN: foundN };
}

// ---- EvalPatchTri (VtMeshDecoder.cpp:3066) ----
function evalPatchTri(meshData, patchIter, baryUv,
                      limitPosData, limitNrmData, cpRedirect,
                      reusableCpData, fullWedgeDstPoints,
                      wedgeNumDstPoints, wedgeStartPoints, irregularLevels) {
  const polyInfo = meshData.m_subdivData.m_polyInfo[patchIter];
  if (polyInfo.m_type !== cc_top_pe.kSubdivPolyType_Triangle) {
    throw new Error('evalPatchTri: poly is not a triangle');
  }

  const { foundCorner, subQuadBary: patchUv } = subQuadAndUvFromTriBaryLevel1(baryUv);
  const cornerGrid = polyInfo.m_triRegionGrids[foundCorner];
  if (polyInfo.m_subdivLevel !== 2) throw new Error('evalPatchTri: subdivLevel != 2');
  if (cornerGrid.m_pointsNum !== 25) throw new Error('evalPatchTri: tri corner pointsNum != 25');

  const edgeSharpness = [
    polyInfo.m_edgeSharpness[(foundCorner + 0) % 3],
    0,
    0,
    polyInfo.m_edgeSharpness[(foundCorner + 2) % 3],
  ];

  return evalSurfaceUvPointHelper(
    meshData, cornerGrid, edgeSharpness, 1,
    limitPosData, limitNrmData, cpRedirect,
    reusableCpData, fullWedgeDstPoints,
    wedgeNumDstPoints, wedgeStartPoints, irregularLevels, patchUv);
}

// ---- EvalPatchQuad (VtMeshDecoder.cpp:3121) ----
function evalPatchQuad(meshData, patchIter, patchUv,
                       limitPosData, limitNrmData, cpRedirect,
                       reusableCpData, fullWedgeDstPoints,
                       wedgeNumDstPoints, wedgeStartPoints, irregularLevels) {
  const polyInfo = meshData.m_subdivData.m_polyInfo[patchIter];
  const t = polyInfo.m_type;
  if (t !== cc_top_pe.kSubdivPolyType_QuadRegular &&
      t !== cc_top_pe.kSubdivPolyType_QuadIrregular1 &&
      t !== cc_top_pe.kSubdivPolyType_QuadIrregular2) {
    throw new Error('evalPatchQuad: bad poly type ' + t);
  }
  return evalSurfaceUvPointHelper(
    meshData, polyInfo.m_quadRegionGrid, polyInfo.m_edgeSharpness, polyInfo.m_subdivLevel,
    limitPosData, limitNrmData, cpRedirect,
    reusableCpData, fullWedgeDstPoints,
    wedgeNumDstPoints, wedgeStartPoints, irregularLevels, patchUv);
}

// ---- EvalPatchQuadOrTri (VtMeshDecoder.cpp:3150) ----
function evalPatchQuadOrTri(meshData, patchIter, baryUv,
                            limitPosData, limitNrmData, cpRedirect,
                            reusableCpData, fullWedgeDstPoints,
                            wedgeNumDstPoints, wedgeStartPoints, irregularLevels) {
  const polyInfo = meshData.m_subdivData.m_polyInfo[patchIter];
  if (polyInfo.m_type === cc_top_pe.kSubdivPolyType_Triangle) {
    return evalPatchTri(meshData, patchIter, baryUv,
                        limitPosData, limitNrmData, cpRedirect,
                        reusableCpData, fullWedgeDstPoints,
                        wedgeNumDstPoints, wedgeStartPoints, irregularLevels);
  }
  return evalPatchQuad(meshData, patchIter, baryUv,
                       limitPosData, limitNrmData, cpRedirect,
                       reusableCpData, fullWedgeDstPoints,
                       wedgeNumDstPoints, wedgeStartPoints, irregularLevels);
}

if (typeof module !== 'undefined' && module.exports) {
  module.exports = {
    ccBeginFrameComputeReusableCpData,
    subQuadBaryFromTriBary, rotateTriBaryForCorner, subQuadAndUvFromTriBaryLevel1,
    regionInfoGetSubRegionRange,
    evalBicubicBSplineSharpness,
    getPointsOffsetsForIrregularRegion,
    evaluateIrregularPointFull,
    evalSurfaceUvPointHelper,
    evalPatchTri, evalPatchQuad, evalPatchQuadOrTri,
  };
}
