Skip to content

Commit 4ab923e

Browse files
arul-trenserarul-trenser
committed
fix: freehand roi in volume viewport for oblique data
1 parent 5cb1086 commit 4ab923e

2 files changed

Lines changed: 176 additions & 88 deletions

File tree

‎packages/tools/src/tools/annotation/PlanarFreehandROITool.ts‎

Lines changed: 108 additions & 85 deletions
Original file line numberDiff line numberDiff line change
@@ -907,25 +907,6 @@ class PlanarFreehandROITool extends ContourSegmentationBaseTool {
907907

908908
const indexPoints = points.map((point) => imageData.worldToIndex(point));
909909

910-
let iMin = Number.MAX_SAFE_INTEGER;
911-
let iMax = Number.MIN_SAFE_INTEGER;
912-
let jMin = Number.MAX_SAFE_INTEGER;
913-
let jMax = Number.MIN_SAFE_INTEGER;
914-
let kMin = Number.MAX_SAFE_INTEGER;
915-
let kMax = Number.MIN_SAFE_INTEGER;
916-
917-
for (let j = 0; j < points.length; j++) {
918-
const worldPosIndex = indexPoints[j].map(Math.floor);
919-
iMin = Math.min(iMin, worldPosIndex[0]);
920-
iMax = Math.max(iMax, worldPosIndex[0]);
921-
922-
jMin = Math.min(jMin, worldPosIndex[1]);
923-
jMax = Math.max(jMax, worldPosIndex[1]);
924-
925-
kMin = Math.min(kMin, worldPosIndex[2]);
926-
kMax = Math.max(kMax, worldPosIndex[2]);
927-
}
928-
929910
// Convert from canvas_pixels ^2 to mm^2
930911
const area = polyline.getArea(canvasCoordinates) * deltaInX * deltaInY;
931912

@@ -935,74 +916,116 @@ class PlanarFreehandROITool extends ContourSegmentationBaseTool {
935916
closed
936917
);
937918

938-
// Expand bounding box
939-
const iDelta = 0.01 * (iMax - iMin);
940-
const jDelta = 0.01 * (jMax - jMin);
941-
const kDelta = 0.01 * (kMax - kMin);
942-
943-
iMin = Math.floor(iMin - iDelta);
944-
iMax = Math.ceil(iMax + iDelta);
945-
jMin = Math.floor(jMin - jDelta);
946-
jMax = Math.ceil(jMax + jDelta);
947-
kMin = Math.floor(kMin - kDelta);
948-
kMax = Math.ceil(kMax + kDelta);
949-
950-
const boundsIJK = [
951-
[iMin, iMax],
952-
[jMin, jMax],
953-
[kMin, kMax],
954-
] as [Types.Point2, Types.Point2, Types.Point2];
955-
956-
const worldPosEnd = imageData.indexToWorld([iMax, jMax, kMax]);
957-
const canvasPosEnd = viewport.worldToCanvas(worldPosEnd);
958-
959-
let curRow = 0;
960-
let intersections = [];
961-
let intersectionCounter = 0;
962-
963-
let pointsInShape;
964-
if (voxelManager) {
965-
pointsInShape = voxelManager.forEach(
966-
this.configuration.statsCalculator.statsCallback,
967-
{
968-
imageData,
969-
isInObject: (pointLPS, _pointIJK) => {
970-
let result = true;
971-
const point = viewport.worldToCanvas(pointLPS);
972-
if (point[1] != curRow) {
973-
intersectionCounter = 0;
974-
curRow = point[1];
975-
intersections = getLineSegmentIntersectionsCoordinates(
976-
canvasCoordinates,
977-
point,
978-
[canvasPosEnd[0], point[1]]
979-
);
980-
intersections.sort(
981-
(function (index) {
982-
return function (a, b) {
983-
return a[index] === b[index]
984-
? 0
985-
: a[index] < b[index]
986-
? -1
987-
: 1;
988-
};
989-
})(0)
990-
);
991-
}
992-
if (intersections.length && point[0] > intersections[0][0]) {
993-
intersections.shift();
994-
intersectionCounter++;
995-
}
996-
if (intersectionCounter % 2 === 0) {
997-
result = false;
998-
}
999-
return result;
1000-
},
1001-
boundsIJK,
1002-
returnPoints: this.configuration.storePointData,
1003-
}
919+
let pointsInShape = [];
920+
const visited = new Set<string>();
921+
922+
// Viewport-space ROI sampling for both orthogonal and oblique MPR views.
923+
//
924+
// Instead of iterating a 3D IJK bounding box, we rasterize the ROI in
925+
// canvas space, project each canvas pixel back to world space, convert
926+
// to IJK, and sample the corresponding voxel.
927+
//
928+
// This avoids the instability of voxel-driven scanline traversal in
929+
// oblique views, where IJK iteration order no longer matches canvas
930+
// row order, leading to missing voxels and invalid statistics.
931+
//
932+
// Complexity scales with projected ROI size rather than IJK volume,
933+
// making it stable and efficient for all viewport orientations.
934+
935+
const {
936+
maxX: canvasMaxX,
937+
maxY: canvasMaxY,
938+
minX: canvasMinX,
939+
minY: canvasMinY,
940+
} = math.polyline.getAABB(canvasCoordinates);
941+
942+
const startX = Math.floor(canvasMinX);
943+
const endX = Math.ceil(canvasMaxX);
944+
const startY = Math.floor(canvasMinY);
945+
const endY = Math.ceil(canvasMaxY);
946+
947+
const dimensions = imageData.getDimensions();
948+
const canvasMaxXPadded = endX + 1;
949+
950+
for (let cy = startY; cy <= endY; cy++) {
951+
// Compute all intersections of the polygon with this scanline row
952+
const intersections = getLineSegmentIntersectionsCoordinates(
953+
canvasCoordinates,
954+
[startX - 1, cy] as Types.Point2,
955+
[canvasMaxXPadded, cy] as Types.Point2
1004956
);
957+
958+
if (!intersections || intersections.length === 0) {
959+
continue;
960+
}
961+
962+
// Sort intersections by X coordinate
963+
intersections.sort((a, b) => a[0] - b[0]);
964+
965+
// Walk through intersection pairs (entry/exit)
966+
for (let i = 0; i + 1 < intersections.length; i += 2) {
967+
const xEnter = Math.ceil(intersections[i][0]);
968+
const xExit = Math.floor(intersections[i + 1][0]);
969+
970+
for (let cx = xEnter; cx <= xExit; cx++) {
971+
const canvasPoint: Types.Point2 = [cx, cy];
972+
const worldPoint = viewport.canvasToWorld(canvasPoint);
973+
const indexPoint = csUtils.transformWorldToIndex(
974+
imageData,
975+
worldPoint
976+
);
977+
978+
// Nearest-neighbor: round to closest voxel
979+
const ijkPoint: Types.Point3 = [
980+
Math.round(indexPoint[0]),
981+
Math.round(indexPoint[1]),
982+
Math.round(indexPoint[2]),
983+
];
984+
985+
// Bounds check
986+
if (
987+
ijkPoint[0] < 0 ||
988+
ijkPoint[0] >= dimensions[0] ||
989+
ijkPoint[1] < 0 ||
990+
ijkPoint[1] >= dimensions[1] ||
991+
ijkPoint[2] < 0 ||
992+
ijkPoint[2] >= dimensions[2]
993+
) {
994+
continue;
995+
}
996+
997+
const value = voxelManager.getAtIJKPoint(ijkPoint);
998+
999+
if (value === undefined || value === null) {
1000+
continue;
1001+
}
1002+
1003+
const key = `${ijkPoint[0]}_${ijkPoint[1]}_${ijkPoint[2]}`;
1004+
1005+
// avoid duplicate voxels
1006+
if (visited.has(key)) {
1007+
continue;
1008+
}
1009+
1010+
visited.add(key);
1011+
1012+
const sample = {
1013+
value: value as number,
1014+
pointLPS: worldPoint as Types.Point3,
1015+
pointIJK: ijkPoint,
1016+
};
1017+
1018+
pointsInShape.push(sample);
1019+
1020+
this.configuration.statsCalculator.statsCallback({
1021+
value: value as number,
1022+
pointLPS: worldPoint as Types.Point3,
1023+
pointIJK: ijkPoint,
1024+
});
1025+
}
1026+
}
10051027
}
1028+
10061029
const stats = this.configuration.statsCalculator.getStatistics();
10071030

10081031
const namedArea: Statistics = {

‎packages/tools/src/utilities/math/polyline/getSubPixelSpacingAndXYDirections.ts‎

Lines changed: 68 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -48,8 +48,8 @@ const getSubPixelSpacingAndXYDirections = (
4848
const kVector = direction.slice(6, 9) as Types.Point3;
4949

5050
const viewRight = vec3.create(); // Get the X direction of the viewport
51-
5251
vec3.cross(viewRight, <vec3>viewUp, <vec3>viewPlaneNormal);
52+
vec3.normalize(viewRight, viewRight);
5353

5454
const absViewRightDotI = Math.abs(vec3.dot(viewRight, iVector));
5555
const absViewRightDotJ = Math.abs(vec3.dot(viewRight, jVector));
@@ -67,7 +67,17 @@ const getSubPixelSpacingAndXYDirections = (
6767
xSpacing = volumeSpacing[2];
6868
xDir = kVector;
6969
} else {
70-
throw new Error('No support yet for oblique plane planar contours');
70+
// Oblique plane: viewRight is not aligned with any single volume axis.
71+
// Compute effective spacing by projecting the view direction onto
72+
// the volume's voxel grid.
73+
xSpacing = _computeEffectiveSpacing(
74+
viewRight as unknown as Types.Point3,
75+
iVector,
76+
jVector,
77+
kVector,
78+
volumeSpacing
79+
);
80+
xDir = Array.from(viewRight) as Types.Point3;
7181
}
7282

7383
const absViewUpDotI = Math.abs(vec3.dot(viewUp, iVector));
@@ -86,7 +96,15 @@ const getSubPixelSpacingAndXYDirections = (
8696
ySpacing = volumeSpacing[2];
8797
yDir = kVector;
8898
} else {
89-
throw new Error('No support yet for oblique plane planar contours');
99+
// Oblique plane: compute effective spacing along viewUp
100+
ySpacing = _computeEffectiveSpacing(
101+
viewUp as Types.Point3,
102+
iVector,
103+
jVector,
104+
kVector,
105+
volumeSpacing
106+
);
107+
yDir = Array.from(viewUp) as Types.Point3;
90108
}
91109

92110
spacing = [xSpacing, ySpacing];
@@ -100,4 +118,51 @@ const getSubPixelSpacingAndXYDirections = (
100118
return { spacing: subPixelSpacing, xDir, yDir };
101119
};
102120

121+
/**
122+
* Computes the effective voxel spacing along an arbitrary world-space direction.
123+
*
124+
* For a unit direction vector `d` and volume axes i, j, k with spacings
125+
* s_i, s_j, s_k, the effective spacing is:
126+
*
127+
* s_eff = 1 / sqrt( (d·i / s_i)² + (d·j / s_j)² + (d·k / s_k)² )
128+
*
129+
* This represents the world-space distance you must travel along `d` to cross
130+
* one voxel-equivalent boundary. It naturally reduces to the orthogonal case:
131+
* when `d` is aligned with axis i, d·j = d·k = 0, so s_eff = s_i.
132+
*
133+
* @param direction - A normalized world-space direction vector.
134+
* @param iVector - The volume's I direction (row).
135+
* @param jVector - The volume's J direction (column).
136+
* @param kVector - The volume's K direction (slice).
137+
* @param volumeSpacing - The [i, j, k] voxel spacings.
138+
* @returns The effective spacing along `direction`.
139+
*/
140+
function _computeEffectiveSpacing(
141+
direction: Types.Point3,
142+
iVector: Types.Point3,
143+
jVector: Types.Point3,
144+
kVector: Types.Point3,
145+
volumeSpacing: number[]
146+
): number {
147+
const dotI = vec3.dot(
148+
direction as unknown as vec3,
149+
iVector as unknown as vec3
150+
);
151+
const dotJ = vec3.dot(
152+
direction as unknown as vec3,
153+
jVector as unknown as vec3
154+
);
155+
const dotK = vec3.dot(
156+
direction as unknown as vec3,
157+
kVector as unknown as vec3
158+
);
159+
160+
const sum =
161+
(dotI * dotI) / (volumeSpacing[0] * volumeSpacing[0]) +
162+
(dotJ * dotJ) / (volumeSpacing[1] * volumeSpacing[1]) +
163+
(dotK * dotK) / (volumeSpacing[2] * volumeSpacing[2]);
164+
165+
return 1.0 / Math.sqrt(sum);
166+
}
167+
103168
export default getSubPixelSpacingAndXYDirections;

0 commit comments

Comments
 (0)