fix: 🐛 Correctly orient non-axial volumes (#69)

* fix: 🐛 Correctly orient non-axial volumes
* Sort data if positions are mixed.

Author: JamesAPetts

examples/VTKLoadImageDataExample.js

@@ -44,6 +44,15 @@ function createActorMapper(imageData) {
   };
 }
 
+function shuffleImageNames(imageNames) {
+  for (let index = imageNames.length - 1; index > 0; index--) {
+    const randomIndex = Math.floor(Math.random() * index);
+    const temporaryValue = imageNames[index];
+    imageNames[index] = imageNames[randomIndex];
+    imageNames[randomIndex] = temporaryValue;
+  }
+}
+
 function getImageIds() {
   const ROOT_URL =
     'https://s3.amazonaws.com/IsomicsPublic/SampleData/QIN-H%2BN-0139/PET-sorted/';
@@ -54,6 +63,9 @@ function getImageIds() {
     imageNames.push('PET_HeadNeck_0-' + i + '.dcm');
   }
 
+  // Shuffle the image names to test that src/lib/data/sortDatasetsByImagePosition.js works.
+  shuffleImageNames(imageNames);
+
   return imageNames.map(name => `dicomweb:${ROOT_URL}${name}`);
 }
 

src/lib/data/getSliceIndex.js

