Compose a series of monotonic piecewise functions and a color transfer function into a single color transfer function

Sources/Common/DataModel/PiecewiseFunction/helpers.d.ts

@@ -0,0 +1,38 @@
+import { vtkPiecewiseFunction } from '../PiecewiseFunction';
+import { vtkColorTransferFunction } from '../../../Rendering/Core/ColorTransferFunction';
+
+/**
+ * Compose a chain of piecewise (value) transform functions and a color
+ * transfer function into a single output color transfer function.
+ *
+ * Collects all x-positions across all transform functions and chains their
+ * outputs through to the color function, producing equivalent break points
+ * in the composed result. h(g(f(x))): g's x-values live in f's output domain
+ * and must be pulled back through f-inverse; h's x-values need g-inverse then
+ * f-inverse, and so on.
+ *
+ * Every transform function in fnList must be monotonic so that its
+ * breakpoints can be pulled back through the chain unambiguously. Without
+ * an errorThreshold the composed result interpolates linearly between the
+ * collected breakpoints, so every transform function must additionally be
+ * piecewise linear (midpoint 0.5, sharpness 0 on every segment); with an
+ * errorThreshold, shaped (non-linear) segments in any stage are captured
+ * through subdivision instead. The final-stage color transfer function is
+ * never restricted: its breakpoint colors are always reproduced exactly,
+ * and its shaped segments are approximated to within errorThreshold when
+ * one is provided (they are linearized between breakpoints otherwise).
+ *
+ * @param {vtkPiecewiseFunction[]} fnList ordered list of value transform functions, e.g. [modalityFn, voiFn, userFn]
+ * @param {vtkColorTransferFunction} colorFn final-stage color transfer function
+ * @param {vtkColorTransferFunction} outputFn function to populate with the composed result
+ * @param {Number} [errorThreshold] per-color-channel tolerance; when given, each output segment is recursively subdivided at its midpoint until linear interpolation matches the true composed color to within the threshold on every channel
+ * @param {Number} [maxSubdivisionDepth] recursion depth cap per segment (default 24). The subdivision is depth-first, so the threshold alone cannot bound the node count; tune both to control the maximum nodes produced (worst case 2^depth - 1 extra nodes per segment)
+ * @returns {Number|null} the maximum measured midpoint deviation between the composed output and the true composed chain (0 means exact; a value above errorThreshold means the depth cap stopped refinement somewhere), or null if any transform function is not monotonic, or is not piecewise linear when no errorThreshold is given (outputFn is left untouched). 0 is falsy: test the result against null for success, never its truthiness
+ */
+export function compose(
+  fnList: vtkPiecewiseFunction[],
+  colorFn: vtkColorTransferFunction,
+  outputFn: vtkColorTransferFunction,
+  errorThreshold?: number,
+  maxSubdivisionDepth?: number
+): number | null;

Sources/Common/DataModel/PiecewiseFunction/helpers.js

