path: root/tools/mq_editor/mq_extract.js

// MQ Extraction Algorithm
// McAulay-Quatieri sinusoidal analysis

// Extract partials from audio buffer
function extractPartials(params, stftCache) {
  const {fftSize, threshold, sampleRate, freqWeight, prominence} = params;
  const numFrames = stftCache.getNumFrames();

  const frames = [];
  for (let i = 0; i < numFrames; ++i) {
    const cachedFrame = stftCache.getFrameAtIndex(i);
    const squaredAmp = cachedFrame.squaredAmplitude;
    const phase = cachedFrame.phase;
    const peaks = detectPeaks(squaredAmp, phase, fftSize, sampleRate, threshold, freqWeight, prominence);
    frames.push({time: cachedFrame.time, peaks});
  }

  const partials = trackPartials(frames, params);

  // Second pass: extend partials leftward to recover onset frames
  expandPartialsLeft(partials, frames);

  for (const partial of partials) {
    partial.freqCurve = fitBezier(partial.times, partial.freqs);
    const ac = fitBezier(partial.times, partial.amps);
    partial.freqCurve.a0 = ac.v0; partial.freqCurve.a1 = ac.v1;
    partial.freqCurve.a2 = ac.v2; partial.freqCurve.a3 = ac.v3;
  }

  return {partials, frames};
}

// Helper to interpolate phase via quadratic formula on unwrapped neighbors.
// This provides a more accurate phase estimate at the sub-bin peak location.
function phaseInterp(p_minus, p_center, p_plus, p_frac) {
  // unwrap neighbors relative to center
  let dp_minus = p_minus - p_center;
  while (dp_minus > Math.PI) dp_minus -= 2 * Math.PI;
  while (dp_minus < -Math.PI) dp_minus += 2 * Math.PI;

  let dp_plus = p_plus - p_center;
  while (dp_plus > Math.PI) dp_plus -= 2 * Math.PI;
  while (dp_plus < -Math.PI) dp_plus += 2 * Math.PI;
  
  const p_interp = p_center + (dp_plus - dp_minus) * p_frac * 0.5 + (dp_plus + dp_minus) * 0.5 * p_frac * p_frac;
  return p_interp;
}

// Detect spectral peaks via local maxima + parabolic interpolation
// squaredAmp: pre-computed re*re+im*im per bin
// phase: pre-computed atan2(im,re) per bin
// freqWeight: if true, weight by f before peak detection (f * Power(f))
function detectPeaks(squaredAmp, phase, fftSize, sampleRate, thresholdDB, freqWeight, prominenceDB = 0) {
  const mag = new Float32Array(fftSize / 2);
  const binHz = sampleRate / fftSize;
  for (let i = 0; i < fftSize / 2; ++i) {
    const w = freqWeight ? (i * binHz) : 1.0;
    mag[i] = 10 * Math.log10(Math.max(squaredAmp[i] * w, 1e-20));
  }

  const peaks = [];
  for (let i = 2; i < mag.length - 2; ++i) {
    if (mag[i] > thresholdDB &&
        mag[i] > mag[i-1] && mag[i] > mag[i-2] &&
        mag[i] > mag[i+1] && mag[i] > mag[i+2]) {

      // Check prominence if requested
      if (prominenceDB > 0) {
        let minLeft = mag[i];
        for (let k = i - 1; k >= 0; --k) {
          if (mag[k] > mag[i]) break; // Found higher peak
          if (mag[k] < minLeft) minLeft = mag[k];
        }

        let minRight = mag[i];
        for (let k = i + 1; k < mag.length; ++k) {
          if (mag[k] > mag[i]) break; // Found higher peak
          if (mag[k] < minRight) minRight = mag[k];
        }

        const valley = Math.max(minLeft, minRight);
        if (mag[i] - valley < prominenceDB) continue;
      }

      // Parabolic interpolation for sub-bin accuracy on frequency, amplitude, and phase
      const alpha = mag[i-1];
      const beta  = mag[i];
      const gamma = mag[i+1];
      const p = 0.5 * (alpha - gamma) / (alpha - 2*beta + gamma);
      
      const p_phase = phaseInterp(phase[i-1], phase[i], phase[i+1], p);

      const freq  = (i + p) * sampleRate / fftSize;
      const ampDB = beta - 0.25 * (alpha - gamma) * p;
      peaks.push({freq, amp: Math.pow(10, ampDB / 20), phase: p_phase});
    }
  }

  return peaks;
}

