IDC2345: add tollerance parameter to generateToolState (#200)

* fix: 🐛 IDC2345: add tollerance parameter to generateToolState and use 1.e-3 as default

* fix: 🐛 implement better floating point rounding checks

* fix: 🐛 move nearlyEqual into utils and add explicit copyrights

Author: Punzo

src/adapters/Cornerstone/Segmentation_4X.js

@@ -11,7 +11,8 @@ import {
     rotateDirectionCosinesInPlane,
     flipImageOrientationPatient as flipIOP,
     flipMatrix2D,
-    rotateMatrix902D
+    rotateMatrix902D,
+    nearlyEqual
 } from "../../utilities/orientation/index.js";
 import {
     encode,
@@ -258,9 +259,10 @@ function _createSegFromImages(images, isMultiframe, options) {
  * generateToolState - Given a set of cornrstoneTools imageIds and a Segmentation buffer,
  * derive cornerstoneTools toolState and brush metadata.
  *
- * @param  {string[]} imageIds    An array of the imageIds.
- * @param  {ArrayBuffer} arrayBuffer The SEG arrayBuffer.
- * @param  {*} metadataProvider
+ * @param  {string[]} imageIds - An array of the imageIds.
+ * @param  {ArrayBuffer} arrayBuffer - The SEG arrayBuffer.
+ * @param  {*} metadataProvider.
+ * @param  {number} tolerance - checks tolerance, default value 1.e-3.
  *
  * @return {[]ArrayBuffer}a list of array buffer for each labelMap
  * @return {Object} an object from which the segment metadata can be derived
@@ -272,7 +274,8 @@ function generateToolState(
     imageIds,
     arrayBuffer,
     metadataProvider,
-    skipOverlapping = false
+    skipOverlapping = false,
+    tolerance = 1e-3
 ) {
     const dicomData = DicomMessage.readFile(arrayBuffer);
     const dataset = DicomMetaDictionary.naturalizeDataset(dicomData.dict);
@@ -340,11 +343,12 @@ function generateToolState(
         }
     }
 
-    const orientation = checkOrientation(multiframe, validOrientations, [
-        imagePlaneModule.rows,
-        imagePlaneModule.columns,
-        imageIds.length
-    ]);
+    const orientation = checkOrientation(
+        multiframe,
+        validOrientations,
+        [imagePlaneModule.rows, imagePlaneModule.columns, imageIds.length],
+        tolerance
+    );
 
     let overlapping = false;
     if (!skipOverlapping) {
@@ -353,7 +357,8 @@ function generateToolState(
             multiframe,
             imageIds,
             validOrientations,
-            metadataProvider
+            metadataProvider,
+            tolerance
         );
     }
 
@@ -401,7 +406,8 @@ function generateToolState(
         multiframe,
         imageIds,
         validOrientations,
-        metadataProvider
+        metadataProvider,
+        tolerance
     );
 
     return {
@@ -617,7 +623,8 @@ function checkSEGsOverlapping(
     multiframe,
     imageIds,
     validOrientations,
-    metadataProvider
+    metadataProvider,
+    tolerance
 ) {
     const {
         SharedFunctionalGroupsSequence,
@@ -678,7 +685,8 @@ function checkSEGsOverlapping(
         const alignedPixelDataI = alignPixelDataWithSourceData(
             pixelDataI2D,
             ImageOrientationPatientI,
-            validOrientations
+            validOrientations,
+            tolerance
         );
 
         if (!alignedPixelDataI) {
@@ -777,7 +785,8 @@ function insertOverlappingPixelDataPlanar(
     multiframe,
     imageIds,
     validOrientations,
-    metadataProvider
+    metadataProvider,
+    tolerance
 ) {
     const {
         SharedFunctionalGroupsSequence,
@@ -850,7 +859,8 @@ function insertOverlappingPixelDataPlanar(
             const alignedPixelDataI = alignPixelDataWithSourceData(
                 pixelDataI2D,
                 ImageOrientationPatientI,
-                validOrientations
+                validOrientations,
+                tolerance
             );
 
             if (!alignedPixelDataI) {
@@ -992,7 +1002,8 @@ function insertPixelDataPlanar(
     multiframe,
     imageIds,
     validOrientations,
-    metadataProvider
+    metadataProvider,
+    tolerance
 ) {
     const {
         SharedFunctionalGroupsSequence,
@@ -1027,7 +1038,8 @@ function insertPixelDataPlanar(
         const alignedPixelDataI = alignPixelDataWithSourceData(
             pixelDataI2D,
             ImageOrientationPatientI,
-            validOrientations
+            validOrientations,
+            tolerance
         );
 
         if (!alignedPixelDataI) {
@@ -1119,7 +1131,12 @@ function insertPixelDataPlanar(
     }
 }
 
-function checkOrientation(multiframe, validOrientations, sourceDataDimensions) {
+function checkOrientation(
+    multiframe,
+    validOrientations,
+    sourceDataDimensions,
+    tolerance
+) {
     const {
         SharedFunctionalGroupsSequence,
         PerFrameFunctionalGroupsSequence
@@ -1139,15 +1156,15 @@ function checkOrientation(multiframe, validOrientations, sourceDataDimensions) {
             .ImageOrientationPatient;
 
     const inPlane = validOrientations.some(operation =>
-        compareIOP(iop, operation)
+        compareIOP(iop, operation, tolerance)
     );
 
     if (inPlane) {
         return "Planar";
     }
 
     if (
-        checkIfPerpendicular(iop, validOrientations[0]) &&
+        checkIfPerpendicular(iop, validOrientations[0], tolerance) &&
         sourceDataDimensions.includes(multiframe.Rows) &&
         sourceDataDimensions.includes(multiframe.Rows)
     ) {
@@ -1159,14 +1176,15 @@ function checkOrientation(multiframe, validOrientations, sourceDataDimensions) {
 }
 
 /**
- * compareIOP - Returns true if iop1 and iop2 are equal
- * within a tollerance, dx.
+ * checkIfPerpendicular - Returns true if iop1 and iop2 are perpendicular
+ * within a tolerance.
  *
  * @param  {Number[6]} iop1 An ImageOrientationPatient array.
  * @param  {Number[6]} iop2 An ImageOrientationPatient array.
- * @return {Boolean}      True if iop1 and iop2 are equal.
+ * @param  {Number} tolerance.
+ * @return {Boolean} True if iop1 and iop2 are equal.
  */
-function checkIfPerpendicular(iop1, iop2) {
+function checkIfPerpendicular(iop1, iop2, tolerance) {
     const absDotColumnCosines = Math.abs(
         iop1[0] * iop2[0] + iop1[1] * iop2[1] + iop1[2] * iop2[2]
     );
@@ -1175,8 +1193,10 @@ function checkIfPerpendicular(iop1, iop2) {
     );
 
     return (
-        (absDotColumnCosines < dx || Math.abs(absDotColumnCosines - 1) < dx) &&
-        (absDotRowCosines < dx || Math.abs(absDotRowCosines - 1) < dx)
+        (absDotColumnCosines < tolerance ||
+            Math.abs(absDotColumnCosines - 1) < tolerance) &&
+        (absDotRowCosines < tolerance ||
+            Math.abs(absDotRowCosines - 1) < tolerance)
     );
 }
 
@@ -1343,44 +1363,50 @@ function getValidOrientations(iop) {
 /**
  * alignPixelDataWithSourceData -
  *
- * @param {Ndarray} pixelData2D The data to align.
- * @param  {Number[6]} iop The orientation of the image slice.
- * @param  {Number[8][6]} orientations   An array of valid imageOrientationPatient values.
- * @return {Ndarray}                         The aligned pixelData.
+ * @param {Ndarray} pixelData2D - The data to align.
+ * @param {Number[6]} iop - The orientation of the image slice.
+ * @param {Number[8][6]} orientations - An array of valid imageOrientationPatient values.
+ * @param {Number} tolerance.
+ * @return {Ndarray} The aligned pixelData.
  */
-function alignPixelDataWithSourceData(pixelData2D, iop, orientations) {
-    if (compareIOP(iop, orientations[0])) {
+function alignPixelDataWithSourceData(
+    pixelData2D,
+    iop,
+    orientations,
+    tolerance
+) {
+    if (compareIOP(iop, orientations[0], tolerance)) {
         return pixelData2D;
-    } else if (compareIOP(iop, orientations[1])) {
+    } else if (compareIOP(iop, orientations[1], tolerance)) {
         // Flipped vertically.
 
         // Undo Flip
         return flipMatrix2D.v(pixelData2D);
-    } else if (compareIOP(iop, orientations[2])) {
+    } else if (compareIOP(iop, orientations[2], tolerance)) {
         // Flipped horizontally.
 
         // Unfo flip
         return flipMatrix2D.h(pixelData2D);
-    } else if (compareIOP(iop, orientations[3])) {
+    } else if (compareIOP(iop, orientations[3], tolerance)) {
         //Rotated 90 degrees
 
         // Rotate back
         return rotateMatrix902D(pixelData2D);
-    } else if (compareIOP(iop, orientations[4])) {
+    } else if (compareIOP(iop, orientations[4], tolerance)) {
         //Rotated 90 degrees and fliped horizontally.
 
         // Undo flip and rotate back.
         return rotateMatrix902D(flipMatrix2D.h(pixelData2D));
-    } else if (compareIOP(iop, orientations[5])) {
+    } else if (compareIOP(iop, orientations[5], tolerance)) {
         // Rotated 90 degrees and fliped vertically
 
         // Unfo flip and rotate back.
         return rotateMatrix902D(flipMatrix2D.v(pixelData2D));
-    } else if (compareIOP(iop, orientations[6])) {
+    } else if (compareIOP(iop, orientations[6], tolerance)) {
         // Rotated 180 degrees. // TODO -> Do this more effeciently, there is a 1:1 mapping like 90 degree rotation.
 
         return rotateMatrix902D(rotateMatrix902D(pixelData2D));
-    } else if (compareIOP(iop, orientations[7])) {
+    } else if (compareIOP(iop, orientations[7], tolerance)) {
         // Rotated 270 degrees
 
         // Rotate back.
@@ -1390,24 +1416,23 @@ function alignPixelDataWithSourceData(pixelData2D, iop, orientations) {
     }
 }
 
-const dx = 1e-5;
-
 /**
  * compareIOP - Returns true if iop1 and iop2 are equal
- * within a tollerance, dx.
+ * within a tolerance.
  *
- * @param  {Number[6]} iop1 An ImageOrientationPatient array.
- * @param  {Number[6]} iop2 An ImageOrientationPatient array.
- * @return {Boolean}      True if iop1 and iop2 are equal.
+ * @param  {Number[6]} iop1 - An ImageOrientationPatient array.
+ * @param  {Number[6]} iop2 - An ImageOrientationPatient array.
+ * @param {Number} tolerance.
+ * @return {Boolean} True if iop1 and iop2 are equal.
  */
-function compareIOP(iop1, iop2) {
+function compareIOP(iop1, iop2, tolerance) {
     return (
-        Math.abs(iop1[0] - iop2[0]) < dx &&
-        Math.abs(iop1[1] - iop2[1]) < dx &&
-        Math.abs(iop1[2] - iop2[2]) < dx &&
-        Math.abs(iop1[3] - iop2[3]) < dx &&
-        Math.abs(iop1[4] - iop2[4]) < dx &&
-        Math.abs(iop1[5] - iop2[5]) < dx
+        nearlyEqual(iop1[0], iop2[0], tolerance) &&
+        nearlyEqual(iop1[1], iop2[1], tolerance) &&
+        nearlyEqual(iop1[2], iop2[2], tolerance) &&
+        nearlyEqual(iop1[3], iop2[3], tolerance) &&
+        nearlyEqual(iop1[4], iop2[4], tolerance) &&
+        nearlyEqual(iop1[5], iop2[5], tolerance)
     );
 }
 

src/utilities/orientation/index.js

@@ -4,12 +4,14 @@ import rotateDirectionCosinesInPlane from "./rotateDirectionCosinesInPlane.js";
 import rotateVectorAroundUnitVector from "./rotateVectorAroundUnitVector.js";
 import { flipMatrix2D } from "./flipMatrix2D.js";
 import rotateMatrix902D from "./rotateMatrix902D.js";
+import nearlyEqual from "./nearlyEqual.js";
 
 export {
     crossProduct3D,
     flipImageOrientationPatient,
     rotateDirectionCosinesInPlane,
     rotateVectorAroundUnitVector,
     flipMatrix2D,
-    rotateMatrix902D
+    rotateMatrix902D,
+    nearlyEqual
 };

src/utilities/orientation/nearlyEqual.js

@@ -0,0 +1,32 @@
+/**
+ * nearlyEqual - Returns true if a and b are nearly equal
+ * within a tolerance.
+ *
+ * This function logic source comes from:
+ * https://floating-point-gui.de/errors/comparison/
+ *
+ * https://floating-point-gui.de is published under
+ * the Creative Commons Attribution License (BY):
+ * http://creativecommons.org/licenses/by/3.0/
+ *
+ * @param  {Number} a
+ * @param  {Number} b
+ * @param {Number} tolerance.
+ * @return {Boolean} True if a and b are nearly equal.
+ */
+export default function nearlyEqual(a, b, epsilon) {
+    const absA = Math.abs(a);
+    const absB = Math.abs(b);
+    const diff = Math.abs(a - b);
+    if (a == b) {
+        // shortcut, handles infinities
+        return true;
+    } else if (a == 0 || b == 0 || absA + absB < Number.EPSILON) {
+        // a or b is zero or both are extremely close to it
+        // relative error is less meaningful here
+        return diff < epsilon * Number.EPSILON;
+    } else {
+        // use relative error
+        return diff / Math.min(absA + absB, Number.MAX_VALUE) < epsilon;
+    }
+}