@@ -0,0 +1,158 @@
+import macro from 'vtk.js/Sources/macros';
+
+const { vtkWarningMacro } = macro;
+
+// Default cap on recursive subdivision in compose(): 2^-24 of a segment
+// is beyond float32 texture precision.
+const DEFAULT_MAX_SUBDIVISION_DEPTH = 24;
+
+function getColorFunctionXValues(cfun) {
+  const size = cfun.getSize();
+  const xValues = new Array(size);
+  cfun.getDataPointer().forEach((n, index) => {
+    xValues[index] = n.x;
+  });
+  return xValues;
+}
+
+function isPiecewiseLinear(fn) {
+  const v = [0, 0, 0, 0]; // [x, y, midpoint, sharpness]
+  // The last node's midpoint and sharpness shape the transition to a
+  // following node that does not exist, so it is exempt.
+  for (let i = 0; i < fn.getSize() - 1; i++) {
+    fn.getNodeValue(i, v);
+    if (v[2] !== 0.5 || v[3] !== 0) {
+      return false;
+    }
+  }
+  return true;
+}
+
+// Compose a chain of piecewise value-transform functions and a final-stage
+// color transfer function into a single output color transfer function.
+// See helpers.d.ts for the full contract.
+export function compose(
+  fnList,
+  colorFn,
+  outputFn,
+  errorThreshold = null,
+  maxSubdivisionDepth = DEFAULT_MAX_SUBDIVISION_DEPTH
+) {
+  const variedIdx = fnList.findIndex((fn) => fn.getType() === 'Varied');
+  if (variedIdx !== -1) {
+    vtkWarningMacro(
+      `compose: transform function at index ${variedIdx} is not ` +
+        'monotonic; aborting composition.'
+    );
+    return null;
+  }
+
+  if (errorThreshold == null) {
+    const nonlinearIdx = fnList.findIndex((fn) => !isPiecewiseLinear(fn));
+    if (nonlinearIdx !== -1) {
+      vtkWarningMacro(
+        `compose: transform function at index ${nonlinearIdx} is not ` +
+          'piecewise linear (midpoint 0.5, sharpness 0) and no ' +
+          'errorThreshold is given to capture its shape through ' +
+          'subdivision; aborting composition.'
+      );
+      return null;
+    }
+  }
+
+  const xSet = new Set();
+
+  // Pull a breakpoint back through fnList[idx - 1] .. fnList[0] into the
+  // original data domain. Returns null when the breakpoint cannot be
+  // inverted (findX returns null for values unattainable by the stage).
+  function pullBack(startIdx, value) {
+    let x = value;
+    for (let j = startIdx; j >= 0; j--) {
+      x = fnList[j].findX(x);
+      if (x == null) {
+        return null;
+      }
+    }
+    return x;
+  }
+
+  // Each function's breakpoint x-values live in that function's input domain.
+  fnList.forEach((fn, idx) => {
+    const data = fn.getDataPointer();
+    if (!data) return;
+    for (let i = 0; i < data.length; i += 2) {
+      const x = pullBack(idx - 1, data[i]);
+      if (x != null) xSet.add(x);
+    }
+  });
+
+  // Also reverse compute from x-values of the final-stage color transfer function,
+  // and add those to our xSet so that we don't miss any break points defined
+  // within the color transfer function.
+  const colorXs = getColorFunctionXValues(colorFn);
+  colorXs.forEach((x) => {
+    const t = pullBack(fnList.length - 1, x);
+    if (t != null) xSet.add(t);
+  });
+
+  // Now use the gathered x values to propogate through the entire set of
+  // functions to determine the final color values for each, and add them
+  // as nodes to our new color transfer function.
+  const xs = Array.from(xSet).sort((a, b) => a - b);
+
+  function composedColor(x) {
+    const finalScalar = fnList.reduce((val, fn) => fn.getValue(val), x);
+    const rgb = [0, 0, 0];
+    colorFn.getColor(finalScalar, rgb);
+    return rgb;
+  }
+
+  // Recursively split [x0, x1] at its midpoint until linear interpolation
+  // of the endpoint colors matches the true composed color on every
+  // channel to within errorThreshold, and return the maximum midpoint
+  // deviation left in the emitted segments. The depth cap guarantees
+  // termination around discontinuities (colorFn sharpness ~1), where the
+  // midpoint test can never pass; a residual above errorThreshold means
+  // the cap stopped refinement.
+  function subdivide(x0, c0, x1, c1, depth) {
+    const xm = 0.5 * (x0 + x1);
+    if (xm <= x0 || xm >= x1) {
+      // Float precision exhausted; nothing left to measure or refine.
+      return 0;
+    }
+    const cm = composedColor(xm);
+    const deviation = cm.reduce(
+      (max, c, i) => Math.max(max, Math.abs(c - 0.5 * (c0[i] + c1[i]))),
+      0
+    );
+    if (depth <= 0 || (errorThreshold != null && deviation <= errorThreshold)) {
+      return deviation;
+    }
+    // The midpoint becomes a node (exact there), so the residual is
+    // whatever the two halves leave behind.
+    const leftError = subdivide(x0, c0, xm, cm, depth - 1);
+    outputFn.addRGBPoint(xm, cm[0], cm[1], cm[2]);
+    const rightError = subdivide(xm, cm, x1, c1, depth - 1);
+    return Math.max(leftError, rightError);
+  }
+
+  const colors = xs.map(composedColor);
+  // Without a threshold, depth 0 measures each segment's midpoint
+  // deviation without subdividing, so the returned error is meaningful
+  // in both modes.
+  const depth = errorThreshold != null ? maxSubdivisionDepth : 0;
+  let maxError = 0;
+  outputFn.removeAllPoints();
+  xs.forEach((x, i) => {
+    const rgb = colors[i];
+    outputFn.addRGBPoint(x, rgb[0], rgb[1], rgb[2]);
+    if (i + 1 < xs.length) {
+      maxError = Math.max(
+        maxError,
+        subdivide(x, rgb, xs[i + 1], colors[i + 1], depth)
+      );
+    }
+  });
+
+  return maxError;
+}

Sources/Common/DataModel/PiecewiseFunction/index.d.ts