// Helper to compute shortest angle difference (e.g., between -pi and pi)
function normalizeAngle(angle) {
  return angle - 2 * Math.PI * Math.floor((angle + Math.PI) / (2 * Math.PI));
}

// Find best matching peak for a predicted freq/phase. Returns {bestIdx, bestCost}.
function findBestPeak(peaks, matched, predictedFreq, predictedPhase, tol, phaseErrorWeight) {
  let bestIdx = -1, bestCost = Infinity;
  for (let i = 0; i < peaks.length; ++i) {
    if (matched.has(i)) continue;
    const pk = peaks[i];
    const freqError = Math.abs(pk.freq - predictedFreq);
    if (freqError > tol) continue;
    const phaseError = Math.abs(normalizeAngle(pk.phase - predictedPhase));
    const cost = freqError + phaseErrorWeight * phaseError * predictedFreq;
    if (cost < bestCost) { bestCost = cost; bestIdx = i; }
  }
  return { bestIdx, bestCost };
}

// Track partials across frames using phase coherence for robust matching.
function trackPartials(frames, params) {
  const {
    sampleRate, hopSize,
    birthPersistence = 3,
    deathAge = 5,
    minLength = 10,
    phaseErrorWeight = 2.0
  } = params;
  const partials        = [];
  const activePartials  = [];
  const candidates      = []; // pre-birth

  const trackingRatio   = 0.05; // 5% frequency tolerance
  const minTrackingHz   = 20;

  for (const frame of frames) {
    const matched = new Set();

    // --- Continue active partials ---
    for (const partial of activePartials) {
      const lastFreq = partial.freqs[partial.freqs.length - 1];
      const lastPhase = partial.phases[partial.phases.length - 1];
      const velocity = partial.velocity || 0;
      const predictedFreq = lastFreq + velocity;
      
      // Predict phase for the current frame based on the last frame's frequency.
      // Multiply by (age+1) to account for frames missed during a gap.
      const phaseAdvance = 2 * Math.PI * lastFreq * (partial.age + 1) * hopSize / sampleRate;
      const predictedPhase = lastPhase + phaseAdvance;
      
      const tol = Math.max(predictedFreq * trackingRatio, minTrackingHz);
      // Find the peak in the new frame with the lowest cost (freq + phase error).
      const { bestIdx } = findBestPeak(frame.peaks, matched, predictedFreq, predictedPhase, tol, phaseErrorWeight);

      if (bestIdx >= 0) {
        const pk = frame.peaks[bestIdx];
        partial.times.push(frame.time);
        partial.freqs.push(pk.freq);
        partial.amps.push(pk.amp);
        partial.phases.push(pk.phase);
        partial.age = 0;
        partial.velocity = pk.freq - lastFreq;
        matched.add(bestIdx);
      } else {
        partial.age++;
      }
    }

    // --- Advance candidates ---
    for (let i = candidates.length - 1; i >= 0; --i) {
      const cand = candidates[i];
      const lastFreq = cand.freqs[cand.freqs.length - 1];
      const lastPhase = cand.phases[cand.phases.length - 1];
      const velocity = cand.velocity || 0;
      const predictedFreq = lastFreq + velocity;

      // Candidates die on first miss so age is always 0 here, but kept consistent.
      const phaseAdvance = 2 * Math.PI * lastFreq * hopSize / sampleRate;
      const predictedPhase = lastPhase + phaseAdvance;

      const tol = Math.max(predictedFreq * trackingRatio, minTrackingHz);
      const { bestIdx } = findBestPeak(frame.peaks, matched, predictedFreq, predictedPhase, tol, phaseErrorWeight);

      if (bestIdx >= 0) {
        const pk = frame.peaks[bestIdx];
        cand.times.push(frame.time);
        cand.freqs.push(pk.freq);
        cand.amps.push(pk.amp);
        cand.phases.push(pk.phase);
        cand.velocity = pk.freq - lastFreq;
        matched.add(bestIdx);
        // "graduate" a candidate to a full partial
        if (cand.times.length >= birthPersistence) {
          activePartials.push(cand);
          candidates.splice(i, 1);
        }
      } else {
        candidates.splice(i, 1); // kill candidate
      }
    }

    // --- Spawn new candidates from unmatched peaks ---
    for (let i = 0; i < frame.peaks.length; ++i) {
      if (matched.has(i)) continue;
      const pk = frame.peaks[i];
      candidates.push({
        times: [frame.time], 
        freqs: [pk.freq], 
        amps: [pk.amp],
        phases: [pk.phase],
        age: 0,
        velocity: 0
      });
    }

    // --- Kill aged-out partials ---
    for (let i = activePartials.length - 1; i >= 0; --i) {
      if (activePartials[i].age > deathAge) {
        if (activePartials[i].times.length >= minLength) partials.push(activePartials[i]);
        activePartials.splice(i, 1);
      }
    }
  }

  // --- Collect remaining active partials ---
  for (const partial of activePartials) {
    if (partial.times.length >= minLength) partials.push(partial);
  }

  return partials;
}

