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 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 | 465x 465x 465x 465x 465x 465x 465x 465x 930x 440x 465x 465x 465x 465x 465x 465x 465x 465x 465x 465x 465x 465x 465x 1395x 465x 465x 465x 465x 60750x 60750x 60285x 465x 465x 465x 465x 465x 465x 465x 23532x 23532x 23532x 84282x 84282x 84282x 84282x 84282x 82164x 82164x 82164x 82164x 82164x 2118x 2118x 2118x 2118x 2118x 2118x 6354x 2118x 4236x 1653x 1653x 2583x 2583x 2583x 2118x | import type { mat3 } from 'gl-matrix';
import { vec3 } from 'gl-matrix';
import type { BoundsIJK, IImageVolume, Point2, Point3 } from '../../types';
import getVoxelThicknessAlongNormal from './getVoxelThicknessAlongNormal';
import type { SlabDepthCoverage } from './slabMembership';
import {
getSlabEpsilon,
getSlabHalfWidth,
isSlabDepthLowInclusive,
resolveReferencePlaneThickness,
} from './slabMembership';
export type VolumeGeometry =
| IImageVolume
| { direction: mat3; spacing: Point3; origin: Point3; dimensions: Point3 };
/**
* The slab of Rule M expressed in index space.
*
* The depth test is exactly linear in the integer voxel indices, which is what
* lets the iterator emit exact runs rather than test every voxel:
*
* ```
* depth(p) = (centre - planePoint) . n = p . g + c0
* g = Mᵀ n (the index space normal, deliberately not normalised)
* c0 = (origin - planePoint) . n
* ```
*
* Iteration nests outer -> row -> column. Derivation and axis roles:
* `docs/docs/concepts/cornerstone-tools/annotation/voxel-statistics.md`.
*/
export interface IndexSpaceSlab {
/** `g`, the unnormalised index space normal. */
g: Point3;
/** `c0`, the world depth of index (0, 0, 0). */
c0: number;
/** The voxel thickness along the normal. */
voxelThickness: number;
/** The resolved reference plane thickness. */
referencePlaneThickness: number;
/**
* A bound on `Math.abs(depth)` for any qualifying voxel. For the open
* interval this is the strict half width, already tightened by the epsilon,
* so `Math.abs(depth) < halfWidth` is the test. For the half-open interval
* this is only a bound, and {@link IndexSpaceSlab.depthLow} and
* {@link IndexSpaceSlab.depthHigh} give the exact interval.
*/
halfWidth: number;
/** The low end of the depth interval. Inclusive when `halfOpen` is set. */
depthLow: number;
/** The high end of the depth interval. Always exclusive. */
depthHigh: number;
/**
* Whether `depthLow` is inclusive, which makes the interval half-open.
*
* An open interval drops both of its endpoints, so a voxel centre that lands
* exactly on a boundary belongs to neither of two consecutive slabs. A
* half-open interval gives that voxel to exactly one of the two, which is
* what makes consecutive slabs tile. Follows from the slab's
* {@link SlabDepthCoverage}.
*/
halfOpen: boolean;
/** `argmax |g|`. Swept outermost. */
outerAxis: 0 | 1 | 2;
/** Swept in the middle loop. */
rowAxis: 0 | 1 | 2;
/** Runs are emitted along this axis, innermost. */
columnAxis: 0 | 1 | 2;
}
/**
* Computes `g = Mᵀ n`, the index space normal.
*
* The transpose, not the inverse: `M⁻¹ n` gives the wrong vector as soon as the
* spacing is anisotropic.
*/
export function getIndexSpaceNormal(
volume: VolumeGeometry,
normal: Point3
): Point3 {
const { direction, spacing } = volume;
const iVector = direction.slice(0, 3) as Point3;
const jVector = direction.slice(3, 6) as Point3;
const kVector = direction.slice(6, 9) as Point3;
return [
spacing[0] * vec3.dot(iVector, normal as vec3),
spacing[1] * vec3.dot(jVector, normal as vec3),
spacing[2] * vec3.dot(kVector, normal as vec3),
];
}
/** `argmax |g|`, with ties resolved to the lowest axis for determinism. */
export function pickOuterAxis(g: Point3): 0 | 1 | 2 {
const magnitudes = [Math.abs(g[0]), Math.abs(g[1]), Math.abs(g[2])];
let best: 0 | 1 | 2 = 0;
for (const axis of [1, 2] as const) {
if (magnitudes[axis] > magnitudes[best]) {
best = axis;
}
}
return best;
}
/**
* Builds the index space form of an annotation's slab.
*
* @param volume - Geometry of the volume being measured.
* @param planePoint - The annotation plane anchor, in world coordinates.
* @param normal - The annotation view plane normal. Must be unit length.
* @param referencePlaneThickness - The thickness in mm, or null/undefined
* to default to one voxel.
* @param options.columnAxis - Force which of the two non-outer axes carries the
* runs. Defaults to the lower-numbered of the two, so an
* acquisition-orientation volume emits runs along i for each j, matching
* row-major memory order.
* @param options.depthCoverage - Which voxels along the normal the slab
* selects: `'overlapping'`, the default, takes every voxel whose box the
* slab reaches, and `'centerInside'` takes only the voxels whose centre the
* slab contains. The half width and the depth interval both follow from it.
* See {@link SlabDepthCoverage}.
*/
export function buildIndexSpaceSlab(
volume: VolumeGeometry,
planePoint: Point3,
normal: Point3,
referencePlaneThickness?: number | null,
options: {
columnAxis?: 0 | 1 | 2;
depthCoverage?: SlabDepthCoverage;
} = {}
): IndexSpaceSlab {
const { origin } = volume;
const g = getIndexSpaceNormal(volume, normal);
const voxelThickness = getVoxelThicknessAlongNormal(volume, normal);
const thickness = resolveReferencePlaneThickness(
referencePlaneThickness,
voxelThickness
);
const c0 =
(origin[0] - planePoint[0]) * normal[0] +
(origin[1] - planePoint[1]) * normal[1] +
(origin[2] - planePoint[2]) * normal[2];
const rawHalfWidth = getSlabHalfWidth(
thickness,
voxelThickness,
options.depthCoverage
);
const epsilon = getSlabEpsilon(voxelThickness);
const halfOpen = isSlabDepthLowInclusive(options.depthCoverage);
// The epsilon shifts the half-open interval rather than narrowing it, so the
// interval keeps a width of exactly `2 * rawHalfWidth` and consecutive slabs
// still tile. Narrowing both ends, as the open interval does, would drop the
// low boundary that the half-open interval must include.
const depthLow = halfOpen ? -rawHalfWidth - epsilon : -rawHalfWidth + epsilon;
const depthHigh = rawHalfWidth - epsilon;
const halfWidth = halfOpen ? rawHalfWidth + epsilon : rawHalfWidth - epsilon;
const outerAxis = pickOuterAxis(g);
const remaining = [0, 1, 2].filter((axis) => axis !== outerAxis) as [
0 | 1 | 2,
0 | 1 | 2,
];
let columnAxis = remaining[0];
let rowAxis = remaining[1];
Iif (options.columnAxis !== undefined) {
if (options.columnAxis === outerAxis) {
throw new Error(
`columnAxis ${options.columnAxis} cannot also be the outer axis`
);
}
columnAxis = options.columnAxis;
rowAxis = remaining[0] === columnAxis ? remaining[1] : remaining[0];
}
return {
g,
c0,
voxelThickness,
referencePlaneThickness: thickness,
halfWidth,
depthLow,
depthHigh,
halfOpen,
outerAxis,
rowAxis,
columnAxis,
};
}
/**
* The world depth of a voxel index, i.e. its signed distance from the
* annotation plane along the normal.
*/
export function depthAtIndex(slab: IndexSpaceSlab, ijk: Point3): number {
const { g, c0 } = slab;
return ijk[0] * g[0] + ijk[1] * g[1] + ijk[2] * g[2] + c0;
}
/**
* The integers `x` satisfying `lo < x * coeff < hi`, as an inclusive range.
*
* Bounds are open, so an endpoint landing exactly on an integer excludes that
* integer - which is what makes a one-voxel-thick annotation select exactly one
* layer. Returns null when no integer qualifies, and `[-Infinity, Infinity]`
* when the constraint is vacuous.
*/
function integersWithProductInOpenInterval(
coeff: number,
lo: number,
hi: number
): Point2 | null {
Iif (lo >= hi) {
return null;
}
if (coeff === 0) {
// The product is always 0, so either every integer qualifies or none does.
return lo < 0 && 0 < hi ? [-Infinity, Infinity] : null;
}
const a = lo / coeff;
const b = hi / coeff;
const lowExclusive = coeff > 0 ? a : b;
const highExclusive = coeff > 0 ? b : a;
const min = Math.floor(lowExclusive) + 1;
const max = Math.ceil(highExclusive) - 1;
return min > max ? null : [min, max];
}
/**
* The integers `x` satisfying `lo <= x * coeff < hi`, as an inclusive range.
*
* The low end is closed and the high end is open, so two intervals that share a
* boundary hold every integer exactly one time between them. This is what makes
* consecutive fills tile the volume. {@link integersWithProductInOpenInterval}
* drops both boundaries instead, which leaves an integer on a boundary out of
* both intervals.
*
* Returns null when no integer qualifies, and `[-Infinity, Infinity]` when
* every integer does.
*/
function integersWithProductInHalfOpenInterval(
coeff: number,
lo: number,
hi: number
): Point2 | null {
Iif (lo >= hi) {
return null;
}
Eif (coeff === 0) {
// The product is always 0, so either every integer qualifies or none does.
return lo <= 0 && 0 < hi ? [-Infinity, Infinity] : null;
}
const a = lo / coeff;
const b = hi / coeff;
// A negative coefficient reverses the two ends, and it also moves the closed
// end from the low side to the high side.
const min = coeff > 0 ? Math.ceil(a) : Math.floor(b) + 1;
const max = coeff > 0 ? Math.ceil(b) - 1 : Math.floor(a);
return min > max ? null : [min, max];
}
/** Intersects an inclusive range with an optional inclusive clamp. */
function clampRange(range: Point2 | null, clampTo?: Point2): Point2 | null {
Iif (!range) {
return null;
}
Iif (!clampTo) {
return Number.isFinite(range[0]) && Number.isFinite(range[1])
? range
: null;
}
const min = Math.max(range[0], clampTo[0]);
const max = Math.min(range[1], clampTo[1]);
return min > max ? null : [min, max];
}
/**
* The inclusive run along `slab.columnAxis` satisfying the depth half of Rule M
* for a fixed outer and row position, or null when the column holds no voxels.
*
* @param clampTo - Inclusive bounds along the column axis, normally the volume
* dimension. Required whenever depth does not vary along the column axis,
* since the run is then unbounded.
*/
export function getDepthRun(
slab: IndexSpaceSlab,
outerIndex: number,
rowIndex: number,
clampTo?: Point2
): Point2 | null {
const {
g,
c0,
halfWidth,
depthLow,
depthHigh,
halfOpen,
outerAxis,
rowAxis,
columnAxis,
} = slab;
Iif (!(halfWidth > 0)) {
return null;
}
const base = outerIndex * g[outerAxis] + rowIndex * g[rowAxis] + c0;
const solve = halfOpen
? integersWithProductInHalfOpenInterval
: integersWithProductInOpenInterval;
return clampRange(
solve(g[columnAxis], depthLow - base, depthHigh - base),
clampTo
);
}
/**
* The inclusive range along `axis` that could contain any qualifying voxel,
* given that the other two axes are confined to `bounds`.
*
* Used to tighten the outer and row loops so they do not sweep the whole volume
* for a thin oblique slab. This is a bound rather than an exact answer: an
* index inside it may still turn out to hold no voxels, but no index outside it
* can hold any.
*/
export function getSlabAxisBound(
slab: IndexSpaceSlab,
axis: 0 | 1 | 2,
bounds: BoundsIJK,
fixed?: { axis: 0 | 1 | 2; index: number }
): Point2 | null {
const { g, c0, halfWidth } = slab;
Iif (!(halfWidth > 0)) {
return null;
}
// Slack is the range the free axes can contribute to the depth.
let slackLow = 0;
let slackHigh = 0;
let base = c0;
for (const other of [0, 1, 2] as const) {
if (other === axis) {
continue;
}
if (fixed && other === fixed.axis) {
base += fixed.index * g[other];
continue;
}
const contributions = [
bounds[other][0] * g[other],
bounds[other][1] * g[other],
];
slackLow += Math.min(contributions[0], contributions[1]);
slackHigh += Math.max(contributions[0], contributions[1]);
}
// Need |base + x * g[axis] + t| < halfWidth for some t in [slackLow, slackHigh].
return clampRange(
integersWithProductInOpenInterval(
g[axis],
-halfWidth - slackHigh - base,
halfWidth - slackLow - base
),
bounds[axis]
);
}
|