All files / core/src/utilities/voxelSlab/shapes createRectangleShape.ts

75.67% Statements 28/37
50% Branches 14/28
80% Functions 4/5
75.67% Lines 28/37

Press n or j to go to the next uncovered block, b, p or k for the previous block.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167                                            137x 6x                                                                                                           3x   3x     3x       3x 3x   3x 3x       3x       3x       3x       3x                                                 288x 288x 288x   288x   288x   465x       465x   465x       465x 140x       148x 148x 148x       3x     3x      
import type { Point2, Point3 } from '../../../types';
import type { IndexSpaceSlab, VolumeGeometry } from '../indexSpaceSlab';
import getVoxelThicknessAlongNormal from '../getVoxelThicknessAlongNormal';
import { projectPointOntoPlane } from '../slabMembership';
import type { PlaneBasis, RealRange, VoxelSlabShape } from './shapeGeometry';
import {
  UNBOUNDED,
  createColumnLineResolver,
  createPlaneBasis,
  SHAPE_BOUNDARY_EPSILON,
  intersectRanges,
  solveAbsLinearLeq,
  toIntegerRun,
} from './shapeGeometry';
 
/**
 * Expands a half extent so a voxel centre on the face counts as inside.
 *
 * The slack scales with the half extent, with a floor of one unit so a
 * sub-millimetre extent keeps a usable tolerance. `createPolylineShape` floors
 * its own slack the same way.
 */
const expand = (halfExtent: number) =>
  halfExtent + Math.max(halfExtent, 1) * SHAPE_BOUNDARY_EPSILON;
 
export interface RectangleShapeOptions {
  /**
   * The volume being measured. Pass an `IImageVolume`, or any object that has
   * `direction`, `spacing`, `origin` and `dimensions`. The structural form
   * lets a caller that holds no cached volume, such as a test, use this code.
   */
  volume: VolumeGeometry;
  /** The annotation plane anchor. Defines the plane's depth. */
  planePoint: Point3;
  /** The annotation view plane normal. Unit length. */
  viewPlaneNormal: Point3;
  /**
   * The rectangle centre in world coordinates. For the planar form it is
   * projected onto the annotation plane.
   */
  centerWorld: Point3;
  /**
   * Direction of the major axis. Need not be unit length, and need not already
   * lie in the plane - its component along the normal is removed.
   */
  majorAxis: Point3;
  /** Half length along `majorAxis`, in mm. */
  majorHalfLength: number;
  /** Half length perpendicular to `majorAxis` within the plane, in mm. */
  minorHalfLength: number;
  /**
   * Half length along the normal, in mm. Omit for a flat rectangle lying in the
   * plane; supply it for a box reaching out of the plane on both sides.
   */
  depthHalfLength?: number;
}
 
/**
 * A rectangle lying in the annotation plane, or a box centred on it.
 *
 * A third half extent promotes the flat rectangle to a solid.
 * `createEllipseShape` takes the same options, so the two are interchangeable.
 * The depth extent is widened by half a voxel; the in-plane extents are not.
 *
 * See `docs/docs/concepts/cornerstone-tools/annotation/voxel-statistics.md`.
 */
export function createRectangleShape(
  options: RectangleShapeOptions
): VoxelSlabShape {
  const {
    volume,
    planePoint,
    viewPlaneNormal: normal,
    majorAxis,
    majorHalfLength,
    minorHalfLength,
    depthHalfLength,
  } = options;
 
  Iif (!(majorHalfLength > 0) || !(minorHalfLength > 0)) {
    throw new Error('Rectangle half lengths must be positive');
  }
  Iif (depthHalfLength !== undefined && !(depthHalfLength > 0)) {
    throw new Error('Rectangle depthHalfLength must be positive when supplied');
  }
 
  const basis: PlaneBasis = createPlaneBasis(normal, majorAxis);
  const isSolid = depthHalfLength !== undefined;
 
  const voxelThickness = getVoxelThicknessAlongNormal(volume, basis.n);
  const depthExtent = isSolid
    ? (depthHalfLength as number) + voxelThickness / 2
    : 0;
 
  const center = isSolid
    ? ([...options.centerWorld] as Point3)
    : projectPointOntoPlane(options.centerWorld, planePoint, basis.n);
 
  const constraints: { axis: Point3; halfExtent: number }[] = [
    { axis: basis.u, halfExtent: expand(majorHalfLength) },
    { axis: basis.v, halfExtent: expand(minorHalfLength) },
  ];
  Iif (isSolid) {
    constraints.push({ axis: basis.n, halfExtent: expand(depthExtent) });
  }
 
  const resolveColumnLine = createColumnLineResolver(volume, basis.n);
 
  function containsPoint(point: Point3): boolean {
    const tested = isSolid
      ? point
      : projectPointOntoPlane(point, planePoint, basis.n);
 
    for (const { axis, halfExtent } of constraints) {
      const offset =
        (tested[0] - center[0]) * axis[0] +
        (tested[1] - center[1]) * axis[1] +
        (tested[2] - center[2]) * axis[2];
      if (Math.abs(offset) > halfExtent) {
        return false;
      }
    }
    return true;
  }
 
  function* getRuns(
    outerIndex: number,
    rowIndex: number,
    _depthRun: Point2,
    slab: IndexSpaceSlab
  ): Generator<Point2, void, undefined> {
    const line = resolveColumnLine(slab, outerIndex, rowIndex);
    const lineBase = isSolid ? line.base : line.projectedBase;
    const lineStep = isSolid ? line.step : line.projectedStep;
 
    let range: RealRange | null = UNBOUNDED;
 
    for (const { axis, halfExtent } of constraints) {
      const offset =
        (lineBase[0] - center[0]) * axis[0] +
        (lineBase[1] - center[1]) * axis[1] +
        (lineBase[2] - center[2]) * axis[2];
      const slope =
        lineStep[0] * axis[0] + lineStep[1] * axis[1] + lineStep[2] * axis[2];
 
      range = intersectRanges(
        range,
        solveAbsLinearLeq(offset, slope, halfExtent)
      );
      if (!range) {
        return;
      }
    }
 
    const run = toIntegerRun(range);
    Eif (run) {
      yield run;
    }
  }
 
  return {
    containsPoint,
    getRuns,
    getRequiredThickness: () => (isSolid ? depthExtent * 2 : 0),
  };
}