diff --git a/Plugins/PluginList.cmake b/Plugins/PluginList.cmake index 5cf41c4112..399f32ed0f 100644 --- a/Plugins/PluginList.cmake +++ b/Plugins/PluginList.cmake @@ -1,95 +1,96 @@ # Plug-ins must be ordered according to their dependencies set(MITK_PLUGINS org.blueberry.core.runtime:ON org.blueberry.core.expressions:OFF org.blueberry.core.commands:OFF org.blueberry.core.jobs:OFF org.blueberry.ui.qt:OFF org.blueberry.ui.qt.help:ON org.blueberry.ui.qt.log:ON org.blueberry.ui.qt.objectinspector:OFF org.mitk.core.services:ON org.mitk.gui.common:ON org.mitk.planarfigure:ON org.mitk.core.ext:OFF org.mitk.core.jobs:OFF org.mitk.gui.qt.application:ON org.mitk.gui.qt.ext:OFF org.mitk.gui.qt.extapplication:OFF org.mitk.gui.qt.common:ON org.mitk.gui.qt.stdmultiwidgeteditor:ON org.mitk.gui.qt.mxnmultiwidgeteditor:OFF org.mitk.gui.qt.common.legacy:OFF org.mitk.gui.qt.cmdlinemodules:OFF org.mitk.gui.qt.chartExample:OFF org.mitk.gui.qt.datamanager:ON org.mitk.gui.qt.datamanagerlight:OFF org.mitk.gui.qt.datastorageviewertest:OFF org.mitk.gui.qt.properties:ON org.mitk.gui.qt.basicimageprocessing:OFF org.mitk.gui.qt.dicom:OFF org.mitk.gui.qt.dicominspector:OFF org.mitk.gui.qt.dosevisualization:OFF org.mitk.gui.qt.geometrytools:OFF org.mitk.gui.qt.igtexamples:OFF org.mitk.gui.qt.igttracking:OFF org.mitk.gui.qt.lasercontrol:OFF org.mitk.gui.qt.openigtlink:OFF org.mitk.gui.qt.imagecropper:OFF org.mitk.gui.qt.imagenavigator:ON org.mitk.gui.qt.viewnavigator:OFF org.mitk.gui.qt.materialeditor:OFF org.mitk.gui.qt.measurementtoolbox:OFF org.mitk.gui.qt.moviemaker:OFF org.mitk.gui.qt.pointsetinteraction:OFF org.mitk.gui.qt.pointsetinteractionmultispectrum:OFF org.mitk.gui.qt.python:OFF org.mitk.gui.qt.remeshing:OFF org.mitk.gui.qt.segmentation:OFF org.mitk.gui.qt.aicpregistration:OFF org.mitk.gui.qt.renderwindowmanager:OFF org.mitk.gui.qt.semanticrelations:OFF org.mitk.gui.qt.toftutorial:OFF org.mitk.gui.qt.tofutil:OFF org.mitk.gui.qt.tubegraph:OFF org.mitk.gui.qt.ugvisualization:OFF org.mitk.gui.qt.photoacoustics.pausviewer:OFF org.mitk.gui.qt.photoacoustics.pausmotioncompensation:OFF org.mitk.gui.qt.photoacoustics.imageprocessing:OFF org.mitk.gui.qt.photoacoustics.simulation:OFF org.mitk.gui.qt.photoacoustics.spectralunmixing:OFF org.mitk.gui.qt.ultrasound:OFF org.mitk.gui.qt.volumevisualization:OFF org.mitk.gui.qt.eventrecorder:OFF org.mitk.gui.qt.xnat:OFF org.mitk.gui.qt.igt.app.ultrasoundtrackingnavigation:OFF org.mitk.gui.qt.spectrocamrecorder:OFF org.mitk.gui.qt.classificationsegmentation:OFF org.mitk.gui.qt.overlaymanager:OFF org.mitk.gui.qt.igt.app.hummelprotocolmeasurements:OFF org.mitk.gui.qt.multilabelsegmentation:OFF org.mitk.matchpoint.core.helper:OFF org.mitk.gui.qt.matchpoint.algorithm.browser:OFF org.mitk.gui.qt.matchpoint.algorithm.control:OFF org.mitk.gui.qt.matchpoint.mapper:OFF org.mitk.gui.qt.matchpoint.framereg:OFF org.mitk.gui.qt.matchpoint.visualizer:OFF org.mitk.gui.qt.matchpoint.evaluator:OFF org.mitk.gui.qt.matchpoint.manipulator:OFF org.mitk.gui.qt.preprocessing.resampling:OFF org.mitk.gui.qt.radiomics:OFF org.mitk.gui.qt.cest:OFF org.mitk.gui.qt.fit.demo:OFF org.mitk.gui.qt.fit.inspector:OFF org.mitk.gui.qt.fit.genericfitting:OFF org.mitk.gui.qt.pharmacokinetics.mri:OFF org.mitk.gui.qt.pharmacokinetics.pet:OFF org.mitk.gui.qt.pharmacokinetics.simulation:OFF org.mitk.gui.qt.pharmacokinetics.curvedescriptor:OFF org.mitk.gui.qt.pharmacokinetics.concentration.mri:OFF org.mitk.gui.qt.flowapplication:OFF org.mitk.gui.qt.flow.segmentation:OFF + org.mitk.gui.qt.standardplanetool:ON ) diff --git a/Plugins/org.mitk.gui.qt.standardplanetool/CMakeLists.txt b/Plugins/org.mitk.gui.qt.standardplanetool/CMakeLists.txt new file mode 100644 index 0000000000..4fe5f72c1d --- /dev/null +++ b/Plugins/org.mitk.gui.qt.standardplanetool/CMakeLists.txt @@ -0,0 +1,9 @@ +project(org_mitk_gui_qt_standardplanetool) + +mitk_create_plugin( + EXPORT_DIRECTIVE MITK_QT_STANDARDPLANETOOL + EXPORTED_INCLUDE_SUFFIXES src + MODULE_DEPENDS MitkRemeshing MitkSegmentation MitkQtWidgetsExt MitkOpenCVVideoSupport MitkSceneSerialization MitkSurfaceInterpolation + PACKAGE_DEPENDS PUBLIC OpenCV VTK ITK +) + diff --git a/Plugins/org.mitk.gui.qt.standardplanetool/documentation/UserManual/StandardPlaneTool.dox b/Plugins/org.mitk.gui.qt.standardplanetool/documentation/UserManual/StandardPlaneTool.dox new file mode 100644 index 0000000000..c93bf54ca8 --- /dev/null +++ b/Plugins/org.mitk.gui.qt.standardplanetool/documentation/UserManual/StandardPlaneTool.dox @@ -0,0 +1,4 @@ +/** +\page org_mitk_gui_qt_standardplanetool Standard Plane Tool + +*/ \ No newline at end of file diff --git a/Plugins/org.mitk.gui.qt.standardplanetool/documentation/doxygen/modules.dox b/Plugins/org.mitk.gui.qt.standardplanetool/documentation/doxygen/modules.dox new file mode 100644 index 0000000000..83ddb2d2a1 --- /dev/null +++ b/Plugins/org.mitk.gui.qt.standardplanetool/documentation/doxygen/modules.dox @@ -0,0 +1,16 @@ +/** + \defgroup org_mitk_gui_qt_imagecropper org.mitk.gui.qt.imagecropper + \ingroup MITKPlugins + + \brief This is the image cropper plugin. It can crop images. + +*/ + +/** + \defgroup org_mitk_gui_qt_imagecropper_internal Internal + \ingroup org_mitk_gui_qt_imagecropper + + \brief This subcategory includes the internal classes of the org.mitk.gui.qt.imagecropper plugin. Other + plugins must not rely on these classes. They contain implementation details and their interface + may change at any time. We mean it. +*/ \ No newline at end of file diff --git a/Plugins/org.mitk.gui.qt.standardplanetool/files.cmake b/Plugins/org.mitk.gui.qt.standardplanetool/files.cmake new file mode 100644 index 0000000000..1102945e3f --- /dev/null +++ b/Plugins/org.mitk.gui.qt.standardplanetool/files.cmake @@ -0,0 +1,42 @@ +set(SRC_CPP_FILES + +) + +set(INTERNAL_CPP_FILES + org_mitk_gui_qt_standardplanetool_Activator.cpp + QmitkStandardPlaneTool.cpp + mitkShapeContourCostFunction.cpp + mitk2D3DShapeRegistration.cpp + vtkSilhouette.cpp +) + +set(UI_FILES + src/internal/QmitkStandardPlaneToolControls.ui +) + +set(MOC_H_FILES + src/internal/org_mitk_gui_qt_standardplanetool_Activator.h + src/internal/QmitkStandardPlaneTool.h + src/internal/mitkShapeContourCostFunction.h + src/internal/vtkSilhouette.h + src/internal/mitk2D3DShapeRegistration.h +) + +set(CACHED_RESOURCE_FILES + resources/icon.png + plugin.xml +) + +#set(QRC_FILES +# resources/imagecropper.qrc +#) + +set(CPP_FILES) + +foreach(file ${SRC_CPP_FILES}) + set(CPP_FILES ${CPP_FILES} src/${file}) +endforeach(file ${SRC_CPP_FILES}) + +foreach(file ${INTERNAL_CPP_FILES}) + set(CPP_FILES ${CPP_FILES} src/internal/${file}) +endforeach(file ${INTERNAL_CPP_FILES}) diff --git a/Plugins/org.mitk.gui.qt.standardplanetool/manifest_headers.cmake b/Plugins/org.mitk.gui.qt.standardplanetool/manifest_headers.cmake new file mode 100644 index 0000000000..a6f7e48634 --- /dev/null +++ b/Plugins/org.mitk.gui.qt.standardplanetool/manifest_headers.cmake @@ -0,0 +1,5 @@ +set(Plugin-Name "MITK Standard Planes") +set(Plugin-Version "1.0.0") +set(Plugin-Vendor "DKFZ, Medical and Biological Informatics") +set(Plugin-ContactAddress "http://www.mitk.org") +set(Require-Plugin org.mitk.gui.qt.common) \ No newline at end of file diff --git a/Plugins/org.mitk.gui.qt.standardplanetool/plugin.xml b/Plugins/org.mitk.gui.qt.standardplanetool/plugin.xml new file mode 100644 index 0000000000..892629096d --- /dev/null +++ b/Plugins/org.mitk.gui.qt.standardplanetool/plugin.xml @@ -0,0 +1,21 @@ + + + + + + + + + Crop images to a given size + + + + + diff --git a/Plugins/org.mitk.gui.qt.standardplanetool/resources/icon.png b/Plugins/org.mitk.gui.qt.standardplanetool/resources/icon.png new file mode 100644 index 0000000000..0fdfa0dfd1 Binary files /dev/null and b/Plugins/org.mitk.gui.qt.standardplanetool/resources/icon.png differ diff --git a/Plugins/org.mitk.gui.qt.standardplanetool/src/internal/2d3dRegistration.h b/Plugins/org.mitk.gui.qt.standardplanetool/src/internal/2d3dRegistration.h new file mode 100644 index 0000000000..545ea8e038 --- /dev/null +++ b/Plugins/org.mitk.gui.qt.standardplanetool/src/internal/2d3dRegistration.h @@ -0,0 +1,108 @@ +/*=================================================================== + +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. + +===================================================================*/ +#ifndef mitk2D3DShapeRegistration_h +#define mitk2D3DShapeRegistration_h + + + +// The following lines instantiate the evolutionary optimizer. + +typedef itk::OnePlusOneEvolutionaryOptimizer OptimizerType; +OptimizerType::Pointer optimizer = OptimizerType::New(); + +// Evolutionary algorithms are based on testing random variations +// of parameters. In order to support the computation of random values, +// ITK provides a family of random number generators. In this example, we +// use the \doxygen{NormalVariateGenerator} which generates values with a +// normal distribution. + +itk::Statistics::NormalVariateGenerator::Pointer generator += itk::Statistics::NormalVariateGenerator::New(); + +// The random number generator must be initialized with a seed. + +generator->Initialize(1234567); + +// The OnePlusOneEvolutionaryOptimizer is initialized by +// specifying the random number generator, the number of samples for the +// initial population and the maximum number of iterations. + +optimizer->SetNormalVariateGenerator(generator); +optimizer->SetCostFunction(costFunction); +optimizer->Initialize(1000); +optimizer->SetInitialRadius(0.5); +optimizer->SetMaximumIteration(10000); + +// scale the components appropriately +itk::Array parametersScale; +parametersScale.SetSize(6); +for (unsigned int i = 0; i<2; i++) +{ + parametersScale[i] = houghScaling; // angle scale +} +MeasurementVectorType::ComponentType maxMV = std::max(mv[2], mv[3]); +parametersScale[2] = 1.; // / maxMV; +parametersScale[3] = 1.; // / maxMV; + +for (unsigned int i = 4; i<6; i++) +{ + parametersScale[i] = 256; // offset scale +} +optimizer->SetScales(parametersScale); + +// Here we instantiate the Command object that will act as an +// observer of the registration method and print out parameters at each +// iteration. Earlier, we defined this command as a class templated over the +// optimizer type. Once it is created with the \code{New()} method, we +// connect the optimizer to the command. + +typedef IterationCallback< OptimizerType > IterationCallbackType; +IterationCallbackType::Pointer callback = IterationCallbackType::New(); +callback->SetOptimizer(optimizer); + +OptimizerType::ParametersType params; +params.SetSize(6); +for (int i = 0; i<4; ++i) + params[i] = mv[i]; +params[4] = minClCand; +params[5] = maxClCand; +optimizer->SetInitialPosition(params); + +std::cout << "Initial Parameters : " << params << std::endl; + +// make optimizer maximize +optimizer->MaximizeOn(); + +// do registration + +try +{ + optimizer->StartOptimization(); + std::cout << "Optimizer stop condition: " + << optimizer->GetStopConditionDescription() + << std::endl; +} +catch (itk::ExceptionObject & exp) +{ + std::cerr << "Exception caught ! " << std::endl; + std::cerr << exp << std::endl; +} + +// get registration result +OptimizerType::ParametersType finalParameters = +optimizer->GetCurrentPosition(); + +std::cout << "Final Solution is : " << finalParameters << std::endl; \ No newline at end of file diff --git a/Plugins/org.mitk.gui.qt.standardplanetool/src/internal/QmitkStandardPlaneTool.cpp b/Plugins/org.mitk.gui.qt.standardplanetool/src/internal/QmitkStandardPlaneTool.cpp new file mode 100644 index 0000000000..3a489b5188 --- /dev/null +++ b/Plugins/org.mitk.gui.qt.standardplanetool/src/internal/QmitkStandardPlaneTool.cpp @@ -0,0 +1,2712 @@ + +/*=================================================================== + +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 "QmitkStandardPlaneTool.h" + +#include "Poco/DirectoryIterator.h" + +#include +#include +#include + +// MITK +#include "mitk2D3DShapeRegistration.h" +#include "mitkConvert2Dto3DImageFilter.h" +#include "mitkImageSliceSelector.h" +#include "mitkNodePredicateDataType.h" +#include "mitkNodePredicateProperty.h" +#include "mitkPlanarCircle.h" +#include "mitkPlanePositionManager.h" +#include "mitkVtkRepresentationProperty.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "org_mitk_gui_qt_standardplanetool_Activator.h" +#include "vtkPoints.h" +#include "vtkSilhouette.h" + +#include + +// Qmitk +#include "QmitkRenderWindow.h" +#include "QmitkStdMultiWidget.h" + +// Qt +#include +#include +#include + +// VTK (for testing evaluation methods) +#include "vtkDoubleArray.h" +#include "vtkFloatArray.h" +#include "vtkPointData.h" +#include "vtkProperty.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include "itkImage.h" +#include "itkOpenCVImageBridge.h" +#include "itkRescaleIntensityImageFilter.h" +#include +#include +#include +#include + +// Boost +#include +#include + +// Poco +#include +#include +#include +#include +#include +#include +#include + +#include + +const std::string QmitkStandardPlaneTool::VIEW_ID = "org.mitk.views.qmitkstandardplanetool"; + +bool ReadPlaneOrientationsFromFile(std::string plane_file_name, + std::string &file_name, + mitk::PlaneGeometry::Pointer &ax, + mitk::PlaneGeometry::Pointer &sag, + mitk::PlaneGeometry::Pointer &cor) +{ + std::vector> result; + Poco::Path filePath(plane_file_name); + + Poco::File inputFile(filePath); + + if (inputFile.exists()) + { + Poco::FileInputStream fis(filePath.toString()); + Poco::JSON::Parser parser; + Poco::Dynamic::Var result; + Poco::JSON::Object::Ptr json_object; + + parser.parse(fis); + result = parser.result(); + json_object = result.extract(); + MITK_INFO << json_object->getObject("FileName"); + Poco::Dynamic::Var v = json_object->get("MPR.Axial: Current.*"); + + // Extract the JSON Object + Poco::DynamicStruct jsonStruct = *result.extract(); + + Poco::JSON::Array::ConstIterator it; + Poco::JSON::Array::Ptr array; + + // auto axialPlane = mitk::PlaneGeometry::New(); + // auto sagittalPlane = mitk::PlaneGeometry::New(); + // auto coronalPlane = mitk::PlaneGeometry::New(); + + for (int i = 0; i < 4; i++) + for (int j = 0; j < 4; j++) + { + jsonStruct["MPR.Axial: Current.*"]["PlaneMatrix"][i][j]; + jsonStruct["MPR.Sagittal: Current.*"]["PlaneMatrix"][i][j]; + jsonStruct["MPR.Coronal: Current.*"]["PlaneMatrix"][i][j]; + ax->GetVtkMatrix()->SetElement(i, j, jsonStruct["MPR.Axial: Current.*"]["PlaneMatrix"][i][j]); + sag->GetVtkMatrix()->SetElement(i, j, jsonStruct["MPR.Sagittal: Current.*"]["PlaneMatrix"][i][j]); + cor->GetVtkMatrix()->SetElement(i, j, jsonStruct["MPR.Coronal: Current.*"]["PlaneMatrix"][i][j]); + } + + // std::string axialString = jsonStruct["MPR.Axial: Current.*"]["PlaneMatrix"][0].toString(); + // std::string sagittalString = jsonStruct["MPR.Sagittal: Current.*"]["PlaneMatrix"].toString(); + // std::string coronalString = jsonStruct["MPR.Coronal: Current.*"]["PlaneMatrix"].toString(); + + MITK_INFO << *ax->GetVtkMatrix() << " " + << jsonStruct["MPR.Axial: Current.*"]["PlaneMatrix"].toString(); // << " " << axialString; + MITK_INFO << jsonStruct["FileName"].toString(); + + size_t pos = std::string::npos; + std::string toErase = "D:\\Ankle_dataset\\"; + file_name = jsonStruct["FileName"].toString(); + + // Search for the substring in string in a loop untill nothing is found + while ((pos = file_name.find(toErase)) != std::string::npos) + { + // If found then erase it from string + file_name.erase(pos, toErase.length()); + } + + // return result; + } + return 0; +} + +QmitkStandardPlaneTool::QmitkStandardPlaneTool(QObject *parent) + : m_OldNode(nullptr), + m_WorkingNode(nullptr), + m_Projection1(nullptr), + m_Projection2(nullptr), + m_Silhouette1(nullptr), + m_Silhouette2(nullptr), + m_Segmentation(nullptr), + m_Surface(nullptr), + m_showSPV(0), + m_showPEV(0), + m_showSV(0), + m_CurrentMortiseID(-1), + m_CurrentLateralID(-1), + m_CurrentAPID(-1) +{ +} + +QmitkStandardPlaneTool::~QmitkStandardPlaneTool() {} + +void QmitkStandardPlaneTool::CreateQtPartControl(QWidget *parent) +{ + // create GUI widgets from the Qt Designer's .ui file + m_Controls.setupUi(parent); + + ////invalidate timer for first measurement of plane-adjustment-time (user tab) + // m_timer.invalidate(); + + // create GUI widgets from the Qt Designer's .ui file + m_Controls.patientImageSelector->SetDataStorage(this->GetDataStorage()); + m_Controls.patientImageSelector->SetPredicate( + mitk::NodePredicateAnd::New(mitk::TNodePredicateDataType::New(), + mitk::NodePredicateNot::New(mitk::NodePredicateProperty::New("helper object")))); + m_WorkingNode = m_Controls.patientImageSelector->GetSelectedNode(); + + //// create signal/slot connections of comboboxes + connect(m_Controls.patientImageSelector, + SIGNAL(OnSelectionChanged(const mitk::DataNode *)), + this, + SLOT(OnComboBoxSelectionChanged(const mitk::DataNode *))); + + connect(m_Controls.pB_StandardPlaneTool, SIGNAL(clicked()), this, SLOT(ToggleStandardPlaneToolView())); + connect(m_Controls.pB_ProjectionExtractionTool, SIGNAL(clicked()), this, SLOT(ToggleProjectionExtractionToolView())); + connect(m_Controls.pB_SilhouetteTool, SIGNAL(clicked()), this, SLOT(ToggleSilhouetteToolView())); + + m_Controls.patientImageSelector->setEnabled(true); + + this->InitProjectionExtractionToolWidget(); + this->InitStandardPlaneToolWidget(); + this->InitSilhouetteToolWidget(); +} + +void QmitkStandardPlaneTool::InitSilhouetteToolWidget() +{ + connect(m_Controls.pB_ExtractSilhouettes, SIGNAL(clicked()), this, SLOT(OnExtractSilhouetteClicked())); + connect(m_Controls.pB_StartRegistration, SIGNAL(clicked()), this, SLOT(OnRegistrationClicked())); + connect(m_Controls.pB_StartManualRegistration, SIGNAL(clicked()), this, SLOT(OnManualRegistrationClicked())); + connect(m_Controls.pB_ExtractContourAndSlices, SIGNAL(clicked()), this, SLOT(OnExtractContourAndSlicesClicked())); + + m_Controls.sB_Tx->setRange(0, 100); + m_Controls.sB_Ty->setRange(0, 100); + m_Controls.sB_Tz->setRange(0, 100); + m_Controls.sB_alpha->setRange(0, 100); + m_Controls.sB_beta->setRange(0, 100); + m_Controls.sB_gamma->setRange(0, 100); + + connect(m_Controls.sB_Tx, SIGNAL(sliderMoved(int)), this, SLOT(OnSliderValueChanged())); + connect(m_Controls.sB_Ty, SIGNAL(sliderMoved(int)), this, SLOT(OnSliderValueChanged())); + connect(m_Controls.sB_Tz, SIGNAL(sliderMoved(int)), this, SLOT(OnSliderValueChanged())); + connect(m_Controls.sB_alpha, SIGNAL(sliderMoved(int)), this, SLOT(OnSliderValueChanged())); + connect(m_Controls.sB_beta, SIGNAL(sliderMoved(int)), this, SLOT(OnSliderValueChanged())); + connect(m_Controls.sB_gamma, SIGNAL(sliderMoved(int)), this, SLOT(OnSliderValueChanged())); + + m_Controls.gB_SilhouetteTool->hide(); +} + +void QmitkStandardPlaneTool::OnSaveLateralSliceClicked() +{ + QmitkRenderWindow *RenderWindowAxial = this->GetRenderWindowPart()->GetQmitkRenderWindow("axial"); + m_CurrentLateralID = RenderWindowAxial->GetSliceNavigationController()->GetSlice()->GetSteps() - + RenderWindowAxial->GetSliceNavigationController()->GetSlice()->GetPos() - 1; + mitk::DataNode *imageNode(m_Controls.patientImageSelector->GetSelectedNode()); + if (imageNode == nullptr) + return; + imageNode->SetIntProperty("StandardProjection.Lateral", m_CurrentLateralID); + imageNode->Modified(); + this->m_Controls.lcd_L->display(m_CurrentLateralID); +} +void QmitkStandardPlaneTool::OnSaveMortiseSliceClicked() +{ + QmitkRenderWindow *RenderWindowAxial = this->GetRenderWindowPart()->GetQmitkRenderWindow("axial"); + m_CurrentMortiseID = RenderWindowAxial->GetSliceNavigationController()->GetSlice()->GetSteps() - + RenderWindowAxial->GetSliceNavigationController()->GetSlice()->GetPos() - 1; + mitk::DataNode *imageNode(m_Controls.patientImageSelector->GetSelectedNode()); + if (imageNode == nullptr) + return; + imageNode->SetIntProperty("StandardProjection.Mortise", m_CurrentMortiseID); + imageNode->Modified(); + this->m_Controls.lcd_M->display(m_CurrentMortiseID); +} +void QmitkStandardPlaneTool::OnSaveAPSliceClicked() +{ + QmitkRenderWindow *RenderWindowAxial = this->GetRenderWindowPart()->GetQmitkRenderWindow("axial"); + m_CurrentAPID = RenderWindowAxial->GetSliceNavigationController()->GetSlice()->GetSteps() - + RenderWindowAxial->GetSliceNavigationController()->GetSlice()->GetPos() - 1; + mitk::DataNode *imageNode(m_Controls.patientImageSelector->GetSelectedNode()); + if (imageNode == nullptr) + return; + imageNode->SetIntProperty("StandardProjection.AP", m_CurrentAPID); + imageNode->Modified(); + this->m_Controls.lcd_AP->display(m_CurrentAPID); +} + +void QmitkStandardPlaneTool::OnShowLateralSliceClicked() +{ + if (m_CurrentLateralID >= 0) + this->ShowAxialSlice(m_CurrentLateralID); +} +void QmitkStandardPlaneTool::OnShowMortiseSliceClicked() +{ + if (m_CurrentMortiseID >= 0) + this->ShowAxialSlice(m_CurrentMortiseID); +} +void QmitkStandardPlaneTool::OnShowAPSliceClicked() +{ + if (m_CurrentAPID >= 0) + this->ShowAxialSlice(m_CurrentAPID); +} + +void QmitkStandardPlaneTool::OnLoadSavedProjections() +{ + mitk::DataNode *imageNode(m_Controls.patientImageSelector->GetSelectedNode()); + if (imageNode == nullptr) + return; + + if (imageNode->GetIntProperty("StandardProjection.AP", m_CurrentAPID)) + { + imageNode->GetIntProperty("StandardProjection.AP", m_CurrentAPID); + this->m_Controls.lcd_AP->display(m_CurrentAPID); + } + else + { + this->m_Controls.lcd_AP->display(0); + m_CurrentAPID = -1; + } + + if (imageNode->GetIntProperty("StandardProjection.Lateral", m_CurrentLateralID)) + { + imageNode->GetIntProperty("StandardProjection.Lateral", m_CurrentLateralID); + this->m_Controls.lcd_L->display(m_CurrentLateralID); + } + else + { + this->m_Controls.lcd_L->display(0); + m_CurrentLateralID = -1; + } + if (imageNode->GetIntProperty("StandardProjection.Mortise", m_CurrentMortiseID)) + { + imageNode->GetIntProperty("StandardProjection.Mortise", m_CurrentMortiseID); + this->m_Controls.lcd_M->display(m_CurrentMortiseID); + } + else + { + m_CurrentMortiseID = -1; + this->m_Controls.lcd_M->display(0); + } +} + +void ExtractSlice(mitk::DataStorage *dataStorage, + mitk::Image::Pointer image, + int slice, + QString fileName, + bool setFirst, + bool save, + bool show = 0) +{ + MITK_INFO << slice; + auto selector = mitk::ImageSliceSelector::New(); + selector->SetSliceNr(slice); + selector->SetInput(image); + selector->Update(); + mitk::Image::Pointer projection1 = selector->GetOutput(); + + if (setFirst) + { + mitk::Point3D origin = projection1->GetGeometry()->GetOrigin(); + origin[2] = 0.0; + projection1->GetGeometry()->SetOrigin(origin[0]); + projection1->Modified(); + } + + auto convert2DTo3DImageFilter = mitk::Convert2Dto3DImageFilter::New(); + convert2DTo3DImageFilter->SetInput(projection1); + convert2DTo3DImageFilter->Update(); + + if (show) + { + auto projection = mitk::DataNode::New(); + projection->SetProperty("name", mitk::StringProperty::New("ImageSource")); + projection->SetData(projection1); + projection->SetColor(1.0, 1.0, 1.0); + projection->SetVisibility(false); + projection->SetName(fileName.toStdString()); + dataStorage->Add(projection); + } + + if (save) + { + if (fileName == nullptr) + return; + mitk::IOUtil::Save(convert2DTo3DImageFilter->GetOutput(), fileName.toStdString() + ".nrrd"); + MITK_INFO << "Image Saved"; + } +} + +void QmitkStandardPlaneTool::OnExportProjections() +{ + mitk::DataNode *imageNode(m_Controls.patientImageSelector->GetSelectedNode()); + if (imageNode == nullptr) + return; + + mitk::Image::Pointer image = dynamic_cast(imageNode->GetData()); + + QString workingDir("C:/Users/thomass/Desktop/StandardProjektionen/"); + QString LateralFileName, MortiseFileName, APFileName; + + ofstream myfile; + myfile.open(workingDir.append(imageNode->GetName().c_str()).append(".csv").toStdString()); + myfile << m_CurrentAPID << "," << m_CurrentLateralID << "," << m_CurrentMortiseID; + myfile.close(); + + // ExtractSlice(this->GetDataStorage(), image, m_CurrentAPID, + // APFileName.append(workingDir).append(imageNode->GetName().c_str()).append("_").append(QString::number(m_CurrentAPID)).append("_AP"), + // true, true); ExtractSlice(this->GetDataStorage(), image, m_CurrentLateralID, + // LateralFileName.append(workingDir).append(imageNode->GetName().c_str()).append("_").append(QString::number(m_CurrentLateralID)).append("_Lateral"), + // true, true); ExtractSlice(this->GetDataStorage(), image, m_CurrentMortiseID, + // MortiseFileName.append(workingDir).append(imageNode->GetName().c_str()).append("_").append(QString::number(m_CurrentMortiseID)).append("_Mortise"), + // true, true); + imageNode->SetName(imageNode->GetName().append("Extracted")); +} + +void QmitkStandardPlaneTool::InitProjectionExtractionToolWidget() +{ + connect(m_Controls.pb_LoadSavedProjections, SIGNAL(clicked()), this, SLOT(OnLoadSavedProjections())); // Load Position + connect(m_Controls.pB_ExportProjections, SIGNAL(clicked()), this, SLOT(OnExportProjections())); // Export Position + + connect(m_Controls.pB_SaveMortise, SIGNAL(clicked()), this, SLOT(OnSaveMortiseSliceClicked())); + connect(m_Controls.pB_SaveAP, SIGNAL(clicked()), this, SLOT(OnSaveAPSliceClicked())); + connect(m_Controls.pB_SaveLateral, SIGNAL(clicked()), this, SLOT(OnSaveLateralSliceClicked())); + + connect(m_Controls.pB_ShowMortise, SIGNAL(clicked()), this, SLOT(OnShowMortiseSliceClicked())); + connect(m_Controls.pB_ShowAP, SIGNAL(clicked()), this, SLOT(OnShowAPSliceClicked())); + connect(m_Controls.pB_ShowLateral, SIGNAL(clicked()), this, SLOT(OnShowLateralSliceClicked())); + + m_Controls.gB_ProjectionExtractionTool->hide(); +} + +void QmitkStandardPlaneTool::InitStandardPlaneToolWidget() +{ + connect(m_Controls.pB_savePosition, SIGNAL(clicked()), this, SLOT(AddPositionNode())); + connect(m_Controls.pB_loadPositions, SIGNAL(clicked()), this, SLOT(LoadPositionFromNode())); + connect(m_Controls.pB_LoadFromFile, SIGNAL(clicked()), this, SLOT(LoadPositionFromFile())); + connect(m_Controls.pB_SaveToFile, SIGNAL(clicked()), this, SLOT(SaveIntersectionFromNodeToFile())); + + m_Controls.gB_StandardPlaneTool->hide(); +} + +void QmitkStandardPlaneTool::ToggleStandardPlaneToolView() +{ + if (!m_showSPV) + m_Controls.gB_StandardPlaneTool->show(); + else + m_Controls.gB_StandardPlaneTool->hide(); + + m_showSPV = !m_showSPV; +} + +void QmitkStandardPlaneTool::ToggleProjectionExtractionToolView() +{ + if (!m_showPEV) + m_Controls.gB_ProjectionExtractionTool->show(); + else + m_Controls.gB_ProjectionExtractionTool->hide(); + + m_showPEV = !m_showPEV; +} + +void QmitkStandardPlaneTool::ToggleSilhouetteToolView() +{ + if (!m_showSV) + m_Controls.gB_SilhouetteTool->show(); + else + m_Controls.gB_SilhouetteTool->hide(); + + m_showSV = !m_showSV; +} + +void QmitkStandardPlaneTool::ShowAxialSlice(int sliceIndex) +{ + QmitkRenderWindow *RenderWindowAxial = this->GetRenderWindowPart()->GetQmitkRenderWindow("axial"); + RenderWindowAxial->GetSliceNavigationController()->GetSlice()->SetPos(sliceIndex); + RenderWindowAxial->GetSliceNavigationController()->SendSlice(); + + mitk::RenderingManager::GetInstance()->RequestUpdateAll(); +} + +void QmitkStandardPlaneTool::OnSelectionChanged(const mitk::DataNode *node) +{ + mitk::DataNode *selectedNode = const_cast(node); + QList nodes; + nodes.push_back(selectedNode); + // TODO because a part is needed a part is given (this should be changed to something else) + + mitk::NodePredicateDataType::Pointer isImage = mitk::NodePredicateDataType::New("Image"); + mitk::DataStorage::SetOfObjects::ConstPointer images = this->GetDataStorage()->GetSubset(isImage); + if (!images->empty()) + { + // iterate over all child nodes with NodePredicateProperty + for (mitk::DataStorage::SetOfObjects::const_iterator iter = images->begin(); iter != images->end(); ++iter) + { + if ((*iter) != node) + (*iter)->SetVisibility(false); + } + } + + this->OnSelectionChanged(this->GetSite()->GetPart(), nodes); + mitk::RenderingManager::GetInstance()->RequestUpdateAll(); + // this->OnLoadSavedProjections(); +} + +void QmitkStandardPlaneTool::OnSelectionChanged(berry::IWorkbenchPart::Pointer /*source*/, + const QList &nodes) +{ + if (!nodes.isEmpty()) + { + mitk::DataNode::Pointer referenceData = nodes.front(); + if (referenceData.IsNull()) + return; + + if (m_OldNode.IsNull()) + { + m_OldNode = referenceData; + } + else if (this->GetDataStorage()->GetNamedNode(m_WorkingNode->GetName())) + { + m_OldNode = m_WorkingNode; + } + else + { + m_OldNode = referenceData; + } + m_WorkingNode = referenceData; + // set the working node to visible + m_WorkingNode->SetProperty("visible", mitk::BoolProperty::New(true)); + + // set comboBox to reference image + disconnect(m_Controls.patientImageSelector, + SIGNAL(OnSelectionChanged(const mitk::DataNode *)), + this, + SLOT(OnComboBoxSelectionChanged(const mitk::DataNode *))); + + m_Controls.patientImageSelector->setCurrentIndex(m_Controls.patientImageSelector->Find(referenceData)); + + connect(m_Controls.patientImageSelector, + SIGNAL(OnSelectionChanged(const mitk::DataNode *)), + this, + SLOT(OnComboBoxSelectionChanged(const mitk::DataNode *))); + + mitk::NodePredicateProperty::Pointer isMarker = + mitk::NodePredicateProperty::New("isPositionNode", mitk::BoolProperty::New(true)); + // load all nodes that are children of the working node in data storage + mitk::DataStorage::SetOfObjects::ConstPointer markers = + this->GetDataStorage()->GetDerivations(m_WorkingNode, isMarker); + if (!markers->empty()) + { + // iterate over all child nodes with NodePredicateProperty + for (mitk::DataStorage::SetOfObjects::const_iterator iter = markers->begin(); iter != markers->end(); ++iter) + { + int markerId; + (*iter)->GetIntProperty("Plane Position ID", markerId); + this->LoadPositionNode((*iter), markerId); + } + } + + ////check if time node exists and if not create node (start timmer) until new node is selected + ////checkbox that is activated by the user when everything is loaded and time recording starts + // if (m_Controls.checkBoxMeasureTime->isChecked()) + //{ + // int time = 0; + // if (m_timer.isValid() && (!m_OldNode->GetIntProperty("time elapsed", time) || time == 0)) + // m_OldNode->SetIntProperty("time elapsed", m_timer.elapsed()); m_timer.restart(); + //} + } +} + +void QmitkStandardPlaneTool::OnComboBoxSelectionChanged(const mitk::DataNode *node) +{ + // if at least one node remains in the data manager + if (node != NULL) + m_Controls.patientImageSelector->setEnabled(true); + else + m_Controls.patientImageSelector->setEnabled(false); + + this->OnSelectionChanged(node); +} + +void QmitkStandardPlaneTool::AddPositionNode() +{ + MITK_INFO << m_Controls.patientImageSelector->currentIndex(); + if (m_Controls.patientImageSelector->currentIndex() < 0) + return; + // we need to save all planes... + for (int i = 1; i <= 3; i++) + { + std::stringstream multiWidgetStream; + multiWidgetStream << "stdmulti.widget"; + multiWidgetStream << i; + std::string multiWidgetString = multiWidgetStream.str(); + + mitk::BaseRenderer *m_Rendereri = mitk::BaseRenderer::GetByName(multiWidgetString.c_str()); + + if (m_Rendereri) + { + const mitk::PlaneGeometry *plane = dynamic_cast( + m_Rendereri->GetSliceNavigationController()->GetCurrentGeometry3D()) + ->GetPlaneGeometry(0); + + // Getting Service + ctkPluginContext *context = mitk::org_mitk_gui_qt_standardplanetool_Activator::GetContext(); + ctkServiceReference ppmRef = context->getServiceReference(); + mitk::PlanePositionManagerService *service = context->getService(ppmRef); + // unsigned int size = service->GetNumberOfPlanePositions(); + // context->ungetService(ppmRef); + + // node predicate property to identify all childnodes for planepositioning + mitk::NodePredicateProperty::Pointer isMarker = + mitk::NodePredicateProperty::New("isPositionNode", mitk::BoolProperty::New(true)); + // removes all nodes that are children of the working node in data storage + mitk::DataStorage::SetOfObjects::ConstPointer markers = + this->GetDataStorage()->GetDerivations(m_WorkingNode, isMarker); + if (!markers->empty() && i == 1) + { + this->GetDataStorage()->Remove(markers); + // remove all entries of these nodes in the service + for (mitk::DataStorage::SetOfObjects::const_iterator iter = markers->begin(); iter != markers->end(); ++iter) + { + int markerId; + (*iter)->GetIntProperty("Plane Position ID", markerId); + // get all position nodes in the data manager + mitk::NodePredicateProperty::Pointer isMarker = + mitk::NodePredicateProperty::New("isPositionNode", mitk::BoolProperty::New(true)); + mitk::DataStorage::SetOfObjects::ConstPointer positionNodes = this->GetDataStorage()->GetSubset(isMarker); + bool duplicate = false; + int duplicateId; + // iterate over the found nodes + for (unsigned int i = 0; i < positionNodes->size(); i++) // TODO make this efficient! + { + positionNodes->at(i)->GetIntProperty("Plane Position ID", duplicateId); + if (markerId == duplicateId) + duplicate = true; + } + if (!duplicate) + service->RemovePlanePosition(markerId); + } + } + + unsigned int id = + service->AddNewPlanePosition(plane, m_Rendereri->GetSliceNavigationController()->GetSlice()->GetPos()); + + // construct node name with some property-coding + std::stringstream nodeNameStream; + nodeNameStream << "Planeposition "; + // insert name extension axial/sagittal/coronal + switch (i) + { + case 1: + nodeNameStream << "axial"; + break; + case 2: + nodeNameStream << "sagittal"; + break; + case 3: + nodeNameStream << "coronal"; + } + + std::string nameString = nodeNameStream.str(); + + mitk::PlanarCircle::Pointer positionMarker = mitk::PlanarCircle::New(); + mitk::Point2D p1; + plane->Map(plane->GetCenter(), p1); + mitk::Point2D p2 = p1; + p2[0] -= plane->GetSpacing()[0]; + p2[1] -= plane->GetSpacing()[1]; + positionMarker->PlaceFigure(p1); + positionMarker->SetCurrentControlPoint(p1); + positionMarker->SetPlaneGeometry(const_cast(plane)); + + // the current selected node in the image selector + mitk::DataNode *imageNode(m_Controls.patientImageSelector->GetSelectedNode()); + + mitk::DataNode::Pointer planePositionNode = mitk::DataNode::New(); + + planePositionNode->SetData(positionMarker); + // TODO workaround for loading the plane positions (filling the service) + planePositionNode->SetIntProperty("sliceIndex", + m_Rendereri->GetSliceNavigationController()->GetSlice()->GetPos()); + // Saving the ID of the plane position as a property of the node + planePositionNode->SetIntProperty("Plane Position ID", id); + planePositionNode->SetProperty("name", mitk::StringProperty::New(nameString.c_str())); + planePositionNode->SetProperty("isPositionNode", mitk::BoolProperty::New(true)); + planePositionNode->SetProperty("includeInBoundingBox", mitk::BoolProperty::New(false)); + planePositionNode->SetBoolProperty("PlanarFigureInitializedWindow", true, m_Rendereri); + // TODO if the position nodes are set to helper object the loading does not work + planePositionNode->SetProperty("helper object", mitk::BoolProperty::New(false)); + // unchecks the planeposition nodes in datamanager + planePositionNode->SetProperty("visible", mitk::BoolProperty::New(false)); + + if (plane) + { + this->GetDataStorage()->Add(planePositionNode, imageNode); + } + } + } +} + +void QmitkStandardPlaneTool::LoadPositionNode(const mitk::DataNode *node, int id) +{ + MITK_INFO << m_Controls.patientImageSelector->currentIndex(); + if (m_Controls.patientImageSelector->currentIndex() < 0) + return; + + QmitkRenderWindow *selectedRenderWindow = 0; + QmitkRenderWindow *RenderWindow1 = this->GetRenderWindowPart()->GetQmitkRenderWindow("axial"); + QmitkRenderWindow *RenderWindow2 = this->GetRenderWindowPart()->GetQmitkRenderWindow("sagittal"); + QmitkRenderWindow *RenderWindow3 = this->GetRenderWindowPart()->GetQmitkRenderWindow("coronal"); + + bool PlanarFigureInitializedWindow = false; + + // find initialized renderwindow + if (node->GetBoolProperty( + "PlanarFigureInitializedWindow", PlanarFigureInitializedWindow, RenderWindow1->GetRenderer())) + { + selectedRenderWindow = RenderWindow1; + } + if (!selectedRenderWindow && + node->GetBoolProperty( + "PlanarFigureInitializedWindow", PlanarFigureInitializedWindow, RenderWindow2->GetRenderer())) + { + selectedRenderWindow = RenderWindow2; + } + if (!selectedRenderWindow && + node->GetBoolProperty( + "PlanarFigureInitializedWindow", PlanarFigureInitializedWindow, RenderWindow3->GetRenderer())) + { + selectedRenderWindow = RenderWindow3; + } + + if (selectedRenderWindow) + { + { + ctkPluginContext *context = mitk::org_mitk_gui_qt_standardplanetool_Activator::GetContext(); + ctkServiceReference ppmRef = context->getServiceReference(); + mitk::PlanePositionManagerService *service = context->getService(ppmRef); + selectedRenderWindow->GetSliceNavigationController()->ExecuteOperation(service->GetPlanePosition(id)); + context->ungetService(ppmRef); + } + mitk::RenderingManager::GetInstance()->RequestUpdateAll(); + } +} + +void QmitkStandardPlaneTool::ResetTimeProperty() +{ + // reset the saved time for plane adjustment of the working node + // m_WorkingNode->SetIntProperty("time elapsed", 0); +} + +void QmitkStandardPlaneTool::LoadPositionsFromNodes() +{ + MITK_INFO << m_Controls.patientImageSelector->currentIndex(); + if (m_Controls.patientImageSelector->currentIndex() < 0) + return; + + // get the planePositionManager Service + ctkPluginContext *context = mitk::org_mitk_gui_qt_standardplanetool_Activator::GetContext(); + ctkServiceReference ppmRef = context->getServiceReference(); + mitk::PlanePositionManagerService *service = context->getService(ppmRef); + service->RemoveAllPlanePositions(); + + int sliceIndex = 0; + mitk::PlanarCircle::Pointer positionMarker = mitk::PlanarCircle::New(); + int id; + // get position nodes in the data manager and add them with the planepositionmanger service + mitk::DataNode::Pointer coronalPosNode = this->GetDataStorage()->GetNamedNode("Planeposition coronal"); + if (coronalPosNode != nullptr) + { + coronalPosNode->GetIntProperty("sliceIndex", sliceIndex); + positionMarker = dynamic_cast(coronalPosNode->GetData()); + const mitk::Geometry2D::ConstPointer corPlane = positionMarker->GetPlaneGeometry(); + id = service->AddNewPlanePosition(corPlane, sliceIndex); + coronalPosNode->SetIntProperty("Plane Position ID", id); + MITK_INFO << "/////////////////Marker Pointer\\\\\\\\\\\\\\\\ "; + MITK_INFO << positionMarker; + } + mitk::DataNode::Pointer axialPosNode = this->GetDataStorage()->GetNamedNode("Planeposition axial"); + if (axialPosNode != nullptr) + { + axialPosNode->GetIntProperty("sliceIndex", sliceIndex); + positionMarker = dynamic_cast(axialPosNode->GetData()); + const mitk::Geometry2D::ConstPointer axiPlane = positionMarker->GetPlaneGeometry(); + id = service->AddNewPlanePosition(axiPlane, sliceIndex); + axialPosNode->SetIntProperty("Plane Position ID", id); + } + mitk::DataNode::Pointer sagittalPosNode = this->GetDataStorage()->GetNamedNode("Planeposition sagittal"); + if (sagittalPosNode != nullptr) + { + sagittalPosNode->GetIntProperty("sliceIndex", sliceIndex); + positionMarker = dynamic_cast(sagittalPosNode->GetData()); + const mitk::Geometry2D::ConstPointer sagPlane = positionMarker->GetPlaneGeometry(); + id = service->AddNewPlanePosition(sagPlane, sliceIndex); + sagittalPosNode->SetIntProperty("Plane Position ID", id); + } + + context->ungetService(ppmRef); +} + +void QmitkStandardPlaneTool::LoadPositionFromNode() +{ + MITK_INFO << m_Controls.patientImageSelector->currentIndex(); + if (m_Controls.patientImageSelector->currentIndex() < 0) + return; + + // get the planePositionManager Service + ctkPluginContext *context = mitk::org_mitk_gui_qt_standardplanetool_Activator::GetContext(); + ctkServiceReference ppmRef = context->getServiceReference(); + mitk::PlanePositionManagerService *service = context->getService(ppmRef); + + int sliceIndex = 0; + mitk::PlanarCircle::Pointer positionMarker = mitk::PlanarCircle::New(); + + mitk::NodePredicateProperty::Pointer isMarker = + mitk::NodePredicateProperty::New("isPositionNode", mitk::BoolProperty::New(true)); + mitk::DataStorage::SetOfObjects::ConstPointer markers; + + // load all nodes that are children of the working node in data storage + mitk::DataNode *imageNode(m_Controls.patientImageSelector->GetSelectedNode()); + markers = this->GetDataStorage()->GetDerivations(imageNode, isMarker); + MITK_INFO << "------------------------------------"; + MITK_INFO << "Process Image: " << imageNode->GetName(); + + if (!markers->empty()) + { + // iterate over all child nodes with NodePredicateProperty + for (mitk::DataStorage::SetOfObjects::const_iterator iter = markers->begin(); iter != markers->end(); ++iter) + { + MITK_INFO << "Marker" << std::endl; + int markerId; + std::string nodeName = (*iter)->GetName(); + // get position nodes in the data manager and add them with the planepositionmanger service + if (nodeName == "Planeposition coronal") + { + mitk::DataNode::Pointer coronalPosNode = (*iter); + coronalPosNode->GetIntProperty("sliceIndex", sliceIndex); + positionMarker = dynamic_cast(coronalPosNode->GetData()); + mitk::PlaneGeometry::ConstPointer corPlane = positionMarker->GetPlaneGeometry(); + markerId = service->AddNewPlanePosition(corPlane, sliceIndex); + coronalPosNode->SetIntProperty("Plane Position ID", markerId); + } + else if (nodeName == "Planeposition axial") + { + mitk::DataNode::Pointer axialPosNode = (*iter); + axialPosNode->GetIntProperty("sliceIndex", sliceIndex); + positionMarker = dynamic_cast(axialPosNode->GetData()); + mitk::PlaneGeometry::ConstPointer axiPlane = positionMarker->GetPlaneGeometry(); + markerId = service->AddNewPlanePosition(axiPlane, sliceIndex); + axialPosNode->SetIntProperty("Plane Position ID", markerId); + } + else if (nodeName == "Planeposition sagittal") + { + mitk::DataNode::Pointer sagittalPosNode = (*iter); + sagittalPosNode->GetIntProperty("sliceIndex", sliceIndex); + positionMarker = dynamic_cast(sagittalPosNode->GetData()); + mitk::PlaneGeometry::ConstPointer sagPlane = positionMarker->GetPlaneGeometry(); + markerId = service->AddNewPlanePosition(sagPlane, sliceIndex); + sagittalPosNode->SetIntProperty("Plane Position ID", markerId); + } + + (*iter)->GetIntProperty("Plane Position ID", markerId); + this->LoadPositionNode((*iter), markerId); + mitk::RenderingManager::GetInstance()->RequestUpdateAll(); + } + } +} + +void QmitkStandardPlaneTool::LoadPositionFromJSON() +{ + for (int i = 188; i < 221; i++) + { + if (i == 62) + continue; + MITK_INFO << "file:" << i; + mitk::Point3D intersection; + std::vector plane_normals; + std::vector imageRot(3); + auto axialPlane = mitk::PlaneGeometry::New(); + auto sagittalPlane = mitk::PlaneGeometry::New(); + auto coronalPlane = mitk::PlaneGeometry::New(); + + std::string prestring = "C:/Users/thomass/Desktop/Data/Ankle_annotations/"; + + std::string poststring = ".json"; + + int n_zero = 3; + std::string midstring = std::string(n_zero - std::to_string(i).length(), '0') + std::to_string(i); + std::string combinedString = prestring + midstring + poststring; + std::string fileName; + ReadPlaneOrientationsFromFile(combinedString, fileName, axialPlane, sagittalPlane, coronalPlane); + MITK_INFO << combinedString; + MITK_INFO << "file: " << fileName; + + std::string imageFolder = "C:/Users/thomass/Desktop/Data/" + fileName + "/"; + mitk::Image::Pointer image = mitk::IOUtil::Load(imageFolder + fileName); + m_WorkingNode = mitk::DataNode::New(); + m_WorkingNode->SetData(image); + m_WorkingNode->Modified(); + this->GetDataStorage()->Add(m_WorkingNode); + OnSelectionChanged(m_WorkingNode); + mitk::RenderingManager::GetInstance()->RequestUpdateAll(); + + std::vector RenderWindows; + RenderWindows.push_back(this->GetRenderWindowPart()->GetQmitkRenderWindow("axial")); + RenderWindows.push_back(this->GetRenderWindowPart()->GetQmitkRenderWindow("sagittal")); + RenderWindows.push_back(this->GetRenderWindowPart()->GetQmitkRenderWindow("coronal")); + + // for (unsigned int i = 0; i < RenderWindows.size(); i++) + mitk::Point3D testPointAx, testPointSag, testPointCor; + mitk::Vector3D testVector2Ax, testVector2Sag, testVector2Cor; + mitk::Vector3D testVector1Ax, testVector1Sag, testVector1Cor; + + for (int i = 0; i < 3; i++) + { + testPointAx[i] = axialPlane->GetVtkMatrix()->GetElement(i, 3); + testVector1Ax[i] = axialPlane->GetVtkMatrix()->GetElement(i, 0); + testVector2Ax[i] = axialPlane->GetVtkMatrix()->GetElement(i, 1); + testPointSag[i] = sagittalPlane->GetVtkMatrix()->GetElement(i, 3); + testVector1Sag[i] = sagittalPlane->GetVtkMatrix()->GetElement(i, 0); + testVector2Sag[i] = sagittalPlane->GetVtkMatrix()->GetElement(i, 1); + testPointCor[i] = coronalPlane->GetVtkMatrix()->GetElement(i, 3); + testVector1Cor[i] = coronalPlane->GetVtkMatrix()->GetElement(i, 0); + testVector2Cor[i] = coronalPlane->GetVtkMatrix()->GetElement(i, 1); + } + MITK_INFO << "Render planes"; + + RenderWindows[0]->GetSliceNavigationController()->ReorientSlices(testPointAx, testVector1Ax, testVector2Ax); + RenderWindows[1]->GetSliceNavigationController()->ReorientSlices(testPointSag, testVector1Sag, testVector2Sag); + RenderWindows[2]->GetSliceNavigationController()->ReorientSlices(testPointCor, testVector1Cor, testVector2Cor); + + mitk::RenderingManager::GetInstance()->RequestUpdateAll(); + // return; + + mitk::PlaneGeometry::ConstPointer sagPlane = + RenderWindows[0]->GetSliceNavigationController()->GetCurrentPlaneGeometry(); + mitk::PlaneGeometry::ConstPointer corPlane = + RenderWindows[1]->GetSliceNavigationController()->GetCurrentPlaneGeometry(); + mitk::PlaneGeometry::ConstPointer axiPlane = + RenderWindows[2]->GetSliceNavigationController()->GetCurrentPlaneGeometry(); + + /* Build coordinate system from all planes */ + mitk::Line3D lineSagAx, lineCorSag, lineCorAx; + mitk::Point3D intersectionPoint; + + axiPlane->IntersectionLine(sagPlane, lineSagAx); + sagPlane->IntersectionLine(corPlane, lineCorSag); + axiPlane->IntersectionLine(corPlane, lineCorAx); + axiPlane->IntersectionPoint(lineCorSag, intersectionPoint); + + auto combinedCoordinateSystem = vtkSmartPointer::New(); + for (int i = 0; i < 4; i++) + { + combinedCoordinateSystem->SetElement(i, 0, lineCorAx.GetDirection()[i]); + combinedCoordinateSystem->SetElement(i, 2, lineSagAx.GetDirection()[i]); + combinedCoordinateSystem->SetElement(i, 1, lineCorSag.GetDirection()[i]); + combinedCoordinateSystem->SetElement(i, 3, intersectionPoint[i]); + } + + auto combinedCoordinateSystem2 = vtkSmartPointer::New(); + for (int i = 0; i < 4; i++) + { + combinedCoordinateSystem2->SetElement(i, 2, lineCorAx.GetDirection()[i]); + combinedCoordinateSystem2->SetElement(i, 1, lineSagAx.GetDirection()[i]); + combinedCoordinateSystem2->SetElement(i, 0, lineCorSag.GetDirection()[i]); + combinedCoordinateSystem2->SetElement(i, 3, intersectionPoint[i]); + } + auto combinedCoordinateSystem3 = vtkSmartPointer::New(); + for (int i = 0; i < 4; i++) + { + combinedCoordinateSystem3->SetElement(i, 1, lineCorAx.GetDirection()[i]); + combinedCoordinateSystem3->SetElement(i, 2, lineSagAx.GetDirection()[i]); + combinedCoordinateSystem3->SetElement(i, 0, lineCorSag.GetDirection()[i]); + combinedCoordinateSystem3->SetElement(i, 3, intersectionPoint[i]); + } + + combinedCoordinateSystem->Modified(); + // MITK_INFO << "combined " << *combinedCoordinateSystem; + + auto vtkMatrix_Img = vtkSmartPointer::New(); + auto vtkMatrix_T1 = vtkSmartPointer::New(); + auto vtkMatrix_meta = vtkSmartPointer::New(); + auto vtkMatrix_new = vtkSmartPointer::New(); + auto vtkMatrix_axial = vtkSmartPointer::New(); + + // auto image = dynamic_cast(m_WorkingNode->GetData()); + // MITK_INFO << *axialPlane->GetVtkTransform()->GetMatrix(); + + // Create a sphere + auto sphereSource = vtkSmartPointer::New(); + sphereSource->SetCenter(intersectionPoint[0], intersectionPoint[1], intersectionPoint[2]); + sphereSource->SetRadius(5.0); + sphereSource->Update(); + + // Create a cube. + auto cubeSourceCor = vtkSmartPointer::New(); + cubeSourceCor->SetXLength(3.0); + cubeSourceCor->SetYLength(3.0); + cubeSourceCor->SetZLength(3.0); + cubeSourceCor->Update(); + + auto cylinderSource = vtkSmartPointer::New(); + cylinderSource->SetCenter(0.0, 0.0, 0.0); + cylinderSource->SetRadius(30.0); + cylinderSource->SetHeight(3.0); + cylinderSource->SetResolution(100); + cylinderSource->Update(); + + auto transform1 = vtkSmartPointer::New(); + transform1->SetMatrix(combinedCoordinateSystem); + transform1->Update(); + + auto transform2 = vtkSmartPointer::New(); + transform2->SetMatrix(combinedCoordinateSystem2); + transform2->Update(); + + auto transform3 = vtkSmartPointer::New(); + transform3->SetMatrix(combinedCoordinateSystem3); + transform3->Update(); + + auto transformFilter1 = vtkSmartPointer::New(); + transformFilter1->SetInputData(cylinderSource->GetOutput()); + transformFilter1->SetTransform(transform1); + transformFilter1->Update(); + + auto transformFilter2 = vtkSmartPointer::New(); + transformFilter2->SetInputData(cylinderSource->GetOutput()); + transformFilter2->SetTransform(transform2); + transformFilter2->Update(); + + auto transformFilter3 = vtkSmartPointer::New(); + transformFilter3->SetInputData(cylinderSource->GetOutput()); + transformFilter3->SetTransform(transform3); + transformFilter3->Update(); + + auto cubeCor = mitk::Surface::New(); + cubeCor->SetVtkPolyData(transformFilter2->GetOutput()); + cubeCor->Modified(); + + auto cubeAx = mitk::Surface::New(); + cubeAx->SetVtkPolyData(transformFilter3->GetOutput()); + cubeAx->Modified(); + + auto cubeSag = mitk::Surface::New(); + cubeSag->SetVtkPolyData(transformFilter1->GetOutput()); + cubeSag->Modified(); + + auto sphere = mitk::Surface::New(); + sphere->SetVtkPolyData(sphereSource->GetOutput()); + sphere->Modified(); + + mitk::Image::Pointer additionalInputImage = mitk::Image::New(); + additionalInputImage->Initialize(mitk::MakeScalarPixelType(), 3, image->GetDimensions()); + additionalInputImage->SetOrigin(image->GetGeometry()->GetOrigin()); + additionalInputImage->GetGeometry()->SetIndexToWorldTransform(image->GetGeometry()->GetIndexToWorldTransform()); + + auto surfaceToImageFilter1 = mitk::SurfaceToImageFilter::New(); + surfaceToImageFilter1->SetInput(cubeSag); + surfaceToImageFilter1->SetImage(additionalInputImage); + surfaceToImageFilter1->MakeOutputBinaryOn(); + surfaceToImageFilter1->Update(); + + auto surfaceToImageFilter2 = mitk::SurfaceToImageFilter::New(); + surfaceToImageFilter2->SetInput(cubeCor); + surfaceToImageFilter2->SetImage(additionalInputImage); + surfaceToImageFilter2->MakeOutputBinaryOn(); + surfaceToImageFilter2->Update(); + + auto surfaceToImageFilter3 = mitk::SurfaceToImageFilter::New(); + surfaceToImageFilter3->SetInput(cubeAx); + surfaceToImageFilter3->SetImage(additionalInputImage); + surfaceToImageFilter3->MakeOutputBinaryOn(); + surfaceToImageFilter3->Update(); + + auto surfaceToImageFilter4 = mitk::SurfaceToImageFilter::New(); + surfaceToImageFilter4->SetInput(sphere); + surfaceToImageFilter4->SetImage(additionalInputImage); + surfaceToImageFilter4->MakeOutputBinaryOn(); + surfaceToImageFilter4->Update(); + + auto segmSag = mitk::Image::New(); + segmSag = surfaceToImageFilter1->GetOutput(); + auto segmCor = mitk::Image::New(); + segmCor = surfaceToImageFilter2->GetOutput(); + auto segmAx = mitk::Image::New(); + segmAx = surfaceToImageFilter3->GetOutput(); + auto segmSphere = mitk::Image::New(); + segmSphere = surfaceToImageFilter4->GetOutput(); + MITK_INFO << "save segmentations"; + + std::string outputfolderSeg = "C:/Users/thomass/Desktop/Data/Ankle_segmentations/"; + std::string outputfolderImg = "C:/Users/thomass/Desktop/Data/Ankle_dataset/"; + + //mitk::IOUtil::Save(segmAx, outputfolderSeg + fileName + "_Ax.nrrd"); + //mitk::IOUtil::Save(segmSag, outputfolderSeg + fileName + "_Sag.nrrd"); + //mitk::IOUtil::Save(segmCor, outputfolderSeg + fileName + "_Cor.nrrd"); + mitk::IOUtil::Save(segmSphere, outputfolderSeg + fileName + "_Sphere.nrrd"); + //mitk::IOUtil::Save(image, outputfolderImg + fileName + ".nrrd"); + } +} + +void QmitkStandardPlaneTool::LoadPositionFromFile() +{ + LoadPositionFromJSON(); + return; + + mitk::Point3D intersection; + std::vector plane_normals; + std::vector imageRot(3); + auto axialPlane = mitk::PlaneGeometry::New(); + auto sagittalPlane = mitk::PlaneGeometry::New(); + auto coronalPlane = mitk::PlaneGeometry::New(); + + std::string prestring = "V:/BGLU/ProjektLisa/Data/SpinSiemens/Ankle_annotations/003.json"; + + int n_zero = 3; + // std::string midstring = std::string(n_zero - std::to_string(i).length(), '0') + std::to_string(i); + // std::string combinedString = prestring + midstring + poststring; + std::string fileName; + ReadPlaneOrientationsFromFile(prestring, fileName, axialPlane, sagittalPlane, coronalPlane); + // MITK_INFO << combinedString; + MITK_INFO << "file: " << fileName; + + std::string imageFolder = "V:/BGLU/ProjektLisa/Data/SpinSiemens/Ankle_dataset/"; + + std::vector RenderWindows; + RenderWindows.push_back(this->GetRenderWindowPart()->GetQmitkRenderWindow("axial")); + RenderWindows.push_back(this->GetRenderWindowPart()->GetQmitkRenderWindow("sagittal")); + RenderWindows.push_back(this->GetRenderWindowPart()->GetQmitkRenderWindow("coronal")); + + // for (unsigned int i = 0; i < RenderWindows.size(); i++) + mitk::Point3D testPointAx, testPointSag, testPointCor; + mitk::Vector3D testVector2Ax, testVector2Sag, testVector2Cor; + mitk::Vector3D testVector1Ax, testVector1Sag, testVector1Cor; + + for (int i = 0; i < 3; i++) + { + testPointAx[i] = axialPlane->GetVtkMatrix()->GetElement(i, 3); + testVector1Ax[i] = axialPlane->GetVtkMatrix()->GetElement(i, 0); + testVector2Ax[i] = axialPlane->GetVtkMatrix()->GetElement(i, 1); + testPointSag[i] = sagittalPlane->GetVtkMatrix()->GetElement(i, 3); + testVector1Sag[i] = sagittalPlane->GetVtkMatrix()->GetElement(i, 0); + testVector2Sag[i] = sagittalPlane->GetVtkMatrix()->GetElement(i, 1); + testPointCor[i] = coronalPlane->GetVtkMatrix()->GetElement(i, 3); + testVector1Cor[i] = coronalPlane->GetVtkMatrix()->GetElement(i, 0); + testVector2Cor[i] = coronalPlane->GetVtkMatrix()->GetElement(i, 1); + } + + RenderWindows[0]->GetSliceNavigationController()->ReorientSlices(testPointAx, testVector1Ax, testVector2Ax); + RenderWindows[1]->GetSliceNavigationController()->ReorientSlices(testPointSag, testVector1Sag, testVector2Sag); + RenderWindows[2]->GetSliceNavigationController()->ReorientSlices(testPointCor, testVector1Cor, testVector2Cor); + + mitk::RenderingManager::GetInstance()->RequestUpdateAll(); + + mitk::PlaneGeometry::ConstPointer sagPlane = + RenderWindows[0]->GetSliceNavigationController()->GetCurrentPlaneGeometry(); + mitk::PlaneGeometry::ConstPointer corPlane = + RenderWindows[1]->GetSliceNavigationController()->GetCurrentPlaneGeometry(); + mitk::PlaneGeometry::ConstPointer axiPlane = + RenderWindows[2]->GetSliceNavigationController()->GetCurrentPlaneGeometry(); + + /* Build coordinate system from all planes */ + mitk::Line3D lineSagAx, lineCorSag, lineCorAx; + mitk::Point3D intersectionPoint; + + axiPlane->IntersectionLine(sagPlane, lineSagAx); + + sagPlane->IntersectionLine(corPlane, lineCorSag); + axiPlane->IntersectionLine(corPlane, lineCorAx); + + axiPlane->IntersectionPoint(lineCorSag, intersectionPoint); + + MITK_INFO << intersectionPoint; + MITK_INFO << lineCorAx.GetDirection(); + MITK_INFO << lineCorSag.GetDirection(); + + auto combinedCoordinateSystem = vtkSmartPointer::New(); + for (int i = 0; i < 4; i++) + { + combinedCoordinateSystem->SetElement(i, 0, lineCorAx.GetDirection()[i]); + combinedCoordinateSystem->SetElement(i, 2, lineSagAx.GetDirection()[i]); + combinedCoordinateSystem->SetElement(i, 1, lineCorSag.GetDirection()[i]); + combinedCoordinateSystem->SetElement(i, 3, intersectionPoint[i]); + } + + auto combinedCoordinateSystem2 = vtkSmartPointer::New(); + for (int i = 0; i < 4; i++) + { + combinedCoordinateSystem2->SetElement(i, 2, lineCorAx.GetDirection()[i]); + combinedCoordinateSystem2->SetElement(i, 1, lineSagAx.GetDirection()[i]); + combinedCoordinateSystem2->SetElement(i, 0, lineCorSag.GetDirection()[i]); + combinedCoordinateSystem2->SetElement(i, 3, intersectionPoint[i]); + } + auto combinedCoordinateSystem3 = vtkSmartPointer::New(); + for (int i = 0; i < 4; i++) + { + combinedCoordinateSystem3->SetElement(i, 1, lineCorAx.GetDirection()[i]); + combinedCoordinateSystem3->SetElement(i, 2, lineSagAx.GetDirection()[i]); + combinedCoordinateSystem3->SetElement(i, 0, lineCorSag.GetDirection()[i]); + combinedCoordinateSystem3->SetElement(i, 3, intersectionPoint[i]); + } + + combinedCoordinateSystem->Modified(); + MITK_INFO << "combined " << *combinedCoordinateSystem; + + auto vtkMatrix_Img = vtkSmartPointer::New(); + auto vtkMatrix_T1 = vtkSmartPointer::New(); + auto vtkMatrix_meta = vtkSmartPointer::New(); + auto vtkMatrix_new = vtkSmartPointer::New(); + auto vtkMatrix_axial = vtkSmartPointer::New(); + + auto image = dynamic_cast(m_WorkingNode->GetData()); + + MITK_INFO << *axialPlane->GetVtkTransform()->GetMatrix(); + + // Create a cube. + auto cubeSourceCor = vtkSmartPointer::New(); + cubeSourceCor->SetXLength(3.0); + cubeSourceCor->SetYLength(3.0); + cubeSourceCor->SetZLength(3.0); + cubeSourceCor->Update(); + + auto cylinderSource = vtkSmartPointer::New(); + cylinderSource->SetCenter(0.0, 0.0, 0.0); + cylinderSource->SetRadius(30.0); + cylinderSource->SetHeight(3.0); + cylinderSource->SetResolution(100); + cylinderSource->Update(); + + auto transform1 = vtkSmartPointer::New(); + transform1->SetMatrix(combinedCoordinateSystem); + transform1->Update(); + + auto transform2 = vtkSmartPointer::New(); + transform2->SetMatrix(combinedCoordinateSystem2); + transform2->Update(); + + auto transform3 = vtkSmartPointer::New(); + transform3->SetMatrix(combinedCoordinateSystem3); + transform3->Update(); + + auto transformFilter1 = vtkSmartPointer::New(); + transformFilter1->SetInputData(cylinderSource->GetOutput()); + transformFilter1->SetTransform(transform1); + transformFilter1->Update(); + + auto transformFilter2 = vtkSmartPointer::New(); + transformFilter2->SetInputData(cylinderSource->GetOutput()); + transformFilter2->SetTransform(transform2); + transformFilter2->Update(); + + auto transformFilter3 = vtkSmartPointer::New(); + transformFilter3->SetInputData(cylinderSource->GetOutput()); + transformFilter3->SetTransform(transform3); + transformFilter3->Update(); + + auto cubeCor = mitk::Surface::New(); + cubeCor->SetVtkPolyData(transformFilter2->GetOutput()); + cubeCor->Modified(); + + auto cubeAx = mitk::Surface::New(); + cubeAx->SetVtkPolyData(transformFilter3->GetOutput()); + cubeAx->Modified(); + + auto cubeSag = mitk::Surface::New(); + cubeSag->SetVtkPolyData(transformFilter1->GetOutput()); + cubeSag->Modified(); + + mitk::Image::Pointer additionalInputImage = mitk::Image::New(); + additionalInputImage->Initialize(mitk::MakeScalarPixelType(), 3, image->GetDimensions()); + additionalInputImage->SetOrigin(image->GetGeometry()->GetOrigin()); + additionalInputImage->GetGeometry()->SetIndexToWorldTransform(image->GetGeometry()->GetIndexToWorldTransform()); + + auto surfaceToImageFilter1 = mitk::SurfaceToImageFilter::New(); + surfaceToImageFilter1->SetInput(cubeSag); + surfaceToImageFilter1->SetImage(additionalInputImage); + surfaceToImageFilter1->MakeOutputBinaryOn(); + surfaceToImageFilter1->Update(); + + auto surfaceToImageFilter2 = mitk::SurfaceToImageFilter::New(); + surfaceToImageFilter2->SetInput(cubeCor); + surfaceToImageFilter2->SetImage(additionalInputImage); + surfaceToImageFilter2->MakeOutputBinaryOn(); + surfaceToImageFilter2->Update(); + + auto surfaceToImageFilter3 = mitk::SurfaceToImageFilter::New(); + surfaceToImageFilter3->SetInput(cubeAx); + surfaceToImageFilter3->SetImage(additionalInputImage); + surfaceToImageFilter3->MakeOutputBinaryOn(); + surfaceToImageFilter3->Update(); + + auto segmSag = mitk::Image::New(); + segmSag = surfaceToImageFilter1->GetOutput(); + auto segmCor = mitk::Image::New(); + segmCor = surfaceToImageFilter2->GetOutput(); + auto segmAx = mitk::Image::New(); + segmAx = surfaceToImageFilter3->GetOutput(); + + //} + // auto segDataNode1 = mitk::DataNode::New(); + // segDataNode1->SetName(m_WorkingNode->GetName() + "_cube_Sag"); + // segDataNode1->SetColor(0.0, 0.0, 1.0); + // segDataNode1->SetData(segmSag); + + // auto segDataNode2 = mitk::DataNode::New(); + // segDataNode2->SetName(m_WorkingNode->GetName() + "_cube_Cor"); + // segDataNode2->SetColor(0.0, 1.0, 0.0); + // segDataNode2->SetData(segmCor); + + // auto segDataNode3 = mitk::DataNode::New(); + // segDataNode3->SetName(m_WorkingNode->GetName() + "_cube_Ax"); + // segDataNode3->SetColor(1.0, 0.0, 0.0); + // segDataNode3->SetData(segmAx); + + // this->GetDataStorage()->Add(segDataNode1); + // this->GetDataStorage()->Add(segDataNode2); + // this->GetDataStorage()->Add(segDataNode3); + + // QString fileName = + // QFileDialog::getOpenFileName(NULL, tr("Open Scene File with Positions"), 0, tr("MITK scene file (*.mitk)")); + // MITK_INFO << fileName; + // mitk::SceneIO::Pointer sceneIO = mitk::SceneIO::New(); + // mitk::DataStorage::Pointer dataStorage = sceneIO->LoadScene(fileName.toStdString(), 0, true); + + MITK_INFO << m_Controls.patientImageSelector->currentIndex(); + if (m_Controls.patientImageSelector->currentIndex() < 0) + return; + + // get the planePositionManager Service + ctkPluginContext *context = mitk::org_mitk_gui_qt_standardplanetool_Activator::GetContext(); + ctkServiceReference ppmRef = context->getServiceReference(); + mitk::PlanePositionManagerService *service = context->getService(ppmRef); + + int sliceIndex = 0; + mitk::PlanarCircle::Pointer positionMarker = mitk::PlanarCircle::New(); + + mitk::NodePredicateProperty::Pointer isMarker = + mitk::NodePredicateProperty::New("isPositionNode", mitk::BoolProperty::New(true)); + mitk::DataStorage::SetOfObjects::ConstPointer markers; + + // load all nodes that are children of the working node in data storage + mitk::DataNode *imageNode(m_Controls.patientImageSelector->GetSelectedNode()); + markers = this->GetDataStorage()->GetSubset(isMarker); + MITK_INFO << "------------------------------------"; + MITK_INFO << "Process Image: " << imageNode->GetName(); + + if (!markers->empty()) + { + // iterate over all child nodes with NodePredicateProperty + for (mitk::DataStorage::SetOfObjects::const_iterator iter = markers->begin(); iter != markers->end(); ++iter) + { + MITK_INFO << "Marker" << std::endl; + int markerId; + std::string nodeName = (*iter)->GetName(); + // get position nodes in the data manager and add them with the planepositionmanger service + if (nodeName == "Planeposition coronal") + { + mitk::DataNode::Pointer coronalPosNode = (*iter); + coronalPosNode->GetIntProperty("sliceIndex", sliceIndex); + positionMarker = dynamic_cast(coronalPosNode->GetData()); + mitk::PlaneGeometry::ConstPointer corPlane = positionMarker->GetPlaneGeometry(); + markerId = service->AddNewPlanePosition(corPlane, sliceIndex); + coronalPosNode->SetIntProperty("Plane Position ID", markerId); + this->GetDataStorage()->Add(coronalPosNode, imageNode); + } + else if (nodeName == "Planeposition axial") + { + mitk::DataNode::Pointer axialPosNode = (*iter); + axialPosNode->GetIntProperty("sliceIndex", sliceIndex); + positionMarker = dynamic_cast(axialPosNode->GetData()); + mitk::PlaneGeometry::ConstPointer axiPlane = positionMarker->GetPlaneGeometry(); + markerId = service->AddNewPlanePosition(axiPlane, sliceIndex); + axialPosNode->SetIntProperty("Plane Position ID", markerId); + this->GetDataStorage()->Add(axialPosNode, imageNode); + } + else if (nodeName == "Planeposition sagittal") + { + mitk::DataNode::Pointer sagittalPosNode = (*iter); + sagittalPosNode->GetIntProperty("sliceIndex", sliceIndex); + positionMarker = dynamic_cast(sagittalPosNode->GetData()); + mitk::PlaneGeometry::ConstPointer sagPlane = positionMarker->GetPlaneGeometry(); + markerId = service->AddNewPlanePosition(sagPlane, sliceIndex); + sagittalPosNode->SetIntProperty("Plane Position ID", markerId); + this->GetDataStorage()->Add(sagittalPosNode, imageNode); + } + + (*iter)->GetIntProperty("Plane Position ID", markerId); + this->LoadPositionNode((*iter), markerId); + mitk::RenderingManager::GetInstance()->RequestUpdateAll(); + } + } +} + +void QmitkStandardPlaneTool::LoadPlaneOrientationFromFile() +{ + // if (m_Controls.patientImageSelector->currentIndex() < 0) + // return; + + // QList selectedNodes = this->GetDataManagerSelection(); + + // MITK_INFO << selectedNodes.size(); + + // QString dir = + // QFileDialog::getExistingDirectory(m_Controls.pB_SaveToFile, tr("Open Directory"), tr(m_working_dir.c_str())); + + // for (mitk::DataNode::Pointer node : selectedNodes) + //{ + // std::string imageKey = node->GetName(); + + // // open file dialog to define folder for reading csv files + + // std::string baseDir = dir.toStdString(); + + // std::string plane_file_name = baseDir + "/" + imageKey; + + // mitk::Point3D intersection; + + // std::vector plane_normals; + + // std::vector imageRot(3); + + // bool success = mitk::ReadPlaneOrientationsFromFile(plane_file_name + ".csv", intersection, plane_normals, + // imageRot); + + // mitk::Image *image = dynamic_cast(node->GetData()); + + // mitk::Vector3D spacing = image->GetGeometry()->GetSpacing(); + + // // adjust normals to world coordinate system + + // for (unsigned int i = 0; i <= plane_normals.size(); i++) + + // { + // plane_normals[i][0] = plane_normals[i][0] / spacing[0]; + + // plane_normals[i][1] = plane_normals[i][1] / spacing[1]; + + // plane_normals[i][2] = plane_normals[i][2] / spacing[2]; + + // image->GetGeometry()->IndexToWorld(plane_normals[i], plane_normals[i]); + // } + + // MITK_INFO << "success: " << success; + + // if (success) + + // { + // // gloabal reinit geometry according to current image node + + // mitk::RenderingManager::GetInstance()->InitializeViews(node->GetData()->GetGeometry()); + + // std::vector RenderWindows; + + // RenderWindows.push_back(this->GetRenderWindowPart()->GetQmitkRenderWindow("axial")); + + // RenderWindows.push_back(this->GetRenderWindowPart()->GetQmitkRenderWindow("sagittal")); + + // RenderWindows.push_back(this->GetRenderWindowPart()->GetQmitkRenderWindow("coronal")); + + // for (unsigned int i = 0; i < RenderWindows.size(); i++) + + // { + // RenderWindows[i]->GetSliceNavigationController()->ReorientSlices(intersection, plane_normals[i]); + // } + + // mitk::RenderingManager::GetInstance()->RequestUpdateAll(); + + // this->AddPositionNode(imageRot, node); + // } + //} +} + +void QmitkStandardPlaneTool::SaveIntersectionFromNodeToFile() +{ + // get BaseRenderer for the three RenderWindows (sagittal, coronal, axial) + mitk::BaseRenderer::Pointer rendererAxi = + mitk::BaseRenderer::GetInstance(mitk::BaseRenderer::GetRenderWindowByName("stdmulti.widget1")); + mitk::BaseRenderer::Pointer rendererSag = + mitk::BaseRenderer::GetInstance(mitk::BaseRenderer::GetRenderWindowByName("stdmulti.widget2")); + mitk::BaseRenderer::Pointer rendererCor = + mitk::BaseRenderer::GetInstance(mitk::BaseRenderer::GetRenderWindowByName("stdmulti.widget3")); + mitk::PlaneGeometry::ConstPointer sagPlane = rendererSag->GetSliceNavigationController()->GetCurrentPlaneGeometry(); + mitk::PlaneGeometry::ConstPointer corPlane = rendererCor->GetSliceNavigationController()->GetCurrentPlaneGeometry(); + mitk::PlaneGeometry::ConstPointer axiPlane = rendererAxi->GetSliceNavigationController()->GetCurrentPlaneGeometry(); + + mitk::DataNode *imageNode(m_Controls.patientImageSelector->GetSelectedNode()); + std::string saveToMpsString = "C:/Users/thomass/Desktop/" + imageNode->GetName(); + + mitk::Line3D intersectionLine; + sagPlane->IntersectionLine(corPlane, intersectionLine); + + mitk::Point3D newPoint, tempPoint1, tempPoint2; + tempPoint1 = intersectionLine.GetPoint1(); + tempPoint2 = intersectionLine.GetPoint2(); + + newPoint[0] = (tempPoint1[0] - tempPoint2[0]) * 300; + newPoint[1] = (tempPoint1[1] - tempPoint2[1]) * 300; + newPoint[2] = (tempPoint1[2] - tempPoint2[2]) * 300; + + tempPoint1[0] = tempPoint1[0] + newPoint[0]; + tempPoint1[1] = tempPoint1[1] + newPoint[1]; + tempPoint1[2] = tempPoint1[2] + newPoint[2]; + + tempPoint2[0] = tempPoint2[0] - newPoint[0]; + tempPoint2[1] = tempPoint2[1] - newPoint[1]; + tempPoint2[2] = tempPoint2[2] - newPoint[2]; + + mitk::PointSet::Pointer intersectionSagCorPointSet = mitk::PointSet::New(); + intersectionSagCorPointSet->InsertPoint(0, tempPoint1); + intersectionSagCorPointSet->InsertPoint(1, tempPoint2); + + mitk::DataNode::Pointer intersectionSagCorPointNode = mitk::DataNode::New(); + intersectionSagCorPointNode->SetData(intersectionSagCorPointSet); + intersectionSagCorPointNode->SetName("intersection sag/cor"); + // activate property for rendering lines between the points + intersectionSagCorPointNode->SetProperty("show contour", mitk::BoolProperty::New(true)); + + mitk::IOUtil::Save(intersectionSagCorPointSet, (saveToMpsString + "_line1.mps")); + + sagPlane->IntersectionLine(axiPlane, intersectionLine); + + tempPoint1 = intersectionLine.GetPoint1(); + tempPoint2 = intersectionLine.GetPoint2(); + + newPoint[0] = (tempPoint1[0] - tempPoint2[0]) * 300; + newPoint[1] = (tempPoint1[1] - tempPoint2[1]) * 300; + newPoint[2] = (tempPoint1[2] - tempPoint2[2]) * 300; + + tempPoint1[0] = tempPoint1[0] + newPoint[0]; + tempPoint1[1] = tempPoint1[1] + newPoint[1]; + tempPoint1[2] = tempPoint1[2] + newPoint[2]; + + tempPoint2[0] = tempPoint2[0] - newPoint[0]; + tempPoint2[1] = tempPoint2[1] - newPoint[1]; + tempPoint2[2] = tempPoint2[2] - newPoint[2]; + + mitk::PointSet::Pointer intersectionSagAxiPointSet = mitk::PointSet::New(); + intersectionSagAxiPointSet->InsertPoint(0, tempPoint1); + intersectionSagAxiPointSet->InsertPoint(1, tempPoint2); + + mitk::DataNode::Pointer intersectionSagAxiPointNode = mitk::DataNode::New(); + intersectionSagAxiPointNode->SetData(intersectionSagAxiPointSet); + intersectionSagAxiPointNode->SetName("intersection sag/axi"); + // activate property for rendering lines between the points + intersectionSagAxiPointNode->SetProperty("show contour", mitk::BoolProperty::New(true)); + + mitk::IOUtil::Save(intersectionSagAxiPointSet, (saveToMpsString + "_line2.mps")); + + corPlane->IntersectionLine(axiPlane, intersectionLine); + + tempPoint1 = intersectionLine.GetPoint1(); + tempPoint2 = intersectionLine.GetPoint2(); + // tempPoint1 = (tempPoint1-tempPoint2)*10; + newPoint[0] = (tempPoint1[0] - tempPoint2[0]) * 300; + newPoint[1] = (tempPoint1[1] - tempPoint2[1]) * 300; + newPoint[2] = (tempPoint1[2] - tempPoint2[2]) * 300; + + tempPoint1[0] = tempPoint1[0] + newPoint[0]; + tempPoint1[1] = tempPoint1[1] + newPoint[1]; + tempPoint1[2] = tempPoint1[2] + newPoint[2]; + + tempPoint2[0] = tempPoint2[0] - newPoint[0]; + tempPoint2[1] = tempPoint2[1] - newPoint[1]; + tempPoint2[2] = tempPoint2[2] - newPoint[2]; + + mitk::PointSet::Pointer intersectionCorAxiPointSet = mitk::PointSet::New(); + intersectionCorAxiPointSet->InsertPoint(0, tempPoint1); + intersectionCorAxiPointSet->InsertPoint(1, tempPoint2); + + mitk::DataNode::Pointer intersectionCorAxiPointNode = mitk::DataNode::New(); + intersectionCorAxiPointNode->SetData(intersectionCorAxiPointSet); + intersectionCorAxiPointNode->SetName("intersection cor/axi"); + // activate property for rendering lines between the points + intersectionCorAxiPointNode->SetProperty("show contour", mitk::BoolProperty::New(true)); + + mitk::IOUtil::Save(intersectionCorAxiPointSet, (saveToMpsString + "_line3.mps")); +} + +void QmitkStandardPlaneTool::ShowSilhouette(int projNum, + vtkSmartPointer projMatrix, + vtkSmartPointer surface, + QString fileName, + int showSilhouette, + double offset, + int sil, + int shape) +{ + MITK_INFO << surface->GetNumberOfPoints(); + auto silhouette = vtkSmartPointer::New(); + auto camera = vtkSmartPointer::New(); + + auto camTransform = vtkSmartPointer::New(); + camTransform->SetMatrix(projMatrix); + camTransform->Update(); + + silhouette->SetProjectionMatrix(camTransform->GetMatrix()); + silhouette->SetInputData(surface); + silhouette->SetCamera(camera); + silhouette->InitializeProjectionGeometry(); + silhouette->InitializeMesh(); + silhouette->SetFeatureAngle(90); + silhouette->Update(); + + auto transform2 = vtkSmartPointer::New(); + transform2->SetMatrix(projMatrix); + transform2->Update(); + + auto transform3 = vtkSmartPointer::New(); + transform3->SetMatrix(projMatrix); + transform3->Update(); + + // auto transformFilter = vtkSmartPointer::New(); + + // if (showSilhouette ==0) + // transformFilter->SetInputData(surface); //silhouette->GetOutput()); //87 + // else + // transformFilter->SetInputData(silhouette->GetOutput()); //87 + + // transformFilter->SetTransform(transform2); + // transformFilter->Update(); + + // MITK_INFO << transformFilter->GetOutput()->GetNumberOfPoints(); + + // auto polydata = vtkSmartPointer::New(); + auto points = vtkSmartPointer::New(); + // for (int i = 0; i < transformFilter->GetOutput()->GetNumberOfPoints(); i++) + //{ + // double* p; + // transformFilter->GetOutput()->GetPoint(i, p); + // p[0] = p[0] / p[2]; + // p[1] = p[1] / p[2]; + // p[2] = p[2] / p[2]; + // points->InsertPoint(i, p); + //} + + // auto polydata = surface;//transformFilter->GetOutput();//silhouette->GetOutput(); //87 + // vtkDataArray* normals = polydata->GetPointData()->GetNormals(); + int idx = 0; + + for (int id = 0; id < silhouette->GetOutput()->GetPoints()->GetNumberOfPoints(); id++) + { + // double normal[3]; + double point[4]; + double *point4D; + // double* trafoNormal; + // double* trafoPoint; + // normals->GetTuple(id, normal); + silhouette->GetOutput()->GetPoint(id, point); + point[0] = point[0]; + point[1] = point[1]; + point[2] = point[2]; + point[3] = 1.0; + + point4D = projMatrix->MultiplyDoublePoint(point); + + point[0] = (point4D[0] / point4D[2]) + offset; + point[1] = point4D[1] / point4D[2]; + point[2] = 0; // point4D[2] / point4D[2]; + // MITK_INFO << "p" << " " << point[0] << " " << point[1] << " " << point[2]; + points->InsertPoint(id, point); + + // point[0] = (point4D[0] / point4D[2]) + offset; + // point[1] = point4D[1] / point4D[2]; + // point[2] = 1; + + // points->InsertPoint(idx++, point); + // point[0] = (point4D[0] / point4D[2]) + offset; + // point[1] = point4D[1] / point4D[2]; + // point[2] = -1; + // points->InsertPoint(idx++, point); + + // trafoNormal = transform2->TransformDoublePoint(point[0] + normal[0], point[1] + normal[1], point[2] + + // normal[2]); trafoPoint = transform3->TransformDoublePoint(polydata->GetPoint(id)); + // normals->SetTuple3(id, trafoPoint[0] - trafoNormal[0], trafoPoint[1] - trafoNormal[1], trafoPoint[2] - + // trafoNormal[2]); + } + silhouette->GetOutput()->SetPoints(points); + silhouette->GetOutput()->Modified(); + auto cleanPolyDataFilter = vtkSmartPointer::New(); + cleanPolyDataFilter->SetInputData(silhouette->GetOutput()); + cleanPolyDataFilter->Update(); + + MITK_INFO << "Silhouette: " << projNum; + typedef itk::Image ImageType; + ImageType::Pointer imageT = ImageType::New(); + ImageType::IndexType start; + start.Fill(0); + ImageType::SizeType size; + size.Fill(1024); + + ImageType::RegionType region(start, size); + imageT->SetRegions(region); + imageT->Allocate(); + imageT->FillBuffer(0); + + mitk::Image::Pointer mitkimageT; + mitk::CastToMitkImage(imageT, mitkimageT); + + // prepare the binary image's voxel grid + vtkSmartPointer whiteImage = vtkSmartPointer::New(); + whiteImage->DeepCopy(mitkimageT->GetVtkImageData()); + + //// fill the image with foreground voxels: + unsigned char outval = 255; + + // sweep polygonal data (this is the important thing with contours!) + auto extruder = vtkSmartPointer::New(); + extruder->SetInputData(cleanPolyDataFilter->GetOutput()); + extruder->SetScaleFactor(2); + extruder->SetExtrusionTypeToVectorExtrusion(); + extruder->SetVector(0, 0, 1); + extruder->Update(); + + vtkSmartPointer pol2stenc = vtkSmartPointer::New(); + pol2stenc->SetTolerance(1); // necessary because of the extruder + pol2stenc->SetInputData(extruder->GetOutput()); + pol2stenc->SetOutputOrigin(whiteImage->GetOrigin()); + pol2stenc->SetOutputSpacing(whiteImage->GetSpacing()); + pol2stenc->SetOutputWholeExtent(whiteImage->GetExtent()); + pol2stenc->Update(); + + // cut the corresponding white image and set the background: + vtkSmartPointer imgstenc = vtkSmartPointer::New(); + + imgstenc->SetInputData(whiteImage); + imgstenc->SetStencilConnection(pol2stenc->GetOutputPort()); + imgstenc->ReverseStencilOn(); + imgstenc->SetBackgroundValue(outval); + imgstenc->Update(); + + mitkimageT->SetVolume(imgstenc->GetOutput()->GetScalarPointer()); // resultImage->GetScalarPointer() + + mitk::CastToItkImage(mitkimageT, imageT); + + typedef itk::BinaryNotImageFilter BinaryNotImageFilterType; + BinaryNotImageFilterType::Pointer binaryNotFilter = BinaryNotImageFilterType::New(); + binaryNotFilter->SetInput(imageT); + binaryNotFilter->Update(); + + typedef itk::ConnectedComponentImageFilter ConnectedComponentImageFilterType; + ConnectedComponentImageFilterType::Pointer connectedFilter = ConnectedComponentImageFilterType::New(); + connectedFilter->SetInput(binaryNotFilter->GetOutput()); + connectedFilter->Update(); + + typedef itk::LabelShapeKeepNObjectsImageFilter LabelShapeKeepNObjectsImageFilterType; + LabelShapeKeepNObjectsImageFilterType::Pointer labelshapeFilter = LabelShapeKeepNObjectsImageFilterType::New(); + labelshapeFilter->SetInput(connectedFilter->GetOutput()); + labelshapeFilter->SetBackgroundValue(0); + labelshapeFilter->SetNumberOfObjects(1); + labelshapeFilter->SetAttribute(LabelShapeKeepNObjectsImageFilterType::LabelObjectType::NUMBER_OF_PIXELS); + labelshapeFilter->Update(); + + typedef itk::RescaleIntensityImageFilter RescaleType; + RescaleType::Pointer rescaleFilter = RescaleType::New(); + rescaleFilter->SetOutputMinimum(0); + rescaleFilter->SetOutputMaximum(itk::NumericTraits::max()); + rescaleFilter->SetInput(labelshapeFilter->GetOutput()); + rescaleFilter->Update(); + + BinaryNotImageFilterType::Pointer binaryNotFilter2 = BinaryNotImageFilterType::New(); + binaryNotFilter2->SetInput(rescaleFilter->GetOutput()); + binaryNotFilter2->Update(); + mitk::CastToMitkImage(binaryNotFilter2->GetOutput(), mitkimageT); // m_Segmentation + + MITK_INFO << "closing done"; + mitkimageT->Modified(); + + // transformFilter->GetOutput()->GetPointData()->SetNormals(normals); + auto projectedTibia = mitk::Surface::New(); + projectedTibia->SetVtkPolyData( + extruder->GetOutput()); // transformFilter->GetOutput());// transformFilter->GetOutput()); + // projectedTibia->GetVtkPolyData()->GetCellData()->SetNormals(normalGenerator->GetOutput()); + projectedTibia->Modified(); + + if (sil == 1) + { + m_Silhouette1->SetName(std::to_string(projNum)); + m_Silhouette1->SetData(projectedTibia); + m_Silhouette1->SetVisibility(true); + m_Silhouette1->Modified(); + } + if (sil == 2) + { + m_Silhouette2->SetName(std::to_string(projNum)); + m_Silhouette2->SetData(projectedTibia); + m_Silhouette2->SetVisibility(true); + m_Silhouette2->Modified(); + } + + if (shape == 2) + { + auto convert2DTo3DImageFilter = mitk::Convert2Dto3DImageFilter::New(); + convert2DTo3DImageFilter->SetInput(mitkimageT); + convert2DTo3DImageFilter->Update(); + + // MITK_INFO << "Segmentation: " << projNum; + // auto segmentation2 = mitk::DataNode::New(); + // segmentation2->SetProperty("name", mitk::StringProperty::New("ImageSource")); + // segmentation2->SetData(mitkimageT); + // segmentation2->SetColor(1.0, 1.0, 1.0); + // segmentation2->SetVisibility(true); + // segmentation2->SetName("segmentation"); + // this->GetDataStorage()->Add(segmentation2); + + if (fileName == nullptr) + return; + mitk::IOUtil::Save(convert2DTo3DImageFilter->GetOutput(), fileName.toStdString() + ".nrrd"); + MITK_INFO << fileName << " Saved"; + } + mitk::RenderingManager::GetInstance()->RequestUpdateAll(); +} + +void QmitkStandardPlaneTool::OnManualRegistrationClicked() +{ + // config file + boost::property_tree::ptree pt; + boost::property_tree::ini_parser::read_ini("C:/Users/thomass.AD/Desktop/config.ini", pt); + std::cout << "Projection1: " << pt.get("Section1.Projection1") << std::endl; + std::cout << "Projection1: " << pt.get("Section1.Projection2") << std::endl; + std::cout << "SurfaceFile: " << pt.get("Section1.SurfaceFile") << std::endl; + std::cout << "ProjectionFile: " << pt.get("Section1.ProjectionFile") << std::endl; + std::cout << "CSVFile: " << pt.get("Section1.CSVFile") << std::endl; + + int showSilhouette = pt.get("Section1.Silhouette"); + m_ProjNum1 = pt.get("Section1.Projection1"); + m_ProjNum2 = pt.get("Section1.Projection2"); + + // load projection matrix + QString csvFileName = QString::fromUtf8(pt.get("Section1.CSVFile").c_str()); + std::ifstream data(csvFileName.toStdString()); + + // load surface + QString surfaceFileName = QString::fromUtf8(pt.get("Section1.SurfaceFile").c_str()); + m_Surface = mitk::IOUtil::Load(surfaceFileName.toStdString()); + + QString imageFileName = QString::fromUtf8(pt.get("Section1.ProjectionFile").c_str()); + mitk::Image::Pointer image = mitk::IOUtil::Load(imageFileName.toStdString()); + + if (csvFileName.isEmpty()) + return; + + std::string line; + for (int i = 0; i < 100; i++) + { + m_ProjectionMat.push_back(vtkSmartPointer::New()); + m_ProjectionMat[i]->Identity(); + } + MITK_INFO << "projection matrix"; + int numProj = 0; + for (int i = 0; i < 100; i++) + { + for (int j = 0; j < 3; j++) + { + std::getline(data, line); + + std::stringstream lineStream(line); + std::string cell; + int num = 0; + + while (std::getline(lineStream, cell, ',')) + { + m_ProjectionMat[numProj]->SetElement(j, num++, std::stod(cell)); + } + num = 0; + } + numProj++; + } + + auto selector = mitk::ImageSliceSelector::New(); + selector->SetSliceNr(m_ProjNum1); + selector->SetInput(image); + selector->Update(); + mitk::Image::Pointer projection1 = selector->GetOutput(); + + mitk::Point3D origin = projection1->GetGeometry()->GetOrigin(); + origin[2] = 0.0; + projection1->GetGeometry()->SetOrigin(origin[0]); + projection1->Modified(); + + auto convert2DTo3DImageFilter = mitk::Convert2Dto3DImageFilter::New(); + convert2DTo3DImageFilter->SetInput(projection1); + convert2DTo3DImageFilter->Update(); + + m_Projection1 = mitk::DataNode::New(); + m_Projection1->SetProperty("name", mitk::StringProperty::New("ImageSource")); + m_Projection1->SetData(projection1); + m_Projection1->SetColor(1.0, 1.0, 1.0); + m_Projection1->SetVisibility(true); + m_Projection1->SetName("Projection1"); + this->GetDataStorage()->Add(m_Projection1); + + auto selector2 = mitk::ImageSliceSelector::New(); + selector2->SetSliceNr(m_ProjNum2); + selector2->SetInput(image); + selector2->Update(); + mitk::Image::Pointer projection2 = selector2->GetOutput(); + + mitk::Point3D origin2 = projection2->GetGeometry()->GetOrigin(); + origin2[2] = 0.0; + origin2[0] = origin2[0] + 1024; + projection2->GetGeometry()->SetOrigin(origin2); + projection2->Modified(); + + auto convert2DTo3DImageFilter2 = mitk::Convert2Dto3DImageFilter::New(); + convert2DTo3DImageFilter2->SetInput(projection2); + convert2DTo3DImageFilter2->Update(); + + m_Projection2 = mitk::DataNode::New(); + m_Projection2->SetProperty("name", mitk::StringProperty::New("ImageSource")); + m_Projection2->SetData(projection2); + m_Projection2->SetColor(1.0, 1.0, 1.0); + m_Projection2->SetVisibility(true); + m_Projection2->SetName("Projection2"); + this->GetDataStorage()->Add(m_Projection2); + + m_Silhouette1 = mitk::DataNode::New(); + m_Silhouette2 = mitk::DataNode::New(); + ShowSilhouette(m_ProjNum1, m_ProjectionMat[m_ProjNum1], m_Surface->GetVtkPolyData(), nullptr, true, 0.0, 1, 2); + ShowSilhouette(m_ProjNum2, m_ProjectionMat[m_ProjNum2], m_Surface->GetVtkPolyData(), nullptr, true, 1024.0, 2, 2); + + this->GetDataStorage()->Add(m_Silhouette1); + this->GetDataStorage()->Add(m_Silhouette2); + + mitk::RenderingManager::GetInstance()->RequestUpdateAll(); +} + +cv::Mat readEdgeMap(const QString csvFileName, mitk::DataStorage *datastorage) +{ + int mySizes[2] = {1024, 1024}; // row, column different from others + cv::Mat edgeMap = cv::Mat(mySizes[0], mySizes[1], CV_64F); + + std::ifstream file(csvFileName.toStdString()); + + for (int row = 0; row < 1024; row++) + { + std::string line; + std::getline(file, line); + if (!file.good()) + break; + + std::stringstream iss(line); + + for (int col = 0; col < 1024; col++) + { + std::string val; + std::getline(iss, val, ','); + if (!iss.good()) + break; + + std::stringstream convertor(val); + float value; + convertor >> value; + edgeMap.at(col, row) = value; + } + } + + typedef itk::Image ImageType; + ImageType::Pointer imageT = ImageType::New(); + ImageType::IndexType start; + start.Fill(0); + ImageType::SizeType size; + size.Fill(1024); + + ImageType::RegionType region(start, size); + imageT->SetRegions(region); + imageT->Allocate(); + imageT->FillBuffer(0); + + for (int n = 0; n < 1024; n++) + { + for (int m = 0; m < 1024; m++) + { + if (edgeMap.at(m, n) > 0) + { + ImageType::IndexType pIndex; + pIndex[0] = m; + pIndex[1] = n; + ImageType::PixelType value = 1000.0; + + imageT->SetPixel(pIndex, value); + ImageType::PixelType pixelValue = imageT->GetPixel(pIndex); + } + } + } + mitk::Image::Pointer mitkimageT; + mitk::CastToMitkImage(imageT, mitkimageT); + + auto datanode2 = mitk::DataNode::New(); + datanode2->SetName("ContourImage"); + // datanode->SetBoolProperty("show contour", true); + // datanode2->SetProperty("contourcolor", ColorProperty::New(0.0, 1.0, 0.0)); + datanode2->SetData(mitkimageT); + datastorage->Add(datanode2); + + return edgeMap; +} + +void QmitkStandardPlaneTool::OnSliderValueChanged() +{ + MITK_INFO << "Change Value"; + typedef itk::Euler3DTransform TransformType; + TransformType::Pointer itkTransform = TransformType::New(); + TransformType::ParametersType rigid_params(6); + rigid_params[0] = m_Controls.sB_Tx->value(); + rigid_params[1] = m_Controls.sB_Ty->value(); + rigid_params[2] = m_Controls.sB_Tz->value(); + rigid_params[3] = m_Controls.sB_alpha->value(); + rigid_params[4] = m_Controls.sB_beta->value(); + rigid_params[5] = m_Controls.sB_gamma->value(); + + itkTransform->SetParameters(rigid_params); + + auto vtkMatrix = vtkSmartPointer::New(); + mitk::TransferItkTransformToVtkMatrix(itkTransform.GetPointer(), vtkMatrix); + auto transformFilter = vtkSmartPointer::New(); + auto vtktransform = vtkSmartPointer::New(); + vtktransform->SetMatrix(vtkMatrix); + transformFilter->SetInputData(m_Surface->GetVtkPolyData()); + transformFilter->SetTransform(vtktransform); + transformFilter->Update(); + + this->ShowSilhouette(m_ProjNum1, m_ProjectionMat[m_ProjNum1], transformFilter->GetOutput(), nullptr, 1, 0.0, 1, 1); + this->ShowSilhouette(m_ProjNum2, m_ProjectionMat[m_ProjNum2], transformFilter->GetOutput(), nullptr, 1, 1024.0, 2, 1); +} + +void ConvertToSurface(mitk::Image::Pointer segmentation, mitk::Surface::Pointer &surface) +{ + // Create a contour from the binary image + auto surfaceFilter = mitk::ManualSegmentationToSurfaceFilter::New(); + surfaceFilter->SetInput(segmentation); + surfaceFilter->SetThreshold(0.5); // expects binary image with zeros and ones + surfaceFilter->SetGaussianStandardDeviation(1.5); + surfaceFilter->SetUseGaussianImageSmooth(true); // apply gaussian to thresholded image ? + surfaceFilter->SetSmooth(true); + surfaceFilter->InterpolationOn(); + surfaceFilter->SetMedianFilter3D(true); // apply median to segmentation before marching cubes ? + surfaceFilter->SetDecimate(mitk::ImageToSurfaceFilter::NoDecimation); + surfaceFilter->UpdateLargestPossibleRegion(); + + ////----------------- REMESHING ---------------------------------- + // auto remesher = mitk::ACVD::RemeshFilter::New(); + // remesher->SetInput(surfaceFilter->GetOutput()); + // remesher->SetTimeStep(0); + // remesher->SetNumVertices(1000); + // remesher->SetGradation(1.0); + //// remesher->SetBoundaryFixing(true) + // remesher->SetSubsampling(50); + // remesher->SetEdgeSplitting(0.0); + // remesher->SetOptimizationLevel(1); + // remesher->SetForceManifold(false); + // remesher->SetBoundaryFixing(false); + // remesher->Update(); + surface = surfaceFilter->GetOutput()->Clone(); +} + +void QmitkStandardPlaneTool::OnExtractContourAndSlicesClicked() +{ + // config file + boost::property_tree::ptree pt; + boost::property_tree::ini_parser::read_ini("C:/Users/thomass.AD/Desktop/configMask.ini", pt); + std::cout << "Folder name2: " << pt.get("Section1.Folder100") << std::endl; + + mitk::NodePredicateDataType::Pointer isImage = mitk::NodePredicateDataType::New("Image"); + mitk::DataStorage::SetOfObjects::ConstPointer images = this->GetDataStorage()->GetSubset(isImage); + if (!images->empty()) + { + // iterate over all child nodes with NodePredicateProperty + for (mitk::DataStorage::SetOfObjects::const_iterator iter = images->begin(); iter != images->end(); ++iter) + { + mitk::DataNode::Pointer node = (*iter); + mitk::Image::Pointer image = dynamic_cast(node->GetData()); + QString nodeName = QString::fromUtf8((*iter)->GetName().c_str()); + + QString folderName = QString::fromUtf8(pt.get("Section1.Folder").c_str()); + QString fibulaFileName = folderName + (QString("/Reorganized/")) + (nodeName + (QString("_Fibula.nrrd"))); + QString talusFileName = folderName + (QString("/Reorganized/")) + (nodeName + (QString("_Talus.nrrd"))); + QString tibiaFileName = folderName + (QString("/Reorganized/")) + (nodeName + (QString("_Tibia.nrrd"))); + QString csvFileName = + folderName + (QString("/XA/")) + (nodeName + (QString("/"))) + (nodeName + (QString(".csv"))); + QString projectionFileName = + folderName + (QString("/XA/")) + (nodeName + (QString("/"))) + (nodeName + (QString(".IMA"))); + QString outputFileName = folderName + (QString("/XA/")) + (nodeName) + (QString("/masks/")); + std::ifstream data(csvFileName.toStdString()); + + // load segmentations + + mitk::Image::Pointer fibula = mitk::IOUtil::Load(fibulaFileName.toStdString()); + mitk::Image::Pointer tibia = mitk::IOUtil::Load(tibiaFileName.toStdString()); + mitk::Image::Pointer talus = mitk::IOUtil::Load(talusFileName.toStdString()); + mitk::Image::Pointer projection = mitk::IOUtil::Load(projectionFileName.toStdString()); + + auto surfaceTibia = mitk::Surface::New(); + auto surfaceFibula = mitk::Surface::New(); + auto surfaceTalus = mitk::Surface::New(); + + ConvertToSurface(fibula, surfaceFibula); + ConvertToSurface(tibia, surfaceTibia); + ConvertToSurface(talus, surfaceTalus); + + if (csvFileName.isEmpty()) + return; + + m_ProjectionMat.empty(); + std::string line; + for (int i = 0; i < 100; i++) + { + m_ProjectionMat.push_back(vtkSmartPointer::New()); + m_ProjectionMat[i]->Identity(); + } + MITK_INFO << "projection matrix"; + int numProj = 0; + for (int i = 0; i < 100; i++) + { + for (int j = 0; j < 3; j++) + { + std::getline(data, line); + + std::stringstream lineStream(line); + std::string cell; + int num = 0; + + while (std::getline(lineStream, cell, ',')) + { + m_ProjectionMat[numProj]->SetElement(j, num++, std::stod(cell)); + } + num = 0; + } + numProj++; + } + + auto selector = mitk::ImageSliceSelector::New(); + selector->SetSliceNr(m_ProjNum1); + selector->SetInput(image); + selector->Update(); + mitk::Image::Pointer projection1 = selector->GetOutput(); + + mitk::Point3D origin = projection1->GetGeometry()->GetOrigin(); + origin[2] = 0.0; + projection1->GetGeometry()->SetOrigin(origin[0]); + projection1->Modified(); + + auto convert2DTo3DImageFilter = mitk::Convert2Dto3DImageFilter::New(); + convert2DTo3DImageFilter->SetInput(projection1); + convert2DTo3DImageFilter->Update(); + + m_Silhouette1 = mitk::DataNode::New(); + m_Silhouette2 = mitk::DataNode::New(); + MITK_INFO << *m_ProjectionMat[0]; + + for (int i = 0; i < 100; i++) + { + mitk::Image::Pointer projection1 = selector->GetOutput(); + ShowSilhouette(i, + m_ProjectionMat[i], + surfaceFibula->GetVtkPolyData(), + outputFileName + (nodeName + ((("_Fibula" + std::to_string(i)).c_str()))), + true, + 0.0, + 1, + 2); + ShowSilhouette(i, + m_ProjectionMat[i], + surfaceTibia->GetVtkPolyData(), + outputFileName + (nodeName + ((("_Tibia" + std::to_string(i)).c_str()))), + true, + 0.0, + 1, + 2); + ShowSilhouette(i, + m_ProjectionMat[i], + surfaceTalus->GetVtkPolyData(), + outputFileName + (nodeName + ((("_Talus" + std::to_string(i)).c_str()))), + true, + 0.0, + 1, + 2); + ExtractSlice(this->GetDataStorage(), + projection, + i, + outputFileName + (nodeName + ((('_' + std::to_string(i)).c_str()))), + true, + true); + } + + /* auto segmentation2 = mitk::DataNode::New(); + segmentation2->SetProperty("name", mitk::StringProperty::New("ImageSource")); + segmentation2->SetData(mask); + segmentation2->SetColor(1.0, 1.0, 1.0); + segmentation2->SetVisibility(true); + segmentation2->SetName("segmentation1"); + this->GetDataStorage()->Add(segmentation2);*/ + } + } +} + +void QmitkStandardPlaneTool::OnRegistrationClicked() +{ + // config file + boost::property_tree::ptree pt; + boost::property_tree::ini_parser::read_ini("C:/Users/thomass.AD/Desktop/config.ini", pt); + std::cout << "Projection1: " << pt.get("Section1.Projection1") << std::endl; + std::cout << "Projection1: " << pt.get("Section1.Projection2") << std::endl; + std::cout << "SurfaceFile: " << pt.get("Section1.SurfaceFile") << std::endl; + std::cout << "EdgeMap1: " << pt.get("Section1.EdgeMapFile1") << std::endl; + std::cout << "EdgeMap2: " << pt.get("Section1.EdgeMapFile2") << std::endl; + std::cout << "ProjectionFile: " << pt.get("Section1.ProjectionFile") << std::endl; + std::cout << "CSVFile: " << pt.get("Section1.CSVFile") << std::endl; + + int showSilhouette = pt.get("Section1.Silhouette"); + int projNum1 = pt.get("Section1.Projection1"); + int projNum2 = pt.get("Section1.Projection2"); + // return; + // load edgemap 1 and 2 + QString csvFileNameEdgeMap1 = + QString::fromUtf8(pt.get("Section1.EdgeMapFile1") + .c_str()); // QFileDialog::getOpenFileName(nullptr, "Select a .csv file", "", "CSV (*.csv)"); + cv::Mat edgeMap1 = readEdgeMap(csvFileNameEdgeMap1, this->GetDataStorage()); + QString csvFileNameEdgeMap2 = + QString::fromUtf8(pt.get("Section1.EdgeMapFile2") + .c_str()); // QFileDialog::getOpenFileName(nullptr, "Select a .csv file", "", "CSV (*.csv)"); + cv::Mat edgeMap2 = readEdgeMap(csvFileNameEdgeMap2, this->GetDataStorage()); + + // load projection matrix + QString csvFileName = + QString::fromUtf8(pt.get("Section1.CSVFile") + .c_str()); // QFileDialog::getOpenFileName(nullptr, "Select a .csv file", "", "CSV (*.csv)"); + std::ifstream data(csvFileName.toStdString()); + + // load surface + QString surfaceFileName = QString::fromUtf8( + pt.get("Section1.SurfaceFile") + .c_str()); // QFileDialog::getOpenFileName(nullptr, "Select a surface file", "", "Surface (*.stl)"); + m_Surface = mitk::IOUtil::Load(surfaceFileName.toStdString()); + + QString imageFileName = QString::fromUtf8(pt.get("Section1.ProjectionFile").c_str()); + mitk::Image::Pointer image = mitk::IOUtil::Load(imageFileName.toStdString()); + + if (csvFileName.isEmpty()) + return; + + std::string line; + std::vector> projectionMat; + for (int i = 0; i < 100; i++) + { + projectionMat.push_back(vtkSmartPointer::New()); + projectionMat[i]->Identity(); + } + MITK_INFO << "projection matrix"; + int numProj = 0; + for (int i = 0; i < 100; i++) + { + for (int j = 0; j < 3; j++) + { + std::getline(data, line); + + std::stringstream lineStream(line); + std::string cell; + int num = 0; + + while (std::getline(lineStream, cell, ',')) + { + projectionMat[numProj]->SetElement(j, num++, std::stod(cell)); + } + num = 0; + } + numProj++; + } + + typedef itk::Euler3DTransform TransformType; + TransformType::Pointer itkTransform = TransformType::New(); + TransformType::ParametersType rigid_params(6); + rigid_params[0] = 0; + rigid_params[1] = 0; + rigid_params[2] = 0; + rigid_params[3] = 0; + rigid_params[4] = 0; + rigid_params[5] = 0; + + itkTransform->SetParameters(rigid_params); + + auto vtkMatrix = vtkSmartPointer::New(); + mitk::TransferItkTransformToVtkMatrix(itkTransform.GetPointer(), vtkMatrix); + auto transformFilter = vtkSmartPointer::New(); + auto vtktransform = vtkSmartPointer::New(); + vtktransform->SetMatrix(vtkMatrix); + transformFilter->SetInputData(m_Surface->GetVtkPolyData()); + transformFilter->SetTransform(vtktransform); + transformFilter->Update(); + m_Surface->SetVtkPolyData(transformFilter->GetOutput()); + m_Surface->Modified(); + + // TODO: Merge to one step ....Add images to be used for gradient calculation + auto selector = mitk::ImageSliceSelector::New(); + selector->SetSliceNr(m_ProjNum1); + selector->SetInput(image); + selector->Update(); + mitk::Image::Pointer projection1 = selector->GetOutput(); + + mitk::Point3D origin = projection1->GetGeometry()->GetOrigin(); + origin[2] = 0.0; + projection1->GetGeometry()->SetOrigin(origin[0]); + projection1->Modified(); + + auto convert2DTo3DImageFilter = mitk::Convert2Dto3DImageFilter::New(); + convert2DTo3DImageFilter->SetInput(projection1); + convert2DTo3DImageFilter->Update(); + + m_Projection1 = mitk::DataNode::New(); + m_Projection1->SetProperty("name", mitk::StringProperty::New("ImageSource")); + m_Projection1->SetData(projection1); + m_Projection1->SetColor(1.0, 1.0, 1.0); + m_Projection1->SetVisibility(true); + m_Projection1->SetName("Projection1"); + this->GetDataStorage()->Add(m_Projection1); + + auto selector2 = mitk::ImageSliceSelector::New(); + selector2->SetSliceNr(m_ProjNum2); + selector2->SetInput(image); + selector2->Update(); + mitk::Image::Pointer projection2 = selector2->GetOutput(); + + mitk::Point3D origin2 = projection2->GetGeometry()->GetOrigin(); + origin2[2] = 0.0; + origin2[0] = origin2[0] + 1024; + projection2->GetGeometry()->SetOrigin(origin2); + projection2->Modified(); + + auto convert2DTo3DImageFilter2 = mitk::Convert2Dto3DImageFilter::New(); + convert2DTo3DImageFilter2->SetInput(projection2); + convert2DTo3DImageFilter2->Update(); + + m_Projection2 = mitk::DataNode::New(); + m_Projection2->SetProperty("name", mitk::StringProperty::New("ImageSource")); + m_Projection2->SetData(projection2); + m_Projection2->SetColor(1.0, 1.0, 1.0); + m_Projection2->SetVisibility(true); + m_Projection2->SetName("Projection2"); + this->GetDataStorage()->Add(m_Projection2); + + const unsigned int Dimension = 2; + typedef unsigned short PixelTypeUShort; + typedef itk::Image ImageTypeUShort; + + ImageTypeUShort::Pointer image1ITK = ImageTypeUShort::New(); + mitk::CastToItkImage(projection1, image1ITK); + cv::Mat image1 = itk::OpenCVImageBridge::ITKImageToCVMat(image1ITK); + ImageTypeUShort::Pointer image2ITK = ImageTypeUShort::New(); + mitk::CastToItkImage(projection2, image2ITK); + cv::Mat image2 = itk::OpenCVImageBridge::ITKImageToCVMat(image2ITK); + + // cv::imshow("orientation", image1); + // return; + + MITK_INFO << "set up registration"; + mitk2D3DShapeRegistration registration; // = mitk2D3DShapeRegistration::New(); + registration.SetSurface(m_Surface); + MITK_INFO << "set up edgemaps"; + registration.SetEdgeMap1(edgeMap1); + registration.SetEdgeMap2(edgeMap2); + registration.SetImage1(image1); + registration.SetImage2(image2); + int mySizes[2] = {1024, 1024}; // row, column different from others + cv::Mat gradImage = cv::Mat(mySizes[0], mySizes[1], CV_32F); + + MITK_INFO << "set up projection matrix"; + registration.SetProjectionMatrix1(projectionMat[projNum1]); + registration.SetProjectionMatrix2(projectionMat[projNum2]); + MITK_INFO << "generate"; + registration.Generate(); + return; + m_Surface->SetVtkPolyData(registration.GetResultSurface()); + m_Surface->Modified(); + m_Silhouette1 = mitk::DataNode::New(); + m_Silhouette2 = mitk::DataNode::New(); + + ShowSilhouette(projNum1, projectionMat[projNum1], m_Surface->GetVtkPolyData(), nullptr, true, 0.0, 1, 1); + ShowSilhouette(projNum2, projectionMat[projNum2], m_Surface->GetVtkPolyData(), nullptr, true, 0.0, 2, 1); + + this->GetDataStorage()->Add(m_Silhouette1); + this->GetDataStorage()->Add(m_Silhouette2); +} + +void QmitkStandardPlaneTool::OnExtractProjectionsClicked() +{ + mitk::NodePredicateDataType::Pointer isImage = mitk::NodePredicateDataType::New("Image"); + mitk::DataStorage::SetOfObjects::ConstPointer images = this->GetDataStorage()->GetSubset(isImage); + if (!images->empty()) + { + // iterate over all child nodes with NodePredicateProperty + for (mitk::DataStorage::SetOfObjects::const_iterator iter = images->begin(); iter != images->end(); ++iter) + { + mitk::DataNode::Pointer node = (*iter); + mitk::Image::Pointer image = dynamic_cast(node->GetData()); + + for (int i = 0; i < 100; i++) + { + QString nodeName = QString::fromUtf8((*iter)->GetName().c_str()); + ExtractSlice(this->GetDataStorage(), + image, + i, + QString("C:/Users/thomass.AD/Desktop/Data/SilhouetteExp/") + + (nodeName + ((('_' + std::to_string(i)).c_str()))), + true, + true); + } + } + } +} + +void QmitkStandardPlaneTool::OnExtractSilhouetteClicked() +{ + QString surfaceFileName = QFileDialog::getOpenFileName(nullptr, "Select a surface file", "", "Surface (*.stl)"); + QString imageFileName = QFileDialog::getOpenFileName(nullptr, "Select an image file", "", "Image (*.nrrd)"); + + if (surfaceFileName.isEmpty()) + return; + + if (imageFileName.isEmpty()) + return; + + mitk::Image::Pointer image = mitk::IOUtil::Load(imageFileName.toStdString()); + // ******************************************************* + // load projection matrices + // ******************************************************* + QString csvFileName = QFileDialog::getOpenFileName(nullptr, "Select a .csv file", "", "CSV (*.csv)"); + std::ifstream data(csvFileName.toStdString()); + + if (csvFileName.isEmpty()) + return; + + std::string line; + std::vector> projectionMat; + for (int i = 0; i < 100; i++) + { + projectionMat.push_back(vtkSmartPointer::New()); + projectionMat[i]->Identity(); + } + MITK_INFO << "projection matrix"; + int numProj = 0; + for (int i = 0; i < 100; i++) + { + for (int j = 0; j < 3; j++) + { + std::getline(data, line); + + std::stringstream lineStream(line); + std::string cell; + int num = 0; + + while (std::getline(lineStream, cell, ',')) + { + projectionMat[numProj]->SetElement(j, num++, std::stod(cell)); + // MITK_INFO << projectionMat[numProj]->GetElement(j, num - 1); + } + num = 0; + } + numProj++; + } + MITK_INFO << "Decompose Matrix"; + // Create a synthetic projection matrix + cv::Mat P(3, 4, cv::DataType::type); + P.at(0, 0) = projectionMat[0]->GetElement(0, 0); + P.at(1, 0) = projectionMat[0]->GetElement(1, 0); + P.at(2, 0) = projectionMat[0]->GetElement(2, 0); + + P.at(0, 1) = projectionMat[0]->GetElement(0, 1); + P.at(1, 1) = projectionMat[0]->GetElement(1, 1); + P.at(2, 1) = projectionMat[0]->GetElement(2, 1); + + P.at(0, 2) = projectionMat[0]->GetElement(0, 2); + P.at(1, 2) = projectionMat[0]->GetElement(1, 2); + P.at(2, 2) = projectionMat[0]->GetElement(2, 2); + + P.at(0, 3) = projectionMat[0]->GetElement(0, 3); + P.at(1, 3) = projectionMat[0]->GetElement(1, 3); + P.at(2, 3) = projectionMat[0]->GetElement(2, 3); + + MITK_INFO << P.at(2, 3); + + // Decompose the projection matrix into: + cv::Mat K(3, 3, cv::DataType::type); // intrinsic parameter matrix + cv::Mat R(3, 3, cv::DataType::type); // rotation matrix + cv::Mat T(4, 1, cv::DataType::type); // translation vector + + // cv::decomposeProjectionMatrix(P, K, R, T); + + // MITK_INFO << "K Matrix"; + // MITK_INFO << K.at(0, 0); + // MITK_INFO << K.at(1, 0); + // MITK_INFO << K.at(2, 0); + // MITK_INFO << K.at(0, 1); + // MITK_INFO << K.at(1, 1); + // MITK_INFO << K.at(2, 1); + // MITK_INFO << K.at(0, 2); + // MITK_INFO << K.at(1, 2); + // MITK_INFO << K.at(2, 2); + + // MITK_INFO << "R Matrix"; + // MITK_INFO << R.at(0, 0); + // MITK_INFO << R.at(1, 0); + // MITK_INFO << R.at(2, 0); + // MITK_INFO << R.at(0, 1); + // MITK_INFO << R.at(1, 1); + // MITK_INFO << R.at(2, 1); + // MITK_INFO << R.at(0, 2); + // MITK_INFO << R.at(1, 2); + // MITK_INFO << R.at(2, 2); + + // MITK_INFO << "T Matrix"; + // MITK_INFO << T.at(0, 0); + // MITK_INFO << T.at(1, 0); + // MITK_INFO << T.at(2, 0); + // MITK_INFO << T.at(3, 0); + + mitk::Surface::Pointer surface = mitk::IOUtil::Load(surfaceFileName.toStdString()); + mitk::Surface::Pointer surface2 = mitk::IOUtil::Load(surfaceFileName.toStdString()); + + for (int i = 0; i < 40; i = i + 39) + { + ShowSilhouette(i, projectionMat[i], surface2->Clone()->GetVtkPolyData(), nullptr, 1, 0.0, 1, 1); + ExtractSlice(this->GetDataStorage(), + image, + i, + QString("C:/Users/thomass.AD/Desktop/Data/SilhouetteExp/test_") + ((std::to_string(i).c_str())), + false, + true); + } + // + // std::cout << "There are " << surface->GetVtkPolyData()->GetNumberOfCells() << " cells." << std::endl; + // std::cout << "There are " << surface->GetVtkPolyData()->GetNumberOfPoints() << " points." << std::endl; + // + // + // + // auto extractEdges = vtkSmartPointer::New(); + // extractEdges->SetInputData(surface->GetVtkPolyData()); + // extractEdges->Update(); + // + // vtkCellArray* lines = extractEdges->GetOutput()->GetLines(); + // vtkPoints* points = extractEdges->GetOutput()->GetPoints(); + // + // std::cout << std::endl << "Edges" << endl << "----------" << std::endl; + // std::cout << "There are " << lines->GetNumberOfCells() << " cells." << std::endl; + // std::cout << "There are " << points->GetNumberOfPoints() << " points." << std::endl; + // + // //// Traverse all of the edges + // //for (vtkIdType i = 0; i < extractEdges->GetOutput()->GetNumberOfCells()/100; i++) + // //{ + // // std::cout << "Type: " << extractEdges->GetOutput()->GetCell(i)->GetClassName() << std::endl; + // // vtkSmartPointer line = vtkLine::SafeDownCast(extractEdges->GetOutput()->GetCell(i)); + // //// std::cout << "Line " << i << " : " << *line << std::endl; + // + // //} + // + // //************************************************************************ + // + // auto triangleFilter = vtkSmartPointer::New(); + // triangleFilter->SetInputData(surface->GetVtkPolyData()); + // triangleFilter->Update(); + // + // vtkIdType cellId = 10; + // + // auto cellPointIds = vtkSmartPointer::New(); + // triangleFilter->GetOutput()->GetCellPoints(cellId, cellPointIds); + // + // + // + // + // + // + // //for (i = 0; iGetVtkPolyData()->GetNumberOfCells(); i++) + // //{ + // // cell = input->GetCell(i); + // // for (j = 0; j<3; j++) + // // { + // // cellEdgeNeighbors->Initialize(); + // // id0 = cell->GetEdge(j)->GetPointIds()->GetId(0); + // // id1 = cell->GetEdge(j)->GetPointIds()->GetId(1); + // // input->GetCellEdgeNeighbors(i, id0, id1, cellEdgeNeighbors); + // // if (cellEdgeNeighbors->GetNumberOfIds() == 0) + // // { + // // boundaryIds->InsertUniqueId(id0); + // // boundaryIds->InsertUniqueId(id1); + // // } + // // } + // //} + // + // + // std::list neighbors; + // + // for (vtkIdType i = 0; i < cellPointIds->GetNumberOfIds(); i++) + // { + // auto idList = vtkSmartPointer::New(); + // + // //add one of the edge points + // idList->InsertNextId(cellPointIds->GetId(i)); + // + // //add the other edge point + // if (i + 1 == cellPointIds->GetNumberOfIds()) + // { + // idList->InsertNextId(cellPointIds->GetId(0)); + // } + // else + // { + // idList->InsertNextId(cellPointIds->GetId(i + 1)); + // } + // + // //get the neighbors of the cell + // auto neighborCellIds = vtkSmartPointer::New(); + // + // triangleFilter->GetOutput()->GetCellNeighbors(cellId, idList, neighborCellIds); + // // triangleFilter->GetOutput()->GetCellEdgeNeighbors(cellId, 1, 0, neighborCellIds); + // + // /* id0 = cell->GetEdge(j)->GetPointIds()->GetId(0); + // id1 = cell->GetEdge(j)->GetPointIds()->GetId(1); + // input->GetCellEdgeNeighbors(i, id0, id1, cellEdgeNeighbors); + //*/ + // + // + // for (vtkIdType j = 0; j < neighborCellIds->GetNumberOfIds(); j++) + // { + // neighbors.push_back(neighborCellIds->GetId(j)); + // } + // + // } + // + // std::cout << "Edge neighbor ids are: " << std::endl; + // + // for (std::list::iterator it1 = neighbors.begin(); it1 != neighbors.end(); it1++) + // { + // std::cout << " " << *it1; + // + // } + // std::cout << std::endl; + // + // + // + // + // // Create a dataset with the cell of interest + // { + // auto ids = vtkSmartPointer::New(); + // ids->SetNumberOfComponents(1); + // ids->InsertNextValue(cellId); + // + // auto selectionNode = vtkSmartPointer::New(); + // selectionNode->SetFieldType(vtkSelectionNode::CELL); + // selectionNode->SetContentType(vtkSelectionNode::INDICES); + // selectionNode->SetSelectionList(ids); + // + // auto selection = vtkSmartPointer::New(); + // selection->AddNode(selectionNode); + // + // auto extractSelection = vtkSmartPointer::New(); + // extractSelection->SetInputData(0, surface->GetVtkPolyData()); + // extractSelection->SetInputData(1, selection); + // extractSelection->Update(); + // + // } + // + // //mainCellActor->GetProperty()->SetColor(1, 0, 0); + // + // // Create a dataset with the neighbor cells + // { + // auto ids = vtkSmartPointer::New(); + // ids->SetNumberOfComponents(1); + // for (std::list::iterator it1 = neighbors.begin(); it1 != neighbors.end(); it1++) + // { + // ids->InsertNextValue(*it1); + // } + // + // vtkSmartPointer selectionNode = + // vtkSmartPointer::New(); + // selectionNode->SetFieldType(vtkSelectionNode::CELL); + // selectionNode->SetContentType(vtkSelectionNode::INDICES); + // selectionNode->SetSelectionList(ids); + // + // vtkSmartPointer selection = + // vtkSmartPointer::New(); + // selection->AddNode(selectionNode); + // + // vtkSmartPointer extractSelection = + // vtkSmartPointer::New(); + // extractSelection->SetInputData(surface->GetVtkPolyData()); + // extractSelection->SetInputData(1, selection); + // extractSelection->Update(); + // + // } + // + // neighborCellsActor->GetProperty()->SetColor(0, 1, 0); + + //************************************************************************* +} + +void QmitkStandardPlaneTool::SetFocus() {} diff --git a/Plugins/org.mitk.gui.qt.standardplanetool/src/internal/QmitkStandardPlaneTool.h b/Plugins/org.mitk.gui.qt.standardplanetool/src/internal/QmitkStandardPlaneTool.h new file mode 100644 index 0000000000..bafb54174c --- /dev/null +++ b/Plugins/org.mitk.gui.qt.standardplanetool/src/internal/QmitkStandardPlaneTool.h @@ -0,0 +1,183 @@ +/*========================================================================= + +Program: Medical Imaging & Interaction Toolkit +Language: C++ +Date: $Date$ +Version: $Revision$ + +Copyright (c) German Cancer Research Center, Division of Medical and +Biological Informatics. All rights reserved. +See MITKCopyright.txt or http://www.mitk.org/copyright.html for details. + +This software is distributed WITHOUT ANY WARRANTY; without even +the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR +PURPOSE. See the above copyright notices for more information. + +=========================================================================*/ + +#ifndef QmitkStandardPlaneTool_h +#define QmitkStandardPlaneTool_h + +#include + +#ifdef WIN32 +#pragma warning( disable : 4250 ) +#endif + +#include +#include "QmitkRegisterClasses.h" + +#include "itkCommand.h" +#include +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include "usServiceRegistration.h" + +#include + +// QT +#include + +#include "ui_QmitkStandardPlaneToolControls.h" + +class QmitkPointListWidget; + +/*! +@brief QmitkImageCropperView +\warning This class is not yet documented. Use "git blame" and ask the author to provide basic documentation. + +\sa QmitkFunctionality +\ingroup ${plugin_target}_internal +*/ +class QmitkStandardPlaneTool : public QmitkAbstractView +{ + // this is needed for all Qt objects that should have a Qt meta-object + // (everything that derives from QObject and wants to have signal/slots) +private: + + Q_OBJECT + +public: + + QmitkStandardPlaneTool(QObject *parent = 0); + + virtual ~QmitkStandardPlaneTool(); + + static const std::string VIEW_ID; + + virtual void CreateQtPartControl(QWidget *parent); + + + ///// \brief Handles selection of data nodes + virtual void OnSelectionChanged(const mitk::DataNode* node); + ///// \brief Handles selection of data nodes + virtual void OnSelectionChanged(berry::IWorkbenchPart::Pointer part, const QList &nodes); + + protected slots: + + /// \brief Handles selection of images through the comboboxes + void OnComboBoxSelectionChanged(const mitk::DataNode* node); + /// \brief Handles selection of images through the comboboxes + //void OnAtlasComboBoxSelectionChanged(const mitk::DataNode* node); + + /// \brief Called for resetting the recorded time in the workingNode (reference data acquisition) + void ResetTimeProperty(); + /// \brief Called for saving the position of the planes (e.g. for saving reference data) + void AddPositionNode(); + + /// \brief Load the position information from the position nodes after loading an MITK project + void LoadPositionsFromNodes(); + /// \brief Load the position information from the selected position node + void LoadPositionFromNode(); + + /// \brief Show or hide tool view + void ToggleStandardPlaneToolView(); + + /// \brief Calculate and save intersection lines between the planes to file + void SaveIntersectionFromNodeToFile(); + + void LoadPositionFromFile(); + + void LoadPositionFromJSON(); + void LoadPlaneOrientationFromFile(); + + void ShowSilhouette(int projNum, vtkSmartPointer projMatrix, vtkSmartPointer surface,QString fileName, int showSilhouette, double offset, int sil, int shape); + + void OnSaveLateralSliceClicked(); + void OnSaveMortiseSliceClicked(); + void OnSaveAPSliceClicked(); + + void OnShowLateralSliceClicked(); + void OnShowMortiseSliceClicked(); + void OnShowAPSliceClicked(); + void OnRegistrationClicked(); + + void ToggleProjectionExtractionToolView(); + void ToggleSilhouetteToolView(); + void InitProjectionExtractionToolWidget(); + void InitStandardPlaneToolWidget(); + + void OnLoadSavedProjections(); + void OnExportProjections(); + void OnExtractProjectionsClicked(); + void OnManualRegistrationClicked(); + + void OnExtractSilhouetteClicked(); + void InitSilhouetteToolWidget(); + + void OnSliderValueChanged(); + void OnExtractContourAndSlicesClicked(); + +protected: + + void ShowAxialSlice(int sliceIndex); + + /// \brief Called for loading the plane position from a data node + void LoadPositionNode(const mitk::DataNode *node, int id); + + /// \brief Old data node selection (used for memory of last selected node) + mitk::DataNode::Pointer m_OldNode; + /// \brief the current selected data node + mitk::DataNode::Pointer m_WorkingNode; + + mitk::DataNode::Pointer m_Projection1; + mitk::DataNode::Pointer m_Projection2; + mitk::DataNode::Pointer m_Silhouette1; + mitk::DataNode::Pointer m_Silhouette2; + mitk::Surface::Pointer m_Surface; + mitk::Image::Pointer m_Segmentation; + + std::vector> m_ProjectionMat; + int m_ProjNum1; + int m_ProjNum2; + + + + int m_CurrentMortiseID; + int m_CurrentLateralID; + int m_CurrentAPID; + + bool m_showSPV; + bool m_showPEV; + bool m_showSV; + virtual void SetFocus(); + + Ui::QmitkStandardPlaneToolControls m_Controls; +}; + + +#endif // QmitkStandardPlaneTool_h \ No newline at end of file diff --git a/Plugins/org.mitk.gui.qt.standardplanetool/src/internal/QmitkStandardPlaneToolControls.ui b/Plugins/org.mitk.gui.qt.standardplanetool/src/internal/QmitkStandardPlaneToolControls.ui new file mode 100644 index 0000000000..ffc718a6d7 --- /dev/null +++ b/Plugins/org.mitk.gui.qt.standardplanetool/src/internal/QmitkStandardPlaneToolControls.ui @@ -0,0 +1,567 @@ + + + QmitkStandardPlaneToolControls + + + + 0 + 0 + 354 + 790 + + + + + 0 + 0 + + + + + 290 + 790 + + + + + 16777215 + 16777215 + + + + QmitkTemplate + + + + + 10 + 10 + 80 + 80 + + + + + 0 + 0 + + + + Standard + Plane + Tool + + + + + + 130 + 10 + 80 + 80 + + + + + 0 + 0 + + + + Projection + Extraction + Tool + + + + + + 10 + 100 + 311 + 31 + + + + + + + false + + + + 0 + 0 + + + + Select an image: + + + Qt::AlignCenter + + + + + + + false + + + + 0 + 0 + + + + + + + + + + 10 + 150 + 311 + 621 + + + + + + + + 0 + 0 + + + + Standard Plane Tool + + + + + + + + Set New Position + + + + + + + Load Current Position + + + + + + + Load Positions From File + + + + + + + false + + + Save All Positionsnodes To File + + + + + + + + + + + + + 0 + 0 + + + + Projection Extraction Tool + + + + + 10 + 30 + 291 + 181 + + + + + + + Load Saved Projections + + + + + + + Qt::Horizontal + + + + + + + 0 + + + + + Lateral + + + + + + + + 0 + 0 + + + + Show L + + + + + + + + 0 + 0 + + + + Save L + + + + + + + + 0 + 0 + + + + QFrame::NoFrame + + + + + + + + + Qt::Horizontal + + + + + + + 0 + + + + + Mortise + + + + + + + + 0 + 0 + + + + Show M + + + + + + + + 0 + 0 + + + + Save M + + + + + + + + 0 + 0 + + + + QFrame::NoFrame + + + + + + + + + Qt::Horizontal + + + + + + + 0 + + + + + AP + + + + + + + + 0 + 0 + + + + Show AP + + + + + + + + 0 + 0 + + + + Save AP + + + + + + + + 0 + 0 + + + + QFrame::NoFrame + + + + + + + + + Qt::Horizontal + + + + + + + Export Saved Projections + + + + + + + Qt::Vertical + + + + 20 + 40 + + + + + + + + + + + + + 0 + 0 + + + + Silhouette Tool + + + + + 10 + 20 + 231 + 232 + + + + + + + Extract Silhouettes + + + + + + + Start Registration + + + + + + + 4 + + + + + Start Manual Registration + + + + + + + Qt::Horizontal + + + + + + + Qt::Horizontal + + + + + + + Qt::Horizontal + + + + + + + Qt::Horizontal + + + + + + + Qt::Horizontal + + + + + + + Qt::Horizontal + + + + + + + Extract Slices and Masks + + + + + + + + + Qt::Vertical + + + + 20 + 40 + + + + + + + + + + + + + + 240 + 10 + 80 + 80 + + + + + 0 + 0 + + + + Silhouette + Tool + + + + + + + QmitkDataStorageComboBox + QComboBox +
QmitkDataStorageComboBox.h
+ + + + + diff --git a/Plugins/org.mitk.gui.qt.standardplanetool/src/internal/QmitkStandardPlaneTool_PC.cpp b/Plugins/org.mitk.gui.qt.standardplanetool/src/internal/QmitkStandardPlaneTool_PC.cpp new file mode 100644 index 0000000000..16ffa322c2 --- /dev/null +++ b/Plugins/org.mitk.gui.qt.standardplanetool/src/internal/QmitkStandardPlaneTool_PC.cpp @@ -0,0 +1,2024 @@ + +/*=================================================================== + +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 "QmitkStandardPlaneTool.h" + +#include "Poco/DirectoryIterator.h" + +#include +#include +#include +#include +#include + +// MITK +#include +#include +#include +#include +#include "mitkNodePredicateProperty.h" +#include "mitkNodePredicateDataType.h" +#include +#include +#include +#include "mitkPlanePositionManager.h" +#include "mitkPlanarCircle.h" +#include +#include "mitkVtkRepresentationProperty.h" +#include "mitkImageSliceSelector.h" +#include "mitkConvert2Dto3DImageFilter.h" +#include "mitk2D3DShapeRegistration.h" +#include +#include +#include +#include +#include +#include +//#include + +#include "org_mitk_gui_qt_standardplanetool_Activator.h" +#include "vtkSilhouette.h" +#include "vtkPoints.h" + +// Qmitk +#include "QmitkStdMultiWidget.h" + +// Qt +#include +#include +#include + + +// VTK (for testing evaluation methods) +#include +#include +#include +#include + +#include +#include "vtkDoubleArray.h" +#include "vtkFloatArray.h" +#include "vtkPointData.h" +#include "vtkProperty.h" +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include "itkImage.h" +#include +#include +#include +#include "itkRescaleIntensityImageFilter.h" + +//Boost +#include +#include + + +const std::string QmitkStandardPlaneTool::VIEW_ID = "org.mitk.views.qmitkstandardplanetool"; + +QmitkStandardPlaneTool::QmitkStandardPlaneTool(QObject *parent) : +m_OldNode(nullptr), +m_WorkingNode(nullptr), +m_Projection1(nullptr), +m_Projection2(nullptr), +m_Silhouette1(nullptr), +m_Silhouette2(nullptr), +m_Segmentation(nullptr), +m_Surface(nullptr), +m_showSPV(0), +m_showPEV(0), +m_showSV(0), +m_CurrentMortiseID(-1), +m_CurrentLateralID(-1), +m_CurrentAPID(-1) +{ +} + +QmitkStandardPlaneTool::~QmitkStandardPlaneTool() +{ +} + +void QmitkStandardPlaneTool::CreateQtPartControl(QWidget *parent) +{ + // create GUI widgets from the Qt Designer's .ui file + m_Controls.setupUi(parent); + + ////invalidate timer for first measurement of plane-adjustment-time (user tab) + //m_timer.invalidate(); + + // create GUI widgets from the Qt Designer's .ui file + m_Controls.patientImageSelector->SetDataStorage(this->GetDataStorage()); + m_Controls.patientImageSelector->SetPredicate(mitk::NodePredicateAnd::New( + mitk::TNodePredicateDataType::New(), + mitk::NodePredicateNot::New(mitk::NodePredicateProperty::New("helper object")))); + m_WorkingNode= m_Controls.patientImageSelector->GetSelectedNode(); + + //// create signal/slot connections of comboboxes + connect(m_Controls.patientImageSelector, SIGNAL(OnSelectionChanged(const mitk::DataNode*)), + this, SLOT(OnComboBoxSelectionChanged(const mitk::DataNode*))); + + + connect(m_Controls.pB_StandardPlaneTool, SIGNAL(clicked()), this, SLOT(ToggleStandardPlaneToolView())); + connect(m_Controls.pB_ProjectionExtractionTool, SIGNAL(clicked()), this, SLOT(ToggleProjectionExtractionToolView())); + connect(m_Controls.pB_SilhouetteTool, SIGNAL(clicked()), this, SLOT(ToggleSilhouetteToolView())); + + + m_Controls.patientImageSelector->setEnabled(true); + + this->InitProjectionExtractionToolWidget(); + this->InitStandardPlaneToolWidget(); + this->InitSilhouetteToolWidget(); +} + + +void QmitkStandardPlaneTool::InitSilhouetteToolWidget() +{ + connect(m_Controls.pB_ExtractSilhouettes, SIGNAL(clicked()), this, SLOT(OnExtractSilhouetteClicked())); + connect(m_Controls.pB_StartRegistration, SIGNAL(clicked()), this, SLOT( OnRegistrationClicked())); + connect(m_Controls.pB_StartManualRegistration, SIGNAL(clicked()), this, SLOT(OnManualRegistrationClicked())); + connect(m_Controls.pB_ExtractContourAndSlices, SIGNAL(clicked()), this, SLOT(OnExtractContourAndSlicesClicked())); + + m_Controls.sB_Tx->setRange(0, 100); + m_Controls.sB_Ty->setRange(0, 100); + m_Controls.sB_Tz->setRange(0, 100); + m_Controls.sB_alpha->setRange(0, 100); + m_Controls.sB_beta->setRange(0, 100); + m_Controls.sB_gamma->setRange(0, 100); + + + connect(m_Controls.sB_Tx, SIGNAL(sliderMoved(int)), this, SLOT(OnSliderValueChanged())); + connect(m_Controls.sB_Ty, SIGNAL(sliderMoved(int)), this, SLOT(OnSliderValueChanged())); + connect(m_Controls.sB_Tz, SIGNAL(sliderMoved(int)), this, SLOT(OnSliderValueChanged())); + connect(m_Controls.sB_alpha, SIGNAL(sliderMoved(int)), this, SLOT(OnSliderValueChanged())); + connect(m_Controls.sB_beta, SIGNAL(sliderMoved(int)), this, SLOT(OnSliderValueChanged())); + connect(m_Controls.sB_gamma, SIGNAL(sliderMoved(int)), this, SLOT(OnSliderValueChanged())); + + + + m_Controls.gB_SilhouetteTool->hide(); +} + +void QmitkStandardPlaneTool::OnSaveLateralSliceClicked() +{ + QmitkRenderWindow* RenderWindowAxial = this->GetRenderWindowPart()->GetQmitkRenderWindow("axial"); + m_CurrentLateralID = RenderWindowAxial->GetSliceNavigationController()->GetSlice()->GetPos(); + mitk::DataNode* imageNode(m_Controls.patientImageSelector->GetSelectedNode()); + if (imageNode == nullptr) + return; + imageNode->SetIntProperty("StandardProjection.Lateral", m_CurrentLateralID); + this->m_Controls.lcd_L->display(m_CurrentLateralID); +} +void QmitkStandardPlaneTool::OnSaveMortiseSliceClicked() +{ + QmitkRenderWindow* RenderWindowAxial = this->GetRenderWindowPart()->GetQmitkRenderWindow("axial"); + m_CurrentMortiseID = RenderWindowAxial->GetSliceNavigationController()->GetSlice()->GetPos(); + mitk::DataNode* imageNode(m_Controls.patientImageSelector->GetSelectedNode()); + if (imageNode == nullptr) + return; + imageNode->SetIntProperty("StandardProjection.Mortise", m_CurrentMortiseID); + this->m_Controls.lcd_M->display(m_CurrentMortiseID); +} +void QmitkStandardPlaneTool::OnSaveAPSliceClicked() +{ + QmitkRenderWindow* RenderWindowAxial = this->GetRenderWindowPart()->GetQmitkRenderWindow("axial"); + m_CurrentAPID = RenderWindowAxial->GetSliceNavigationController()->GetSlice()->GetPos(); + mitk::DataNode* imageNode(m_Controls.patientImageSelector->GetSelectedNode()); + if (imageNode == nullptr) + return; + imageNode->SetIntProperty("StandardProjection.AP", m_CurrentAPID); + this->m_Controls.lcd_AP->display(m_CurrentAPID); +} + + +void QmitkStandardPlaneTool::OnShowLateralSliceClicked() +{ + if (m_CurrentLateralID >= 0) + this->ShowAxialSlice(m_CurrentLateralID); +} +void QmitkStandardPlaneTool::OnShowMortiseSliceClicked() +{ + if (m_CurrentMortiseID>=0) + this->ShowAxialSlice(m_CurrentMortiseID); +} +void QmitkStandardPlaneTool::OnShowAPSliceClicked() +{ + if (m_CurrentAPID >= 0) + this->ShowAxialSlice(m_CurrentAPID); +} + +void QmitkStandardPlaneTool::OnLoadSavedProjections() +{ + mitk::DataNode* imageNode(m_Controls.patientImageSelector->GetSelectedNode()); + if (imageNode == nullptr) + return; + + if (imageNode->GetIntProperty("StandardProjection.AP", m_CurrentAPID)) + { + imageNode->GetIntProperty("StandardProjection.AP", m_CurrentAPID); + this->m_Controls.lcd_AP->display(m_CurrentAPID); + } + else + { + this->m_Controls.lcd_AP->display(0); + m_CurrentAPID = -1; + } + + if (imageNode->GetIntProperty("StandardProjection.Lateral", m_CurrentLateralID)) + { + imageNode->GetIntProperty("StandardProjection.Lateral", m_CurrentLateralID); + this->m_Controls.lcd_L->display(m_CurrentLateralID); + } + else + { + this->m_Controls.lcd_L->display(0); + m_CurrentLateralID = -1; + } + if (imageNode->GetIntProperty("StandardProjection.Mortise", m_CurrentMortiseID)) + { + imageNode->GetIntProperty("StandardProjection.Mortise", m_CurrentMortiseID); + this->m_Controls.lcd_M->display(m_CurrentMortiseID); + } + else + { + m_CurrentMortiseID = -1; + this->m_Controls.lcd_M->display(0); + } + + +} + + + + +void ExtractSlice(mitk::DataStorage* dataStorage, mitk::Image::Pointer image, int slice, QString fileName, bool setFirst, bool save, bool show = 0) +{ + + MITK_INFO << slice; + auto selector = mitk::ImageSliceSelector::New(); + selector->SetSliceNr(slice); + selector->SetInput(image); + selector->Update(); + mitk::Image::Pointer projection1 = selector->GetOutput(); + + if (setFirst) + { + mitk::Point3D origin = projection1->GetGeometry()->GetOrigin(); + origin[2] = 0.0; + projection1->GetGeometry()->SetOrigin(origin[0]); + projection1->Modified(); + } + + auto convert2DTo3DImageFilter = mitk::Convert2Dto3DImageFilter::New(); + convert2DTo3DImageFilter->SetInput(projection1); + convert2DTo3DImageFilter->Update(); + + if (show) + { + auto projection = mitk::DataNode::New(); + projection->SetProperty("name", mitk::StringProperty::New("ImageSource")); + projection->SetData(projection1); + projection->SetColor(1.0, 1.0, 1.0); + projection->SetVisibility(false); + projection->SetName(fileName.toStdString()); + dataStorage->Add(projection); + } + + if (save) + { + if (fileName == nullptr) + return; + mitk::IOUtil::Save(convert2DTo3DImageFilter->GetOutput(), fileName.toStdString() + ".nrrd"); + MITK_INFO << "Image Saved"; + } +} + +void QmitkStandardPlaneTool::OnExportProjections() +{ + mitk::DataNode* imageNode(m_Controls.patientImageSelector->GetSelectedNode()); + if (imageNode == nullptr) + return; + + mitk::Image::Pointer image = dynamic_cast(imageNode->GetData()); + + QString workingDir("C:/Users/thomass.AD/Desktop/StandardProjektionen/"); + QString LateralFileName, MortiseFileName, APFileName; + + ExtractSlice(this->GetDataStorage(), image, m_CurrentAPID, APFileName.append(workingDir).append(imageNode->GetName().c_str()).append("_").append(QString::number(m_CurrentAPID)).append("_AP"), true, true); + ExtractSlice(this->GetDataStorage(), image, m_CurrentLateralID, LateralFileName.append(workingDir).append(imageNode->GetName().c_str()).append("_").append(QString::number(m_CurrentLateralID)).append("_Lateral"), true, true); + ExtractSlice(this->GetDataStorage(), image, m_CurrentMortiseID, MortiseFileName.append(workingDir).append(imageNode->GetName().c_str()).append("_").append(QString::number(m_CurrentMortiseID)).append("_Mortise"), true, true); + imageNode->SetName(imageNode->GetName().append("Extracted")); +} + + +void QmitkStandardPlaneTool::InitProjectionExtractionToolWidget() +{ + connect(m_Controls.pb_LoadSavedProjections, SIGNAL(clicked()), this, SLOT(OnLoadSavedProjections())); //Load Position + connect(m_Controls.pB_ExportProjections, SIGNAL(clicked()), this, SLOT(OnExportProjections())); //Export Position + + connect(m_Controls.pB_SaveMortise, SIGNAL(clicked()), this, SLOT(OnSaveMortiseSliceClicked())); + connect(m_Controls.pB_SaveAP, SIGNAL(clicked()), this, SLOT(OnSaveAPSliceClicked())); + connect(m_Controls.pB_SaveLateral, SIGNAL(clicked()), this, SLOT(OnSaveLateralSliceClicked())); + + connect(m_Controls.pB_ShowMortise, SIGNAL(clicked()), this, SLOT(OnShowMortiseSliceClicked())); + connect(m_Controls.pB_ShowAP, SIGNAL(clicked()), this, SLOT(OnShowAPSliceClicked())); + connect(m_Controls.pB_ShowLateral, SIGNAL(clicked()), this, SLOT(OnShowLateralSliceClicked())); + + m_Controls.gB_ProjectionExtractionTool->hide(); +} + +void QmitkStandardPlaneTool::InitStandardPlaneToolWidget() +{ + connect(m_Controls.pB_savePosition, SIGNAL(clicked()), this, SLOT(AddPositionNode())); + connect(m_Controls.pB_loadPositions, SIGNAL(clicked()), this, SLOT(LoadPositionFromNode())); + connect(m_Controls.pB_LoadFromFile, SIGNAL(clicked()), this, SLOT(LoadPositionFromFile())); + connect(m_Controls.pB_SaveToFile, SIGNAL(clicked()), this, SLOT(SaveIntersectionFromNodeToFile())); + + m_Controls.gB_StandardPlaneTool->hide(); +} + +void QmitkStandardPlaneTool::ToggleStandardPlaneToolView() +{ + if (!m_showSPV) + m_Controls.gB_StandardPlaneTool->show(); + else + m_Controls.gB_StandardPlaneTool->hide(); + + m_showSPV = !m_showSPV; +} + +void QmitkStandardPlaneTool::ToggleProjectionExtractionToolView() +{ + if (!m_showPEV) + m_Controls.gB_ProjectionExtractionTool->show(); + else + m_Controls.gB_ProjectionExtractionTool->hide(); + + m_showPEV = !m_showPEV; +} + +void QmitkStandardPlaneTool::ToggleSilhouetteToolView() +{ + if (!m_showSV) + m_Controls.gB_SilhouetteTool->show(); + else + m_Controls.gB_SilhouetteTool->hide(); + + m_showSV = !m_showSV; +} + + +void QmitkStandardPlaneTool::ShowAxialSlice(int sliceIndex) +{ + QmitkRenderWindow* RenderWindowAxial = this->GetRenderWindowPart()->GetQmitkRenderWindow("axial"); + RenderWindowAxial->GetSliceNavigationController()->GetSlice()->SetPos(sliceIndex); + RenderWindowAxial->GetSliceNavigationController()->SendSlice(); + + mitk::RenderingManager::GetInstance()->RequestUpdateAll(); +} + +void QmitkStandardPlaneTool::OnSelectionChanged(const mitk::DataNode* node) +{ + mitk::DataNode* selectedNode = const_cast(node); + QList nodes; + nodes.push_back(selectedNode); + //TODO because a part is needed a part is given (this should be changed to something else) + + + mitk::NodePredicateDataType::Pointer isImage = mitk::NodePredicateDataType::New("Image"); + mitk::DataStorage::SetOfObjects::ConstPointer images = this->GetDataStorage()->GetSubset(isImage); + if (!images->empty()) + { + //iterate over all child nodes with NodePredicateProperty + for (mitk::DataStorage::SetOfObjects::const_iterator iter = images->begin(); iter != images->end(); ++iter) + { + if ((*iter) != node) + (*iter)->SetVisibility(false); + } + } + + this->OnSelectionChanged(this->GetSite()->GetPart(), nodes); + mitk::RenderingManager::GetInstance()->RequestUpdateAll(); + this->OnLoadSavedProjections(); +} + +void QmitkStandardPlaneTool::OnSelectionChanged(berry::IWorkbenchPart::Pointer /*source*/, + const QList &nodes) +{ + + if (!nodes.isEmpty()) + { + mitk::DataNode::Pointer referenceData = nodes.front(); + if (referenceData.IsNull()) return; + + if (m_OldNode.IsNull()) + { + m_OldNode = referenceData; + } + else if (this->GetDataStorage()->GetNamedNode(m_WorkingNode->GetName())) + { + m_OldNode = m_WorkingNode; + } + else + { + m_OldNode = referenceData; + } + m_WorkingNode = referenceData; + //set the working node to visible + m_WorkingNode->SetProperty("visible", mitk::BoolProperty::New(true)); + + //set comboBox to reference image + disconnect(m_Controls.patientImageSelector, SIGNAL(OnSelectionChanged(const mitk::DataNode*)), + this, SLOT(OnComboBoxSelectionChanged(const mitk::DataNode*))); + + m_Controls.patientImageSelector->setCurrentIndex(m_Controls.patientImageSelector->Find(referenceData)); + + connect(m_Controls.patientImageSelector, SIGNAL(OnSelectionChanged(const mitk::DataNode*)), + this, SLOT(OnComboBoxSelectionChanged(const mitk::DataNode*))); + + + mitk::NodePredicateProperty::Pointer isMarker = mitk::NodePredicateProperty::New("isPositionNode", mitk::BoolProperty::New(true)); + //load all nodes that are children of the working node in data storage + mitk::DataStorage::SetOfObjects::ConstPointer markers = this->GetDataStorage()->GetDerivations(m_WorkingNode, isMarker); + if (!markers->empty()) + { + //iterate over all child nodes with NodePredicateProperty + for (mitk::DataStorage::SetOfObjects::const_iterator iter = markers->begin(); iter != markers->end(); ++iter) + { + int markerId; + (*iter)->GetIntProperty("Plane Position ID", markerId); + this->LoadPositionNode((*iter), markerId); + } + } + + + + ////check if time node exists and if not create node (start timmer) until new node is selected + ////checkbox that is activated by the user when everything is loaded and time recording starts + //if (m_Controls.checkBoxMeasureTime->isChecked()) + //{ + // int time = 0; + // if (m_timer.isValid() && (!m_OldNode->GetIntProperty("time elapsed", time) || time == 0)) m_OldNode->SetIntProperty("time elapsed", m_timer.elapsed()); + // m_timer.restart(); + //} + + } +} + +void QmitkStandardPlaneTool::OnComboBoxSelectionChanged(const mitk::DataNode* node) +{ + //if at least one node remains in the data manager + if (node != NULL) + m_Controls.patientImageSelector->setEnabled(true); + else + m_Controls.patientImageSelector->setEnabled(false); + + this->OnSelectionChanged(node); +} + + +void QmitkStandardPlaneTool::AddPositionNode() +{ + MITK_INFO << m_Controls.patientImageSelector->currentIndex(); + if (m_Controls.patientImageSelector->currentIndex() <0) + return; + //we need to save all planes... + for (int i = 1; i <= 3; i++) + { + + std::stringstream multiWidgetStream; + multiWidgetStream << "stdmulti.widget"; + multiWidgetStream << i; + std::string multiWidgetString = multiWidgetStream.str(); + + mitk::BaseRenderer *m_Rendereri = mitk::BaseRenderer::GetByName(multiWidgetString.c_str()); + + if (m_Rendereri) + { + const mitk::PlaneGeometry* plane = dynamic_cast< const mitk::SlicedGeometry3D*>(m_Rendereri->GetSliceNavigationController()->GetCurrentGeometry3D())->GetPlaneGeometry(0); + + //Getting Service + ctkPluginContext* context = mitk::org_mitk_gui_qt_standardplanetool_Activator::GetContext(); + ctkServiceReference ppmRef = context->getServiceReference(); + mitk::PlanePositionManagerService* service = context->getService(ppmRef); + //unsigned int size = service->GetNumberOfPlanePositions(); + //context->ungetService(ppmRef); + + //node predicate property to identify all childnodes for planepositioning + mitk::NodePredicateProperty::Pointer isMarker = mitk::NodePredicateProperty::New("isPositionNode", mitk::BoolProperty::New(true)); + //removes all nodes that are children of the working node in data storage + mitk::DataStorage::SetOfObjects::ConstPointer markers = this->GetDataStorage()->GetDerivations(m_WorkingNode, isMarker); + if (!markers->empty() && i == 1) + { + this->GetDataStorage()->Remove(markers); + //remove all entries of these nodes in the service + for (mitk::DataStorage::SetOfObjects::const_iterator iter = markers->begin(); iter != markers->end(); ++iter) + { + int markerId; + (*iter)->GetIntProperty("Plane Position ID", markerId); + //get all position nodes in the data manager + mitk::NodePredicateProperty::Pointer isMarker = mitk::NodePredicateProperty::New("isPositionNode", mitk::BoolProperty::New(true)); + mitk::DataStorage::SetOfObjects::ConstPointer positionNodes = this->GetDataStorage()->GetSubset(isMarker); + bool duplicate = false; + int duplicateId; + //iterate over the found nodes + for (unsigned int i = 0; i < positionNodes->size(); i++) //TODO make this efficient! + { + positionNodes->at(i)->GetIntProperty("Plane Position ID", duplicateId); + if (markerId == duplicateId) duplicate = true; + } + if (!duplicate) service->RemovePlanePosition(markerId); + } + } + + unsigned int id = service->AddNewPlanePosition(plane, m_Rendereri->GetSliceNavigationController()->GetSlice()->GetPos()); + + //construct node name with some property-coding + std::stringstream nodeNameStream; + nodeNameStream << "Planeposition "; + //insert name extension axial/sagittal/coronal + switch (i) + { + case 1: nodeNameStream << "axial"; + break; + case 2: nodeNameStream << "sagittal"; + break; + case 3: nodeNameStream << "coronal"; + } + + std::string nameString = nodeNameStream.str(); + + mitk::PlanarCircle::Pointer positionMarker = mitk::PlanarCircle::New(); + mitk::Point2D p1; + plane->Map(plane->GetCenter(), p1); + mitk::Point2D p2 = p1; + p2[0] -= plane->GetSpacing()[0]; + p2[1] -= plane->GetSpacing()[1]; + positionMarker->PlaceFigure(p1); + positionMarker->SetCurrentControlPoint(p1); + positionMarker->SetPlaneGeometry(const_cast(plane)); + + // the current selected node in the image selector + mitk::DataNode* imageNode(m_Controls.patientImageSelector->GetSelectedNode()); + + mitk::DataNode::Pointer planePositionNode = mitk::DataNode::New(); + + planePositionNode->SetData(positionMarker); + //TODO workaround for loading the plane positions (filling the service) + planePositionNode->SetIntProperty("sliceIndex", m_Rendereri->GetSliceNavigationController()->GetSlice()->GetPos()); + //Saving the ID of the plane position as a property of the node + planePositionNode->SetIntProperty("Plane Position ID", id); + planePositionNode->SetProperty("name", mitk::StringProperty::New(nameString.c_str())); + planePositionNode->SetProperty("isPositionNode", mitk::BoolProperty::New(true)); + planePositionNode->SetProperty("includeInBoundingBox", mitk::BoolProperty::New(false)); + planePositionNode->SetBoolProperty("PlanarFigureInitializedWindow", true, m_Rendereri); + //TODO if the position nodes are set to helper object the loading does not work + planePositionNode->SetProperty("helper object", mitk::BoolProperty::New(false)); + //unchecks the planeposition nodes in datamanager + planePositionNode->SetProperty("visible", mitk::BoolProperty::New(false)); + + if (plane) + { + this->GetDataStorage()->Add(planePositionNode, imageNode); + } + } + } +} + +void QmitkStandardPlaneTool::LoadPositionNode(const mitk::DataNode *node, int id) +{ + MITK_INFO << m_Controls.patientImageSelector->currentIndex(); + if (m_Controls.patientImageSelector->currentIndex() <0) + return; + + QmitkRenderWindow* selectedRenderWindow = 0; + QmitkRenderWindow* RenderWindow1 = this->GetRenderWindowPart()->GetQmitkRenderWindow("axial"); + QmitkRenderWindow* RenderWindow2 = this->GetRenderWindowPart()->GetQmitkRenderWindow("sagittal"); + QmitkRenderWindow* RenderWindow3 = this->GetRenderWindowPart()->GetQmitkRenderWindow("coronal"); + + bool PlanarFigureInitializedWindow = false; + + // find initialized renderwindow + if (node->GetBoolProperty("PlanarFigureInitializedWindow", + PlanarFigureInitializedWindow, RenderWindow1->GetRenderer())) + { + selectedRenderWindow = RenderWindow1; + } + if (!selectedRenderWindow && node->GetBoolProperty( + "PlanarFigureInitializedWindow", PlanarFigureInitializedWindow, + RenderWindow2->GetRenderer())) + { + selectedRenderWindow = RenderWindow2; + } + if (!selectedRenderWindow && node->GetBoolProperty( + "PlanarFigureInitializedWindow", PlanarFigureInitializedWindow, + RenderWindow3->GetRenderer())) + { + selectedRenderWindow = RenderWindow3; + } + + if (selectedRenderWindow) + { + { + ctkPluginContext* context = mitk::org_mitk_gui_qt_standardplanetool_Activator::GetContext(); + ctkServiceReference ppmRef = context->getServiceReference(); + mitk::PlanePositionManagerService* service = context->getService(ppmRef); + selectedRenderWindow->GetSliceNavigationController()->ExecuteOperation(service->GetPlanePosition(id)); + context->ungetService(ppmRef); + } + mitk::RenderingManager::GetInstance()->RequestUpdateAll(); + } +} + +void QmitkStandardPlaneTool::ResetTimeProperty() +{ + //reset the saved time for plane adjustment of the working node + // m_WorkingNode->SetIntProperty("time elapsed", 0); +} + + +void QmitkStandardPlaneTool::LoadPositionsFromNodes() +{ + MITK_INFO << m_Controls.patientImageSelector->currentIndex(); + if (m_Controls.patientImageSelector->currentIndex() <0) + return; + + //get the planePositionManager Service + ctkPluginContext* context = mitk::org_mitk_gui_qt_standardplanetool_Activator::GetContext(); + ctkServiceReference ppmRef = context->getServiceReference(); + mitk::PlanePositionManagerService* service = context->getService(ppmRef); + service->RemoveAllPlanePositions(); + + int sliceIndex = 0; + mitk::PlanarCircle::Pointer positionMarker = mitk::PlanarCircle::New(); + int id; + //get position nodes in the data manager and add them with the planepositionmanger service + mitk::DataNode::Pointer coronalPosNode = this->GetDataStorage()->GetNamedNode("Planeposition coronal"); + if (coronalPosNode != nullptr){ + coronalPosNode->GetIntProperty("sliceIndex", sliceIndex); + positionMarker = dynamic_cast (coronalPosNode->GetData()); + const mitk::Geometry2D::ConstPointer corPlane = positionMarker->GetPlaneGeometry(); + id = service->AddNewPlanePosition(corPlane, sliceIndex); + coronalPosNode->SetIntProperty("Plane Position ID", id); + MITK_INFO << "/////////////////Marker Pointer\\\\\\\\\\\\\\\\ "; + MITK_INFO << positionMarker; + } + mitk::DataNode::Pointer axialPosNode = this->GetDataStorage()->GetNamedNode("Planeposition axial"); + if (axialPosNode != nullptr){ + + axialPosNode->GetIntProperty("sliceIndex", sliceIndex); + positionMarker = dynamic_cast (axialPosNode->GetData()); + const mitk::Geometry2D::ConstPointer axiPlane = positionMarker->GetPlaneGeometry(); + id = service->AddNewPlanePosition(axiPlane, sliceIndex); + axialPosNode->SetIntProperty("Plane Position ID", id); + } + mitk::DataNode::Pointer sagittalPosNode = this->GetDataStorage()->GetNamedNode("Planeposition sagittal"); + if (sagittalPosNode != nullptr){ + sagittalPosNode->GetIntProperty("sliceIndex", sliceIndex); + positionMarker = dynamic_cast (sagittalPosNode->GetData()); + const mitk::Geometry2D::ConstPointer sagPlane = positionMarker->GetPlaneGeometry(); + id = service->AddNewPlanePosition(sagPlane, sliceIndex); + sagittalPosNode->SetIntProperty("Plane Position ID", id); + } + + context->ungetService(ppmRef); +} + +void QmitkStandardPlaneTool::LoadPositionFromNode() +{ + MITK_INFO << m_Controls.patientImageSelector->currentIndex(); + if (m_Controls.patientImageSelector->currentIndex() <0) + return; + + //get the planePositionManager Service + ctkPluginContext* context = mitk::org_mitk_gui_qt_standardplanetool_Activator::GetContext(); + ctkServiceReference ppmRef = context->getServiceReference(); + mitk::PlanePositionManagerService* service = context->getService(ppmRef); + + int sliceIndex = 0; + mitk::PlanarCircle::Pointer positionMarker = mitk::PlanarCircle::New(); + + mitk::NodePredicateProperty::Pointer isMarker = mitk::NodePredicateProperty::New("isPositionNode", mitk::BoolProperty::New(true)); + mitk::DataStorage::SetOfObjects::ConstPointer markers; + + //load all nodes that are children of the working node in data storage + mitk::DataNode* imageNode(m_Controls.patientImageSelector->GetSelectedNode()); + markers = this->GetDataStorage()->GetDerivations(imageNode, isMarker); + MITK_INFO << "------------------------------------"; + MITK_INFO << "Process Image: " << imageNode->GetName(); + + if (!markers->empty()) + { + //iterate over all child nodes with NodePredicateProperty + for (mitk::DataStorage::SetOfObjects::const_iterator iter = markers->begin(); iter != markers->end(); ++iter) + { + MITK_INFO << "Marker" << std::endl; + int markerId; + std::string nodeName = (*iter)->GetName(); + //get position nodes in the data manager and add them with the planepositionmanger service + if (nodeName == "Planeposition coronal") + { + mitk::DataNode::Pointer coronalPosNode = (*iter); + coronalPosNode->GetIntProperty("sliceIndex", sliceIndex); + positionMarker = dynamic_cast (coronalPosNode->GetData()); + mitk::PlaneGeometry::ConstPointer corPlane = positionMarker->GetPlaneGeometry(); + markerId = service->AddNewPlanePosition(corPlane, sliceIndex); + coronalPosNode->SetIntProperty("Plane Position ID", markerId); + } + else if (nodeName == "Planeposition axial") + { + mitk::DataNode::Pointer axialPosNode = (*iter); + axialPosNode->GetIntProperty("sliceIndex", sliceIndex); + positionMarker = dynamic_cast (axialPosNode->GetData()); + mitk::PlaneGeometry::ConstPointer axiPlane = positionMarker->GetPlaneGeometry(); + markerId = service->AddNewPlanePosition(axiPlane, sliceIndex); + axialPosNode->SetIntProperty("Plane Position ID", markerId); + } + else if (nodeName == "Planeposition sagittal") + { + mitk::DataNode::Pointer sagittalPosNode = (*iter); + sagittalPosNode->GetIntProperty("sliceIndex", sliceIndex); + positionMarker = dynamic_cast (sagittalPosNode->GetData()); + mitk::PlaneGeometry::ConstPointer sagPlane = positionMarker->GetPlaneGeometry(); + markerId = service->AddNewPlanePosition(sagPlane, sliceIndex); + sagittalPosNode->SetIntProperty("Plane Position ID", markerId); + } + + (*iter)->GetIntProperty("Plane Position ID", markerId); + this->LoadPositionNode((*iter), markerId); + mitk::RenderingManager::GetInstance()->RequestUpdateAll(); + } + } + +} + + +void QmitkStandardPlaneTool::LoadPositionFromFile() +{ + + QString fileName = QFileDialog::getOpenFileName(NULL, + tr("Open Scene File with Positions"),0, tr("MITK scene file (*.mitk)")); + + MITK_INFO << fileName; + + mitk::SceneIO::Pointer sceneIO = mitk::SceneIO::New(); + mitk::DataStorage::Pointer dataStorage = sceneIO->LoadScene(fileName.toStdString(), 0, true); + + + MITK_INFO << m_Controls.patientImageSelector->currentIndex(); + if (m_Controls.patientImageSelector->currentIndex() <0) + return; + + //get the planePositionManager Service + ctkPluginContext* context = mitk::org_mitk_gui_qt_standardplanetool_Activator::GetContext(); + ctkServiceReference ppmRef = context->getServiceReference(); + mitk::PlanePositionManagerService* service = context->getService(ppmRef); + + int sliceIndex = 0; + mitk::PlanarCircle::Pointer positionMarker = mitk::PlanarCircle::New(); + + mitk::NodePredicateProperty::Pointer isMarker = mitk::NodePredicateProperty::New("isPositionNode", mitk::BoolProperty::New(true)); + mitk::DataStorage::SetOfObjects::ConstPointer markers; + + //load all nodes that are children of the working node in data storage + mitk::DataNode* imageNode(m_Controls.patientImageSelector->GetSelectedNode()); + markers = dataStorage->GetSubset(isMarker); + MITK_INFO << "------------------------------------"; + MITK_INFO << "Process Image: " << imageNode->GetName(); + + if (!markers->empty()) + { + //iterate over all child nodes with NodePredicateProperty + for (mitk::DataStorage::SetOfObjects::const_iterator iter = markers->begin(); iter != markers->end(); ++iter) + { + MITK_INFO << "Marker" << std::endl; + int markerId; + std::string nodeName = (*iter)->GetName(); + //get position nodes in the data manager and add them with the planepositionmanger service + if (nodeName == "Planeposition coronal") + { + mitk::DataNode::Pointer coronalPosNode = (*iter); + coronalPosNode->GetIntProperty("sliceIndex", sliceIndex); + positionMarker = dynamic_cast (coronalPosNode->GetData()); + mitk::PlaneGeometry::ConstPointer corPlane = positionMarker->GetPlaneGeometry(); + markerId = service->AddNewPlanePosition(corPlane, sliceIndex); + coronalPosNode->SetIntProperty("Plane Position ID", markerId); + this->GetDataStorage()->Add(coronalPosNode, imageNode); + } + else if (nodeName == "Planeposition axial") + { + mitk::DataNode::Pointer axialPosNode = (*iter); + axialPosNode->GetIntProperty("sliceIndex", sliceIndex); + positionMarker = dynamic_cast (axialPosNode->GetData()); + mitk::PlaneGeometry::ConstPointer axiPlane = positionMarker->GetPlaneGeometry(); + markerId = service->AddNewPlanePosition(axiPlane, sliceIndex); + axialPosNode->SetIntProperty("Plane Position ID", markerId); + this->GetDataStorage()->Add(axialPosNode, imageNode); + } + else if (nodeName == "Planeposition sagittal") + { + mitk::DataNode::Pointer sagittalPosNode = (*iter); + sagittalPosNode->GetIntProperty("sliceIndex", sliceIndex); + positionMarker = dynamic_cast (sagittalPosNode->GetData()); + mitk::PlaneGeometry::ConstPointer sagPlane = positionMarker->GetPlaneGeometry(); + markerId = service->AddNewPlanePosition(sagPlane, sliceIndex); + sagittalPosNode->SetIntProperty("Plane Position ID", markerId); + this->GetDataStorage()->Add(sagittalPosNode, imageNode); + } + + (*iter)->GetIntProperty("Plane Position ID", markerId); + this->LoadPositionNode((*iter), markerId); + mitk::RenderingManager::GetInstance()->RequestUpdateAll(); + } + } +} + +void QmitkStandardPlaneTool::SaveIntersectionFromNodeToFile() +{ + //get BaseRenderer for the three RenderWindows (sagittal, coronal, axial) + mitk::BaseRenderer::Pointer rendererAxi = mitk::BaseRenderer::GetInstance(mitk::BaseRenderer::GetRenderWindowByName("stdmulti.widget1")); + mitk::BaseRenderer::Pointer rendererSag = mitk::BaseRenderer::GetInstance(mitk::BaseRenderer::GetRenderWindowByName("stdmulti.widget2")); + mitk::BaseRenderer::Pointer rendererCor = mitk::BaseRenderer::GetInstance(mitk::BaseRenderer::GetRenderWindowByName("stdmulti.widget3")); + mitk::PlaneGeometry::ConstPointer sagPlane = rendererSag->GetSliceNavigationController()->GetCurrentPlaneGeometry(); + mitk::PlaneGeometry::ConstPointer corPlane = rendererCor->GetSliceNavigationController()->GetCurrentPlaneGeometry(); + mitk::PlaneGeometry::ConstPointer axiPlane = rendererAxi->GetSliceNavigationController()->GetCurrentPlaneGeometry(); + + mitk::DataNode* imageNode(m_Controls.patientImageSelector->GetSelectedNode()); + std::string saveToMpsString = "C:/Users/thomass/Desktop/" + imageNode->GetName(); + + mitk::Line3D intersectionLine; + sagPlane->IntersectionLine(corPlane, intersectionLine); + + mitk::Point3D newPoint, tempPoint1, tempPoint2; + tempPoint1 = intersectionLine.GetPoint1(); + tempPoint2 = intersectionLine.GetPoint2(); + + newPoint[0] = (tempPoint1[0] - tempPoint2[0]) * 300; + newPoint[1] = (tempPoint1[1] - tempPoint2[1]) * 300; + newPoint[2] = (tempPoint1[2] - tempPoint2[2]) * 300; + + tempPoint1[0] = tempPoint1[0] + newPoint[0]; + tempPoint1[1] = tempPoint1[1] + newPoint[1]; + tempPoint1[2] = tempPoint1[2] + newPoint[2]; + + tempPoint2[0] = tempPoint2[0] - newPoint[0]; + tempPoint2[1] = tempPoint2[1] - newPoint[1]; + tempPoint2[2] = tempPoint2[2] - newPoint[2]; + + mitk::PointSet::Pointer intersectionSagCorPointSet = mitk::PointSet::New(); + intersectionSagCorPointSet->InsertPoint(0, tempPoint1); + intersectionSagCorPointSet->InsertPoint(1, tempPoint2); + + mitk::DataNode::Pointer intersectionSagCorPointNode = mitk::DataNode::New(); + intersectionSagCorPointNode->SetData(intersectionSagCorPointSet); + intersectionSagCorPointNode->SetName("intersection sag/cor"); + // activate property for rendering lines between the points + intersectionSagCorPointNode->SetProperty("show contour", mitk::BoolProperty::New(true)); + + mitk::IOUtil::SavePointSet(intersectionSagCorPointSet, (saveToMpsString + "_line1.mps")); + + + + sagPlane->IntersectionLine(axiPlane, intersectionLine); + + tempPoint1 = intersectionLine.GetPoint1(); + tempPoint2 = intersectionLine.GetPoint2(); + + newPoint[0] = (tempPoint1[0] - tempPoint2[0]) * 300; + newPoint[1] = (tempPoint1[1] - tempPoint2[1]) * 300; + newPoint[2] = (tempPoint1[2] - tempPoint2[2]) * 300; + + tempPoint1[0] = tempPoint1[0] + newPoint[0]; + tempPoint1[1] = tempPoint1[1] + newPoint[1]; + tempPoint1[2] = tempPoint1[2] + newPoint[2]; + + tempPoint2[0] = tempPoint2[0] - newPoint[0]; + tempPoint2[1] = tempPoint2[1] - newPoint[1]; + tempPoint2[2] = tempPoint2[2] - newPoint[2]; + + mitk::PointSet::Pointer intersectionSagAxiPointSet = mitk::PointSet::New(); + intersectionSagAxiPointSet->InsertPoint(0, tempPoint1); + intersectionSagAxiPointSet->InsertPoint(1, tempPoint2); + + mitk::DataNode::Pointer intersectionSagAxiPointNode = mitk::DataNode::New(); + intersectionSagAxiPointNode->SetData(intersectionSagAxiPointSet); + intersectionSagAxiPointNode->SetName("intersection sag/axi"); + // activate property for rendering lines between the points + intersectionSagAxiPointNode->SetProperty("show contour", mitk::BoolProperty::New(true)); + + mitk::IOUtil::SavePointSet(intersectionSagAxiPointSet, (saveToMpsString + "_line2.mps")); + + + + corPlane->IntersectionLine(axiPlane, intersectionLine); + + tempPoint1 = intersectionLine.GetPoint1(); + tempPoint2 = intersectionLine.GetPoint2(); + //tempPoint1 = (tempPoint1-tempPoint2)*10; + newPoint[0] = (tempPoint1[0] - tempPoint2[0]) * 300; + newPoint[1] = (tempPoint1[1] - tempPoint2[1]) * 300; + newPoint[2] = (tempPoint1[2] - tempPoint2[2]) * 300; + + tempPoint1[0] = tempPoint1[0] + newPoint[0]; + tempPoint1[1] = tempPoint1[1] + newPoint[1]; + tempPoint1[2] = tempPoint1[2] + newPoint[2]; + + tempPoint2[0] = tempPoint2[0] - newPoint[0]; + tempPoint2[1] = tempPoint2[1] - newPoint[1]; + tempPoint2[2] = tempPoint2[2] - newPoint[2]; + + mitk::PointSet::Pointer intersectionCorAxiPointSet = mitk::PointSet::New(); + intersectionCorAxiPointSet->InsertPoint(0, tempPoint1); + intersectionCorAxiPointSet->InsertPoint(1, tempPoint2); + + mitk::DataNode::Pointer intersectionCorAxiPointNode = mitk::DataNode::New(); + intersectionCorAxiPointNode->SetData(intersectionCorAxiPointSet); + intersectionCorAxiPointNode->SetName("intersection cor/axi"); + // activate property for rendering lines between the points + intersectionCorAxiPointNode->SetProperty("show contour", mitk::BoolProperty::New(true)); + + mitk::IOUtil::SavePointSet(intersectionCorAxiPointSet, (saveToMpsString + "_line3.mps")); +} + + +void QmitkStandardPlaneTool::ShowSilhouette(int projNum, vtkSmartPointer projMatrix, vtkSmartPointer surface, QString fileName, int showSilhouette, double offset, int sil, int shape) +{ + + MITK_INFO << surface->GetNumberOfPoints(); + auto silhouette = vtkSmartPointer::New(); + auto camera = vtkSmartPointer::New(); + + auto camTransform = vtkSmartPointer::New(); + camTransform->SetMatrix(projMatrix); + camTransform->Update(); + + silhouette->SetProjectionMatrix(camTransform->GetMatrix()); + silhouette->SetInputData(surface); + silhouette->SetCamera(camera); + silhouette->SetFeatureAngle(90); + silhouette->Update(); + + auto transform2 = vtkSmartPointer::New(); + transform2->SetMatrix(projMatrix); + transform2->Update(); + + auto transform3 = vtkSmartPointer::New(); + transform3->SetMatrix(projMatrix); + transform3->Update(); + + //auto transformFilter = vtkSmartPointer::New(); + + //if (showSilhouette ==0) + //transformFilter->SetInputData(surface); //silhouette->GetOutput()); //87 + //else + // transformFilter->SetInputData(silhouette->GetOutput()); //87 + + //transformFilter->SetTransform(transform2); + //transformFilter->Update(); + + //MITK_INFO << transformFilter->GetOutput()->GetNumberOfPoints(); + + // auto polydata = vtkSmartPointer::New(); + auto points = vtkSmartPointer::New(); + //for (int i = 0; i < transformFilter->GetOutput()->GetNumberOfPoints(); i++) + //{ + // double* p; + // transformFilter->GetOutput()->GetPoint(i, p); + // p[0] = p[0] / p[2]; + // p[1] = p[1] / p[2]; + // p[2] = p[2] / p[2]; + // points->InsertPoint(i, p); + //} + + //auto polydata = surface;//transformFilter->GetOutput();//silhouette->GetOutput(); //87 + // vtkDataArray* normals = polydata->GetPointData()->GetNormals(); + int idx = 0; + + for (int id = 0; id < silhouette->GetOutput()->GetPoints()->GetNumberOfPoints(); id++) + { + double normal[3]; + double point[4]; + double* point4D; + // double* trafoNormal; + // double* trafoPoint; + //normals->GetTuple(id, normal); + silhouette->GetOutput()->GetPoint(id, point); + point[0] = point[0]; + point[1] = point[1]; + point[2] = point[2]; + point[3] = 1.0; + + point4D = projMatrix->MultiplyDoublePoint(point); + + point[0] = (point4D[0] / point4D[2]) + offset; + point[1] = point4D[1] / point4D[2]; + point[2] = 0;// point4D[2] / point4D[2]; + // MITK_INFO << "p" << " " << point[0] << " " << point[1] << " " << point[2]; + points->InsertPoint(id, point); + + //point[0] = (point4D[0] / point4D[2]) + offset; + //point[1] = point4D[1] / point4D[2]; + //point[2] = 1; + + //points->InsertPoint(idx++, point); + //point[0] = (point4D[0] / point4D[2]) + offset; + //point[1] = point4D[1] / point4D[2]; + //point[2] = -1; + //points->InsertPoint(idx++, point); + + + // trafoNormal = transform2->TransformDoublePoint(point[0] + normal[0], point[1] + normal[1], point[2] + normal[2]); + // trafoPoint = transform3->TransformDoublePoint(polydata->GetPoint(id)); + // normals->SetTuple3(id, trafoPoint[0] - trafoNormal[0], trafoPoint[1] - trafoNormal[1], trafoPoint[2] - trafoNormal[2]); + } + silhouette->GetOutput()->SetPoints(points); + silhouette->GetOutput()->Modified(); + auto cleanPolyDataFilter = vtkSmartPointer::New(); + cleanPolyDataFilter->SetInputData(silhouette->GetOutput()); + cleanPolyDataFilter->Update(); + + + MITK_INFO << "Silhouette: " << projNum; + typedef itk::Image< unsigned char, 2> ImageType; + ImageType::Pointer imageT = ImageType::New(); + ImageType::IndexType start; + start.Fill(0); + ImageType::SizeType size; + size.Fill(1024); + + ImageType::RegionType region(start, size); + imageT->SetRegions(region); + imageT->Allocate(); + imageT->FillBuffer(0); + + mitk::Image::Pointer mitkimageT; + mitk::CastToMitkImage(imageT, mitkimageT); + + // prepare the binary image's voxel grid + vtkSmartPointer whiteImage = vtkSmartPointer::New(); + whiteImage->DeepCopy(mitkimageT->GetVtkImageData()); + + //// fill the image with foreground voxels: + unsigned char outval = 255; + + + // sweep polygonal data (this is the important thing with contours!) + auto extruder = vtkSmartPointer::New(); + extruder->SetInputData(cleanPolyDataFilter->GetOutput()); + extruder->SetScaleFactor(2); + extruder->SetExtrusionTypeToVectorExtrusion(); + extruder->SetVector(0, 0, 1); + extruder->Update(); + + vtkSmartPointer pol2stenc = vtkSmartPointer::New(); + pol2stenc->SetTolerance(1); //necessary because of the extruder + pol2stenc->SetInputData(extruder->GetOutput()); + pol2stenc->SetOutputOrigin(whiteImage->GetOrigin()); + pol2stenc->SetOutputSpacing(whiteImage->GetSpacing()); + pol2stenc->SetOutputWholeExtent(whiteImage->GetExtent()); + pol2stenc->Update(); + + // cut the corresponding white image and set the background: + vtkSmartPointer imgstenc = vtkSmartPointer::New(); + + imgstenc->SetInputData(whiteImage); + imgstenc->SetStencilConnection(pol2stenc->GetOutputPort()); + imgstenc->ReverseStencilOn(); + imgstenc->SetBackgroundValue(outval); + imgstenc->Update(); + + mitkimageT->SetVolume(imgstenc->GetOutput()->GetScalarPointer());//resultImage->GetScalarPointer() + + mitk::CastToItkImage(mitkimageT,imageT); + + typedef itk::BinaryNotImageFilter BinaryNotImageFilterType; + BinaryNotImageFilterType::Pointer binaryNotFilter = BinaryNotImageFilterType::New(); + binaryNotFilter->SetInput(imageT); + binaryNotFilter->Update(); + + typedef itk::ConnectedComponentImageFilter ConnectedComponentImageFilterType; + ConnectedComponentImageFilterType::Pointer connectedFilter = ConnectedComponentImageFilterType::New(); + connectedFilter->SetInput(binaryNotFilter->GetOutput()); + connectedFilter->Update(); + + typedef itk::LabelShapeKeepNObjectsImageFilter LabelShapeKeepNObjectsImageFilterType; + LabelShapeKeepNObjectsImageFilterType::Pointer labelshapeFilter = LabelShapeKeepNObjectsImageFilterType::New(); + labelshapeFilter->SetInput(connectedFilter->GetOutput()); + labelshapeFilter->SetBackgroundValue(0); + labelshapeFilter->SetNumberOfObjects(1); + labelshapeFilter->SetAttribute(LabelShapeKeepNObjectsImageFilterType::LabelObjectType::NUMBER_OF_PIXELS); + labelshapeFilter->Update(); + + typedef itk::RescaleIntensityImageFilter RescaleType; + RescaleType::Pointer rescaleFilter = RescaleType::New(); + rescaleFilter->SetOutputMinimum(0); + rescaleFilter->SetOutputMaximum(itk::NumericTraits::max()); + rescaleFilter->SetInput(labelshapeFilter->GetOutput()); + rescaleFilter->Update(); + + BinaryNotImageFilterType::Pointer binaryNotFilter2 = BinaryNotImageFilterType::New(); + binaryNotFilter2->SetInput(rescaleFilter->GetOutput()); + binaryNotFilter2->Update(); + mitk::CastToMitkImage(binaryNotFilter2->GetOutput(), mitkimageT);//m_Segmentation + + MITK_INFO << "closing done"; + mitkimageT->Modified(); + + + //transformFilter->GetOutput()->GetPointData()->SetNormals(normals); + auto projectedTibia = mitk::Surface::New(); + projectedTibia->SetVtkPolyData(extruder->GetOutput());// transformFilter->GetOutput());// transformFilter->GetOutput()); + //projectedTibia->GetVtkPolyData()->GetCellData()->SetNormals(normalGenerator->GetOutput()); + projectedTibia->Modified(); + + + if (sil == 1) + { + m_Silhouette1->SetName(std::to_string(projNum)); + m_Silhouette1->SetData(projectedTibia); + m_Silhouette1->SetVisibility(true); + m_Silhouette1->Modified(); + } + if (sil ==2) + { + m_Silhouette2->SetName(std::to_string(projNum)); + m_Silhouette2->SetData(projectedTibia); + m_Silhouette2->SetVisibility(true); + m_Silhouette2->Modified(); + } + + + if (shape == 2) + { + auto convert2DTo3DImageFilter = mitk::Convert2Dto3DImageFilter::New(); + convert2DTo3DImageFilter->SetInput(mitkimageT); + convert2DTo3DImageFilter->Update(); + + //MITK_INFO << "Segmentation: " << projNum; + //auto segmentation2 = mitk::DataNode::New(); + //segmentation2->SetProperty("name", mitk::StringProperty::New("ImageSource")); + //segmentation2->SetData(mitkimageT); + //segmentation2->SetColor(1.0, 1.0, 1.0); + //segmentation2->SetVisibility(true); + //segmentation2->SetName("segmentation"); + //this->GetDataStorage()->Add(segmentation2); + + + + if (fileName == nullptr) + return; + mitk::IOUtil::Save(convert2DTo3DImageFilter->GetOutput(), fileName.toStdString() + ".nrrd"); + MITK_INFO << fileName << " Saved"; + + } + mitk::RenderingManager::GetInstance()->RequestUpdateAll(); + + +} + +void QmitkStandardPlaneTool::OnManualRegistrationClicked() +{ + //config file + boost::property_tree::ptree pt; + boost::property_tree::ini_parser::read_ini("C:/Users/thomass.AD/Desktop/config.ini", pt); + std::cout << "Projection1: " << pt.get("Section1.Projection1") << std::endl; + std::cout << "Projection1: " << pt.get("Section1.Projection2") << std::endl; + std::cout << "SurfaceFile: " << pt.get("Section1.SurfaceFile") << std::endl; + std::cout << "ProjectionFile: " << pt.get("Section1.ProjectionFile") << std::endl; + std::cout << "CSVFile: " << pt.get("Section1.CSVFile") << std::endl; + + int showSilhouette = pt.get("Section1.Silhouette"); + m_ProjNum1 = pt.get("Section1.Projection1"); + m_ProjNum2 = pt.get("Section1.Projection2"); + + // load projection matrix + QString csvFileName = QString::fromUtf8(pt.get("Section1.CSVFile").c_str()); + std::ifstream data(csvFileName.toStdString()); + + // load surface + QString surfaceFileName = QString::fromUtf8(pt.get("Section1.SurfaceFile").c_str()); + m_Surface = mitk::IOUtil::LoadSurface(surfaceFileName.toStdString()); + + QString imageFileName = QString::fromUtf8(pt.get("Section1.ProjectionFile").c_str()); + mitk::Image::Pointer image = mitk::IOUtil::LoadImageA(imageFileName.toStdString()); + + if (csvFileName.isEmpty()) + return; + + std::string line; + for (int i = 0; i < 100; i++) + { + m_ProjectionMat.push_back(vtkSmartPointer::New()); + m_ProjectionMat[i]->Identity(); + } + MITK_INFO << "projection matrix"; + int numProj = 0; + for (int i = 0; i < 100; i++) + { + + for (int j = 0; j < 3; j++) + { + std::getline(data, line); + + std::stringstream lineStream(line); + std::string cell; + int num = 0; + + while (std::getline(lineStream, cell, ',')) + { + m_ProjectionMat[numProj]->SetElement(j, num++, std::stod(cell)); + } + num = 0; + } + numProj++; + } + + auto selector = mitk::ImageSliceSelector::New(); + selector->SetSliceNr(m_ProjNum1); + selector->SetInput(image); + selector->Update(); + mitk::Image::Pointer projection1 = selector->GetOutput(); + + mitk::Point3D origin = projection1->GetGeometry()->GetOrigin(); + origin[2] = 0.0; + projection1->GetGeometry()->SetOrigin(origin[0]); + projection1->Modified(); + + auto convert2DTo3DImageFilter = mitk::Convert2Dto3DImageFilter::New(); + convert2DTo3DImageFilter->SetInput(projection1); + convert2DTo3DImageFilter->Update(); + + m_Projection1 = mitk::DataNode::New(); + m_Projection1->SetProperty("name", mitk::StringProperty::New("ImageSource")); + m_Projection1->SetData(projection1); + m_Projection1->SetColor(1.0, 1.0, 1.0); + m_Projection1->SetVisibility(true); + m_Projection1->SetName("Projection1"); + this->GetDataStorage()->Add(m_Projection1); + + + auto selector2 = mitk::ImageSliceSelector::New(); + selector2->SetSliceNr(m_ProjNum2); + selector2->SetInput(image); + selector2->Update(); + mitk::Image::Pointer projection2 = selector2->GetOutput(); + + mitk::Point3D origin2 = projection2->GetGeometry()->GetOrigin(); + origin2[2] = 0.0; + origin2[0] = origin2[0] + 1024; + projection2->GetGeometry()->SetOrigin(origin2); + projection2->Modified(); + + auto convert2DTo3DImageFilter2 = mitk::Convert2Dto3DImageFilter::New(); + convert2DTo3DImageFilter2->SetInput(projection2); + convert2DTo3DImageFilter2->Update(); + + m_Projection2 = mitk::DataNode::New(); + m_Projection2->SetProperty("name", mitk::StringProperty::New("ImageSource")); + m_Projection2->SetData(projection2); + m_Projection2->SetColor(1.0, 1.0, 1.0); + m_Projection2->SetVisibility(true); + m_Projection2->SetName("Projection2"); + this->GetDataStorage()->Add(m_Projection2); + + m_Silhouette1 = mitk::DataNode::New(); + m_Silhouette2 = mitk::DataNode::New(); + ShowSilhouette(m_ProjNum1, m_ProjectionMat[m_ProjNum1], m_Surface->GetVtkPolyData(),nullptr, true, 0.0,1,2); + ShowSilhouette(m_ProjNum2, m_ProjectionMat[m_ProjNum2], m_Surface->GetVtkPolyData(),nullptr,true, 1024.0,2,2); + + this->GetDataStorage()->Add(m_Silhouette1); + this->GetDataStorage()->Add(m_Silhouette2); + + mitk::RenderingManager::GetInstance()->RequestUpdateAll(); + +} + + + +cv::Mat readEdgeMap(const QString csvFileName, mitk::DataStorage* datastorage) +{ + int mySizes[2] = { 1024, 1024 }; // row, column different from others + cv::Mat edgeMap = cv::Mat(mySizes[0],mySizes[1], CV_64F); + + std::ifstream file(csvFileName.toStdString()); + + for (int row = 0; row < 1024; row++) + { + std::string line; + std::getline(file, line); + if (!file.good()) + break; + + std::stringstream iss(line); + + for (int col = 0; col < 1024; col++) + { + std::string val; + std::getline(iss, val, ','); + if (!iss.good()) + break; + + std::stringstream convertor(val); + float value; + convertor >> value; + edgeMap.at(col, row) = value; + } + } + + typedef itk::Image< unsigned char, 2 > ImageType; + ImageType::Pointer imageT = ImageType::New(); + ImageType::IndexType start; + start.Fill(0); + ImageType::SizeType size; + size.Fill(1024); + + ImageType::RegionType region(start, size); + imageT->SetRegions(region); + imageT->Allocate(); + imageT->FillBuffer(0); + + + for (int n = 0; n < 1024; n++) + { + for (int m = 0; m < 1024; m++) + { + if (edgeMap.at(m, n)>0) + { + ImageType::IndexType pIndex; + pIndex[0] = m; + pIndex[1] = n; + ImageType::PixelType value = 1000.0; + + imageT->SetPixel(pIndex, value); + ImageType::PixelType pixelValue = imageT->GetPixel(pIndex); + } + } + } + mitk::Image::Pointer mitkimageT; + mitk::CastToMitkImage(imageT, mitkimageT); + + + auto datanode2 = mitk::DataNode::New(); + datanode2->SetName("ContourImage"); + // datanode->SetBoolProperty("show contour", true); + //datanode2->SetProperty("contourcolor", ColorProperty::New(0.0, 1.0, 0.0)); + datanode2->SetData(mitkimageT); + datastorage->Add(datanode2); + + return edgeMap; +} + + + +void QmitkStandardPlaneTool::OnSliderValueChanged() +{ + MITK_INFO << "Change Value"; + typedef itk::Euler3DTransform< double > TransformType; + TransformType::Pointer itkTransform = TransformType::New(); + TransformType::ParametersType rigid_params(6); + rigid_params[0] = m_Controls.sB_Tx->value(); + rigid_params[1] = m_Controls.sB_Ty->value(); + rigid_params[2] = m_Controls.sB_Tz->value(); + rigid_params[3] = m_Controls.sB_alpha->value(); + rigid_params[4] = m_Controls.sB_beta->value(); + rigid_params[5] = m_Controls.sB_gamma->value(); + + itkTransform->SetParameters(rigid_params); + + auto vtkMatrix = vtkSmartPointer::New(); + mitk::TransferItkTransformToVtkMatrix(itkTransform.GetPointer(), vtkMatrix); + auto transformFilter = vtkSmartPointer::New(); + auto vtktransform = vtkSmartPointer::New(); + vtktransform->SetMatrix(vtkMatrix); + transformFilter->SetInputData(m_Surface->GetVtkPolyData()); + transformFilter->SetTransform(vtktransform); + transformFilter->Update(); + + this->ShowSilhouette(m_ProjNum1, m_ProjectionMat[m_ProjNum1], transformFilter->GetOutput(),nullptr, 1, 0.0, 1,1); + this->ShowSilhouette(m_ProjNum2, m_ProjectionMat[m_ProjNum2], transformFilter->GetOutput(),nullptr, 1, 1024.0, 2,1); + +} + +void ConvertToSurface(mitk::Image::Pointer segmentation, mitk::Surface::Pointer& surface) +{ + // Create a contour from the binary image + auto surfaceFilter = mitk::ManualSegmentationToSurfaceFilter::New(); + surfaceFilter->SetInput(segmentation); + surfaceFilter->SetThreshold(0.5); //expects binary image with zeros and ones + surfaceFilter->SetGaussianStandardDeviation(1.5); + surfaceFilter->SetUseGaussianImageSmooth(true); // apply gaussian to thresholded image ? + surfaceFilter->SetSmooth(true); + surfaceFilter->InterpolationOn(); + surfaceFilter->SetMedianFilter3D(true); // apply median to segmentation before marching cubes ? + surfaceFilter->SetDecimate(mitk::ImageToSurfaceFilter::NoDecimation); + surfaceFilter->UpdateLargestPossibleRegion(); + + ////----------------- REMESHING ---------------------------------- + //auto remesher = mitk::ACVD::RemeshFilter::New(); + //remesher->SetInput(surfaceFilter->GetOutput()); + //remesher->SetTimeStep(0); + //remesher->SetNumVertices(1000); + //remesher->SetGradation(1.0); + //// remesher->SetBoundaryFixing(true) + //remesher->SetSubsampling(50); + //remesher->SetEdgeSplitting(0.0); + //remesher->SetOptimizationLevel(1); + //remesher->SetForceManifold(false); + //remesher->SetBoundaryFixing(false); + //remesher->Update(); + surface = surfaceFilter->GetOutput()->Clone(); +} + +void QmitkStandardPlaneTool::OnExtractContourAndSlicesClicked() +{ + //config file + boost::property_tree::ptree pt; + boost::property_tree::ini_parser::read_ini("C:/Users/thomass.AD/Desktop/configMask.ini", pt); + std::cout << "Folder name: " << pt.get("Section1.Folder") << std::endl; + + mitk::NodePredicateDataType::Pointer isImage = mitk::NodePredicateDataType::New("Image"); + mitk::DataStorage::SetOfObjects::ConstPointer images = this->GetDataStorage()->GetSubset(isImage); + if (!images->empty()) + { + //iterate over all child nodes with NodePredicateProperty + for (mitk::DataStorage::SetOfObjects::const_iterator iter = images->begin(); iter != images->end(); ++iter) + { + + + mitk::DataNode::Pointer node = (*iter); + mitk::Image::Pointer image = dynamic_cast(node->GetData()); + QString nodeName = QString::fromUtf8((*iter)->GetName().c_str()); + + QString folderName = QString::fromUtf8(pt.get("Section1.Folder").c_str()); + QString fibulaFileName = folderName+(QString("/Reorganized/"))+(nodeName+(QString("_Fibula.nrrd"))); + QString talusFileName = folderName+(QString("/Reorganized/"))+(nodeName+(QString("_Talus.nrrd"))); + QString tibiaFileName = folderName+(QString("/Reorganized/"))+(nodeName+(QString("_Tibia.nrrd"))); + QString csvFileName = folderName+(QString("/Reorganized/"))+(nodeName+(QString("/")))+(nodeName+(QString(".csv"))); + QString projectionFileName = folderName+(QString("/XA/"))+(nodeName+(QString("/")))+(nodeName+(QString(".IMA"))); + QString outputFileName = folderName + (QString("/XA/")) + (nodeName + (QString("/"))); + std::ifstream data(csvFileName.toStdString()); + + // load segmentations + + mitk::Image::Pointer fibula = mitk::IOUtil::LoadImageA(fibulaFileName.toStdString()); + mitk::Image::Pointer tibia = mitk::IOUtil::LoadImageA(tibiaFileName.toStdString()); + mitk::Image::Pointer talus = mitk::IOUtil::LoadImageA(talusFileName.toStdString()); + mitk::Image::Pointer projection = mitk::IOUtil::LoadImageA(projectionFileName.toStdString()); + + auto surfaceTibia = mitk::Surface::New(); + auto surfaceFibula = mitk::Surface::New(); + auto surfaceTalus = mitk::Surface::New(); + + ConvertToSurface(fibula, surfaceFibula); + ConvertToSurface(tibia, surfaceTibia); + ConvertToSurface(talus, surfaceTalus); + + if (csvFileName.isEmpty()) + return; + + m_ProjectionMat.empty(); + std::string line; + for (int i = 0; i < 100; i++) + { + m_ProjectionMat.push_back(vtkSmartPointer::New()); + m_ProjectionMat[i]->Identity(); + } + MITK_INFO << "projection matrix"; + int numProj = 0; + for (int i = 0; i < 100; i++) + { + + for (int j = 0; j < 3; j++) + { + std::getline(data, line); + + std::stringstream lineStream(line); + std::string cell; + int num = 0; + + while (std::getline(lineStream, cell, ',')) + { + m_ProjectionMat[numProj]->SetElement(j, num++, std::stod(cell)); + } + num = 0; + } + numProj++; + } + + auto selector = mitk::ImageSliceSelector::New(); + selector->SetSliceNr(m_ProjNum1); + selector->SetInput(image); + selector->Update(); + mitk::Image::Pointer projection1 = selector->GetOutput(); + + mitk::Point3D origin = projection1->GetGeometry()->GetOrigin(); + origin[2] = 0.0; + projection1->GetGeometry()->SetOrigin(origin[0]); + projection1->Modified(); + + auto convert2DTo3DImageFilter = mitk::Convert2Dto3DImageFilter::New(); + convert2DTo3DImageFilter->SetInput(projection1); + convert2DTo3DImageFilter->Update(); + + + m_Silhouette1 = mitk::DataNode::New(); + m_Silhouette2 = mitk::DataNode::New(); + + //const char* path = (outputFileName + nodeName + QString("/")).toStdString().c_str(); + + //if (_mkdir(path)) + //{ + // return; + //} + //else + //{ + // std::cerr << "Directory Created: " << path << std::endl; + + //} + + for (int i = 0; i < 100; i++) + { + ShowSilhouette(i, m_ProjectionMat[i], surfaceFibula->GetVtkPolyData(), outputFileName+(nodeName+((("_Fibula" + std::to_string(i)).c_str()))), true, 0.0, 1, 2); + ShowSilhouette(i, m_ProjectionMat[i], surfaceTibia->GetVtkPolyData(), outputFileName+(nodeName+((("_Tibia" + std::to_string(i)).c_str()))), true, 0.0, 1, 2); + ShowSilhouette(i, m_ProjectionMat[i], surfaceTalus->GetVtkPolyData(), outputFileName+(nodeName+((("_Talus" + std::to_string(i)).c_str()))), true, 0.0, 1, 2); + ExtractSlice(this->GetDataStorage(), projection, i, outputFileName+(nodeName+((('_' + std::to_string(i)).c_str()))), true, true); + } + + + /* auto segmentation2 = mitk::DataNode::New(); + segmentation2->SetProperty("name", mitk::StringProperty::New("ImageSource")); + segmentation2->SetData(mask); + segmentation2->SetColor(1.0, 1.0, 1.0); + segmentation2->SetVisibility(true); + segmentation2->SetName("segmentation1"); + this->GetDataStorage()->Add(segmentation2);*/ + } + } + + +} + + + +void QmitkStandardPlaneTool::OnRegistrationClicked() +{ + + //config file + boost::property_tree::ptree pt; + boost::property_tree::ini_parser::read_ini("C:/Users/thomass.AD/Desktop/config.ini", pt); + std::cout << "Projection1: " << pt.get("Section1.Projection1") << std::endl; + std::cout << "Projection1: " << pt.get("Section1.Projection2") << std::endl; + std::cout << "SurfaceFile: " << pt.get("Section1.SurfaceFile") << std::endl; + std::cout << "EdgeMap1: " << pt.get("Section1.EdgeMapFile1") << std::endl; + std::cout << "EdgeMap2: " << pt.get("Section1.EdgeMapFile2") << std::endl; + std::cout << "ProjectionFile: " << pt.get("Section1.ProjectionFile") << std::endl; + std::cout << "CSVFile: " << pt.get("Section1.CSVFile") << std::endl; + + int showSilhouette = pt.get("Section1.Silhouette"); + int projNum1 = pt.get("Section1.Projection1"); + int projNum2 = pt.get("Section1.Projection2"); + //return; + // load edgemap 1 and 2 + QString csvFileNameEdgeMap1 = QString::fromUtf8(pt.get("Section1.EdgeMapFile1").c_str()); //QFileDialog::getOpenFileName(nullptr, "Select a .csv file", "", "CSV (*.csv)"); + cv::Mat edgeMap1 = readEdgeMap(csvFileNameEdgeMap1, this->GetDataStorage()); + QString csvFileNameEdgeMap2 = QString::fromUtf8(pt.get("Section1.EdgeMapFile2").c_str()); //QFileDialog::getOpenFileName(nullptr, "Select a .csv file", "", "CSV (*.csv)"); + cv::Mat edgeMap2 = readEdgeMap(csvFileNameEdgeMap2, this->GetDataStorage()); + + // load projection matrix + QString csvFileName = QString::fromUtf8(pt.get("Section1.CSVFile").c_str());//QFileDialog::getOpenFileName(nullptr, "Select a .csv file", "", "CSV (*.csv)"); + std::ifstream data(csvFileName.toStdString()); + + // load surface + QString surfaceFileName = QString::fromUtf8(pt.get("Section1.SurfaceFile").c_str());//QFileDialog::getOpenFileName(nullptr, "Select a surface file", "", "Surface (*.stl)"); + m_Surface = mitk::IOUtil::LoadSurface(surfaceFileName.toStdString()); + + QString imageFileName = QString::fromUtf8(pt.get("Section1.ProjectionFile").c_str()); + mitk::Image::Pointer image = mitk::IOUtil::LoadImageA(imageFileName.toStdString()); + + if (csvFileName.isEmpty()) + return; + + std::string line; + std::vector> projectionMat; + for (int i = 0; i < 100; i++) + { + projectionMat.push_back(vtkSmartPointer::New()); + projectionMat[i]->Identity(); + } + MITK_INFO << "projection matrix"; + int numProj = 0; + for (int i = 0; i < 100; i++) + { + + for (int j = 0; j < 3; j++) + { + std::getline(data, line); + + std::stringstream lineStream(line); + std::string cell; + int num = 0; + + while (std::getline(lineStream, cell, ',')) + { + projectionMat[numProj]->SetElement(j, num++, std::stod(cell)); + } + num = 0; + } + numProj++; + } + + typedef itk::Euler3DTransform< double > TransformType; + TransformType::Pointer itkTransform = TransformType::New(); + TransformType::ParametersType rigid_params(6); + rigid_params[0] = 0; + rigid_params[1] = 0; + rigid_params[2] = 0; + rigid_params[3] = 0; + rigid_params[4] = 0; + rigid_params[5] = 0; + + itkTransform->SetParameters(rigid_params); + + auto vtkMatrix = vtkSmartPointer::New(); + mitk::TransferItkTransformToVtkMatrix(itkTransform.GetPointer(), vtkMatrix); + auto transformFilter = vtkSmartPointer::New(); + auto vtktransform = vtkSmartPointer::New(); + vtktransform->SetMatrix(vtkMatrix); + transformFilter->SetInputData(m_Surface->GetVtkPolyData()); + transformFilter->SetTransform(vtktransform); + transformFilter->Update(); + m_Surface->SetVtkPolyData(transformFilter->GetOutput()); + m_Surface->Modified(); + + MITK_INFO << "set up registration"; + mitk2D3DShapeRegistration registration;// = mitk2D3DShapeRegistration::New(); + registration.SetSurface(m_Surface); + MITK_INFO << "set up edgemaps"; + registration.SetEdgeMap1(edgeMap1); + registration.SetEdgeMap2(edgeMap2); + MITK_INFO << "set up projection matrix"; + registration.SetProjectionMatrix1(projectionMat[projNum1]); + registration.SetProjectionMatrix2(projectionMat[projNum2]); + MITK_INFO << "generate"; + registration.Generate(); + + m_Surface->SetVtkPolyData(registration.GetResultSurface()); + m_Surface->Modified(); + m_Silhouette1 = mitk::DataNode::New(); + m_Silhouette2 = mitk::DataNode::New(); + + ShowSilhouette(projNum1, projectionMat[projNum1], m_Surface->GetVtkPolyData(), nullptr, true, 0.0, 1,1); + ShowSilhouette(projNum2, projectionMat[projNum2], m_Surface->GetVtkPolyData(), nullptr, true, 0.0, 2,1); + + this->GetDataStorage()->Add(m_Silhouette1); + this->GetDataStorage()->Add(m_Silhouette2); + + +} + +void QmitkStandardPlaneTool::OnExtractProjectionsClicked() +{ + + mitk::NodePredicateDataType::Pointer isImage = mitk::NodePredicateDataType::New("Image"); + mitk::DataStorage::SetOfObjects::ConstPointer images = this->GetDataStorage()->GetSubset(isImage); + if (!images->empty()) + { + //iterate over all child nodes with NodePredicateProperty + for (mitk::DataStorage::SetOfObjects::const_iterator iter = images->begin(); iter != images->end(); ++iter) + { + + + mitk::DataNode::Pointer node= (*iter); + mitk::Image::Pointer image = dynamic_cast(node->GetData()); + + for (int i = 0; i < 100; i++) + { + QString nodeName = QString::fromUtf8((*iter)->GetName().c_str()); + ExtractSlice(this->GetDataStorage(), image, i, QString("C:/Users/thomass.AD/Desktop/Data/SilhouetteExp/")+(nodeName+((('_'+std::to_string(i)).c_str()))), true, true); + } + + } + } + +} + +void QmitkStandardPlaneTool::OnExtractSilhouetteClicked() +{ + QString surfaceFileName = QFileDialog::getOpenFileName(nullptr, "Select a surface file", "", "Surface (*.stl)"); + QString imageFileName = QFileDialog::getOpenFileName(nullptr, "Select an image file", "", "Image (*.nrrd)"); + + if (surfaceFileName.isEmpty()) + return; + + if (imageFileName.isEmpty()) + return; + + mitk::Image::Pointer image = mitk::IOUtil::LoadImageA(imageFileName.toStdString()); + // ******************************************************* + // load projection matrices + // ******************************************************* + QString csvFileName = QFileDialog::getOpenFileName(nullptr, "Select a .csv file", "", "CSV (*.csv)"); + std::ifstream data(csvFileName.toStdString()); + + if (csvFileName.isEmpty()) + return; + + + std::string line; + std::vector> projectionMat; + for (int i = 0; i < 100; i++) + { + projectionMat.push_back(vtkSmartPointer::New()); + projectionMat[i]->Identity(); + } + MITK_INFO << "projection matrix"; + int numProj = 0; + for (int i = 0; i < 100; i++) + { + + for (int j = 0; j < 3; j++) + { + std::getline(data, line); + + std::stringstream lineStream(line); + std::string cell; + int num = 0; + + while (std::getline(lineStream, cell, ',')) + { + projectionMat[numProj]->SetElement(j, num++, std::stod(cell)); + // MITK_INFO << projectionMat[numProj]->GetElement(j, num - 1); + } + num = 0; + } + numProj++; + } + MITK_INFO << "Decompose Matrix"; + // Create a synthetic projection matrix + cv::Mat P(3, 4, cv::DataType::type); + P.at(0, 0) = projectionMat[0]->GetElement(0, 0); + P.at(1, 0) = projectionMat[0]->GetElement(1, 0); + P.at(2, 0) = projectionMat[0]->GetElement(2, 0); + + P.at(0, 1) = projectionMat[0]->GetElement(0, 1); + P.at(1, 1) = projectionMat[0]->GetElement(1, 1); + P.at(2, 1) = projectionMat[0]->GetElement(2, 1); + + P.at(0, 2) = projectionMat[0]->GetElement(0, 2); + P.at(1, 2) = projectionMat[0]->GetElement(1, 2); + P.at(2, 2) = projectionMat[0]->GetElement(2, 2); + + P.at(0, 3) = projectionMat[0]->GetElement(0, 3); + P.at(1, 3) = projectionMat[0]->GetElement(1, 3); + P.at(2, 3) = projectionMat[0]->GetElement(2, 3); + + MITK_INFO << P.at(2, 3); + + // Decompose the projection matrix into: + cv::Mat K(3, 3, cv::DataType::type); // intrinsic parameter matrix + cv::Mat R(3, 3, cv::DataType::type); // rotation matrix + cv::Mat T(4, 1, cv::DataType::type); // translation vector + + cv::decomposeProjectionMatrix(P, K, R, T); + + MITK_INFO << "K Matrix"; + MITK_INFO << K.at(0, 0); + MITK_INFO << K.at(1, 0); + MITK_INFO << K.at(2, 0); + MITK_INFO << K.at(0, 1); + MITK_INFO << K.at(1, 1); + MITK_INFO << K.at(2, 1); + MITK_INFO << K.at(0, 2); + MITK_INFO << K.at(1, 2); + MITK_INFO << K.at(2, 2); + + MITK_INFO << "R Matrix"; + MITK_INFO << R.at(0, 0); + MITK_INFO << R.at(1, 0); + MITK_INFO << R.at(2, 0); + MITK_INFO << R.at(0, 1); + MITK_INFO << R.at(1, 1); + MITK_INFO << R.at(2, 1); + MITK_INFO << R.at(0, 2); + MITK_INFO << R.at(1, 2); + MITK_INFO << R.at(2, 2); + + MITK_INFO << "T Matrix"; + MITK_INFO << T.at(0, 0); + MITK_INFO << T.at(1, 0); + MITK_INFO << T.at(2, 0); + MITK_INFO << T.at(3, 0); + + mitk::Surface::Pointer surface = mitk::IOUtil::LoadSurface(surfaceFileName.toStdString()); + mitk::Surface::Pointer surface2 = mitk::IOUtil::LoadSurface(surfaceFileName.toStdString()); + + for (int i = 0; i < 40; i=i+39) + { + ShowSilhouette(i, projectionMat[i], surface2->Clone()->GetVtkPolyData(),nullptr,1,0.0,1,1); + ExtractSlice(this->GetDataStorage(), image, i, QString("C:/Users/thomass.AD/Desktop/Data/SilhouetteExp/test_")+((std::to_string(i).c_str())), false, true); + } +// +// std::cout << "There are " << surface->GetVtkPolyData()->GetNumberOfCells() << " cells." << std::endl; +// std::cout << "There are " << surface->GetVtkPolyData()->GetNumberOfPoints() << " points." << std::endl; +// +// +// +// auto extractEdges = vtkSmartPointer::New(); +// extractEdges->SetInputData(surface->GetVtkPolyData()); +// extractEdges->Update(); +// +// vtkCellArray* lines = extractEdges->GetOutput()->GetLines(); +// vtkPoints* points = extractEdges->GetOutput()->GetPoints(); +// +// std::cout << std::endl << "Edges" << endl << "----------" << std::endl; +// std::cout << "There are " << lines->GetNumberOfCells() << " cells." << std::endl; +// std::cout << "There are " << points->GetNumberOfPoints() << " points." << std::endl; +// +// //// Traverse all of the edges +// //for (vtkIdType i = 0; i < extractEdges->GetOutput()->GetNumberOfCells()/100; i++) +// //{ +// // std::cout << "Type: " << extractEdges->GetOutput()->GetCell(i)->GetClassName() << std::endl; +// // vtkSmartPointer line = vtkLine::SafeDownCast(extractEdges->GetOutput()->GetCell(i)); +// //// std::cout << "Line " << i << " : " << *line << std::endl; +// +// //} +// +// //************************************************************************ +// +// auto triangleFilter = vtkSmartPointer::New(); +// triangleFilter->SetInputData(surface->GetVtkPolyData()); +// triangleFilter->Update(); +// +// vtkIdType cellId = 10; +// +// auto cellPointIds = vtkSmartPointer::New(); +// triangleFilter->GetOutput()->GetCellPoints(cellId, cellPointIds); +// +// +// +// +// +// +// //for (i = 0; iGetVtkPolyData()->GetNumberOfCells(); i++) +// //{ +// // cell = input->GetCell(i); +// // for (j = 0; j<3; j++) +// // { +// // cellEdgeNeighbors->Initialize(); +// // id0 = cell->GetEdge(j)->GetPointIds()->GetId(0); +// // id1 = cell->GetEdge(j)->GetPointIds()->GetId(1); +// // input->GetCellEdgeNeighbors(i, id0, id1, cellEdgeNeighbors); +// // if (cellEdgeNeighbors->GetNumberOfIds() == 0) +// // { +// // boundaryIds->InsertUniqueId(id0); +// // boundaryIds->InsertUniqueId(id1); +// // } +// // } +// //} +// +// +// std::list neighbors; +// +// for (vtkIdType i = 0; i < cellPointIds->GetNumberOfIds(); i++) +// { +// auto idList = vtkSmartPointer::New(); +// +// //add one of the edge points +// idList->InsertNextId(cellPointIds->GetId(i)); +// +// //add the other edge point +// if (i + 1 == cellPointIds->GetNumberOfIds()) +// { +// idList->InsertNextId(cellPointIds->GetId(0)); +// } +// else +// { +// idList->InsertNextId(cellPointIds->GetId(i + 1)); +// } +// +// //get the neighbors of the cell +// auto neighborCellIds = vtkSmartPointer::New(); +// +// triangleFilter->GetOutput()->GetCellNeighbors(cellId, idList, neighborCellIds); +// // triangleFilter->GetOutput()->GetCellEdgeNeighbors(cellId, 1, 0, neighborCellIds); +// +// /* id0 = cell->GetEdge(j)->GetPointIds()->GetId(0); +// id1 = cell->GetEdge(j)->GetPointIds()->GetId(1); +// input->GetCellEdgeNeighbors(i, id0, id1, cellEdgeNeighbors); +//*/ +// +// +// for (vtkIdType j = 0; j < neighborCellIds->GetNumberOfIds(); j++) +// { +// neighbors.push_back(neighborCellIds->GetId(j)); +// } +// +// } +// +// std::cout << "Edge neighbor ids are: " << std::endl; +// +// for (std::list::iterator it1 = neighbors.begin(); it1 != neighbors.end(); it1++) +// { +// std::cout << " " << *it1; +// +// } +// std::cout << std::endl; +// +// +// +// +// // Create a dataset with the cell of interest +// { +// auto ids = vtkSmartPointer::New(); +// ids->SetNumberOfComponents(1); +// ids->InsertNextValue(cellId); +// +// auto selectionNode = vtkSmartPointer::New(); +// selectionNode->SetFieldType(vtkSelectionNode::CELL); +// selectionNode->SetContentType(vtkSelectionNode::INDICES); +// selectionNode->SetSelectionList(ids); +// +// auto selection = vtkSmartPointer::New(); +// selection->AddNode(selectionNode); +// +// auto extractSelection = vtkSmartPointer::New(); +// extractSelection->SetInputData(0, surface->GetVtkPolyData()); +// extractSelection->SetInputData(1, selection); +// extractSelection->Update(); +// +// } +// +// //mainCellActor->GetProperty()->SetColor(1, 0, 0); +// +// // Create a dataset with the neighbor cells +// { +// auto ids = vtkSmartPointer::New(); +// ids->SetNumberOfComponents(1); +// for (std::list::iterator it1 = neighbors.begin(); it1 != neighbors.end(); it1++) +// { +// ids->InsertNextValue(*it1); +// } +// +// vtkSmartPointer selectionNode = +// vtkSmartPointer::New(); +// selectionNode->SetFieldType(vtkSelectionNode::CELL); +// selectionNode->SetContentType(vtkSelectionNode::INDICES); +// selectionNode->SetSelectionList(ids); +// +// vtkSmartPointer selection = +// vtkSmartPointer::New(); +// selection->AddNode(selectionNode); +// +// vtkSmartPointer extractSelection = +// vtkSmartPointer::New(); +// extractSelection->SetInputData(surface->GetVtkPolyData()); +// extractSelection->SetInputData(1, selection); +// extractSelection->Update(); +// +// } +// +// neighborCellsActor->GetProperty()->SetColor(0, 1, 0); + + + //************************************************************************* + +} + + +void QmitkStandardPlaneTool::SetFocus() +{ +} + + diff --git a/Plugins/org.mitk.gui.qt.standardplanetool/src/internal/StandardPlaneToolControls.ui b/Plugins/org.mitk.gui.qt.standardplanetool/src/internal/StandardPlaneToolControls.ui new file mode 100644 index 0000000000..01b11a14e2 --- /dev/null +++ b/Plugins/org.mitk.gui.qt.standardplanetool/src/internal/StandardPlaneToolControls.ui @@ -0,0 +1,76 @@ + + + StandardPlaneToolControls + + + + 0 + 0 + 365 + 863 + + + + + 0 + 0 + + + + QmitkTemplate + + + + + + Manual adjustment mode + + + + + + + LoadPosition + + + + + + + Save and Load next Dataset + + + + + + + SavePosition + + + + + + + Show Slice + + + + + + + Qt::Vertical + + + + 20 + 40 + + + + + + + + + + diff --git a/Plugins/org.mitk.gui.qt.standardplanetool/src/internal/itkShapeContourCostFunctionBckup.hxx b/Plugins/org.mitk.gui.qt.standardplanetool/src/internal/itkShapeContourCostFunctionBckup.hxx new file mode 100644 index 0000000000..b518b9196e --- /dev/null +++ b/Plugins/org.mitk.gui.qt.standardplanetool/src/internal/itkShapeContourCostFunctionBckup.hxx @@ -0,0 +1,121 @@ +/*========================================================================= + * + * Copyright Insight Software Consortium + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#ifndef itkShapeContourCostFunction_hxx +#define itkShapeContourCostFunction_hxx + +#include "itkShapeContourCostFunction.h" + +namespace itk +{ +/** + * Constructor + */ +template< typename TFeatureImage, typename TOutputPixel > +ShapeContourCostFunction< TFeatureImage, TOutputPixel > +::ShapeContourCostFunction() +{ + +} + +/** + * PrintSelf + */ +template< typename TFeatureImage, typename TOutputPixel > +void +ShapeContourCostFunction< TFeatureImage, TOutputPixel > +::PrintSelf(std::ostream & os, Indent indent) const +{ + Superclass::PrintSelf(os, indent); + os << indent << "Contours: " << std::endl; + +} + + + +/** + * + */ +template< typename TFeatureImage, typename TOutputPixel > +typename ShapeContourCostFunction< TFeatureImage, TOutputPixel > +::MeasureType +ShapeContourCostFunction< TFeatureImage, TOutputPixel > +::ComputeDistanceTerm(const ParametersType & parameters) const +{ + +} + + + + +/** + * + */ +template< typename TFeatureImage, typename TOutputPixel > +void +ShapeContourCostFunction< TFeatureImage, TOutputPixel > +::Initialize() throw ( ExceptionObject ) +{ + this->Superclass::Initialize(); + + // check if the mean and variances array are of the right size + if ( m_ShapeParameterMeans.Size() < + this->m_ShapeFunction->GetNumberOfShapeParameters() ) + { + itkExceptionMacro(<< "ShapeParameterMeans does not have at least " + << this->m_ShapeFunction->GetNumberOfShapeParameters() + << " number of elements."); + } + + if ( m_ShapeParameterStandardDeviations.Size() < + this->m_ShapeFunction->GetNumberOfShapeParameters() ) + { + itkExceptionMacro(<< "ShapeParameterStandardDeviations does not have at least " + << this->m_ShapeFunction->GetNumberOfShapeParameters() + << " number of elements."); + } +} + + +/** + * + */ +template< typename TFeatureImage, typename TOutputPixel > +typename ShapeContourCostFunction< TFeatureImage, TOutputPixel > +::MeasureType +ShapeContourCostFunction< TFeatureImage, TOutputPixel > +::GetValue(const ParametersType & parameters) const +{ + return (this->ComputeDistanceTerm(parameters) ); +} + +/** + * + */ +template< typename TFeatureImage, typename TOutputPixel > +void +ShapeContourCostFunction< TFeatureImage, TOutputPixel > +::Initialize(void) +throw ( ExceptionObject ) +{ + +} + + +} // end namespace itk + +#endif \ No newline at end of file diff --git a/Plugins/org.mitk.gui.qt.standardplanetool/src/internal/mitk2D3DShapeRegistration.cpp b/Plugins/org.mitk.gui.qt.standardplanetool/src/internal/mitk2D3DShapeRegistration.cpp new file mode 100644 index 0000000000..0ce0a72410 --- /dev/null +++ b/Plugins/org.mitk.gui.qt.standardplanetool/src/internal/mitk2D3DShapeRegistration.cpp @@ -0,0 +1,264 @@ + +#include "mitk2D3DShapeRegistration.h" +#include +#include +#include +#include +#include +#include "itkEuler3DTransform.h" +#include +#include "mitkShapeContourCostFunction.h" +#include +#include +#include + +#include +#include +#include + +#include + +mitk2D3DShapeRegistration::mitk2D3DShapeRegistration(): +m_surface(mitk::Surface::New()), +m_projectionmatrix1(vtkSmartPointer::New()), +m_projectionmatrix2(vtkSmartPointer::New()) +{ + + int mySizes[2] = { 1024, 1024 }; // row, column different from others + m_edgeMap1 = cv::Mat(2, mySizes, CV_32F); + m_edgeMap2 = cv::Mat(2, mySizes, CV_32F); + m_gradImage1 = cv::Mat(2, mySizes, CV_32F); + m_gradImage2 = cv::Mat(2, mySizes, CV_32F); +} + + +mitk2D3DShapeRegistration::~mitk2D3DShapeRegistration() +{ + return; +} + + +void mitk2D3DShapeRegistration::SetSurface(mitk::Surface::Pointer surface) +{ + m_surface = surface; + m_surface->Modified(); +} +void mitk2D3DShapeRegistration::SetEdgeMap1(cv::Mat edgemap) +{ + m_edgeMap1 = edgemap; +} +void mitk2D3DShapeRegistration::SetEdgeMap2(cv::Mat edgemap) +{ + m_edgeMap2 = edgemap; +} + +void mitk2D3DShapeRegistration::SetImage1(cv::Mat image) +{ + m_image1 = image; +} +void mitk2D3DShapeRegistration::SetImage2(cv::Mat image) +{ + m_image2 = image; +} + + +void mitk2D3DShapeRegistration::CalcGradientDirectionImage(cv::Mat image, cv::Mat &gradImage, int kernelsize) +{ + m_image1 = image; + GaussianBlur(m_image1, m_image1, cv::Size(kernelsize, kernelsize), 0, 0, cv::BORDER_DEFAULT); + cv::Mat Sx; + Sobel(m_image1, Sx, CV_32F, 1, 0, 3); + + cv::Mat Sy; + Sobel(m_image1, Sy, CV_32F, 0, 1, 3); + + cv::Mat magnitudeImage, orientationImage; + cv::magnitude(Sx, Sy, magnitudeImage); + cv::phase(Sx, Sy, orientationImage, true); + // normalize(orientationImage, orientationImage, 0.0, 255.0, cv::NORM_MINMAX, CV_8U); + + + gradImage = orientationImage;//.convertTo(gradImage,CV_32F); + //imshow("rg", gradImage); +} + +void mitk2D3DShapeRegistration::SetProjectionMatrix1(vtkSmartPointer mat) +{ + m_projectionmatrix1 = mat; +} + +void mitk2D3DShapeRegistration::SetProjectionMatrix2(vtkSmartPointer mat) +{ + m_projectionmatrix2 = mat; +} + + +vtkSmartPointer mitk2D3DShapeRegistration::GetResultSurface() +{ + return m_surface->GetVtkPolyData(); +} + +vtkSmartPointer mitk2D3DShapeRegistration::Generate() +{ + + if (m_image1.empty() || m_image2.empty()) + return 0; + + this->CalcGradientDirectionImage(m_image1, m_gradImage1, 3); + this->CalcGradientDirectionImage(m_image2, m_gradImage2, 3); + //return nullptr; + mitk::ShapeContourCostFunction::Pointer costFunction = mitk::ShapeContourCostFunction::New(); + costFunction->SetProjectionMatrix1(m_projectionmatrix1); + costFunction->SetProjectionMatrix2(m_projectionmatrix2); + costFunction->SetEdgeMap1(m_edgeMap1); + costFunction->SetEdgeMap2(m_edgeMap2); + costFunction->SetGradientImage1(m_gradImage1); + costFunction->SetGradientImage2(m_gradImage2); + costFunction->SetSurface(m_surface->GetVtkPolyData()); + + + // The following lines instantiate the evolutionary optimizer. + + typedef itk::PowellOptimizer OptimizerType;// + OptimizerType::Pointer optimizer = OptimizerType::New(); + + // Evolutionary algorithms are based on testing random variations + // of parameters. In order to support the computation of random values, + // ITK provides a family of random number generators. In this example, we + // use the \doxygen{NormalVariateGenerator} which generates values with a + // normal distribution. + + itk::Statistics::NormalVariateGenerator::Pointer generator + = itk::Statistics::NormalVariateGenerator::New(); + + // The random number generator must be initialized with a seed. + + generator->Initialize(1234567); + + // The OnePlusOneEvolutionaryOptimizer is initialized by + // specifying the random number generator, the number of samples for the + // initial population and the maximum number of iterations. + // optimizer->SetNormalVariateGenerator(generator); + optimizer->SetCostFunction(costFunction); + optimizer->SetMaximize(false); + optimizer->SetStepLength(0.5); + optimizer->SetStepTolerance(0.1); + optimizer->SetValueTolerance(0.01); + optimizer->SetMaximumIteration(100); + + + + // optimizer->SetGradientMagnitudeTolerance(0.01); + // optimizer->SetMinimumStepLength(0.001); + // optimizer->SetMaximumStepLength(1); + + // optimizer->Initialize(10); + // optimizer->SetInitialRadius(0.1); + // optimizer->SetMaximumIteration(100); + + // scale the components appropriately + itk::Array parametersScale; + parametersScale.SetSize(6); + parametersScale[0] = 10; // / maxMV; + parametersScale[1] = 10; // / maxMV; + parametersScale[2] = 10; // / maxMV; + parametersScale[3] = 0.1; // / maxMV; + parametersScale[4] = 0.1; // / maxMV; + parametersScale[5] = 0.1; // / maxMV; + + //for (unsigned int i = 4; i<6; i++) + //{ + // parametersScale[i] = 256; // offset scale + //} + optimizer->SetScales(parametersScale); + + // Here we instantiate the Command object that will act as an + // observer of the registration method and print out parameters at each + // iteration. Earlier, we defined this command as a class templated over the + // optimizer type. Once it is created with the \code{New()} method, we + // connect the optimizer to the command. + + typedef IterationCallback< OptimizerType > IterationCallbackType; + IterationCallbackType::Pointer callback = IterationCallbackType::New(); + callback->SetOptimizer(optimizer); + + OptimizerType::ParametersType params; + params.SetSize(6); + params[0] = 0.0; + params[1] = 0.0; + params[2] = 0.0; + params[3] = 0.0; + params[4] = 0.0; + params[5] = 0.0; + optimizer->SetInitialPosition(params); + + std::cout << "Initial Parameters : " << params << std::endl; + + //// make optimizer maximize + //optimizer->MinimizeOn(); + + // do registration + + try + { + optimizer->StartOptimization(); + std::cout << "Optimizer stop condition: " + << optimizer->GetStopConditionDescription() + << std::endl; + } + catch (itk::ExceptionObject & exp) + { + std::cerr << "Exception caught ! " << std::endl; + std::cerr << exp << std::endl; + } + + // get registration result + OptimizerType::ParametersType finalParameters = + optimizer->GetCurrentPosition(); + + std::cout << "Final Solution is : " << finalParameters << std::endl; + + + float tx = static_cast(finalParameters[0]); + float ty = static_cast(finalParameters[1]); + float tz = static_cast(finalParameters[2]); + int rx = static_cast(finalParameters[3]); + int ry = static_cast(finalParameters[4]); + int rz = static_cast(finalParameters[5]); + + //MITK_INFO << tx << " " << ty << " " << tz; + + + typedef itk::Euler3DTransform< double > TransformType; + TransformType::Pointer itkTransform = TransformType::New(); + TransformType::ParametersType rigid_params(6); + rigid_params[0] = tx; + rigid_params[1] = ty; + rigid_params[2] = tz; + rigid_params[3] = rx; + rigid_params[4] = ry; + rigid_params[5] = rz; + + itkTransform->SetParameters(rigid_params); + + auto vtkMatrix = vtkSmartPointer::New(); + mitk::TransferItkTransformToVtkMatrix(itkTransform.GetPointer(), vtkMatrix); + //vtkMatrix->Identity(); + //vtkMatrix->SetElement(0, 3, 20); + //vtkMatrix->SetElement(1, 3, 10); + //vtkMatrix->SetElement(2, 3, 10); + auto transformFilter = vtkSmartPointer::New(); + auto vtktransform = vtkSmartPointer::New(); + vtktransform->SetMatrix(vtkMatrix); + // vtktransform->Inverse(); + transformFilter->SetTransform(vtktransform); + transformFilter->SetInputData(m_surface->GetVtkPolyData()); + transformFilter->Update(); + + m_surface->SetVtkPolyData(transformFilter->GetOutput()); + m_surface->Modified(); + + // MITK_INFO << "Num: "<GetOutput()->GetNumberOfPoints(); + return transformFilter->GetOutput(); //m_surface->GetVtkPolyData(); // + +} \ No newline at end of file diff --git a/Plugins/org.mitk.gui.qt.standardplanetool/src/internal/mitk2D3DShapeRegistration.h b/Plugins/org.mitk.gui.qt.standardplanetool/src/internal/mitk2D3DShapeRegistration.h new file mode 100644 index 0000000000..0a508db32b --- /dev/null +++ b/Plugins/org.mitk.gui.qt.standardplanetool/src/internal/mitk2D3DShapeRegistration.h @@ -0,0 +1,58 @@ +/*=================================================================== + +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. + +===================================================================*/ +#ifndef mitk2D3DShapeRegistration_h +#define mitk2D3DShapeRegistration_h + +#include +#include +#include +#include + +class mitk2D3DShapeRegistration +{ +public: + mitk2D3DShapeRegistration(); + + ~mitk2D3DShapeRegistration(); + +// static Pointer New(void); + + void SetSurface(mitk::Surface::Pointer surface); + void SetEdgeMap1(cv::Mat edgemap); + void SetEdgeMap2(cv::Mat edgemap); + void SetProjectionMatrix1(vtkSmartPointer mat); + void SetProjectionMatrix2(vtkSmartPointer mat); + void SetImage1(cv::Mat image); + void SetImage2(cv::Mat image); + void CalcGradientDirectionImage(cv::Mat image, cv::Mat& gradImage, int kernelsize = 3); + + + vtkSmartPointer Generate(); + vtkSmartPointer GetResultSurface(); + +private: + mitk::Surface::Pointer m_surface; + vtkSmartPointer m_projectionmatrix1; + vtkSmartPointer m_projectionmatrix2; + cv::Mat m_edgeMap1; + cv::Mat m_edgeMap2; + cv::Mat m_gradImage1; + cv::Mat m_gradImage2; + cv::Mat m_image1; + cv::Mat m_image2; +}; + +#endif //mitk2D3DShapeRegistration \ No newline at end of file diff --git a/Plugins/org.mitk.gui.qt.standardplanetool/src/internal/mitkShapeContourCostFunction.cpp b/Plugins/org.mitk.gui.qt.standardplanetool/src/internal/mitkShapeContourCostFunction.cpp new file mode 100644 index 0000000000..6ea45c2034 --- /dev/null +++ b/Plugins/org.mitk.gui.qt.standardplanetool/src/internal/mitkShapeContourCostFunction.cpp @@ -0,0 +1,449 @@ +/*=================================================================== + +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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +#include "mitkShapeContourCostFunction.h" +#include "vtkSilhouette.h" +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "itkEuler3DTransform.h" +#include "itkTimeProbe.h" + +namespace mitk +{ + + // set projection matrix + class ShapeContourCostFunction::Impl + { + public: + Impl(): + m_projectionMatrix(vtkSmartPointer::New()), + m_silhouette(vtkSmartPointer::New()), + m_kDTreeProjection( vtkSmartPointer::New()), + m_edgeMap(vtkSmartPointer::New()), + m_camTransform(vtkSmartPointer::New()), + m_pointTransform(vtkSmartPointer::New()), + m_normalTransform(vtkSmartPointer::New()) + { + int mySizes[2] = { 1024, 1024 }; // row, column different from others + m_gradImage = cv::Mat(mySizes[0], mySizes[1], CV_32F); + + }; + + void SetUpKDTree() + { + m_kDTreeProjection->SetDataSet(m_edgeMap); + m_kDTreeProjection->BuildLocator(); + + auto idxToWorldTransformInv = vtkSmartPointer::New(); + idxToWorldTransformInv->Identity(); + + m_camTransform->SetMatrix(m_projectionMatrix); + m_camTransform->Concatenate(idxToWorldTransformInv); + m_camTransform->Update(); + + m_pointTransform->SetMatrix(m_projectionMatrix); + m_pointTransform->Concatenate(idxToWorldTransformInv); + // m_pointTransform->PostMultiply(); + m_pointTransform->Update(); + + m_normalTransform->SetMatrix(m_projectionMatrix); + m_normalTransform->Concatenate(idxToWorldTransformInv); + // m_normalTransform->PostMultiply(); + m_normalTransform->Update(); + + } + + itk::SingleValuedCostFunction::MeasureType GetValue(vtkSmartPointer surface) + { + double distance = 0.0; + auto camera = vtkSmartPointer::New(); + auto silhouette = vtkSmartPointer::New(); + auto poly_silhouette = vtkSmartPointer::New(); + // MITK_INFO << "s " << surface->GetNumberOfPoints(); + silhouette->SetProjectionMatrix(m_camTransform->GetMatrix()); + silhouette->SetInputData(surface); + silhouette->SetCamera(camera); + silhouette->InitializeMesh(); + silhouette->InitializeProjectionGeometry(); + silhouette->SetFeatureAngle(90); + silhouette->Update(); + + auto points = vtkSmartPointer::New(); + + vtkSmartPointer polydata = silhouette->GetOutput(); + vtkDataArray* normals = polydata->GetPointData()->GetNormals(); + for (int id = 0; id < normals->GetNumberOfTuples(); id++) + { + double normal[3]; + double point[4]; + double* point4D; + double* trafoNormal; + double* trafoPoint; + normals->GetTuple(id, normal); + polydata->GetPoint(id, point); + + point[0] = point[0]; + point[1] = point[1]; + point[2] = point[2]; + point[3] = 1.0; + + point4D = m_projectionMatrix->MultiplyDoublePoint(point); + + point[0] = point4D[0] / point4D[2]; + point[1] = point4D[1] / point4D[2]; + point[2] = point4D[2] / point4D[2]; + // MITK_INFO << "p" << " " << point[0] << " " << point[1] << " " << point[2]; + points->InsertPoint(id, point); + + trafoNormal = m_pointTransform->TransformDoublePoint(point[0] + normal[0], point[1] + normal[1], point[2] + normal[2]); + trafoPoint = m_normalTransform->TransformDoublePoint(polydata->GetPoint(id)); + normals->SetTuple3(id, trafoPoint[0] - trafoNormal[0], trafoPoint[1] - trafoNormal[1], trafoPoint[2] - trafoNormal[2]); + } + // MITK_INFO << "b "<GetNumberOfPoints(); + poly_silhouette->SetPoints(points); + poly_silhouette->GetPointData()->SetNormals(normals); + //MITK_INFO << "a " << m_silhouette->GetNumberOfPoints(); + vtkSmartPointer pointNormals = vtkDataArray::SafeDownCast(poly_silhouette->GetPointData()->GetNormals()); + + int numValidPoints = 0; + // kd tree locator + for (int i = 0; i < poly_silhouette->GetNumberOfPoints(); i++) + { + + // point to closest point distance + vtkIdType iD = m_kDTreeProjection->FindClosestPoint(poly_silhouette->GetPoint(i)); + double closestPoint[3]; + m_kDTreeProjection->GetDataSet()->GetPoint(iD, closestPoint); + + + // normal eval + // TODO optimize cv calc? + double pN[3]; + pointNormals->GetTuple(i, pN); + float angle; + angle = cv::fastAtan2((float) pN[0],(float) pN[1]); + + + + double anglediff =( (double)angle)- m_gradImage.at((int)closestPoint[1], (int)closestPoint[0]); // row and column order different from vtk + + //MITK_INFO << m_gradImage.at((int)closestPoint[1], (int)closestPoint[0]) << " " << angle << " " << anglediff; + if (abs(anglediff) < 30) + { + distance = distance + sqrt(vtkMath::Distance2BetweenPoints(poly_silhouette->GetPoint(i), closestPoint)); + + numValidPoints++; + } + + } + MITK_INFO <<"Points "<< numValidPoints; + if (numValidPoints>2) + return distance /= (numValidPoints-1); + else + return 100000; //TODO + }; + + + //Input + vtkSmartPointer m_projectionMatrix; + vtkSmartPointer m_silhouette; + vtkSmartPointer m_kDTreeProjection; + vtkSmartPointer m_edgeMap; + vtkSmartPointer m_camTransform; + vtkSmartPointer m_pointTransform; + vtkSmartPointer m_normalTransform; + cv::Mat m_gradImage; + + bool m_Active; + }; + + + ShapeContourCostFunction::ShapeContourCostFunction() + //(const CvImage& drrImage, const CvImage& carmImage, const cv::Mat& _Mask = cv::Mat()) + : m_Impl1(new Impl), + m_Impl2(new Impl) + { + //set projection matrix + //set surface + + } + + ShapeContourCostFunction::~ShapeContourCostFunction() + { + delete m_Impl1; + delete m_Impl2; + } + + void ShapeContourCostFunction::Initialize() + { + // tree from edge map + // check all edges of the projection surface + // distance + + + } + + void ShapeContourCostFunction::SetEdgeMap2(cv::Mat edgeMap) + { + auto points = vtkSmartPointer::New(); + + for (int row = 0; row < 1024; ++row) + for (int col = 0; col < 1024; ++col) + if (edgeMap.at(col, row)>0) + points->InsertNextPoint(col, row, 1); + + m_Impl2->m_edgeMap->SetPoints(points); + +// auto polydata = vtkSmartPointer::New(); +// polydata->SetPoints(points); +// +// vtkSmartPointer vertexFilter = +// vtkSmartPointer::New(); +//#if VTK_MAJOR_VERSION <= 5 +// vertexFilter->SetInputConnection(polydata->GetProducerPort()); +//#else +// vertexFilter->SetInputData(polydata); +//#endif +// vertexFilter->Update(); +// +// vtkSmartPointer polydata2 = +// vtkSmartPointer::New(); +// polydata2->ShallowCopy(vertexFilter->GetOutput()); +// +// +// // Visualization +// vtkSmartPointer mapper = +// vtkSmartPointer::New(); +//#if VTK_MAJOR_VERSION <= 5 +// mapper->SetInputConnection(polydata2->GetProducerPort()); +//#else +// mapper->SetInputData(polydata2); +//#endif +// +// vtkSmartPointer actor = +// vtkSmartPointer::New(); +// actor->SetMapper(mapper); +// actor->GetProperty()->SetPointSize(5); +// +// +// +// vtkSmartPointer renderer = +// vtkSmartPointer::New(); +// vtkSmartPointer renderWindow = +// vtkSmartPointer::New(); +// renderWindow->AddRenderer(renderer); +// vtkSmartPointer renderWindowInteractor = +// vtkSmartPointer::New(); +// renderWindowInteractor->SetRenderWindow(renderWindow); +// +// renderer->AddActor(actor); +// +// renderWindow->Render(); +// renderWindowInteractor->Start(); + + + m_Impl2->SetUpKDTree(); + + + + } + + void ShapeContourCostFunction::SetEdgeMap1(cv::Mat edgeMap) + { + auto points = vtkSmartPointer::New(); + + for (int row = 0; row < 1024; ++row) + for (int col = 0; col < 1024; ++col) + if (edgeMap.at(col,row)>0) + points->InsertNextPoint(col,row, 1); + + m_Impl1->m_edgeMap->SetPoints(points); + m_Impl1->SetUpKDTree(); + } + + void ShapeContourCostFunction::SetGradientImage1(cv::Mat gradientImage) + { + if (gradientImage.empty()) + return; + this->m_Impl1->m_gradImage = gradientImage; + } + + void ShapeContourCostFunction::SetGradientImage2(cv::Mat gradientImage) + { + if (gradientImage.empty()) + return; + this->m_Impl2->m_gradImage = gradientImage; + } + + + void ShapeContourCostFunction::SetProjectionMatrix1(vtkSmartPointer projMatrix) + { m_Impl1->m_projectionMatrix = projMatrix; } + + void ShapeContourCostFunction::SetProjectionMatrix2(vtkSmartPointer projMatrix) + { m_Impl2->m_projectionMatrix = projMatrix; } + + + itk::SingleValuedCostFunction::MeasureType ShapeContourCostFunction::GetValue(const itk::SingleValuedCostFunction::ParametersType &x) const + { + double tx = static_cast(x[0]); + double ty = static_cast(x[1]); + double tz = static_cast(x[2]); + double rx = static_cast(x[3]); + double ry = static_cast(x[4]); + double rz = static_cast(x[5]); + + MITK_INFO << tx << " " << ty << " " << tz <<" "<< rx << " " << ry << " " << rz; + + + typedef itk::Euler3DTransform< double > TransformType; + TransformType::Pointer itkTransform = TransformType::New(); + TransformType::ParametersType rigid_params(6); + rigid_params[0] = tx; + rigid_params[1] = ty; + rigid_params[2] = tz; + rigid_params[3] = rx; + rigid_params[4] = ry; + rigid_params[5] = rz; + + itkTransform->SetParameters(rigid_params); + + auto vtkMatrix = vtkSmartPointer::New(); + mitk::TransferItkTransformToVtkMatrix(itkTransform.GetPointer(), vtkMatrix); + auto transformFilter = vtkSmartPointer::New(); + auto vtktransform = vtkSmartPointer::New(); + vtktransform->SetMatrix(vtkMatrix); + transformFilter->SetTransform(vtktransform); + auto surface = vtkSmartPointer::New(); + surface->DeepCopy(m_surface); + transformFilter->SetInputData(surface); + transformFilter->Update(); + + //double *bound = transformFilter->GetOutput()->GetBounds(); + //MITK_INFO << "b: " << bound[0] << " " << bound[1] << " " << bound[2] << " " << bound[3]; + // itk::TimeProbe timer; + //timer.Start(); + + itk::SingleValuedCostFunction::MeasureType val1 = m_Impl1->GetValue(transformFilter->GetOutput()); + itk::SingleValuedCostFunction::MeasureType val2 = m_Impl2->GetValue(transformFilter->GetOutput()); + + MITK_INFO << "Metric value is: " << val1 <<" "< ImageType; +//ImageType::Pointer imageT = ImageType::New(); +//ImageType::IndexType start; +//start.Fill(0); +//ImageType::SizeType size; +//size.Fill(1024); +//ImageType::RegionType region(start, size); +//imageT->SetRegions(region); +//imageT->Allocate(); +//imageT->FillBuffer(1000); + +//for (int n = 0; n < transformFilter->GetOutput()->GetNumberOfPoints(); n++) +//{ +// double* point; +// point = projectedTibia->GetVtkPolyData()->GetPoint(n); +// ImageType::IndexType pIndex; +// pIndex[0] = ceil(point[0]); +// pIndex[1] = ceil(point[1]); +// ImageType::PixelType value = 1.0; + +// imageT->SetPixel(pIndex, value); +// ImageType::PixelType pixelValue = imageT->GetPixel(pIndex); +//} +//mitk::Image::Pointer mitkimageT; +//mitk::CastToMitkImage(imageT, mitkimageT); + +// +//endodebug("converting to gray image") +//int channels = _2dImage.channels(); +//endodebugvar(channels) +//cv::Mat _2dImageGray = _2dImage; +//if (channels > 1) +//{ +// cv::cvtColor(_2dImage, _2dImageGray, CV_BGR2GRAY); +//} +// +//cv::Mat _8bit2dImageGray = _2dImageGray; +//int type2dImage = _2dImageGray.type(); +//endodebugvar(type2dImage) +//if (type2dImage == CV_16UC1) +//{ +// endodebug("converting to 8 bit") +// _2dImageGray.convertTo(_8bit2dImageGray, CV_8U, 255. / 65535.); +//} +///* +//else if( type2dImage == CV_32S ) +//{ +//endodebug("converting to 8 bit from 32 bit") +//_2dImage.convertTo(eightBit2dImage, CV_8U, 255./2147483647.); +//} +//*/ +// +//assert(_8bit2dImageGray.type() == CV_8U); +//endodebugimg(_8bit2dImageGray.clone()) +// +//cv::equalizeHist(_8bit2dImageGray, _8bit2dImageGray); +//cv::threshold(_8bit2dImageGray, d->thresholdedImage, d->m_Threshold, 255, cv::THRESH_BINARY_INV); +//endodebugvar(d->thresholdedImage.type()) +//endodebugimg(d->thresholdedImage.clone()) +// +//assert(d->m_Mask.type() == CV_8UC1); +// +//d->m_Active = true; +//d->m_CvImageFromRenderWindow = new QmitkCvImageFromRenderWindow(&d->m_Active, _RenderWindow, &d->renderedImage); + + +// namespace mitk diff --git a/Plugins/org.mitk.gui.qt.standardplanetool/src/internal/mitkShapeContourCostFunction.h b/Plugins/org.mitk.gui.qt.standardplanetool/src/internal/mitkShapeContourCostFunction.h new file mode 100644 index 0000000000..0a39bfb007 --- /dev/null +++ b/Plugins/org.mitk.gui.qt.standardplanetool/src/internal/mitkShapeContourCostFunction.h @@ -0,0 +1,203 @@ +/*=================================================================== + +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. + +===================================================================*/ +#ifndef mitkShapeContourCostFunction_h +#define mitkShapeContourCostFunction_h + +//#include +//#include +//#include +#include +#include +#include +#include +#include "vtkSilhouette.h" + +#include +#include +#include + +#include +#include +#include +#include +#include "itkOnePlusOneEvolutionaryOptimizer.h" +#include "itkNormalVariateGenerator.h" +#include "itkCastImageFilter.h" +#include "itkRescaleIntensityImageFilter.h" +#include "itkImageToSpatialObjectMetric.h" +#include "itkCommand.h" + +#include "cv.h" + +/** +// As in previous registration examples, it is important to +// track the evolution of the optimizer as it progresses through the parameter +// space. This is done by using the Command/Observer paradigm. The +// following lines of code implement the \doxygen{Command} observer that +// monitors the progress of the registration. The code is quite +// similar to what we have used in previous registration examples. +// +// \index{Model to Image Registration!Observer} +*/ + +template < class TOptimizer > +class IterationCallback : public itk::Command +{ +public: + typedef IterationCallback Self; + typedef itk::Command Superclass; + typedef itk::SmartPointer Pointer; + typedef itk::SmartPointer ConstPointer; + + itkTypeMacro(IterationCallback, Superclass); + itkNewMacro(Self); + + /** Type defining the optimizer. */ + typedef TOptimizer OptimizerType; + + /** Method to specify the optimizer. */ + void SetOptimizer(OptimizerType * optimizer) + { + m_Optimizer = optimizer; + m_Optimizer->AddObserver(itk::IterationEvent(), this); + } + + /** Execute method will print data at each iteration */ + void Execute(itk::Object *caller, const itk::EventObject & event) + { + Execute((const itk::Object *)caller, event); + } + + void Execute(const itk::Object *, const itk::EventObject & event) + { + if (typeid(event) == typeid(itk::StartEvent)) + { + std::cout << std::endl << "Position Value"; + std::cout << std::endl << std::endl; + } + else if (typeid(event) == typeid(itk::IterationEvent)) + { + return; + std::cout << m_Optimizer->GetCurrentIteration() << " "; + std::cout << m_Optimizer->GetValue() << " "; + std::cout << m_Optimizer->GetCurrentPosition() << std::endl; + } + else if (typeid(event) == typeid(itk::EndEvent)) + { + std::cout << std::endl << std::endl; + std::cout << "After " << m_Optimizer->GetCurrentIteration(); + std::cout << " iterations " << std::endl; + std::cout << "Solution is = " << m_Optimizer->GetCurrentPosition(); + std::cout << std::endl; + } + } + +protected: + IterationCallback() {}; + itk::WeakPointer m_Optimizer; +}; + +// cost function +/* +- perturbate shape and pose parameter +- project shape and normals +- calculate metric: +- nearest neighbor graph +1 - for each drr contour point --> distance nearest neighbor +2 - for each drr contour point --> difference gradient to surface normal + + + +Input +Params (Pose 1 Shape 2 Appearance 3) +Surface (1) +ProjectionMatrix (1,2,3) +ShapeParams (2,3) +ProjectionImage (Edge / Gradientmap) (1,2,3) + +Output +Value for optimizer (1,2,3) + +*/ + + +namespace mitk +{ + + class ShapeContourCostFunction : + public itk::SingleValuedCostFunction + { + public: + /// + /// init default values and save references + /// + ShapeContourCostFunction();// const CvImage& drrImage, const CvImage& carmImage, const cv::Mat& _Mask = cv::Mat()); + + mitkClassMacroItkParent(ShapeContourCostFunction, itk::SingleValuedCostFunction); + itkNewMacro(Self); + + // get shape metric term + // virtual itk::SingleValuedCostFunction::MeasureType ComputeShapeDistanceTerm(const ParametersType & parameters) const ITK_OVERRIDE; + // get gradient metric term + // virtual itk::SingleValuedCostFunction::MeasureType ComputeGradientDirectionTerm(const ParametersType & parameters) const ITK_OVERRIDE; + // get projection + // virtual itk::SingleValuedCostFunction::MeasureType ComputeCurrentShape(const ParametersType & parameters) const ITK_OVERRIDE; + + /// + /// executes the cost function + /// + itk::SingleValuedCostFunction::MeasureType GetValue(const itk::SingleValuedCostFunction::ParametersType& x) const; + + void SetSurface(vtkSmartPointer surface){ m_surface = surface; }; + + void SetProjectionMatrix1(vtkSmartPointer projMatrix); + void SetProjectionMatrix2(vtkSmartPointer projMatrix); + + void SetEdgeMap2(cv::Mat edgeMap); + void SetEdgeMap1(cv::Mat edgeMap); + + void SetGradientImage2(cv::Mat gradientImage); + void SetGradientImage1(cv::Mat gradientImage); + + void Initialize(); + + + virtual void GetDerivative( + const itk::SingleValuedCostFunction::ParametersType& parameters, + itk::SingleValuedCostFunction::DerivativeType& derivative) const + { + return; + //Superclass::GetDerivative(parameters, derivative); + // throw itk::ExceptionObject(__FILE__, __LINE__, "No derivate available"); + } + + unsigned int GetNumberOfParameters() const; + /// + /// delete d pointer + /// + virtual ~ShapeContourCostFunction(); + + private: + + class Impl; + Impl *m_Impl1; + Impl *m_Impl2; + + vtkSmartPointer m_surface; + }; +} // namespace mitk + +#endif // mitkShapeContourCostFunction_h diff --git a/Plugins/org.mitk.gui.qt.standardplanetool/src/internal/org_mitk_gui_qt_standardplanetool_Activator.cpp b/Plugins/org.mitk.gui.qt.standardplanetool/src/internal/org_mitk_gui_qt_standardplanetool_Activator.cpp new file mode 100644 index 0000000000..6a2b51a4e3 --- /dev/null +++ b/Plugins/org.mitk.gui.qt.standardplanetool/src/internal/org_mitk_gui_qt_standardplanetool_Activator.cpp @@ -0,0 +1,51 @@ +/*========================================================================= + +Program: Medical Imaging & Interaction Toolkit +Language: C++ +Date: $Date$ +Version: $Revision$ + +Copyright (c) German Cancer Research Center, Division of Medical and +Biological Informatics. All rights reserved. +See MITKCopyright.txt or http://www.mitk.org/copyright.html for details. + +This software is distributed WITHOUT ANY WARRANTY; without even +the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR +PURPOSE. See the above copyright notices for more information. + +=========================================================================*/ + +#include "org_mitk_gui_qt_standardplanetool_Activator.h" + +#include + +#include "QmitkStandardPlaneTool.h" +#include +#include + +US_INITIALIZE_MODULE + +namespace mitk { + ctkPluginContext* org_mitk_gui_qt_standardplanetool_Activator::m_Context = 0; + + void org_mitk_gui_qt_standardplanetool_Activator::start(ctkPluginContext* context) + { +// RegisterBoundingShapeObjectFactory(); + BERRY_REGISTER_EXTENSION_CLASS(QmitkStandardPlaneTool, context) + m_Context = context; + } + + void org_mitk_gui_qt_standardplanetool_Activator::stop(ctkPluginContext* context) + { + Q_UNUSED(context) + } + + ctkPluginContext* org_mitk_gui_qt_standardplanetool_Activator::GetContext() + { + return m_Context; + } +} + +#if QT_VERSION < QT_VERSION_CHECK(5, 0, 0) +Q_EXPORT_PLUGIN2(org_mitk_gui_qt_standardplanetool, mitk::PluginActivator) +#endif \ No newline at end of file diff --git a/Plugins/org.mitk.gui.qt.standardplanetool/src/internal/org_mitk_gui_qt_standardplanetool_Activator.h b/Plugins/org.mitk.gui.qt.standardplanetool/src/internal/org_mitk_gui_qt_standardplanetool_Activator.h new file mode 100644 index 0000000000..bb6bda0032 --- /dev/null +++ b/Plugins/org.mitk.gui.qt.standardplanetool/src/internal/org_mitk_gui_qt_standardplanetool_Activator.h @@ -0,0 +1,48 @@ +/*=================================================================== + +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. + +===================================================================*/ +#ifndef org_mitk_gui_qt_standardplanetool_Activator_h +#define org_mitk_gui_qt_standardplanetool_Activator_h + +#include + +namespace mitk { + + class org_mitk_gui_qt_standardplanetool_Activator : + public QObject, public ctkPluginActivator +{ + Q_OBJECT + +#if QT_VERSION >= QT_VERSION_CHECK(5, 0, 0) + Q_PLUGIN_METADATA(IID "org_mitk_gui_qt_standardplanetool") +#endif + Q_INTERFACES(ctkPluginActivator) + +public: + + void start(ctkPluginContext* context) override; + void stop(ctkPluginContext* context) override; + + static ctkPluginContext* GetContext(); + +private: + + static ctkPluginContext* m_Context; + +}; // PluginActivator + +} + +#endif // org_mitk_gui_qt_standardplanetool_Activator_h diff --git a/Plugins/org.mitk.gui.qt.standardplanetool/src/internal/vtkSilhouette.cpp b/Plugins/org.mitk.gui.qt.standardplanetool/src/internal/vtkSilhouette.cpp new file mode 100644 index 0000000000..31f0f74402 --- /dev/null +++ b/Plugins/org.mitk.gui.qt.standardplanetool/src/internal/vtkSilhouette.cpp @@ -0,0 +1,487 @@ +/*========================================================================= + + Program: Visualization Toolkit + Module: vtkSilhouette.cxx + + Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen + All rights reserved. + See Copyright.txt or http://www.kitware.com/Copyright.htm for details. + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notice for more information. + +=========================================================================*/ +// .SECTION Thanks +// Contribution by Thierry Carrard
+// CEA/DIF - Commissariat a l'Energie Atomique, Centre DAM Ile-De-France
+// BP12, F-91297 Arpajon, France.
+ +#include "vtkSilhouette.h" + +#include "vtkCamera.h" +#include "vtkCellData.h" +#include "vtkGenericCell.h" +#include "vtkMath.h" +#include "vtkInformation.h" +#include "vtkInformationVector.h" +#include "vtkObjectFactory.h" +#include "vtkPointData.h" +#include "vtkPolyData.h" +#include "vtkProp3D.h" +#include "vtkTransform.h" +#include "vtkUnsignedIntArray.h" +#include "vtkIdTypeArray.h" +#include "vtkPoints.h" +#include "vtkCellArray.h" +#include "vtkPolygon.h" +#include "vtkSmartPointer.h" +#include "vtkCleanPolyData.h" +#include "vtkPolyDataNormals.h" +#include "vtkFloatArray.h" +#include + +vtkStandardNewMacro(vtkSilhouette); + +vtkCxxSetObjectMacro(vtkSilhouette,Camera,vtkCamera); + +struct vtkOrderedEdge +{ + inline vtkOrderedEdge(vtkIdType a, vtkIdType b) + { + if(a<=b) { p1=a; p2=b; } + else { p1=b; p2=a; } + } + inline bool operator < (const vtkOrderedEdge& oe) const + { + return (p1 edges; + bool * edgeFlag; + vtkCellArray* lines; + inline vtkPolyDataEdges() : edgeFlag(0), lines(0) { vec[0]=vec[1]=vec[2]=0.0; } +}; + +vtkSilhouette::vtkSilhouette() +{ + this->Camera = NULL; + this->Prop3D = NULL; + this->Direction = VTK_DIRECTION_CAMERA_ORIGIN; + this->Vector[0] = this->Vector[1] = this->Vector[2] = 0.0; + this->Origin[0] = this->Origin[1] = this->Origin[2] = 0.0; + this->Transform = vtkTransform::New(); + this->ProjectionMatrix = vtkMatrix4x4::New(); + this->EnableFeatureAngle = 1; + this->FeatureAngle = 60; + this->BorderEdges = 0; + this->PieceInvariant = 1; + this->PreComp = new vtkPolyDataEdges(); +} + +vtkSilhouette::~vtkSilhouette() +{ + //if (this->Transform) + //{ + // MITK_INFO << "delete"; + // this->Transform->Delete(); + //} +// this->ProjectionMatrix->Delete(); + + if ( this->Camera ) + { + this->Camera->Delete(); + } + + delete [] this->PreComp->edgeFlag; + if (this->PreComp->lines) + { + this->PreComp->lines->Delete(); + } + delete this->PreComp; + + //Note: vtkProp3D is not deleted to avoid reference count cycle +} + + +void vtkSilhouette::InitializeProjectionGeometry() +{ + double CMat[3][3]; + double CMatInv[3][3]; + for (int m = 0; m < 3; m++) + for (int n = 0; n < 3; n++) + CMat[m][n] = this->ProjectionMatrix->GetElement(m, n); + + double imageCenter[3]; + imageCenter[0] = this->ProjectionMatrix->GetElement(0, 3); + imageCenter[1] = this->ProjectionMatrix->GetElement(1, 3); + imageCenter[2] = this->ProjectionMatrix->GetElement(2, 3); + + double source[3]; + vtkMath::Invert3x3(CMat, CMatInv); + vtkMath::Multiply3x3(CMatInv, imageCenter, source); + + double isoCenter[3]; + isoCenter[0] = 0.0; + isoCenter[1] = 0.0; + isoCenter[2] = 0.0; + + this->ProjectionMatrix->Invert(); + auto inverseTransform = vtkSmartPointer::New(); + inverseTransform->SetMatrix(this->ProjectionMatrix); + + + + bool vectorMode = true; + double normVec = 1; + double* origin = inverseTransform->TransformDoublePoint(isoCenter); + Origin[0] = origin[0]; + Origin[1] = origin[1]; + Origin[2] = origin[2]; + + // Use specified direction vector from projection matrix origin + CameraVector[0] = (source[0] - Origin[0]); + CameraVector[1] = (source[1] - Origin[1]); + CameraVector[2] = (source[2] - Origin[2]); + + normVec = vtkMath::Norm(CameraVector); + CameraVector[0] = CameraVector[0] / normVec; + CameraVector[1] = CameraVector[1] / normVec; + CameraVector[2] = CameraVector[2] / normVec; + +} + + +void vtkSilhouette::InitializeMesh() +{ + // get the input and ouptut + vtkSmartPointer input = vtkPolyData::SafeDownCast(this->GetPolyDataInput(0)); + auto helperPolyData = vtkSmartPointer::New(); + + if (input == nullptr) + { + vtkErrorMacro(<< "Need correct connections"); + } + vtkDebugMacro(<< "RequestData\n"); + + vtkIdType nPolys = input->GetNumberOfPolys(); + vtkSmartPointer polysArray = input->GetPolys()->GetData(); + vtkSmartPointer inPoints = input->GetPoints(); + vtkSmartPointer inNormals = nullptr; + vtkSmartPointeroutNormals = vtkFloatArray::New(); + outNormals->SetName("Normals"); + outNormals->SetNumberOfComponents(3); + + // determine point_data normals for the poly data points. + auto normalsGenerator = vtkSmartPointer::New(); + normalsGenerator->SetInputData(input); + normalsGenerator->FlipNormalsOn(); + normalsGenerator->Update(); + inNormals = normalsGenerator->GetOutput()->GetPointData()->GetNormals(); + input->GetPointData()->SetNormals(inNormals); + + auto outPoints = vtkSmartPointer::New(); + // auto outNormals = vtkSmartPointer::New(); + + if (input->GetMTime() > this->PreComp->mtime.GetMTime()) + { + vtkDebugMacro(<< "Compute edge-face connectivity and face normals\n"); + + this->PreComp->mtime.Modified(); + this->PreComp->edges.clear(); + + vtkIdType* polys = polysArray->GetPointer(0); + + for (vtkIdType i = 0; i < nPolys; i++) + { + int np = *(polys); ++polys; + double normal[3]; + vtkPolygon::ComputeNormal(inPoints, np, polys, normal); + + for (int j = 0; j < np; j++) + { + vtkIdType p1 = j, p2 = (j + 1) % np; + vtkOrderedEdge oe(polys[p1], polys[p2]); + vtkTwoNormals& tn = this->PreComp->edges[oe]; + if (polys[p1] < polys[p2]) + { +#ifdef DEBUG + if (vtkMath::Dot(tn.leftNormal, tn.leftNormal) > 0) + { + vtkDebugMacro(<< "Warning: vtkSilhouette: non-manifold mesh: edge-L (" << polys[p1] << "," << polys[p2] << ") counted more than once\n"); + } +#endif + tn.leftNormal[0] = normal[0]; + tn.leftNormal[1] = normal[1]; + tn.leftNormal[2] = normal[2]; + } + else + { +#ifdef DEBUG + if (vtkMath::Dot(tn.rightNormal, tn.rightNormal) > 0) + { + vtkDebugMacro(<< "Warning: vtkSilhouette: non-manifold mesh: edge-R (" << polys[p1] << "," << polys[p2] << ") counted more than once\n"); + } +#endif + tn.rightNormal[0] = normal[0]; + tn.rightNormal[1] = normal[1]; + tn.rightNormal[2] = normal[2]; + } + } + polys += np; + } + if (this->PreComp->edgeFlag != 0) delete[] this->PreComp->edgeFlag; + this->PreComp->edgeFlag = new bool[this->PreComp->edges.size()]; + } +} + +// Don't reference count to avoid nasty cycle +void vtkSilhouette::SetProp3D(vtkProp3D *prop3d) +{ + if ( this->Prop3D != prop3d ) + { + this->Prop3D = prop3d; + this->Modified(); + } +} + +vtkProp3D *vtkSilhouette::GetProp3D() +{ + return this->Prop3D; +} + +int vtkSilhouette::RequestData( + vtkInformation *vtkNotUsed(request), + vtkInformationVector **inputVector, + vtkInformationVector *outputVector) +{ + // get the info objects + vtkInformation *inInfo = inputVector[0]->GetInformationObject(0); + vtkInformation *outInfo = outputVector->GetInformationObject(0); + + // get the input and ouptut + vtkSmartPointer input = vtkPolyData::SafeDownCast( + inInfo->Get(vtkDataObject::DATA_OBJECT())); + vtkSmartPointeroutput = vtkPolyData::SafeDownCast( + outInfo->Get(vtkDataObject::DATA_OBJECT())); + auto helperPolyData = vtkSmartPointer::New(); + + if (input == 0 || output == 0) + { + vtkErrorMacro(<< "Need correct connections"); + return 0; + } + + vtkDebugMacro(<< "RequestData\n"); + + const double featureAngleCos = cos(vtkMath::RadiansFromDegrees(this->FeatureAngle)); + + + vtkIdType nPolys = input->GetNumberOfPolys(); + vtkSmartPointer polysArray = input->GetPolys()->GetData(); + vtkSmartPointer inPoints = input->GetPoints(); + vtkSmartPointer inNormals = input->GetPointData()->GetNormals(); + auto outNormals = vtkFloatArray::New(); + outNormals->SetName("Normals"); + outNormals->SetNumberOfComponents(3); + auto outPoints = vtkSmartPointer::New(); + bool vecChanged = false; + + + for(int d=0;d<3;d++) + { + vecChanged = vecChanged || this->PreComp->vec[d]!=CameraVector[d] ; + } + + if( ( this->PreComp->mtime.GetMTime() > output->GetMTime() ) || + ( this->Camera && this->Camera->GetMTime() > output->GetMTime() ) || + ( this->Prop3D && this->Prop3D->GetMTime() > output->GetMTime() ) || + vecChanged ) + { + vtkDebugMacro(<<"Extract edges\n"); + + vtkIdType i=0, silhouetteEdges=0; + + for(std::map::iterator it=this->PreComp->edges.begin(); it!=this->PreComp->edges.end(); ++it) + { + double d1,d2; + + // does this edge have two co-faces ? + bool winged = vtkMath::Norm(it->second.leftNormal)>0.5 && vtkMath::Norm(it->second.rightNormal)>0.5 ; + + // cosine of feature angle, to be compared with scalar product of two co-faces normals + double edgeAngleCos = vtkMath::Dot( it->second.leftNormal, it->second.rightNormal ); + + double p1[3]; + double p2[3]; + double vec[3]; + inPoints->GetPoint( it->first.p1 , p1 ); + inPoints->GetPoint( it->first.p2 , p2 ); + vec[0] = Origin[0] - ( (p1[0]+p2[0])*0.5 ); + vec[1] = Origin[1] - ( (p1[1]+p2[1])*0.5 ); + vec[2] = Origin[2] - ( (p1[2]+p2[2])*0.5 ); + d1 = vtkMath::Dot( vec, it->second.leftNormal ); + d2 = vtkMath::Dot( vec, it->second.rightNormal ); + + // shall we output this edge ? + bool outputEdge = + ( winged && (d1*d2)<0 ) + || ( this->EnableFeatureAngle && edgeAngleCosBorderEdges && !winged ); + + if( outputEdge ) // add this edge + { + this->PreComp->edgeFlag[i] = true ; + ++silhouetteEdges; + } + else // skip this edge + { + this->PreComp->edgeFlag[i] = false ; + } + ++i; + } + + // build output data set (lines) + vtkSmartPointer la = vtkIdTypeArray::New(); + la->SetNumberOfValues( 3*silhouetteEdges ); + vtkIdType* laPtr = la->WritePointer(0,3*silhouetteEdges); + + i=0; + silhouetteEdges=0; + for(std::map::iterator it=this->PreComp->edges.begin(); it!=this->PreComp->edges.end(); ++it) + { + if( this->PreComp->edgeFlag[i] ) + { + laPtr[ silhouetteEdges*3+0 ] = 2 ; + laPtr[ silhouetteEdges*3+1 ] = it->first.p1 ; + laPtr[ silhouetteEdges*3+2 ] = it->first.p2 ; + ++silhouetteEdges; + outPoints->InsertPoint(it->first.p1,inPoints->GetPoint(it->first.p1)); + outPoints->InsertPoint(it->first.p2,inPoints->GetPoint(it->first.p2)); + double normal[3]; + inNormals->GetTuple(it->first.p1, normal); + outNormals->InsertTuple3(it->first.p1, normal[0], normal[1], normal[2]); + inNormals->GetTuple(it->first.p2, normal); + outNormals->InsertTuple3(it->first.p2, normal[0], normal[1], normal[2]); + } + ++i; + } + + if( this->PreComp->lines == 0 ) + { + this->PreComp->lines = vtkCellArray::New(); + } + this->PreComp->lines->SetCells( silhouetteEdges, la ); + la->Delete(); + } + + helperPolyData->Initialize(); + helperPolyData->SetPoints(inPoints); + helperPolyData->SetLines(this->PreComp->lines); + helperPolyData->GetPointData()->SetNormals(outNormals); + auto cleanPolyDataFilter = vtkSmartPointer::New(); + cleanPolyDataFilter->SetInputData(helperPolyData); + cleanPolyDataFilter->Update(); + output->Initialize(); + output->DeepCopy(cleanPolyDataFilter->GetOutput()); + + return 1; +} + + +//unsigned long int vtkSilhouette::GetMTime() +//{ +// unsigned long mTime=this->Superclass::GetMTime(); +// +// if ( this->Direction != VTK_DIRECTION_SPECIFIED_VECTOR ) +// { +// unsigned long time; +// if ( this->Camera != NULL ) +// { +// time = this->Camera->GetMTime(); +// mTime = ( time > mTime ? time : mTime ); +// } +// +// if ( this->Prop3D != NULL ) +// { +// time = this->Prop3D->GetMTime(); +// mTime = ( time > mTime ? time : mTime ); +// } +// } +// +// return mTime; +//} + +void vtkSilhouette::PrintSelf(ostream& os, vtkIndent indent) +{ + + this->Superclass::PrintSelf(os,indent); + + if ( this->Camera ) + { + os << indent << "Camera:\n"; + this->Camera->PrintSelf(os,indent.GetNextIndent()); + } + else + { + os << indent << "Camera: (none)\n"; + } + + if ( this->Prop3D ) + { + os << indent << "Prop3D:\n"; + this->Prop3D->PrintSelf(os,indent.GetNextIndent()); + } + else + { + os << indent << "Prop3D: (none)\n"; + } + + os << indent << "Direction: "; +#define DIRECTION_CASE(name) case VTK_DIRECTION_##name :os << "VTK_DIRECTION_" << #name <<"\n"; break + switch(this->Direction) + { + DIRECTION_CASE( SPECIFIED_ORIGIN ); + DIRECTION_CASE( SPECIFIED_VECTOR ); + DIRECTION_CASE( CAMERA_ORIGIN ); + DIRECTION_CASE( CAMERA_VECTOR ); + } +#undef DIRECTION_CASE + + if( this->Direction == VTK_DIRECTION_SPECIFIED_VECTOR ) + { + os << "Specified Vector: (" << this->Vector[0] << ", " << this->Vector[1] << ", " << this->Vector[2] << ")\n"; + } + if( this->Direction == VTK_DIRECTION_SPECIFIED_ORIGIN ) + { + os << "Specified Origin: (" << this->Origin[0] << ", " << this->Origin[1] << ", " << this->Origin[2] << ")\n"; + } + + os << indent << "PieceInvariant: " << this->PieceInvariant << "\n"; + os << indent << "FeatureAngle: " << this->FeatureAngle << "\n"; + os << indent << "EnableFeatureAngle: " << this->EnableFeatureAngle << "\n"; + os << indent << "BorderEdges: " << this->BorderEdges << "\n"; +} + diff --git a/Plugins/org.mitk.gui.qt.standardplanetool/src/internal/vtkSilhouette.h b/Plugins/org.mitk.gui.qt.standardplanetool/src/internal/vtkSilhouette.h new file mode 100644 index 0000000000..b88b22ec24 --- /dev/null +++ b/Plugins/org.mitk.gui.qt.standardplanetool/src/internal/vtkSilhouette.h @@ -0,0 +1,172 @@ +/*========================================================================= + + Program: Visualization Toolkit + Module: vtkPolyDataSilhouette.h + + Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen + All rights reserved. + See Copyright.txt or http://www.kitware.com/Copyright.htm for details. + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notice for more information. + +=========================================================================*/ +// .NAME vtkPolyDataSilhouette - sort polydata along camera view direction +// +// .SECTION Description +// vtkPolyDataSilhouette extracts a subset of a polygonal mesh edges to +// generate an outline (silhouette) of the corresponding 3D object. In +// addition, this filter can also extracts sharp edges (aka feature angles). +// In order to use this filter you must specify the a point of view (origin) or +// a direction (vector). given this direction or origin, a silhouette is +// generated wherever the surface's normal is orthogonal to the view direction. +// +// .SECTION Caveats +// when the active camera is used, almost everything is recomputed for each +// frame, keep this in mind when dealing with extremely large surface data +// sets. +// +// .SECTION Thanks +// Contribution by Thierry Carrard
+// CEA/DIF - Commissariat a l'Energie Atomique, Centre DAM Ile-De-France
+// BP12, F-91297 Arpajon, France.
+ +#ifndef vtkSilhouette_h +#define vtkSilhouette_h + +#include "vtkPolyDataAlgorithm.h" + +#define VTK_DIRECTION_SPECIFIED_VECTOR 0 +#define VTK_DIRECTION_SPECIFIED_ORIGIN 1 +#define VTK_DIRECTION_CAMERA_ORIGIN 2 +#define VTK_DIRECTION_CAMERA_VECTOR 3 + +class vtkCamera; +class vtkProp3D; +class vtkTransform; +class vtkPolyDataEdges; +class vtkMatrix4x4; + +class vtkSilhouette : public vtkPolyDataAlgorithm +{ +public: + // Description: + // Instantiate object. + static vtkSilhouette *New(); + + vtkTypeMacro(vtkSilhouette,vtkPolyDataAlgorithm); + void PrintSelf(ostream& os, vtkIndent indent); + + // Description: + // Enables or Disables generation of silhouette edges along sharp edges + vtkSetMacro(EnableFeatureAngle,int); + vtkGetMacro(EnableFeatureAngle,int); + + // Description: + // Sets/Gets minimal angle for sharp edges detection. Default is 60 + vtkSetMacro(FeatureAngle,double); + vtkGetMacro(FeatureAngle,double); + + // Description: + // Enables or Disables generation of border edges. Note: borders exist only + // in case of non closed surface + vtkSetMacro(BorderEdges,int); + vtkGetMacro(BorderEdges,int); + vtkBooleanMacro(BorderEdges,int); + + // Description: + // Enables or Disables piece invariance. This is useful when dealing with + // multi-block data sets. Note: requires one level of ghost cells + vtkSetMacro(PieceInvariant,int); + vtkGetMacro(PieceInvariant,int); + vtkBooleanMacro(PieceInvariant,int); + + // Description: + // Specify how view direction is computed. By default, the + // camera origin (eye) is used. + vtkSetMacro(Direction,int); + vtkGetMacro(Direction,int); + void SetDirectionToSpecifiedVector() + {this->SetDirection( VTK_DIRECTION_SPECIFIED_VECTOR ); } + void SetDirectionToSpecifiedOrigin() + {this->SetDirection( VTK_DIRECTION_SPECIFIED_ORIGIN ); } + void SetDirectionToCameraVector() + {this->SetDirection( VTK_DIRECTION_CAMERA_VECTOR ); } + void SetDirectionToCameraOrigin() + {this->SetDirection( VTK_DIRECTION_CAMERA_ORIGIN ); } + + // Description: + // Specify a camera that is used to define the view direction. This ivar + // only has effect if the direction is set to VTK_DIRECTION_CAMERA_ORIGIN or + // VTK_DIRECTION_CAMERA_VECTOR, and a camera is specified. + virtual void SetCamera(vtkCamera VTK_WRAP_EXTERN*); + vtkGetObjectMacro(Camera,vtkCamera VTK_WRAP_EXTERN); + + + void SetProjectionMatrix(vtkMatrix4x4* matrix) + { + this->ProjectionMatrix = matrix; + }; + vtkGetObjectMacro(ProjectionMatrix, vtkMatrix4x4); + + void InitializeProjectionGeometry(); + void InitializeMesh(); + + // Description: + // Specify a transformation matrix (via the vtkProp3D::GetMatrix() method) + // that is used to include the effects of transformation. This ivar only has + // effect if the direction is set to VTK_DIRECTION_CAMERA_ORIGIN or + // VTK_DIRECTION_CAMERA_VECTOR, and a camera is specified. Specifying the + // vtkProp3D is optional. + void SetProp3D(vtkProp3D VTK_WRAP_EXTERN*); + vtkProp3D VTK_WRAP_EXTERN*GetProp3D(); + + // Description: + // Set/Get the sort direction. This ivar only has effect if the sort + // direction is set to SetDirectionToSpecifiedVector(). The edge detection + // occurs in the direction of the vector. + vtkSetVector3Macro(Vector,double); + vtkGetVectorMacro(Vector,double,3); + + // Description: + // Set/Get the sort origin. This ivar only has effect if the sort direction + // is set to SetDirectionToSpecifiedOrigin(). The edge detection occurs in + // the direction of the origin to each edge's center. + vtkSetVector3Macro(Origin,double); + vtkGetVectorMacro(Origin,double,3); + + // Description: + // Return MTime also considering the dependent objects: the camera + // and/or the prop3D. +// unsigned long GetMTime(); + +protected: + vtkSilhouette(); + ~vtkSilhouette(); + + int RequestData(vtkInformation *, vtkInformationVector **, vtkInformationVector *); + + int Direction; + vtkCamera *Camera; + vtkProp3D *Prop3D; + vtkTransform *Transform; + vtkMatrix4x4 *ProjectionMatrix; + double Vector[3]; + double Origin[3]; + double CameraVector[3]; + + int EnableFeatureAngle; + double FeatureAngle; + + int BorderEdges; + int PieceInvariant; + + vtkPolyDataEdges* PreComp; // precomputed data for a given point of view + +private: + vtkSilhouette(const vtkSilhouette&); // Not implemented. + void operator=(const vtkSilhouette&); // Not implemented. +}; + +#endif