Merge pull request #1 from sct-pipeline/orthogonal-plane-select-in-pr…

…ogress

First version of SCT orthogonal plane selection algorithm

README.md

@@ -1,2 +1,36 @@
 # venus-integration
 Code to process data for integrating acquisition planning with VENUS
+
+## How?
+### Step 1. Set up your directory structure
+```
+input/
+	2022-11-16-Scene.mrml
+	input-pointNormal-Plane-markup.json
+	input-anatomical-image.nii.gz (e.g. t2.nii.gz)
+output/
+preprocessing.sh
+slice_select.py
+write_slicer_markup_json.py
+```
+
+### Step 2. Preprocessing your data
+Label the spinal cord, vertebrae and vertebral boundaries within which you want to compute your slices. \
+Usage: `./preprocessing.sh anatomical_image.nii.gz contrast upper_vertebra lower_vertebra` \
+Labels (integer values) corresponding to each vertebra and disc can be found [here](https://spinalcordtoolbox.com/user_section/tutorials/registration-to-template/vertebral-labeling/labeling-conventions.html).
+```
+./preprocessing.sh t2.nii.gz t2 2 5 # 2 = mid C2; 5 = mid C5
+```
+
+### Step 3. Slice selection and orthogonal plane generation
+Find the indices of N slices (N = 5 in this example) that are equidistant along the centerline. \
+At each slice, compute a plane that is orthogonal to the centerline. 
+```
+python slice_select.py t2.nii.gz t2_seg.nii.gz t2_boundary.nii.gz t2 5
+```
+
+## Input
+* `input/t2.nii.gz` was downloaded from the [SCT t2 single subject tutorial](https://spinalcordtoolbox.com/user_section/tutorials/segmentation/before-starting.html).
+* `input/2022-11-16-Scene.mrml`: necessary to generate the planes as a markup file that can be read by slicer.
+* `input/input-pointNormal-Plane-markup.json`: necessary to generate the planes as a markup file that can be read by slicer.
+

fix_slicer_markup.py

@@ -0,0 +1,59 @@
+# Load Packages
+import sys
+import os
+import json
+
+class pointNormalPlane(object):
+    """
+    Class to represent point-Normal plane (that can be written to and from a .json file)
+
+    :param origin: list [x,y,z] with origin coordinates (float or int)
+    :param normal: list [a,b,c] with normal components (float or int)
+    :param orientation: 3-character string (e.g. 'RAS','LPS','LPI')
+    :param filename: string ending in .json
+    :param space: 'SCT_image' or 'anatomical'
+    """
+    def __init__(self,origin,normal,orientation,space='SCT_image'):
+        self.origin = origin
+        self.normal = normal
+        self.orientation = orientation
+        self.space = space
+
+    @classmethod
+    def fromJsonFile(cls,filename):
+        with open(filename, "r") as f:
+            data=f.read()
+        dict_tmp=json.loads(data)
+        return(pointNormalPlane(origin=dict_tmp['origin'],normal=dict_tmp['normal'],orientation=dict_tmp['orientation'],space=dict_tmp['space']))
+
+    def write_plane_json(self,filename):
+        new_plane_dict_json = json.dumps(self.__dict__)
+        if (filename.endswith('json')):
+            with open(filename, "w") as f:
+                f.write(new_plane_dict_json)
+        else:
+            with open(f'{filename}.json', "w") as f:
+                f.write(new_plane_dict_json)
+
+fname_in = sys.argv[1]
+plane=pointNormalPlane.fromJsonFile(sys.argv[1])
+print(plane.normal)
+print(plane.origin)
+slicer.util.loadScene('../../../input/2022-11-16-Scene.mrml')
+input_Node = getNode("input-pointNormal-Plane-markup")
+
+print(f'input plane normal: {input_Node.GetNormal()}')
+print(f'input plane normal: {input_Node.GetOrigin()}')
+print('\n')
+input_Node.SetNormal(plane.normal)
+input_Node.SetOrigin(plane.origin)
+print('New plane values:')
+print(f'new plane normal: {input_Node.GetNormal()}')
+print(f'new plane origin: {input_Node.GetOrigin()}')
+myStorageNode = input_Node.CreateDefaultStorageNode()
+
+filename=f'markup_{fname_in}'
+myStorageNode.SetFileName(filename)
+myStorageNode.WriteData(input_Node)
+myStorageNode.UnRegister(None)
+sys.exit(0)

fix_slicer_markup_memory_leak.sh

@@ -0,0 +1,14 @@
+#!/bin/bash
+
+for dir in output/*/test*slices
+do
+	cd $dir
+	for file in plane*RAS_*json
+	do 
+		if [ ! -f markup_$(basename $file) ]; then
+			echo "markup_$(basename $file) does not exist!"
+			/Applications/Slicer.app/Contents/MacOS/Slicer --no-splash --no-main-window --python-script ../../../fix_slicer_markup.py $file
+		fi
+	done
+	cd ../../..
+done

input/2022-11-16-Scene.mrml

@@ -0,0 +1,42 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<MRML version="Slicer 5.1.0 31234" userTags="">
+<Crosshair
+  id="vtkMRMLCrosshairNodedefault" name="Crosshair" hideFromEditors="true" selectable="true" selected="false" singletonTag="default" crosshairMode="NoCrosshair" crosshairBehavior="OffsetJumpSlice" crosshairThickness="Fine" crosshairRAS="0 0 0"></Crosshair>
+ <Selection
+  id="vtkMRMLSelectionNodeSingleton" name="Selection" hideFromEditors="true" selectable="true" selected="false" singletonTag="Singleton" references="ActivePlaceNode:vtkMRMLMarkupsPlaneNode1;unit/area:vtkMRMLUnitNodeApplicationArea;unit/frequency:vtkMRMLUnitNodeApplicationFrequency;unit/intensity:vtkMRMLUnitNodeApplicationIntensity;unit/length:vtkMRMLUnitNodeApplicationLength;unit/time:vtkMRMLUnitNodeApplicationTime;unit/velocity:vtkMRMLUnitNodeApplicationVelocity;unit/volume:vtkMRMLUnitNodeApplicationVolume;" activePlaceNodeClassName="vtkMRMLMarkupsPlaneNode" ></Selection>
+ <Interaction
+  id="vtkMRMLInteractionNodeSingleton" name="Interaction" hideFromEditors="true" selectable="true" selected="false" singletonTag="Singleton" currentInteractionMode="ViewTransform" placeModePersistence="false" lastInteractionMode="ViewTransform" ></Interaction>
+ <View
+  id="vtkMRMLViewNode1" name="View1" hideFromEditors="false" selectable="true" selected="false" singletonTag="1" attributes="MappedInLayout:1" layoutLabel="1" layoutName="1" active="false" visibility="true" backgroundColor="0.756863 0.764706 0.909804" backgroundColor2="0.454902 0.470588 0.745098" layoutColor="0.454902 0.513725 0.913725" orientationMarkerType="none" orientationMarkerSize="medium" rulerType="none" rulerColor="white" AxisLabels="L;R;P;A;I;S" fieldOfView="200" letterSize="0.05" boxVisible="true" boxColor="1 0 1" fiducialsVisible="true" fiducialLabelsVisible="true" axisLabelsVisible="true" axisLabelsCameraDependent="true" animationMode="Off" viewAxisMode="LookFrom" spinDegrees="2" spinMs="5" spinDirection="YawLeft" rotateDegrees="5" rockLength="200" rockCount="0" stereoType="NoStereo" renderMode="Perspective" useDepthPeeling="1" gpuMemorySize="0" autoReleaseGraphicsResources="false" expectedFPS="8" volumeRenderingQuality="Normal" raycastTechnique="Composite" volumeRenderingSurfaceSmoothing="0" volumeRenderingOversamplingFactor="2" linkedControl="0" ></View>
+ <Slice
+  id="vtkMRMLSliceNodeRed" name="Red" hideFromEditors="false" selectable="true" selected="false" singletonTag="Red" attributes="MappedInLayout:1" layoutLabel="R" layoutName="Red" active="false" visibility="true" backgroundColor="0 0 0" backgroundColor2="0 0 0" layoutColor="0.952941 0.290196 0.2" orientationMarkerType="none" orientationMarkerSize="medium" rulerType="none" rulerColor="white" AxisLabels="L;R;P;A;I;S" fieldOfView="299.842 182.353 1" dimensions="490 298 1" xyzOrigin="0 0 0" sliceResolutionMode="1" uvwExtents="299.842 182.353 1" uvwDimensions="256 256 1" uvwOrigin="0 0 0" activeSlice="0" layoutGridRows="1" layoutGridColumns="1" sliceToRAS="-1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1" orientationMatrixAxial="-1 0 0 0 1 0 0 0 1" orientationMatrixSagittal="0 0 -1 -1 0 0 0 1 0" orientationMatrixCoronal="-1 0 0 0 0 1 0 1 0" orientation="Axial" defaultOrientation="Axial" orientationReference="Axial" jumpMode="1" sliceVisibility="false" widgetVisibility="false" widgetOutlineVisibility="true" useLabelOutline="false" sliceSpacingMode="0" prescribedSliceSpacing="1 1 1" ></Slice>
+ <Slice
+  id="vtkMRMLSliceNodeGreen" name="Green" hideFromEditors="false" selectable="true" selected="false" singletonTag="Green" attributes="MappedInLayout:1" layoutLabel="G" layoutName="Green" active="false" visibility="true" backgroundColor="0 0 0" backgroundColor2="0 0 0" layoutColor="0.431373 0.690196 0.294118" orientationMarkerType="none" orientationMarkerSize="medium" rulerType="none" rulerColor="white" AxisLabels="L;R;P;A;I;S" fieldOfView="301.868 182.353 1" dimensions="490 296 1" xyzOrigin="0 0 0" sliceResolutionMode="1" uvwExtents="301.868 182.353 1" uvwDimensions="256 256 1" uvwOrigin="0 0 0" activeSlice="0" layoutGridRows="1" layoutGridColumns="1" sliceToRAS="-1 0 0 0 0 0 1 0 0 1 0 0 0 0 0 1" orientationMatrixAxial="-1 0 0 0 1 0 0 0 1" orientationMatrixSagittal="0 0 -1 -1 0 0 0 1 0" orientationMatrixCoronal="-1 0 0 0 0 1 0 1 0" orientation="Coronal" defaultOrientation="Coronal" orientationReference="Coronal" jumpMode="1" sliceVisibility="false" widgetVisibility="false" widgetOutlineVisibility="true" useLabelOutline="false" sliceSpacingMode="0" prescribedSliceSpacing="1 1 1" ></Slice>
+ <Slice
+  id="vtkMRMLSliceNodeYellow" name="Yellow" hideFromEditors="false" selectable="true" selected="false" singletonTag="Yellow" attributes="MappedInLayout:1" layoutLabel="Y" layoutName="Yellow" active="false" visibility="true" backgroundColor="0 0 0" backgroundColor2="0 0 0" layoutColor="0.929412 0.835294 0.298039" orientationMarkerType="none" orientationMarkerSize="medium" rulerType="none" rulerColor="white" AxisLabels="L;R;P;A;I;S" fieldOfView="300.103 181.287 1" dimensions="490 296 1" xyzOrigin="0 0 0" sliceResolutionMode="1" uvwExtents="300.103 181.287 1" uvwDimensions="256 256 1" uvwOrigin="0 0 0" activeSlice="0" layoutGridRows="1" layoutGridColumns="1" sliceToRAS="0 0 -1 0 -1 0 0 0 0 1 0 0 0 0 0 1" orientationMatrixAxial="-1 0 0 0 1 0 0 0 1" orientationMatrixSagittal="0 0 -1 -1 0 0 0 1 0" orientationMatrixCoronal="-1 0 0 0 0 1 0 1 0" orientation="Sagittal" defaultOrientation="Sagittal" orientationReference="Sagittal" jumpMode="1" sliceVisibility="false" widgetVisibility="false" widgetOutlineVisibility="true" useLabelOutline="false" sliceSpacingMode="0" prescribedSliceSpacing="1 1 1" ></Slice>
+ <Layout
+  id="vtkMRMLLayoutNodevtkMRMLLayoutNode" name="Layout" hideFromEditors="true" selectable="true" selected="false" singletonTag="vtkMRMLLayoutNode" currentViewArrangement="3" guiPanelVisibility="1" bottomPanelVisibility ="1" guiPanelLR="0" collapseSliceControllers="0"
+ numberOfCompareViewRows="1" numberOfCompareViewColumns="1" numberOfLightboxRows="6" numberOfLightboxColumns="6" mainPanelSize="400" secondaryPanelSize="400" ></Layout>
+ <SliceComposite
+  id="vtkMRMLSliceCompositeNodeRed" name="SliceComposite" hideFromEditors="true" selectable="true" selected="false" singletonTag="Red" compositing="0" foregroundOpacity="0" labelOpacity="1" linkedControl="0" hotLinkedControl="0" fiducialVisibility="1" fiducialLabelVisibility="1" layoutName="Red" annotationSpace="IJKAndRAS" annotationMode="All" doPropagateVolumeSelection="1" ></SliceComposite>
+ <SliceComposite
+  id="vtkMRMLSliceCompositeNodeGreen" name="SliceComposite_1" hideFromEditors="true" selectable="true" selected="false" singletonTag="Green" compositing="0" foregroundOpacity="0" labelOpacity="1" linkedControl="0" hotLinkedControl="0" fiducialVisibility="1" fiducialLabelVisibility="1" layoutName="Green" annotationSpace="IJKAndRAS" annotationMode="All" doPropagateVolumeSelection="1" ></SliceComposite>
+ <SubjectHierarchy
+  id="vtkMRMLSubjectHierarchyNode1" name="SubjectHierarchy" hideFromEditors="false" selectable="true" selected="false" attributes="SubjectHierarchyVersion:2" >
+   <SubjectHierarchyItem id="3" name="Scene" parent="0" type="" expanded="true" attributes="Level^Scene|">
+     <SubjectHierarchyItem id="7" dataNode="vtkMRMLMarkupsPlaneNode1" parent="3" type="Markups" expanded="true"></SubjectHierarchyItem></SubjectHierarchyItem></SubjectHierarchy>
+ <SliceComposite
+  id="vtkMRMLSliceCompositeNodeYellow" name="SliceComposite_2" hideFromEditors="true" selectable="true" selected="false" singletonTag="Yellow" compositing="0" foregroundOpacity="0" labelOpacity="1" linkedControl="0" hotLinkedControl="0" fiducialVisibility="1" fiducialLabelVisibility="1" layoutName="Yellow" annotationSpace="IJKAndRAS" annotationMode="All" doPropagateVolumeSelection="1" ></SliceComposite>
+ <Camera
+  id="vtkMRMLCameraNode1" name="Camera" description="Default Scene Camera" hideFromEditors="false" selectable="true" selected="false" singletonTag="1" userTags="" position="0 500 0" focalPoint="0 0 0" viewUp="0 0 1" parallelProjection="false" parallelScale="1" viewAngle="30" appliedTransform="1 0 0 0  0 1 0 0  0 0 1 0  0 0 0 1" ></Camera>
+ <ClipModels
+  id="vtkMRMLClipModelsNodevtkMRMLClipModelsNode" name="ClipModels" hideFromEditors="true" selectable="true" selected="false" singletonTag="vtkMRMLClipModelsNode" clipType="0" redSliceClipState="0" yellowSliceClipState="0" greenSliceClipState="0" ></ClipModels>
+ <ScriptedModule
+  id="vtkMRMLScriptedModuleNodeDataProbe" name="ScriptedModule" hideFromEditors="true" selectable="true" selected="false" singletonTag="DataProbe" ModuleName ="DataProbe" ></ScriptedModule>
+ <MarkupsPlaneJsonStorage
+  id="vtkMRMLMarkupsPlaneJsonStorageNode1" name="MarkupsPlaneJsonStorage" hideFromEditors="true" selectable="true" selected="false" fileName="input-pointNormal-Plane-markup.json" useCompression="1" defaultWriteFileExtension="mrk.json" readState="0" writeState="0" coordinateSystem="LPS" ></MarkupsPlaneJsonStorage>
+ <MarkupsPlane
+  id="vtkMRMLMarkupsPlaneNode1" name="input-pointNormal-Plane-markup" description="" hideFromEditors="false" selectable="true" selected="false" references="display:vtkMRMLMarkupsPlaneDisplayNode1;storage:vtkMRMLMarkupsPlaneJsonStorageNode1;" userTags="" locked="false" controlPointLabelFormat="%N-%d" interactionHandleToWorldMatrix="1 0 0 4.19558  0 1 0 -18.3518  0 0 1 -17.5416  0 0 0 1" maximumNumberOfControlPoints="1" requiredNumberOfControlPoints="1" ></MarkupsPlane>
+ <MarkupsPlaneDisplay
+  id="vtkMRMLMarkupsPlaneDisplayNode1" name="MarkupsPlaneDisplay" hideFromEditors="true" selectable="true" selected="false" color="0.4 1 0" edgeColor="0 0 0" selectedColor="1 0.500008 0.500008" selectedAmbient="0.4" ambient="0" diffuse="1" selectedSpecular="0.5" specular="0" power="1" metallic="0" roughness="0.5" opacity="1" sliceIntersectionOpacity="1" pointSize="1" lineWidth="1" representation="2" lighting="true" interpolation="1" shading="true" visibility="true" visibility2D="true" visibility3D="true" edgeVisibility="false" clipping="false" sliceIntersectionThickness="1" frontfaceCulling="false" backfaceCulling="false" scalarVisibility="false" vectorVisibility="false" tensorVisibility="false" interpolateTexture="false" scalarRangeFlag="UseData" scalarRange="0 100" activeAttributeLocation="point" viewNodeRef="" folderDisplayOverrideAllowed="true" propertiesLabelVisibility="true" pointLabelsVisibility="false" textScale="3" glyphScale="3" glyphSize="5" useGlyphScale="true" glyphType="Sphere3D" snapMode="toVisibleSurface" sliceProjection="false" sliceProjectionUseFiducialColor="true" sliceProjectionOutlinedBehindSlicePlane="false" sliceProjectionColor="1 1 1" sliceProjectionOpacity="0.6" curveLineSizeMode="UseLineThickness" lineThickness="0.2" lineDiameter="1" lineColorFadingStart="1" lineColorFadingEnd="10" lineColorFadingSaturation="1" lineColorFadingHueOffset="0" handlesInteractive="true" translationHandleVisibility="false" rotationHandleVisibility="false" scaleHandleVisibility="true" interactionHandleScale="3" fillVisibility="true" outlineVisibility="true" fillOpacity="0.5" outlineOpacity="1" occludedVisibility="false" occludedOpacity="0.3" textProperty="font-family:Arial;font-size:5px;font-style:normal;font-weight:bold;color:rgba(255,255,255,1);background-color:rgba(0,0,0,0);border-width:1px;border-color:rgba(255,255,255,0.0);text-shadow:1px -1px 2px rgba(0,0,0,1.0);" activeColor="0.4 1 0" ></MarkupsPlaneDisplay>
+</MRML>

input/input-pointNormal-Plane-markup.json

@@ -0,0 +1,38 @@
+{
+    "@schema": "https://raw.githubusercontent.com/slicer/slicer/master/Modules/Loadable/Markups/Resources/Schema/markups-schema-v1.0.3.json#",
+    "markups": [
+        {
+            "type": "Plane",
+            "coordinateSystem": "LPS",
+            "coordinateUnits": "mm",
+            "locked": false,
+            "fixedNumberOfControlPoints": false,
+            "labelFormat": "%N-%d",
+            "lastUsedControlPointNumber": 1,
+            "planeType": "pointNormal",
+            "sizeMode": "auto",
+            "autoScalingFactor": 1.0,
+            "center": [-4.195577, 18.3518, -17.54164],
+            "normal": [-0.0, -0.0, 1.0],
+            "objectToBase": [-1.0, -0.0, -0.0, -0.0, -0.0, -1.0, -0.0, -0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0],
+            "baseToNode": [-1.0, -0.0, -0.0, -4.195577, -0.0, -1.0, -0.0, 18.3518, 0.0, 0.0, 1.0, -17.54164, 0.0, 0.0, 0.0, 1.0],
+            "orientation": [-1.0, -0.0, -0.0, -0.0, -1.0, -0.0, 0.0, 0.0, 1.0],
+            "size": [100.0, 100.0, 0.0],
+            "planeBounds": [-50.0, 50.0, -50.0, 50.0],
+            "controlPoints": [
+                {
+                    "id": "1",
+                    "label": "Plane_RAS-1",
+                    "description": "",
+                    "associatedNodeID": "",
+                    "position": [-4.195577, 18.3518, -17.54164],
+                    "orientation": [-1.0, -0.0, -0.0, -0.0, -1.0, -0.0, 0.0, 0.0, 1.0],
+                    "selected": true,
+                    "locked": false,
+                    "visibility": true,
+                    "positionStatus": "defined"
+                }
+            ]
+        }
+    ]
+}

preprocessing.sh

@@ -0,0 +1,23 @@
+#!/bin/bash
+
+if [[ $# -ne 4 ]]; then
+	echo 'Usage: ./sct_plane_select.sh img.nii.gz contrast upper_vert lower_vert'
+    exit
+fi
+
+input_img=$1
+contrast=$2
+upper_vert=$3
+lower_vert=$4
+
+cd output
+
+# Label SC (w/ deep learning, instead of sct_prop_seg)
+sct_deepseg_sc -i ../input/$input_img -c $contrast -o $(basename $input_img .nii.gz)_seg.nii.gz
+
+# label vertebrae & disks
+sct_label_vertebrae -i ../input/$input_img -s $(basename $input_img .nii.gz)_seg.nii.gz -c $contrast
+
+# select area we want to work with (ex: from mid-vertebrum C2 to mid-vertebrum C5)
+sct_label_utils -i $(basename $input_img .nii.gz)_seg_labeled.nii.gz -vert-body $upper_vert,$lower_vert -o $(basename $input_img .nii.gz)_boundary.nii.gz
+