NRRDLoader handling of coronal and sagittal oriented scans. (#21962)

* NRRDLoader.js now identifies anchor points from within the "space directions" field within an NRRD file to orient and transform the volume structures properly in 3D space, fixing the previously broken sagittal and coronal oriented scans.

* Improved initial solution. Retained most of the core-logic and changed only how the spacing (and thus, the size) is calculated by changing its read target. This is done by storing axis order as a new member of the Volume objects

* Better semantics for spacing.

* Removed extra variables.

Co-authored-by: Geoffrey Lee <[email protected]>

examples/js/loaders/NRRDLoader.js

@@ -359,40 +359,81 @@
 			volume.windowLow = min;
 			volume.windowHigh = max; // get the image dimensions
 
+
+			// get the image dimensions
 			volume.dimensions = [ headerObject.sizes[ 0 ], headerObject.sizes[ 1 ], headerObject.sizes[ 2 ] ];
 			volume.xLength = volume.dimensions[ 0 ];
 			volume.yLength = volume.dimensions[ 1 ];
-			volume.zLength = volume.dimensions[ 2 ]; // spacing
+			volume.zLength = volume.dimensions[ 2 ];
+
+			// Identify axis order in the space-directions matrix from the header if possible.
+			if (headerObject.vectors) {
+				const xIndex = headerObject.vectors.findIndex(vector => vector[0] !== 0);
+				const yIndex = headerObject.vectors.findIndex(vector => vector[1] !== 0);
+				const zIndex = headerObject.vectors.findIndex(vector => vector[2] !== 0);
+
+				const axisOrder = [];
+				axisOrder[xIndex] = 'x';
+				axisOrder[yIndex] = 'y';
+				axisOrder[zIndex] = 'z';
+				volume.axisOrder = axisOrder;
+			}
+			else {
+				volume.axisOrder = ['x', 'y', 'z'];
+			}
 
-			const spacingX = new THREE.Vector3( headerObject.vectors[ 0 ][ 0 ], headerObject.vectors[ 0 ][ 1 ], headerObject.vectors[ 0 ][ 2 ] ).length();
-			const spacingY = new THREE.Vector3( headerObject.vectors[ 1 ][ 0 ], headerObject.vectors[ 1 ][ 1 ], headerObject.vectors[ 1 ][ 2 ] ).length();
-			const spacingZ = new THREE.Vector3( headerObject.vectors[ 2 ][ 0 ], headerObject.vectors[ 2 ][ 1 ], headerObject.vectors[ 2 ][ 2 ] ).length();
-			volume.spacing = [ spacingX, spacingY, spacingZ ]; // Create IJKtoRAS matrix
+			// spacing
+			const spacingX = new THREE.Vector3().fromArray( headerObject.vectors[ 0 ] ).length();
+			const spacingY = new THREE.Vector3().fromArray( headerObject.vectors[ 1 ] ).length();
+			const spacingZ = new THREE.Vector3().fromArray( headerObject.vectors[ 2 ] ).length();
+			volume.spacing = [ spacingX, spacingY, spacingZ ];
 
+			// Create IJKtoRAS matrix
 			volume.matrix = new THREE.Matrix4();
-			let _spaceX = 1;
-			let _spaceY = 1;
-			const _spaceZ = 1;
 
-			if ( headerObject.space == 'left-posterior-superior' ) {
+			const transitionMatrix = new THREE.Matrix4();
+
+			if ( headerObject.space === 'left-posterior-superior' ) {
 
-				_spaceX = - 1;
-				_spaceY = - 1;
+				transitionMatrix.set(
+					- 1, 0, 0, 0,
+					0, -1, 0, 0,
+					0, 0, 1, 0,
+					0, 0, 0, 1 );
 
 			} else if ( headerObject.space === 'left-anterior-superior' ) {
 
-				_spaceX = - 1;
+				transitionMatrix.set(
+					1, 0, 0, 0,
+					0, 1, 0, 0,
+					0, 0, -1, 0,
+					0, 0, 0, 1 );
 
 			}
 
+
 			if ( ! headerObject.vectors ) {
 
-				volume.matrix.set( _spaceX, 0, 0, 0, 0, _spaceY, 0, 0, 0, 0, _spaceZ, 0, 0, 0, 0, 1 );
+				volume.matrix.set(
+					1, 0, 0, 0,
+					0, 1, 0, 0,
+					0, 0, 1, 0,
+					0, 0, 0, 1 );
 
 			} else {
 
 				const v = headerObject.vectors;
-				volume.matrix.set( _spaceX * v[ 0 ][ 0 ], _spaceX * v[ 1 ][ 0 ], _spaceX * v[ 2 ][ 0 ], 0, _spaceY * v[ 0 ][ 1 ], _spaceY * v[ 1 ][ 1 ], _spaceY * v[ 2 ][ 1 ], 0, _spaceZ * v[ 0 ][ 2 ], _spaceZ * v[ 1 ][ 2 ], _spaceZ * v[ 2 ][ 2 ], 0, 0, 0, 0, 1 );
+
+				const ijk_to_transition = ( new THREE.Matrix4() ).set(
+					v[ 0 ][ 0 ], v[ 1 ][ 0 ], v[ 2 ][ 0 ], 0,
+					v[ 0 ][ 1 ], v[ 1 ][ 1 ], v[ 2 ][ 1 ], 0,
+					v[ 0 ][ 2 ], v[ 1 ][ 2 ], v[ 2 ][ 2 ], 0,
+					0, 0, 0, 1
+				)
+
+				const transition_to_ras = (new THREE.Matrix4()).multiplyMatrices( ijk_to_transition, transitionMatrix );
+
+				volume.matrix = transition_to_ras;
 
 			}
 

examples/js/misc/Volume.js

@@ -31,6 +31,13 @@
      */
 
 			this.zLength = Number( zLength ) || 1;
+
+			/**
+			 * @member {Array<string>} The order of the Axis dictated by the NRRD header
+			 */
+
+			this.axisOrder = [ 'x', 'y', 'z' ];
+
 			/**
      * @member {TypedArray} data Data of the volume
      */
@@ -306,8 +313,8 @@
 					axisInIJK.set( 1, 0, 0 );
 					firstDirection.set( 0, 0, - 1 );
 					secondDirection.set( 0, - 1, 0 );
-					firstSpacing = this.spacing[ 2 ];
-					secondSpacing = this.spacing[ 1 ];
+					firstSpacing = this.spacing[ this.axisOrder.indexOf('z') ];
+					secondSpacing = this.spacing[ this.axisOrder.indexOf('y') ];
 					IJKIndex = new THREE.Vector3( RASIndex, 0, 0 );
 					planeMatrix.multiply( new THREE.Matrix4().makeRotationY( Math.PI / 2 ) );
 					positionOffset = ( volume.RASDimensions[ 0 ] - 1 ) / 2;
@@ -318,8 +325,8 @@
 					axisInIJK.set( 0, 1, 0 );
 					firstDirection.set( 1, 0, 0 );
 					secondDirection.set( 0, 0, 1 );
-					firstSpacing = this.spacing[ 0 ];
-					secondSpacing = this.spacing[ 2 ];
+					firstSpacing = this.spacing[ this.axisOrder.indexOf('x') ];
+					secondSpacing = this.spacing[ this.axisOrder.indexOf('z') ];
 					IJKIndex = new THREE.Vector3( 0, RASIndex, 0 );
 					planeMatrix.multiply( new THREE.Matrix4().makeRotationX( - Math.PI / 2 ) );
 					positionOffset = ( volume.RASDimensions[ 1 ] - 1 ) / 2;
@@ -331,8 +338,8 @@
 					axisInIJK.set( 0, 0, 1 );
 					firstDirection.set( 1, 0, 0 );
 					secondDirection.set( 0, - 1, 0 );
-					firstSpacing = this.spacing[ 0 ];
-					secondSpacing = this.spacing[ 1 ];
+					firstSpacing = this.spacing[ this.axisOrder.indexOf('x') ];
+					secondSpacing = this.spacing[ this.axisOrder.indexOf('y') ];
 					IJKIndex = new THREE.Vector3( 0, 0, RASIndex );
 					positionOffset = ( volume.RASDimensions[ 2 ] - 1 ) / 2;
 					planeMatrix.setPosition( new THREE.Vector3( 0, 0, RASIndex - positionOffset ) );

examples/jsm/loaders/NRRDLoader.js

@@ -83,7 +83,7 @@ class NRRDLoader extends Loader {
 				case 'schar':
 					_array_type = Int8Array;
 					break;
-					// 2 byte data types
+				// 2 byte data types
 				case 'ushort':
 					_array_type = Uint16Array;
 					_chunkSize = 2;
@@ -92,7 +92,7 @@ class NRRDLoader extends Loader {
 					_array_type = Int16Array;
 					_chunkSize = 2;
 					break;
-					// 4 byte data types
+				// 4 byte data types
 				case 'uint':
 					_array_type = Uint32Array;
 					_chunkSize = 4;
@@ -363,31 +363,50 @@ class NRRDLoader extends Loader {
 		volume.xLength = volume.dimensions[ 0 ];
 		volume.yLength = volume.dimensions[ 1 ];
 		volume.zLength = volume.dimensions[ 2 ];
+
+		// Identify axis order in the space-directions matrix from the header if possible.
+		if (headerObject.vectors) {
+			const xIndex = headerObject.vectors.findIndex(vector => vector[0] !== 0);
+			const yIndex = headerObject.vectors.findIndex(vector => vector[1] !== 0);
+			const zIndex = headerObject.vectors.findIndex(vector => vector[2] !== 0);
+
+			const axisOrder = [];
+			axisOrder[xIndex] = 'x';
+			axisOrder[yIndex] = 'y';
+			axisOrder[zIndex] = 'z';
+			volume.axisOrder = axisOrder;
+		}
+		else {
+			volume.axisOrder = ['x', 'y', 'z'];
+		}
+
 		// spacing
-		const spacingX = ( new Vector3( headerObject.vectors[ 0 ][ 0 ], headerObject.vectors[ 0 ][ 1 ],
-			headerObject.vectors[ 0 ][ 2 ] ) ).length();
-		const spacingY = ( new Vector3( headerObject.vectors[ 1 ][ 0 ], headerObject.vectors[ 1 ][ 1 ],
-			headerObject.vectors[ 1 ][ 2 ] ) ).length();
-		const spacingZ = ( new Vector3( headerObject.vectors[ 2 ][ 0 ], headerObject.vectors[ 2 ][ 1 ],
-			headerObject.vectors[ 2 ][ 2 ] ) ).length();
+		const spacingX = new Vector3().fromArray( headerObject.vectors[ 0 ] ).length();
+		const spacingY = new Vector3().fromArray( headerObject.vectors[ 1 ] ).length();
+		const spacingZ = new Vector3().fromArray( headerObject.vectors[ 2 ] ).length();
 		volume.spacing = [ spacingX, spacingY, spacingZ ];
 
 
 		// Create IJKtoRAS matrix
 		volume.matrix = new Matrix4();
 
-		let _spaceX = 1;
-		let _spaceY = 1;
-		const _spaceZ = 1;
+		const transitionMatrix = new Matrix4();
 
-		if ( headerObject.space == 'left-posterior-superior' ) {
+		if ( headerObject.space === 'left-posterior-superior' ) {
 
-			_spaceX = - 1;
-			_spaceY = - 1;
+			transitionMatrix.set(
+				- 1, 0, 0, 0,
+				0, -1, 0, 0,
+				0, 0, 1, 0,
+				0, 0, 0, 1 );
 
 		} else if ( headerObject.space === 'left-anterior-superior' ) {
 
-			_spaceX = - 1;
+			transitionMatrix.set(
+				1, 0, 0, 0,
+				0, 1, 0, 0,
+				0, 0, -1, 0,
+				0, 0, 0, 1 );
 
 		}
 
@@ -395,20 +414,25 @@ class NRRDLoader extends Loader {
 		if ( ! headerObject.vectors ) {
 
 			volume.matrix.set(
-				_spaceX, 0, 0, 0,
-				0, _spaceY, 0, 0,
-				0, 0, _spaceZ, 0,
+				1, 0, 0, 0,
+				0, 1, 0, 0,
+				0, 0, 1, 0,
 				0, 0, 0, 1 );
 
 		} else {
 
 			const v = headerObject.vectors;
 
-			volume.matrix.set(
-				_spaceX * v[ 0 ][ 0 ], _spaceX * v[ 1 ][ 0 ], _spaceX * v[ 2 ][ 0 ], 0,
-				_spaceY * v[ 0 ][ 1 ], _spaceY * v[ 1 ][ 1 ], _spaceY * v[ 2 ][ 1 ], 0,
-				_spaceZ * v[ 0 ][ 2 ], _spaceZ * v[ 1 ][ 2 ], _spaceZ * v[ 2 ][ 2 ], 0,
-				0, 0, 0, 1 );
+			const ijk_to_transition = ( new Matrix4() ).set(
+				v[ 0 ][ 0 ], v[ 1 ][ 0 ], v[ 2 ][ 0 ], 0,
+				v[ 0 ][ 1 ], v[ 1 ][ 1 ], v[ 2 ][ 1 ], 0,
+				v[ 0 ][ 2 ], v[ 1 ][ 2 ], v[ 2 ][ 2 ], 0,
+				0, 0, 0, 1
+			)
+
+			const transition_to_ras = (new Matrix4()).multiplyMatrices( ijk_to_transition, transitionMatrix );
+
+			volume.matrix = transition_to_ras;
 
 		}
 

examples/jsm/misc/Volume.js

@@ -33,7 +33,10 @@ var Volume = function ( xLength, yLength, zLength, type, arrayBuffer ) {
 		 * @member {number} zLength Depth of the volume in the IJK coordinate system
 		 */
 		this.zLength = Number( zLength ) || 1;
-
+		/**
+		 * @member {Array<string>} The order of the Axis dictated by the NRRD header
+		 */
+		this.axisOrder = [ 'x', 'y', 'z' ];
 		/**
 		 * @member {TypedArray} data Data of the volume
 		 */
@@ -301,8 +304,8 @@ Volume.prototype = {
 				axisInIJK.set( 1, 0, 0 );
 				firstDirection.set( 0, 0, - 1 );
 				secondDirection.set( 0, - 1, 0 );
-				firstSpacing = this.spacing[ 2 ];
-				secondSpacing = this.spacing[ 1 ];
+				firstSpacing = this.spacing[ this.axisOrder.indexOf('z') ];
+				secondSpacing = this.spacing[ this.axisOrder.indexOf('y') ];
 				IJKIndex = new Vector3( RASIndex, 0, 0 );
 
 				planeMatrix.multiply( ( new Matrix4() ).makeRotationY( Math.PI / 2 ) );
@@ -313,8 +316,8 @@ Volume.prototype = {
 				axisInIJK.set( 0, 1, 0 );
 				firstDirection.set( 1, 0, 0 );
 				secondDirection.set( 0, 0, 1 );
-				firstSpacing = this.spacing[ 0 ];
-				secondSpacing = this.spacing[ 2 ];
+				firstSpacing = this.spacing[ this.axisOrder.indexOf('x') ];
+				secondSpacing = this.spacing[ this.axisOrder.indexOf('z') ];
 				IJKIndex = new Vector3( 0, RASIndex, 0 );
 
 				planeMatrix.multiply( ( new Matrix4() ).makeRotationX( - Math.PI / 2 ) );
@@ -326,8 +329,8 @@ Volume.prototype = {
 				axisInIJK.set( 0, 0, 1 );
 				firstDirection.set( 1, 0, 0 );
 				secondDirection.set( 0, - 1, 0 );
-				firstSpacing = this.spacing[ 0 ];
-				secondSpacing = this.spacing[ 1 ];
+				firstSpacing = this.spacing[ this.axisOrder.indexOf('x') ];
+				secondSpacing = this.spacing[ this.axisOrder.indexOf('y') ];
 				IJKIndex = new Vector3( 0, 0, RASIndex );
 
 				positionOffset = ( volume.RASDimensions[ 2 ] - 1 ) / 2;