@@ -85,12 +85,26 @@ export interface vtkPiecewiseFunction extends vtkObject {
   getFirstNonZeroValue(): number;
 
   /**
-   * For the node specified by index, set/get the location (X), value (Y),
+   * Inverse of getValue(): given a value y, returns an x such that
+   * getValue(x) === y, bisecting the first segment whose endpoint values
+   * bracket y using forward evaluation, so midpoint and sharpness are
+   * honored. For a discontinuous segment (sharpness ~1) the jump location
+   * is returned. When y is outside the function's output range: with
+   * clamping on, returns the endpoint x whose node holds the nearest
+   * extremum; returns null if that extremum is only attained at an
+   * interior node (no x clamps to y) or if clamping is off.
+   * @param {Number} y
+   */
+  findX(y: number): number | null;
+
+  /**
+   * For the node specified by index, get the location (X), value (Y),
    * midpoint, and sharpness values at the node.
    * @param {Number} index
    * @param val
+   * @returns -1 if index is out of range, returns 1 otherwise.
    */
-  getNodeValue(index: number, val: any[]): void;
+  getNodeValue(index: number, val: number[]): -1 | 1;
 
   /**
    * Returns the min and max node locations of the function.
@@ -169,11 +183,13 @@ export interface vtkPiecewiseFunction extends vtkObject {
   setClamping(clamping: boolean): boolean;
 
   /**
-   *
+   * For the node specified by index, set the location (X), value (Y),
+   * midpoint, and sharpness values at the node.
    * @param {Number} index
    * @param val
+   * @returns -1 if index is out of range, returns 1 otherwise.
    */
-  setNodeValue(index: number, val: any[]): number;
+  setNodeValue(index: number, val: number[]): -1 | 1;
 
   /**
    *

Sources/Common/DataModel/PiecewiseFunction/index.js

@@ -19,8 +19,9 @@ function vtkPiecewiseFunction(publicAPI, model) {
   //    1 : NonDecreasing   (Always increasing or zero slope)
   //    2 : NonIncreasing   (Always decreasing or zero slope)
   //    3 : Varied          (Contains both decreasing and increasing slopes)
-  //    4 : Unknown         (Error condition)
   //
+  // Unlike VTK/C++ there is no 'Unknown' error state: functionType can
+  // only ever be one of the four values above.
   publicAPI.getType = () => {
     let value;
     let prevValue = 0.0;
@@ -338,6 +339,75 @@ function vtkPiecewiseFunction(publicAPI, model) {
     return table[0];
   };
 
+  // Inverse of getValue(). See index.d.ts for the full contract.
+  publicAPI.findX = (y) => {
+    const { nodes } = model;
+    for (let i = 0; i < nodes.length - 1; i++) {
+      const { x: x0, y: y0 } = nodes[i];
+      const { x: x1, y: y1 } = nodes[i + 1];
+      if (y === y0) {
+        return x0;
+      }
+      if (y === y1) {
+        return x1;
+      }
+      const minY = Math.min(y0, y1);
+      const maxY = Math.max(y0, y1);
+      if (y > minY && y < maxY) {
+        // Within a segment, forward evaluation is monotonic from y0 to
+        // y1 regardless of midpoint and sharpness, so bisect until
+        // floating-point convergence.
+        const increasing = y1 > y0;
+        let lo = x0;
+        let hi = x1;
+        let mid = 0.5 * (lo + hi);
+        while (lo < mid && mid < hi) {
+          const value = publicAPI.getValue(mid);
+          if (value === y) {
+            return mid;
+          }
+          if (value < y === increasing) {
+            lo = mid;
+          } else {
+            hi = mid;
+          }
+          mid = 0.5 * (lo + hi);
+        }
+        return mid;
+      }
+    }
+    if (model.clamping && nodes.length > 0) {
+      let minY = Infinity;
+      let maxY = -Infinity;
+      nodes.forEach((node) => {
+        minY = Math.min(minY, node.y);
+        maxY = Math.max(maxY, node.y);
+      });
+      const first = nodes[0];
+      const last = nodes[nodes.length - 1];
+      if (y >= maxY) {
+        if (first.y === maxY) {
+          return first.x;
+        }
+        if (last.y === maxY) {
+          return last.x;
+        }
+        return null;
+      }
+      if (y <= minY) {
+        if (first.y === minY) {
+          return first.x;
+        }
+        if (last.y === minY) {
+          return last.x;
+        }
+        return null;
+      }
+    }
+
+    return null;
+  };
+
   // Remove all points outside the range, and make sure a point
   // exists at each end of the range. Used as a convenience method
   // for transfer function editors