@@ -1,32 +1,159 @@
-import cornerstone from 'cornerstone-core';
-
-// insert the slice at the z index location.
+/**
+ *
+ * @param {Object} imageData - The vtkImageData
+ * @param {*} sliceIndex - The index of the slice you are inserting.
+ * @param {*} acquistionDirection - The acquistion direction of the slice.
+ * @param {*} image The cornerstone image to pull pixel data data from.
+ * @param {*} modality The modality of the image.
+ * @param {*} modalitySpecificScalingParameters Specific scaling paramaters for this modality. E.g. Patient weight.
+ */
 export default function insertSlice(
   imageData,
   sliceIndex,
+  acquistionDirection,
   image,
   modality,
   modalitySpecificScalingParameters
 ) {
-  const pixels = image.getPixelData();
-  const { rows, columns } = image;
   const scalars = imageData.getPointData().getScalars();
   const scalarData = scalars.getData();
-  const sliceLength = pixels.length;
-  let pixelIndex = 0;
 
   const scalingFunction = _getScalingFunction(
     modality,
     image,
     modalitySpecificScalingParameters
   );
 
+  const vtkImageDimensions = imageData.getDimensions();
+
+  let minAndMax;
+
+  switch (acquistionDirection) {
+    case 'coronal':
+      minAndMax = insertCoronalSlice(
+        image,
+        sliceIndex,
+        vtkImageDimensions,
+        scalarData,
+        scalingFunction
+      );
+      break;
+    case 'sagittal':
+      minAndMax = insertSagittalSlice(
+        image,
+        sliceIndex,
+        vtkImageDimensions,
+        scalarData,
+        scalingFunction
+      );
+      break;
+    case 'axial':
+      minAndMax = insertAxialSlice(
+        image,
+        sliceIndex,
+        scalarData,
+        scalingFunction
+      );
+      break;
+  }
+
+  return minAndMax;
+}
+
+/**
+ *
+ * @param {object} image The cornerstone image to pull pixel data data from.
+ * @param {number} xIndex The x index of axially oriented vtk volume to put the sagital slice.
+ * @param {number[]} vtkImageDimensions The dimensions of the axially oriented vtk volume.
+ * @param {number[]} scalarData The data array for the axially oriented vtk volume.
+ * @param {function} scalingFunction The modality specific scaling function.
+ *
+ * @returns {object} The min and max pixel values in the inserted slice.
+ */
+function insertSagittalSlice(
+  image,
+  xIndex,
+  vtkImageDimensions,
+  scalarData,
+  scalingFunction
+) {
+  const pixels = image.getPixelData();
+  const { rows, columns } = image;
+
+  let pixelIndex = 0;
+  let max = scalingFunction(pixels[pixelIndex]);
+  let min = max;
+
+  const vtkImageDimensionsX = vtkImageDimensions[0];
+  const vtkImageDimensionsY = vtkImageDimensions[1];
+  const vtkImageDimensionsZ = vtkImageDimensions[2];
+
+  const axialSliceLength = vtkImageDimensionsX * vtkImageDimensionsY;
+
+  for (let row = 0; row < rows; row++) {
+    for (let col = 0; col < columns; col++) {
+      const yPos = vtkImageDimensionsY - col;
+      const zPos = vtkImageDimensionsZ - row;
+
+      const destIdx =
+        zPos * axialSliceLength + yPos * vtkImageDimensionsX + xIndex;
+
+      const pixel = pixels[pixelIndex];
+      const pixelValue = scalingFunction(pixel);
+
+      if (pixelValue > max) {
+        max = pixelValue;
+      } else if (pixelValue < min) {
+        min = pixelValue;
+      }
+
+      scalarData[destIdx] = pixelValue;
+      pixelIndex++;
+    }
+  }
+
+  return { min, max };
+}
+
+/**
+ *
+ * @param {object} image The cornerstone image to pull pixel data data from.
+ * @param {number} yIndex The y index of axially oriented vtk volume to put the coronal slice.
+ * @param {number[]} vtkImageDimensions The dimensions of the axially oriented vtk volume.
+ * @param {number[]} scalarData The data array for the axially oriented vtk volume.
+ * @param {function} scalingFunction The modality specific scaling function.
+ *
+ * @returns {object} The min and max pixel values in the inserted slice.
+ */
+function insertCoronalSlice(
+  image,
+  yIndex,
+  vtkImageDimensions,
+  scalarData,
+  scalingFunction
+) {
+  const pixels = image.getPixelData();
+  const { rows, columns } = image;
+
+  let pixelIndex = 0;
   let max = scalingFunction(pixels[pixelIndex]);
   let min = max;
 
-  for (let row = 0; row <= rows; row++) {
-    for (let col = 0; col <= columns; col++) {
-      const destIdx = pixelIndex + sliceIndex * sliceLength;
+  const vtkImageDimensionsX = vtkImageDimensions[0];
+  const vtkImageDimensionsY = vtkImageDimensions[1];
+  const vtkImageDimensionsZ = vtkImageDimensions[2];
+
+  const axialSliceLength = vtkImageDimensionsX * vtkImageDimensionsY;
+
+  for (let row = 0; row < rows; row++) {
+    for (let col = 0; col < columns; col++) {
+      const xPos = col;
+      const yPos = yIndex;
+      const zPos = vtkImageDimensionsZ - row;
+
+      const destIdx =
+        zPos * axialSliceLength + yPos * vtkImageDimensionsX + xPos;
+
       const pixel = pixels[pixelIndex];
       const pixelValue = scalingFunction(pixel);
 
@@ -44,6 +171,40 @@ export default function insertSlice(
   return { min, max };
 }
 
+/**
+ *
+ * @param {object} image The cornerstone image to pull pixel data data from.
+ * @param {number} zIndex The z index of axially oriented vtk volume to put the axial slice.
+ * @param {number[]} scalarData The data array for the axially oriented vtk volume.
+ * @param {function} scalingFunction The modality specific scaling function.
+ *
+ * @returns {object} The min and max pixel values in the inserted slice.
+ */
+function insertAxialSlice(image, zIndex, scalarData, scalingFunction) {
+  const pixels = image.getPixelData();
+  const sliceLength = pixels.length;
+
+  let pixelIndex = 0;
+  let max = scalingFunction(pixels[pixelIndex]);
+  let min = max;
+
+  for (let pixelIndex = 0; pixelIndex < pixels.length; pixelIndex++) {
+    const destIdx = pixelIndex + zIndex * sliceLength;
+    const pixel = pixels[pixelIndex];
+    const pixelValue = scalingFunction(pixel);
+
+    if (pixelValue > max) {
+      max = pixelValue;
+    } else if (pixelValue < min) {
+      min = pixelValue;
+    }
+
+    scalarData[destIdx] = pixelValue;
+  }
+
+  return { min, max };
+}
+
 function _getScalingFunction(
   modality,
   image,

src/lib/data/insertSlice.js

@@ -0,0 +1,44 @@
+import { Vector3 } from 'cornerstone-math';
+
+export default function sortDatasetsByImagePosition(
+  scanAxisNormal,
+  imageMetaDataMap
+) {
+  // See https://github.com/dcmjs-org/dcmjs/blob/4849ed50db8788741c2773b3d9c75cc52441dbcb/src/normalizers.js#L167
+  // TODO: Find a way to make this code generic?
+
+  const datasets = Array.from(imageMetaDataMap.values());
+  const referenceDataset = datasets[0];
+
+  const refIppVec = new Vector3(...referenceDataset.imagePositionPatient);
+
+  const distanceDatasetPairs = datasets.map(function(dataset) {
+    const ippVec = new Vector3(...dataset.imagePositionPatient);
+    const positionVector = refIppVec.clone().sub(ippVec);
+    const distance = positionVector.dot(scanAxisNormal);
+
+    return {
+      distance,
+      dataset,
+    };
+  });
+
+  distanceDatasetPairs.sort(function(a, b) {
+    return b.distance - a.distance;
+  });
+
+  const sortedDatasets = distanceDatasetPairs.map(a => a.dataset);
+  const distances = distanceDatasetPairs.map(a => a.distance);
+
+  // TODO: The way we calculate spacing determines how the volume shows up if
+  // we have missing slices.
+  // - Should we just bail out for now if missing slices are present?
+  //const spacing = mean(diff(distances));
+  const spacing = Math.abs(distances[1] - distances[0]);
+
+  return {
+    spacing,
+    origin: distanceDatasetPairs[0].dataset.imagePositionPatient,
+    sortedDatasets,
+  };
+}