diff --git a/Modules/IGTUI/Qmitk/QmitkInteractiveTransformationWidgetControls.ui b/Modules/IGTUI/Qmitk/QmitkInteractiveTransformationWidgetControls.ui index 6daf485458..659f57055c 100644 --- a/Modules/IGTUI/Qmitk/QmitkInteractiveTransformationWidgetControls.ui +++ b/Modules/IGTUI/Qmitk/QmitkInteractiveTransformationWidgetControls.ui @@ -1,463 +1,463 @@ QmitkInteractiveTransformationWidgetControls 0 0 435 381 Form Interactive Translation 50 false x-Direction (Frontal): false 0 0 ArrowCursor -500 500 1 Qt::Horizontal false false QSlider::TicksAbove 100 - 1 + 4 -500.000000000000000 500.000000000000000 50 false y-Direction (Sagittal): false 0 0 -500 500 Qt::Horizontal false QSlider::TicksAbove 100 - 1 + 4 -500.000000000000000 500.000000000000000 50 false z-Direction (Transversal): false 0 0 -500 500 0 Qt::Horizontal false QSlider::TicksAbove 100 - 1 + 4 -500.000000000000000 500.000000000000000 Interactive Rotation 50 false x-Axis (Frontal): false 0 0 ArrowCursor -180 180 Qt::Horizontal QSlider::TicksAbove 45 - 1 + 4 -180.000000000000000 180.000000000000000 50 false y-Axis (Sagittal): false 0 0 -180 180 Qt::Horizontal QSlider::TicksAbove 45 - 1 + 4 -180.000000000000000 180.000000000000000 50 false z-Axis (Transversal): false 0 0 -180 180 Qt::Horizontal QSlider::TicksAbove 45 - 1 + 4 -180.000000000000000 180.000000000000000 Use Manipulated ToolTip Reset To Identity Revert Changes Cancel diff --git a/Plugins/org.mitk.gui.qt.igt.app.echotrack/src/internal/NavigationStepWidgets/QmitkUSNavigationStepCtUsRegistration.cpp b/Plugins/org.mitk.gui.qt.igt.app.echotrack/src/internal/NavigationStepWidgets/QmitkUSNavigationStepCtUsRegistration.cpp index cdfe37e46f..6f4da18cab 100644 --- a/Plugins/org.mitk.gui.qt.igt.app.echotrack/src/internal/NavigationStepWidgets/QmitkUSNavigationStepCtUsRegistration.cpp +++ b/Plugins/org.mitk.gui.qt.igt.app.echotrack/src/internal/NavigationStepWidgets/QmitkUSNavigationStepCtUsRegistration.cpp @@ -1,2331 +1,2369 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "QmitkUSNavigationStepCtUsRegistration.h" #include "ui_QmitkUSNavigationStepCtUsRegistration.h" #include #include "mitkNodeDisplacementFilter.h" #include "../QmitkUSNavigationMarkerPlacement.h" #include #include #include #include "mitkProperties.h" #include #include #include #include #include #include #include #include #include #include #include #include static const int NUMBER_FIDUCIALS_NEEDED = 8; QmitkUSNavigationStepCtUsRegistration::QmitkUSNavigationStepCtUsRegistration(QWidget *parent) : QmitkUSAbstractNavigationStep(parent), ui(new Ui::QmitkUSNavigationStepCtUsRegistration), m_PerformingGroundTruthProtocolEvaluation(false) { this->UnsetFloatingImageGeometry(); this->DefineDataStorageImageFilter(); this->CreateQtPartControl(this); + this->InitializeTransformationSensorCSToMarkerCS(); } QmitkUSNavigationStepCtUsRegistration::~QmitkUSNavigationStepCtUsRegistration() { delete ui; } bool QmitkUSNavigationStepCtUsRegistration::OnStartStep() { MITK_INFO << "OnStartStep()"; return true; } bool QmitkUSNavigationStepCtUsRegistration::OnStopStep() { MITK_INFO << "OnStopStep()"; return true; } bool QmitkUSNavigationStepCtUsRegistration::OnFinishStep() { MITK_INFO << "OnFinishStep()"; return true; } bool QmitkUSNavigationStepCtUsRegistration::OnActivateStep() { MITK_INFO << "OnActivateStep()"; ui->floatingImageComboBox->SetDataStorage(this->GetDataStorage()); ui->groundTruthImageComboBox->SetDataStorage(this->GetDataStorage()); ui->ctImagesToChooseComboBox->SetDataStorage(this->GetDataStorage()); return true; } bool QmitkUSNavigationStepCtUsRegistration::OnDeactivateStep() { MITK_INFO << "OnDeactivateStep()"; return true; } void QmitkUSNavigationStepCtUsRegistration::OnUpdate() { if (m_NavigationDataSource.IsNull()) { return; } m_NavigationDataSource->Update(); } void QmitkUSNavigationStepCtUsRegistration::OnSettingsChanged(const itk::SmartPointer settingsNode) { } QString QmitkUSNavigationStepCtUsRegistration::GetTitle() { return "CT-to-US registration"; } QmitkUSAbstractNavigationStep::FilterVector QmitkUSNavigationStepCtUsRegistration::GetFilter() { return FilterVector(); } void QmitkUSNavigationStepCtUsRegistration::OnSetCombinedModality() { mitk::AbstractUltrasoundTrackerDevice::Pointer combinedModality = this->GetCombinedModality(false); if (combinedModality.IsNotNull()) { m_NavigationDataSource = combinedModality->GetNavigationDataSource(); } } void QmitkUSNavigationStepCtUsRegistration::UnsetFloatingImageGeometry() { m_ImageDimension[0] = 0; m_ImageDimension[1] = 0; m_ImageDimension[2] = 0; m_ImageSpacing[0] = 1; m_ImageSpacing[1] = 1; m_ImageSpacing[2] = 1; } void QmitkUSNavigationStepCtUsRegistration::SetFloatingImageGeometryInformation(mitk::Image * image) { m_ImageDimension[0] = image->GetDimension(0); m_ImageDimension[1] = image->GetDimension(1); m_ImageDimension[2] = image->GetDimension(2); m_ImageSpacing[0] = image->GetGeometry()->GetSpacing()[0]; m_ImageSpacing[1] = image->GetGeometry()->GetSpacing()[1]; m_ImageSpacing[2] = image->GetGeometry()->GetSpacing()[2]; } double QmitkUSNavigationStepCtUsRegistration::GetVoxelVolume() { if (m_FloatingImage.IsNull()) { return 0.0; } MITK_INFO << "ImageSpacing = " << m_ImageSpacing; return m_ImageSpacing[0] * m_ImageSpacing[1] * m_ImageSpacing[2]; } double QmitkUSNavigationStepCtUsRegistration::GetFiducialVolume(double radius) { return 1.333333333 * 3.141592 * (radius * radius * radius); } bool QmitkUSNavigationStepCtUsRegistration::FilterFloatingImage() { if (m_FloatingImage.IsNull()) { return false; } ImageType::Pointer itkImage1 = ImageType::New(); mitk::CastToItkImage(m_FloatingImage, itkImage1); this->InitializeImageFilters(); m_ThresholdFilter->SetInput(itkImage1); m_LaplacianFilter1->SetInput(m_ThresholdFilter->GetOutput()); m_LaplacianFilter2->SetInput(m_LaplacianFilter1->GetOutput()); m_BinaryThresholdFilter->SetInput(m_LaplacianFilter2->GetOutput()); m_HoleFillingFilter->SetInput(m_BinaryThresholdFilter->GetOutput()); m_BinaryImageToShapeLabelMapFilter->SetInput(m_HoleFillingFilter->GetOutput()); m_BinaryImageToShapeLabelMapFilter->Update(); ImageType::Pointer binaryImage = ImageType::New(); binaryImage = m_HoleFillingFilter->GetOutput(); this->EliminateTooSmallLabeledObjects(binaryImage); //mitk::CastToMitkImage(binaryImage, m_FloatingImage); return true; } void QmitkUSNavigationStepCtUsRegistration::InitializeImageFilters() { //Initialize threshold filters m_ThresholdFilter = itk::ThresholdImageFilter::New(); m_ThresholdFilter->SetOutsideValue(0); m_ThresholdFilter->SetLower(500); m_ThresholdFilter->SetUpper(3200); //Initialize binary threshold filter 1 m_BinaryThresholdFilter = BinaryThresholdImageFilterType::New(); m_BinaryThresholdFilter->SetOutsideValue(0); m_BinaryThresholdFilter->SetInsideValue(1); m_BinaryThresholdFilter->SetLowerThreshold(350); m_BinaryThresholdFilter->SetUpperThreshold(10000); //Initialize laplacian recursive gaussian image filter m_LaplacianFilter1 = LaplacianRecursiveGaussianImageFilterType::New(); m_LaplacianFilter2 = LaplacianRecursiveGaussianImageFilterType::New(); //Initialize binary hole filling filter m_HoleFillingFilter = VotingBinaryIterativeHoleFillingImageFilterType::New(); VotingBinaryIterativeHoleFillingImageFilterType::InputSizeType radius; radius.Fill(1); m_HoleFillingFilter->SetRadius(radius); m_HoleFillingFilter->SetBackgroundValue(0); m_HoleFillingFilter->SetForegroundValue(1); m_HoleFillingFilter->SetMaximumNumberOfIterations(5); //Initialize binary image to shape label map filter m_BinaryImageToShapeLabelMapFilter = BinaryImageToShapeLabelMapFilterType::New(); m_BinaryImageToShapeLabelMapFilter->SetInputForegroundValue(1); } double QmitkUSNavigationStepCtUsRegistration::GetCharacteristicDistanceAWithUpperMargin() { switch (ui->fiducialMarkerConfigurationComboBox->currentIndex()) { // case 0 is equal to fiducial marker configuration A (10mm distance) case 0: return 12.07; // case 1 is equal to fiducial marker configuration B (15mm distance) case 1: return 18.105; // case 2 is equal to fiducial marker configuration C (20mm distance) case 2: return 24.14; } return 0.0; } double QmitkUSNavigationStepCtUsRegistration::GetCharacteristicDistanceBWithLowerMargin() { switch (ui->fiducialMarkerConfigurationComboBox->currentIndex()) { // case 0 is equal to fiducial marker configuration A (10mm distance) case 0: return 12.07; // case 1 is equal to fiducial marker configuration B (15mm distance) case 1: return 18.105; // case 2 is equal to fiducial marker configuration C (20mm distance) case 2: return 24.14; } return 0.0; } double QmitkUSNavigationStepCtUsRegistration::GetCharacteristicDistanceBWithUpperMargin() { switch (ui->fiducialMarkerConfigurationComboBox->currentIndex()) { // case 0 is equal to fiducial marker configuration A (10mm distance) case 0: return 15.73; // case 1 is equal to fiducial marker configuration B (15mm distance) case 1: return 23.595; // case 2 is equal to fiducial marker configuration C (20mm distance) case 2: return 31.46; } return 0.0; } double QmitkUSNavigationStepCtUsRegistration::GetMinimalFiducialConfigurationDistance() { switch (ui->fiducialMarkerConfigurationComboBox->currentIndex()) { // case 0 is equal to fiducial marker configuration A (10mm distance) case 0: return 10.0; // case 1 is equal to fiducial marker configuration B (15mm distance) case 1: return 15.0; // case 2 is equal to fiducial marker configuration C (20mm distance) case 2: return 20.0; } return 0.0; } void QmitkUSNavigationStepCtUsRegistration::CreateMarkerModelCoordinateSystemPointSet() { if (m_MarkerModelCoordinateSystemPointSet.IsNull()) { m_MarkerModelCoordinateSystemPointSet = mitk::PointSet::New(); } else { m_MarkerModelCoordinateSystemPointSet->Clear(); } mitk::Point3D fiducial1; mitk::Point3D fiducial2; mitk::Point3D fiducial3; mitk::Point3D fiducial4; mitk::Point3D fiducial5; mitk::Point3D fiducial6; mitk::Point3D fiducial7; mitk::Point3D fiducial8; switch (ui->fiducialMarkerConfigurationComboBox->currentIndex()) { // case 0 is equal to fiducial marker configuration A (10mm distance) case 0: fiducial1[0] = 0; fiducial1[1] = 0; fiducial1[2] = 0; fiducial2[0] = 0; fiducial2[1] = 10; fiducial2[2] = 0; fiducial3[0] = 10; fiducial3[1] = 0; fiducial3[2] = 0; fiducial4[0] = 20; fiducial4[1] = 20; fiducial4[2] = 0; fiducial5[0] = 0; fiducial5[1] = 20; fiducial5[2] = 10; fiducial6[0] = 10; fiducial6[1] = 20; fiducial6[2] = 10; fiducial7[0] = 20; fiducial7[1] = 10; fiducial7[2] = 10; fiducial8[0] = 20; fiducial8[1] = 0; fiducial8[2] = 10; break; // case 1 is equal to fiducial marker configuration B (15mm distance) case 1: fiducial1[0] = 0; fiducial1[1] = 0; fiducial1[2] = 0; fiducial2[0] = 0; fiducial2[1] = 15; fiducial2[2] = 0; fiducial3[0] = 15; fiducial3[1] = 0; fiducial3[2] = 0; fiducial4[0] = 30; fiducial4[1] = 30; fiducial4[2] = 0; fiducial5[0] = 0; fiducial5[1] = 30; fiducial5[2] = 15; fiducial6[0] = 15; fiducial6[1] = 30; fiducial6[2] = 15; fiducial7[0] = 30; fiducial7[1] = 15; fiducial7[2] = 15; fiducial8[0] = 30; fiducial8[1] = 0; fiducial8[2] = 15; break; // case 2 is equal to fiducial marker configuration C (20mm distance) case 2: fiducial1[0] = 0; fiducial1[1] = 0; fiducial1[2] = 0; fiducial2[0] = 0; fiducial2[1] = 20; fiducial2[2] = 0; fiducial3[0] = 20; fiducial3[1] = 0; fiducial3[2] = 0; fiducial4[0] = 40; fiducial4[1] = 40; fiducial4[2] = 0; fiducial5[0] = 0; fiducial5[1] = 40; fiducial5[2] = 20; fiducial6[0] = 20; fiducial6[1] = 40; fiducial6[2] = 20; fiducial7[0] = 40; fiducial7[1] = 20; fiducial7[2] = 20; fiducial8[0] = 40; fiducial8[1] = 0; fiducial8[2] = 20; break; } m_MarkerModelCoordinateSystemPointSet->InsertPoint(0, fiducial1); m_MarkerModelCoordinateSystemPointSet->InsertPoint(1, fiducial2); m_MarkerModelCoordinateSystemPointSet->InsertPoint(2, fiducial3); m_MarkerModelCoordinateSystemPointSet->InsertPoint(3, fiducial4); m_MarkerModelCoordinateSystemPointSet->InsertPoint(4, fiducial5); m_MarkerModelCoordinateSystemPointSet->InsertPoint(5, fiducial6); m_MarkerModelCoordinateSystemPointSet->InsertPoint(6, fiducial7); m_MarkerModelCoordinateSystemPointSet->InsertPoint(7, fiducial8); /*mitk::DataNode::Pointer node = this->GetDataStorage()->GetNamedNode("Marker Model Coordinate System Point Set"); if (node == nullptr) { node = mitk::DataNode::New(); node->SetName("Marker Model Coordinate System Point Set"); node->SetData(m_MarkerModelCoordinateSystemPointSet); this->GetDataStorage()->Add(node); } else { node->SetData(m_MarkerModelCoordinateSystemPointSet); this->GetDataStorage()->Modified(); }*/ } void QmitkUSNavigationStepCtUsRegistration::InitializeTransformationSensorCSToMarkerCS() { - //At first calculate the inverse transformation (marker-CS to sensor-CS) + //The following calculations are related to the 3mm | 15mm fiducial configuration + m_TransformSensorCSToMarkerCS = mitk::AffineTransform3D::New(); + mitk::Vector3D translation; + translation[0] = -18.175; + translation[1] = 15.000; + translation[2] = 8.001; // considering a distance from the base plate of 0.315 inch and not 0.313 inch + + m_TransformSensorCSToMarkerCS->SetOffset(translation); + + // Quaternion (x, y, z, r) --> n = (1,0,0) --> q(sin(90°),0,0,cos(90°)) + mitk::Quaternion q1(1, 0, 0, 0); // corresponding to a rotation of 180° around the normal x-axis. + // .transpose() is needed for changing the rows and the columns of the returned rotation_matrix_transpose + vnl_matrix_fixed vnl_rotation = q1.rotation_matrix_transpose().transpose(); // :-) + mitk::Matrix3D rotationMatrix; + + for (int i = 0; i < 3; ++i) { + for (int j = 0; j < 3; ++j) { + rotationMatrix[i][j] = vnl_rotation[i][j]; + } + } + + m_TransformSensorCSToMarkerCS->SetMatrix(rotationMatrix); + //The transformation from the sensor-CS to the marker-CS is calculated now. + MITK_INFO << "TransformSensorCSToMarkerCS = " << m_TransformSensorCSToMarkerCS; + + /*mitk::NavigationData::Pointer navData = mitk::NavigationData::New(); + navData->SetOrientation(q1); + navData->SetPosition(translation); + navData->SetHasOrientation(true); + navData->SetHasPosition(true);*/ } void QmitkUSNavigationStepCtUsRegistration::InitializePointsToTransformForGroundTruthProtocol() { m_PointsToTransformGroundTruthProtocol.clear(); mitk::Point3D point0mm; mitk::Point3D point20mm; mitk::Point3D point40mm; mitk::Point3D point60mm; mitk::Point3D point80mm; mitk::Point3D point100mm; point0mm[0] = 0.0; point0mm[1] = 0.0; point0mm[2] = 0.0; point20mm[0] = 0.0; point20mm[1] = 0.0; point20mm[2] = 0.0; point40mm[0] = 0.0; point40mm[1] = 0.0; point40mm[2] = 0.0; point60mm[0] = 0.0; point60mm[1] = 0.0; point60mm[2] = 0.0; point80mm[0] = 0.0; point80mm[1] = 0.0; point80mm[2] = 0.0; point100mm[0] = 0.0; point100mm[1] = 0.0; point100mm[2] = 0.0; m_PointsToTransformGroundTruthProtocol.insert(std::pair(0, point0mm)); m_PointsToTransformGroundTruthProtocol.insert(std::pair(20, point20mm)); m_PointsToTransformGroundTruthProtocol.insert(std::pair(40, point40mm)); m_PointsToTransformGroundTruthProtocol.insert(std::pair(60, point60mm)); m_PointsToTransformGroundTruthProtocol.insert(std::pair(80, point80mm)); m_PointsToTransformGroundTruthProtocol.insert(std::pair(100, point100mm)); } void QmitkUSNavigationStepCtUsRegistration::CreatePointsToTransformForGroundTruthProtocol() { this->InitializePointsToTransformForGroundTruthProtocol(); switch (ui->fiducialMarkerConfigurationComboBox->currentIndex()) { // case 0 is equal to fiducial marker configuration A (10mm distance) case 0: MITK_WARN << "For this marker configuration (10mm) there does not exist a point to transform."; break; // case 1 is equal to fiducial marker configuration B (15mm distance) case 1: m_PointsToTransformGroundTruthProtocol.at(0)[0] = 130; // = 30mm to end of clipping plate + 100 mm to middle axis of measurement plate m_PointsToTransformGroundTruthProtocol.at(0)[1] = 15; m_PointsToTransformGroundTruthProtocol.at(0)[2] = -7; // = 5mm distance to clipping plate + 2mm to base m_PointsToTransformGroundTruthProtocol.at(20)[0] = 130; m_PointsToTransformGroundTruthProtocol.at(20)[1] = 15; m_PointsToTransformGroundTruthProtocol.at(20)[2] = -27; // = 5mm distance to clipping plate + 2mm to base + 20mm depth m_PointsToTransformGroundTruthProtocol.at(40)[0] = 130; m_PointsToTransformGroundTruthProtocol.at(40)[1] = 15; m_PointsToTransformGroundTruthProtocol.at(40)[2] = -47; // = 5mm distance to clipping plate + 2mm to base + 40mm depth m_PointsToTransformGroundTruthProtocol.at(60)[0] = 130; m_PointsToTransformGroundTruthProtocol.at(60)[1] = 15; m_PointsToTransformGroundTruthProtocol.at(60)[2] = -67; // = 5mm distance to clipping plate + 2mm to base + 60mm depth m_PointsToTransformGroundTruthProtocol.at(80)[0] = 130; m_PointsToTransformGroundTruthProtocol.at(80)[1] = 15; m_PointsToTransformGroundTruthProtocol.at(80)[2] = -87; // = 5mm distance to clipping plate + 2mm to base + 80mm depth m_PointsToTransformGroundTruthProtocol.at(100)[0] = 130; m_PointsToTransformGroundTruthProtocol.at(100)[1] = 15; m_PointsToTransformGroundTruthProtocol.at(100)[2] = -107; // = 5mm distance to clipping plate + 2mm to base + 100mm depth break; // case 2 is equal to fiducial marker configuration C (20mm distance) case 2: m_PointsToTransformGroundTruthProtocol.at(0)[0] = 135; // = 20 + 15mm to end of clipping plate + 100 mm to middle axis of measurement plate m_PointsToTransformGroundTruthProtocol.at(0)[1] = 20; m_PointsToTransformGroundTruthProtocol.at(0)[2] = -9; // = 7mm distance to clipping plate + 2mm to base m_PointsToTransformGroundTruthProtocol.at(20)[0] = 135; m_PointsToTransformGroundTruthProtocol.at(20)[1] = 20; m_PointsToTransformGroundTruthProtocol.at(20)[2] = -29; // = 7mm distance to clipping plate + 2mm to base + 20mm depth m_PointsToTransformGroundTruthProtocol.at(40)[0] = 135; m_PointsToTransformGroundTruthProtocol.at(40)[1] = 20; m_PointsToTransformGroundTruthProtocol.at(40)[2] = -49; // = 7mm distance to clipping plate + 2mm to base + 40mm depth m_PointsToTransformGroundTruthProtocol.at(60)[0] = 135; m_PointsToTransformGroundTruthProtocol.at(60)[1] = 20; m_PointsToTransformGroundTruthProtocol.at(60)[2] = -69; // = 7mm distance to clipping plate + 2mm to base + 60mm depth m_PointsToTransformGroundTruthProtocol.at(80)[0] = 135; m_PointsToTransformGroundTruthProtocol.at(80)[1] = 20; m_PointsToTransformGroundTruthProtocol.at(80)[2] = -89; // = 7mm distance to clipping plate + 2mm to base + 80mm depth m_PointsToTransformGroundTruthProtocol.at(100)[0] = 135; m_PointsToTransformGroundTruthProtocol.at(100)[1] = 20; m_PointsToTransformGroundTruthProtocol.at(100)[2] = -109; // = 7mm distance to clipping plate + 2mm to base + 100mm depth break; } } void QmitkUSNavigationStepCtUsRegistration::TransformPointsGroundTruthProtocol() { if (m_GroundTruthProtocolTransformedPoints.find(0) == m_GroundTruthProtocolTransformedPoints.end()) { mitk::PointSet::Pointer pointSet = mitk::PointSet::New(); pointSet->InsertPoint(m_TransformMarkerCSToFloatingImageCS->TransformPoint(m_PointsToTransformGroundTruthProtocol.at(0))); m_GroundTruthProtocolTransformedPoints.insert(std::pair(0, pointSet)); } else { m_GroundTruthProtocolTransformedPoints.at(0)->InsertPoint(m_TransformMarkerCSToFloatingImageCS->TransformPoint(m_PointsToTransformGroundTruthProtocol.at(0))); } if (m_GroundTruthProtocolTransformedPoints.find(20) == m_GroundTruthProtocolTransformedPoints.end()) { mitk::PointSet::Pointer pointSet = mitk::PointSet::New(); pointSet->InsertPoint(m_TransformMarkerCSToFloatingImageCS->TransformPoint(m_PointsToTransformGroundTruthProtocol.at(20))); m_GroundTruthProtocolTransformedPoints.insert(std::pair(20, pointSet)); } else { m_GroundTruthProtocolTransformedPoints.at(20)->InsertPoint(m_TransformMarkerCSToFloatingImageCS->TransformPoint(m_PointsToTransformGroundTruthProtocol.at(20))); } if (m_GroundTruthProtocolTransformedPoints.find(40) == m_GroundTruthProtocolTransformedPoints.end()) { mitk::PointSet::Pointer pointSet = mitk::PointSet::New(); pointSet->InsertPoint(m_TransformMarkerCSToFloatingImageCS->TransformPoint(m_PointsToTransformGroundTruthProtocol.at(40))); m_GroundTruthProtocolTransformedPoints.insert(std::pair(40, pointSet)); } else { m_GroundTruthProtocolTransformedPoints.at(40)->InsertPoint(m_TransformMarkerCSToFloatingImageCS->TransformPoint(m_PointsToTransformGroundTruthProtocol.at(40))); } if (m_GroundTruthProtocolTransformedPoints.find(60) == m_GroundTruthProtocolTransformedPoints.end()) { mitk::PointSet::Pointer pointSet = mitk::PointSet::New(); pointSet->InsertPoint(m_TransformMarkerCSToFloatingImageCS->TransformPoint(m_PointsToTransformGroundTruthProtocol.at(60))); m_GroundTruthProtocolTransformedPoints.insert(std::pair(60, pointSet)); } else { m_GroundTruthProtocolTransformedPoints.at(60)->InsertPoint(m_TransformMarkerCSToFloatingImageCS->TransformPoint(m_PointsToTransformGroundTruthProtocol.at(60))); } if (m_GroundTruthProtocolTransformedPoints.find(80) == m_GroundTruthProtocolTransformedPoints.end()) { mitk::PointSet::Pointer pointSet = mitk::PointSet::New(); pointSet->InsertPoint(m_TransformMarkerCSToFloatingImageCS->TransformPoint(m_PointsToTransformGroundTruthProtocol.at(80))); m_GroundTruthProtocolTransformedPoints.insert(std::pair(80, pointSet)); } else { m_GroundTruthProtocolTransformedPoints.at(80)->InsertPoint(m_TransformMarkerCSToFloatingImageCS->TransformPoint(m_PointsToTransformGroundTruthProtocol.at(80))); } if (m_GroundTruthProtocolTransformedPoints.find(100) == m_GroundTruthProtocolTransformedPoints.end()) { mitk::PointSet::Pointer pointSet = mitk::PointSet::New(); pointSet->InsertPoint(m_TransformMarkerCSToFloatingImageCS->TransformPoint(m_PointsToTransformGroundTruthProtocol.at(100))); m_GroundTruthProtocolTransformedPoints.insert(std::pair(100, pointSet)); } else { m_GroundTruthProtocolTransformedPoints.at(100)->InsertPoint(m_TransformMarkerCSToFloatingImageCS->TransformPoint(m_PointsToTransformGroundTruthProtocol.at(100))); } } void QmitkUSNavigationStepCtUsRegistration::AddTransformedPointsToDataStorage() { if (m_GroundTruthProtocolTransformedPoints.find(0) == m_GroundTruthProtocolTransformedPoints.end() || m_GroundTruthProtocolTransformedPoints.find(20) == m_GroundTruthProtocolTransformedPoints.end() || m_GroundTruthProtocolTransformedPoints.find(40) == m_GroundTruthProtocolTransformedPoints.end() || m_GroundTruthProtocolTransformedPoints.find(60) == m_GroundTruthProtocolTransformedPoints.end() || m_GroundTruthProtocolTransformedPoints.find(80) == m_GroundTruthProtocolTransformedPoints.end() || m_GroundTruthProtocolTransformedPoints.find(100) == m_GroundTruthProtocolTransformedPoints.end()) { QMessageBox msgBox; msgBox.setText("Cannot add transformed Points to DataStorage because they do not exist.\ Stopping evaluation the protocol."); msgBox.exec(); return; } std::string nameNode0mm = "GroundTruthProt-Depth0mm"; std::string nameNode20mm = "GroundTruthProt-Depth20mm"; std::string nameNode40mm = "GroundTruthProt-Depth40mm"; std::string nameNode60mm = "GroundTruthProt-Depth60mm"; std::string nameNode80mm = "GroundTruthProt-Depth80mm"; std::string nameNode100mm = "GroundTruthProt-Depth100mm"; //Add transformed points of depth 0mm to the data storage mitk::DataNode::Pointer node0mm = this->GetDataStorage()->GetNamedNode(nameNode0mm); if (node0mm.IsNull()) { node0mm = mitk::DataNode::New(); node0mm->SetName(nameNode0mm); node0mm->SetData(m_GroundTruthProtocolTransformedPoints.at(0)); this->GetDataStorage()->Add(node0mm); } else { node0mm->SetData(m_GroundTruthProtocolTransformedPoints.at(0)); this->GetDataStorage()->Modified(); } if(ui->protocolEvaluationTypeComboBox->currentText().compare("PLANE") == 0 ) { //Add transformed points of depth 20mm to the data storage mitk::DataNode::Pointer node20mm = this->GetDataStorage()->GetNamedNode(nameNode20mm); if (node20mm.IsNull()) { node20mm = mitk::DataNode::New(); node20mm->SetName(nameNode20mm); node20mm->SetData(m_GroundTruthProtocolTransformedPoints.at(20)); this->GetDataStorage()->Add(node20mm); } else { node20mm->SetData(m_GroundTruthProtocolTransformedPoints.at(20)); this->GetDataStorage()->Modified(); } //Add transformed points of depth 40mm to the data storage mitk::DataNode::Pointer node40mm = this->GetDataStorage()->GetNamedNode(nameNode40mm); if (node40mm.IsNull()) { node40mm = mitk::DataNode::New(); node40mm->SetName(nameNode40mm); node40mm->SetData(m_GroundTruthProtocolTransformedPoints.at(40)); this->GetDataStorage()->Add(node40mm); } else { node40mm->SetData(m_GroundTruthProtocolTransformedPoints.at(40)); this->GetDataStorage()->Modified(); } //Add transformed points of depth 60mm to the data storage mitk::DataNode::Pointer node60mm = this->GetDataStorage()->GetNamedNode(nameNode60mm); if (node60mm.IsNull()) { node60mm = mitk::DataNode::New(); node60mm->SetName(nameNode60mm); node60mm->SetData(m_GroundTruthProtocolTransformedPoints.at(60)); this->GetDataStorage()->Add(node60mm); } else { node60mm->SetData(m_GroundTruthProtocolTransformedPoints.at(60)); this->GetDataStorage()->Modified(); } //Add transformed points of depth 80mm to the data storage mitk::DataNode::Pointer node80mm = this->GetDataStorage()->GetNamedNode(nameNode80mm); if (node80mm.IsNull()) { node80mm = mitk::DataNode::New(); node80mm->SetName(nameNode80mm); node80mm->SetData(m_GroundTruthProtocolTransformedPoints.at(80)); this->GetDataStorage()->Add(node80mm); } else { node80mm->SetData(m_GroundTruthProtocolTransformedPoints.at(80)); this->GetDataStorage()->Modified(); } //Add transformed points of depth 100mm to the data storage mitk::DataNode::Pointer node100mm = this->GetDataStorage()->GetNamedNode(nameNode100mm); if (node100mm.IsNull()) { node100mm = mitk::DataNode::New(); node100mm->SetName(nameNode100mm); node100mm->SetData(m_GroundTruthProtocolTransformedPoints.at(100)); this->GetDataStorage()->Add(node100mm); } else { node100mm->SetData(m_GroundTruthProtocolTransformedPoints.at(100)); this->GetDataStorage()->Modified(); } } //Do a global reinit mitk::RenderingManager::GetInstance()->InitializeViewsByBoundingObjects(this->GetDataStorage()); } double QmitkUSNavigationStepCtUsRegistration::CalculateMeanFRE() { double meanFRE = 0.0; for (int counter = 0; counter < m_GroundTruthProtocolFRE.size(); ++counter) { meanFRE += m_GroundTruthProtocolFRE[counter]; } return meanFRE / m_GroundTruthProtocolFRE.size(); } double QmitkUSNavigationStepCtUsRegistration::CalculateStandardDeviationOfFRE(double meanFRE) { double variance = 0.0; for (int counter = 0; counter < m_GroundTruthProtocolFRE.size(); ++counter) { variance += ((meanFRE - m_GroundTruthProtocolFRE[counter]) * (meanFRE - m_GroundTruthProtocolFRE[counter])); } variance /= m_GroundTruthProtocolFRE.size(); // calculate the empirical variance (n) and not the sampling variance (n-1) return sqrt(variance); } void QmitkUSNavigationStepCtUsRegistration::CalculateGroundTruthProtocolTRE() { if (m_GroundTruthProtocolTransformedPoints.find(0) == m_GroundTruthProtocolTransformedPoints.end() || m_GroundTruthProtocolTransformedPoints.find(20) == m_GroundTruthProtocolTransformedPoints.end() || m_GroundTruthProtocolTransformedPoints.find(40) == m_GroundTruthProtocolTransformedPoints.end() || m_GroundTruthProtocolTransformedPoints.find(60) == m_GroundTruthProtocolTransformedPoints.end() || m_GroundTruthProtocolTransformedPoints.find(80) == m_GroundTruthProtocolTransformedPoints.end() || m_GroundTruthProtocolTransformedPoints.find(100) == m_GroundTruthProtocolTransformedPoints.end()) { QMessageBox msgBox; msgBox.setText("Cannot calculate TRE of Ground-Truth-Protocol because points were not transformed."); msgBox.exec(); return; } // clear the std::map containing possibly data of earlier TRE calculations m_GroundTruthProtocolTRE.clear(); // loop through all existing point sets containing the transformed points for (int counter = 0; m_GroundTruthProtocolTransformedPoints.find(counter) != m_GroundTruthProtocolTransformedPoints.end(); counter += 20) { //calculate the middle point of the point set mitk::PointSet::Pointer pointSet = m_GroundTruthProtocolTransformedPoints.at(counter); mitk::Point3D middlePoint; middlePoint[0] = 0.0; middlePoint[1] = 0.0; middlePoint[2] = 0.0; for (int position = 0; position < pointSet->GetSize(); ++position) { middlePoint[0] += pointSet->GetPoint(position)[0]; middlePoint[1] += pointSet->GetPoint(position)[1]; middlePoint[2] += pointSet->GetPoint(position)[2]; } middlePoint[0] /= pointSet->GetSize(); middlePoint[1] /= pointSet->GetSize(); middlePoint[2] /= pointSet->GetSize(); MITK_INFO << "Calculated MiddlePoint: " << middlePoint; //sum up the euclidean distances between the middle point and each transformed point double meanDistance = 0.0; for (int position = 0; position < pointSet->GetSize(); ++position) { meanDistance += middlePoint.SquaredEuclideanDistanceTo(pointSet->GetPoint(position)); MITK_INFO << "SquaredEuclideanDistance: " << middlePoint.SquaredEuclideanDistanceTo(pointSet->GetPoint(position)); } meanDistance /= pointSet->GetSize(); // this can be interpreted as empirical variance // the root of the empirical variance can be interpreted as the protocols registration TRE m_GroundTruthProtocolTRE.insert(std::pair(counter, sqrt(meanDistance))); MITK_INFO << "Ground-Truth-Protocol TRE: " << sqrt(meanDistance); } } void QmitkUSNavigationStepCtUsRegistration::EliminateTooSmallLabeledObjects( ImageType::Pointer binaryImage) { BinaryImageToShapeLabelMapFilterType::OutputImageType::Pointer labelMap = m_BinaryImageToShapeLabelMapFilter->GetOutput(); double voxelVolume = this->GetVoxelVolume(); double fiducialVolume; unsigned int numberOfPixels; if (ui->fiducialDiameter3mmRadioButton->isChecked()) { fiducialVolume = this->GetFiducialVolume(1.5); numberOfPixels = ceil(fiducialVolume / voxelVolume); } else { fiducialVolume = this->GetFiducialVolume(2.5); numberOfPixels = ceil(fiducialVolume / voxelVolume); } MITK_INFO << "Voxel Volume = " << voxelVolume << "; Fiducial Volume = " << fiducialVolume; MITK_INFO << "Number of pixels = " << numberOfPixels; labelMap = m_BinaryImageToShapeLabelMapFilter->GetOutput(); // The output of this filter is an itk::LabelMap, which contains itk::LabelObject's MITK_INFO << "There are " << labelMap->GetNumberOfLabelObjects() << " objects."; // Loop over each region for (int i = labelMap->GetNumberOfLabelObjects() - 1; i >= 0; --i) { // Get the ith region BinaryImageToShapeLabelMapFilterType::OutputImageType::LabelObjectType* labelObject = labelMap->GetNthLabelObject(i); MITK_INFO << "Object " << i << " contains " << labelObject->Size() << " pixel"; //TODO: Threshold-Wert evtl. experimentell besser abstimmen, // um zu verhindern, dass durch Threshold wahre Fiducial-Kandidaten elimiert werden. if (labelObject->Size() < numberOfPixels * 0.8) { for (unsigned int pixelId = 0; pixelId < labelObject->Size(); pixelId++) { binaryImage->SetPixel(labelObject->GetIndex(pixelId), 0); } labelMap->RemoveLabelObject(labelObject); } } } bool QmitkUSNavigationStepCtUsRegistration::EliminateFiducialCandidatesByEuclideanDistances() { if (m_CentroidsOfFiducialCandidates.size() < NUMBER_FIDUCIALS_NEEDED) { return false; } for (int counter = 0; counter < m_CentroidsOfFiducialCandidates.size(); ++counter) { int amountOfAcceptedFiducials = 0; mitk::Point3D fiducialCentroid(m_CentroidsOfFiducialCandidates.at(counter)); //Loop through all fiducial candidates and calculate the distance between the chosen fiducial // candidate and the other candidates. For each candidate with a right distance between // Configuration A: 7.93mm and 31.0mm (10 mm distance between fiducial centers) or // Configuration B: 11.895mm and 45.0mm (15 mm distance between fiducial centers) or // Configuration C: 15.86mm and 59.0mm (20 mm distance between fiducial centers) // // increase the amountOfAcceptedFiducials. for (int position = 0; position < m_CentroidsOfFiducialCandidates.size(); ++position) { if (position == counter) { continue; } mitk::Point3D otherCentroid(m_CentroidsOfFiducialCandidates.at(position)); double distance = fiducialCentroid.EuclideanDistanceTo(otherCentroid); switch (ui->fiducialMarkerConfigurationComboBox->currentIndex()) { // case 0 is equal to fiducial marker configuration A (10mm distance) case 0: if (distance > 7.93 && distance < 31.0) { ++amountOfAcceptedFiducials; } break; // case 1 is equal to fiducial marker configuration B (15mm distance) case 1: if (distance > 11.895 && distance < 45.0) { ++amountOfAcceptedFiducials; } break; // case 2 is equal to fiducial marker configuration C (20mm distance) case 2: if (distance > 15.86 && distance < 59.0) { ++amountOfAcceptedFiducials; } break; } } //The amountOfAcceptedFiducials must be at least 7. Otherwise delete the fiducial candidate // from the list of candidates. if (amountOfAcceptedFiducials < NUMBER_FIDUCIALS_NEEDED - 1) { MITK_INFO << "Deleting fiducial candidate at position: " << m_CentroidsOfFiducialCandidates.at(counter); m_CentroidsOfFiducialCandidates.erase(m_CentroidsOfFiducialCandidates.begin() + counter); if (m_CentroidsOfFiducialCandidates.size() < NUMBER_FIDUCIALS_NEEDED ) { return false; } counter = -1; } } //Classify the rested fiducial candidates by its characteristic Euclidean distances // between the canidates and remove all candidates with a false distance configuration: this->ClassifyFiducialCandidates(); return true; } void QmitkUSNavigationStepCtUsRegistration::ClassifyFiducialCandidates() { MITK_INFO << "ClassifyFiducialCandidates()"; std::vector fiducialCandidatesToBeRemoved; std::vector> distanceVectorsFiducials; this->CalculateDistancesBetweenFiducials(distanceVectorsFiducials); for (int counter = 0; counter < distanceVectorsFiducials.size(); ++counter) { int distanceA = 0; // => 10,00mm distance int distanceB = 0; // => 14,14mm distance int distanceC = 0; // => 17,32mm distance int distanceD = 0; // => 22,36mm distance int distanceE = 0; // => 24,49mm distance int distanceF = 0; // => 28,28mm distance std::vector &distances = distanceVectorsFiducials.at(counter); for (int number = 0; number < distances.size(); ++number) { double &distance = distances.at(number); switch (ui->fiducialMarkerConfigurationComboBox->currentIndex()) { // case 0 is equal to fiducial marker configuration A (10mm distance) case 0: if (distance > 7.93 && distance <= 12.07) { ++distanceA; } else if (distance > 12.07 && distance <= 15.73) { ++distanceB; } else if (distance > 15.73 && distance <= 19.84) { ++distanceC; } else if (distance > 19.84 && distance <= 23.425) { ++distanceD; } else if (distance > 23.425 && distance <= 26.385) { ++distanceE; } else if (distance > 26.385 && distance <= 31.00) { ++distanceF; } break; // case 1 is equal to fiducial marker configuration B (15mm distance) case 1: if (distance > 11.895 && distance <= 18.105) { ++distanceA; } else if (distance > 18.105 && distance <= 23.595) { ++distanceB; } else if (distance > 23.595 && distance <= 29.76) { ++distanceC; } else if (distance > 29.76 && distance <= 35.1375) { ++distanceD; + if (distance > 33.54) + { + ++distanceE; + } } else if (distance > 35.1375 && distance <= 39.5775) { ++distanceE; + if (distance < 36.735) + { + ++distanceD; + } } else if (distance > 39.5775 && distance <= 45.00) { ++distanceF; } break; // case 2 is equal to fiducial marker configuration C (20mm distance) case 2: if (distance > 15.86 && distance <= 24.14) { ++distanceA; } else if (distance > 24.14 && distance <= 31.46) { ++distanceB; } else if (distance > 31.46 && distance <= 39.68) { ++distanceC; } else if (distance > 39.68 && distance <= 46.85) { ++distanceD; } else if (distance > 46.85 && distance <= 52.77) { ++distanceE; } else if (distance > 52.77 && distance <= 59.00) { ++distanceF; } break; } }// End for-loop distances-vector //Now, having looped through all distances of one fiducial candidate, check // if the combination of different distances is known. The >= is due to the // possible occurrence of other fiducial candidates that have an distance equal to // one of the distances A - E. However, false fiducial candidates outside // the fiducial marker does not have the right distance configuration: if ((distanceA >= 2 && distanceD >= 2 && distanceE >= 2 && distanceF >= 1 || distanceA >= 1 && distanceB >= 2 && distanceC >= 1 && distanceD >= 2 && distanceE >= 1 || distanceB >= 2 && distanceD >= 4 && distanceF >= 1 || distanceA >= 1 && distanceB >= 1 && distanceD >= 3 && distanceE >= 1 && distanceF >= 1) == false) { MITK_INFO << "Detected fiducial candidate with unknown distance configuration."; fiducialCandidatesToBeRemoved.push_back(counter); } } for (int count = fiducialCandidatesToBeRemoved.size() - 1; count >= 0; --count) { MITK_INFO << "Removing fiducial candidate " << fiducialCandidatesToBeRemoved.at(count); m_CentroidsOfFiducialCandidates.erase(m_CentroidsOfFiducialCandidates.begin() + fiducialCandidatesToBeRemoved.at(count)); } } void QmitkUSNavigationStepCtUsRegistration::GetCentroidsOfLabeledObjects() { MITK_INFO << "GetCentroidsOfLabeledObjects()"; BinaryImageToShapeLabelMapFilterType::OutputImageType::Pointer labelMap = m_BinaryImageToShapeLabelMapFilter->GetOutput(); for (int i = labelMap->GetNumberOfLabelObjects() - 1; i >= 0; --i) { // Get the ith region BinaryImageToShapeLabelMapFilterType::OutputImageType::LabelObjectType* labelObject = labelMap->GetNthLabelObject(i); MITK_INFO << "Object " << i << " contains " << labelObject->Size() << " pixel"; mitk::Vector3D centroid; centroid[0] = labelObject->GetCentroid()[0]; centroid[1] = labelObject->GetCentroid()[1]; centroid[2] = labelObject->GetCentroid()[2]; m_CentroidsOfFiducialCandidates.push_back(centroid); } //evtl. for later: itk::LabelMapOverlayImageFilter } void QmitkUSNavigationStepCtUsRegistration::CalculatePCA() { //Step 1: Construct data matrix int columnSize = m_CentroidsOfFiducialCandidates.size(); if (columnSize == 0) { MITK_INFO << "Cannot calculate PCA. There are no fiducial candidates."; return; } vnl_matrix pointSetMatrix(3, columnSize, 0.0); for (int counter = 0; counter < columnSize; ++counter) { pointSetMatrix[0][counter] = m_CentroidsOfFiducialCandidates.at(counter)[0]; pointSetMatrix[1][counter] = m_CentroidsOfFiducialCandidates.at(counter)[1]; pointSetMatrix[2][counter] = m_CentroidsOfFiducialCandidates.at(counter)[2]; } //Step 2: Remove average for each row (Mittelwertbefreiung) for (int counter = 0; counter < columnSize; ++counter) { m_MeanCentroidFiducialCandidates += mitk::Vector3D(pointSetMatrix.get_column(counter)); } //TODO: für später überprüfen, ob Division durch integer nicht zu Rechenproblemen führt. m_MeanCentroidFiducialCandidates /= columnSize; for (int counter = 0; counter < columnSize; ++counter) { pointSetMatrix.get_column(counter) -= m_MeanCentroidFiducialCandidates; } //Step 3: Compute covariance matrix vnl_matrix covarianceMatrix = (1.0 / (columnSize - 1.0)) * pointSetMatrix * pointSetMatrix.transpose(); //Step 4: Singular value composition vnl_svd svd(covarianceMatrix); //Storing results: for (int counter = 0; counter < 3; ++counter) { mitk::Vector3D eigenVector = svd.U().get_column(counter); double eigenValue = sqrt(svd.W(counter)); m_EigenVectorsFiducialCandidates[eigenValue] = eigenVector; m_EigenValuesFiducialCandidates.push_back(eigenValue); } std::sort( m_EigenValuesFiducialCandidates.begin(), m_EigenValuesFiducialCandidates.end() ); mitk::DataNode::Pointer axis1Node = mitk::DataNode::New(); axis1Node->SetName("Eigenvector 1"); mitk::PointSet::Pointer axis1 = mitk::PointSet::New(); axis1->InsertPoint(0, m_MeanCentroidFiducialCandidates); axis1->InsertPoint(1, (m_MeanCentroidFiducialCandidates + m_EigenVectorsFiducialCandidates.at(m_EigenValuesFiducialCandidates.at(2))*m_EigenValuesFiducialCandidates.at(2))); axis1Node->SetData(axis1); axis1Node->SetBoolProperty("show contour", true); axis1Node->SetColor(1, 0, 0); this->GetDataStorage()->Add(axis1Node); mitk::DataNode::Pointer axis2Node = mitk::DataNode::New(); axis2Node->SetName("Eigenvector 2"); mitk::PointSet::Pointer axis2 = mitk::PointSet::New(); axis2->InsertPoint(0, m_MeanCentroidFiducialCandidates); axis2->InsertPoint(1, (m_MeanCentroidFiducialCandidates + m_EigenVectorsFiducialCandidates.at(m_EigenValuesFiducialCandidates.at(1))*m_EigenValuesFiducialCandidates.at(1))); axis2Node->SetData(axis2); axis2Node->SetBoolProperty("show contour", true); axis2Node->SetColor(2, 0, 0); this->GetDataStorage()->Add(axis2Node); mitk::DataNode::Pointer axis3Node = mitk::DataNode::New(); axis3Node->SetName("Eigenvector 3"); mitk::PointSet::Pointer axis3 = mitk::PointSet::New(); axis3->InsertPoint(0, m_MeanCentroidFiducialCandidates); axis3->InsertPoint(1, (m_MeanCentroidFiducialCandidates + m_EigenVectorsFiducialCandidates.at(m_EigenValuesFiducialCandidates.at(0))*m_EigenValuesFiducialCandidates.at(0))); axis3Node->SetData(axis3); axis3Node->SetBoolProperty("show contour", true); axis3Node->SetColor(3, 0, 0); this->GetDataStorage()->Add(axis3Node); MITK_INFO << "Mean: " << m_MeanCentroidFiducialCandidates; MITK_INFO << "Eigenvektor 1: " << m_EigenVectorsFiducialCandidates.at(m_EigenValuesFiducialCandidates.at(2)); MITK_INFO << "Eigenvektor 2: " << m_EigenVectorsFiducialCandidates.at(m_EigenValuesFiducialCandidates.at(1)); MITK_INFO << "Eigenvektor 3: " << m_EigenVectorsFiducialCandidates.at(m_EigenValuesFiducialCandidates.at(0)); MITK_INFO << "Eigenwert 1: " << m_EigenValuesFiducialCandidates.at(2); MITK_INFO << "Eigenwert 2: " << m_EigenValuesFiducialCandidates.at(1); MITK_INFO << "Eigenwert 3: " << m_EigenValuesFiducialCandidates.at(0); } void QmitkUSNavigationStepCtUsRegistration::NumerateFiducialMarks() { MITK_INFO << "NumerateFiducialMarks()"; bool successFiducialNo1; bool successFiducialNo4; bool successFiducialNo2And3; bool successFiducialNo5; bool successFiducialNo8; bool successFiducialNo6; bool successFiducialNo7; std::vector> distanceVectorsFiducials; this->CalculateDistancesBetweenFiducials(distanceVectorsFiducials); successFiducialNo1 = this->FindFiducialNo1(distanceVectorsFiducials); successFiducialNo4 = this->FindFiducialNo4(distanceVectorsFiducials); successFiducialNo2And3 = this->FindFiducialNo2And3(); successFiducialNo5 = this->FindFiducialNo5(); successFiducialNo8 = this->FindFiducialNo8(); successFiducialNo6 = this->FindFiducialNo6(); successFiducialNo7 = this->FindFiducialNo7(); if (!successFiducialNo1 || !successFiducialNo4 || !successFiducialNo2And3 || !successFiducialNo5 || !successFiducialNo8 || !successFiducialNo6 || !successFiducialNo7) { QMessageBox msgBox; msgBox.setText("Cannot numerate/localize all fiducials successfully."); msgBox.exec(); return; } if( ui->optimizeFiducialPositionsCheckBox->isChecked()) { this->OptimizeFiducialPositions(); } if (m_MarkerFloatingImageCoordinateSystemPointSet.IsNull()) { m_MarkerFloatingImageCoordinateSystemPointSet = mitk::PointSet::New(); } else if (m_MarkerFloatingImageCoordinateSystemPointSet->GetSize() != 0) { m_MarkerFloatingImageCoordinateSystemPointSet->Clear(); } for (int counter = 1; counter <= m_FiducialMarkerCentroids.size(); ++counter) { m_MarkerFloatingImageCoordinateSystemPointSet->InsertPoint(counter - 1, m_FiducialMarkerCentroids.at(counter)); } if( !m_PerformingGroundTruthProtocolEvaluation ) { mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData(m_MarkerFloatingImageCoordinateSystemPointSet); node->SetName("MarkerFloatingImageCSPointSet"); //node->SetFloatProperty("pointsize", 5.0); this->GetDataStorage()->Add(node); } } void QmitkUSNavigationStepCtUsRegistration::ShowGroundTruthMarkerEdges() { mitk::Point3D edge1; mitk::Point3D edge2; mitk::Point3D edge3; mitk::Point3D edge4; mitk::Point3D edge5; mitk::Point3D edge6; mitk::Point3D edge7; mitk::Point3D edge8; edge1[0] = -5; edge1[1] = -5; edge1[2] = -7.5; edge2[0] = -5; edge2[1] = 25; edge2[2] = -7.5; edge3[0] = 25; edge3[1] = 25; edge3[2] = -7.5; edge4[0] = 25; edge4[1] = -5; edge4[2] = -7.5; edge5[0] = -5; edge5[1] = -5; edge5[2] = 15; edge6[0] = -5; edge6[1] = 25; edge6[2] = 15; edge7[0] = 25; edge7[1] = 25; edge7[2] = 15; edge8[0] = 25; edge8[1] = -5; edge8[2] = 15; mitk::PointSet::Pointer edgePointSet = mitk::PointSet::New(); edgePointSet->InsertPoint(0, edge1); edgePointSet->InsertPoint(1, edge2); edgePointSet->InsertPoint(2, edge3); edgePointSet->InsertPoint(3, edge4); edgePointSet->InsertPoint(4, edge1); edgePointSet->InsertPoint(5, edge5); edgePointSet->InsertPoint(6, edge6); edgePointSet->InsertPoint(7, edge7); edgePointSet->InsertPoint(8, edge8); edgePointSet->InsertPoint(9, edge5); edgePointSet->InsertPoint(10, edge6); edgePointSet->InsertPoint(11, edge2); edgePointSet->InsertPoint(12, edge6); edgePointSet->InsertPoint(13, edge7); edgePointSet->InsertPoint(14, edge3); edgePointSet->InsertPoint(15, edge7); edgePointSet->InsertPoint(16, edge8); edgePointSet->InsertPoint(17, edge4); mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData(edgePointSet); node->SetName("Marker Edges PointSet"); node->SetBoolProperty("show contour", true); node->SetFloatProperty("pointsize", 0.25); this->GetDataStorage()->Add(node); mitk::PointSet* pointSet_orig = edgePointSet; mitk::PointSet::Pointer pointSet_moved = mitk::PointSet::New(); for (int i = 0; i < pointSet_orig->GetSize(); i++) { pointSet_moved->InsertPoint(m_TransformMarkerCSToFloatingImageCS->TransformPoint(pointSet_orig->GetPoint(i))); } pointSet_orig->Clear(); for (int i = 0; i < pointSet_moved->GetSize(); i++) pointSet_orig->InsertPoint(pointSet_moved->GetPoint(i)); //Do a global reinit mitk::RenderingManager::GetInstance()->InitializeViewsByBoundingObjects(this->GetDataStorage()); } void QmitkUSNavigationStepCtUsRegistration::CalculateDistancesBetweenFiducials(std::vector>& distanceVectorsFiducials) { std::vector distancesBetweenFiducials; for (int i = 0; i < m_CentroidsOfFiducialCandidates.size(); ++i) { distancesBetweenFiducials.clear(); mitk::Point3D fiducialCentroid(m_CentroidsOfFiducialCandidates.at(i)); for (int n = 0; n < m_CentroidsOfFiducialCandidates.size(); ++n) { mitk::Point3D otherCentroid(m_CentroidsOfFiducialCandidates.at(n)); distancesBetweenFiducials.push_back(fiducialCentroid.EuclideanDistanceTo(otherCentroid)); } //Sort the distances from low to big numbers std::sort(distancesBetweenFiducials.begin(), distancesBetweenFiducials.end()); //First entry of the distance vector must be 0, so erase it if (distancesBetweenFiducials.at(0) == 0.0) { distancesBetweenFiducials.erase(distancesBetweenFiducials.begin()); } //Add the distance vector to the collecting distances vector distanceVectorsFiducials.push_back(distancesBetweenFiducials); } for (int i = 0; i < distanceVectorsFiducials.size(); ++i) { MITK_INFO << "Vector " << i << ":"; for (int k = 0; k < distanceVectorsFiducials.at(i).size(); ++k) { MITK_INFO << distanceVectorsFiducials.at(i).at(k); } } } bool QmitkUSNavigationStepCtUsRegistration::FindFiducialNo1(std::vector>& distanceVectorsFiducials) { for (int i = 0; i < distanceVectorsFiducials.size(); ++i) { std::vector &distances = distanceVectorsFiducials.at(i); if (distances.size() < NUMBER_FIDUCIALS_NEEDED - 1 ) { MITK_WARN << "Cannot find fiducial 1, there aren't found enough fiducial candidates."; return false; } double characteristicDistanceAWithUpperMargin = this->GetCharacteristicDistanceAWithUpperMargin(); if (distances.at(0) <= characteristicDistanceAWithUpperMargin && distances.at(1) <= characteristicDistanceAWithUpperMargin) { MITK_INFO << "Found Fiducial 1 (PointSet number " << i << ")"; m_FiducialMarkerCentroids.insert( std::pair(1, m_CentroidsOfFiducialCandidates.at(i))); distanceVectorsFiducials.erase(distanceVectorsFiducials.begin() + i); m_CentroidsOfFiducialCandidates.erase(m_CentroidsOfFiducialCandidates.begin() + i); return true; } } return false; } bool QmitkUSNavigationStepCtUsRegistration::FindFiducialNo2And3() { if (m_FiducialMarkerCentroids.find(1) == m_FiducialMarkerCentroids.end() ) { MITK_WARN << "Cannot find fiducial No 2 and 3. Before must be found fiducial No 1."; return false; } mitk::Point3D fiducialNo1(m_FiducialMarkerCentroids.at(1)); mitk::Vector3D fiducialVectorA; mitk::Vector3D fiducialVectorB; mitk::Point3D fiducialPointA; mitk::Point3D fiducialPointB; bool foundFiducialA = false; bool foundFiducialB = false; mitk::Vector3D vectorFiducial1ToFiducialA; mitk::Vector3D vectorFiducial1ToFiducialB; for (int i = 0; i < m_CentroidsOfFiducialCandidates.size(); ++i) { mitk::Point3D fiducialCentroid(m_CentroidsOfFiducialCandidates.at(i)); double distance = fiducialNo1.EuclideanDistanceTo(fiducialCentroid); if (distance <= this->GetCharacteristicDistanceAWithUpperMargin()) { fiducialVectorA = m_CentroidsOfFiducialCandidates.at(i); fiducialPointA = fiducialCentroid; m_CentroidsOfFiducialCandidates.erase(m_CentroidsOfFiducialCandidates.begin() + i); foundFiducialA = true; break; } } for (int i = 0; i < m_CentroidsOfFiducialCandidates.size(); ++i) { mitk::Point3D fiducialCentroid(m_CentroidsOfFiducialCandidates.at(i)); double distance = fiducialNo1.EuclideanDistanceTo(fiducialCentroid); if (distance <= this->GetCharacteristicDistanceAWithUpperMargin()) { fiducialVectorB = m_CentroidsOfFiducialCandidates.at(i); fiducialPointB = fiducialCentroid; m_CentroidsOfFiducialCandidates.erase(m_CentroidsOfFiducialCandidates.begin() + i); foundFiducialB = true; break; } } if (!foundFiducialA || !foundFiducialB) { MITK_WARN << "Cannot identify fiducial candidates 2 and 3"; return false; } else if (m_CentroidsOfFiducialCandidates.size() == 0) { MITK_WARN << "Too less fiducials detected. Cannot identify fiducial candidates 2 and 3"; return false; } vectorFiducial1ToFiducialA = fiducialVectorA - m_FiducialMarkerCentroids.at(1); vectorFiducial1ToFiducialB = fiducialVectorB - m_FiducialMarkerCentroids.at(1); vnl_vector crossProductVnl = vnl_cross_3d(vectorFiducial1ToFiducialA.GetVnlVector(), vectorFiducial1ToFiducialB.GetVnlVector()); mitk::Vector3D crossProduct; crossProduct.SetVnlVector(crossProductVnl); mitk::Vector3D vectorFiducial1ToRandomLeftFiducial = m_CentroidsOfFiducialCandidates.at(0) - m_FiducialMarkerCentroids.at(1); double scalarProduct = (crossProduct * vectorFiducial1ToRandomLeftFiducial) / (crossProduct.GetNorm() * vectorFiducial1ToRandomLeftFiducial.GetNorm()); double alpha = acos(scalarProduct) * 57.29578; //Transform into degree MITK_INFO << "Scalar Product = " << alpha; if (alpha <= 90) { m_FiducialMarkerCentroids[3] = fiducialVectorA; m_FiducialMarkerCentroids[2] = fiducialVectorB; } else { m_FiducialMarkerCentroids[2] = fiducialVectorA; m_FiducialMarkerCentroids[3] = fiducialVectorB; } MITK_INFO << "Found Fiducial 2, PointSet: " << m_FiducialMarkerCentroids.at(2); MITK_INFO << "Found Fiducial 3, PointSet: " << m_FiducialMarkerCentroids.at(3); return true; } bool QmitkUSNavigationStepCtUsRegistration::FindFiducialNo4(std::vector>& distanceVectorsFiducials) { double characteristicDistanceBWithLowerMargin = this->GetCharacteristicDistanceBWithLowerMargin(); double characteristicDistanceBWithUpperMargin = this->GetCharacteristicDistanceBWithUpperMargin(); for (int i = 0; i < distanceVectorsFiducials.size(); ++i) { std::vector &distances = distanceVectorsFiducials.at(i); if (distances.size() < NUMBER_FIDUCIALS_NEEDED - 1) { MITK_WARN << "Cannot find fiducial 4, there aren't found enough fiducial candidates."; return false; } if (distances.at(0) > characteristicDistanceBWithLowerMargin && distances.at(0) <= characteristicDistanceBWithUpperMargin && distances.at(1) > characteristicDistanceBWithLowerMargin && distances.at(1) <= characteristicDistanceBWithUpperMargin) { MITK_INFO << "Found Fiducial 4 (PointSet number " << i << ")"; m_FiducialMarkerCentroids.insert(std::pair(4, m_CentroidsOfFiducialCandidates.at(i))); distanceVectorsFiducials.erase(distanceVectorsFiducials.begin() + i); m_CentroidsOfFiducialCandidates.erase(m_CentroidsOfFiducialCandidates.begin() + i); return true; } } return false; } bool QmitkUSNavigationStepCtUsRegistration::FindFiducialNo5() { if (m_FiducialMarkerCentroids.find(2) == m_FiducialMarkerCentroids.end()) { MITK_WARN << "To find fiducial No 5, fiducial No 2 has to be found before."; return false; } double characteristicDistanceBWithUpperMargin = this->GetCharacteristicDistanceBWithUpperMargin(); mitk::Point3D fiducialNo2(m_FiducialMarkerCentroids.at(2)); for (int counter = 0; counter < m_CentroidsOfFiducialCandidates.size(); ++counter) { mitk::Point3D fiducialCentroid(m_CentroidsOfFiducialCandidates.at(counter)); double distance = fiducialNo2.EuclideanDistanceTo(fiducialCentroid); if (distance <= characteristicDistanceBWithUpperMargin) { m_FiducialMarkerCentroids[5] = m_CentroidsOfFiducialCandidates.at(counter); m_CentroidsOfFiducialCandidates.erase(m_CentroidsOfFiducialCandidates.begin() + counter); MITK_INFO << "Found Fiducial No 5, PointSet: " << m_FiducialMarkerCentroids[5]; return true; } } MITK_WARN << "Cannot find fiducial No 5."; return false; } bool QmitkUSNavigationStepCtUsRegistration::FindFiducialNo6() { if (m_FiducialMarkerCentroids.find(5) == m_FiducialMarkerCentroids.end()) { MITK_WARN << "To find fiducial No 6, fiducial No 5 has to be found before."; return false; } double characteristicDistanceAWithUpperMargin = this->GetCharacteristicDistanceAWithUpperMargin(); mitk::Point3D fiducialNo5(m_FiducialMarkerCentroids.at(5)); for (int counter = 0; counter < m_CentroidsOfFiducialCandidates.size(); ++counter) { mitk::Point3D fiducialCentroid(m_CentroidsOfFiducialCandidates.at(counter)); double distance = fiducialNo5.EuclideanDistanceTo(fiducialCentroid); if (distance <= characteristicDistanceAWithUpperMargin) { m_FiducialMarkerCentroids[6] = m_CentroidsOfFiducialCandidates.at(counter); m_CentroidsOfFiducialCandidates.erase(m_CentroidsOfFiducialCandidates.begin() + counter); MITK_INFO << "Found Fiducial No 6, PointSet: " << m_FiducialMarkerCentroids[6]; return true; } } MITK_WARN << "Cannot find fiducial No 6."; return false; } bool QmitkUSNavigationStepCtUsRegistration::FindFiducialNo7() { if (m_FiducialMarkerCentroids.find(8) == m_FiducialMarkerCentroids.end()) { MITK_WARN << "To find fiducial No 7, fiducial No 8 has to be found before."; return false; } double characteristicDistanceAWithUpperMargin = this->GetCharacteristicDistanceAWithUpperMargin(); mitk::Point3D fiducialNo8(m_FiducialMarkerCentroids.at(8)); for (int counter = 0; counter < m_CentroidsOfFiducialCandidates.size(); ++counter) { mitk::Point3D fiducialCentroid(m_CentroidsOfFiducialCandidates.at(counter)); double distance = fiducialNo8.EuclideanDistanceTo(fiducialCentroid); if (distance <= characteristicDistanceAWithUpperMargin) { m_FiducialMarkerCentroids[7] = m_CentroidsOfFiducialCandidates.at(counter); m_CentroidsOfFiducialCandidates.erase(m_CentroidsOfFiducialCandidates.begin() + counter); MITK_INFO << "Found Fiducial No 7, PointSet: " << m_FiducialMarkerCentroids[7]; return true; } } MITK_WARN << "Cannot find fiducial No 7."; return false; } bool QmitkUSNavigationStepCtUsRegistration::FindFiducialNo8() { if (m_FiducialMarkerCentroids.find(3) == m_FiducialMarkerCentroids.end()) { MITK_WARN << "To find fiducial No 8, fiducial No 3 has to be found before."; return false; } double characteristicDistanceBWithUpperMargin = this->GetCharacteristicDistanceBWithUpperMargin(); mitk::Point3D fiducialNo3(m_FiducialMarkerCentroids.at(3)); for (int counter = 0; counter < m_CentroidsOfFiducialCandidates.size(); ++counter) { mitk::Point3D fiducialCentroid(m_CentroidsOfFiducialCandidates.at(counter)); double distance = fiducialNo3.EuclideanDistanceTo(fiducialCentroid); if (distance <= characteristicDistanceBWithUpperMargin) { m_FiducialMarkerCentroids[8] = m_CentroidsOfFiducialCandidates.at(counter); m_CentroidsOfFiducialCandidates.erase(m_CentroidsOfFiducialCandidates.begin() + counter); MITK_INFO << "Found Fiducial No 8, PointSet: " << m_FiducialMarkerCentroids[8]; return true; } } MITK_WARN << "Cannot find fiducial No 8."; return false; } void QmitkUSNavigationStepCtUsRegistration::OptimizeFiducialPositions() { MITK_INFO << "OptimizeFiducialPositions()"; //Initialization for planes 1 - 2 - 3 - 4 and 5 - 6 - 7 - 8 mitk::PlaneFit::Pointer planeFit1234 = mitk::PlaneFit::New(); mitk::PlaneFit::Pointer planeFit5678 = mitk::PlaneFit::New(); mitk::PointSet::Pointer pointSet1234 = mitk::PointSet::New(); mitk::PointSet::Pointer pointSet5678 = mitk::PointSet::New(); mitk::PlaneGeometry::Pointer planeGeometry1234 = mitk::PlaneGeometry::New(); mitk::PlaneGeometry::Pointer planeGeometry5678 = mitk::PlaneGeometry::New(); for (int counter = 1; counter <= 4; ++counter) { pointSet1234->InsertPoint(counter - 1, m_FiducialMarkerCentroids.at(counter)); } for (int counter = 5; counter <= 8; ++counter) { pointSet5678->InsertPoint(counter - 5, m_FiducialMarkerCentroids.at(counter)); } //Make planes parallel this->CreateParallelPlanes(planeFit1234, planeFit5678, pointSet1234, pointSet5678, planeGeometry1234, planeGeometry5678, true); this->MovePlanes(planeGeometry1234, planeGeometry5678, this->GetMinimalFiducialConfigurationDistance()); //Move the points into the parallel planes for (int counter = 1; counter <= 4; ++counter) { this->MovePoint(planeGeometry1234, counter); } for (int counter = 5; counter <= 8; ++counter) { this->MovePoint(planeGeometry5678, counter); } MITK_INFO << "NormalVector plane 1234 = " << planeGeometry1234->GetNormal(); MITK_INFO << "NormalVector plane 5678 = " << planeGeometry5678->GetNormal(); MITK_INFO << "Are parallel: " << planeGeometry1234->IsParallel(planeGeometry5678); //Optimize now the positions of Fiducials 1 - 2 - 5 and 4 - 7 - 8 //Initialization for planes 1 - 2 - 5 and 4 - 7 - 8 mitk::PlaneFit::Pointer planeFit125 = mitk::PlaneFit::New(); mitk::PlaneFit::Pointer planeFit478 = mitk::PlaneFit::New(); mitk::PointSet::Pointer pointSet125 = mitk::PointSet::New(); mitk::PointSet::Pointer pointSet478 = mitk::PointSet::New(); mitk::PlaneGeometry::Pointer planeGeometry125 = mitk::PlaneGeometry::New(); mitk::PlaneGeometry::Pointer planeGeometry478 = mitk::PlaneGeometry::New(); pointSet125->InsertPoint(0, m_FiducialMarkerCentroids.at(1)); pointSet125->InsertPoint(1, m_FiducialMarkerCentroids.at(2)); pointSet125->InsertPoint(2, m_FiducialMarkerCentroids.at(5)); //Add the points projected onto the opposite (parallel) plane to the pointset // By means of these projected points the calculated plane // is orthogonal to the plane 1234 and 5678. pointSet125->InsertPoint(3, planeGeometry5678->ProjectPointOntoPlane(m_FiducialMarkerCentroids.at(1))); pointSet125->InsertPoint(4, planeGeometry5678->ProjectPointOntoPlane(m_FiducialMarkerCentroids.at(2))); pointSet125->InsertPoint(5, planeGeometry1234->ProjectPointOntoPlane(m_FiducialMarkerCentroids.at(5))); pointSet478->InsertPoint(0, m_FiducialMarkerCentroids.at(4)); pointSet478->InsertPoint(1, m_FiducialMarkerCentroids.at(7)); pointSet478->InsertPoint(2, m_FiducialMarkerCentroids.at(8)); //Add the points projected onto the opposite (parallel) plane to the pointset // By means of these projected points the calculated plane // is orthogonal to the plane 1234 and 5678. pointSet478->InsertPoint(3, planeGeometry5678->ProjectPointOntoPlane(m_FiducialMarkerCentroids.at(4))); pointSet478->InsertPoint(4, planeGeometry1234->ProjectPointOntoPlane(m_FiducialMarkerCentroids.at(7))); pointSet478->InsertPoint(5, planeGeometry1234->ProjectPointOntoPlane(m_FiducialMarkerCentroids.at(8))); //Make planes parallel this->CreateParallelPlanes(planeFit125, planeFit478, pointSet125, pointSet478, planeGeometry125, planeGeometry478, false); this->MovePlanes(planeGeometry125, planeGeometry478, (2 * this->GetMinimalFiducialConfigurationDistance())); //Move the points into the parallel planes this->MovePoint(planeGeometry125, 1); this->MovePoint(planeGeometry125, 2); this->MovePoint(planeGeometry125, 5); this->MovePoint(planeGeometry478, 4); this->MovePoint(planeGeometry478, 7); this->MovePoint(planeGeometry478, 8); MITK_INFO << "NormalVector plane 125 = " << planeGeometry125->GetNormal(); MITK_INFO << "NormalVector plane 478 = " << planeGeometry478->GetNormal(); MITK_INFO << "Are parallel: " << planeGeometry125->IsParallel(planeGeometry478); //Optimize now the positions of Fiducials 1 - 3 - 8 and 4 - 5 - 6 //Initialization for planes 1 - 3 - 8 and 4 - 5 - 6 mitk::PlaneFit::Pointer planeFit138 = mitk::PlaneFit::New(); mitk::PlaneFit::Pointer planeFit456 = mitk::PlaneFit::New(); mitk::PointSet::Pointer pointSet138 = mitk::PointSet::New(); mitk::PointSet::Pointer pointSet456 = mitk::PointSet::New(); mitk::PlaneGeometry::Pointer planeGeometry138 = mitk::PlaneGeometry::New(); mitk::PlaneGeometry::Pointer planeGeometry456 = mitk::PlaneGeometry::New(); pointSet138->InsertPoint(0, m_FiducialMarkerCentroids.at(1)); pointSet138->InsertPoint(1, m_FiducialMarkerCentroids.at(3)); pointSet138->InsertPoint(2, m_FiducialMarkerCentroids.at(8)); //Add the points projected onto the opposite (parallel) plane to the pointset // By means of these projected points the calculated plane // is orthogonal to the plane 1234 and 5678. pointSet138->InsertPoint(3, planeGeometry5678->ProjectPointOntoPlane(m_FiducialMarkerCentroids.at(1))); pointSet138->InsertPoint(4, planeGeometry5678->ProjectPointOntoPlane(m_FiducialMarkerCentroids.at(3))); pointSet138->InsertPoint(5, planeGeometry1234->ProjectPointOntoPlane(m_FiducialMarkerCentroids.at(8))); pointSet456->InsertPoint(0, m_FiducialMarkerCentroids.at(4)); pointSet456->InsertPoint(1, m_FiducialMarkerCentroids.at(5)); pointSet456->InsertPoint(2, m_FiducialMarkerCentroids.at(6)); //Add the points projected onto the opposite (parallel) plane to the pointset // By means of these projected points the calculated plane // is orthogonal to the plane 1234 and 5678. pointSet456->InsertPoint(3, planeGeometry5678->ProjectPointOntoPlane(m_FiducialMarkerCentroids.at(4))); pointSet456->InsertPoint(4, planeGeometry1234->ProjectPointOntoPlane(m_FiducialMarkerCentroids.at(5))); pointSet456->InsertPoint(5, planeGeometry1234->ProjectPointOntoPlane(m_FiducialMarkerCentroids.at(6))); //Make planes parallel this->CreateParallelPlanes(planeFit138, planeFit456, pointSet138, pointSet456, planeGeometry138, planeGeometry456, false); this->MovePlanes(planeGeometry138, planeGeometry456, (2 * this->GetMinimalFiducialConfigurationDistance())); //Move the points into the parallel planes this->MovePoint(planeGeometry138, 1); this->MovePoint(planeGeometry138, 3); this->MovePoint(planeGeometry138, 8); this->MovePoint(planeGeometry456, 4); this->MovePoint(planeGeometry456, 5); this->MovePoint(planeGeometry456, 6); MITK_INFO << "NormalVector plane 138 = " << planeGeometry138->GetNormal(); MITK_INFO << "NormalVector plane 456 = " << planeGeometry456->GetNormal(); MITK_INFO << "Are parallel: " << planeGeometry138->IsParallel(planeGeometry456); } void QmitkUSNavigationStepCtUsRegistration::CreateParallelPlanes( mitk::PlaneFit::Pointer planeA, mitk::PlaneFit::Pointer planeB, mitk::PointSet::Pointer pointSetA, mitk::PointSet::Pointer pointSetB, mitk::PlaneGeometry::Pointer planeGeometryA, mitk::PlaneGeometry::Pointer planeGeometryB, bool minimizeInfluenceOutliers ) { planeA->SetInput(pointSetA); planeA->Update(); mitk::PlaneGeometry::Pointer geometryA = dynamic_cast(planeA->GetOutput()->GetGeometry()); planeB->SetInput(pointSetB); planeB->Update(); mitk::PlaneGeometry::Pointer geometryB = dynamic_cast(planeB->GetOutput()->GetGeometry()); mitk::DataNode::Pointer node1 = mitk::DataNode::New(); MITK_INFO << "Angle before = " << geometryA->Angle(geometryB) * 57.29578; //Transform into degree //Minimize influence of outliers concerning the inclination angle of the plane: if (minimizeInfluenceOutliers) { bool minimizeA = true; bool minimizeB = true; while(minimizeA) { MITK_INFO << "Minimize A"; std::vector> distancesA; distancesA.push_back(std::pair(geometryA->Distance(m_FiducialMarkerCentroids.at(1)), 1)); distancesA.push_back(std::pair(geometryA->Distance(m_FiducialMarkerCentroids.at(2)), 2)); distancesA.push_back(std::pair(geometryA->Distance(m_FiducialMarkerCentroids.at(3)), 3)); distancesA.push_back(std::pair(geometryA->Distance(m_FiducialMarkerCentroids.at(4)), 4)); std::sort(distancesA.begin(), distancesA.end()); for (std::vector>::iterator it = distancesA.begin(); it != distancesA.end(); ++it) { MITK_INFO << "First = " << (*it).first << " Second = " << (*it).second; } if ((*distancesA.rbegin()).first < 0.005) { minimizeA = false; break; } std::vector>::reverse_iterator it = distancesA.rbegin(); int fiducialNoToBeMovedToPlane1 = (*it).second; ++it; int fiducialNoToBeMovedToPlane2 = (*it).second; this->MovePoint(geometryA, fiducialNoToBeMovedToPlane1); this->MovePoint(geometryA, fiducialNoToBeMovedToPlane2); pointSetA->Clear(); for (int counter = 1; counter <= 4; ++counter) { pointSetA->InsertPoint(counter - 1, m_FiducialMarkerCentroids.at(counter)); } planeA = mitk::PlaneFit::New(); planeA->SetInput(pointSetA); planeA->Update(); geometryA = dynamic_cast(planeA->GetOutput()->GetGeometry()); } while (minimizeB) { MITK_INFO << "Minimize B"; std::vector> distancesB; distancesB.push_back(std::pair(geometryB->Distance(m_FiducialMarkerCentroids.at(5)), 5)); distancesB.push_back(std::pair(geometryB->Distance(m_FiducialMarkerCentroids.at(6)), 6)); distancesB.push_back(std::pair(geometryB->Distance(m_FiducialMarkerCentroids.at(7)), 7)); distancesB.push_back(std::pair(geometryB->Distance(m_FiducialMarkerCentroids.at(8)), 8)); std::sort(distancesB.begin(), distancesB.end()); for (std::vector>::iterator it = distancesB.begin(); it != distancesB.end(); ++it) { MITK_INFO << "First = " << (*it).first << " Second = " << (*it).second; } if ((*distancesB.rbegin()).first < 0.005) { minimizeB = false; break; } std::vector>::reverse_iterator it = distancesB.rbegin(); int fiducialNoToBeMovedToPlane1 = (*it).second; ++it; int fiducialNoToBeMovedToPlane2 = (*it).second; this->MovePoint(geometryB, fiducialNoToBeMovedToPlane1); this->MovePoint(geometryB, fiducialNoToBeMovedToPlane2); pointSetB->Clear(); for (int counter = 5; counter <= 8; ++counter) { pointSetB->InsertPoint(counter - 5, m_FiducialMarkerCentroids.at(counter)); } planeB = mitk::PlaneFit::New(); planeB->SetInput(pointSetB); planeB->Update(); geometryB = dynamic_cast(planeB->GetOutput()->GetGeometry()); } MITK_INFO << "Angle after = " << geometryA->Angle(geometryB) * 57.29578; //Transform into degree } // End outliers minimization mitk::Point3D originPlaneA = geometryA->GetOrigin(); mitk::Point3D originPlaneB = geometryB->GetOrigin(); mitk::Vector3D normalPlaneA = geometryA->GetNormal(); mitk::Vector3D normalPlaneB = geometryB->GetNormal(); double scalarProduct = (normalPlaneA * normalPlaneB) / (normalPlaneA.GetNorm() * normalPlaneB.GetNorm()); double alpha = acos(scalarProduct) * 57.29578; //Transform into degree if (alpha > 90) { normalPlaneB *= -1; } mitk::Vector3D combinedNormalPlaneA_B = 0.5 * (normalPlaneA + normalPlaneB); combinedNormalPlaneA_B.Normalize(); planeGeometryA->InitializePlane(originPlaneA, combinedNormalPlaneA_B); planeGeometryB->InitializePlane(originPlaneB, combinedNormalPlaneA_B); } void QmitkUSNavigationStepCtUsRegistration::MovePlanes(mitk::PlaneGeometry::Pointer planeA, mitk::PlaneGeometry::Pointer planeB, double referenceDistance) { const mitk::PlaneGeometry *constPlaneB = planeB; double distanceBetweenPlanes = planeA->DistanceFromPlane(constPlaneB); MITK_INFO << "Distance between planes before moving = " << distanceBetweenPlanes; //If distanceBetweenPlanes is > referenceDistance, the result will be negative, // otherwise it will be positive: double distanceToMove = 0.5 * (referenceDistance - distanceBetweenPlanes); const mitk::Vector3D movingVectorPositive = distanceToMove * planeA->GetNormal(); const mitk::Vector3D movingVectorNegative = distanceToMove * planeA->GetNormal() * -1; //Check, whether the planeB is above planeA if (planeA->IsAbove(planeB->GetOrigin())) { //planeB is above planeA planeA->Translate(movingVectorNegative); planeB->Translate(movingVectorPositive); } else { //planeB is below planeA planeA->Translate(movingVectorPositive); planeB->Translate(movingVectorNegative); } const mitk::PlaneGeometry *constPlaneBMoved = planeB; MITK_INFO << "Distance between planes after moving = " << planeA->DistanceFromPlane(constPlaneBMoved); } void QmitkUSNavigationStepCtUsRegistration::MovePoint( mitk::PlaneGeometry::Pointer planeGeometry, int fiducialNo) { MITK_INFO << "Moved Point from " << m_FiducialMarkerCentroids.at(fiducialNo); //Move all points outside the plane to the plane (finally all points lie into the plane): mitk::Point3D movedPoint = planeGeometry-> ProjectPointOntoPlane(m_FiducialMarkerCentroids.at(fiducialNo)); m_FiducialMarkerCentroids.at(fiducialNo)[0] = movedPoint[0]; m_FiducialMarkerCentroids.at(fiducialNo)[1] = movedPoint[1]; m_FiducialMarkerCentroids.at(fiducialNo)[2] = movedPoint[2]; MITK_INFO << " to " << m_FiducialMarkerCentroids.at(fiducialNo); } void QmitkUSNavigationStepCtUsRegistration::DefineDataStorageImageFilter() { m_IsAPointSetPredicate = mitk::TNodePredicateDataType::New(); mitk::TNodePredicateDataType::Pointer isImage = mitk::TNodePredicateDataType::New(); auto isSegmentation = mitk::NodePredicateDataType::New("Segment"); mitk::NodePredicateOr::Pointer validImages = mitk::NodePredicateOr::New(); validImages->AddPredicate(mitk::NodePredicateAnd::New(isImage, mitk::NodePredicateNot::New(isSegmentation))); mitk::NodePredicateNot::Pointer isNotAHelperObject = mitk::NodePredicateNot::New(mitk::NodePredicateProperty::New("helper object", mitk::BoolProperty::New(true))); m_IsOfTypeImagePredicate = mitk::NodePredicateAnd::New(validImages, isNotAHelperObject); mitk::NodePredicateProperty::Pointer isBinaryPredicate = mitk::NodePredicateProperty::New("binary", mitk::BoolProperty::New(true)); mitk::NodePredicateNot::Pointer isNotBinaryPredicate = mitk::NodePredicateNot::New(isBinaryPredicate); mitk::NodePredicateAnd::Pointer isABinaryImagePredicate = mitk::NodePredicateAnd::New(m_IsOfTypeImagePredicate, isBinaryPredicate); mitk::NodePredicateAnd::Pointer isNotABinaryImagePredicate = mitk::NodePredicateAnd::New(m_IsOfTypeImagePredicate, isNotBinaryPredicate); m_IsASegmentationImagePredicate = mitk::NodePredicateOr::New(isABinaryImagePredicate, mitk::TNodePredicateDataType::New()); m_IsAPatientImagePredicate = mitk::NodePredicateAnd::New(isNotABinaryImagePredicate, mitk::NodePredicateNot::New(mitk::TNodePredicateDataType::New())); } void QmitkUSNavigationStepCtUsRegistration::CreateQtPartControl(QWidget *parent) { ui->setupUi(parent); ui->floatingImageComboBox->SetPredicate(m_IsAPatientImagePredicate); ui->groundTruthImageComboBox->SetPredicate(m_IsAPatientImagePredicate); ui->ctImagesToChooseComboBox->SetPredicate(m_IsAPatientImagePredicate); // create signal/slot connections connect(ui->floatingImageComboBox, SIGNAL(OnSelectionChanged(const mitk::DataNode*)), this, SLOT(OnFloatingImageComboBoxSelectionChanged(const mitk::DataNode*))); connect(ui->groundTruthImageComboBox, SIGNAL(OnSelectionChanged(const mitk::DataNode*)), this, SLOT(OnGroundTruthImageComboBoxSelectionChanged(const mitk::DataNode*))); connect(ui->doRegistrationMarkerToImagePushButton, SIGNAL(clicked()), this, SLOT(OnRegisterMarkerToFloatingImageCS())); connect(ui->localizeFiducialMarkerPushButton, SIGNAL(clicked()), this, SLOT(OnLocalizeFiducials())); connect(ui->filterImageGroundTruthEvaluationPushButton, SIGNAL(clicked()), this, SLOT(OnFilterGroundTruthImage())); connect(ui->addCtImagePushButton, SIGNAL(clicked()), this, SLOT(OnAddCtImageClicked())); connect(ui->removeCtImagePushButton, SIGNAL(clicked()), this, SLOT(OnRemoveCtImageClicked())); connect(ui->evaluateProtocolPushButton, SIGNAL(clicked()), this, SLOT(OnEvaluateGroundTruthFiducialLocalizationProtocol())); } void QmitkUSNavigationStepCtUsRegistration::OnFloatingImageComboBoxSelectionChanged(const mitk::DataNode* node) { MITK_INFO << "OnFloatingImageComboBoxSelectionChanged()"; if (m_FloatingImage.IsNotNull()) { //TODO: Define, what will happen if the imageCT is not null... } if (node == nullptr) { this->UnsetFloatingImageGeometry(); m_FloatingImage = nullptr; return; } mitk::DataNode* selectedFloatingImage = ui->floatingImageComboBox->GetSelectedNode(); if (selectedFloatingImage == nullptr) { this->UnsetFloatingImageGeometry(); m_FloatingImage = nullptr; return; } mitk::Image::Pointer floatingImage = dynamic_cast(selectedFloatingImage->GetData()); if (floatingImage.IsNull()) { MITK_WARN << "Failed to cast selected segmentation node to mitk::Image*"; this->UnsetFloatingImageGeometry(); m_FloatingImage = nullptr; return; } m_FloatingImage = floatingImage; this->SetFloatingImageGeometryInformation(floatingImage.GetPointer()); } void QmitkUSNavigationStepCtUsRegistration::OnGroundTruthImageComboBoxSelectionChanged(const mitk::DataNode* node) { MITK_INFO << "OnGroundTruthImageComboBoxSelectionChanged()"; if (m_GroundTruthImage.IsNotNull()) { //TODO: Define, what will happen if the imageCT is not null... } if (node == nullptr) { m_GroundTruthImage = nullptr; return; } mitk::DataNode* selectedGroundTruthImage = ui->groundTruthImageComboBox->GetSelectedNode(); if (selectedGroundTruthImage == nullptr) { m_GroundTruthImage = nullptr; return; } mitk::Image::Pointer groundTruthImage = dynamic_cast(selectedGroundTruthImage->GetData()); if (groundTruthImage.IsNull()) { MITK_WARN << "Failed to cast selected groundTruth image node to mitk::Image*"; m_GroundTruthImage = nullptr; return; } m_GroundTruthImage = groundTruthImage; } void QmitkUSNavigationStepCtUsRegistration::OnRegisterMarkerToFloatingImageCS() { this->CreateMarkerModelCoordinateSystemPointSet(); //Check for initialization if( m_MarkerModelCoordinateSystemPointSet.IsNull() || m_MarkerFloatingImageCoordinateSystemPointSet.IsNull() ) { MITK_WARN << "Fiducial Landmarks are not initialized yet, cannot register"; return; } //Retrieve fiducials if (m_MarkerFloatingImageCoordinateSystemPointSet->GetSize() != m_MarkerModelCoordinateSystemPointSet->GetSize()) { MITK_WARN << "Not the same number of fiducials, cannot register"; return; } else if (m_MarkerFloatingImageCoordinateSystemPointSet->GetSize() < 3) { MITK_WARN << "Need at least 3 fiducials, cannot register"; return; } //############### conversion to vtk data types (we will use the vtk landmark based transform) ########################## //convert point sets to vtk poly data vtkSmartPointer sourcePoints = vtkSmartPointer::New(); vtkSmartPointer targetPoints = vtkSmartPointer::New(); for (int i = 0; iGetSize(); i++) { double point[3] = { m_MarkerModelCoordinateSystemPointSet->GetPoint(i)[0], m_MarkerModelCoordinateSystemPointSet->GetPoint(i)[1], m_MarkerModelCoordinateSystemPointSet->GetPoint(i)[2] }; sourcePoints->InsertNextPoint(point); double point_targets[3] = { m_MarkerFloatingImageCoordinateSystemPointSet->GetPoint(i)[0], m_MarkerFloatingImageCoordinateSystemPointSet->GetPoint(i)[1], m_MarkerFloatingImageCoordinateSystemPointSet->GetPoint(i)[2] }; targetPoints->InsertNextPoint(point_targets); } //########################### here, the actual transform is computed ########################## //compute transform vtkSmartPointer transform = vtkSmartPointer::New(); transform->SetSourceLandmarks(sourcePoints); transform->SetTargetLandmarks(targetPoints); transform->SetModeToRigidBody(); transform->Modified(); transform->Update(); //compute FRE of transform double FRE = mitk::StaticIGTHelperFunctions::ComputeFRE(m_MarkerModelCoordinateSystemPointSet, m_MarkerFloatingImageCoordinateSystemPointSet, transform); MITK_INFO << "FRE: " << FRE << " mm"; if (m_PerformingGroundTruthProtocolEvaluation) { m_GroundTruthProtocolFRE.push_back(FRE); } //############################################################################################# //############### conversion back to itk/mitk data types ########################## //convert from vtk to itk data types itk::Matrix rotationFloat = itk::Matrix(); itk::Vector translationFloat = itk::Vector(); itk::Matrix rotationDouble = itk::Matrix(); itk::Vector translationDouble = itk::Vector(); vtkSmartPointer m = transform->GetMatrix(); for (int k = 0; k<3; k++) for (int l = 0; l<3; l++) { rotationFloat[k][l] = m->GetElement(k, l); rotationDouble[k][l] = m->GetElement(k, l); } for (int k = 0; k<3; k++) { translationFloat[k] = m->GetElement(k, 3); translationDouble[k] = m->GetElement(k, 3); } //create mitk affine transform 3D and save it to the class member m_TransformMarkerCSToFloatingImageCS = mitk::AffineTransform3D::New(); m_TransformMarkerCSToFloatingImageCS->SetMatrix(rotationDouble); m_TransformMarkerCSToFloatingImageCS->SetOffset(translationDouble); MITK_INFO << m_TransformMarkerCSToFloatingImageCS; //################################################################ //############### object is transformed ########################## //transform surface/image //only move image if we have one. Sometimes, this widget is used just to register point sets without images. /*if (m_ImageNode.IsNotNull()) { //first we have to store the original ct image transform to compose it with the new transform later mitk::AffineTransform3D::Pointer imageTransform = m_ImageNode->GetData()->GetGeometry()->GetIndexToWorldTransform(); imageTransform->Compose(mitkTransform); mitk::AffineTransform3D::Pointer newImageTransform = mitk::AffineTransform3D::New(); //create new image transform... setting the composed directly leads to an error itk::Matrix rotationFloatNew = imageTransform->GetMatrix(); itk::Vector translationFloatNew = imageTransform->GetOffset(); newImageTransform->SetMatrix(rotationFloatNew); newImageTransform->SetOffset(translationFloatNew); m_ImageNode->GetData()->GetGeometry()->SetIndexToWorldTransform(newImageTransform); }*/ //If this option is set, each point will be transformed and the acutal coordinates of the points change. if( !m_PerformingGroundTruthProtocolEvaluation ) { mitk::PointSet* pointSet_orig = m_MarkerModelCoordinateSystemPointSet; mitk::PointSet::Pointer pointSet_moved = mitk::PointSet::New(); for (int i = 0; i < pointSet_orig->GetSize(); i++) { pointSet_moved->InsertPoint(m_TransformMarkerCSToFloatingImageCS->TransformPoint(pointSet_orig->GetPoint(i))); } pointSet_orig->Clear(); for (int i = 0; i < pointSet_moved->GetSize(); i++) pointSet_orig->InsertPoint(pointSet_moved->GetPoint(i)); //Do a global reinit mitk::RenderingManager::GetInstance()->InitializeViewsByBoundingObjects(this->GetDataStorage()); } } void QmitkUSNavigationStepCtUsRegistration::OnLocalizeFiducials() { m_FiducialMarkerCentroids.clear(); m_CentroidsOfFiducialCandidates.clear(); if (m_MarkerFloatingImageCoordinateSystemPointSet.IsNotNull()) { m_MarkerFloatingImageCoordinateSystemPointSet->Clear(); } if (!this->FilterFloatingImage()) { QMessageBox msgBox; msgBox.setText("Cannot perform filtering of the image. The floating image = nullptr."); msgBox.exec(); return; } this->GetCentroidsOfLabeledObjects(); if (!this->EliminateFiducialCandidatesByEuclideanDistances() || m_CentroidsOfFiducialCandidates.size() != NUMBER_FIDUCIALS_NEEDED) { QMessageBox msgBox; QString text = QString("Have found %1 instead of 8 fiducial candidates.\ Cannot perform fiducial localization procedure.").arg(m_CentroidsOfFiducialCandidates.size()); msgBox.setText(text); msgBox.exec(); return; } //Before calling NumerateFiducialMarks it must be sure, // that there rested only 8 fiducial candidates. this->NumerateFiducialMarks(); } void QmitkUSNavigationStepCtUsRegistration::OnFilterGroundTruthImage() { /*if (m_GroundTruthImage.IsNull()) { return; }*/ this->ShowGroundTruthMarkerEdges(); } void QmitkUSNavigationStepCtUsRegistration::OnAddCtImageClicked() { mitk::DataNode* selectedCtImage = ui->ctImagesToChooseComboBox->GetSelectedNode(); if (selectedCtImage == nullptr) { return; } mitk::Image::Pointer ctImage = dynamic_cast(selectedCtImage->GetData()); if (ctImage.IsNull()) { MITK_WARN << "Failed to cast selected segmentation node to mitk::Image*"; return; } QString name = QString::fromStdString(selectedCtImage->GetName()); for( int counter = 0; counter < ui->chosenCtImagesListWidget->count(); ++counter) { MITK_INFO << ui->chosenCtImagesListWidget->item(counter)->text() << " - " << counter; MITK_INFO << m_ImagesGroundTruthProtocol.at(counter).GetPointer(); if (ui->chosenCtImagesListWidget->item(counter)->text().compare(name) == 0) { MITK_INFO << "CT image already exist in list of chosen CT images. Do not add the image."; return; } } ui->chosenCtImagesListWidget->addItem(name); m_ImagesGroundTruthProtocol.push_back(ctImage); } void QmitkUSNavigationStepCtUsRegistration::OnRemoveCtImageClicked() { int position = ui->chosenCtImagesListWidget->currentRow(); if (ui->chosenCtImagesListWidget->count() == 0 || position < 0) { return; } m_ImagesGroundTruthProtocol.erase(m_ImagesGroundTruthProtocol.begin() + position); QListWidgetItem *item = ui->chosenCtImagesListWidget->currentItem(); ui->chosenCtImagesListWidget->removeItemWidget(item); delete item; } void QmitkUSNavigationStepCtUsRegistration::OnEvaluateGroundTruthFiducialLocalizationProtocol() { m_GroundTruthProtocolFRE.clear(); if (m_ImagesGroundTruthProtocol.size() != 6) { QMessageBox msgBox; msgBox.setText("For evaluating the Ground-Truth-Fiducial-Localization-Protocol there must be loaded 6 different CT images."); msgBox.exec(); return; } m_PerformingGroundTruthProtocolEvaluation = true; this->CreatePointsToTransformForGroundTruthProtocol(); m_GroundTruthProtocolTransformedPoints.clear(); for (int cycleNo = 0; cycleNo < m_ImagesGroundTruthProtocol.size(); ++cycleNo) { m_FloatingImage = m_ImagesGroundTruthProtocol.at(cycleNo); this->SetFloatingImageGeometryInformation(m_FloatingImage.GetPointer()); this->OnLocalizeFiducials(); this->OnRegisterMarkerToFloatingImageCS(); this->TransformPointsGroundTruthProtocol(); } this->AddTransformedPointsToDataStorage(); double meanFRE = this->CalculateMeanFRE(); double sdOfFRE = this->CalculateStandardDeviationOfFRE(meanFRE); this->CalculateGroundTruthProtocolTRE(); ui->meanFREValue->setText(QString("%1").arg(meanFRE)); ui->sdFREValue->setText(QString("%1").arg(sdOfFRE)); if (ui->protocolEvaluationTypeComboBox->currentText().compare("ANGLE") == 0) { if (m_GroundTruthProtocolTRE.find(0) != m_GroundTruthProtocolTRE.end()) { ui->TREValue->setText(QString("%1").arg(m_GroundTruthProtocolTRE.at(0))); } } else if (ui->protocolEvaluationTypeComboBox->currentText().compare("PLANE") == 0) { if (m_GroundTruthProtocolTRE.find(0) != m_GroundTruthProtocolTRE.end() && m_GroundTruthProtocolTRE.find(20) != m_GroundTruthProtocolTRE.end() && m_GroundTruthProtocolTRE.find(40) != m_GroundTruthProtocolTRE.end() && m_GroundTruthProtocolTRE.find(60) != m_GroundTruthProtocolTRE.end() && m_GroundTruthProtocolTRE.find(80) != m_GroundTruthProtocolTRE.end() && m_GroundTruthProtocolTRE.find(100) != m_GroundTruthProtocolTRE.end()) { ui->TREValue->setText(QString("Depth 0mm: %1\nDepth 20mm: %2\nDepth 40mm: %3\ \nDepth 60mm: %4\nDepth 80mm: %5\nDepth 100mm: %6") .arg(m_GroundTruthProtocolTRE.at(0)) .arg(m_GroundTruthProtocolTRE.at(20)) .arg(m_GroundTruthProtocolTRE.at(40)) .arg(m_GroundTruthProtocolTRE.at(60)) .arg(m_GroundTruthProtocolTRE.at(80)) .arg(m_GroundTruthProtocolTRE.at(100))); } } m_PerformingGroundTruthProtocolEvaluation = false; } diff --git a/Plugins/org.mitk.gui.qt.igt.app.echotrack/src/internal/QmitkUltrasoundCalibration.cpp b/Plugins/org.mitk.gui.qt.igt.app.echotrack/src/internal/QmitkUltrasoundCalibration.cpp index 019828609b..cf6efdaeef 100644 --- a/Plugins/org.mitk.gui.qt.igt.app.echotrack/src/internal/QmitkUltrasoundCalibration.cpp +++ b/Plugins/org.mitk.gui.qt.igt.app.echotrack/src/internal/QmitkUltrasoundCalibration.cpp @@ -1,1158 +1,1164 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ // Blueberry #include #include // Qmitk #include "QmitkUltrasoundCalibration.h" #include // Qt #include #include #include #include // MITK #include //#include #include #include #include #include #include #include "mitkIRenderingManager.h" // us #include //VTK #include #include #include #include #include #include #include "internal/org_mbi_gui_qt_usnavigation_Activator.h" //sleep headers #include #include const std::string QmitkUltrasoundCalibration::VIEW_ID = "org.mitk.views.ultrasoundcalibration"; QmitkUltrasoundCalibration::QmitkUltrasoundCalibration() : m_USDeviceChanged(this, &QmitkUltrasoundCalibration::OnUSDepthChanged) { ctkPluginContext* pluginContext = mitk::PluginActivator::GetContext(); if (pluginContext) { // to be notified about service event of an USDevice pluginContext->connectServiceListener(this, "OnDeviceServiceEvent", QString::fromStdString("(" + us::ServiceConstants::OBJECTCLASS() + "=" + us_service_interface_iid() + ")")); } } QmitkUltrasoundCalibration::~QmitkUltrasoundCalibration() { m_Controls.m_CombinedModalityManagerWidget->blockSignals(true); mitk::AbstractUltrasoundTrackerDevice::Pointer combinedModality; combinedModality = m_Controls.m_CombinedModalityManagerWidget->GetSelectedCombinedModality(); if (combinedModality.IsNotNull()) { combinedModality->GetUltrasoundDevice()->RemovePropertyChangedListener(m_USDeviceChanged); } m_Timer->stop(); // Sleep(500); //This might be problematic... seems like sometimes some ressources are still in use at calling time. this->OnStopCalibrationProcess(); this->OnStopPlusCalibration(); /*mitk::DataNode::Pointer node = this->GetDataStorage()->GetNamedNode("Tool Calibration Points"); if (node.IsNotNull())this->GetDataStorage()->Remove(node); node = this->GetDataStorage()->GetNamedNode("Image Calibration Points"); if (node.IsNotNull())this->GetDataStorage()->Remove(node); node = this->GetDataStorage()->GetNamedNode("US Image Stream"); if (node.IsNotNull())this->GetDataStorage()->Remove(node);*/ mitk::DataNode::Pointer node = this->GetDataStorage()->GetNamedNode("Needle Path"); if (node.IsNotNull())this->GetDataStorage()->Remove(node); this->GetDataStorage()->Remove(m_VerificationReferencePointsDataNode); delete m_Timer; } void QmitkUltrasoundCalibration::SetFocus() { m_Controls.m_ToolBox->setFocus(); } void QmitkUltrasoundCalibration::CreateQtPartControl(QWidget *parent) { // create GUI widgets from the Qt Designer's .ui file m_Controls.setupUi(parent); m_Controls.m_CombinedModalityManagerWidget->SetCalibrationLoadedNecessary(false); m_Timer = new QTimer(this); m_StreamingTimer = new QTimer(this); m_Controls.m_SpacingBtnFreeze->setEnabled(true); m_Controls.m_SpacingAddPoint->setEnabled(false); m_Controls.m_CalculateSpacing->setEnabled(false); m_SpacingPointsCount = 0; m_SpacingPoints = mitk::PointSet::New(); m_SpacingNode = mitk::DataNode::New(); m_SpacingNode->SetName("Spacing Points"); m_SpacingNode->SetData(this->m_SpacingPoints); this->GetDataStorage()->Add(m_SpacingNode); // Pointset for Calibration Points m_CalibPointsTool = mitk::PointSet::New(); // Pointset for Worldpoints m_CalibPointsImage = mitk::PointSet::New(); m_CalibPointsCount = 0; // Evaluation Pointsets (Non-Visualized) m_EvalPointsImage = mitk::PointSet::New(); m_EvalPointsTool = mitk::PointSet::New(); m_EvalPointsProjected = mitk::PointSet::New(); // Neelde Projection Filter m_NeedleProjectionFilter = mitk::NeedleProjectionFilter::New(); // Tracking Status Widgets m_Controls.m_CalibTrackingStatus->ShowStatusLabels(); m_Controls.m_EvalTrackingStatus->ShowStatusLabels(); // General & Device Selection connect(m_Timer, SIGNAL(timeout()), this, SLOT(Update())); //connect(m_Controls.m_ToolBox, SIGNAL(currentChanged(int)), this, SLOT(OnTabSwitch(int))); // Calibration connect(m_Controls.m_CalibBtnFreeze, SIGNAL(clicked()), this, SLOT(SwitchFreeze())); // Freeze connect(m_Controls.m_CalibBtnAddPoint, SIGNAL(clicked()), this, SLOT(OnAddCalibPoint())); // Tracking & Image Points (Calibration) connect(m_Controls.m_CalibBtnCalibrate, SIGNAL(clicked()), this, SLOT(OnCalibration())); // Perform Calibration // Evaluation connect(m_Controls.m_EvalBtnStep1, SIGNAL(clicked()), this, SLOT(OnAddEvalProjectedPoint())); // Needle Projection connect(m_Controls.m_EvalBtnStep2, SIGNAL(clicked()), this, SLOT(SwitchFreeze())); // Freeze connect(m_Controls.m_EvalBtnStep3, SIGNAL(clicked()), this, SLOT(OnAddEvalTargetPoint())); // Tracking & Image Points (Evaluation) connect(m_Controls.m_EvalBtnSave, SIGNAL(clicked()), this, SLOT(OnSaveEvaluation())); // Save Evaluation Results connect(m_Controls.m_CalibBtnSaveCalibration, SIGNAL(clicked()), this, SLOT(OnSaveCalibration())); // Save Evaluation Results connect(m_Controls.m_BtnReset, SIGNAL(clicked()), this, SLOT(OnReset())); // Reset Pointsets // PLUS Calibration connect(m_Controls.m_GetCalibrationFromPLUS, SIGNAL(clicked()), this, SLOT(OnGetPlusCalibration())); connect(m_Controls.m_StartStreaming, SIGNAL(clicked()), this, SLOT(OnStartStreaming())); connect(m_StreamingTimer, SIGNAL(timeout()), this, SLOT(OnStreamingTimerTimeout())); connect(m_Controls.m_StopPlusCalibration, SIGNAL(clicked()), this, SLOT(OnStopPlusCalibration())); connect(m_Controls.m_SavePlusCalibration, SIGNAL(clicked()), this, SLOT(OnSaveCalibration())); connect(this, SIGNAL(NewConnectionSignal()), this, SLOT(OnNewConnection())); //Determine Spacing for Calibration of USVideoDevice connect(m_Controls.m_SpacingBtnFreeze, SIGNAL(clicked()), this, SLOT(OnFreezeClicked())); connect(m_Controls.m_SpacingAddPoint, SIGNAL(clicked()), this, SLOT(OnAddSpacingPoint())); connect(m_Controls.m_CalculateSpacing, SIGNAL(clicked()), this, SLOT(OnCalculateSpacing())); //connect( m_Controls.m_CombinedModalityManagerWidget, SIGNAL(SignalCombinedModalitySelected(mitk::USCombinedModality::Pointer)), // this, SLOT(OnSelectDevice(mitk::USCombinedModality::Pointer)) ); connect(m_Controls.m_CombinedModalityManagerWidget, SIGNAL(SignalReadyForNextStep()), this, SLOT(OnDeviceSelected())); connect(m_Controls.m_CombinedModalityManagerWidget, SIGNAL(SignalNoLongerReadyForNextStep()), this, SLOT(OnDeviceDeselected())); connect(m_Controls.m_StartCalibrationButton, SIGNAL(clicked()), this, SLOT(OnStartCalibrationProcess())); connect(m_Controls.m_StartPlusCalibrationButton, SIGNAL(clicked()), this, SLOT(OnStartPlusCalibration())); connect(m_Controls.m_CalibBtnRestartCalibration, SIGNAL(clicked()), this, SLOT(OnReset())); connect(m_Controls.m_CalibBtnStopCalibration, SIGNAL(clicked()), this, SLOT(OnStopCalibrationProcess())); connect(m_Controls.m_AddReferencePoints, SIGNAL(clicked()), this, SLOT(OnAddCurrentTipPositionToReferencePoints())); connect(m_Controls.m_AddCurrentPointerTipForVerification, SIGNAL(clicked()), this, SLOT(OnAddCurrentTipPositionForVerification())); connect(m_Controls.m_StartVerification, SIGNAL(clicked()), this, SLOT(OnStartVerification())); //initialize data storage combo box m_Controls.m_ReferencePointsComboBox->SetDataStorage(this->GetDataStorage()); m_Controls.m_ReferencePointsComboBox->SetAutoSelectNewItems(true); m_Controls.m_ReferencePointsComboBox->SetPredicate(mitk::NodePredicateDataType::New("PointSet")); //initialize point list widget if (m_VerificationReferencePoints.IsNull()) { m_VerificationReferencePoints = mitk::PointSet::New(); } if (m_VerificationReferencePointsDataNode.IsNull()) { m_VerificationReferencePointsDataNode = mitk::DataNode::New(); m_VerificationReferencePointsDataNode->SetName("US Verification Reference Points"); m_VerificationReferencePointsDataNode->SetData(m_VerificationReferencePoints); this->GetDataStorage()->Add(m_VerificationReferencePointsDataNode); } m_Controls.m_ReferencePointsPointListWidget->SetPointSetNode(m_VerificationReferencePointsDataNode); m_Controls.m_ToolBox->setCurrentIndex(0); } void QmitkUltrasoundCalibration::OnSelectionChanged(berry::IWorkbenchPart::Pointer /*source*/, const QList& /*nodes*/) { } void QmitkUltrasoundCalibration::OnTabSwitch(int index) { switch (index) { case 0: if (m_Controls.m_ToolBox->isItemEnabled(1) || m_Controls.m_ToolBox->isItemEnabled(2)) { this->OnStopCalibrationProcess(); } break; default: ; } } //void QmitkUltrasoundCalibration::OnSelectDevice(mitk::USCombinedModality::Pointer combinedModality) void QmitkUltrasoundCalibration::OnDeviceSelected() { mitk::AbstractUltrasoundTrackerDevice::Pointer combinedModality; combinedModality = m_Controls.m_CombinedModalityManagerWidget->GetSelectedCombinedModality(); if (combinedModality.IsNotNull()) { //m_Tracker = m_CombinedModality->GetNavigationDataSource(); // Construct Pipeline //this->m_NeedleProjectionFilter->SetInput(0, m_Tracker->GetOutput(0)); combinedModality->GetUltrasoundDevice()->AddPropertyChangedListener(m_USDeviceChanged); m_Controls.m_StartCalibrationButton->setEnabled(true); m_Controls.m_StartPlusCalibrationButton->setEnabled(true); m_Controls.m_ToolBox->setItemEnabled(1, true); m_Controls.m_ToolBox->setItemEnabled(2, true); } } void QmitkUltrasoundCalibration::OnDeviceDeselected() { mitk::AbstractUltrasoundTrackerDevice::Pointer combinedModality; combinedModality = m_Controls.m_CombinedModalityManagerWidget->GetSelectedCombinedModality(); if (combinedModality.IsNotNull()) { combinedModality->GetUltrasoundDevice()->RemovePropertyChangedListener(m_USDeviceChanged); } m_Controls.m_StartCalibrationButton->setEnabled(false); m_Controls.m_StartPlusCalibrationButton->setEnabled(false); m_Controls.m_ToolBox->setCurrentIndex(0); m_Controls.m_ToolBox->setItemEnabled(1, false); m_Controls.m_ToolBox->setItemEnabled(2, false); } void QmitkUltrasoundCalibration::OnAddCurrentTipPositionToReferencePoints() { if (m_Controls.m_VerificationPointerChoser->GetSelectedNavigationDataSource().IsNull() || (m_Controls.m_VerificationPointerChoser->GetSelectedToolID() == -1)) { MITK_WARN << "No tool selected, aborting"; return; } mitk::NavigationData::Pointer currentPointerData = m_Controls.m_VerificationPointerChoser->GetSelectedNavigationDataSource()->GetOutput(m_Controls.m_VerificationPointerChoser->GetSelectedToolID()); mitk::Point3D currentTipPosition = currentPointerData->GetPosition(); m_VerificationReferencePoints->InsertPoint(m_VerificationReferencePoints->GetSize(), currentTipPosition); } void QmitkUltrasoundCalibration::OnStartVerification() { m_currentPoint = 0; mitk::PointSet::Pointer selectedPointSet = dynamic_cast(m_Controls.m_ReferencePointsComboBox->GetSelectedNode()->GetData()); m_Controls.m_CurrentPointLabel->setText("Point " + QString::number(m_currentPoint) + " of " + QString::number(selectedPointSet->GetSize())); m_allErrors = std::vector(); m_allReferencePoints = std::vector(); for (int i = 0; i < selectedPointSet->GetSize(); i++) { m_allReferencePoints.push_back(selectedPointSet->GetPoint(i)); } } void QmitkUltrasoundCalibration::OnAddCurrentTipPositionForVerification() { if (m_currentPoint == -1) { MITK_WARN << "Cannot add point"; return; } if (m_Controls.m_VerificationPointerChoser->GetSelectedNavigationDataSource().IsNull() || (m_Controls.m_VerificationPointerChoser->GetSelectedToolID() == -1)) { MITK_WARN << "No tool selected, aborting"; return; } mitk::NavigationData::Pointer currentPointerData = m_Controls.m_VerificationPointerChoser->GetSelectedNavigationDataSource()->GetOutput(m_Controls.m_VerificationPointerChoser->GetSelectedToolID()); mitk::Point3D currentTipPosition = currentPointerData->GetPosition(); double currentError = m_allReferencePoints.at(m_currentPoint).EuclideanDistanceTo(currentTipPosition); MITK_INFO << "Current Error: " << currentError << " mm"; m_allErrors.push_back(currentError); if (++m_currentPoint < static_cast(m_allReferencePoints.size())) { m_Controls.m_CurrentPointLabel->setText("Point " + QString::number(m_currentPoint) + " of " + QString::number(m_allReferencePoints.size())); } else { m_currentPoint = -1; double meanError = 0; for (std::size_t i = 0; i < m_allErrors.size(); ++i) { meanError += m_allErrors[i]; } meanError /= m_allErrors.size(); QString result = "Finished verification! \n Verification of " + QString::number(m_allErrors.size()) + " points, mean error: " + QString::number(meanError) + " mm"; m_Controls.m_ResultsTextEdit->setText(result); MITK_INFO << result.toStdString(); } } void QmitkUltrasoundCalibration::OnStartCalibrationProcess() { // US Image Stream m_Node = mitk::DataNode::New(); m_Node->SetName("US Calibration Viewing Stream"); //create a dummy image (gray values 0..255) for correct initialization of level window, etc. mitk::Image::Pointer dummyImage = mitk::ImageGenerator::GenerateRandomImage(100, 100, 1, 1, 1, 1, 1, 255, 0); m_Node->SetData(dummyImage); this->GetDataStorage()->Add(m_Node); // data node for calibration point set m_CalibNode = mitk::DataNode::New(); m_CalibNode->SetName("Tool Calibration Points"); m_CalibNode->SetData(this->m_CalibPointsImage); this->GetDataStorage()->Add(m_CalibNode); // data node for world point set m_WorldNode = mitk::DataNode::New(); m_WorldNode->SetName("Image Calibration Points"); m_WorldNode->SetData(this->m_CalibPointsTool); this->GetDataStorage()->Add(m_WorldNode); m_CombinedModality = m_Controls.m_CombinedModalityManagerWidget->GetSelectedCombinedModality(); m_CombinedModality->SetCalibration(mitk::AffineTransform3D::New()); //dummy calibration because without a calibration the comined modality was laggy (maybe a bug?) if (m_CombinedModality.IsNull()) { return; } m_Tracker = m_CombinedModality->GetNavigationDataSource(); //QString curDepth = service.getProperty(QString::fromStdString(mitk::USDevice::US_PROPKEY_BMODE_DEPTH)).toString(); // Construct Pipeline this->m_NeedleProjectionFilter->SetInput(0, m_Tracker->GetOutput(0)); QApplication::setOverrideCursor(Qt::WaitCursor); // make sure that the combined modality is in connected state before using it if (m_CombinedModality->GetUltrasoundDevice()->GetDeviceState() < mitk::USDevice::State_Connected) { m_CombinedModality->GetUltrasoundDevice()->Connect(); } if (m_CombinedModality->GetUltrasoundDevice()->GetDeviceState() < mitk::USDevice::State_Activated) { m_CombinedModality->GetUltrasoundDevice()->Activate(); } QApplication::restoreOverrideCursor(); this->SwitchFreeze(); //Trigger the ProbeChanged method for initializing/updating the spacing of the ultrasound image correctly std::string probeName = m_CombinedModality->GetUltrasoundDevice()->GetCurrentProbe()->GetName(); m_CombinedModality->GetUltrasoundDevice()->ProbeChanged(probeName); + mitk::DataNode::Pointer usNode = this->GetDataStorage()->GetNamedNode("US Viewing Stream - Image 0"); + if (usNode.IsNotNull()) + { + this->GetDataStorage()->Remove(usNode); + } + // Todo: Maybe display this elsewhere this->ShowNeedlePath(); // Switch active tab to Calibration page m_Controls.m_ToolBox->setItemEnabled(1, true); m_Controls.m_ToolBox->setCurrentIndex(1); } void QmitkUltrasoundCalibration::OnStartPlusCalibration() { if (m_CombinedModality.IsNull()){ m_CombinedModality = m_Controls.m_CombinedModalityManagerWidget->GetSelectedCombinedModality(); if (m_CombinedModality.IsNull()) { return; } //something went wrong, there is no combined modality } //setup server to send UltrasoundImages to PLUS mitk::IGTLServer::Pointer m_USServer = mitk::IGTLServer::New(true); m_USServer->SetName("EchoTrack Image Source"); m_USServer->SetHostname("127.0.0.1"); m_USServer->SetPortNumber(18944); m_USMessageProvider = mitk::IGTLMessageProvider::New(); m_USMessageProvider->SetIGTLDevice(m_USServer); m_USMessageProvider->SetFPS(5); m_USImageToIGTLMessageFilter = mitk::ImageToIGTLMessageFilter::New(); m_USImageToIGTLMessageFilter->ConnectTo(m_CombinedModality->GetUltrasoundDevice()); m_USImageToIGTLMessageFilter->SetName("USImage Filter"); //setup server to send TrackingData to PLUS m_TrackingServer = mitk::IGTLServer::New(true); m_TrackingServer->SetName("EchoTrack Tracking Source"); m_TrackingServer->SetHostname("127.0.0.1"); m_TrackingServer->SetPortNumber(18945); m_TrackingMessageProvider = mitk::IGTLMessageProvider::New(); m_TrackingMessageProvider->SetIGTLDevice(m_TrackingServer); m_TrackingMessageProvider->SetFPS(5); m_TrackingToIGTLMessageFilter = mitk::NavigationDataToIGTLMessageFilter::New(); m_TrackingToIGTLMessageFilter->ConnectTo(m_CombinedModality->GetTrackingDeviceDataSource()); m_TrackingToIGTLMessageFilter->SetName("Tracker Filter"); typedef itk::SimpleMemberCommand< QmitkUltrasoundCalibration > CurCommandType; CurCommandType::Pointer newConnectionCommand = CurCommandType::New(); newConnectionCommand->SetCallbackFunction( this, &QmitkUltrasoundCalibration::OnPlusConnected); this->m_NewConnectionObserverTag = this->m_TrackingServer->AddObserver( mitk::NewClientConnectionEvent(), newConnectionCommand); //Open connections of both servers if (m_USServer->OpenConnection()) { MITK_INFO << "US Server opened its connection successfully"; m_USServer->StartCommunication(); } else { MITK_INFO << "US Server could not open its connection"; } if (m_TrackingServer->OpenConnection()) { MITK_INFO << "Tracking Server opened its connection successfully"; m_TrackingServer->StartCommunication(); } else { MITK_INFO << "Tracking Server could not open its connection"; } if (m_USMessageProvider->IsCommunicating() && m_TrackingMessageProvider->IsCommunicating()) { m_Controls.m_StartPlusCalibrationButton->setEnabled(false); m_Controls.m_GetCalibrationFromPLUS->setEnabled(true); m_Controls.m_StartStreaming->setEnabled(false); m_Controls.m_SavePlusCalibration->setEnabled(false); m_Controls.m_SetupStatus->setStyleSheet("QLabel { color : green; }"); m_Controls.m_SetupStatus->setText("Setup successfull you can now connect PLUS"); } else { m_Controls.m_SetupStatus->setStyleSheet("QLabel { color : red; }"); m_Controls.m_SetupStatus->setText("Something went wrong. Please try again"); } } void QmitkUltrasoundCalibration::OnStopPlusCalibration() { //closing all server and clients when PlusCalibration is finished if (m_USMessageProvider.IsNotNull()) { if (m_USMessageProvider->IsStreaming()) { m_USMessageProvider->StopStreamingOfSource(m_USImageToIGTLMessageFilter); } } if (m_TrackingMessageProvider.IsNotNull()) { if (m_TrackingMessageProvider->IsStreaming()) { m_TrackingMessageProvider->StopStreamingOfSource(m_TrackingToIGTLMessageFilter); } } if (m_USServer.IsNotNull()) { m_USServer->CloseConnection(); } if (m_TrackingServer.IsNotNull()) { m_TrackingServer->CloseConnection(); } if (m_TransformClient.IsNotNull()) { m_TransformClient->CloseConnection(); } m_Controls.m_GotCalibrationLabel->setText(""); m_Controls.m_ConnectionStatus->setText(""); m_Controls.m_SetupStatus->setText(""); m_Controls.m_StartPlusCalibrationButton->setEnabled(true); m_StreamingTimer->stop(); delete m_StreamingTimer; } void QmitkUltrasoundCalibration::OnPlusConnected() { emit NewConnectionSignal(); } void QmitkUltrasoundCalibration::OnNewConnection() { m_Controls.m_StartStreaming->setEnabled(true); m_Controls.m_ConnectionStatus->setStyleSheet("QLabel { color : green; }"); m_Controls.m_ConnectionStatus->setText("Connection successfull you can now start streaming"); } void QmitkUltrasoundCalibration::OnStreamingTimerTimeout() { m_USMessageProvider->Update(); m_TrackingMessageProvider->Update(); } void QmitkUltrasoundCalibration::OnStartStreaming() { m_USMessageProvider->StartStreamingOfSource(m_USImageToIGTLMessageFilter, 5); m_TrackingMessageProvider->StartStreamingOfSource(m_TrackingToIGTLMessageFilter, 5); m_Controls.m_StartStreaming->setEnabled(false); m_Controls.m_ConnectionStatus->setText(""); m_StreamingTimer->start((1.0 / 5.0 * 1000.0)); } void QmitkUltrasoundCalibration::OnGetPlusCalibration() { m_TransformClient = mitk::IGTLClient::New(true); m_TransformClient->SetHostname("127.0.0.1"); m_TransformClient->SetPortNumber(18946); m_TransformDeviceSource = mitk::IGTLDeviceSource::New(); m_TransformDeviceSource->SetIGTLDevice(m_TransformClient); m_TransformDeviceSource->Connect(); if (m_TransformDeviceSource->IsConnected()) { MITK_INFO << "successfully connected"; m_TransformDeviceSource->StartCommunication(); if (m_TransformDeviceSource->IsCommunicating()) { MITK_INFO << "communication started"; mitk::IGTLMessage::Pointer receivedMessage; bool condition = false; igtl::Matrix4x4 transformPLUS; while (!(receivedMessage.IsNotNull() && receivedMessage->IsDataValid())) { std::this_thread::sleep_for(std::chrono::milliseconds(50)); m_TransformDeviceSource->Update(); receivedMessage = m_TransformDeviceSource->GetOutput(); igtl::TransformMessage::Pointer msg = dynamic_cast(m_TransformDeviceSource->GetOutput()->GetMessage().GetPointer()); if (msg == nullptr || msg.IsNull()) { MITK_INFO << "Received message could not be casted to TransformMessage. Skipping.."; continue; } else { if (std::strcmp(msg->GetDeviceName(), "ImageToTracker") != 0) { MITK_INFO << "Was not Image to Tracker Transform. Skipping..."; continue; } else { msg->GetMatrix(transformPLUS); condition = true; break; } } } if (condition) { this->ProcessPlusCalibration(transformPLUS); } else { m_Controls.m_GotCalibrationLabel->setStyleSheet("QLabel { color : red; }"); m_Controls.m_GotCalibrationLabel->setText("Something went wrong. Please try again"); } } else { MITK_INFO << " no connection"; m_Controls.m_GotCalibrationLabel->setStyleSheet("QLabel { color : red; }"); m_Controls.m_GotCalibrationLabel->setText("Something went wrong. Please try again"); } } else { m_Controls.m_GotCalibrationLabel->setStyleSheet("QLabel { color : red; }"); m_Controls.m_GotCalibrationLabel->setText("Something went wrong. Please try again"); } } void QmitkUltrasoundCalibration::ProcessPlusCalibration(igtl::Matrix4x4& imageToTracker) { mitk::AffineTransform3D::Pointer imageToTrackerTransform = mitk::AffineTransform3D::New(); itk::Matrix rotationFloat = itk::Matrix(); itk::Vector translationFloat = itk::Vector(); rotationFloat[0][0] = imageToTracker[0][0]; rotationFloat[0][1] = imageToTracker[0][1]; rotationFloat[0][2] = imageToTracker[0][2]; rotationFloat[1][0] = imageToTracker[1][0]; rotationFloat[1][1] = imageToTracker[1][1]; rotationFloat[1][2] = imageToTracker[1][2]; rotationFloat[2][0] = imageToTracker[2][0]; rotationFloat[2][1] = imageToTracker[2][1]; rotationFloat[2][2] = imageToTracker[2][2]; translationFloat[0] = imageToTracker[0][3]; translationFloat[1] = imageToTracker[1][3]; translationFloat[2] = imageToTracker[2][3]; imageToTrackerTransform->SetTranslation(translationFloat); imageToTrackerTransform->SetMatrix(rotationFloat); m_CombinedModality->SetCalibration(imageToTrackerTransform); m_Controls.m_ToolBox->setItemEnabled(2, true); m_Controls.m_SavePlusCalibration->setEnabled(true); m_Controls.m_GotCalibrationLabel->setStyleSheet("QLabel { color : green; }"); m_Controls.m_GotCalibrationLabel->setText("Recieved Calibration from PLUS you can now save it"); } void QmitkUltrasoundCalibration::OnStopCalibrationProcess() { this->ClearTemporaryMembers(); m_Timer->stop(); this->GetDataStorage()->Remove(m_Node); m_Node = 0; this->GetDataStorage()->Remove(m_CalibNode); m_CalibNode = 0; this->GetDataStorage()->Remove(m_WorldNode); m_WorldNode = 0; m_Controls.m_ToolBox->setCurrentIndex(0); } void QmitkUltrasoundCalibration::OnDeviceServiceEvent(const ctkServiceEvent event) { if (m_CombinedModality.IsNull() || event.getType() != ctkServiceEvent::MODIFIED) { return; } ctkServiceReference service = event.getServiceReference(); QString curDepth = service.getProperty(QString::fromStdString(mitk::USDevice::GetPropertyKeys().US_PROPKEY_BMODE_DEPTH)).toString(); if (m_CurrentDepth != curDepth) { m_CurrentDepth = curDepth; this->OnReset(); } } void QmitkUltrasoundCalibration::OnAddCalibPoint() { mitk::Point3D world = this->GetRenderWindowPart()->GetSelectedPosition(); this->m_CalibPointsImage->InsertPoint(m_CalibPointsCount, world); this->m_CalibPointsTool->InsertPoint(m_CalibPointsCount, this->m_FreezePoint); QString text = text.number(m_CalibPointsCount + 1); text = "Point " + text; this->m_Controls.m_CalibPointList->addItem(text); m_CalibPointsCount++; SwitchFreeze(); } void QmitkUltrasoundCalibration::OnCalibration() { // Compute transformation vtkSmartPointer transform = vtkSmartPointer::New(); transform->SetSourceLandmarks(this->ConvertPointSetToVtkPolyData(m_CalibPointsImage)->GetPoints()); transform->SetTargetLandmarks(this->ConvertPointSetToVtkPolyData(m_CalibPointsTool)->GetPoints()); if( !m_CombinedModality->GetIsTrackedUltrasoundActive() ) { if (m_Controls.m_ScaleTransform->isChecked()) { transform->SetModeToSimilarity(); } //use affine transform else { transform->SetModeToRigidBody(); } //use similarity transform: scaling is not touched MITK_INFO << "TEST"; } else { transform->SetModeToRigidBody();//use similarity transform: scaling is not touched } transform->Modified(); transform->Update(); // Convert from vtk to itk data types itk::Matrix rotationFloat = itk::Matrix(); itk::Vector translationFloat = itk::Vector(); vtkSmartPointer m = transform->GetMatrix(); rotationFloat[0][0] = m->GetElement(0, 0); rotationFloat[0][1] = m->GetElement(0, 1); rotationFloat[0][2] = m->GetElement(0, 2); rotationFloat[1][0] = m->GetElement(1, 0); rotationFloat[1][1] = m->GetElement(1, 1); rotationFloat[1][2] = m->GetElement(1, 2); rotationFloat[2][0] = m->GetElement(2, 0); rotationFloat[2][1] = m->GetElement(2, 1); rotationFloat[2][2] = m->GetElement(2, 2); translationFloat[0] = m->GetElement(0, 3); translationFloat[1] = m->GetElement(1, 3); translationFloat[2] = m->GetElement(2, 3); mitk::DataNode::Pointer CalibPointsImage = mitk::DataNode::New(); CalibPointsImage->SetName("Calibration Points Image"); CalibPointsImage->SetData(m_CalibPointsImage); this->GetDataStorage()->Add(CalibPointsImage); mitk::DataNode::Pointer CalibPointsTracking = mitk::DataNode::New(); CalibPointsTracking->SetName("Calibration Points Tracking"); CalibPointsTracking->SetData(m_CalibPointsTool); this->GetDataStorage()->Add(CalibPointsTracking); mitk::PointSet::Pointer ImagePointsTransformed = m_CalibPointsImage->Clone(); this->ApplyTransformToPointSet(ImagePointsTransformed, transform); mitk::DataNode::Pointer CalibPointsImageTransformed = mitk::DataNode::New(); CalibPointsImageTransformed->SetName("Calibration Points Image (Transformed)"); CalibPointsImageTransformed->SetData(ImagePointsTransformed); this->GetDataStorage()->Add(CalibPointsImageTransformed); // Set output variable mitk::AffineTransform3D::Pointer oldUSImageTransform = m_CombinedModality->GetUltrasoundDevice()->GetOutput()->GetGeometry()->GetIndexToWorldTransform(); //including spacing! MITK_INFO << "Old US Image transform: " << oldUSImageTransform; mitk::AffineTransform3D::Pointer calibTransform = mitk::AffineTransform3D::New(); calibTransform->SetTranslation(translationFloat); calibTransform->SetMatrix(rotationFloat); MITK_INFO << "Calibration transform: " << calibTransform; m_Transformation = mitk::AffineTransform3D::New(); if( !m_CombinedModality->GetIsTrackedUltrasoundActive() ) { if( !m_Controls.m_ScaleTransform->isChecked() ) { m_Transformation->Compose(oldUSImageTransform); } MITK_INFO << "Used old USImageTransform"; } m_Transformation->Compose(calibTransform); MITK_INFO << "New combined transform: " << m_Transformation; mitk::SlicedGeometry3D::Pointer sliced3d = dynamic_cast (m_Node->GetData()->GetGeometry()); mitk::PlaneGeometry::Pointer plane = const_cast (sliced3d->GetPlaneGeometry(0)); plane->SetIndexToWorldTransform(m_Transformation); // Save to US-Device m_CombinedModality->SetCalibration(m_Transformation); m_Controls.m_ToolBox->setItemEnabled(2, true); // Save to NeedleProjectionFilter m_NeedleProjectionFilter->SetTargetPlane(m_Transformation); // Update Calibration FRE m_CalibrationStatistics = mitk::PointSetDifferenceStatisticsCalculator::New(); mitk::PointSet::Pointer p1 = this->m_CalibPointsTool->Clone(); // We use clones to calculate statistics to avoid concurrency Problems // Create point set with transformed image calibration points for // calculating the difference of image calibration and tool // calibration points in one geometry space mitk::PointSet::Pointer p2 = mitk::PointSet::New(); int n = 0; for (mitk::PointSet::PointsConstIterator it = m_CalibPointsImage->Begin(); it != m_CalibPointsImage->End(); ++it, ++n) { p2->InsertPoint(n, m_Transformation->TransformPoint(it->Value())); } m_CalibrationStatistics->SetPointSets(p1, p2); //QString text = text.number(m_CalibrationStatistics->GetRMS()); QString text = QString::number(ComputeFRE(m_CalibPointsImage, m_CalibPointsTool, transform)); MITK_INFO << "Calibration FRE: " << text.toStdString().c_str(); m_Controls.m_EvalLblCalibrationFRE->setText(text); m_Node->SetStringProperty("Calibration FRE", text.toStdString().c_str()); // Enable Button to save Calibration m_Controls.m_CalibBtnSaveCalibration->setEnabled(true); } void QmitkUltrasoundCalibration::OnAddEvalTargetPoint() { mitk::Point3D world = this->GetRenderWindowPart()->GetSelectedPosition(); this->m_EvalPointsImage->InsertPoint(m_EvalPointsImage->GetSize(), world); this->m_EvalPointsTool->InsertPoint(m_EvalPointsTool->GetSize(), this->m_FreezePoint); QString text = text.number(this->m_EvalPointsTool->GetSize()); this->m_Controls.m_EvalLblNumTargetPoints->setText(text); // Update FREs // Update Evaluation FRE, but only if it contains more than one point (will crash otherwise) if ((m_EvalPointsProjected->GetSize() > 1) && (m_EvalPointsTool->GetSize() > 1)) { m_EvaluationStatistics = mitk::PointSetDifferenceStatisticsCalculator::New(); m_ProjectionStatistics = mitk::PointSetDifferenceStatisticsCalculator::New(); mitk::PointSet::Pointer p1 = this->m_EvalPointsTool->Clone(); // We use clones to calculate statistics to avoid concurrency Problems mitk::PointSet::Pointer p2 = this->m_EvalPointsImage->Clone(); mitk::PointSet::Pointer p3 = this->m_EvalPointsProjected->Clone(); m_EvaluationStatistics->SetPointSets(p1, p2); m_ProjectionStatistics->SetPointSets(p1, p3); QString evalText = evalText.number(m_EvaluationStatistics->GetRMS()); QString projText = projText.number(m_ProjectionStatistics->GetRMS()); m_Controls.m_EvalLblEvaluationFRE->setText(evalText); m_Controls.m_EvalLblProjectionFRE->setText(projText); } SwitchFreeze(); } void QmitkUltrasoundCalibration::OnAddEvalProjectedPoint() { MITK_WARN << "Projection Evaluation may currently be inaccurate."; // TODO: Verify correct Evaluation. Is the Point that is added really current? mitk::Point3D projection = this->m_NeedleProjectionFilter->GetProjection()->GetPoint(1); m_EvalPointsProjected->InsertPoint(m_EvalPointsProjected->GetSize(), projection); QString text = text.number(this->m_EvalPointsProjected->GetSize()); this->m_Controls.m_EvalLblNumProjectionPoints->setText(text); } void QmitkUltrasoundCalibration::OnSaveEvaluation() { //Filename without suffix QString filename = m_Controls.m_EvalFilePath->text() + "//" + m_Controls.m_EvalFilePrefix->text(); MITK_WARN << "CANNOT SAVE, ABORTING!"; /* not working any more TODO! mitk::PointSetWriter::Pointer psWriter = mitk::PointSetWriter::New(); psWriter->SetInput(0, m_CalibPointsImage); psWriter->SetInput(1, m_CalibPointsTool); psWriter->SetInput(2, m_EvalPointsImage); psWriter->SetInput(3, m_EvalPointsTool); psWriter->SetInput(4, m_EvalPointsProjected); psWriter->SetFileName(filename.toStdString() + ".xml"); psWriter->Write(); */ // TODO: New writer for transformations must be implemented. /* mitk::TransformationFileWriter::Pointer tWriter = mitk::TransformationFileWriter::New(); tWriter->SetInput(0, m_CalibPointsImage); tWriter->SetInput(1, m_CalibPointsTool); tWriter->SetInput(2, m_EvalPointsImage); tWriter->SetInput(3, m_EvalPointsTool); tWriter->SetInput(4, m_EvalPointsProjected); tWriter->SetOutputFilename(filename.toStdString() + ".txt"); tWriter->DoWrite(this->m_Transformation); */ } void QmitkUltrasoundCalibration::OnSaveCalibration() { m_Controls.m_GotCalibrationLabel->setText(""); QString filename = QFileDialog::getSaveFileName(QApplication::activeWindow(), "Save Calibration", "", "Calibration files *.cal"); QFile file(filename); if (!file.open(QIODevice::WriteOnly | QIODevice::Text | QIODevice::Truncate)) { MITK_WARN << "Cannot open file '" << filename.toStdString() << "' for writing."; return; } std::string calibrationSerialization = m_CombinedModality->SerializeCalibration(); QTextStream outStream(&file); outStream << QString::fromStdString(calibrationSerialization); //save additional information if (m_Controls.m_saveAdditionalCalibrationLog->isChecked()) { mitk::SceneIO::Pointer mySceneIO = mitk::SceneIO::New(); QString filenameScene = filename + "_mitkScene.mitk"; mitk::NodePredicateNot::Pointer isNotHelperObject = mitk::NodePredicateNot::New(mitk::NodePredicateProperty::New("helper object", mitk::BoolProperty::New(true))); mitk::DataStorage::SetOfObjects::ConstPointer nodesToBeSaved = this->GetDataStorage()->GetSubset(isNotHelperObject); mySceneIO->SaveScene(nodesToBeSaved, this->GetDataStorage(), filenameScene.toStdString().c_str()); } } void QmitkUltrasoundCalibration::OnReset() { this->ClearTemporaryMembers(); if (m_Transformation.IsNull()) { m_Transformation = mitk::AffineTransform3D::New(); } m_Transformation->SetIdentity(); if (m_Node.IsNotNull() && (m_Node->GetData() != nullptr) && (m_Node->GetData()->GetGeometry() != nullptr)) { mitk::SlicedGeometry3D::Pointer sliced3d = dynamic_cast (m_Node->GetData()->GetGeometry()); mitk::PlaneGeometry::Pointer plane = const_cast (sliced3d->GetPlaneGeometry(0)); plane->SetIndexToWorldTransform(m_Transformation); } QString text1 = text1.number(this->m_EvalPointsTool->GetSize()); this->m_Controls.m_EvalLblNumTargetPoints->setText(text1); QString text2 = text2.number(this->m_EvalPointsProjected->GetSize()); this->m_Controls.m_EvalLblNumProjectionPoints->setText(text2); } void QmitkUltrasoundCalibration::Update() { //QList nodes = this->GetDataManagerSelection(); // if (nodes.empty()) return; // Update Tracking Data std::vector* datas = new std::vector(); datas->push_back(m_Tracker->GetOutput()); m_Controls.m_CalibTrackingStatus->SetNavigationDatas(datas); m_Controls.m_CalibTrackingStatus->Refresh(); m_Controls.m_EvalTrackingStatus->SetNavigationDatas(datas); m_Controls.m_EvalTrackingStatus->Refresh(); /* if (m_Image.IsNotNull() && m_Image->IsInitialized()) { m_Node->SetData(m_Image); } else { m_Image = m_CombinedModality->GetOutput(); m_Node->SetData(m_Image); }*/ m_CombinedModality->Modified(); m_CombinedModality->Update(); // Update US Image mitk::Image::Pointer image = m_CombinedModality->GetOutput(); // make sure that always the current image is set to the data node if (image.IsNotNull() && m_Node->GetData() != image.GetPointer() && image->IsInitialized()) { m_Node->SetData(image); } // Update Needle Projection m_NeedleProjectionFilter->Update(); //only update 2d window because it is faster this->RequestRenderWindowUpdate(mitk::RenderingManager::REQUEST_UPDATE_2DWINDOWS); } void QmitkUltrasoundCalibration::SwitchFreeze() { m_Controls.m_CalibBtnAddPoint->setEnabled(false); // generally deactivate // We use the activity state of the timer to determine whether we are currently viewing images if (!m_Timer->isActive()) // Activate Imaging { // if (m_Node) m_Node->ReleaseData(); if (m_CombinedModality.IsNull()){ m_Timer->stop(); return; } m_CombinedModality->Update(); m_Image = m_CombinedModality->GetOutput(); if (m_Image.IsNotNull() && m_Image->IsInitialized()) { m_Node->SetData(m_Image); } std::vector datas; datas.push_back(m_Tracker->GetOutput()); m_Controls.m_CalibTrackingStatus->SetNavigationDatas(&datas); m_Controls.m_CalibTrackingStatus->ShowStatusLabels(); m_Controls.m_CalibTrackingStatus->Refresh(); m_Controls.m_EvalTrackingStatus->SetNavigationDatas(&datas); m_Controls.m_EvalTrackingStatus->ShowStatusLabels(); m_Controls.m_EvalTrackingStatus->Refresh(); int interval = 40; m_Timer->setInterval(interval); m_Timer->start(); m_CombinedModality->SetIsFreezed(false); } else if (this->m_Tracker->GetOutput(0)->IsDataValid()) { //deactivate Imaging m_Timer->stop(); // Remember last tracking coordinates m_FreezePoint = this->m_Tracker->GetOutput(0)->GetPosition(); m_Controls.m_CalibBtnAddPoint->setEnabled(true); // activate only, if valid point is set m_CombinedModality->SetIsFreezed(true); } } void QmitkUltrasoundCalibration::ShowNeedlePath() { // Init Filter this->m_NeedleProjectionFilter->SelectInput(0); // Create Node for Pointset mitk::DataNode::Pointer node = this->GetDataStorage()->GetNamedNode("Needle Path"); if (node.IsNull()) { node = mitk::DataNode::New(); node->SetName("Needle Path"); node->SetData(m_NeedleProjectionFilter->GetProjection()); node->SetBoolProperty("show contour", true); this->GetDataStorage()->Add(node); } } void QmitkUltrasoundCalibration::ClearTemporaryMembers() { m_CalibPointsTool->Clear(); m_CalibPointsImage->Clear(); m_CalibPointsCount = 0; m_EvalPointsImage->Clear(); m_EvalPointsTool->Clear(); m_EvalPointsProjected->Clear(); this->m_Controls.m_CalibPointList->clear(); m_SpacingPoints->Clear(); m_Controls.m_SpacingPointsList->clear(); m_SpacingPointsCount = 0; } vtkSmartPointer QmitkUltrasoundCalibration::ConvertPointSetToVtkPolyData(mitk::PointSet::Pointer PointSet) { vtkSmartPointer returnValue = vtkSmartPointer::New(); vtkSmartPointer points = vtkSmartPointer::New(); for (int i = 0; i < PointSet->GetSize(); i++) { double point[3] = { PointSet->GetPoint(i)[0], PointSet->GetPoint(i)[1], PointSet->GetPoint(i)[2] }; points->InsertNextPoint(point); } vtkSmartPointer temp = vtkSmartPointer::New(); temp->SetPoints(points); vtkSmartPointer vertexFilter = vtkSmartPointer::New(); vertexFilter->SetInputData(temp); vertexFilter->Update(); returnValue->ShallowCopy(vertexFilter->GetOutput()); return returnValue; } double QmitkUltrasoundCalibration::ComputeFRE(mitk::PointSet::Pointer imageFiducials, mitk::PointSet::Pointer realWorldFiducials, vtkSmartPointer transform) { if (imageFiducials->GetSize() != realWorldFiducials->GetSize()) return -1; double FRE = 0; for (int i = 0; i < imageFiducials->GetSize(); ++i) { itk::Point current_image_fiducial_point = imageFiducials->GetPoint(i); if (transform != nullptr) { current_image_fiducial_point = transform->TransformPoint(imageFiducials->GetPoint(i)[0], imageFiducials->GetPoint(i)[1], imageFiducials->GetPoint(i)[2]); } double cur_error_squared = current_image_fiducial_point.SquaredEuclideanDistanceTo(realWorldFiducials->GetPoint(i)); FRE += cur_error_squared; } FRE = sqrt(FRE / (double)imageFiducials->GetSize()); return FRE; } void QmitkUltrasoundCalibration::ApplyTransformToPointSet(mitk::PointSet::Pointer pointSet, vtkSmartPointer transform) { for (int i = 0; i < pointSet->GetSize(); ++i) { itk::Point current_point_transformed = itk::Point(); current_point_transformed = transform->TransformPoint(pointSet->GetPoint(i)[0], pointSet->GetPoint(i)[1], pointSet->GetPoint(i)[2]); pointSet->SetPoint(i, current_point_transformed); } } void QmitkUltrasoundCalibration::OnFreezeClicked() { if (m_CombinedModality->GetIsFreezed()) { //device was already frozen so we need to delete all Spacing points because they need to be collected all at once // no need to check if all four points are already collected, because if thats the case you can no longer click the Freeze Button m_SpacingPoints->Clear(); m_Controls.m_SpacingPointsList->clear(); m_SpacingPointsCount = 0; m_Controls.m_SpacingAddPoint->setEnabled(false); } else { m_Controls.m_SpacingAddPoint->setEnabled(true); } SwitchFreeze(); } void QmitkUltrasoundCalibration::OnAddSpacingPoint() { mitk::Point3D point = this->GetRenderWindowPart()->GetSelectedPosition(); this->m_SpacingPoints->InsertPoint(m_SpacingPointsCount, point); QString text = text.number(m_SpacingPointsCount + 1); text = "Point " + text; this->m_Controls.m_SpacingPointsList->addItem(text); m_SpacingPointsCount++; if (m_SpacingPointsCount == 4) //now we have all 4 points needed { m_Controls.m_SpacingAddPoint->setEnabled(false); m_Controls.m_CalculateSpacing->setEnabled(true); m_Controls.m_SpacingBtnFreeze->setEnabled(false); } } void QmitkUltrasoundCalibration::OnCalculateSpacing() { mitk::Point3D horizontalOne = m_SpacingPoints->GetPoint(0); mitk::Point3D horizontalTwo = m_SpacingPoints->GetPoint(1); mitk::Point3D verticalOne = m_SpacingPoints->GetPoint(2); mitk::Point3D verticalTwo = m_SpacingPoints->GetPoint(3); //Get the distances between the points in the image double xDistance = horizontalOne.EuclideanDistanceTo(horizontalTwo); double yDistance = verticalOne.EuclideanDistanceTo(verticalTwo); //Calculate the spacing of the image and fill a vector with it double xSpacing = 30 / xDistance; double ySpacing = 20 / yDistance; m_CombinedModality->GetUltrasoundDevice()->SetSpacing(xSpacing, ySpacing); //Now that the spacing is set clear all stuff and return to Calibration m_SpacingPoints->Clear(); m_Controls.m_SpacingPointsList->clear(); m_SpacingPointsCount = 0; m_CombinedModality->SetIsFreezed(false); } void QmitkUltrasoundCalibration::OnUSDepthChanged(const std::string& key, const std::string&) { //whenever depth of USImage is changed the spacing should no longer be overwritten if (key == mitk::USDevice::GetPropertyKeys().US_PROPKEY_BMODE_DEPTH) { } }