1#![forbid(unsafe_code)]
10
11use openbnct_core::GridGeometry;
12use thiserror::Error;
13
14const ALIGNMENT_TOLERANCE: f64 = 1.0e-6;
15const IDENTITY_DIRECTION: [f64; 9] = [1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0];
16
17#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
18pub enum AnatomicalPlane {
19 Axial,
20 Coronal,
21 Sagittal,
22}
23
24impl AnatomicalPlane {
25 pub const ALL: [Self; 3] = [Self::Axial, Self::Coronal, Self::Sagittal];
26
27 #[must_use]
28 pub const fn name(self) -> &'static str {
29 match self {
30 Self::Axial => "Axial",
31 Self::Coronal => "Coronal",
32 Self::Sagittal => "Sagittal",
33 }
34 }
35
36 #[must_use]
37 pub const fn edge_labels(self) -> EdgeLabels {
38 match self {
39 Self::Axial => EdgeLabels {
40 left: "R",
41 right: "L",
42 top: "A",
43 bottom: "P",
44 },
45 Self::Coronal => EdgeLabels {
46 left: "R",
47 right: "L",
48 top: "S",
49 bottom: "I",
50 },
51 Self::Sagittal => EdgeLabels {
52 left: "A",
53 right: "P",
54 top: "S",
55 bottom: "I",
56 },
57 }
58 }
59
60 const fn fixed_axis(self) -> usize {
61 match self {
62 Self::Axial => 2,
63 Self::Coronal => 1,
64 Self::Sagittal => 0,
65 }
66 }
67}
68
69#[derive(Debug, Clone, Copy, PartialEq, Eq)]
70pub struct EdgeLabels {
71 pub left: &'static str,
72 pub right: &'static str,
73 pub top: &'static str,
74 pub bottom: &'static str,
75}
76
77#[derive(Debug, Clone, Copy, PartialEq, Eq)]
78pub struct Crosshair {
79 voxel: [u32; 3],
80}
81
82impl Crosshair {
83 pub fn centered(grid: &PatientAlignedGrid) -> Self {
84 Self {
85 voxel: grid.geometry.shape.map(|extent| extent / 2),
86 }
87 }
88
89 pub fn new(grid: &PatientAlignedGrid, voxel: [u32; 3]) -> Result<Self, ViewError> {
90 validate_voxel(&grid.geometry, voxel)?;
91 Ok(Self { voxel })
92 }
93
94 #[must_use]
95 pub const fn voxel(self) -> [u32; 3] {
96 self.voxel
97 }
98
99 pub fn set_voxel(
100 &mut self,
101 grid: &PatientAlignedGrid,
102 voxel: [u32; 3],
103 ) -> Result<(), ViewError> {
104 validate_voxel(&grid.geometry, voxel)?;
105 self.voxel = voxel;
106 Ok(())
107 }
108
109 pub fn select_pixel(
110 &mut self,
111 grid: &PatientAlignedGrid,
112 view: &SliceView,
113 pixel: [u32; 2],
114 ) -> Result<(), ViewError> {
115 self.set_voxel(grid, view.voxel_at(pixel)?)
116 }
117
118 pub fn world_lps_mm(self, grid: &PatientAlignedGrid) -> Result<[f64; 3], ViewError> {
119 grid.geometry
120 .voxel_center_lps_mm(self.voxel)
121 .map_err(|error| ViewError::InvalidGeometry(error.to_string()))
122 }
123}
124
125#[derive(Debug, Clone, PartialEq)]
126pub struct PatientAlignedGrid {
127 geometry: GridGeometry,
128 voxel_count: usize,
129}
130
131impl PatientAlignedGrid {
132 pub fn new(geometry: &GridGeometry) -> Result<Self, ViewError> {
133 let voxel_count = geometry
134 .voxel_count()
135 .map_err(|error| ViewError::InvalidGeometry(error.to_string()))?;
136 if geometry
137 .direction
138 .into_iter()
139 .zip(IDENTITY_DIRECTION)
140 .any(|(observed, expected)| (observed - expected).abs() > ALIGNMENT_TOLERANCE)
141 {
142 return Err(ViewError::NotPatientAligned(geometry.direction));
143 }
144 Ok(Self {
145 geometry: geometry.clone(),
146 voxel_count,
147 })
148 }
149
150 #[must_use]
151 pub fn geometry(&self) -> &GridGeometry {
152 &self.geometry
153 }
154
155 #[must_use]
156 pub const fn voxel_count(&self) -> usize {
157 self.voxel_count
158 }
159
160 pub fn linear_index(&self, voxel: [u32; 3]) -> Result<usize, ViewError> {
161 validate_voxel(&self.geometry, voxel)?;
162 Ok(linear_index(self.geometry.shape, voxel))
163 }
164
165 pub fn slice(
166 &self,
167 plane: AnatomicalPlane,
168 crosshair: Crosshair,
169 ) -> Result<SliceView, ViewError> {
170 validate_voxel(&self.geometry, crosshair.voxel)?;
171 let dimensions = match plane {
172 AnatomicalPlane::Axial => [self.geometry.shape[0], self.geometry.shape[1]],
173 AnatomicalPlane::Coronal => [self.geometry.shape[0], self.geometry.shape[2]],
174 AnatomicalPlane::Sagittal => [self.geometry.shape[1], self.geometry.shape[2]],
175 };
176 Ok(SliceView {
177 plane,
178 fixed_index: crosshair.voxel[plane.fixed_axis()],
179 dimensions,
180 volume_shape: self.geometry.shape,
181 pixel_spacing_mm: match plane {
182 AnatomicalPlane::Axial => {
183 [self.geometry.spacing_mm[0], self.geometry.spacing_mm[1]]
184 }
185 AnatomicalPlane::Coronal => {
186 [self.geometry.spacing_mm[0], self.geometry.spacing_mm[2]]
187 }
188 AnatomicalPlane::Sagittal => {
189 [self.geometry.spacing_mm[1], self.geometry.spacing_mm[2]]
190 }
191 },
192 })
193 }
194}
195
196#[derive(Debug, Clone, Copy, PartialEq)]
197pub struct SliceView {
198 plane: AnatomicalPlane,
199 fixed_index: u32,
200 dimensions: [u32; 2],
201 volume_shape: [u32; 3],
202 pixel_spacing_mm: [f64; 2],
203}
204
205impl SliceView {
206 #[must_use]
207 pub const fn plane(self) -> AnatomicalPlane {
208 self.plane
209 }
210
211 #[must_use]
212 pub const fn fixed_index(self) -> u32 {
213 self.fixed_index
214 }
215
216 #[must_use]
217 pub const fn dimensions(self) -> [u32; 2] {
218 self.dimensions
219 }
220
221 #[must_use]
222 pub const fn pixel_spacing_mm(self) -> [f64; 2] {
223 self.pixel_spacing_mm
224 }
225
226 #[must_use]
227 pub const fn edge_labels(self) -> EdgeLabels {
228 self.plane.edge_labels()
229 }
230
231 pub fn voxel_at(self, pixel: [u32; 2]) -> Result<[u32; 3], ViewError> {
232 validate_pixel(self.dimensions, pixel)?;
233 let inverted_vertical = self.dimensions[1] - 1 - pixel[1];
234 Ok(match self.plane {
235 AnatomicalPlane::Axial => [pixel[0], pixel[1], self.fixed_index],
236 AnatomicalPlane::Coronal => [pixel[0], self.fixed_index, inverted_vertical],
237 AnatomicalPlane::Sagittal => [self.fixed_index, pixel[0], inverted_vertical],
238 })
239 }
240
241 pub fn voxel_at_fraction(self, fraction: [f32; 2]) -> Result<[u32; 3], ViewError> {
242 if fraction.iter().any(|value| !value.is_finite()) {
243 return Err(ViewError::NonFiniteScreenCoordinate);
244 }
245 let pixel = [
246 fraction_to_pixel(fraction[0], self.dimensions[0]),
247 fraction_to_pixel(fraction[1], self.dimensions[1]),
248 ];
249 self.voxel_at(pixel)
250 }
251
252 pub fn pixel_for_voxel(self, voxel: [u32; 3]) -> Result<[u32; 2], ViewError> {
253 if voxel
254 .into_iter()
255 .zip(self.volume_shape)
256 .any(|(index, extent)| index >= extent)
257 {
258 return Err(ViewError::VoxelOutOfBounds {
259 voxel,
260 shape: self.volume_shape,
261 });
262 }
263 if voxel[self.plane.fixed_axis()] != self.fixed_index {
264 return Err(ViewError::VoxelNotOnSlice {
265 voxel,
266 plane: self.plane,
267 fixed_index: self.fixed_index,
268 });
269 }
270 Ok(match self.plane {
271 AnatomicalPlane::Axial => [voxel[0], voxel[1]],
272 AnatomicalPlane::Coronal => [voxel[0], self.dimensions[1] - 1 - voxel[2]],
273 AnatomicalPlane::Sagittal => [voxel[1], self.dimensions[1] - 1 - voxel[2]],
274 })
275 }
276
277 pub fn linear_index_at(self, pixel: [u32; 2]) -> Result<usize, ViewError> {
278 let voxel = self.voxel_at(pixel)?;
279 Ok(linear_index(self.volume_shape, voxel))
280 }
281
282 pub fn extract<T: Copy>(self, values: &[T]) -> Result<Vec<T>, ViewError> {
283 let expected = self
284 .volume_shape
285 .into_iter()
286 .try_fold(1_usize, |count, extent| count.checked_mul(extent as usize));
287 if expected != Some(values.len()) {
288 return Err(ViewError::VolumeLength {
289 expected: expected.unwrap_or(usize::MAX),
290 actual: values.len(),
291 });
292 }
293 let capacity = self.dimensions[0] as usize * self.dimensions[1] as usize;
294 let mut output = Vec::with_capacity(capacity);
295 for vertical in 0..self.dimensions[1] {
296 for horizontal in 0..self.dimensions[0] {
297 output.push(values[self.linear_index_at([horizontal, vertical])?]);
298 }
299 }
300 Ok(output)
301 }
302}
303
304#[derive(Debug, Error, PartialEq)]
305pub enum ViewError {
306 #[error("invalid grid geometry: {0}")]
307 InvalidGeometry(String),
308 #[error("grid direction {0:?} is not aligned to canonical DICOM LPS axes")]
309 NotPatientAligned([f64; 9]),
310 #[error("voxel {voxel:?} is outside grid shape {shape:?}")]
311 VoxelOutOfBounds { voxel: [u32; 3], shape: [u32; 3] },
312 #[error("pixel {pixel:?} is outside slice dimensions {dimensions:?}")]
313 PixelOutOfBounds {
314 pixel: [u32; 2],
315 dimensions: [u32; 2],
316 },
317 #[error("voxel {voxel:?} is not on {plane:?} slice with fixed grid index {fixed_index}")]
318 VoxelNotOnSlice {
319 voxel: [u32; 3],
320 plane: AnatomicalPlane,
321 fixed_index: u32,
322 },
323 #[error("screen coordinate contains NaN or infinity")]
324 NonFiniteScreenCoordinate,
325 #[error("volume contains {actual} values; expected {expected}")]
326 VolumeLength { expected: usize, actual: usize },
327}
328
329fn validate_voxel(geometry: &GridGeometry, voxel: [u32; 3]) -> Result<(), ViewError> {
330 if voxel
331 .into_iter()
332 .zip(geometry.shape)
333 .any(|(index, extent)| index >= extent)
334 {
335 return Err(ViewError::VoxelOutOfBounds {
336 voxel,
337 shape: geometry.shape,
338 });
339 }
340 Ok(())
341}
342
343fn validate_pixel(dimensions: [u32; 2], pixel: [u32; 2]) -> Result<(), ViewError> {
344 if pixel
345 .into_iter()
346 .zip(dimensions)
347 .any(|(index, extent)| index >= extent)
348 {
349 return Err(ViewError::PixelOutOfBounds { pixel, dimensions });
350 }
351 Ok(())
352}
353
354fn fraction_to_pixel(fraction: f32, extent: u32) -> u32 {
355 let scaled = fraction.clamp(0.0, 1.0) * extent as f32;
356 (scaled.floor() as u32).min(extent - 1)
357}
358
359fn linear_index(shape: [u32; 3], voxel: [u32; 3]) -> usize {
360 (voxel[2] as usize * shape[1] as usize + voxel[1] as usize) * shape[0] as usize
361 + voxel[0] as usize
362}
363
364#[cfg(test)]
365mod tests {
366 use super::*;
367
368 fn geometry() -> GridGeometry {
369 GridGeometry {
370 shape: [4, 6, 8],
371 spacing_mm: [2.0, 3.0, 4.0],
372 origin_mm: [-3.0, -7.5, -14.0],
373 direction: IDENTITY_DIRECTION,
374 }
375 }
376
377 #[test]
378 fn maps_patient_aligned_view_edges_independently() {
379 let grid = PatientAlignedGrid::new(&geometry()).expect("patient-aligned grid");
380 let crosshair = Crosshair::new(&grid, [1, 2, 3]).expect("crosshair");
381
382 let axial = grid
383 .slice(AnatomicalPlane::Axial, crosshair)
384 .expect("axial");
385 assert_eq!(axial.dimensions(), [4, 6]);
386 assert_eq!(axial.pixel_spacing_mm(), [2.0, 3.0]);
387 assert_eq!(axial.voxel_at([0, 0]), Ok([0, 0, 3]));
388 assert_eq!(axial.voxel_at([3, 5]), Ok([3, 5, 3]));
389 assert_eq!(
390 axial.edge_labels(),
391 EdgeLabels {
392 left: "R",
393 right: "L",
394 top: "A",
395 bottom: "P"
396 }
397 );
398
399 let coronal = grid
400 .slice(AnatomicalPlane::Coronal, crosshair)
401 .expect("coronal");
402 assert_eq!(coronal.dimensions(), [4, 8]);
403 assert_eq!(coronal.pixel_spacing_mm(), [2.0, 4.0]);
404 assert_eq!(coronal.voxel_at([0, 0]), Ok([0, 2, 7]));
405 assert_eq!(coronal.voxel_at([3, 7]), Ok([3, 2, 0]));
406 assert_eq!(
407 coronal.edge_labels(),
408 EdgeLabels {
409 left: "R",
410 right: "L",
411 top: "S",
412 bottom: "I"
413 }
414 );
415
416 let sagittal = grid
417 .slice(AnatomicalPlane::Sagittal, crosshair)
418 .expect("sagittal");
419 assert_eq!(sagittal.dimensions(), [6, 8]);
420 assert_eq!(sagittal.pixel_spacing_mm(), [3.0, 4.0]);
421 assert_eq!(sagittal.voxel_at([0, 0]), Ok([1, 0, 7]));
422 assert_eq!(sagittal.voxel_at([5, 7]), Ok([1, 5, 0]));
423 assert_eq!(
424 sagittal.edge_labels(),
425 EdgeLabels {
426 left: "A",
427 right: "P",
428 top: "S",
429 bottom: "I"
430 }
431 );
432 }
433
434 #[test]
435 fn crosshair_round_trips_through_every_view() {
436 let grid = PatientAlignedGrid::new(&geometry()).expect("patient-aligned grid");
437 let crosshair = Crosshair::new(&grid, [1, 2, 3]).expect("crosshair");
438 for plane in AnatomicalPlane::ALL {
439 let view = grid.slice(plane, crosshair).expect("view");
440 let pixel = view
441 .pixel_for_voxel(crosshair.voxel())
442 .expect("crosshair pixel");
443 assert_eq!(view.voxel_at(pixel), Ok(crosshair.voxel()), "{plane:?}");
444 }
445 }
446
447 #[test]
448 fn extraction_uses_columns_fastest_and_superior_at_top() {
449 let grid = PatientAlignedGrid::new(&geometry()).expect("patient-aligned grid");
450 let crosshair = Crosshair::new(&grid, [1, 2, 3]).expect("crosshair");
451 let values: Vec<_> = (0..8)
452 .flat_map(|slice| {
453 (0..6)
454 .flat_map(move |row| (0..4).map(move |column| 100 * slice + 10 * row + column))
455 })
456 .collect();
457
458 let coronal = grid
459 .slice(AnatomicalPlane::Coronal, crosshair)
460 .expect("coronal");
461 let extracted = coronal.extract(&values).expect("extract");
462 assert_eq!(extracted[0], 720);
463 assert_eq!(extracted[3], 723);
464 assert_eq!(extracted[7 * 4], 20);
465 assert_eq!(extracted[7 * 4 + 3], 23);
466 }
467
468 #[test]
469 fn click_selection_updates_all_linked_coordinates() {
470 let grid = PatientAlignedGrid::new(&geometry()).expect("patient-aligned grid");
471 let mut crosshair = Crosshair::centered(&grid);
472 assert_eq!(crosshair.voxel(), [2, 3, 4]);
473 let sagittal = grid
474 .slice(AnatomicalPlane::Sagittal, crosshair)
475 .expect("sagittal");
476 crosshair
477 .select_pixel(&grid, &sagittal, [1, 2])
478 .expect("select pixel");
479 assert_eq!(crosshair.voxel(), [2, 1, 5]);
480 assert_eq!(crosshair.world_lps_mm(&grid), Ok([1.0, -4.5, 6.0]));
481 }
482
483 #[test]
484 fn rejects_anatomical_labels_for_misaligned_grid() {
485 let mut geometry = geometry();
486 geometry.direction = [0.0, 1.0, 0.0, -1.0, 0.0, 0.0, 0.0, 0.0, 1.0];
487 assert!(matches!(
488 PatientAlignedGrid::new(&geometry),
489 Err(ViewError::NotPatientAligned(_))
490 ));
491 }
492
493 #[test]
494 fn normalized_edges_map_inside_the_last_pixel() {
495 let grid = PatientAlignedGrid::new(&geometry()).expect("patient-aligned grid");
496 let crosshair = Crosshair::new(&grid, [1, 2, 3]).expect("crosshair");
497 let axial = grid
498 .slice(AnatomicalPlane::Axial, crosshair)
499 .expect("axial");
500 assert_eq!(axial.voxel_at_fraction([0.0, 0.0]), Ok([0, 0, 3]));
501 assert_eq!(axial.voxel_at_fraction([1.0, 1.0]), Ok([3, 5, 3]));
502 }
503}