// Second pass: extend each partial leftward to recover onset frames missed
// by the birthPersistence requirement in the forward pass.
function expandPartialsLeft(partials, frames) {
  const trackingRatio = 0.05;
  const minTrackingHz = 20;

  // Build time → frame index map
  const timeToIdx = new Map();
  for (let i = 0; i < frames.length; ++i) timeToIdx.set(frames[i].time, i);

  for (const partial of partials) {
    if (!partial.phases) partial.phases = []; // Ensure old partials have phase array

    let startIdx = timeToIdx.get(partial.times[0]);
    if (startIdx == null || startIdx === 0) continue;

    for (let i = startIdx - 1; i >= 0; --i) {
      const frame = frames[i];
      const refFreq = partial.freqs[0];
      const tol = Math.max(refFreq * trackingRatio, minTrackingHz);

      let bestIdx = -1, bestDist = Infinity;
      for (let j = 0; j < frame.peaks.length; ++j) {
        const dist = Math.abs(frame.peaks[j].freq - refFreq);
        if (dist < tol && dist < bestDist) { bestDist = dist; bestIdx = j; }
      }

      if (bestIdx < 0) break;

      const pk = frame.peaks[bestIdx];
      partial.times.unshift(frame.time);
      partial.freqs.unshift(pk.freq);
      partial.amps.unshift(pk.amp);
      partial.phases.unshift(pk.phase);
    }
  }
}

// Autodetect spread_above / spread_below from the spectrogram.
// For each (subsampled) STFT frame within the partial, measures the
// half-power (-3dB) width of the spectral peak above and below the center.
// spread = half_bandwidth / f0  (fractional).
function autodetectSpread(partial, stftCache, fftSize, sampleRate) {
  const curve = partial.freqCurve;
  if (!curve || !stftCache) return {spread_above: 0.02, spread_below: 0.02};

  const numFrames = stftCache.getNumFrames();
  const binHz    = sampleRate / fftSize;
  const halfBins = fftSize / 2;

  let sumAbove = 0, sumBelow = 0, count = 0;

  const STEP = 4;
  for (let fi = 0; fi < numFrames; fi += STEP) {
    const frame = stftCache.getFrameAtIndex(fi);
    if (!frame) continue;
    const t = frame.time;
    if (t < curve.t0 || t > curve.t3) continue;

    const f0 = evalBezier(curve, t);
    if (f0 <= 0) continue;

    const sq = frame.squaredAmplitude;
    if (!sq) continue;

    // Find peak bin in ±10% window
    const binCenter = f0 / binHz;
    const searchBins = Math.max(3, Math.round(f0 * 0.10 / binHz));
    const binLo = Math.max(1, Math.floor(binCenter - searchBins));
    const binHi = Math.min(halfBins - 2, Math.ceil(binCenter + searchBins));

    let peakBin = binLo, peakVal = sq[binLo];
    for (let b = binLo + 1; b <= binHi; ++b) {
      if (sq[b] > peakVal) { peakVal = sq[b]; peakBin = b; }
    }

    const halfPower = peakVal * 0.5;  // -3dB in power

    // Walk above peak until half-power, interpolate crossing
    let aboveBin = peakBin;
    while (aboveBin < halfBins - 1 && sq[aboveBin] > halfPower) ++aboveBin;
    const tA = aboveBin > peakBin && sq[aboveBin - 1] !== sq[aboveBin]
      ? (halfPower - sq[aboveBin - 1]) / (sq[aboveBin] - sq[aboveBin - 1])
      : 0;
    const widthAbove = (aboveBin - 1 + tA - peakBin) * binHz;

    // Walk below peak until half-power, interpolate crossing
    let belowBin = peakBin;
    while (belowBin > 1 && sq[belowBin] > halfPower) --belowBin;
    const tB = belowBin < peakBin && sq[belowBin + 1] !== sq[belowBin]
      ? (halfPower - sq[belowBin + 1]) / (sq[belowBin] - sq[belowBin + 1])
      : 0;
    const widthBelow = (peakBin - belowBin - 1 + tB) * binHz;

    sumAbove += (widthAbove / f0) * (widthAbove / f0);
    sumBelow += (widthBelow / f0) * (widthBelow / f0);
    ++count;
  }

  const spread_above = count > 0 ? Math.sqrt(sumAbove / count) : 0.01;
  const spread_below = count > 0 ? Math.sqrt(sumBelow / count) : 0.01;
  return {spread_above, spread_below};
}

