From a7e86dd3fcc3535f616554d1e5d09e4f4246de94 Mon Sep 17 00:00:00 2001 From: Mikhail Polkovnikov Date: Mon, 6 Nov 2023 15:57:10 +0300 Subject: [PATCH] ENH: Initial support of orientated CT volumes If IJKToRASDirection matrix _A_ isn't LPS transform, | a11 a12 a13 a14 | | -1 0 0 0 | _A_ = | a21 a22 a23 a24 | != | 0 -1 0 0 | | a31 a32 a33 a34 | | 0 0 1 0 | | a41 a42 a43 a44 | | 0 0 0 1 | then, the volume direction to LPS transform matrix _R_ can be applied to the CT volume to make IJKToRASDirection matrix of the CT volume to be LPS orientated after hardening the transform. the matrix transformation: | -1 0 0 0 | _R_ = _A_^-1 * | 0 -1 0 0 | | 0 0 1 0 | | 0 0 0 1 | the matrix widget shows elements of the transformation matrix _R_, if the matrix _A_ is LPS orientated => the matrix _R_ is an identity matrix. --- Docs/user_guide/modules/drr.md | 74 ++++-- .../vtkSlicerDrrImageComputationLogic.cxx | 249 +++++++++++++++++- .../Logic/vtkSlicerDrrImageComputationLogic.h | 14 +- .../MRML/vtkMRMLDrrImageComputationNode.cxx | 42 ++- .../MRML/vtkMRMLDrrImageComputationNode.h | 2 + .../qSlicerDrrImageComputationModuleWidget.ui | 138 +++++++--- ...qSlicerDrrImageComputationModuleWidget.cxx | 105 +++++++- .../qSlicerDrrImageComputationModuleWidget.h | 1 + PlmDrr/plastimatch_slicer_drr.cxx | 88 ++++++- 9 files changed, 635 insertions(+), 78 deletions(-) diff --git a/Docs/user_guide/modules/drr.md b/Docs/user_guide/modules/drr.md index b8058222..e5ef28df 100644 --- a/Docs/user_guide/modules/drr.md +++ b/Docs/user_guide/modules/drr.md @@ -50,24 +50,56 @@ based on [plastimatch drr](http://www.plastimatch.org/drr.html) ones. ### Graphical User Interface loadable module (GUI) -![image](https://user-images.githubusercontent.com/3785912/118298680-726b9e80-b4e8-11eb-8a39-c5e607459e62.png) +![image](https://private-user-images.githubusercontent.com/3785912/294439297-c315f0db-9ab7-4a91-9eaa-10a4698d7cd8.png?jwt=eyJhbGciOiJIUzI1NiIsInR5cCI6IkpXVCJ9.eyJpc3MiOiJnaXRodWIuY29tIiwiYXVkIjoicmF3LmdpdGh1YnVzZXJjb250ZW50LmNvbSIsImtleSI6ImtleTUiLCJleHAiOjE3MDQ0NDQ3NzMsIm5iZiI6MTcwNDQ0NDQ3MywicGF0aCI6Ii8zNzg1OTEyLzI5NDQzOTI5Ny1jMzE1ZjBkYi05YWI3LTRhOTEtOWVhYS0xMGE0Njk4ZDdjZDgucG5nP1gtQW16LUFsZ29yaXRobT1BV1M0LUhNQUMtU0hBMjU2JlgtQW16LUNyZWRlbnRpYWw9QUtJQVZDT0RZTFNBNTNQUUs0WkElMkYyMDI0MDEwNSUyRnVzLWVhc3QtMSUyRnMzJTJGYXdzNF9yZXF1ZXN0JlgtQW16LURhdGU9MjAyNDAxMDVUMDg0NzUzWiZYLUFtei1FeHBpcmVzPTMwMCZYLUFtei1TaWduYXR1cmU9MmYyNGIzYmNmYjM2ZDI0YWU0NGFiM2RiYjg0ZGE4ZDQ3NGIyMjFiYWNkMmQ2Y2Y4MzlmY2I2OWI0MGUzZDI3YyZYLUFtei1TaWduZWRIZWFkZXJzPWhvc3QmYWN0b3JfaWQ9MCZrZXlfaWQ9MCZyZXBvX2lkPTAifQ.J_kchHVITj00Qsd6Z84Wujp8kEXvknmQYZ9SoiWgqbA) Loadable GUI module "DRR Image Computation" uses CLI module's logic and nodes data to calculate and visualize DRR image. It also shows basic detector elements such as: detector boundary, detector normal vector, detector view-up vector, detector image origin (0,0) pixel, image subwindow boundary as markups data on a slice and 3D views. -Markups data is only for WYSIWYG purpose. +This module _partially_ supports calculation of DRR for orientated CT volumes, +if IJKToRASDirection matrix of the volume isn't as below. + +| -1 0 0 0 | +| 0 -1 0 0 | +| 0 0 1 0 | +| 0 0 0 1 | + +The markups data is only for WYSIWYG purpose. #### Reference input nodes -![image](https://user-images.githubusercontent.com/3785912/118298709-7e576080-b4e8-11eb-8c2b-b1a4f8222eba.png) +![image](https://private-user-images.githubusercontent.com/3785912/294439383-a0718f87-81bf-4a9f-880b-34876a6152ed.png?jwt=eyJhbGciOiJIUzI1NiIsInR5cCI6IkpXVCJ9.eyJpc3MiOiJnaXRodWIuY29tIiwiYXVkIjoicmF3LmdpdGh1YnVzZXJjb250ZW50LmNvbSIsImtleSI6ImtleTUiLCJleHAiOjE3MDQ0NDQ3NzMsIm5iZiI6MTcwNDQ0NDQ3MywicGF0aCI6Ii8zNzg1OTEyLzI5NDQzOTM4My1hMDcxOGY4Ny04MWJmLTRhOWYtODgwYi0zNDg3NmE2MTUyZWQucG5nP1gtQW16LUFsZ29yaXRobT1BV1M0LUhNQUMtU0hBMjU2JlgtQW16LUNyZWRlbnRpYWw9QUtJQVZDT0RZTFNBNTNQUUs0WkElMkYyMDI0MDEwNSUyRnVzLWVhc3QtMSUyRnMzJTJGYXdzNF9yZXF1ZXN0JlgtQW16LURhdGU9MjAyNDAxMDVUMDg0NzUzWiZYLUFtei1FeHBpcmVzPTMwMCZYLUFtei1TaWduYXR1cmU9NzNlYWI1MDk2M2Q0MzQ3NWQ5MzZhMWNiMzlmNmVmNmQ0NzFkNDQ5Y2U1YjBmNDE5NzZlNTYzM2E1NGEyYjY3NSZYLUFtei1TaWduZWRIZWFkZXJzPWhvc3QmYWN0b3JfaWQ9MCZrZXlfaWQ9MCZyZXBvX2lkPTAifQ.ZrUVivDARk6G1kLMgdk36E3r-z0UzOV7SLK2OUBZzbg) 1. **CT volume**: Input CT data 2. **RT beam**: Input RT beam (vtkMRMLRTBeamNode) for source and detector orientation parameters 2. **Camera**: Camera node (vtkMRMLCameraNode) to update if needed beam geometry and transformation 3. **Update beam**: Update beam geometry and transformation using camera node data -4. **Show DRR markups**: Show or hide detector markups +4. **Volume direction to LPS transform matrix**: The matrix widget shows elements of the transformation matrix R + The volume direction to LPS transform matrix _R_ + can be applied to the volume, to make IJKToRASDirectionMatrix of the volume a LPS orientated + after hardening the transform. + + If IJKToRASDirection matrix _A_ isn't LPS transform + + | a11 a12 a13 a14 | | -1 0 0 0 | + _A_ = | a21 a22 a23 a24 | != | 0 -1 0 0 | + | a31 a32 a33 a34 | | 0 0 1 0 | + | a41 a42 a43 a44 | | 0 0 0 1 | + + the volume direction to LPS transform matrix _R_ can be applied to the volume, + by clicking on `Apply transform to volume` button, to make + the IJKToRASDirection matrix a LPS orientated after hardening the transform. + the matrix transformation: + + | -1 0 0 0 | + _R_ = _A_^-1 * | 0 -1 0 0 | + | 0 0 1 0 | + | 0 0 0 1 | + + if matrix _A_ is LPS orientated => matrix _R_ is an identity matrix. + +5. **Show DRR markups**: Show or hide detector markups CT data is represented by vtkMRMLScalarVolumeNode object. RT beam data is represented by vtkMRMLRTBeamNode object. Camera data is represented by vtkMRMLCameraNode object. @@ -76,27 +108,33 @@ represented by vtkMRMLRTBeamNode object. Camera data is represented by vtkMRMLCa ![image](https://user-images.githubusercontent.com/3785912/107616499-40d00680-6c5f-11eb-9d45-c6cccc9e12cd.png) -4. **Isocenter to imager distance**: Distance from isocenter to detector center in mm -5. **Imager resolution (columns, rows)**: Detector resolution in pixels -6. **Imager spacing in mm (columns, rows)**: Detector spacing in mm -7. **Image window parameters**: Use and setup image subwindow or a whole detector +6. **Isocenter to imager distance**: Distance from isocenter to detector center in mm +7. **Imager resolution (columns, rows)**: Detector resolution in pixels +8. **Imager spacing in mm (columns, rows)**: Detector spacing in mm +9. **Image window parameters**: Use and setup image subwindow or a whole detector #### Image Window Parameters -8. **Columns**: Number of image columns in subwindow -9. **Rows**: Number of image rows in subwindow +10. **Columns**: Number of image columns in subwindow +11. **Rows**: Number of image rows in subwindow + +#### Markups perspective projection on the imager plane + +![image](https://private-user-images.githubusercontent.com/3785912/294439414-d1bbb6ac-4838-47a7-88fc-9ec76bfe19b2.png?jwt=eyJhbGciOiJIUzI1NiIsInR5cCI6IkpXVCJ9.eyJpc3MiOiJnaXRodWIuY29tIiwiYXVkIjoicmF3LmdpdGh1YnVzZXJjb250ZW50LmNvbSIsImtleSI6ImtleTUiLCJleHAiOjE3MDQ0NDQ3NzMsIm5iZiI6MTcwNDQ0NDQ3MywicGF0aCI6Ii8zNzg1OTEyLzI5NDQzOTQxNC1kMWJiYjZhYy00ODM4LTQ3YTctODhmYy05ZWM3NmJmZTE5YjIucG5nP1gtQW16LUFsZ29yaXRobT1BV1M0LUhNQUMtU0hBMjU2JlgtQW16LUNyZWRlbnRpYWw9QUtJQVZDT0RZTFNBNTNQUUs0WkElMkYyMDI0MDEwNSUyRnVzLWVhc3QtMSUyRnMzJTJGYXdzNF9yZXF1ZXN0JlgtQW16LURhdGU9MjAyNDAxMDVUMDg0NzUzWiZYLUFtei1FeHBpcmVzPTMwMCZYLUFtei1TaWduYXR1cmU9M2I4NzQ3ZWU0ZWFlOTExZTllZTM1YTMzZjE0OTI2N2VkOWQ0MWEwOTdmYzQ0MGFhMTU0MWFiY2Y0NjdmYTRkNSZYLUFtei1TaWduZWRIZWFkZXJzPWhvc3QmYWN0b3JfaWQ9MCZrZXlfaWQ9MCZyZXBvX2lkPTAifQ.yUYlwYdj9JdqD7CUyS-zFfZ4tV0Yezux__UI35jEUlg) + +Calculates projection of fiducial markups points in imager plane. #### Plastimatch DRR image processing ![image](https://user-images.githubusercontent.com/3785912/107617306-b38db180-6c60-11eb-9dd1-b2751b23a314.png) -10. **Use exponential mapping**: Apply exponential mapping of output values -11. **Autoscale**: Automatically rescale intensity -12. **Invert**: Invert image intensity -13. **Range**: Range intensity in form (min, max) -14. **Reconstruction algorithm**: Type of reconstruction algorithm (Type of exposure algorithm in CLI module) -15. **Hounsfield units conversion**: Type of Hounsfield Units conversion during computation -16. **Threading**: Type of parallelism of computation -17. **Hounsfield units threshold**: HU values below threshold are counted as -1000 (Air value) +12. **Use exponential mapping**: Apply exponential mapping of output values +13. **Autoscale**: Automatically rescale intensity +14. **Invert**: Invert image intensity +15. **Range**: Range intensity in form (min, max) +16. **Reconstruction algorithm**: Type of reconstruction algorithm (Type of exposure algorithm in CLI module) +17. **Hounsfield units conversion**: Type of Hounsfield Units conversion during computation +18. **Threading**: Type of parallelism of computation +19. **Hounsfield units threshold**: HU values below threshold are counted as -1000 (Air value) #### Plastimatch DRR command arguments (read only) diff --git a/DrrImageComputation/Logic/vtkSlicerDrrImageComputationLogic.cxx b/DrrImageComputation/Logic/vtkSlicerDrrImageComputationLogic.cxx index 82f85bd8..31c6229d 100644 --- a/DrrImageComputation/Logic/vtkSlicerDrrImageComputationLogic.cxx +++ b/DrrImageComputation/Logic/vtkSlicerDrrImageComputationLogic.cxx @@ -79,6 +79,38 @@ const char* vtkSlicerDrrImageComputationLogic::NORMAL_VECTOR_MARKUPS_NODE_NAME = const char* vtkSlicerDrrImageComputationLogic::VUP_VECTOR_MARKUPS_NODE_NAME = "VupVector"; // line const char* vtkSlicerDrrImageComputationLogic::RTIMAGE_TRANSFORM_NODE_NAME = "DrrImageComputationTransform"; +const char* vtkSlicerDrrImageComputationLogic::VOLUME_TO_LPS_TRANSFORM_NODE_NAME = "DrrVolumeToLpsTransform"; + +namespace { + +bool IsMatrixTransformsRasToLps(vtkMatrix4x4* mat) +{ + if (!mat) + { + return false; + } +// constexpr double local_epsilon = std::numeric_limits< double >::epsilon(); + constexpr double local_epsilon = EPSILON; + bool result = bool(fabs(mat->GetElement(0, 0) + 1.) <= local_epsilon); + result &= bool(fabs(mat->GetElement(1, 1) + 1.) <= local_epsilon); + result &= bool(fabs(mat->GetElement(2, 2) - 1.) <= local_epsilon); + result &= bool(fabs(mat->GetElement(3, 3) - 1.) <= local_epsilon); + result &= bool(fabs(mat->GetElement(1, 0)) <= local_epsilon); + result &= bool(fabs(mat->GetElement(2, 0)) <= local_epsilon); + result &= bool(fabs(mat->GetElement(2, 1)) <= local_epsilon); + result &= bool(fabs(mat->GetElement(3, 0)) <= local_epsilon); + result &= bool(fabs(mat->GetElement(3, 1)) <= local_epsilon); + result &= bool(fabs(mat->GetElement(3, 2)) <= local_epsilon); + result &= bool(fabs(mat->GetElement(0, 1)) <= local_epsilon); + result &= bool(fabs(mat->GetElement(0, 2)) <= local_epsilon); + result &= bool(fabs(mat->GetElement(1, 2)) <= local_epsilon); + result &= bool(fabs(mat->GetElement(0, 3)) <= local_epsilon); + result &= bool(fabs(mat->GetElement(1, 3)) <= local_epsilon); + result &= bool(fabs(mat->GetElement(2, 3)) <= local_epsilon); + return result; +} + +} //---------------------------------------------------------------------------- vtkStandardNewMacro(vtkSlicerDrrImageComputationLogic); @@ -1004,7 +1036,7 @@ vtkMRMLScalarVolumeNode* vtkSlicerDrrImageComputationLogic::ComputePlastimatchDR std::stringstream isocenterStream; double isocenter[3] = {}; - parameterNode->GetIsocenterPositionLPS(isocenter); + parameterNode->GetIsocenterPositionRAS(isocenter); isocenterStream << isocenter[0] << "," << isocenter[1] << "," << isocenter[2]; cmdNode->SetParameterAsString( "isocenterPosition", isocenterStream.str()); @@ -1234,8 +1266,18 @@ bool vtkSlicerDrrImageComputationLogic::SetupGeometry( vtkMRMLDrrImageComputatio vtkNew rtImageCenterToCornerTransform; rtImageCenterToCornerTransform->Identity(); - rtImageCenterToCornerTransform->Translate( -1. * rtImagePosition[0], 0.0, rtImagePosition[1]); + vtkMRMLScalarVolumeNode* ctVolumeNode = parameterNode->GetReferenceVolumeNode(); + vtkNew< vtkMatrix4x4 > ijkToRasDirectionMatrix; + ctVolumeNode->GetIJKToRASDirectionMatrix(ijkToRasDirectionMatrix); + if (!IsMatrixTransformsRasToLps(ijkToRasDirectionMatrix)) + { + rtImageCenterToCornerTransform->Translate( rtImagePosition[0], 0.0, rtImagePosition[1]); + } + else + { + rtImageCenterToCornerTransform->Translate( -1. * rtImagePosition[0], 0.0, rtImagePosition[1]); + } // Create isocenter to RAS transform // The transformation below is based section C.8.8 in DICOM standard volume 3: // "Note: IEC document 62C/269/CDV 'Amendment to IEC 61217: Radiotherapy Equipment - @@ -1659,6 +1701,8 @@ vtkMRMLLinearTransformNode* vtkSlicerDrrImageComputationLogic::UpdateImageTransf // Update transforms in IEC logic from beam node parameters iecLogic->UpdateIECTransformsFromBeam(beamNode); + // (a BUG?) For RT Image correct orientation PatientSupport -> Fixed Reference MUST have a negative sign + iecLogic->UpdatePatientSupportRotationToFixedReferenceTransform(-1. * beamNode->GetCouchAngle()); // Dynamic transform from Gantry to RAS // Transformation path: @@ -1701,6 +1745,8 @@ bool vtkSlicerDrrImageComputationLogic::GetRtImageTransformMatrixFromBeam(vtkMRM // Update transforms in IEC logic from beam node parameters iecLogic->UpdateIECTransformsFromBeam(beamNode); + // (a BUG?) For RT Image correct orientation PatientSupport -> Fixed Reference MUST have a negative sign + iecLogic->UpdatePatientSupportRotationToFixedReferenceTransform(-1. * beamNode->GetCouchAngle()); // Dynamic transform from Gantry to RAS // Transformation path: @@ -1747,12 +1793,30 @@ void vtkSlicerDrrImageComputationLogic::UpdateNormalAndVupVectors(vtkMRMLDrrImag return; } + vtkMRMLLinearTransformNode* volumeToLpsTransformNode = this->UpdateVolumeToLPSTransform(parameterNode); + if (volumeToLpsTransformNode) + { + vtkTransform* volumeToLpsTransform = vtkTransform::SafeDownCast(volumeToLpsTransformNode->GetTransformToParent()); + if (!volumeToLpsTransform) + { + return; + } + } + else + { + return; + } + vtkMRMLTransformNode* beamTransformNode = nullptr; if (vtkMRMLNode* node = scene->GetFirstNodeByName(RTIMAGE_TRANSFORM_NODE_NAME)) { beamTransformNode = vtkMRMLLinearTransformNode::SafeDownCast(node); } + vtkMRMLScalarVolumeNode* volumeNode = parameterNode->GetReferenceVolumeNode(); + vtkNew ijkToRasDirectionMatrix; + volumeNode->GetIJKToRASDirectionMatrix(ijkToRasDirectionMatrix); + vtkTransform* beamTransform = nullptr; vtkNew mat; // DICOM beam transform matrix mat->Identity(); @@ -1760,18 +1824,37 @@ void vtkSlicerDrrImageComputationLogic::UpdateNormalAndVupVectors(vtkMRMLDrrImag if (beamTransformNode) { beamTransform = vtkTransform::SafeDownCast(beamTransformNode->GetTransformToParent()); + ijkToRasDirectionMatrix->Invert(); vtkNew rasToLpsTransform; rasToLpsTransform->Identity(); rasToLpsTransform->RotateZ(180.0); - - vtkNew dicomBeamTransform; - dicomBeamTransform->Identity(); - dicomBeamTransform->PreMultiply(); - dicomBeamTransform->Concatenate(rasToLpsTransform); - dicomBeamTransform->Concatenate(beamTransform); - dicomBeamTransform->GetMatrix(mat); + vtkNew ijkToLpsDirectionMatrix; + ijkToLpsDirectionMatrix->Identity(); + ijkToLpsDirectionMatrix->PreMultiply(); + ijkToLpsDirectionMatrix->Concatenate(ijkToRasDirectionMatrix); + ijkToLpsDirectionMatrix->Concatenate(rasToLpsTransform); + ijkToLpsDirectionMatrix->Update(); + + if (!IsMatrixTransformsRasToLps(ijkToRasDirectionMatrix)) + { + vtkNew dicomBeamTransform; + dicomBeamTransform->Identity(); + dicomBeamTransform->PreMultiply(); + dicomBeamTransform->Concatenate(ijkToLpsDirectionMatrix); + dicomBeamTransform->Concatenate(beamTransform); + dicomBeamTransform->GetMatrix(mat); + } + else + { + vtkNew dicomBeamTransform; + dicomBeamTransform->Identity(); + dicomBeamTransform->PreMultiply(); + dicomBeamTransform->Concatenate(rasToLpsTransform); + dicomBeamTransform->Concatenate(beamTransform); + dicomBeamTransform->GetMatrix(mat); + } } else { @@ -1779,12 +1862,22 @@ void vtkSlicerDrrImageComputationLogic::UpdateNormalAndVupVectors(vtkMRMLDrrImag } double n[4], vup[4]; - const double normalVector[4] = { 0., 0., 1., 0. }; // beam positive Z-axis - const double viewUpVector[4] = { -1., 0., 0., 0. }; // beam negative X-axis + if (!IsMatrixTransformsRasToLps(ijkToRasDirectionMatrix)) + { + const double normalVector[4] = { 0., 0., -1., 0. }; // beam negative Z-axis + const double viewUpVector[4] = { -1., 0., 0., 0. }; // beam negative X-axis - mat->MultiplyPoint( normalVector, n); - mat->MultiplyPoint( viewUpVector, vup); + mat->MultiplyPoint( normalVector, n); + mat->MultiplyPoint( viewUpVector, vup); + } + else + { + const double normalVector[4] = { 0., 0., 1., 0. }; // beam positive Z-axis + const double viewUpVector[4] = { -1., 0., 0., 0. }; // beam negative X-axis + mat->MultiplyPoint( normalVector, n); + mat->MultiplyPoint( viewUpVector, vup); + } parameterNode->SetNormalVector(n); parameterNode->SetViewUpVector(vup); } @@ -2063,6 +2156,136 @@ vtkMRMLTableNode* vtkSlicerDrrImageComputationLogic::CreateProjectionsTableNode( return tableNode; } +//------------------------------------------------------------------------------ +vtkMRMLLinearTransformNode* vtkSlicerDrrImageComputationLogic::UpdateVolumeToLPSTransform(vtkMRMLDrrImageComputationNode* parameterNode) +{ + if (!parameterNode) + { + vtkErrorMacro("UpdateVolumeToLPSTransform: Invalid parameter node"); + return nullptr; + } + + vtkMRMLRTBeamNode* beamNode = parameterNode->GetBeamNode(); + if (!beamNode) + { + vtkErrorMacro("UpdateVolumeToLPSTransform: Invalid RT Beam node"); + return nullptr; + } + // Get RT plan for beam + vtkMRMLRTPlanNode *planNode = beamNode->GetParentPlanNode(); + if (!planNode) + { + vtkErrorMacro("UpdateVolumeToLPSTransform: Failed to retrieve valid plan node for beam '" << beamNode->GetName() << "'"); + return nullptr; + } + vtkMRMLScalarVolumeNode* refNode = planNode->GetReferenceVolumeNode(); + if (!refNode) + { + vtkErrorMacro("UpdateVolumeToLPSTransform: Failed to retrieve reference volume node from plan '" << planNode->GetName() << "'"); + return nullptr; + } + return this->UpdateVolumeToLPSTransform(refNode); +} + +//------------------------------------------------------------------------------ +vtkMRMLLinearTransformNode* vtkSlicerDrrImageComputationLogic::UpdateVolumeToLPSTransform(vtkMRMLScalarVolumeNode* ctInputVolume) +{ + vtkMRMLScene* scene = this->GetMRMLScene(); + if (!scene) + { + vtkErrorMacro("UpdateVolumeToLPSTransform: Invalid MRML scene"); + return nullptr; + } + + vtkSmartPointer transformNode; + if (!scene->GetFirstNodeByName(VOLUME_TO_LPS_TRANSFORM_NODE_NAME)) + { + transformNode = vtkSmartPointer::New(); + transformNode->SetName(VOLUME_TO_LPS_TRANSFORM_NODE_NAME); + transformNode->SetHideFromEditors(1); + transformNode->SetSingletonTag("VolumeToLPS_Transform"); + scene->AddNode(transformNode); + } + else + { + transformNode = vtkMRMLLinearTransformNode::SafeDownCast( + scene->GetFirstNodeByName(VOLUME_TO_LPS_TRANSFORM_NODE_NAME)); + } + + vtkNew matrix; + if (transformNode) + { + ctInputVolume->GetIJKToRASDirectionMatrix(matrix); + vtkNew transform; + vtkNew rasToLpsTransform; + rasToLpsTransform->Identity(); + rasToLpsTransform->RotateZ(180.0); + + matrix->Invert(); + transform->Identity(); + transform->PostMultiply(); + transform->Concatenate(matrix); // RAS -> IJK + transform->Concatenate(rasToLpsTransform); // RAS -> LPS + transform->Update(); + + transformNode->SetAndObserveTransformToParent(transform); + return transformNode; + } + return nullptr; +} + +//------------------------------------------------------------------------------ +void vtkSlicerDrrImageComputationLogic::GetVolumeToLPSTransformMatrix(vtkMRMLDrrImageComputationNode* parameterNode, vtkMatrix4x4* matrix) +{ + if (!parameterNode) + { + vtkErrorMacro("GetVolumeToLPSTransform: Invalid parameter node"); + return; + } + if (!matrix) + { + vtkErrorMacro("GetVolumeToLPSTransform: Invalid matrix"); + return; + } + vtkMRMLTransformNode* transformNode = this->UpdateVolumeToLPSTransform(parameterNode); + if (!transformNode) + { + vtkErrorMacro("GetVolumeToLPSTransform: Invalid volume to LPS transfrom node"); + return; + } + vtkLinearTransform* linearTransform = vtkLinearTransform::SafeDownCast(transformNode->GetTransformToParent()); + if (linearTransform) + { + linearTransform->GetMatrix(matrix); + } +} + +//------------------------------------------------------------------------------ +void vtkSlicerDrrImageComputationLogic::GetVolumeToLPSTransformMatrix(vtkMRMLScalarVolumeNode* ctVolumeNode, vtkMatrix4x4* matrix) +{ + if (!ctVolumeNode) + { + vtkErrorMacro("GetVolumeToLPSTransform: Invalid CT volume node"); + return; + } + if (!matrix) + { + vtkErrorMacro("GetVolumeToLPSTransform: Invalid matrix"); + return; + } + vtkMRMLTransformNode* transformNode = this->UpdateVolumeToLPSTransform(ctVolumeNode); + if (!transformNode) + { + vtkErrorMacro("GetVolumeToLPSTransform: Invalid volume to LPS transfrom node"); + return; + } + vtkLinearTransform* linearTransform = vtkLinearTransform::SafeDownCast(transformNode->GetTransformToParent()); + if (linearTransform) + { + linearTransform->GetMatrix(matrix); + } +} + //------------------------------------------------------------------------------ void vtkSlicerDrrImageComputationLogic::SetPlanarImageLogic(vtkSlicerPlanarImageModuleLogic* planarImageLogic) { diff --git a/DrrImageComputation/Logic/vtkSlicerDrrImageComputationLogic.h b/DrrImageComputation/Logic/vtkSlicerDrrImageComputationLogic.h index 8e19d72e..c10d2911 100644 --- a/DrrImageComputation/Logic/vtkSlicerDrrImageComputationLogic.h +++ b/DrrImageComputation/Logic/vtkSlicerDrrImageComputationLogic.h @@ -50,6 +50,7 @@ class vtkSlicerBeamsModuleLogic; class vtkSlicerPlanarImageModuleLogic; class vtkSlicerCLIModuleLogic; +class vtkTransform; class vtkMatrix4x4; /// \ingroup Slicer_QtModules_DrrImageComputation @@ -62,7 +63,8 @@ class VTK_SLICER_DRRIMAGECOMPUTATION_MODULE_LOGIC_EXPORT vtkSlicerDrrImageComput static const char* FIDUCIALS_MARKUPS_NODE_NAME; // fiducial static const char* NORMAL_VECTOR_MARKUPS_NODE_NAME; // line static const char* VUP_VECTOR_MARKUPS_NODE_NAME; // line - static const char* RTIMAGE_TRANSFORM_NODE_NAME; + static const char* RTIMAGE_TRANSFORM_NODE_NAME; // RtImage transform node name + static const char* VOLUME_TO_LPS_TRANSFORM_NODE_NAME; // CT Volume to LPS transform node name static vtkSlicerDrrImageComputationLogic *New(); vtkTypeMacro(vtkSlicerDrrImageComputationLogic, vtkSlicerModuleLogic); @@ -71,6 +73,16 @@ class VTK_SLICER_DRRIMAGECOMPUTATION_MODULE_LOGIC_EXPORT vtkSlicerDrrImageComput /// Update normal and view up vectors of RT Image void UpdateNormalAndVupVectors(vtkMRMLDrrImageComputationNode* parameterNode); + /// Update volume RAS to LPS transform node + /// @param parameterNode - input CT volume node + vtkMRMLLinearTransformNode* UpdateVolumeToLPSTransform(vtkMRMLScalarVolumeNode* ctVolumeNode); + /// Update volume RAS to LPS transform node + /// @param parameterNode - input parameter node + vtkMRMLLinearTransformNode* UpdateVolumeToLPSTransform(vtkMRMLDrrImageComputationNode* parameterNode); + /// Get volume RAS to LPS transform matrix + void GetVolumeToLPSTransformMatrix(vtkMRMLDrrImageComputationNode* parameterNode, vtkMatrix4x4* matrix); + void GetVolumeToLPSTransformMatrix(vtkMRMLScalarVolumeNode* ctVolumeNode, vtkMatrix4x4* matrix); + /// Create markups nodes for visualization void CreateMarkupsNodes(vtkMRMLDrrImageComputationNode* parameterNode); /// Update markups nodes using parameter node data diff --git a/DrrImageComputation/MRML/vtkMRMLDrrImageComputationNode.cxx b/DrrImageComputation/MRML/vtkMRMLDrrImageComputationNode.cxx index c295d817..842cbd02 100644 --- a/DrrImageComputation/MRML/vtkMRMLDrrImageComputationNode.cxx +++ b/DrrImageComputation/MRML/vtkMRMLDrrImageComputationNode.cxx @@ -18,6 +18,7 @@ // Beams includes #include +#include // MRML includes #include @@ -84,7 +85,7 @@ vtkMRMLDrrImageComputationNode::vtkMRMLDrrImageComputationNode() InvertIntensityFlag = true; IsocenterImagerDistance = 300.; - HUThresholdBelow = -1000; + HUThresholdBelow = 150; // Observe RTBeam node events (like change of transform or geometry) vtkNew nodeEvents; @@ -394,6 +395,20 @@ void vtkMRMLDrrImageComputationNode::SetThreading(int threading) }; } +//---------------------------------------------------------------------------- +void vtkMRMLDrrImageComputationNode::GetIsocenterPositionRAS(double isocenterPosition[3]) +{ + vtkMRMLRTBeamNode* beamNode = this->GetBeamNode(); + if (!beamNode) + { + vtkErrorMacro("GetIsocenterPositionLPS: RT Beam node is invalid"); + return; + } + + // Isocenter RAS position, for plastimatch isocenter MUST BE in LPS system + beamNode->GetPlanIsocenterPosition(isocenterPosition); +} + //---------------------------------------------------------------------------- void vtkMRMLDrrImageComputationNode::GetIsocenterPositionLPS(double isocenterPosition[3]) { @@ -436,3 +451,28 @@ void vtkMRMLDrrImageComputationNode::GetImageCenter(double imageCenter[2]) image_center[1] = static_cast(imageCenter[1]); this->SetImageCenter(image_center); } + +//---------------------------------------------------------------------------- +vtkMRMLScalarVolumeNode* vtkMRMLDrrImageComputationNode::GetReferenceVolumeNode() +{ + vtkMRMLRTBeamNode* beamNode = this->GetBeamNode(); + if (!beamNode) + { + vtkErrorMacro("GetReferenceVolumeNode: Invalid RT Beam node"); + return nullptr; + } + // Get RT plan for beam + vtkMRMLRTPlanNode *planNode = beamNode->GetParentPlanNode(); + if (!planNode) + { + vtkErrorMacro("GetReferenceVolumeNode: Failed to retrieve valid plan node for beam '" << beamNode->GetName() << "'"); + return nullptr; + } + vtkMRMLScalarVolumeNode* refNode = planNode->GetReferenceVolumeNode(); + if (!refNode) + { + vtkErrorMacro("GetReferenceVolumeNode: Failed to retrieve reference volume node from plan '" << planNode->GetName() << "'"); + return nullptr; + } + return refNode; +} diff --git a/DrrImageComputation/MRML/vtkMRMLDrrImageComputationNode.h b/DrrImageComputation/MRML/vtkMRMLDrrImageComputationNode.h index ad990e69..4de8dc6b 100644 --- a/DrrImageComputation/MRML/vtkMRMLDrrImageComputationNode.h +++ b/DrrImageComputation/MRML/vtkMRMLDrrImageComputationNode.h @@ -76,9 +76,11 @@ class VTK_SLICER_DRRIMAGECOMPUTATION_MODULE_MRML_EXPORT vtkMRMLDrrImageComputati void GetRTImagePosition(double position[2]); void GetIsocenterPositionLPS(double position[3]); + void GetIsocenterPositionRAS(double position[3]); /// Handles events registered in the observer manager void ProcessMRMLEvents(vtkObject *caller, unsigned long eventID, void *callData) override; + vtkMRMLScalarVolumeNode* GetReferenceVolumeNode(); public: /// Get beam node diff --git a/DrrImageComputation/Resources/UI/qSlicerDrrImageComputationModuleWidget.ui b/DrrImageComputation/Resources/UI/qSlicerDrrImageComputationModuleWidget.ui index 56e8f90b..f335d66f 100644 --- a/DrrImageComputation/Resources/UI/qSlicerDrrImageComputationModuleWidget.ui +++ b/DrrImageComputation/Resources/UI/qSlicerDrrImageComputationModuleWidget.ui @@ -6,8 +6,8 @@ 0 0 - 458 - 859 + 535 + 951 @@ -96,13 +96,6 @@ - - - - CT volume: - - - @@ -127,37 +120,6 @@ - - - - - vtkMRMLRTBeamNode - - - - false - - - true - - - false - - - false - - - Please select a RTBeam node - - - - - - - Camera: - - - @@ -196,6 +158,86 @@ + + + + Volume direction to LPS transform matrix + + + true + + + true + + + + + + + + -1.000000000000000 + + + 1.000000000000000 + + + 4 + + + 0.000100000000000 + + + + + + + + + Apply transform to volume + + + + + + + + + + + vtkMRMLRTBeamNode + + + + false + + + true + + + false + + + false + + + Please select a RTBeam node + + + + + + + CT volume: + + + + + + + Camera: + + + @@ -443,10 +485,10 @@ - + false - + false @@ -565,11 +607,22 @@
ctkCollapsibleButton.h
1 + + ctkCollapsibleGroupBox + QGroupBox +
ctkCollapsibleGroupBox.h
+ 1 + ctkCoordinatesWidget QWidget
ctkCoordinatesWidget.h
 + + ctkMatrixWidget + QWidget +
ctkMatrixWidget.h
+ ctkRangeWidget QWidget @@ -596,6 +649,7 @@ qSlicerSimpleMarkupsWidget qSlicerWidget
qSlicerSimpleMarkupsWidget.h
+ 1 qSlicerDrrImageComputationPlastimatchParametersWidget diff --git a/DrrImageComputation/qSlicerDrrImageComputationModuleWidget.cxx b/DrrImageComputation/qSlicerDrrImageComputationModuleWidget.cxx index e6e4097f..820e1592 100644 --- a/DrrImageComputation/qSlicerDrrImageComputationModuleWidget.cxx +++ b/DrrImageComputation/qSlicerDrrImageComputationModuleWidget.cxx @@ -30,11 +30,15 @@ #include #include #include +#include // SlicerRT MRML Beams includes #include #include +// SlicerRT includes +#include + // SlicerRT MRML DrrImageComputation includes #include "vtkMRMLDrrImageComputationNode.h" @@ -42,6 +46,7 @@ #include "vtkSlicerDrrImageComputationLogic.h" // VTK includes +#include #include #include #include @@ -59,6 +64,33 @@ constexpr int PROJECTED_POINT_ROW_COLUMN = 7; constexpr int PROJECTED_POINT_STATUS_COLUMN = 8; constexpr int PROJECTED_POINT_COLUMNS = 9; +bool IsIdentityMatrix(vtkMatrix4x4* mat) +{ + if (!mat) + { + return false; + } +// constexpr double local_epsilon = std::numeric_limits< double >::epsilon(); + constexpr double local_epsilon = EPSILON; + bool result = bool(fabs(mat->GetElement(0, 0) - 1.) <= local_epsilon); + result &= bool(fabs(mat->GetElement(1, 1) - 1.) <= local_epsilon); + result &= bool(fabs(mat->GetElement(2, 2) - 1.) <= local_epsilon); + result &= bool(fabs(mat->GetElement(3, 3) - 1.) <= local_epsilon); + result &= bool(fabs(mat->GetElement(1, 0)) <= local_epsilon); + result &= bool(fabs(mat->GetElement(2, 0)) <= local_epsilon); + result &= bool(fabs(mat->GetElement(2, 1)) <= local_epsilon); + result &= bool(fabs(mat->GetElement(3, 0)) <= local_epsilon); + result &= bool(fabs(mat->GetElement(3, 1)) <= local_epsilon); + result &= bool(fabs(mat->GetElement(3, 2)) <= local_epsilon); + result &= bool(fabs(mat->GetElement(0, 1)) <= local_epsilon); + result &= bool(fabs(mat->GetElement(0, 2)) <= local_epsilon); + result &= bool(fabs(mat->GetElement(1, 2)) <= local_epsilon); + result &= bool(fabs(mat->GetElement(0, 3)) <= local_epsilon); + result &= bool(fabs(mat->GetElement(1, 3)) <= local_epsilon); + result &= bool(fabs(mat->GetElement(2, 3)) <= local_epsilon); + return result; +} + }; //----------------------------------------------------------------------------- @@ -217,6 +249,7 @@ void qSlicerDrrImageComputationModuleWidget::setup() connect( d->PushButton_UpdateBeamFromCamera, SIGNAL(clicked()), this, SLOT(onUpdateBeamFromCameraClicked())); connect( d->PushButton_ProjectControlPoints, SIGNAL(clicked()), this, SLOT(onProjectMarkupsControlPointsClicked())); connect( d->PushButton_ClearProjectedTableWidget, SIGNAL(clicked()), this, SLOT(onClearProjectedTableClicked())); + connect( d->PushButton_ApplyToVolumeTransform, SIGNAL(clicked()), this, SLOT(onApplyVolumeToLpsTransformToVolumeClicked())); // Handle scene change event if occurs qvtkConnect( d->logic(), vtkCommand::ModifiedEvent, this, SLOT(onLogicModified())); @@ -373,6 +406,19 @@ void qSlicerDrrImageComputationModuleWidget::updateWidgetFromMRML() res = d->logic()->GetPlastimatchProjectionMatrix(parameterNode, mat); d->PlastimatchParametersWidget->setEnabled(res); d->PlastimatchParametersWidget->setProjectionMatrix(mat); + // Volume IJKtoRAS to LPS direction matrix transform + mat->Identity(); + d->logic()->GetVolumeToLPSTransformMatrix(parameterNode, mat); + + for (int i = 0; i < 4; i++) + { + for (int j = 0; j < 4; j++) + { + double v = mat->GetElement( i, j ); + d->MatrixWidget_VolumeDirectionToLpsTransformMatrix->setValue( i, j, v ); + } + } + d->PushButton_ApplyToVolumeTransform->setEnabled(!IsIdentityMatrix(mat)); } //----------------------------------------------------------------------------- @@ -454,6 +500,21 @@ void qSlicerDrrImageComputationModuleWidget::onCtVolumeNodeChanged(vtkMRMLNode* qCritical() << Q_FUNC_INFO << ": Invalid reference CT volume node"; return; } + + // Volume IJKtoRAS to LPS direction matrix transform + vtkNew mat; + mat->Identity(); + d->logic()->GetVolumeToLPSTransformMatrix(volumeNode, mat); + + for (int i = 0; i < 4; i++) + { + for (int j = 0; j < 4; j++) + { + double v = mat->GetElement( i, j ); + d->MatrixWidget_VolumeDirectionToLpsTransformMatrix->setValue( i, j, v ); + } + } + d->PushButton_ApplyToVolumeTransform->setEnabled(!IsIdentityMatrix(mat)); } //----------------------------------------------------------------------------- @@ -695,6 +756,41 @@ void qSlicerDrrImageComputationModuleWidget::onUseImageWindowToggled(bool value) parameterNode->SetImageWindowFlag(value); // Update imager and image markups, DRR arguments } +//----------------------------------------------------------------------------- +void qSlicerDrrImageComputationModuleWidget::onApplyVolumeToLpsTransformToVolumeClicked() +{ + Q_D(qSlicerDrrImageComputationModuleWidget); + vtkMRMLDrrImageComputationNode* parameterNode = vtkMRMLDrrImageComputationNode::SafeDownCast(d->MRMLNodeComboBox_ParameterSet->currentNode()); + vtkMRMLScalarVolumeNode* ctVolumeNode = vtkMRMLScalarVolumeNode::SafeDownCast(d->MRMLNodeComboBox_CtVolume->currentNode()); + + if (!parameterNode) + { + qCritical() << Q_FUNC_INFO << ": Invalid parameter node"; + return; + } + + if (!ctVolumeNode) + { + qCritical() << Q_FUNC_INFO << ": Invalid CT volume node"; + return; + } + vtkMRMLLinearTransformNode* transformNode = d->logic()->UpdateVolumeToLPSTransform(ctVolumeNode); + if (transformNode) + { + vtkNew< vtkMatrix4x4 > mat; + mat->Identity(); + vtkLinearTransform* linearTransform = vtkLinearTransform::SafeDownCast(transformNode->GetTransformToParent()); + if (linearTransform) + { + linearTransform->GetMatrix(mat); + } + if (!mat->IsIdentity()) + { + ctVolumeNode->SetAndObserveTransformNodeID(transformNode->GetID()); + } + } +} + //----------------------------------------------------------------------------- void qSlicerDrrImageComputationModuleWidget::onComputeDrrClicked() { @@ -721,7 +817,13 @@ void qSlicerDrrImageComputationModuleWidget::onComputeDrrClicked() } QApplication::setOverrideCursor(Qt::WaitCursor); - +/* + vtkNew transform; + if (d->logic()->GetVolumeToLPSTransform(parameterNode, transform, ctVolumeNode)) + { + qDebug() << Q_FUNC_INFO << "Harden the transform"; + } +*/ vtkMRMLScalarVolumeNode* drrImageNode = d->logic()->ComputePlastimatchDRR( parameterNode, ctVolumeNode); if (drrImageNode) { @@ -729,6 +831,7 @@ void qSlicerDrrImageComputationModuleWidget::onComputeDrrClicked() QApplication::restoreOverrideCursor(); return; } + QApplication::restoreOverrideCursor(); } diff --git a/DrrImageComputation/qSlicerDrrImageComputationModuleWidget.h b/DrrImageComputation/qSlicerDrrImageComputationModuleWidget.h index 413f4dd9..d6fe6407 100644 --- a/DrrImageComputation/qSlicerDrrImageComputationModuleWidget.h +++ b/DrrImageComputation/qSlicerDrrImageComputationModuleWidget.h @@ -70,6 +70,7 @@ public slots: void onProjectMarkupsControlPointsClicked(); void onClearProjectedTableClicked(); + void onApplyVolumeToLpsTransformToVolumeClicked(); /// Update widget GUI from RT Image parameters node void updateWidgetFromMRML(); diff --git a/PlmDrr/plastimatch_slicer_drr.cxx b/PlmDrr/plastimatch_slicer_drr.cxx index a1ceda98..0a4d2d89 100644 --- a/PlmDrr/plastimatch_slicer_drr.cxx +++ b/PlmDrr/plastimatch_slicer_drr.cxx @@ -22,10 +22,15 @@ #include #include #include +#include +#include #include #include #include #include +#include +#include +#include // Plastimatch includes #include @@ -183,6 +188,45 @@ int DoIt( int argc, char * argv[], Drr_options& options, TPixel ) // Apply threshold filter if HU thresholdBelow is higher than -1000 typename InputImageType::Pointer inputImagePointer = inputReader->GetOutput(); + typename InputImageType::PointType origin = inputImagePointer->GetOrigin(); + typename InputImageType::DirectionType originalDirection = inputImagePointer->GetDirection(); + + // Transform origin and orientation into LPS + using FilterType = itk::ChangeInformationImageFilter; + auto filter = FilterType::New(); + filter->SetInput(inputImagePointer); + + typename InputImageType::DirectionType direction = inputImagePointer->GetDirection().GetInverse(); + + using AffineTransformType = itk::AffineTransform; + auto transform = AffineTransformType::New(); + transform->SetMatrix(direction); + + typename InputImageType::PointType originLPS = transform->TransformPoint(origin); + options.isocenter[0] *= -1.; + options.isocenter[1] *= -1.; + typename InputImageType::PointType isocenterLPS = transform->TransformPoint(options.isocenter); + options.isocenter[0] = isocenterLPS[0]; + options.isocenter[1] = isocenterLPS[1]; + options.isocenter[2] = isocenterLPS[2]; + + transform->SetIdentity(); + typename InputImageType::DirectionType identity = transform->GetMatrix(); // RAI orientation, unity matrix + + filter->SetOutputOrigin(originLPS); + filter->ChangeOriginOn(); + filter->SetOutputDirection(identity); + filter->ChangeDirectionOn(); + try + { + filter->UpdateOutputInformation(); + inputImagePointer = filter->GetOutput(); + } + catch ( itk::ExceptionObject& excep ) + { + throw; + } + using ThresholdFilterType = itk::ThresholdImageFilter< InputImageType >; typename ThresholdFilterType::Pointer thresholdFilter = ThresholdFilterType::New(); @@ -305,6 +349,32 @@ int DoIt( int argc, char * argv[], Drr_options& options, TPixel ) rescale->SetOutputMaximum(options.autoscale_range[1]); rescale->SetInput(drrReader->GetOutput()); + // flip image data in case of original image direction (orientation) isn't LPS + typename InputImageType::DirectionType identity; + identity.SetIdentity(); + + using FlipImageFilterType = itk::FlipImageFilter< PlmDrrImageType >; + FlipImageFilterType::FlipAxesArrayType flipAxes; + auto flipFilter = FlipImageFilterType::New(); + + if (originalDirection != identity) + { + flipAxes[0] = false; + flipAxes[1] = true; + flipFilter->SetFlipAxes(flipAxes); + flipFilter->SetInput(rescale->GetOutput()); + std::cout << "Original volume orientation isn't LPS" << std::endl; + + try + { + flipFilter->Update(); + } + catch ( itk::ExceptionObject& excep ) + { + throw; + } + } + // write data into Slicer using WriterType = itk::ImageFileWriter< PlmDrrImageType >; WriterType::Pointer writer = WriterType::New(); @@ -316,7 +386,14 @@ int DoIt( int argc, char * argv[], Drr_options& options, TPixel ) // invert using InvertFilterType = itk::InvertIntensityImageFilter< PlmDrrImageType, PlmDrrImageType >; InvertFilterType::Pointer invert = InvertFilterType::New(); - invert->SetInput(rescale->GetOutput()); + if (originalDirection != identity) // Flip image filter if original volume orientation isn't LPS + { + invert->SetInput( flipFilter->GetOutput() ); + } + else + { + invert->SetInput( rescale->GetOutput() ); + } invert->SetMaximum(options.autoscale_range[0]); // inverted input @@ -333,7 +410,14 @@ int DoIt( int argc, char * argv[], Drr_options& options, TPixel ) } else { - writer->SetInput( rescale->GetOutput() ); + if (originalDirection != identity) // Flip image filter if original volume orientation isn't LPS + { + writer->SetInput( flipFilter->GetOutput() ); + } + else + { + writer->SetInput( rescale->GetOutput() ); + } try {