// Track a single partial starting from a (time, freq) seed position.
// Snaps to nearest spectral peak, then tracks forward and backward.
// Returns a partial object (with freqCurve), or null if no peak found near seed.
function trackFromSeed(frames, seedTime, seedFreq, params) {
  if (!frames || frames.length === 0) return null;

  // Find nearest frame to seedTime
  let seedFrameIdx = 0;
  let bestDt = Infinity;
  for (let i = 0; i < frames.length; ++i) {
    const dt = Math.abs(frames[i].time - seedTime);
    if (dt < bestDt) { bestDt = dt; seedFrameIdx = i; }
  }

  // Snap to nearest spectral peak within 10% freq tolerance
  const seedFrame = frames[seedFrameIdx];
  const snapTol = Math.max(seedFreq * 0.10, 50);
  let seedPeak = null, bestDist = snapTol;
  for (const pk of seedFrame.peaks) {
    const d = Math.abs(pk.freq - seedFreq);
    if (d < bestDist) { bestDist = d; seedPeak = pk; }
  }
  if (!seedPeak) return null;

  const { hopSize, sampleRate, deathAge = 5, phaseErrorWeight = 2.0 } = params;
  const trackingRatio = 0.05;
  const minTrackingHz = 20;

  // Forward pass from seed frame
  const times  = [seedFrame.time];
  const freqs  = [seedPeak.freq];
  const amps   = [seedPeak.amp];
  const phases = [seedPeak.phase];

  let fwdFreq = seedPeak.freq, fwdPhase = seedPeak.phase, fwdVel = 0, fwdAge = 0;
  for (let i = seedFrameIdx + 1; i < frames.length; ++i) {
    const predicted = fwdFreq + fwdVel;
    const predPhase = fwdPhase + 2 * Math.PI * fwdFreq * (fwdAge + 1) * hopSize / sampleRate;
    const tol = Math.max(predicted * trackingRatio, minTrackingHz);
    const { bestIdx } = findBestPeak(frames[i].peaks, new Set(), predicted, predPhase, tol, phaseErrorWeight);
    if (bestIdx >= 0) {
      const pk = frames[i].peaks[bestIdx];
      times.push(frames[i].time);
      freqs.push(pk.freq);
      amps.push(pk.amp);
      phases.push(pk.phase);
      fwdVel = pk.freq - fwdFreq;
      fwdFreq = pk.freq; fwdPhase = pk.phase; fwdAge = 0;
    } else {
      fwdAge++;
      if (fwdAge > deathAge) break;
    }
  }

  // Backward pass from seed frame
  const bwdTimes = [], bwdFreqs = [], bwdAmps = [], bwdPhases = [];
  let bwdFreq = seedPeak.freq, bwdAge = 0;
  for (let i = seedFrameIdx - 1; i >= 0; --i) {
    const tol = Math.max(bwdFreq * trackingRatio, minTrackingHz);
    let bestIdx = -1, bDist = tol;
    for (let j = 0; j < frames[i].peaks.length; ++j) {
      const d = Math.abs(frames[i].peaks[j].freq - bwdFreq);
      if (d < bDist) { bDist = d; bestIdx = j; }
    }
    if (bestIdx >= 0) {
      const pk = frames[i].peaks[bestIdx];
      bwdTimes.unshift(frames[i].time);
      bwdFreqs.unshift(pk.freq);
      bwdAmps.unshift(pk.amp);
      bwdPhases.unshift(pk.phase);
      bwdFreq = pk.freq; bwdAge = 0;
    } else {
      bwdAge++;
      if (bwdAge > deathAge) break;
    }
  }

  const allTimes  = [...bwdTimes,  ...times];
  const allFreqs  = [...bwdFreqs,  ...freqs];
  const allAmps   = [...bwdAmps,   ...amps];
  const allPhases = [...bwdPhases, ...phases];

  if (allTimes.length < 2) return null;

  const freqCurve = fitBezier(allTimes, allFreqs);
  const ac = fitBezier(allTimes, allAmps);
  freqCurve.a0 = ac.v0; freqCurve.a1 = ac.v1;
  freqCurve.a2 = ac.v2; freqCurve.a3 = ac.v3;

  return {
    times: allTimes, freqs: allFreqs, amps: allAmps, phases: allPhases,
    muted: false, freqCurve,
    replicas: { decay_alpha: 0.1, jitter: 0.05, spread_above: 0.02, spread_below: 0.02 },
  };
}

// Track an iso-energy contour starting from (seedTime, seedFreq).
// Instead of following spectral peaks, follows where energy ≈ seedEnergy.
// Useful for broad/diffuse bass regions with no detectable peaks.
// Returns a partial with large default spread, or null if seed energy is zero.
function trackIsoContour(stftCache, seedTime, seedFreq, params) {
  const { sampleRate, deathAge = 8 } = params;
  const numFrames = stftCache.getNumFrames();
  const fftSize   = stftCache.fftSize;
  const binHz     = sampleRate / fftSize;
  const halfBins  = fftSize / 2;

  // Find seed frame
  let seedFrameIdx = 0, bestDt = Infinity;
  for (let i = 0; i < numFrames; ++i) {
    const dt = Math.abs(stftCache.getFrameAtIndex(i).time - seedTime);
    if (dt < bestDt) { bestDt = dt; seedFrameIdx = i; }
  }

  const seedFrame = stftCache.getFrameAtIndex(seedFrameIdx);
  const seedBin   = Math.max(1, Math.min(halfBins - 2, Math.round(seedFreq / binHz)));
  const targetSq  = seedFrame.squaredAmplitude[seedBin];
  if (targetSq <= 0) return null;
  const targetDB  = 10 * Math.log10(targetSq);

  const trackingRatio = 0.15; // larger search window than peak tracker
  const minTrackHz    = 30;
  const maxDbDev      = 15;   // dB: declare miss if nothing within this range

  // Find bin minimizing |dB(b) - targetDB| near refBin, with mild position bias.
  function findContourBin(sq, refBin) {
    const tol     = Math.max(refBin * binHz * trackingRatio, minTrackHz);
    const tolBins = Math.ceil(tol / binHz);
    const lo = Math.max(1, refBin - tolBins);
    const hi = Math.min(halfBins - 2, refBin + tolBins);
    let bestBin = -1, bestCost = Infinity;
    for (let b = lo; b <= hi; ++b) {
      const dE = Math.abs(10 * Math.log10(Math.max(sq[b], 1e-20)) - targetDB);
      if (dE > maxDbDev) continue;
      const dPos = Math.abs(b - refBin) / Math.max(1, tolBins);
      const cost = dE + 3 * dPos; // energy match dominates, position breaks ties
      if (cost < bestCost) { bestCost = cost; bestBin = b; }
    }
    return bestBin;
  }

  const times = [seedFrame.time];
  const freqs = [seedBin * binHz];
  const amps  = [Math.sqrt(Math.max(0, targetSq))];

  // Forward pass
  let fwdBin = seedBin, fwdAge = 0;
  for (let i = seedFrameIdx + 1; i < numFrames; ++i) {
    const frame = stftCache.getFrameAtIndex(i);
    const b = findContourBin(frame.squaredAmplitude, fwdBin);
    if (b >= 0) {
      times.push(frame.time);
      freqs.push(b * binHz);
      amps.push(Math.sqrt(Math.max(0, frame.squaredAmplitude[b])));
      fwdBin = b; fwdAge = 0;
    } else { if (++fwdAge > deathAge) break; }
  }

  // Backward pass
  const bwdTimes = [], bwdFreqs = [], bwdAmps = [];
  let bwdBin = seedBin, bwdAge = 0;
  for (let i = seedFrameIdx - 1; i >= 0; --i) {
    const frame = stftCache.getFrameAtIndex(i);
    const b = findContourBin(frame.squaredAmplitude, bwdBin);
    if (b >= 0) {
      bwdTimes.unshift(frame.time);
      bwdFreqs.unshift(b * binHz);
      bwdAmps.unshift(Math.sqrt(Math.max(0, frame.squaredAmplitude[b])));
      bwdBin = b; bwdAge = 0;
    } else { if (++bwdAge > deathAge) break; }
  }

  const allTimes = [...bwdTimes, ...times];
  const allFreqs = [...bwdFreqs, ...freqs];
  const allAmps  = [...bwdAmps,  ...amps];
  if (allTimes.length < 2) return null;

  const freqCurve = fitBezier(allTimes, allFreqs);
  const ac = fitBezier(allTimes, allAmps);
  freqCurve.a0 = ac.v0; freqCurve.a1 = ac.v1;
  freqCurve.a2 = ac.v2; freqCurve.a3 = ac.v3;

  return {
    times: allTimes, freqs: allFreqs, amps: allAmps,
    phases: new Array(allTimes.length).fill(0),
    muted: false, freqCurve,
    replicas: { decay_alpha: 0.1, jitter: 0.05, spread_above: 0.15, spread_below: 0.15 },
  };
}

// Fit interpolating curve to trajectory via least-squares for inner control point values.
// Inner knots fixed at u=1/3 and u=2/3 (t = t0+dt/3, t0+2*dt/3).
// The curve passes through all 4 control points (Lagrange interpolation).
// TODO: support arbitrary number of inner control points
function fitBezier(times, values) {
  const n = times.length - 1;
  const t0 = times[0], v0 = values[0];
  const t3 = times[n], v3 = values[n];
  const dt = t3 - t0;

  if (dt <= 1e-9 || n < 2) {
    return {t0, v0, t1: t0 + dt / 3, v1: v0 + (v3 - v0) / 3, t2: t0 + 2 * dt / 3, v2: v0 + 2 * (v3 - v0) / 3, t3, v3};
  }

  // Lagrange basis with inner knots at u1=1/3, u2=2/3
  // l1(u) = u*(u-2/3)*(u-1) / ((1/3)*(1/3-2/3)*(1/3-1))  = 13.5*u*(u-2/3)*(u-1)
  // l2(u) = u*(u-1/3)*(u-1) / ((2/3)*(2/3-1/3)*(2/3-1))  = -13.5*u*(u-1/3)*(u-1)
  // l0(u) = (u-1/3)*(u-2/3)*(u-1) / ((-1/3)*(-2/3)*(-1)) = -4.5*(u-1/3)*(u-2/3)*(u-1)
  // l3(u) = u*(u-1/3)*(u-2/3) / ((2/3)*(1/3))             = 4.5*u*(u-1/3)*(u-2/3)
  // Least-squares: minimize Σ(l1*v1 + l2*v2 - target_i)^2
  // target_i = values[i] - l0*v0 - l3*v3

  let sA2 = 0, sB2 = 0, sAB = 0, sAT = 0, sBT = 0;

  for (let i = 0; i <= n; ++i) {
    const u  = (times[i] - t0) / dt;
    const l0 = -4.5 * (u - 1/3) * (u - 2/3) * (u - 1);
    const l1 =  13.5 * u * (u - 2/3) * (u - 1);
    const l2 = -13.5 * u * (u - 1/3) * (u - 1);
    const l3 =   4.5 * u * (u - 1/3) * (u - 2/3);
    const A = l1, B = l2;
    const target = values[i] - l0 * v0 - l3 * v3;
    sA2 += A * A; sB2 += B * B; sAB += A * B;
    sAT += A * target; sBT += B * target;
  }

  const det = sA2 * sB2 - sAB * sAB;
  let v1, v2;

  if (Math.abs(det) < 1e-9) {
    const idx1 = Math.round(n / 3);
    const idx2 = Math.round(2 * n / 3);
    v1 = values[idx1];
    v2 = values[idx2];
  } else {
    v1 = (sB2 * sAT - sAB * sBT) / det;
    v2 = (sA2 * sBT - sAB * sAT) / det;
  }

  return {t0, v0, t1: t0 + dt / 3, v1, t2: t0 + 2 * dt / 3, v2, t3, v3};
}