diff --git a/CMake/BuildConfigurations/DiffusionRelease.cmake b/CMake/BuildConfigurations/DiffusionRelease.cmake index c1e855e2ec..67ceff6737 100644 --- a/CMake/BuildConfigurations/DiffusionRelease.cmake +++ b/CMake/BuildConfigurations/DiffusionRelease.cmake @@ -1,75 +1,76 @@ message(STATUS "Configuring MITK Diffusion Release Build") # Enable non-optional external dependencies set(MITK_USE_Vigra ON CACHE BOOL "MITK Use Vigra Library" FORCE) set(MITK_USE_HDF5 ON CACHE BOOL "MITK Use HDF5 Library" FORCE) # Disable all apps but MITK Diffusion set(MITK_BUILD_ALL_APPS OFF CACHE BOOL "Build all MITK applications" FORCE) set(MITK_BUILD_APP_CoreApp OFF CACHE BOOL "Build the MITK CoreApp" FORCE) set(MITK_BUILD_APP_Workbench OFF CACHE BOOL "Build the MITK Workbench" FORCE) set(MITK_BUILD_APP_Fiberfox OFF CACHE BOOL "Build MITK Fiberfox" FORCE) set(MITK_BUILD_APP_Diffusion ON CACHE BOOL "Build MITK Diffusion" FORCE) # Activate Diffusion Mini Apps set(BUILD_DiffusionCoreCmdApps OFF CACHE BOOL "Build commandline tools for diffusion" FORCE) set(BUILD_DiffusionFiberProcessingCmdApps ON CACHE BOOL "Build commandline tools for diffusion fiber processing" FORCE) set(BUILD_DiffusionFiberfoxCmdApps ON CACHE BOOL "Build commandline tools for diffusion data simulation (Fiberfox)" FORCE) set(BUILD_DiffusionMiscCmdApps OFF CACHE BOOL "Build miscellaneous commandline tools for diffusion" FORCE) set(BUILD_DiffusionQuantificationCmdApps OFF CACHE BOOL "Build commandline tools for diffusion quantification (IVIM, ADC, ...)" FORCE) set(BUILD_DiffusionTractographyCmdApps ON CACHE BOOL "Build commandline tools for diffusion fiber tractography" FORCE) set(BUILD_DiffusionTractographyEvaluationCmdApps OFF CACHE BOOL "Build commandline tools for diffusion fiber tractography evaluation" FORCE) -set(BUILD_DiffusionConnectomicsCmdApps OFF CACHE BOOL "Build commandline tools for diffusion connectomics" FORCE) +set(BUILD_DiffusionConnectomicsCmdApps ON CACHE BOOL "Build commandline tools for diffusion connectomics" FORCE) # Build neither all plugins nor examples set(MITK_BUILD_ALL_PLUGINS OFF CACHE BOOL "Build all MITK plugins" FORCE) set(MITK_BUILD_EXAMPLES OFF CACHE BOOL "Build the MITK examples" FORCE) set(BUILD_TESTING OFF CACHE BOOL "Build the MITK tests" FORCE) # Activate in-application help generation set(MITK_DOXYGEN_GENERATE_QCH_FILES ON CACHE BOOL "Use doxygen to generate Qt compressed help files for MITK docs" FORCE) set(BLUEBERRY_USE_QT_HELP ON CACHE BOOL "Enable support for integrating bundle documentation into Qt Help" FORCE) # Disable console window set(MITK_SHOW_CONSOLE_WINDOW ON CACHE BOOL "Use this to enable or disable the console window when starting MITK GUI Applications" FORCE) # Activate required plugins # should be identical to the list in /Applications/Diffusion/CMakeLists.txt set(_plugins org.commontk.configadmin org.commontk.eventadmin org.blueberry.compat org.blueberry.core.runtime org.blueberry.core.expressions org.blueberry.core.commands org.blueberry.ui.qt org.blueberry.ui.qt.log org.blueberry.ui.qt.help org.mitk.core.services org.mitk.gui.common org.mitk.planarfigure org.mitk.core.ext org.mitk.gui.qt.application org.mitk.gui.qt.ext org.mitk.gui.qt.diffusionimagingapp org.mitk.gui.qt.common org.mitk.gui.qt.stdmultiwidgeteditor org.mitk.gui.qt.common.legacy org.mitk.gui.qt.datamanager org.mitk.gui.qt.measurementtoolbox org.mitk.gui.qt.segmentation org.mitk.gui.qt.volumevisualization org.mitk.gui.qt.diffusionimaging org.mitk.gui.qt.diffusionimaging.fiberfox org.mitk.gui.qt.diffusionimaging.fiberprocessing org.mitk.gui.qt.diffusionimaging.ivim org.mitk.gui.qt.diffusionimaging.odfpeaks org.mitk.gui.qt.diffusionimaging.preprocessing org.mitk.gui.qt.diffusionimaging.reconstruction org.mitk.gui.qt.diffusionimaging.tractography + org.mitk.gui.qt.diffusionimaging.connectomics org.mitk.gui.qt.imagenavigator org.mitk.gui.qt.moviemaker org.mitk.gui.qt.basicimageprocessing org.mitk.gui.qt.properties org.mitk.gui.qt.viewnavigator ) diff --git a/Modules/DiffusionImaging/Connectomics/Algorithms/itkConnectomicsNetworkToConnectivityMatrixImageFilter.cpp b/Modules/DiffusionImaging/Connectomics/Algorithms/itkConnectomicsNetworkToConnectivityMatrixImageFilter.cpp index 3198e34caf..33e6285f00 100644 --- a/Modules/DiffusionImaging/Connectomics/Algorithms/itkConnectomicsNetworkToConnectivityMatrixImageFilter.cpp +++ b/Modules/DiffusionImaging/Connectomics/Algorithms/itkConnectomicsNetworkToConnectivityMatrixImageFilter.cpp @@ -1,186 +1,186 @@ /*=================================================================== 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 ITK_ConnectomicsNetworkToConnectivityMatrixImageFilter_CPP #define ITK_ConnectomicsNetworkToConnectivityMatrixImageFilter_CPP #include #include itk::ConnectomicsNetworkToConnectivityMatrixImageFilter::ConnectomicsNetworkToConnectivityMatrixImageFilter() : m_BinaryConnectivity(false) , m_RescaleConnectivity(false) , m_InputNetwork(nullptr) { } itk::ConnectomicsNetworkToConnectivityMatrixImageFilter::~ConnectomicsNetworkToConnectivityMatrixImageFilter() { } void itk::ConnectomicsNetworkToConnectivityMatrixImageFilter::GenerateData() { this->AllocateOutputs(); OutputImageType::Pointer output = this->GetOutput(); if(m_InputNetwork.IsNull()) { MITK_ERROR << "Failed to generate data, no valid network was set."; return; } typedef mitk::ConnectomicsNetwork::NetworkType NetworkType; typedef boost::graph_traits< NetworkType >::vertex_descriptor DescriptorType; typedef boost::graph_traits< NetworkType >::vertex_iterator IteratorType; // prepare connectivity matrix for data std::vector< std::vector< int > > connectivityMatrix; int numberOfVertices = m_InputNetwork->GetNumberOfVertices(); connectivityMatrix.resize( numberOfVertices ); for( int index(0); index < numberOfVertices; index++ ) { connectivityMatrix[ index ].resize( numberOfVertices ); } // Create LabelToIndex translation std::map< std::string, DescriptorType > labelToIndexMap; NetworkType boostGraph = *(m_InputNetwork->GetBoostGraph()); // translate index to label IteratorType iterator, end; boost::tie(iterator, end) = boost::vertices( boostGraph ); for ( ; iterator != end; ++iterator) { labelToIndexMap.insert( std::pair< std::string, DescriptorType >( boostGraph[ *iterator ].label, *iterator ) ); } std::vector< std::string > indexToLabel; // translate index to label indexToLabel.resize( numberOfVertices ); boost::tie(iterator, end) = boost::vertices( boostGraph ); for ( ; iterator != end; ++iterator) { indexToLabel[ *iterator ] = boostGraph[ *iterator ].label; } //translate label to position std::vector< std::string > positionToLabelVector; positionToLabelVector = indexToLabel; std::sort ( positionToLabelVector.begin(), positionToLabelVector.end() ); for( int outerLoop( 0 ); outerLoop < numberOfVertices; outerLoop++ ) { DescriptorType fromVertexDescriptor = labelToIndexMap.find( positionToLabelVector[ outerLoop ] )->second; for( int innerLoop( outerLoop + 1 ); innerLoop < numberOfVertices; innerLoop++ ) { DescriptorType toVertexDescriptor = labelToIndexMap.find( positionToLabelVector[ innerLoop ] )->second; int weight( 0 ); if( boost::edge(toVertexDescriptor, fromVertexDescriptor, boostGraph ).second ) { weight = m_InputNetwork->GetEdge( fromVertexDescriptor , toVertexDescriptor ).weight;; } connectivityMatrix[outerLoop][innerLoop] = weight; connectivityMatrix[innerLoop][outerLoop] = weight; } } OutputImageType::SpacingType spacing; spacing[0] = 1.0; spacing[1] = 1.0; OutputImageType::PointType origin; origin[0] = 0.0; origin[1] = 0.0; OutputImageType::IndexType index; index[0] = 0.0; index[1] = 0.0; OutputImageType::SizeType size; size[0] = numberOfVertices; size[1] = numberOfVertices; OutputImageType::RegionType region; region.SetIndex( index ); region.SetSize( size ); output->SetSpacing( spacing ); // Set the image spacing output->SetOrigin( origin ); // Set the image origin output->SetRegions( region ); output->Allocate(); output->FillBuffer(0); itk::ImageRegionIterator< OutputImageType > it_connect(output, output->GetLargestPossibleRegion()); int counter( 0 ); double rescaleFactor( 1.0 ); double rescaleMax( 255.0 ); if( m_RescaleConnectivity ) { rescaleFactor = rescaleMax / m_InputNetwork->GetMaximumWeight(); } // Colour pixels according to connectivity if( m_BinaryConnectivity ) { // binary connectivity while( !it_connect.IsAtEnd() ) { if( connectivityMatrix[ ( counter - counter % numberOfVertices ) / numberOfVertices][ counter % numberOfVertices ] ) { it_connect.Set( 1 ); } else { it_connect.Set( 0 ); } ++it_connect; counter++; } } else { // if desired rescale to the 0-255 range while( !it_connect.IsAtEnd() ) { - it_connect.Set( ( unsigned short ) rescaleFactor * - connectivityMatrix[ ( counter - counter % numberOfVertices ) / numberOfVertices][ counter % numberOfVertices ] + it_connect.Set( ( unsigned short ) (rescaleFactor * + connectivityMatrix[ ( counter - counter % numberOfVertices ) / numberOfVertices][ counter % numberOfVertices ]) ); ++it_connect; counter++; } } } #endif diff --git a/Modules/DiffusionImaging/DiffusionIO/mitkFiberBundleVtkWriter.cpp b/Modules/DiffusionImaging/DiffusionIO/mitkFiberBundleVtkWriter.cpp index af04976873..f0599994c3 100644 --- a/Modules/DiffusionImaging/DiffusionIO/mitkFiberBundleVtkWriter.cpp +++ b/Modules/DiffusionImaging/DiffusionIO/mitkFiberBundleVtkWriter.cpp @@ -1,155 +1,172 @@ /*=================================================================== 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 "mitkFiberBundleVtkWriter.h" #include #include #include #include #include #include #include #include #include #include #include #include "mitkDiffusionIOMimeTypes.h" mitk::FiberBundleVtkWriter::FiberBundleVtkWriter() - : mitk::AbstractFileWriter(mitk::FiberBundle::GetStaticNameOfClass(), mitk::DiffusionIOMimeTypes::FIBERBUNDLE_VTK_MIMETYPE_NAME(), "VTK Fiber Bundle Writer") + : mitk::AbstractFileWriter(mitk::FiberBundle::GetStaticNameOfClass(), mitk::DiffusionIOMimeTypes::FIBERBUNDLE_VTK_MIMETYPE_NAME(), "VTK Fiber Bundle Writer") { - Options defaultOptions; - defaultOptions["Save as binary file"] = true; - defaultOptions["Save as xml file (vtp style)"] = false; - defaultOptions["Save color information"] = false; - defaultOptions["Save fiber weights"] = true; - this->SetDefaultOptions(defaultOptions); - RegisterService(); + Options defaultOptions; + defaultOptions["Save as binary file"] = true; + defaultOptions["Save as xml file (vtp style)"] = false; + defaultOptions["Save color information"] = false; + defaultOptions["Save fiber weights"] = true; + this->SetDefaultOptions(defaultOptions); + RegisterService(); } mitk::FiberBundleVtkWriter::FiberBundleVtkWriter(const mitk::FiberBundleVtkWriter & other) - :mitk::AbstractFileWriter(other) + :mitk::AbstractFileWriter(other) {} mitk::FiberBundleVtkWriter::~FiberBundleVtkWriter() {} mitk::FiberBundleVtkWriter * mitk::FiberBundleVtkWriter::Clone() const { - return new mitk::FiberBundleVtkWriter(*this); + return new mitk::FiberBundleVtkWriter(*this); } void mitk::FiberBundleVtkWriter::Write() { - std::ostream* out; - std::ofstream outStream; + std::ostream* out; + std::ofstream outStream; - if( this->GetOutputStream() ) + if( this->GetOutputStream() ) + { + out = this->GetOutputStream(); + }else{ + outStream.open( this->GetOutputLocation().c_str() ); + out = &outStream; + } + + if ( !out->good() ) + { + mitkThrow() << "Stream not good."; + } + + try + { + const std::string& locale = "C"; + const std::string& currLocale = setlocale( LC_ALL, nullptr ); + setlocale(LC_ALL, locale.c_str()); + + std::locale previousLocale(out->getloc()); + std::locale I("C"); + out->imbue(I); + + std::string filename = this->GetOutputLocation().c_str(); + + mitk::FiberBundle::ConstPointer input = dynamic_cast(this->GetInput()); + std::string ext = itksys::SystemTools::GetFilenameLastExtension(this->GetOutputLocation().c_str()); + Options options = this->GetOptions(); + + vtkSmartPointer weights = nullptr; + vtkSmartPointer fiberColors = nullptr; + + vtkSmartPointer fibPoly = input->GetFiberPolyData(); + if (us::any_cast(options["Save fiber weights"])) + { + MITK_INFO << "Adding fiber weight information"; + fibPoly->GetCellData()->AddArray(input->GetFiberWeights()); + } + else if (fibPoly->GetCellData()->HasArray("FIBER_WEIGHTS")) { - out = this->GetOutputStream(); - }else{ - outStream.open( this->GetOutputLocation().c_str() ); - out = &outStream; + weights = input->GetFiberWeights(); + fibPoly->GetCellData()->RemoveArray("FIBER_WEIGHTS"); } - if ( !out->good() ) + if (us::any_cast(options["Save color information"])) { - mitkThrow() << "Stream not good."; + MITK_INFO << "Adding color information"; + fibPoly->GetPointData()->AddArray(input->GetFiberColors()); + } + else if (fibPoly->GetPointData()->HasArray("FIBER_COLORS")) + { + fiberColors = input->GetFiberColors(); + fibPoly->GetPointData()->RemoveArray("FIBER_COLORS"); + } + + // default extension is .fib + if(ext == "") + { + ext = ".fib"; + this->SetOutputLocation(this->GetOutputLocation() + ext); } - try + if (us::any_cast(options["Save as xml file (vtp style)"])) { - const std::string& locale = "C"; - const std::string& currLocale = setlocale( LC_ALL, nullptr ); - setlocale(LC_ALL, locale.c_str()); - - std::locale previousLocale(out->getloc()); - std::locale I("C"); - out->imbue(I); - - std::string filename = this->GetOutputLocation().c_str(); - - mitk::FiberBundle::ConstPointer input = dynamic_cast(this->GetInput()); - std::string ext = itksys::SystemTools::GetFilenameLastExtension(this->GetOutputLocation().c_str()); - Options options = this->GetOptions(); - - vtkSmartPointer fibPoly = input->GetFiberPolyData(); - if (us::any_cast(options["Save fiber weights"])) - { - MITK_INFO << "Adding fiber weight information"; - fibPoly->GetCellData()->AddArray(input->GetFiberWeights()); - } - else if (fibPoly->GetCellData()->HasArray("FIBER_WEIGHTS")) - fibPoly->GetCellData()->RemoveArray("FIBER_WEIGHTS"); - if (us::any_cast(options["Save color information"])) - { - MITK_INFO << "Adding color information"; - fibPoly->GetPointData()->AddArray(input->GetFiberColors()); - } - else if (fibPoly->GetPointData()->HasArray("FIBER_COLORS")) - fibPoly->GetPointData()->RemoveArray("FIBER_COLORS"); - - // default extension is .fib - if(ext == "") - { - ext = ".fib"; - this->SetOutputLocation(this->GetOutputLocation() + ext); - } - - if (us::any_cast(options["Save as xml file (vtp style)"])) - { - vtkSmartPointer writer = vtkSmartPointer::New(); - writer->SetInputData(fibPoly); - writer->SetFileName(filename.c_str()); - if (us::any_cast(options["Save as binary file"])) - { - MITK_INFO << "Writing fiber bundle as vtk binary file"; - writer->SetDataModeToBinary(); - } - else - { - MITK_INFO << "Writing fiber bundle as vtk ascii file"; - writer->SetDataModeToAscii(); - } - writer->Write(); - } - else - { - vtkSmartPointer writer = vtkSmartPointer::New(); - writer->SetInputData(fibPoly); - writer->SetFileName(filename.c_str()); - if (us::any_cast(options["Save as binary file"])) - { - MITK_INFO << "Writing fiber bundle as vtk binary file"; - writer->SetFileTypeToBinary(); - } - else - { - MITK_INFO << "Writing fiber bundle as vtk ascii file"; - writer->SetFileTypeToASCII(); - } - writer->Write(); - } - - setlocale(LC_ALL, currLocale.c_str()); - MITK_INFO << "Fiber bundle written"; + vtkSmartPointer writer = vtkSmartPointer::New(); + writer->SetInputData(fibPoly); + writer->SetFileName(filename.c_str()); + if (us::any_cast(options["Save as binary file"])) + { + MITK_INFO << "Writing fiber bundle as vtk binary file"; + writer->SetDataModeToBinary(); + } + else + { + MITK_INFO << "Writing fiber bundle as vtk ascii file"; + writer->SetDataModeToAscii(); + } + writer->Write(); } - catch(...) + else { - throw; + vtkSmartPointer writer = vtkSmartPointer::New(); + writer->SetInputData(fibPoly); + writer->SetFileName(filename.c_str()); + if (us::any_cast(options["Save as binary file"])) + { + MITK_INFO << "Writing fiber bundle as vtk binary file"; + writer->SetFileTypeToBinary(); + } + else + { + MITK_INFO << "Writing fiber bundle as vtk ascii file"; + writer->SetFileTypeToASCII(); + } + writer->Write(); } + + // restore arrays + if ( weights!=nullptr ) + fibPoly->GetPointData()->AddArray(weights); + + if (fiberColors != nullptr) + fibPoly->GetPointData()->AddArray(fiberColors); + + setlocale(LC_ALL, currLocale.c_str()); + MITK_INFO << "Fiber bundle written"; + } + catch(...) + { + throw; + } } diff --git a/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundle/mitkFiberBundle.cpp b/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundle/mitkFiberBundle.cpp index 7ce5700324..af15024329 100755 --- a/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundle/mitkFiberBundle.cpp +++ b/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundle/mitkFiberBundle.cpp @@ -1,2224 +1,2283 @@ /*=================================================================== 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. ===================================================================*/ #define _USE_MATH_DEFINES #include "mitkFiberBundle.h" #include #include #include #include "mitkImagePixelReadAccessor.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include const char* mitk::FiberBundle::FIBER_ID_ARRAY = "Fiber_IDs"; using namespace std; mitk::FiberBundle::FiberBundle( vtkPolyData* fiberPolyData ) - : m_NumFibers(0) + : m_NumFibers(0) { - m_FiberWeights = vtkSmartPointer::New(); - m_FiberWeights->SetName("FIBER_WEIGHTS"); + m_FiberWeights = vtkSmartPointer::New(); + m_FiberWeights->SetName("FIBER_WEIGHTS"); - m_FiberPolyData = vtkSmartPointer::New(); - if (fiberPolyData != nullptr) - m_FiberPolyData = fiberPolyData; + m_FiberPolyData = vtkSmartPointer::New(); + if (fiberPolyData != nullptr) + m_FiberPolyData = fiberPolyData; - this->UpdateFiberGeometry(); - this->GenerateFiberIds(); - this->ColorFibersByOrientation(); + this->UpdateFiberGeometry(); + this->GenerateFiberIds(); + this->ColorFibersByOrientation(); } mitk::FiberBundle::~FiberBundle() { } mitk::FiberBundle::Pointer mitk::FiberBundle::GetDeepCopy() { - mitk::FiberBundle::Pointer newFib = mitk::FiberBundle::New(m_FiberPolyData); - newFib->SetFiberColors(this->m_FiberColors); - newFib->SetFiberWeights(this->m_FiberWeights); - return newFib; + mitk::FiberBundle::Pointer newFib = mitk::FiberBundle::New(m_FiberPolyData); + newFib->SetFiberColors(this->m_FiberColors); + newFib->SetFiberWeights(this->m_FiberWeights); + return newFib; } vtkSmartPointer mitk::FiberBundle::GeneratePolyDataByIds(std::vector fiberIds) { - vtkSmartPointer newFiberPolyData = vtkSmartPointer::New(); - vtkSmartPointer newLineSet = vtkSmartPointer::New(); - vtkSmartPointer newPointSet = vtkSmartPointer::New(); - - auto finIt = fiberIds.begin(); - while ( finIt != fiberIds.end() ) - { - if (*finIt < 0 || *finIt>GetNumFibers()){ - MITK_INFO << "FiberID can not be negative or >NumFibers!!! check id Extraction!" << *finIt; - break; - } - - vtkSmartPointer fiber = m_FiberIdDataSet->GetCell(*finIt);//->DeepCopy(fiber); - vtkSmartPointer fibPoints = fiber->GetPoints(); - vtkSmartPointer newFiber = vtkSmartPointer::New(); - newFiber->GetPointIds()->SetNumberOfIds( fibPoints->GetNumberOfPoints() ); + vtkSmartPointer newFiberPolyData = vtkSmartPointer::New(); + vtkSmartPointer newLineSet = vtkSmartPointer::New(); + vtkSmartPointer newPointSet = vtkSmartPointer::New(); + + auto finIt = fiberIds.begin(); + while ( finIt != fiberIds.end() ) + { + if (*finIt < 0 || *finIt>GetNumFibers()){ + MITK_INFO << "FiberID can not be negative or >NumFibers!!! check id Extraction!" << *finIt; + break; + } - for(int i=0; iGetNumberOfPoints(); i++) - { - newFiber->GetPointIds()->SetId(i, newPointSet->GetNumberOfPoints()); - newPointSet->InsertNextPoint(fibPoints->GetPoint(i)[0], fibPoints->GetPoint(i)[1], fibPoints->GetPoint(i)[2]); - } + vtkSmartPointer fiber = m_FiberIdDataSet->GetCell(*finIt);//->DeepCopy(fiber); + vtkSmartPointer fibPoints = fiber->GetPoints(); + vtkSmartPointer newFiber = vtkSmartPointer::New(); + newFiber->GetPointIds()->SetNumberOfIds( fibPoints->GetNumberOfPoints() ); - newLineSet->InsertNextCell(newFiber); - ++finIt; + for(int i=0; iGetNumberOfPoints(); i++) + { + newFiber->GetPointIds()->SetId(i, newPointSet->GetNumberOfPoints()); + newPointSet->InsertNextPoint(fibPoints->GetPoint(i)[0], fibPoints->GetPoint(i)[1], fibPoints->GetPoint(i)[2]); } - newFiberPolyData->SetPoints(newPointSet); - newFiberPolyData->SetLines(newLineSet); - return newFiberPolyData; + newLineSet->InsertNextCell(newFiber); + ++finIt; + } + + newFiberPolyData->SetPoints(newPointSet); + newFiberPolyData->SetLines(newLineSet); + return newFiberPolyData; } // merge two fiber bundles mitk::FiberBundle::Pointer mitk::FiberBundle::AddBundle(mitk::FiberBundle* fib) { - if (fib==nullptr) - { - MITK_WARN << "trying to call AddBundle with nullptr argument"; - return nullptr; - } - MITK_INFO << "Adding fibers"; - - vtkSmartPointer vNewPolyData = vtkSmartPointer::New(); - vtkSmartPointer vNewLines = vtkSmartPointer::New(); - vtkSmartPointer vNewPoints = vtkSmartPointer::New(); - - // add current fiber bundle - vtkSmartPointer weights = vtkSmartPointer::New(); - weights->SetNumberOfValues(this->GetNumFibers()+fib->GetNumFibers()); - - unsigned int counter = 0; - for (int i=0; iGetNumberOfCells(); i++) + if (fib==nullptr) + { + MITK_WARN << "trying to call AddBundle with nullptr argument"; + return nullptr; + } + MITK_INFO << "Adding fibers"; + + vtkSmartPointer vNewPolyData = vtkSmartPointer::New(); + vtkSmartPointer vNewLines = vtkSmartPointer::New(); + vtkSmartPointer vNewPoints = vtkSmartPointer::New(); + + // add current fiber bundle + vtkSmartPointer weights = vtkSmartPointer::New(); + weights->SetNumberOfValues(this->GetNumFibers()+fib->GetNumFibers()); + + unsigned int counter = 0; + for (int i=0; iGetNumberOfCells(); i++) + { + vtkCell* cell = m_FiberPolyData->GetCell(i); + int numPoints = cell->GetNumberOfPoints(); + vtkPoints* points = cell->GetPoints(); + + vtkSmartPointer container = vtkSmartPointer::New(); + for (int j=0; jGetCell(i); - int numPoints = cell->GetNumberOfPoints(); - vtkPoints* points = cell->GetPoints(); - - vtkSmartPointer container = vtkSmartPointer::New(); - for (int j=0; jGetPoint(j, p); + double p[3]; + points->GetPoint(j, p); - vtkIdType id = vNewPoints->InsertNextPoint(p); - container->GetPointIds()->InsertNextId(id); - } - weights->InsertValue(counter, this->GetFiberWeight(i)); - vNewLines->InsertNextCell(container); - counter++; + vtkIdType id = vNewPoints->InsertNextPoint(p); + container->GetPointIds()->InsertNextId(id); } - - // add new fiber bundle - for (int i=0; iGetFiberPolyData()->GetNumberOfCells(); i++) + weights->InsertValue(counter, this->GetFiberWeight(i)); + vNewLines->InsertNextCell(container); + counter++; + } + + // add new fiber bundle + for (int i=0; iGetFiberPolyData()->GetNumberOfCells(); i++) + { + vtkCell* cell = fib->GetFiberPolyData()->GetCell(i); + int numPoints = cell->GetNumberOfPoints(); + vtkPoints* points = cell->GetPoints(); + + vtkSmartPointer container = vtkSmartPointer::New(); + for (int j=0; jGetFiberPolyData()->GetCell(i); - int numPoints = cell->GetNumberOfPoints(); - vtkPoints* points = cell->GetPoints(); + double p[3]; + points->GetPoint(j, p); - vtkSmartPointer container = vtkSmartPointer::New(); - for (int j=0; jGetPoint(j, p); - - vtkIdType id = vNewPoints->InsertNextPoint(p); - container->GetPointIds()->InsertNextId(id); - } - weights->InsertValue(counter, fib->GetFiberWeight(i)); - vNewLines->InsertNextCell(container); - counter++; + vtkIdType id = vNewPoints->InsertNextPoint(p); + container->GetPointIds()->InsertNextId(id); } - - // initialize PolyData - vNewPolyData->SetPoints(vNewPoints); - vNewPolyData->SetLines(vNewLines); - - // initialize fiber bundle - mitk::FiberBundle::Pointer newFib = mitk::FiberBundle::New(vNewPolyData); - newFib->SetFiberWeights(weights); - return newFib; + weights->InsertValue(counter, fib->GetFiberWeight(i)); + vNewLines->InsertNextCell(container); + counter++; + } + + // initialize PolyData + vNewPolyData->SetPoints(vNewPoints); + vNewPolyData->SetLines(vNewLines); + + // initialize fiber bundle + mitk::FiberBundle::Pointer newFib = mitk::FiberBundle::New(vNewPolyData); + newFib->SetFiberWeights(weights); + return newFib; } // subtract two fiber bundles mitk::FiberBundle::Pointer mitk::FiberBundle::SubtractBundle(mitk::FiberBundle* fib) { - MITK_INFO << "Subtracting fibers"; + MITK_INFO << "Subtracting fibers"; - vtkSmartPointer vNewPolyData = vtkSmartPointer::New(); - vtkSmartPointer vNewLines = vtkSmartPointer::New(); - vtkSmartPointer vNewPoints = vtkSmartPointer::New(); + vtkSmartPointer vNewPolyData = vtkSmartPointer::New(); + vtkSmartPointer vNewLines = vtkSmartPointer::New(); + vtkSmartPointer vNewPoints = vtkSmartPointer::New(); - std::vector< std::vector< itk::Point > > points1; - for( int i=0; iGetCell(i); - int numPoints = cell->GetNumberOfPoints(); - vtkPoints* points = cell->GetPoints(); + std::vector< std::vector< itk::Point > > points1; + for( int i=0; iGetCell(i); + int numPoints = cell->GetNumberOfPoints(); + vtkPoints* points = cell->GetPoints(); - if (points==nullptr || numPoints<=0) - continue; + if (points==nullptr || numPoints<=0) + continue; - itk::Point start = GetItkPoint(points->GetPoint(0)); - itk::Point end = GetItkPoint(points->GetPoint(numPoints-1)); + itk::Point start = GetItkPoint(points->GetPoint(0)); + itk::Point end = GetItkPoint(points->GetPoint(numPoints-1)); - points1.push_back( {start, end} ); - } + points1.push_back( {start, end} ); + } - std::vector< std::vector< itk::Point > > points2; - for( int i=0; iGetNumFibers(); i++ ) - { - vtkCell* cell = fib->GetFiberPolyData()->GetCell(i); - int numPoints = cell->GetNumberOfPoints(); - vtkPoints* points = cell->GetPoints(); + std::vector< std::vector< itk::Point > > points2; + for( int i=0; iGetNumFibers(); i++ ) + { + vtkCell* cell = fib->GetFiberPolyData()->GetCell(i); + int numPoints = cell->GetNumberOfPoints(); + vtkPoints* points = cell->GetPoints(); - if (points==nullptr || numPoints<=0) - continue; + if (points==nullptr || numPoints<=0) + continue; - itk::Point start = GetItkPoint(points->GetPoint(0)); - itk::Point end = GetItkPoint(points->GetPoint(numPoints-1)); + itk::Point start = GetItkPoint(points->GetPoint(0)); + itk::Point end = GetItkPoint(points->GetPoint(numPoints-1)); - points2.push_back( {start, end} ); - } + points2.push_back( {start, end} ); + } - int progress = 0; - std::vector< int > ids; + int progress = 0; + std::vector< int > ids; #pragma omp parallel for - for (int i=0; iGetCell(i); - int numPoints = cell->GetNumberOfPoints(); - vtkPoints* points = cell->GetPoints(); + for( int i : ids ) + { + vtkCell* cell = m_FiberPolyData->GetCell(i); + int numPoints = cell->GetNumberOfPoints(); + vtkPoints* points = cell->GetPoints(); - if (points==nullptr || numPoints<=0) - continue; + if (points==nullptr || numPoints<=0) + continue; - vtkSmartPointer container = vtkSmartPointer::New(); - for( int j=0; jInsertNextPoint(points->GetPoint(j)); - container->GetPointIds()->InsertNextId(id); - } - vNewLines->InsertNextCell(container); + vtkSmartPointer container = vtkSmartPointer::New(); + for( int j=0; jInsertNextPoint(points->GetPoint(j)); + container->GetPointIds()->InsertNextId(id); } - if(vNewLines->GetNumberOfCells()==0) - return nullptr; - // initialize PolyData - vNewPolyData->SetPoints(vNewPoints); - vNewPolyData->SetLines(vNewLines); - - // initialize fiber bundle - return mitk::FiberBundle::New(vNewPolyData); + vNewLines->InsertNextCell(container); + } + if(vNewLines->GetNumberOfCells()==0) + return nullptr; + // initialize PolyData + vNewPolyData->SetPoints(vNewPoints); + vNewPolyData->SetLines(vNewLines); + + // initialize fiber bundle + return mitk::FiberBundle::New(vNewPolyData); } itk::Point mitk::FiberBundle::GetItkPoint(double point[3]) { - itk::Point itkPoint; - itkPoint[0] = point[0]; - itkPoint[1] = point[1]; - itkPoint[2] = point[2]; - return itkPoint; + itk::Point itkPoint; + itkPoint[0] = point[0]; + itkPoint[1] = point[1]; + itkPoint[2] = point[2]; + return itkPoint; } /* * set PolyData (additional flag to recompute fiber geometry, default = true) */ void mitk::FiberBundle::SetFiberPolyData(vtkSmartPointer fiberPD, bool updateGeometry) { - if (fiberPD == nullptr) - this->m_FiberPolyData = vtkSmartPointer::New(); - else - m_FiberPolyData->DeepCopy(fiberPD); + if (fiberPD == nullptr) + this->m_FiberPolyData = vtkSmartPointer::New(); + else + m_FiberPolyData->DeepCopy(fiberPD); - m_NumFibers = m_FiberPolyData->GetNumberOfLines(); + m_NumFibers = m_FiberPolyData->GetNumberOfLines(); - if (updateGeometry) - UpdateFiberGeometry(); - GenerateFiberIds(); - ColorFibersByOrientation(); + if (updateGeometry) + UpdateFiberGeometry(); + GenerateFiberIds(); + ColorFibersByOrientation(); } /* * return vtkPolyData */ vtkSmartPointer mitk::FiberBundle::GetFiberPolyData() const { - return m_FiberPolyData; + return m_FiberPolyData; } void mitk::FiberBundle::ColorFibersByOrientation() { - //===== FOR WRITING A TEST ======================== - // colorT size == tupelComponents * tupelElements - // compare color results - // to cover this code 100% also PolyData needed, where colorarray already exists - // + one fiber with exactly 1 point - // + one fiber with 0 points - //================================================= - - vtkPoints* extrPoints = nullptr; - extrPoints = m_FiberPolyData->GetPoints(); - int numOfPoints = 0; - if (extrPoints!=nullptr) - numOfPoints = extrPoints->GetNumberOfPoints(); - - //colors and alpha value for each single point, RGBA = 4 components - unsigned char rgba[4] = {0,0,0,0}; - int componentSize = 4; - m_FiberColors = vtkSmartPointer::New(); - m_FiberColors->Allocate(numOfPoints * componentSize); - m_FiberColors->SetNumberOfComponents(componentSize); - m_FiberColors->SetName("FIBER_COLORS"); - - int numOfFibers = m_FiberPolyData->GetNumberOfLines(); - if (numOfFibers < 1) - return; - - /* extract single fibers of fiberBundle */ - vtkCellArray* fiberList = m_FiberPolyData->GetLines(); - fiberList->InitTraversal(); - for (int fi=0; fiGetNextCell(pointsPerFiber, idList); - - /* single fiber checkpoints: is number of points valid */ - if (pointsPerFiber > 1) + //===== FOR WRITING A TEST ======================== + // colorT size == tupelComponents * tupelElements + // compare color results + // to cover this code 100% also PolyData needed, where colorarray already exists + // + one fiber with exactly 1 point + // + one fiber with 0 points + //================================================= + + vtkPoints* extrPoints = nullptr; + extrPoints = m_FiberPolyData->GetPoints(); + int numOfPoints = 0; + if (extrPoints!=nullptr) + numOfPoints = extrPoints->GetNumberOfPoints(); + + //colors and alpha value for each single point, RGBA = 4 components + unsigned char rgba[4] = {0,0,0,0}; + int componentSize = 4; + m_FiberColors = vtkSmartPointer::New(); + m_FiberColors->Allocate(numOfPoints * componentSize); + m_FiberColors->SetNumberOfComponents(componentSize); + m_FiberColors->SetName("FIBER_COLORS"); + + int numOfFibers = m_FiberPolyData->GetNumberOfLines(); + if (numOfFibers < 1) + return; + + /* extract single fibers of fiberBundle */ + vtkCellArray* fiberList = m_FiberPolyData->GetLines(); + fiberList->InitTraversal(); + for (int fi=0; fiGetNextCell(pointsPerFiber, idList); + + /* single fiber checkpoints: is number of points valid */ + if (pointsPerFiber > 1) + { + /* operate on points of single fiber */ + for (int i=0; i 0) { - /* operate on points of single fiber */ - for (int i=0; i 0) - { - /* The color value of the current point is influenced by the previous point and next point. */ - vnl_vector_fixed< double, 3 > currentPntvtk(extrPoints->GetPoint(idList[i])[0], extrPoints->GetPoint(idList[i])[1],extrPoints->GetPoint(idList[i])[2]); - vnl_vector_fixed< double, 3 > nextPntvtk(extrPoints->GetPoint(idList[i+1])[0], extrPoints->GetPoint(idList[i+1])[1], extrPoints->GetPoint(idList[i+1])[2]); - vnl_vector_fixed< double, 3 > prevPntvtk(extrPoints->GetPoint(idList[i-1])[0], extrPoints->GetPoint(idList[i-1])[1], extrPoints->GetPoint(idList[i-1])[2]); - - vnl_vector_fixed< double, 3 > diff1; - diff1 = currentPntvtk - nextPntvtk; - - vnl_vector_fixed< double, 3 > diff2; - diff2 = currentPntvtk - prevPntvtk; - - vnl_vector_fixed< double, 3 > diff; - diff = (diff1 - diff2) / 2.0; - diff.normalize(); - - rgba[0] = (unsigned char) (255.0 * std::fabs(diff[0])); - rgba[1] = (unsigned char) (255.0 * std::fabs(diff[1])); - rgba[2] = (unsigned char) (255.0 * std::fabs(diff[2])); - rgba[3] = (unsigned char) (255.0); - } - else if (i==0) - { - /* First point has no previous point, therefore only diff1 is taken */ - - vnl_vector_fixed< double, 3 > currentPntvtk(extrPoints->GetPoint(idList[i])[0], extrPoints->GetPoint(idList[i])[1],extrPoints->GetPoint(idList[i])[2]); - vnl_vector_fixed< double, 3 > nextPntvtk(extrPoints->GetPoint(idList[i+1])[0], extrPoints->GetPoint(idList[i+1])[1], extrPoints->GetPoint(idList[i+1])[2]); - - vnl_vector_fixed< double, 3 > diff1; - diff1 = currentPntvtk - nextPntvtk; - diff1.normalize(); - - rgba[0] = (unsigned char) (255.0 * std::fabs(diff1[0])); - rgba[1] = (unsigned char) (255.0 * std::fabs(diff1[1])); - rgba[2] = (unsigned char) (255.0 * std::fabs(diff1[2])); - rgba[3] = (unsigned char) (255.0); - } - else if (i==pointsPerFiber-1) - { - /* Last point has no next point, therefore only diff2 is taken */ - vnl_vector_fixed< double, 3 > currentPntvtk(extrPoints->GetPoint(idList[i])[0], extrPoints->GetPoint(idList[i])[1],extrPoints->GetPoint(idList[i])[2]); - vnl_vector_fixed< double, 3 > prevPntvtk(extrPoints->GetPoint(idList[i-1])[0], extrPoints->GetPoint(idList[i-1])[1], extrPoints->GetPoint(idList[i-1])[2]); - - vnl_vector_fixed< double, 3 > diff2; - diff2 = currentPntvtk - prevPntvtk; - diff2.normalize(); - - rgba[0] = (unsigned char) (255.0 * std::fabs(diff2[0])); - rgba[1] = (unsigned char) (255.0 * std::fabs(diff2[1])); - rgba[2] = (unsigned char) (255.0 * std::fabs(diff2[2])); - rgba[3] = (unsigned char) (255.0); - } - m_FiberColors->InsertTupleValue(idList[i], rgba); - } + /* The color value of the current point is influenced by the previous point and next point. */ + vnl_vector_fixed< double, 3 > currentPntvtk(extrPoints->GetPoint(idList[i])[0], extrPoints->GetPoint(idList[i])[1],extrPoints->GetPoint(idList[i])[2]); + vnl_vector_fixed< double, 3 > nextPntvtk(extrPoints->GetPoint(idList[i+1])[0], extrPoints->GetPoint(idList[i+1])[1], extrPoints->GetPoint(idList[i+1])[2]); + vnl_vector_fixed< double, 3 > prevPntvtk(extrPoints->GetPoint(idList[i-1])[0], extrPoints->GetPoint(idList[i-1])[1], extrPoints->GetPoint(idList[i-1])[2]); + + vnl_vector_fixed< double, 3 > diff1; + diff1 = currentPntvtk - nextPntvtk; + + vnl_vector_fixed< double, 3 > diff2; + diff2 = currentPntvtk - prevPntvtk; + + vnl_vector_fixed< double, 3 > diff; + diff = (diff1 - diff2) / 2.0; + diff.normalize(); + + rgba[0] = (unsigned char) (255.0 * std::fabs(diff[0])); + rgba[1] = (unsigned char) (255.0 * std::fabs(diff[1])); + rgba[2] = (unsigned char) (255.0 * std::fabs(diff[2])); + rgba[3] = (unsigned char) (255.0); } - else if (pointsPerFiber == 1) + else if (i==0) { - /* a single point does not define a fiber (use vertex mechanisms instead */ - continue; + /* First point has no previous point, therefore only diff1 is taken */ + + vnl_vector_fixed< double, 3 > currentPntvtk(extrPoints->GetPoint(idList[i])[0], extrPoints->GetPoint(idList[i])[1],extrPoints->GetPoint(idList[i])[2]); + vnl_vector_fixed< double, 3 > nextPntvtk(extrPoints->GetPoint(idList[i+1])[0], extrPoints->GetPoint(idList[i+1])[1], extrPoints->GetPoint(idList[i+1])[2]); + + vnl_vector_fixed< double, 3 > diff1; + diff1 = currentPntvtk - nextPntvtk; + diff1.normalize(); + + rgba[0] = (unsigned char) (255.0 * std::fabs(diff1[0])); + rgba[1] = (unsigned char) (255.0 * std::fabs(diff1[1])); + rgba[2] = (unsigned char) (255.0 * std::fabs(diff1[2])); + rgba[3] = (unsigned char) (255.0); } - else + else if (i==pointsPerFiber-1) { - MITK_DEBUG << "Fiber with 0 points detected... please check your tractography algorithm!" ; - continue; + /* Last point has no next point, therefore only diff2 is taken */ + vnl_vector_fixed< double, 3 > currentPntvtk(extrPoints->GetPoint(idList[i])[0], extrPoints->GetPoint(idList[i])[1],extrPoints->GetPoint(idList[i])[2]); + vnl_vector_fixed< double, 3 > prevPntvtk(extrPoints->GetPoint(idList[i-1])[0], extrPoints->GetPoint(idList[i-1])[1], extrPoints->GetPoint(idList[i-1])[2]); + + vnl_vector_fixed< double, 3 > diff2; + diff2 = currentPntvtk - prevPntvtk; + diff2.normalize(); + + rgba[0] = (unsigned char) (255.0 * std::fabs(diff2[0])); + rgba[1] = (unsigned char) (255.0 * std::fabs(diff2[1])); + rgba[2] = (unsigned char) (255.0 * std::fabs(diff2[2])); + rgba[3] = (unsigned char) (255.0); } + m_FiberColors->InsertTupleValue(idList[i], rgba); + } + } + else if (pointsPerFiber == 1) + { + /* a single point does not define a fiber (use vertex mechanisms instead */ + continue; + } + else + { + MITK_DEBUG << "Fiber with 0 points detected... please check your tractography algorithm!" ; + continue; } - m_UpdateTime3D.Modified(); - m_UpdateTime2D.Modified(); + } + m_UpdateTime3D.Modified(); + m_UpdateTime2D.Modified(); } void mitk::FiberBundle::ColorFibersByCurvature(bool minMaxNorm) { - double window = 5; - - //colors and alpha value for each single point, RGBA = 4 components - unsigned char rgba[4] = {0,0,0,0}; - int componentSize = 4; - m_FiberColors = vtkSmartPointer::New(); - m_FiberColors->Allocate(m_FiberPolyData->GetNumberOfPoints() * componentSize); - m_FiberColors->SetNumberOfComponents(componentSize); - m_FiberColors->SetName("FIBER_COLORS"); - - mitk::LookupTable::Pointer mitkLookup = mitk::LookupTable::New(); - vtkSmartPointer lookupTable = vtkSmartPointer::New(); - lookupTable->SetTableRange(0.0, 0.8); - lookupTable->Build(); - mitkLookup->SetVtkLookupTable(lookupTable); - mitkLookup->SetType(mitk::LookupTable::JET); - - vector< double > values; - double min = 1; - double max = 0; - MITK_INFO << "Coloring fibers by curvature"; - boost::progress_display disp(m_FiberPolyData->GetNumberOfCells()); - for (int i=0; iGetNumberOfCells(); i++) + double window = 5; + + //colors and alpha value for each single point, RGBA = 4 components + unsigned char rgba[4] = {0,0,0,0}; + int componentSize = 4; + m_FiberColors = vtkSmartPointer::New(); + m_FiberColors->Allocate(m_FiberPolyData->GetNumberOfPoints() * componentSize); + m_FiberColors->SetNumberOfComponents(componentSize); + m_FiberColors->SetName("FIBER_COLORS"); + + mitk::LookupTable::Pointer mitkLookup = mitk::LookupTable::New(); + vtkSmartPointer lookupTable = vtkSmartPointer::New(); + lookupTable->SetTableRange(0.0, 0.8); + lookupTable->Build(); + mitkLookup->SetVtkLookupTable(lookupTable); + mitkLookup->SetType(mitk::LookupTable::JET); + + vector< double > values; + double min = 1; + double max = 0; + MITK_INFO << "Coloring fibers by curvature"; + boost::progress_display disp(m_FiberPolyData->GetNumberOfCells()); + for (int i=0; iGetNumberOfCells(); i++) + { + ++disp; + vtkCell* cell = m_FiberPolyData->GetCell(i); + int numPoints = cell->GetNumberOfPoints(); + vtkPoints* points = cell->GetPoints(); + + // calculate curvatures + for (int j=0; jGetCell(i); - int numPoints = cell->GetNumberOfPoints(); - vtkPoints* points = cell->GetPoints(); - - // calculate curvatures - for (int j=0; j > vectors; - vnl_vector_fixed< float, 3 > meanV; meanV.fill(0.0); - while(dist1) - { - double p1[3]; - points->GetPoint(c-1, p1); - double p2[3]; - points->GetPoint(c, p2); - - vnl_vector_fixed< float, 3 > v; - v[0] = p2[0]-p1[0]; - v[1] = p2[1]-p1[1]; - v[2] = p2[2]-p1[2]; - dist += v.magnitude(); - v.normalize(); - vectors.push_back(v); - if (c==j) - meanV += v; - c--; - } - c = j; - dist = 0; - while(distGetPoint(c, p1); - double p2[3]; - points->GetPoint(c+1, p2); - - vnl_vector_fixed< float, 3 > v; - v[0] = p2[0]-p1[0]; - v[1] = p2[1]-p1[1]; - v[2] = p2[2]-p1[2]; - dist += v.magnitude(); - v.normalize(); - vectors.push_back(v); - if (c==j) - meanV += v; - c++; - } - meanV.normalize(); - - double dev = 0; - for (unsigned int c=0; c1.0) - angle = 1.0; - if (angle<-1.0) - angle = -1.0; - dev += acos(angle)*180/M_PI; - } - if (vectors.size()>0) - dev /= vectors.size(); - - dev = 1.0-dev/180.0; - values.push_back(dev); - if (devmax) - max = dev; - } + double dist = 0; + int c = j; + std::vector< vnl_vector_fixed< float, 3 > > vectors; + vnl_vector_fixed< float, 3 > meanV; meanV.fill(0.0); + while(dist1) + { + double p1[3]; + points->GetPoint(c-1, p1); + double p2[3]; + points->GetPoint(c, p2); + + vnl_vector_fixed< float, 3 > v; + v[0] = p2[0]-p1[0]; + v[1] = p2[1]-p1[1]; + v[2] = p2[2]-p1[2]; + dist += v.magnitude(); + v.normalize(); + vectors.push_back(v); + if (c==j) + meanV += v; + c--; + } + c = j; + dist = 0; + while(distGetPoint(c, p1); + double p2[3]; + points->GetPoint(c+1, p2); + + vnl_vector_fixed< float, 3 > v; + v[0] = p2[0]-p1[0]; + v[1] = p2[1]-p1[1]; + v[2] = p2[2]-p1[2]; + dist += v.magnitude(); + v.normalize(); + vectors.push_back(v); + if (c==j) + meanV += v; + c++; + } + meanV.normalize(); + + double dev = 0; + for (unsigned int c=0; c1.0) + angle = 1.0; + if (angle<-1.0) + angle = -1.0; + dev += acos(angle)*180/M_PI; + } + if (vectors.size()>0) + dev /= vectors.size(); + + dev = 1.0-dev/180.0; + values.push_back(dev); + if (devmax) + max = dev; } - unsigned int count = 0; - for (int i=0; iGetNumberOfCells(); i++) + } + unsigned int count = 0; + for (int i=0; iGetNumberOfCells(); i++) + { + vtkCell* cell = m_FiberPolyData->GetCell(i); + int numPoints = cell->GetNumberOfPoints(); + for (int j=0; jGetCell(i); - int numPoints = cell->GetNumberOfPoints(); - for (int j=0; jGetColor(dev, color); - - rgba[0] = (unsigned char) (255.0 * color[0]); - rgba[1] = (unsigned char) (255.0 * color[1]); - rgba[2] = (unsigned char) (255.0 * color[2]); - rgba[3] = (unsigned char) (255.0); - m_FiberColors->InsertTupleValue(cell->GetPointId(j), rgba); - count++; - } + double color[3]; + double dev = values.at(count); + if (minMaxNorm) + dev = (dev-min)/(max-min); + // double dev = values.at(count)*values.at(count); + lookupTable->GetColor(dev, color); + + rgba[0] = (unsigned char) (255.0 * color[0]); + rgba[1] = (unsigned char) (255.0 * color[1]); + rgba[2] = (unsigned char) (255.0 * color[2]); + rgba[3] = (unsigned char) (255.0); + m_FiberColors->InsertTupleValue(cell->GetPointId(j), rgba); + count++; } - m_UpdateTime3D.Modified(); - m_UpdateTime2D.Modified(); + } + m_UpdateTime3D.Modified(); + m_UpdateTime2D.Modified(); } void mitk::FiberBundle::SetFiberOpacity(vtkDoubleArray* FAValArray) { - for(long i=0; iGetNumberOfTuples(); i++) - { - double faValue = FAValArray->GetValue(i); - faValue = faValue * 255.0; - m_FiberColors->SetComponent(i,3, (unsigned char) faValue ); - } - m_UpdateTime3D.Modified(); - m_UpdateTime2D.Modified(); + for(long i=0; iGetNumberOfTuples(); i++) + { + double faValue = FAValArray->GetValue(i); + faValue = faValue * 255.0; + m_FiberColors->SetComponent(i,3, (unsigned char) faValue ); + } + m_UpdateTime3D.Modified(); + m_UpdateTime2D.Modified(); } void mitk::FiberBundle::ResetFiberOpacity() { - for(long i=0; iGetNumberOfTuples(); i++) - m_FiberColors->SetComponent(i,3, 255.0 ); - m_UpdateTime3D.Modified(); - m_UpdateTime2D.Modified(); + for(long i=0; iGetNumberOfTuples(); i++) + m_FiberColors->SetComponent(i,3, 255.0 ); + m_UpdateTime3D.Modified(); + m_UpdateTime2D.Modified(); } void mitk::FiberBundle::ColorFibersByScalarMap(mitk::Image::Pointer FAimage, bool opacity) { - mitkPixelTypeMultiplex2( ColorFibersByScalarMap, FAimage->GetPixelType(), FAimage, opacity ); - m_UpdateTime3D.Modified(); - m_UpdateTime2D.Modified(); + mitkPixelTypeMultiplex2( ColorFibersByScalarMap, FAimage->GetPixelType(), FAimage, opacity ); + m_UpdateTime3D.Modified(); + m_UpdateTime2D.Modified(); } template void mitk::FiberBundle::ColorFibersByScalarMap(const mitk::PixelType, mitk::Image::Pointer image, bool opacity) { - m_FiberColors = vtkSmartPointer::New(); - m_FiberColors->Allocate(m_FiberPolyData->GetNumberOfPoints() * 4); - m_FiberColors->SetNumberOfComponents(4); - m_FiberColors->SetName("FIBER_COLORS"); + m_FiberColors = vtkSmartPointer::New(); + m_FiberColors->Allocate(m_FiberPolyData->GetNumberOfPoints() * 4); + m_FiberColors->SetNumberOfComponents(4); + m_FiberColors->SetName("FIBER_COLORS"); + + mitk::ImagePixelReadAccessor readimage(image, image->GetVolumeData(0)); + + unsigned char rgba[4] = {0,0,0,0}; + vtkPoints* pointSet = m_FiberPolyData->GetPoints(); + + mitk::LookupTable::Pointer mitkLookup = mitk::LookupTable::New(); + vtkSmartPointer lookupTable = vtkSmartPointer::New(); + lookupTable->SetTableRange(0.0, 0.8); + lookupTable->Build(); + mitkLookup->SetVtkLookupTable(lookupTable); + mitkLookup->SetType(mitk::LookupTable::JET); + + for(long i=0; iGetNumberOfPoints(); ++i) + { + Point3D px; + px[0] = pointSet->GetPoint(i)[0]; + px[1] = pointSet->GetPoint(i)[1]; + px[2] = pointSet->GetPoint(i)[2]; + double pixelValue = readimage.GetPixelByWorldCoordinates(px); + + double color[3]; + lookupTable->GetColor(1-pixelValue, color); + + rgba[0] = (unsigned char) (255.0 * color[0]); + rgba[1] = (unsigned char) (255.0 * color[1]); + rgba[2] = (unsigned char) (255.0 * color[2]); + if (opacity) + rgba[3] = (unsigned char) (255.0 * pixelValue); + else + rgba[3] = (unsigned char) (255.0); + m_FiberColors->InsertTupleValue(i, rgba); + } + m_UpdateTime3D.Modified(); + m_UpdateTime2D.Modified(); +} - mitk::ImagePixelReadAccessor readimage(image, image->GetVolumeData(0)); - unsigned char rgba[4] = {0,0,0,0}; - vtkPoints* pointSet = m_FiberPolyData->GetPoints(); +void mitk::FiberBundle::ColorFibersByFiberWeights() +{ + m_FiberColors = vtkSmartPointer::New(); + m_FiberColors->Allocate(m_FiberPolyData->GetNumberOfPoints() * 4); + m_FiberColors->SetNumberOfComponents(4); + m_FiberColors->SetName("FIBER_COLORS"); + + mitk::LookupTable::Pointer mitkLookup = mitk::LookupTable::New(); + vtkSmartPointer lookupTable = vtkSmartPointer::New(); + lookupTable->SetTableRange(0.0, 0.8); + lookupTable->Build(); + mitkLookup->SetVtkLookupTable(lookupTable); + mitkLookup->SetType(mitk::LookupTable::JET); + + unsigned char rgba[4] = {0,0,0,0}; + unsigned int counter = 0; + + float max = -999999; + float min = 999999; + for (int i=0; iGetFiberWeight(i); + if (weight>max) + max = weight; + if (weightGetCell(i); + int numPoints = cell->GetNumberOfPoints(); + double weight = this->GetFiberWeight(i); + + for (int j=0; jGetColor(1-(weight-min)/(max-min), color); - mitk::LookupTable::Pointer mitkLookup = mitk::LookupTable::New(); - vtkSmartPointer lookupTable = vtkSmartPointer::New(); - lookupTable->SetTableRange(0.0, 0.8); - lookupTable->Build(); - mitkLookup->SetVtkLookupTable(lookupTable); - mitkLookup->SetType(mitk::LookupTable::JET); + rgba[0] = (unsigned char) (255.0 * color[0]); + rgba[1] = (unsigned char) (255.0 * color[1]); + rgba[2] = (unsigned char) (255.0 * color[2]); + rgba[3] = (unsigned char) (255.0); - for(long i=0; iGetNumberOfPoints(); ++i) - { - Point3D px; - px[0] = pointSet->GetPoint(i)[0]; - px[1] = pointSet->GetPoint(i)[1]; - px[2] = pointSet->GetPoint(i)[2]; - double pixelValue = readimage.GetPixelByWorldCoordinates(px); - - double color[3]; - lookupTable->GetColor(1-pixelValue, color); - - rgba[0] = (unsigned char) (255.0 * color[0]); - rgba[1] = (unsigned char) (255.0 * color[1]); - rgba[2] = (unsigned char) (255.0 * color[2]); - if (opacity) - rgba[3] = (unsigned char) (255.0 * pixelValue); - else - rgba[3] = (unsigned char) (255.0); - m_FiberColors->InsertTupleValue(i, rgba); + m_FiberColors->InsertTupleValue(counter, rgba); + counter++; } - m_UpdateTime3D.Modified(); - m_UpdateTime2D.Modified(); + } + + m_UpdateTime3D.Modified(); + m_UpdateTime2D.Modified(); } void mitk::FiberBundle::SetFiberColors(float r, float g, float b, float alpha) { - m_FiberColors = vtkSmartPointer::New(); - m_FiberColors->Allocate(m_FiberPolyData->GetNumberOfPoints() * 4); - m_FiberColors->SetNumberOfComponents(4); - m_FiberColors->SetName("FIBER_COLORS"); - - unsigned char rgba[4] = {0,0,0,0}; - for(long i=0; iGetNumberOfPoints(); ++i) - { - rgba[0] = (unsigned char) r; - rgba[1] = (unsigned char) g; - rgba[2] = (unsigned char) b; - rgba[3] = (unsigned char) alpha; - m_FiberColors->InsertTupleValue(i, rgba); - } - m_UpdateTime3D.Modified(); - m_UpdateTime2D.Modified(); + m_FiberColors = vtkSmartPointer::New(); + m_FiberColors->Allocate(m_FiberPolyData->GetNumberOfPoints() * 4); + m_FiberColors->SetNumberOfComponents(4); + m_FiberColors->SetName("FIBER_COLORS"); + + unsigned char rgba[4] = {0,0,0,0}; + for(long i=0; iGetNumberOfPoints(); ++i) + { + rgba[0] = (unsigned char) r; + rgba[1] = (unsigned char) g; + rgba[2] = (unsigned char) b; + rgba[3] = (unsigned char) alpha; + m_FiberColors->InsertTupleValue(i, rgba); + } + m_UpdateTime3D.Modified(); + m_UpdateTime2D.Modified(); } void mitk::FiberBundle::GenerateFiberIds() { - if (m_FiberPolyData == nullptr) - return; + if (m_FiberPolyData == nullptr) + return; - vtkSmartPointer idFiberFilter = vtkSmartPointer::New(); - idFiberFilter->SetInputData(m_FiberPolyData); - idFiberFilter->CellIdsOn(); - // idFiberFilter->PointIdsOn(); // point id's are not needed - idFiberFilter->SetIdsArrayName(FIBER_ID_ARRAY); - idFiberFilter->FieldDataOn(); - idFiberFilter->Update(); + vtkSmartPointer idFiberFilter = vtkSmartPointer::New(); + idFiberFilter->SetInputData(m_FiberPolyData); + idFiberFilter->CellIdsOn(); + // idFiberFilter->PointIdsOn(); // point id's are not needed + idFiberFilter->SetIdsArrayName(FIBER_ID_ARRAY); + idFiberFilter->FieldDataOn(); + idFiberFilter->Update(); - m_FiberIdDataSet = idFiberFilter->GetOutput(); + m_FiberIdDataSet = idFiberFilter->GetOutput(); } mitk::FiberBundle::Pointer mitk::FiberBundle::ExtractFiberSubset(ItkUcharImgType* mask, bool anyPoint, bool invert, bool bothEnds, float fraction) { - vtkSmartPointer PolyData = m_FiberPolyData; - if (anyPoint) - { - float minSpacing = 1; - if(mask->GetSpacing()[0]GetSpacing()[1] && mask->GetSpacing()[0]GetSpacing()[2]) - minSpacing = mask->GetSpacing()[0]; - else if (mask->GetSpacing()[1] < mask->GetSpacing()[2]) - minSpacing = mask->GetSpacing()[1]; - else - minSpacing = mask->GetSpacing()[2]; + vtkSmartPointer PolyData = m_FiberPolyData; + if (anyPoint) + { + float minSpacing = 1; + if(mask->GetSpacing()[0]GetSpacing()[1] && mask->GetSpacing()[0]GetSpacing()[2]) + minSpacing = mask->GetSpacing()[0]; + else if (mask->GetSpacing()[1] < mask->GetSpacing()[2]) + minSpacing = mask->GetSpacing()[1]; + else + minSpacing = mask->GetSpacing()[2]; - mitk::FiberBundle::Pointer fibCopy = this->GetDeepCopy(); - fibCopy->ResampleLinear(minSpacing/5); - PolyData = fibCopy->GetFiberPolyData(); - } - vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); - vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); + mitk::FiberBundle::Pointer fibCopy = this->GetDeepCopy(); + fibCopy->ResampleLinear(minSpacing/5); + PolyData = fibCopy->GetFiberPolyData(); + } + vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); + vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); - MITK_INFO << "Extracting fibers"; - boost::progress_display disp(m_NumFibers); - for (int i=0; iGetCell(i); - int numPoints = cell->GetNumberOfPoints(); - vtkPoints* points = cell->GetPoints(); + vtkCell* cell = PolyData->GetCell(i); + int numPoints = cell->GetNumberOfPoints(); + vtkPoints* points = cell->GetPoints(); - vtkCell* cellOriginal = m_FiberPolyData->GetCell(i); - int numPointsOriginal = cellOriginal->GetNumberOfPoints(); - vtkPoints* pointsOriginal = cellOriginal->GetPoints(); + vtkCell* cellOriginal = m_FiberPolyData->GetCell(i); + int numPointsOriginal = cellOriginal->GetNumberOfPoints(); + vtkPoints* pointsOriginal = cellOriginal->GetPoints(); - vtkSmartPointer container = vtkSmartPointer::New(); + vtkSmartPointer container = vtkSmartPointer::New(); - if (numPoints>1 && numPointsOriginal) + if (numPoints>1 && numPointsOriginal) + { + if (anyPoint) + { + int inside = 0; + int outside = 0; + if (!invert) { - if (anyPoint) + for (int j=0; jGetPoint(j); + + itk::Point itkP; + itkP[0] = p[0]; itkP[1] = p[1]; itkP[2] = p[2]; + itk::Index<3> idx; + mask->TransformPhysicalPointToIndex(itkP, idx); + + if ( mask->GetLargestPossibleRegion().IsInside(idx) && mask->GetPixel(idx) != 0 ) { - int inside = 0; - int outside = 0; - if (!invert) - { - for (int j=0; jGetPoint(j); - - itk::Point itkP; - itkP[0] = p[0]; itkP[1] = p[1]; itkP[2] = p[2]; - itk::Index<3> idx; - mask->TransformPhysicalPointToIndex(itkP, idx); - - if ( mask->GetLargestPossibleRegion().IsInside(idx) && mask->GetPixel(idx) != 0 ) - { - inside++; - if (fraction==0) - break; - } - else - outside++; - } - - float current_fraction = 0.0; - if (inside+outside>0) - current_fraction = (float)inside/(inside+outside); - - if (current_fraction>fraction) - { - for (int k=0; kGetPoint(k); - vtkIdType id = vtkNewPoints->InsertNextPoint(p); - container->GetPointIds()->InsertNextId(id); - } - } - } - else - { - bool includeFiber = true; - for (int j=0; jGetPoint(j); - - itk::Point itkP; - itkP[0] = p[0]; itkP[1] = p[1]; itkP[2] = p[2]; - itk::Index<3> idx; - mask->TransformPhysicalPointToIndex(itkP, idx); - - if ( mask->GetPixel(idx) != 0 && mask->GetLargestPossibleRegion().IsInside(idx) ) - { - inside++; - includeFiber = false; - break; - } - else - outside++; - } - if (includeFiber) - { - - for (int k=0; kGetPoint(k); - vtkIdType id = vtkNewPoints->InsertNextPoint(p); - container->GetPointIds()->InsertNextId(id); - } - } - } + inside++; + if (fraction==0) + break; } else + outside++; + } + + float current_fraction = 0.0; + if (inside+outside>0) + current_fraction = (float)inside/(inside+outside); + + if (current_fraction>fraction) + { + for (int k=0; kGetPoint(0); - itk::Point itkStart; - itkStart[0] = start[0]; itkStart[1] = start[1]; itkStart[2] = start[2]; - itk::Index<3> idxStart; - mask->TransformPhysicalPointToIndex(itkStart, idxStart); - - double* end = pointsOriginal->GetPoint(numPointsOriginal-1); - itk::Point itkEnd; - itkEnd[0] = end[0]; itkEnd[1] = end[1]; itkEnd[2] = end[2]; - itk::Index<3> idxEnd; - mask->TransformPhysicalPointToIndex(itkEnd, idxEnd); - - if (invert) - { - if (bothEnds) - { - if ( mask->GetPixel(idxStart) == 0 && mask->GetPixel(idxEnd) == 0 ) - { - for (int j=0; jGetPoint(j); - vtkIdType id = vtkNewPoints->InsertNextPoint(p); - container->GetPointIds()->InsertNextId(id); - } - } - } - else if ( mask->GetPixel(idxStart) == 0 || mask->GetPixel(idxEnd) == 0 ) - { - for (int j=0; jGetPoint(j); - vtkIdType id = vtkNewPoints->InsertNextPoint(p); - container->GetPointIds()->InsertNextId(id); - } - } - } - else - { - if (bothEnds) - { - if ( mask->GetPixel(idxStart) != 0 && mask->GetPixel(idxEnd) != 0 && mask->GetLargestPossibleRegion().IsInside(idxStart) && mask->GetLargestPossibleRegion().IsInside(idxEnd) ) - { - for (int j=0; jGetPoint(j); - vtkIdType id = vtkNewPoints->InsertNextPoint(p); - container->GetPointIds()->InsertNextId(id); - } - } - } - else if ( (mask->GetPixel(idxStart) != 0 && mask->GetLargestPossibleRegion().IsInside(idxStart)) || (mask->GetPixel(idxEnd) != 0 && mask->GetLargestPossibleRegion().IsInside(idxEnd)) ) - { - for (int j=0; jGetPoint(j); - vtkIdType id = vtkNewPoints->InsertNextPoint(p); - container->GetPointIds()->InsertNextId(id); - } - } - } + double* p = points->GetPoint(k); + vtkIdType id = vtkNewPoints->InsertNextPoint(p); + container->GetPointIds()->InsertNextId(id); } + } } + else + { + bool includeFiber = true; + for (int j=0; jGetPoint(j); - vtkNewCells->InsertNextCell(container); + itk::Point itkP; + itkP[0] = p[0]; itkP[1] = p[1]; itkP[2] = p[2]; + itk::Index<3> idx; + mask->TransformPhysicalPointToIndex(itkP, idx); + + if ( mask->GetPixel(idx) != 0 && mask->GetLargestPossibleRegion().IsInside(idx) ) + { + inside++; + includeFiber = false; + break; + } + else + outside++; + } + if (includeFiber) + { + + for (int k=0; kGetPoint(k); + vtkIdType id = vtkNewPoints->InsertNextPoint(p); + container->GetPointIds()->InsertNextId(id); + } + } + } + } + else + { + double* start = pointsOriginal->GetPoint(0); + itk::Point itkStart; + itkStart[0] = start[0]; itkStart[1] = start[1]; itkStart[2] = start[2]; + itk::Index<3> idxStart; + mask->TransformPhysicalPointToIndex(itkStart, idxStart); + + double* end = pointsOriginal->GetPoint(numPointsOriginal-1); + itk::Point itkEnd; + itkEnd[0] = end[0]; itkEnd[1] = end[1]; itkEnd[2] = end[2]; + itk::Index<3> idxEnd; + mask->TransformPhysicalPointToIndex(itkEnd, idxEnd); + + if (invert) + { + if (bothEnds) + { + if ( mask->GetPixel(idxStart) == 0 && mask->GetPixel(idxEnd) == 0 ) + { + for (int j=0; jGetPoint(j); + vtkIdType id = vtkNewPoints->InsertNextPoint(p); + container->GetPointIds()->InsertNextId(id); + } + } + } + else if ( mask->GetPixel(idxStart) == 0 || mask->GetPixel(idxEnd) == 0 ) + { + for (int j=0; jGetPoint(j); + vtkIdType id = vtkNewPoints->InsertNextPoint(p); + container->GetPointIds()->InsertNextId(id); + } + } + } + else + { + if (bothEnds) + { + if ( mask->GetPixel(idxStart) != 0 && mask->GetPixel(idxEnd) != 0 && mask->GetLargestPossibleRegion().IsInside(idxStart) && mask->GetLargestPossibleRegion().IsInside(idxEnd) ) + { + for (int j=0; jGetPoint(j); + vtkIdType id = vtkNewPoints->InsertNextPoint(p); + container->GetPointIds()->InsertNextId(id); + } + } + } + else if ( (mask->GetPixel(idxStart) != 0 && mask->GetLargestPossibleRegion().IsInside(idxStart)) || (mask->GetPixel(idxEnd) != 0 && mask->GetLargestPossibleRegion().IsInside(idxEnd)) ) + { + for (int j=0; jGetPoint(j); + vtkIdType id = vtkNewPoints->InsertNextPoint(p); + container->GetPointIds()->InsertNextId(id); + } + } + } + } } - if (vtkNewCells->GetNumberOfCells()<=0) - return nullptr; + vtkNewCells->InsertNextCell(container); + } - vtkSmartPointer newPolyData = vtkSmartPointer::New(); - newPolyData->SetPoints(vtkNewPoints); - newPolyData->SetLines(vtkNewCells); - return mitk::FiberBundle::New(newPolyData); + if (vtkNewCells->GetNumberOfCells()<=0) + return nullptr; + + vtkSmartPointer newPolyData = vtkSmartPointer::New(); + newPolyData->SetPoints(vtkNewPoints); + newPolyData->SetLines(vtkNewCells); + return mitk::FiberBundle::New(newPolyData); } mitk::FiberBundle::Pointer mitk::FiberBundle::RemoveFibersOutside(ItkUcharImgType* mask, bool invert) { - float minSpacing = 1; - if(mask->GetSpacing()[0]GetSpacing()[1] && mask->GetSpacing()[0]GetSpacing()[2]) - minSpacing = mask->GetSpacing()[0]; - else if (mask->GetSpacing()[1] < mask->GetSpacing()[2]) - minSpacing = mask->GetSpacing()[1]; - else - minSpacing = mask->GetSpacing()[2]; - - mitk::FiberBundle::Pointer fibCopy = this->GetDeepCopy(); - fibCopy->ResampleLinear(minSpacing/10); - vtkSmartPointer PolyData =fibCopy->GetFiberPolyData(); - - vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); - vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); - - MITK_INFO << "Cutting fibers"; - boost::progress_display disp(m_NumFibers); - for (int i=0; iGetSpacing()[0]GetSpacing()[1] && mask->GetSpacing()[0]GetSpacing()[2]) + minSpacing = mask->GetSpacing()[0]; + else if (mask->GetSpacing()[1] < mask->GetSpacing()[2]) + minSpacing = mask->GetSpacing()[1]; + else + minSpacing = mask->GetSpacing()[2]; + + mitk::FiberBundle::Pointer fibCopy = this->GetDeepCopy(); + fibCopy->ResampleLinear(minSpacing/10); + vtkSmartPointer PolyData =fibCopy->GetFiberPolyData(); + + vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); + vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); + + MITK_INFO << "Cutting fibers"; + boost::progress_display disp(m_NumFibers); + for (int i=0; iGetCell(i); + int numPoints = cell->GetNumberOfPoints(); + vtkPoints* points = cell->GetPoints(); + + vtkSmartPointer container = vtkSmartPointer::New(); + if (numPoints>1) { - ++disp; + int newNumPoints = 0; + for (int j=0; jGetPoint(j); - vtkCell* cell = PolyData->GetCell(i); - int numPoints = cell->GetNumberOfPoints(); - vtkPoints* points = cell->GetPoints(); + itk::Point itkP; + itkP[0] = p[0]; itkP[1] = p[1]; itkP[2] = p[2]; + itk::Index<3> idx; + mask->TransformPhysicalPointToIndex(itkP, idx); - vtkSmartPointer container = vtkSmartPointer::New(); - if (numPoints>1) + if ( mask->GetPixel(idx) != 0 && mask->GetLargestPossibleRegion().IsInside(idx) && !invert ) { - int newNumPoints = 0; - for (int j=0; jGetPoint(j); - - itk::Point itkP; - itkP[0] = p[0]; itkP[1] = p[1]; itkP[2] = p[2]; - itk::Index<3> idx; - mask->TransformPhysicalPointToIndex(itkP, idx); - - if ( mask->GetPixel(idx) != 0 && mask->GetLargestPossibleRegion().IsInside(idx) && !invert ) - { - vtkIdType id = vtkNewPoints->InsertNextPoint(p); - container->GetPointIds()->InsertNextId(id); - newNumPoints++; - } - else if ( (mask->GetPixel(idx) == 0 || !mask->GetLargestPossibleRegion().IsInside(idx)) && invert ) - { - vtkIdType id = vtkNewPoints->InsertNextPoint(p); - container->GetPointIds()->InsertNextId(id); - newNumPoints++; - } - else if (newNumPoints>0) - { - vtkNewCells->InsertNextCell(container); - - newNumPoints = 0; - container = vtkSmartPointer::New(); - } - } + vtkIdType id = vtkNewPoints->InsertNextPoint(p); + container->GetPointIds()->InsertNextId(id); + newNumPoints++; + } + else if ( (mask->GetPixel(idx) == 0 || !mask->GetLargestPossibleRegion().IsInside(idx)) && invert ) + { + vtkIdType id = vtkNewPoints->InsertNextPoint(p); + container->GetPointIds()->InsertNextId(id); + newNumPoints++; + } + else if (newNumPoints>0) + { + vtkNewCells->InsertNextCell(container); - if (newNumPoints>0) - vtkNewCells->InsertNextCell(container); + newNumPoints = 0; + container = vtkSmartPointer::New(); } + } + if (newNumPoints>0) + vtkNewCells->InsertNextCell(container); } - if (vtkNewCells->GetNumberOfCells()<=0) - return nullptr; + } - vtkSmartPointer newPolyData = vtkSmartPointer::New(); - newPolyData->SetPoints(vtkNewPoints); - newPolyData->SetLines(vtkNewCells); - mitk::FiberBundle::Pointer newFib = mitk::FiberBundle::New(newPolyData); - newFib->Compress(0.1); - return newFib; + if (vtkNewCells->GetNumberOfCells()<=0) + return nullptr; + + vtkSmartPointer newPolyData = vtkSmartPointer::New(); + newPolyData->SetPoints(vtkNewPoints); + newPolyData->SetLines(vtkNewCells); + mitk::FiberBundle::Pointer newFib = mitk::FiberBundle::New(newPolyData); + newFib->Compress(0.1); + return newFib; } mitk::FiberBundle::Pointer mitk::FiberBundle::ExtractFiberSubset(DataNode* roi, DataStorage* storage) { - if (roi==nullptr || !(dynamic_cast(roi->GetData()) || dynamic_cast(roi->GetData())) ) - return nullptr; + if (roi==nullptr || !(dynamic_cast(roi->GetData()) || dynamic_cast(roi->GetData())) ) + return nullptr; - std::vector tmp = ExtractFiberIdSubset(roi, storage); + std::vector tmp = ExtractFiberIdSubset(roi, storage); - if (tmp.size()<=0) - return mitk::FiberBundle::New(); - vtkSmartPointer pTmp = GeneratePolyDataByIds(tmp); - return mitk::FiberBundle::New(pTmp); + if (tmp.size()<=0) + return mitk::FiberBundle::New(); + vtkSmartPointer pTmp = GeneratePolyDataByIds(tmp); + return mitk::FiberBundle::New(pTmp); } std::vector mitk::FiberBundle::ExtractFiberIdSubset(DataNode *roi, DataStorage* storage) { - std::vector result; - if (roi==nullptr || roi->GetData()==nullptr) - return result; - - mitk::PlanarFigureComposite::Pointer pfc = dynamic_cast(roi->GetData()); - if (!pfc.IsNull()) // handle composite - { - DataStorage::SetOfObjects::ConstPointer children = storage->GetDerivations(roi); - if (children->size()==0) - return result; + std::vector result; + if (roi==nullptr || roi->GetData()==nullptr) + return result; - switch (pfc->getOperationType()) - { - case 0: // AND - { - MITK_INFO << "AND"; - result = this->ExtractFiberIdSubset(children->ElementAt(0), storage); - std::vector::iterator it; - for (unsigned int i=1; iSize(); ++i) - { - std::vector inRoi = this->ExtractFiberIdSubset(children->ElementAt(i), storage); + mitk::PlanarFigureComposite::Pointer pfc = dynamic_cast(roi->GetData()); + if (!pfc.IsNull()) // handle composite + { + DataStorage::SetOfObjects::ConstPointer children = storage->GetDerivations(roi); + if (children->size()==0) + return result; - std::vector rest(std::min(result.size(),inRoi.size())); - it = std::set_intersection(result.begin(), result.end(), inRoi.begin(), inRoi.end(), rest.begin() ); - rest.resize( it - rest.begin() ); - result = rest; - } - break; - } - case 1: // OR - { - MITK_INFO << "OR"; - result = ExtractFiberIdSubset(children->ElementAt(0), storage); - std::vector::iterator it; - for (unsigned int i=1; iSize(); ++i) - { - it = result.end(); - std::vector inRoi = ExtractFiberIdSubset(children->ElementAt(i), storage); - result.insert(it, inRoi.begin(), inRoi.end()); - } + switch (pfc->getOperationType()) + { + case 0: // AND + { + MITK_INFO << "AND"; + result = this->ExtractFiberIdSubset(children->ElementAt(0), storage); + std::vector::iterator it; + for (unsigned int i=1; iSize(); ++i) + { + std::vector inRoi = this->ExtractFiberIdSubset(children->ElementAt(i), storage); + + std::vector rest(std::min(result.size(),inRoi.size())); + it = std::set_intersection(result.begin(), result.end(), inRoi.begin(), inRoi.end(), rest.begin() ); + rest.resize( it - rest.begin() ); + result = rest; + } + break; + } + case 1: // OR + { + MITK_INFO << "OR"; + result = ExtractFiberIdSubset(children->ElementAt(0), storage); + std::vector::iterator it; + for (unsigned int i=1; iSize(); ++i) + { + it = result.end(); + std::vector inRoi = ExtractFiberIdSubset(children->ElementAt(i), storage); + result.insert(it, inRoi.begin(), inRoi.end()); + } + + // remove duplicates + sort(result.begin(), result.end()); + it = unique(result.begin(), result.end()); + result.resize( it - result.begin() ); + break; + } + case 2: // NOT + { + MITK_INFO << "NOT"; + for(long i=0; iGetNumFibers(); i++) + result.push_back(i); + + std::vector::iterator it; + for (unsigned int i=0; iSize(); ++i) + { + std::vector inRoi = ExtractFiberIdSubset(children->ElementAt(i), storage); + + std::vector rest(result.size()-inRoi.size()); + it = std::set_difference(result.begin(), result.end(), inRoi.begin(), inRoi.end(), rest.begin() ); + rest.resize( it - rest.begin() ); + result = rest; + } + break; + } + } + } + else if ( dynamic_cast(roi->GetData()) ) // actual extraction + { + if ( dynamic_cast(roi->GetData()) ) + { + mitk::PlanarFigure::Pointer planarPoly = dynamic_cast(roi->GetData()); + + //create vtkPolygon using controlpoints from planarFigure polygon + vtkSmartPointer polygonVtk = vtkSmartPointer::New(); + for (unsigned int i=0; iGetNumberOfControlPoints(); ++i) + { + itk::Point p = planarPoly->GetWorldControlPoint(i); + vtkIdType id = polygonVtk->GetPoints()->InsertNextPoint(p[0], p[1], p[2] ); + polygonVtk->GetPointIds()->InsertNextId(id); + } + + MITK_INFO << "Extracting with polygon"; + boost::progress_display disp(m_NumFibers); + for (int i=0; iGetCell(i); + int numPoints = cell->GetNumberOfPoints(); + vtkPoints* points = cell->GetPoints(); - // remove duplicates - sort(result.begin(), result.end()); - it = unique(result.begin(), result.end()); - result.resize( it - result.begin() ); - break; - } - case 2: // NOT + for (int j=0; jGetNumFibers(); i++) - result.push_back(i); - - std::vector::iterator it; - for (unsigned int i=0; iSize(); ++i) - { - std::vector inRoi = ExtractFiberIdSubset(children->ElementAt(i), storage); - - std::vector rest(result.size()-inRoi.size()); - it = std::set_difference(result.begin(), result.end(), inRoi.begin(), inRoi.end(), rest.begin() ); - rest.resize( it - rest.begin() ); - result = rest; - } + // Inputs + double p1[3] = {0,0,0}; + points->GetPoint(j, p1); + double p2[3] = {0,0,0}; + points->GetPoint(j+1, p2); + double tolerance = 0.001; + + // Outputs + double t = 0; // Parametric coordinate of intersection (0 (corresponding to p1) to 1 (corresponding to p2)) + double x[3] = {0,0,0}; // The coordinate of the intersection + double pcoords[3] = {0,0,0}; + int subId = 0; + + int iD = polygonVtk->IntersectWithLine(p1, p2, tolerance, t, x, pcoords, subId); + if (iD!=0) + { + result.push_back(i); break; + } } - } + } } - else if ( dynamic_cast(roi->GetData()) ) // actual extraction + else if ( dynamic_cast(roi->GetData()) ) { - if ( dynamic_cast(roi->GetData()) ) - { - mitk::PlanarFigure::Pointer planarPoly = dynamic_cast(roi->GetData()); - - //create vtkPolygon using controlpoints from planarFigure polygon - vtkSmartPointer polygonVtk = vtkSmartPointer::New(); - for (unsigned int i=0; iGetNumberOfControlPoints(); ++i) - { - itk::Point p = planarPoly->GetWorldControlPoint(i); - vtkIdType id = polygonVtk->GetPoints()->InsertNextPoint(p[0], p[1], p[2] ); - polygonVtk->GetPointIds()->InsertNextId(id); - } + mitk::PlanarFigure::Pointer planarFigure = dynamic_cast(roi->GetData()); + Vector3D planeNormal = planarFigure->GetPlaneGeometry()->GetNormal(); + planeNormal.Normalize(); - MITK_INFO << "Extracting with polygon"; - boost::progress_display disp(m_NumFibers); - for (int i=0; iGetCell(i); - int numPoints = cell->GetNumberOfPoints(); - vtkPoints* points = cell->GetPoints(); - - for (int j=0; jGetPoint(j, p1); - double p2[3] = {0,0,0}; - points->GetPoint(j+1, p2); - double tolerance = 0.001; - - // Outputs - double t = 0; // Parametric coordinate of intersection (0 (corresponding to p1) to 1 (corresponding to p2)) - double x[3] = {0,0,0}; // The coordinate of the intersection - double pcoords[3] = {0,0,0}; - int subId = 0; - - int iD = polygonVtk->IntersectWithLine(p1, p2, tolerance, t, x, pcoords, subId); - if (iD!=0) - { - result.push_back(i); - break; - } - } - } - } - else if ( dynamic_cast(roi->GetData()) ) - { - mitk::PlanarFigure::Pointer planarFigure = dynamic_cast(roi->GetData()); - Vector3D planeNormal = planarFigure->GetPlaneGeometry()->GetNormal(); - planeNormal.Normalize(); + //calculate circle radius + mitk::Point3D V1w = planarFigure->GetWorldControlPoint(0); //centerPoint + mitk::Point3D V2w = planarFigure->GetWorldControlPoint(1); //radiusPoint - //calculate circle radius - mitk::Point3D V1w = planarFigure->GetWorldControlPoint(0); //centerPoint - mitk::Point3D V2w = planarFigure->GetWorldControlPoint(1); //radiusPoint + double radius = V1w.EuclideanDistanceTo(V2w); + radius *= radius; - double radius = V1w.EuclideanDistanceTo(V2w); - radius *= radius; + MITK_INFO << "Extracting with circle"; + boost::progress_display disp(m_NumFibers); + for (int i=0; iGetCell(i); + int numPoints = cell->GetNumberOfPoints(); + vtkPoints* points = cell->GetPoints(); - MITK_INFO << "Extracting with circle"; - boost::progress_display disp(m_NumFibers); - for (int i=0; iGetPoint(j, p1); + double p2[3] = {0,0,0}; + points->GetPoint(j+1, p2); + + // Outputs + double t = 0; // Parametric coordinate of intersection (0 (corresponding to p1) to 1 (corresponding to p2)) + double x[3] = {0,0,0}; // The coordinate of the intersection + + int iD = vtkPlane::IntersectWithLine(p1,p2,planeNormal.GetDataPointer(),V1w.GetDataPointer(),t,x); + + if (iD!=0) + { + double dist = (x[0]-V1w[0])*(x[0]-V1w[0])+(x[1]-V1w[1])*(x[1]-V1w[1])+(x[2]-V1w[2])*(x[2]-V1w[2]); + if( dist <= radius) { - ++disp ; - vtkCell* cell = m_FiberPolyData->GetCell(i); - int numPoints = cell->GetNumberOfPoints(); - vtkPoints* points = cell->GetPoints(); - - for (int j=0; jGetPoint(j, p1); - double p2[3] = {0,0,0}; - points->GetPoint(j+1, p2); - - // Outputs - double t = 0; // Parametric coordinate of intersection (0 (corresponding to p1) to 1 (corresponding to p2)) - double x[3] = {0,0,0}; // The coordinate of the intersection - - int iD = vtkPlane::IntersectWithLine(p1,p2,planeNormal.GetDataPointer(),V1w.GetDataPointer(),t,x); - - if (iD!=0) - { - double dist = (x[0]-V1w[0])*(x[0]-V1w[0])+(x[1]-V1w[1])*(x[1]-V1w[1])+(x[2]-V1w[2])*(x[2]-V1w[2]); - if( dist <= radius) - { - result.push_back(i); - break; - } - } - } + result.push_back(i); + break; } + } } - return result; + } } - return result; + } + + return result; } void mitk::FiberBundle::UpdateFiberGeometry() { - vtkSmartPointer cleaner = vtkSmartPointer::New(); - cleaner->SetInputData(m_FiberPolyData); - cleaner->PointMergingOff(); - cleaner->Update(); - m_FiberPolyData = cleaner->GetOutput(); - - m_FiberLengths.clear(); - m_MeanFiberLength = 0; - m_MedianFiberLength = 0; - m_LengthStDev = 0; - m_NumFibers = m_FiberPolyData->GetNumberOfCells(); - - if (m_FiberColors==nullptr || m_FiberColors->GetNumberOfTuples()!=m_FiberPolyData->GetNumberOfPoints()) - this->ColorFibersByOrientation(); + vtkSmartPointer cleaner = vtkSmartPointer::New(); + cleaner->SetInputData(m_FiberPolyData); + cleaner->PointMergingOff(); + cleaner->Update(); + m_FiberPolyData = cleaner->GetOutput(); + + m_FiberLengths.clear(); + m_MeanFiberLength = 0; + m_MedianFiberLength = 0; + m_LengthStDev = 0; + m_NumFibers = m_FiberPolyData->GetNumberOfCells(); + + if (m_FiberColors==nullptr || m_FiberColors->GetNumberOfTuples()!=m_FiberPolyData->GetNumberOfPoints()) + this->ColorFibersByOrientation(); - if (m_FiberWeights->GetSize()!=m_NumFibers) + if (m_FiberWeights->GetSize()!=m_NumFibers) + { + m_FiberWeights = vtkSmartPointer::New(); + m_FiberWeights->SetName("FIBER_WEIGHTS"); + m_FiberWeights->SetNumberOfValues(m_NumFibers); + this->SetFiberWeights(1); + } + + if (m_NumFibers<=0) // no fibers present; apply default geometry + { + m_MinFiberLength = 0; + m_MaxFiberLength = 0; + mitk::Geometry3D::Pointer geometry = mitk::Geometry3D::New(); + geometry->SetImageGeometry(false); + float b[] = {0, 1, 0, 1, 0, 1}; + geometry->SetFloatBounds(b); + SetGeometry(geometry); + return; + } + double b[6]; + m_FiberPolyData->GetBounds(b); + + // calculate statistics + for (int i=0; iGetNumberOfCells(); i++) + { + vtkCell* cell = m_FiberPolyData->GetCell(i); + int p = cell->GetNumberOfPoints(); + vtkPoints* points = cell->GetPoints(); + float length = 0; + for (int j=0; j::New(); - m_FiberWeights->SetName("FIBER_WEIGHTS"); - m_FiberWeights->SetNumberOfValues(m_NumFibers); - this->SetFiberWeights(1); - } + double p1[3]; + points->GetPoint(j, p1); + double p2[3]; + points->GetPoint(j+1, p2); - if (m_NumFibers<=0) // no fibers present; apply default geometry - { - m_MinFiberLength = 0; - m_MaxFiberLength = 0; - mitk::Geometry3D::Pointer geometry = mitk::Geometry3D::New(); - geometry->SetImageGeometry(false); - float b[] = {0, 1, 0, 1, 0, 1}; - geometry->SetFloatBounds(b); - SetGeometry(geometry); - return; + float dist = std::sqrt((p1[0]-p2[0])*(p1[0]-p2[0])+(p1[1]-p2[1])*(p1[1]-p2[1])+(p1[2]-p2[2])*(p1[2]-p2[2])); + length += dist; } - double b[6]; - m_FiberPolyData->GetBounds(b); - - // calculate statistics - for (int i=0; iGetNumberOfCells(); i++) + m_FiberLengths.push_back(length); + m_MeanFiberLength += length; + if (i==0) { - vtkCell* cell = m_FiberPolyData->GetCell(i); - int p = cell->GetNumberOfPoints(); - vtkPoints* points = cell->GetPoints(); - float length = 0; - for (int j=0; jGetPoint(j, p1); - double p2[3]; - points->GetPoint(j+1, p2); - - float dist = std::sqrt((p1[0]-p2[0])*(p1[0]-p2[0])+(p1[1]-p2[1])*(p1[1]-p2[1])+(p1[2]-p2[2])*(p1[2]-p2[2])); - length += dist; - } - m_FiberLengths.push_back(length); - m_MeanFiberLength += length; - if (i==0) - { - m_MinFiberLength = length; - m_MaxFiberLength = length; - } - else - { - if (lengthm_MaxFiberLength) - m_MaxFiberLength = length; - } + m_MinFiberLength = length; + m_MaxFiberLength = length; } - m_MeanFiberLength /= m_NumFibers; - - std::vector< float > sortedLengths = m_FiberLengths; - std::sort(sortedLengths.begin(), sortedLengths.end()); - for (int i=0; i1) - m_LengthStDev /= (m_NumFibers-1); else - m_LengthStDev = 0; - m_LengthStDev = std::sqrt(m_LengthStDev); - m_MedianFiberLength = sortedLengths.at(m_NumFibers/2); + { + if (lengthm_MaxFiberLength) + m_MaxFiberLength = length; + } + } + m_MeanFiberLength /= m_NumFibers; + + std::vector< float > sortedLengths = m_FiberLengths; + std::sort(sortedLengths.begin(), sortedLengths.end()); + for (int i=0; i1) + m_LengthStDev /= (m_NumFibers-1); + else + m_LengthStDev = 0; + m_LengthStDev = std::sqrt(m_LengthStDev); + m_MedianFiberLength = sortedLengths.at(m_NumFibers/2); - mitk::Geometry3D::Pointer geometry = mitk::Geometry3D::New(); - geometry->SetFloatBounds(b); - this->SetGeometry(geometry); + mitk::Geometry3D::Pointer geometry = mitk::Geometry3D::New(); + geometry->SetFloatBounds(b); + this->SetGeometry(geometry); - m_UpdateTime3D.Modified(); - m_UpdateTime2D.Modified(); + m_UpdateTime3D.Modified(); + m_UpdateTime2D.Modified(); } float mitk::FiberBundle::GetFiberWeight(unsigned int fiber) const { - return m_FiberWeights->GetValue(fiber); + return m_FiberWeights->GetValue(fiber); } void mitk::FiberBundle::SetFiberWeights(float newWeight) { - for (int i=0; iGetSize(); i++) - m_FiberWeights->SetValue(i, newWeight); + for (int i=0; iGetSize(); i++) + m_FiberWeights->SetValue(i, newWeight); } void mitk::FiberBundle::SetFiberWeights(vtkSmartPointer weights) { - if (m_NumFibers!=weights->GetSize()) - { - MITK_INFO << "Weights array not equal to number of fibers!"; - return; - } + if (m_NumFibers!=weights->GetSize()) + { + MITK_INFO << "Weights array not equal to number of fibers!"; + return; + } - for (int i=0; iGetSize(); i++) - m_FiberWeights->SetValue(i, weights->GetValue(i)); + for (int i=0; iGetSize(); i++) + m_FiberWeights->SetValue(i, weights->GetValue(i)); - m_FiberWeights->SetName("FIBER_WEIGHTS"); + m_FiberWeights->SetName("FIBER_WEIGHTS"); } void mitk::FiberBundle::SetFiberWeight(unsigned int fiber, float weight) { - m_FiberWeights->SetValue(fiber, weight); + m_FiberWeights->SetValue(fiber, weight); } void mitk::FiberBundle::SetFiberColors(vtkSmartPointer fiberColors) { - for(long i=0; iGetNumberOfPoints(); ++i) - { - unsigned char source[4] = {0,0,0,0}; - fiberColors->GetTupleValue(i, source); - - unsigned char target[4] = {0,0,0,0}; - target[0] = source[0]; - target[1] = source[1]; - target[2] = source[2]; - target[3] = source[3]; - m_FiberColors->InsertTupleValue(i, target); - } - m_UpdateTime3D.Modified(); - m_UpdateTime2D.Modified(); + for(long i=0; iGetNumberOfPoints(); ++i) + { + unsigned char source[4] = {0,0,0,0}; + fiberColors->GetTupleValue(i, source); + + unsigned char target[4] = {0,0,0,0}; + target[0] = source[0]; + target[1] = source[1]; + target[2] = source[2]; + target[3] = source[3]; + m_FiberColors->InsertTupleValue(i, target); + } + m_UpdateTime3D.Modified(); + m_UpdateTime2D.Modified(); } itk::Matrix< double, 3, 3 > mitk::FiberBundle::TransformMatrix(itk::Matrix< double, 3, 3 > m, double rx, double ry, double rz) { - rx = rx*M_PI/180; - ry = ry*M_PI/180; - rz = rz*M_PI/180; + rx = rx*M_PI/180; + ry = ry*M_PI/180; + rz = rz*M_PI/180; - itk::Matrix< double, 3, 3 > rotX; rotX.SetIdentity(); - rotX[1][1] = cos(rx); - rotX[2][2] = rotX[1][1]; - rotX[1][2] = -sin(rx); - rotX[2][1] = -rotX[1][2]; + itk::Matrix< double, 3, 3 > rotX; rotX.SetIdentity(); + rotX[1][1] = cos(rx); + rotX[2][2] = rotX[1][1]; + rotX[1][2] = -sin(rx); + rotX[2][1] = -rotX[1][2]; - itk::Matrix< double, 3, 3 > rotY; rotY.SetIdentity(); - rotY[0][0] = cos(ry); - rotY[2][2] = rotY[0][0]; - rotY[0][2] = sin(ry); - rotY[2][0] = -rotY[0][2]; + itk::Matrix< double, 3, 3 > rotY; rotY.SetIdentity(); + rotY[0][0] = cos(ry); + rotY[2][2] = rotY[0][0]; + rotY[0][2] = sin(ry); + rotY[2][0] = -rotY[0][2]; - itk::Matrix< double, 3, 3 > rotZ; rotZ.SetIdentity(); - rotZ[0][0] = cos(rz); - rotZ[1][1] = rotZ[0][0]; - rotZ[0][1] = -sin(rz); - rotZ[1][0] = -rotZ[0][1]; + itk::Matrix< double, 3, 3 > rotZ; rotZ.SetIdentity(); + rotZ[0][0] = cos(rz); + rotZ[1][1] = rotZ[0][0]; + rotZ[0][1] = -sin(rz); + rotZ[1][0] = -rotZ[0][1]; - itk::Matrix< double, 3, 3 > rot = rotZ*rotY*rotX; + itk::Matrix< double, 3, 3 > rot = rotZ*rotY*rotX; - m = rot*m; + m = rot*m; - return m; + return m; } itk::Point mitk::FiberBundle::TransformPoint(vnl_vector_fixed< double, 3 > point, double rx, double ry, double rz, double tx, double ty, double tz) { - rx = rx*M_PI/180; - ry = ry*M_PI/180; - rz = rz*M_PI/180; - - vnl_matrix_fixed< double, 3, 3 > rotX; rotX.set_identity(); - rotX[1][1] = cos(rx); - rotX[2][2] = rotX[1][1]; - rotX[1][2] = -sin(rx); - rotX[2][1] = -rotX[1][2]; - - vnl_matrix_fixed< double, 3, 3 > rotY; rotY.set_identity(); - rotY[0][0] = cos(ry); - rotY[2][2] = rotY[0][0]; - rotY[0][2] = sin(ry); - rotY[2][0] = -rotY[0][2]; - - vnl_matrix_fixed< double, 3, 3 > rotZ; rotZ.set_identity(); - rotZ[0][0] = cos(rz); - rotZ[1][1] = rotZ[0][0]; - rotZ[0][1] = -sin(rz); - rotZ[1][0] = -rotZ[0][1]; - - vnl_matrix_fixed< double, 3, 3 > rot = rotZ*rotY*rotX; - - mitk::BaseGeometry::Pointer geom = this->GetGeometry(); - mitk::Point3D center = geom->GetCenter(); - - point[0] -= center[0]; - point[1] -= center[1]; - point[2] -= center[2]; - point = rot*point; - point[0] += center[0]+tx; - point[1] += center[1]+ty; - point[2] += center[2]+tz; - itk::Point out; out[0] = point[0]; out[1] = point[1]; out[2] = point[2]; - return out; + rx = rx*M_PI/180; + ry = ry*M_PI/180; + rz = rz*M_PI/180; + + vnl_matrix_fixed< double, 3, 3 > rotX; rotX.set_identity(); + rotX[1][1] = cos(rx); + rotX[2][2] = rotX[1][1]; + rotX[1][2] = -sin(rx); + rotX[2][1] = -rotX[1][2]; + + vnl_matrix_fixed< double, 3, 3 > rotY; rotY.set_identity(); + rotY[0][0] = cos(ry); + rotY[2][2] = rotY[0][0]; + rotY[0][2] = sin(ry); + rotY[2][0] = -rotY[0][2]; + + vnl_matrix_fixed< double, 3, 3 > rotZ; rotZ.set_identity(); + rotZ[0][0] = cos(rz); + rotZ[1][1] = rotZ[0][0]; + rotZ[0][1] = -sin(rz); + rotZ[1][0] = -rotZ[0][1]; + + vnl_matrix_fixed< double, 3, 3 > rot = rotZ*rotY*rotX; + + mitk::BaseGeometry::Pointer geom = this->GetGeometry(); + mitk::Point3D center = geom->GetCenter(); + + point[0] -= center[0]; + point[1] -= center[1]; + point[2] -= center[2]; + point = rot*point; + point[0] += center[0]+tx; + point[1] += center[1]+ty; + point[2] += center[2]+tz; + itk::Point out; out[0] = point[0]; out[1] = point[1]; out[2] = point[2]; + return out; } void mitk::FiberBundle::TransformFibers(double rx, double ry, double rz, double tx, double ty, double tz) { - rx = rx*M_PI/180; - ry = ry*M_PI/180; - rz = rz*M_PI/180; - - vnl_matrix_fixed< double, 3, 3 > rotX; rotX.set_identity(); - rotX[1][1] = cos(rx); - rotX[2][2] = rotX[1][1]; - rotX[1][2] = -sin(rx); - rotX[2][1] = -rotX[1][2]; - - vnl_matrix_fixed< double, 3, 3 > rotY; rotY.set_identity(); - rotY[0][0] = cos(ry); - rotY[2][2] = rotY[0][0]; - rotY[0][2] = sin(ry); - rotY[2][0] = -rotY[0][2]; - - vnl_matrix_fixed< double, 3, 3 > rotZ; rotZ.set_identity(); - rotZ[0][0] = cos(rz); - rotZ[1][1] = rotZ[0][0]; - rotZ[0][1] = -sin(rz); - rotZ[1][0] = -rotZ[0][1]; - - vnl_matrix_fixed< double, 3, 3 > rot = rotZ*rotY*rotX; - - mitk::BaseGeometry::Pointer geom = this->GetGeometry(); - mitk::Point3D center = geom->GetCenter(); - - vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); - vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); - - for (int i=0; i rotX; rotX.set_identity(); + rotX[1][1] = cos(rx); + rotX[2][2] = rotX[1][1]; + rotX[1][2] = -sin(rx); + rotX[2][1] = -rotX[1][2]; + + vnl_matrix_fixed< double, 3, 3 > rotY; rotY.set_identity(); + rotY[0][0] = cos(ry); + rotY[2][2] = rotY[0][0]; + rotY[0][2] = sin(ry); + rotY[2][0] = -rotY[0][2]; + + vnl_matrix_fixed< double, 3, 3 > rotZ; rotZ.set_identity(); + rotZ[0][0] = cos(rz); + rotZ[1][1] = rotZ[0][0]; + rotZ[0][1] = -sin(rz); + rotZ[1][0] = -rotZ[0][1]; + + vnl_matrix_fixed< double, 3, 3 > rot = rotZ*rotY*rotX; + + mitk::BaseGeometry::Pointer geom = this->GetGeometry(); + mitk::Point3D center = geom->GetCenter(); + + vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); + vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); + + for (int i=0; iGetCell(i); + int numPoints = cell->GetNumberOfPoints(); + vtkPoints* points = cell->GetPoints(); + + vtkSmartPointer container = vtkSmartPointer::New(); + for (int j=0; jGetCell(i); - int numPoints = cell->GetNumberOfPoints(); - vtkPoints* points = cell->GetPoints(); - - vtkSmartPointer container = vtkSmartPointer::New(); - for (int j=0; jGetPoint(j); - vnl_vector_fixed< double, 3 > dir; - dir[0] = p[0]-center[0]; - dir[1] = p[1]-center[1]; - dir[2] = p[2]-center[2]; - dir = rot*dir; - dir[0] += center[0]+tx; - dir[1] += center[1]+ty; - dir[2] += center[2]+tz; - vtkIdType id = vtkNewPoints->InsertNextPoint(dir.data_block()); - container->GetPointIds()->InsertNextId(id); - } - vtkNewCells->InsertNextCell(container); + double* p = points->GetPoint(j); + vnl_vector_fixed< double, 3 > dir; + dir[0] = p[0]-center[0]; + dir[1] = p[1]-center[1]; + dir[2] = p[2]-center[2]; + dir = rot*dir; + dir[0] += center[0]+tx; + dir[1] += center[1]+ty; + dir[2] += center[2]+tz; + vtkIdType id = vtkNewPoints->InsertNextPoint(dir.data_block()); + container->GetPointIds()->InsertNextId(id); } + vtkNewCells->InsertNextCell(container); + } - m_FiberPolyData = vtkSmartPointer::New(); - m_FiberPolyData->SetPoints(vtkNewPoints); - m_FiberPolyData->SetLines(vtkNewCells); - this->SetFiberPolyData(m_FiberPolyData, true); + m_FiberPolyData = vtkSmartPointer::New(); + m_FiberPolyData->SetPoints(vtkNewPoints); + m_FiberPolyData->SetLines(vtkNewCells); + this->SetFiberPolyData(m_FiberPolyData, true); } void mitk::FiberBundle::RotateAroundAxis(double x, double y, double z) { - x = x*M_PI/180; - y = y*M_PI/180; - z = z*M_PI/180; - - vnl_matrix_fixed< double, 3, 3 > rotX; rotX.set_identity(); - rotX[1][1] = cos(x); - rotX[2][2] = rotX[1][1]; - rotX[1][2] = -sin(x); - rotX[2][1] = -rotX[1][2]; - - vnl_matrix_fixed< double, 3, 3 > rotY; rotY.set_identity(); - rotY[0][0] = cos(y); - rotY[2][2] = rotY[0][0]; - rotY[0][2] = sin(y); - rotY[2][0] = -rotY[0][2]; - - vnl_matrix_fixed< double, 3, 3 > rotZ; rotZ.set_identity(); - rotZ[0][0] = cos(z); - rotZ[1][1] = rotZ[0][0]; - rotZ[0][1] = -sin(z); - rotZ[1][0] = -rotZ[0][1]; - - mitk::BaseGeometry::Pointer geom = this->GetGeometry(); - mitk::Point3D center = geom->GetCenter(); - - vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); - vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); - - for (int i=0; i rotX; rotX.set_identity(); + rotX[1][1] = cos(x); + rotX[2][2] = rotX[1][1]; + rotX[1][2] = -sin(x); + rotX[2][1] = -rotX[1][2]; + + vnl_matrix_fixed< double, 3, 3 > rotY; rotY.set_identity(); + rotY[0][0] = cos(y); + rotY[2][2] = rotY[0][0]; + rotY[0][2] = sin(y); + rotY[2][0] = -rotY[0][2]; + + vnl_matrix_fixed< double, 3, 3 > rotZ; rotZ.set_identity(); + rotZ[0][0] = cos(z); + rotZ[1][1] = rotZ[0][0]; + rotZ[0][1] = -sin(z); + rotZ[1][0] = -rotZ[0][1]; + + mitk::BaseGeometry::Pointer geom = this->GetGeometry(); + mitk::Point3D center = geom->GetCenter(); + + vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); + vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); + + for (int i=0; iGetCell(i); + int numPoints = cell->GetNumberOfPoints(); + vtkPoints* points = cell->GetPoints(); + + vtkSmartPointer container = vtkSmartPointer::New(); + for (int j=0; jGetCell(i); - int numPoints = cell->GetNumberOfPoints(); - vtkPoints* points = cell->GetPoints(); - - vtkSmartPointer container = vtkSmartPointer::New(); - for (int j=0; jGetPoint(j); - vnl_vector_fixed< double, 3 > dir; - dir[0] = p[0]-center[0]; - dir[1] = p[1]-center[1]; - dir[2] = p[2]-center[2]; - dir = rotZ*rotY*rotX*dir; - dir[0] += center[0]; - dir[1] += center[1]; - dir[2] += center[2]; - vtkIdType id = vtkNewPoints->InsertNextPoint(dir.data_block()); - container->GetPointIds()->InsertNextId(id); - } - vtkNewCells->InsertNextCell(container); + double* p = points->GetPoint(j); + vnl_vector_fixed< double, 3 > dir; + dir[0] = p[0]-center[0]; + dir[1] = p[1]-center[1]; + dir[2] = p[2]-center[2]; + dir = rotZ*rotY*rotX*dir; + dir[0] += center[0]; + dir[1] += center[1]; + dir[2] += center[2]; + vtkIdType id = vtkNewPoints->InsertNextPoint(dir.data_block()); + container->GetPointIds()->InsertNextId(id); } + vtkNewCells->InsertNextCell(container); + } - m_FiberPolyData = vtkSmartPointer::New(); - m_FiberPolyData->SetPoints(vtkNewPoints); - m_FiberPolyData->SetLines(vtkNewCells); - this->SetFiberPolyData(m_FiberPolyData, true); + m_FiberPolyData = vtkSmartPointer::New(); + m_FiberPolyData->SetPoints(vtkNewPoints); + m_FiberPolyData->SetLines(vtkNewCells); + this->SetFiberPolyData(m_FiberPolyData, true); } void mitk::FiberBundle::ScaleFibers(double x, double y, double z, bool subtractCenter) { - MITK_INFO << "Scaling fibers"; - boost::progress_display disp(m_NumFibers); + MITK_INFO << "Scaling fibers"; + boost::progress_display disp(m_NumFibers); - mitk::BaseGeometry* geom = this->GetGeometry(); - mitk::Point3D c = geom->GetCenter(); + mitk::BaseGeometry* geom = this->GetGeometry(); + mitk::Point3D c = geom->GetCenter(); - vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); - vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); + vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); + vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); - for (int i=0; iGetCell(i); - int numPoints = cell->GetNumberOfPoints(); - vtkPoints* points = cell->GetPoints(); + for (int i=0; iGetCell(i); + int numPoints = cell->GetNumberOfPoints(); + vtkPoints* points = cell->GetPoints(); - vtkSmartPointer container = vtkSmartPointer::New(); - for (int j=0; jGetPoint(j); - if (subtractCenter) - { - p[0] -= c[0]; p[1] -= c[1]; p[2] -= c[2]; - } - p[0] *= x; - p[1] *= y; - p[2] *= z; - if (subtractCenter) - { - p[0] += c[0]; p[1] += c[1]; p[2] += c[2]; - } - vtkIdType id = vtkNewPoints->InsertNextPoint(p); - container->GetPointIds()->InsertNextId(id); - } - vtkNewCells->InsertNextCell(container); + vtkSmartPointer container = vtkSmartPointer::New(); + for (int j=0; jGetPoint(j); + if (subtractCenter) + { + p[0] -= c[0]; p[1] -= c[1]; p[2] -= c[2]; + } + p[0] *= x; + p[1] *= y; + p[2] *= z; + if (subtractCenter) + { + p[0] += c[0]; p[1] += c[1]; p[2] += c[2]; + } + vtkIdType id = vtkNewPoints->InsertNextPoint(p); + container->GetPointIds()->InsertNextId(id); } + vtkNewCells->InsertNextCell(container); + } - m_FiberPolyData = vtkSmartPointer::New(); - m_FiberPolyData->SetPoints(vtkNewPoints); - m_FiberPolyData->SetLines(vtkNewCells); - this->SetFiberPolyData(m_FiberPolyData, true); + m_FiberPolyData = vtkSmartPointer::New(); + m_FiberPolyData->SetPoints(vtkNewPoints); + m_FiberPolyData->SetLines(vtkNewCells); + this->SetFiberPolyData(m_FiberPolyData, true); } void mitk::FiberBundle::TranslateFibers(double x, double y, double z) { - vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); - vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); + vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); + vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); - for (int i=0; iGetCell(i); - int numPoints = cell->GetNumberOfPoints(); - vtkPoints* points = cell->GetPoints(); + for (int i=0; iGetCell(i); + int numPoints = cell->GetNumberOfPoints(); + vtkPoints* points = cell->GetPoints(); - vtkSmartPointer container = vtkSmartPointer::New(); - for (int j=0; jGetPoint(j); - p[0] += x; - p[1] += y; - p[2] += z; - vtkIdType id = vtkNewPoints->InsertNextPoint(p); - container->GetPointIds()->InsertNextId(id); - } - vtkNewCells->InsertNextCell(container); + vtkSmartPointer container = vtkSmartPointer::New(); + for (int j=0; jGetPoint(j); + p[0] += x; + p[1] += y; + p[2] += z; + vtkIdType id = vtkNewPoints->InsertNextPoint(p); + container->GetPointIds()->InsertNextId(id); } + vtkNewCells->InsertNextCell(container); + } - m_FiberPolyData = vtkSmartPointer::New(); - m_FiberPolyData->SetPoints(vtkNewPoints); - m_FiberPolyData->SetLines(vtkNewCells); - this->SetFiberPolyData(m_FiberPolyData, true); + m_FiberPolyData = vtkSmartPointer::New(); + m_FiberPolyData->SetPoints(vtkNewPoints); + m_FiberPolyData->SetLines(vtkNewCells); + this->SetFiberPolyData(m_FiberPolyData, true); } void mitk::FiberBundle::MirrorFibers(unsigned int axis) { - if (axis>2) - return; + if (axis>2) + return; - MITK_INFO << "Mirroring fibers"; - boost::progress_display disp(m_NumFibers); + MITK_INFO << "Mirroring fibers"; + boost::progress_display disp(m_NumFibers); - vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); - vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); + vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); + vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); - for (int i=0; iGetCell(i); - int numPoints = cell->GetNumberOfPoints(); - vtkPoints* points = cell->GetPoints(); + for (int i=0; iGetCell(i); + int numPoints = cell->GetNumberOfPoints(); + vtkPoints* points = cell->GetPoints(); - vtkSmartPointer container = vtkSmartPointer::New(); - for (int j=0; jGetPoint(j); - p[axis] = -p[axis]; - vtkIdType id = vtkNewPoints->InsertNextPoint(p); - container->GetPointIds()->InsertNextId(id); - } - vtkNewCells->InsertNextCell(container); + vtkSmartPointer container = vtkSmartPointer::New(); + for (int j=0; jGetPoint(j); + p[axis] = -p[axis]; + vtkIdType id = vtkNewPoints->InsertNextPoint(p); + container->GetPointIds()->InsertNextId(id); } + vtkNewCells->InsertNextCell(container); + } - m_FiberPolyData = vtkSmartPointer::New(); - m_FiberPolyData->SetPoints(vtkNewPoints); - m_FiberPolyData->SetLines(vtkNewCells); - this->SetFiberPolyData(m_FiberPolyData, true); + m_FiberPolyData = vtkSmartPointer::New(); + m_FiberPolyData->SetPoints(vtkNewPoints); + m_FiberPolyData->SetLines(vtkNewCells); + this->SetFiberPolyData(m_FiberPolyData, true); } void mitk::FiberBundle::RemoveDir(vnl_vector_fixed dir, double threshold) { - dir.normalize(); - vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); - vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); - - boost::progress_display disp(m_FiberPolyData->GetNumberOfCells()); - for (int i=0; iGetNumberOfCells(); i++) + dir.normalize(); + vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); + vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); + + boost::progress_display disp(m_FiberPolyData->GetNumberOfCells()); + for (int i=0; iGetNumberOfCells(); i++) + { + ++disp ; + vtkCell* cell = m_FiberPolyData->GetCell(i); + int numPoints = cell->GetNumberOfPoints(); + vtkPoints* points = cell->GetPoints(); + + // calculate curvatures + vtkSmartPointer container = vtkSmartPointer::New(); + bool discard = false; + for (int j=0; jGetCell(i); - int numPoints = cell->GetNumberOfPoints(); - vtkPoints* points = cell->GetPoints(); - - // calculate curvatures - vtkSmartPointer container = vtkSmartPointer::New(); - bool discard = false; - for (int j=0; jGetPoint(j, p1); + double p2[3]; + points->GetPoint(j+1, p2); + + vnl_vector_fixed< double, 3 > v1; + v1[0] = p2[0]-p1[0]; + v1[1] = p2[1]-p1[1]; + v1[2] = p2[2]-p1[2]; + if (v1.magnitude()>0.001) + { + v1.normalize(); + + if (fabs(dot_product(v1,dir))>threshold) { - double p1[3]; - points->GetPoint(j, p1); - double p2[3]; - points->GetPoint(j+1, p2); - - vnl_vector_fixed< double, 3 > v1; - v1[0] = p2[0]-p1[0]; - v1[1] = p2[1]-p1[1]; - v1[2] = p2[2]-p1[2]; - if (v1.magnitude()>0.001) - { - v1.normalize(); - - if (fabs(dot_product(v1,dir))>threshold) - { - discard = true; - break; - } - } - } - if (!discard) - { - for (int j=0; jGetPoint(j, p1); - - vtkIdType id = vtkNewPoints->InsertNextPoint(p1); - container->GetPointIds()->InsertNextId(id); - } - vtkNewCells->InsertNextCell(container); + discard = true; + break; } + } } + if (!discard) + { + for (int j=0; jGetPoint(j, p1); + + vtkIdType id = vtkNewPoints->InsertNextPoint(p1); + container->GetPointIds()->InsertNextId(id); + } + vtkNewCells->InsertNextCell(container); + } + } - m_FiberPolyData = vtkSmartPointer::New(); - m_FiberPolyData->SetPoints(vtkNewPoints); - m_FiberPolyData->SetLines(vtkNewCells); + m_FiberPolyData = vtkSmartPointer::New(); + m_FiberPolyData->SetPoints(vtkNewPoints); + m_FiberPolyData->SetLines(vtkNewCells); - this->SetFiberPolyData(m_FiberPolyData, true); + this->SetFiberPolyData(m_FiberPolyData, true); - // UpdateColorCoding(); - // UpdateFiberGeometry(); + // UpdateColorCoding(); + // UpdateFiberGeometry(); } bool mitk::FiberBundle::ApplyCurvatureThreshold(float minRadius, bool deleteFibers) { - if (minRadius<0) - return true; - - vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); - vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); + if (minRadius<0) + return true; - MITK_INFO << "Applying curvature threshold"; - boost::progress_display disp(m_FiberPolyData->GetNumberOfCells()); - for (int i=0; iGetNumberOfCells(); i++) + vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); + vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); + + MITK_INFO << "Applying curvature threshold"; + boost::progress_display disp(m_FiberPolyData->GetNumberOfCells()); + for (int i=0; iGetNumberOfCells(); i++) + { + ++disp ; + vtkCell* cell = m_FiberPolyData->GetCell(i); + int numPoints = cell->GetNumberOfPoints(); + vtkPoints* points = cell->GetPoints(); + + // calculate curvatures + vtkSmartPointer container = vtkSmartPointer::New(); + for (int j=0; jGetCell(i); - int numPoints = cell->GetNumberOfPoints(); - vtkPoints* points = cell->GetPoints(); + double p1[3]; + points->GetPoint(j, p1); + double p2[3]; + points->GetPoint(j+1, p2); + double p3[3]; + points->GetPoint(j+2, p3); - // calculate curvatures - vtkSmartPointer container = vtkSmartPointer::New(); - for (int j=0; jGetPoint(j, p1); - double p2[3]; - points->GetPoint(j+1, p2); - double p3[3]; - points->GetPoint(j+2, p3); + vnl_vector_fixed< float, 3 > v1, v2, v3; - vnl_vector_fixed< float, 3 > v1, v2, v3; + v1[0] = p2[0]-p1[0]; + v1[1] = p2[1]-p1[1]; + v1[2] = p2[2]-p1[2]; - v1[0] = p2[0]-p1[0]; - v1[1] = p2[1]-p1[1]; - v1[2] = p2[2]-p1[2]; + v2[0] = p3[0]-p2[0]; + v2[1] = p3[1]-p2[1]; + v2[2] = p3[2]-p2[2]; - v2[0] = p3[0]-p2[0]; - v2[1] = p3[1]-p2[1]; - v2[2] = p3[2]-p2[2]; + v3[0] = p1[0]-p3[0]; + v3[1] = p1[1]-p3[1]; + v3[2] = p1[2]-p3[2]; - v3[0] = p1[0]-p3[0]; - v3[1] = p1[1]-p3[1]; - v3[2] = p1[2]-p3[2]; + float a = v1.magnitude(); + float b = v2.magnitude(); + float c = v3.magnitude(); + float r = a*b*c/std::sqrt((a+b+c)*(a+b-c)*(b+c-a)*(a-b+c)); // radius of triangle via Heron's formula (area of triangle) - float a = v1.magnitude(); - float b = v2.magnitude(); - float c = v3.magnitude(); - float r = a*b*c/std::sqrt((a+b+c)*(a+b-c)*(b+c-a)*(a-b+c)); // radius of triangle via Heron's formula (area of triangle) + vtkIdType id = vtkNewPoints->InsertNextPoint(p1); + container->GetPointIds()->InsertNextId(id); - vtkIdType id = vtkNewPoints->InsertNextPoint(p1); - container->GetPointIds()->InsertNextId(id); + if (deleteFibers && rInsertNextCell(container); - container = vtkSmartPointer::New(); - } - else if (j==numPoints-3) - { - id = vtkNewPoints->InsertNextPoint(p2); - container->GetPointIds()->InsertNextId(id); - id = vtkNewPoints->InsertNextPoint(p3); - container->GetPointIds()->InsertNextId(id); - vtkNewCells->InsertNextCell(container); - } - } + if (rInsertNextCell(container); + container = vtkSmartPointer::New(); + } + else if (j==numPoints-3) + { + id = vtkNewPoints->InsertNextPoint(p2); + container->GetPointIds()->InsertNextId(id); + id = vtkNewPoints->InsertNextPoint(p3); + container->GetPointIds()->InsertNextId(id); + vtkNewCells->InsertNextCell(container); + } } + } - if (vtkNewCells->GetNumberOfCells()<=0) - return false; + if (vtkNewCells->GetNumberOfCells()<=0) + return false; - m_FiberPolyData = vtkSmartPointer::New(); - m_FiberPolyData->SetPoints(vtkNewPoints); - m_FiberPolyData->SetLines(vtkNewCells); - this->SetFiberPolyData(m_FiberPolyData, true); - return true; + m_FiberPolyData = vtkSmartPointer::New(); + m_FiberPolyData->SetPoints(vtkNewPoints); + m_FiberPolyData->SetLines(vtkNewCells); + this->SetFiberPolyData(m_FiberPolyData, true); + return true; } bool mitk::FiberBundle::RemoveShortFibers(float lengthInMM) { - MITK_INFO << "Removing short fibers"; - if (lengthInMM<=0 || lengthInMMm_MaxFiberLength) // can't remove all fibers - { - MITK_WARN << "Process aborted. No fibers would be left!"; - return false; - } + if (lengthInMM>m_MaxFiberLength) // can't remove all fibers + { + MITK_WARN << "Process aborted. No fibers would be left!"; + return false; + } - vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); - vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); - float min = m_MaxFiberLength; + vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); + vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); + float min = m_MaxFiberLength; - boost::progress_display disp(m_NumFibers); - for (int i=0; iGetCell(i); - int numPoints = cell->GetNumberOfPoints(); - vtkPoints* points = cell->GetPoints(); + boost::progress_display disp(m_NumFibers); + for (int i=0; iGetCell(i); + int numPoints = cell->GetNumberOfPoints(); + vtkPoints* points = cell->GetPoints(); - if (m_FiberLengths.at(i)>=lengthInMM) - { - vtkSmartPointer container = vtkSmartPointer::New(); - for (int j=0; jGetPoint(j); - vtkIdType id = vtkNewPoints->InsertNextPoint(p); - container->GetPointIds()->InsertNextId(id); - } - vtkNewCells->InsertNextCell(container); - if (m_FiberLengths.at(i)=lengthInMM) + { + vtkSmartPointer container = vtkSmartPointer::New(); + for (int j=0; jGetPoint(j); + vtkIdType id = vtkNewPoints->InsertNextPoint(p); + container->GetPointIds()->InsertNextId(id); + } + vtkNewCells->InsertNextCell(container); + if (m_FiberLengths.at(i)GetNumberOfCells()<=0) - return false; + if (vtkNewCells->GetNumberOfCells()<=0) + return false; - m_FiberPolyData = vtkSmartPointer::New(); - m_FiberPolyData->SetPoints(vtkNewPoints); - m_FiberPolyData->SetLines(vtkNewCells); - this->SetFiberPolyData(m_FiberPolyData, true); - return true; + m_FiberPolyData = vtkSmartPointer::New(); + m_FiberPolyData->SetPoints(vtkNewPoints); + m_FiberPolyData->SetLines(vtkNewCells); + this->SetFiberPolyData(m_FiberPolyData, true); + return true; } bool mitk::FiberBundle::RemoveLongFibers(float lengthInMM) { - if (lengthInMM<=0 || lengthInMM>m_MaxFiberLength) - return true; + if (lengthInMM<=0 || lengthInMM>m_MaxFiberLength) + return true; - if (lengthInMM vtkNewPoints = vtkSmartPointer::New(); - vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); + vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); + vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); - MITK_INFO << "Removing long fibers"; - boost::progress_display disp(m_NumFibers); - for (int i=0; iGetCell(i); - int numPoints = cell->GetNumberOfPoints(); - vtkPoints* points = cell->GetPoints(); + MITK_INFO << "Removing long fibers"; + boost::progress_display disp(m_NumFibers); + for (int i=0; iGetCell(i); + int numPoints = cell->GetNumberOfPoints(); + vtkPoints* points = cell->GetPoints(); - if (m_FiberLengths.at(i)<=lengthInMM) - { - vtkSmartPointer container = vtkSmartPointer::New(); - for (int j=0; jGetPoint(j); - vtkIdType id = vtkNewPoints->InsertNextPoint(p); - container->GetPointIds()->InsertNextId(id); - } - vtkNewCells->InsertNextCell(container); - } + if (m_FiberLengths.at(i)<=lengthInMM) + { + vtkSmartPointer container = vtkSmartPointer::New(); + for (int j=0; jGetPoint(j); + vtkIdType id = vtkNewPoints->InsertNextPoint(p); + container->GetPointIds()->InsertNextId(id); + } + vtkNewCells->InsertNextCell(container); } + } - if (vtkNewCells->GetNumberOfCells()<=0) - return false; + if (vtkNewCells->GetNumberOfCells()<=0) + return false; - m_FiberPolyData = vtkSmartPointer::New(); - m_FiberPolyData->SetPoints(vtkNewPoints); - m_FiberPolyData->SetLines(vtkNewCells); - this->SetFiberPolyData(m_FiberPolyData, true); - return true; + m_FiberPolyData = vtkSmartPointer::New(); + m_FiberPolyData->SetPoints(vtkNewPoints); + m_FiberPolyData->SetLines(vtkNewCells); + this->SetFiberPolyData(m_FiberPolyData, true); + return true; } void mitk::FiberBundle::ResampleSpline(float pointDistance, double tension, double continuity, double bias ) { - if (pointDistance<=0) - return; + if (pointDistance<=0) + return; - vtkSmartPointer vtkSmoothPoints = vtkSmartPointer::New(); //in smoothpoints the interpolated points representing a fiber are stored. + vtkSmartPointer vtkSmoothPoints = vtkSmartPointer::New(); //in smoothpoints the interpolated points representing a fiber are stored. - //in vtkcells all polylines are stored, actually all id's of them are stored - vtkSmartPointer vtkSmoothCells = vtkSmartPointer::New(); //cellcontainer for smoothed lines - vtkIdType pointHelperCnt = 0; + //in vtkcells all polylines are stored, actually all id's of them are stored + vtkSmartPointer vtkSmoothCells = vtkSmartPointer::New(); //cellcontainer for smoothed lines + vtkIdType pointHelperCnt = 0; - MITK_INFO << "Smoothing fibers"; - vtkSmartPointer newFiberWeights = vtkSmartPointer::New(); - newFiberWeights->SetName("FIBER_WEIGHTS"); - newFiberWeights->SetNumberOfValues(m_NumFibers); + MITK_INFO << "Smoothing fibers"; + vtkSmartPointer newFiberWeights = vtkSmartPointer::New(); + newFiberWeights->SetName("FIBER_WEIGHTS"); + newFiberWeights->SetNumberOfValues(m_NumFibers); - boost::progress_display disp(m_NumFibers); + boost::progress_display disp(m_NumFibers); #pragma omp parallel for - for (int i=0; i newPoints = vtkSmartPointer::New(); - float length = 0; - float weight = 1; + for (int i=0; i newPoints = vtkSmartPointer::New(); + float length = 0; + float weight = 1; #pragma omp critical - { - length = m_FiberLengths.at(i); - weight = m_FiberWeights->GetValue(i); - ++disp; - vtkCell* cell = m_FiberPolyData->GetCell(i); - int numPoints = cell->GetNumberOfPoints(); - vtkPoints* points = cell->GetPoints(); - for (int j=0; jInsertNextPoint(points->GetPoint(j)); - } + { + length = m_FiberLengths.at(i); + weight = m_FiberWeights->GetValue(i); + ++disp; + vtkCell* cell = m_FiberPolyData->GetCell(i); + int numPoints = cell->GetNumberOfPoints(); + vtkPoints* points = cell->GetPoints(); + for (int j=0; jInsertNextPoint(points->GetPoint(j)); + } - int sampling = std::ceil(length/pointDistance); + int sampling = std::ceil(length/pointDistance); - vtkSmartPointer xSpline = vtkSmartPointer::New(); - vtkSmartPointer ySpline = vtkSmartPointer::New(); - vtkSmartPointer zSpline = vtkSmartPointer::New(); - xSpline->SetDefaultBias(bias); xSpline->SetDefaultTension(tension); xSpline->SetDefaultContinuity(continuity); - ySpline->SetDefaultBias(bias); ySpline->SetDefaultTension(tension); ySpline->SetDefaultContinuity(continuity); - zSpline->SetDefaultBias(bias); zSpline->SetDefaultTension(tension); zSpline->SetDefaultContinuity(continuity); + vtkSmartPointer xSpline = vtkSmartPointer::New(); + vtkSmartPointer ySpline = vtkSmartPointer::New(); + vtkSmartPointer zSpline = vtkSmartPointer::New(); + xSpline->SetDefaultBias(bias); xSpline->SetDefaultTension(tension); xSpline->SetDefaultContinuity(continuity); + ySpline->SetDefaultBias(bias); ySpline->SetDefaultTension(tension); ySpline->SetDefaultContinuity(continuity); + zSpline->SetDefaultBias(bias); zSpline->SetDefaultTension(tension); zSpline->SetDefaultContinuity(continuity); - vtkSmartPointer spline = vtkSmartPointer::New(); - spline->SetXSpline(xSpline); - spline->SetYSpline(ySpline); - spline->SetZSpline(zSpline); - spline->SetPoints(newPoints); + vtkSmartPointer spline = vtkSmartPointer::New(); + spline->SetXSpline(xSpline); + spline->SetYSpline(ySpline); + spline->SetZSpline(zSpline); + spline->SetPoints(newPoints); - vtkSmartPointer functionSource = vtkSmartPointer::New(); - functionSource->SetParametricFunction(spline); - functionSource->SetUResolution(sampling); - functionSource->SetVResolution(sampling); - functionSource->SetWResolution(sampling); - functionSource->Update(); + vtkSmartPointer functionSource = vtkSmartPointer::New(); + functionSource->SetParametricFunction(spline); + functionSource->SetUResolution(sampling); + functionSource->SetVResolution(sampling); + functionSource->SetWResolution(sampling); + functionSource->Update(); - vtkPolyData* outputFunction = functionSource->GetOutput(); - vtkPoints* tmpSmoothPnts = outputFunction->GetPoints(); //smoothPoints of current fiber + vtkPolyData* outputFunction = functionSource->GetOutput(); + vtkPoints* tmpSmoothPnts = outputFunction->GetPoints(); //smoothPoints of current fiber - vtkSmartPointer smoothLine = vtkSmartPointer::New(); - smoothLine->GetPointIds()->SetNumberOfIds(tmpSmoothPnts->GetNumberOfPoints()); + vtkSmartPointer smoothLine = vtkSmartPointer::New(); + smoothLine->GetPointIds()->SetNumberOfIds(tmpSmoothPnts->GetNumberOfPoints()); #pragma omp critical - { - for (int j=0; jGetNumberOfPoints(); j++) - { - smoothLine->GetPointIds()->SetId(j, j+pointHelperCnt); - vtkSmoothPoints->InsertNextPoint(tmpSmoothPnts->GetPoint(j)); - } - - newFiberWeights->SetValue(vtkSmoothCells->GetNumberOfCells(), weight); - vtkSmoothCells->InsertNextCell(smoothLine); - pointHelperCnt += tmpSmoothPnts->GetNumberOfPoints(); - } + { + for (int j=0; jGetNumberOfPoints(); j++) + { + smoothLine->GetPointIds()->SetId(j, j+pointHelperCnt); + vtkSmoothPoints->InsertNextPoint(tmpSmoothPnts->GetPoint(j)); + } + + newFiberWeights->SetValue(vtkSmoothCells->GetNumberOfCells(), weight); + vtkSmoothCells->InsertNextCell(smoothLine); + pointHelperCnt += tmpSmoothPnts->GetNumberOfPoints(); } + } - SetFiberWeights(newFiberWeights); - m_FiberPolyData = vtkSmartPointer::New(); - m_FiberPolyData->SetPoints(vtkSmoothPoints); - m_FiberPolyData->SetLines(vtkSmoothCells); - this->SetFiberPolyData(m_FiberPolyData, true); + SetFiberWeights(newFiberWeights); + m_FiberPolyData = vtkSmartPointer::New(); + m_FiberPolyData->SetPoints(vtkSmoothPoints); + m_FiberPolyData->SetLines(vtkSmoothCells); + this->SetFiberPolyData(m_FiberPolyData, true); } void mitk::FiberBundle::ResampleSpline(float pointDistance) { - ResampleSpline(pointDistance, 0, 0, 0 ); + ResampleSpline(pointDistance, 0, 0, 0 ); } unsigned long mitk::FiberBundle::GetNumberOfPoints() const { - unsigned long points = 0; - for (int i=0; iGetNumberOfCells(); i++) - { - vtkCell* cell = m_FiberPolyData->GetCell(i); - points += cell->GetNumberOfPoints(); - } - return points; + unsigned long points = 0; + for (int i=0; iGetNumberOfCells(); i++) + { + vtkCell* cell = m_FiberPolyData->GetCell(i); + points += cell->GetNumberOfPoints(); + } + return points; } void mitk::FiberBundle::Compress(float error) { - vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); - vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); + vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); + vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); - MITK_INFO << "Compressing fibers"; - unsigned long numRemovedPoints = 0; - boost::progress_display disp(m_FiberPolyData->GetNumberOfCells()); - vtkSmartPointer newFiberWeights = vtkSmartPointer::New(); - newFiberWeights->SetName("FIBER_WEIGHTS"); - newFiberWeights->SetNumberOfValues(m_NumFibers); + MITK_INFO << "Compressing fibers"; + unsigned long numRemovedPoints = 0; + boost::progress_display disp(m_FiberPolyData->GetNumberOfCells()); + vtkSmartPointer newFiberWeights = vtkSmartPointer::New(); + newFiberWeights->SetName("FIBER_WEIGHTS"); + newFiberWeights->SetNumberOfValues(m_NumFibers); #pragma omp parallel for - for (int i=0; iGetNumberOfCells(); i++) - { + for (int i=0; iGetNumberOfCells(); i++) + { - std::vector< vnl_vector_fixed< double, 3 > > vertices; - float weight = 1; + std::vector< vnl_vector_fixed< double, 3 > > vertices; + float weight = 1; #pragma omp critical - { - ++disp; - weight = m_FiberWeights->GetValue(i); - vtkCell* cell = m_FiberPolyData->GetCell(i); - int numPoints = cell->GetNumberOfPoints(); - vtkPoints* points = cell->GetPoints(); + { + ++disp; + weight = m_FiberWeights->GetValue(i); + vtkCell* cell = m_FiberPolyData->GetCell(i); + int numPoints = cell->GetNumberOfPoints(); + vtkPoints* points = cell->GetPoints(); + + for (int j=0; jGetPoint(j, cand); + vnl_vector_fixed< double, 3 > candV; + candV[0]=cand[0]; candV[1]=cand[1]; candV[2]=cand[2]; + vertices.push_back(candV); + } + } - for (int j=0; jGetPoint(j, cand); - vnl_vector_fixed< double, 3 > candV; - candV[0]=cand[0]; candV[1]=cand[1]; candV[2]=cand[2]; - vertices.push_back(candV); - } - } + // calculate curvatures + int numPoints = vertices.size(); + std::vector< int > removedPoints; removedPoints.resize(numPoints, 0); + removedPoints[0]=-1; removedPoints[numPoints-1]=-1; - // calculate curvatures - int numPoints = vertices.size(); - std::vector< int > removedPoints; removedPoints.resize(numPoints, 0); - removedPoints[0]=-1; removedPoints[numPoints-1]=-1; + vtkSmartPointer container = vtkSmartPointer::New(); + int remCounter = 0; - vtkSmartPointer container = vtkSmartPointer::New(); - int remCounter = 0; + bool pointFound = true; + while (pointFound) + { + pointFound = false; + double minError = error; + int removeIndex = -1; - bool pointFound = true; - while (pointFound) + for (int j=0; j candV = vertices.at(j); - for (int j=0; j pred; + for (int k=j-1; k>=0; k--) + if (removedPoints[k]<=0) { - if (removedPoints[j]==0) - { - vnl_vector_fixed< double, 3 > candV = vertices.at(j); - - int validP = -1; - vnl_vector_fixed< double, 3 > pred; - for (int k=j-1; k>=0; k--) - if (removedPoints[k]<=0) - { - pred = vertices.at(k); - validP = k; - break; - } - int validS = -1; - vnl_vector_fixed< double, 3 > succ; - for (int k=j+1; k=0 && validS>=0) - { - double a = (candV-pred).magnitude(); - double b = (candV-succ).magnitude(); - double c = (pred-succ).magnitude(); - double s=0.5*(a+b+c); - double hc=(2.0/c)*sqrt(fabs(s*(s-a)*(s-b)*(s-c))); - - if (hc succ; + for (int k=j+1; k=0 && validS>=0) + { + double a = (candV-pred).magnitude(); + double b = (candV-succ).magnitude(); + double c = (pred-succ).magnitude(); + double s=0.5*(a+b+c); + double hc=(2.0/c)*sqrt(fabs(s*(s-a)*(s-b)*(s-c))); + + if (hcInsertNextPoint(vertices.at(j).data_block()); - container->GetPointIds()->InsertNextId(id); - } + removeIndex = j; + minError = hc; + pointFound = true; } + } } + } + if (pointFound) + { + removedPoints[removeIndex] = 1; + remCounter++; + } + } + + for (int j=0; jSetValue(vtkNewCells->GetNumberOfCells(), weight); - numRemovedPoints += remCounter; - vtkNewCells->InsertNextCell(container); + vtkIdType id = vtkNewPoints->InsertNextPoint(vertices.at(j).data_block()); + container->GetPointIds()->InsertNextId(id); } + } } - if (vtkNewCells->GetNumberOfCells()>0) +#pragma omp critical { - MITK_INFO << "Removed points: " << numRemovedPoints; - SetFiberWeights(newFiberWeights); - m_FiberPolyData = vtkSmartPointer::New(); - m_FiberPolyData->SetPoints(vtkNewPoints); - m_FiberPolyData->SetLines(vtkNewCells); - this->SetFiberPolyData(m_FiberPolyData, true); + newFiberWeights->SetValue(vtkNewCells->GetNumberOfCells(), weight); + numRemovedPoints += remCounter; + vtkNewCells->InsertNextCell(container); } + } + + if (vtkNewCells->GetNumberOfCells()>0) + { + MITK_INFO << "Removed points: " << numRemovedPoints; + SetFiberWeights(newFiberWeights); + m_FiberPolyData = vtkSmartPointer::New(); + m_FiberPolyData->SetPoints(vtkNewPoints); + m_FiberPolyData->SetLines(vtkNewCells); + this->SetFiberPolyData(m_FiberPolyData, true); + } } void mitk::FiberBundle::ResampleLinear(double pointDistance) { - vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); - vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); + vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); + vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); - MITK_INFO << "Resampling fibers (linear)"; - boost::progress_display disp(m_FiberPolyData->GetNumberOfCells()); - vtkSmartPointer newFiberWeights = vtkSmartPointer::New(); - newFiberWeights->SetName("FIBER_WEIGHTS"); - newFiberWeights->SetNumberOfValues(m_NumFibers); + MITK_INFO << "Resampling fibers (linear)"; + boost::progress_display disp(m_FiberPolyData->GetNumberOfCells()); + vtkSmartPointer newFiberWeights = vtkSmartPointer::New(); + newFiberWeights->SetName("FIBER_WEIGHTS"); + newFiberWeights->SetNumberOfValues(m_NumFibers); #pragma omp parallel for - for (int i=0; iGetNumberOfCells(); i++) - { + for (int i=0; iGetNumberOfCells(); i++) + { - std::vector< vnl_vector_fixed< double, 3 > > vertices; - float weight = 1; + std::vector< vnl_vector_fixed< double, 3 > > vertices; + float weight = 1; #pragma omp critical - { - ++disp; - weight = m_FiberWeights->GetValue(i); - vtkCell* cell = m_FiberPolyData->GetCell(i); - int numPoints = cell->GetNumberOfPoints(); - vtkPoints* points = cell->GetPoints(); - - for (int j=0; jGetPoint(j, cand); - vnl_vector_fixed< double, 3 > candV; - candV[0]=cand[0]; candV[1]=cand[1]; candV[2]=cand[2]; - vertices.push_back(candV); - } - } + { + ++disp; + weight = m_FiberWeights->GetValue(i); + vtkCell* cell = m_FiberPolyData->GetCell(i); + int numPoints = cell->GetNumberOfPoints(); + vtkPoints* points = cell->GetPoints(); + + for (int j=0; jGetPoint(j, cand); + vnl_vector_fixed< double, 3 > candV; + candV[0]=cand[0]; candV[1]=cand[1]; candV[2]=cand[2]; + vertices.push_back(candV); + } + } - vtkSmartPointer container = vtkSmartPointer::New(); - vnl_vector_fixed< double, 3 > lastV = vertices.at(0); + vtkSmartPointer container = vtkSmartPointer::New(); + vnl_vector_fixed< double, 3 > lastV = vertices.at(0); #pragma omp critical + { + vtkIdType id = vtkNewPoints->InsertNextPoint(lastV.data_block()); + container->GetPointIds()->InsertNextId(id); + } + for (unsigned int j=1; j vec = vertices.at(j) - lastV; + double new_dist = vec.magnitude(); + + if (new_dist >= pointDistance) + { + vnl_vector_fixed< double, 3 > newV = lastV; + if ( new_dist-pointDistance <= mitk::eps ) { - vtkIdType id = vtkNewPoints->InsertNextPoint(lastV.data_block()); - container->GetPointIds()->InsertNextId(id); + vec.normalize(); + newV += vec * pointDistance; } - for (unsigned int j=1; j vec = vertices.at(j) - lastV; - double new_dist = vec.magnitude(); + // intersection between sphere (radius 'pointDistance', center 'lastV') and line (direction 'd' and point 'p') + vnl_vector_fixed< double, 3 > p = vertices.at(j-1); + vnl_vector_fixed< double, 3 > d = vertices.at(j) - p; - if (new_dist >= pointDistance) - { - vnl_vector_fixed< double, 3 > newV = lastV; - if ( new_dist-pointDistance <= mitk::eps ) - { - vec.normalize(); - newV += vec * pointDistance; - } - else - { - // intersection between sphere (radius 'pointDistance', center 'lastV') and line (direction 'd' and point 'p') - vnl_vector_fixed< double, 3 > p = vertices.at(j-1); - vnl_vector_fixed< double, 3 > d = vertices.at(j) - p; - - double a = d[0]*d[0] + d[1]*d[1] + d[2]*d[2]; - double b = 2 * (d[0] * (p[0] - lastV[0]) + d[1] * (p[1] - lastV[1]) + d[2] * (p[2] - lastV[2])); - double c = (p[0] - lastV[0])*(p[0] - lastV[0]) + (p[1] - lastV[1])*(p[1] - lastV[1]) + (p[2] - lastV[2])*(p[2] - lastV[2]) - pointDistance*pointDistance; - - double v1 =(-b + std::sqrt(b*b-4*a*c))/(2*a); - double v2 =(-b - std::sqrt(b*b-4*a*c))/(2*a); - - if (v1>0) - newV = p + d * v1; - else if (v2>0) - newV = p + d * v2; - else - MITK_INFO << "ERROR1 - linear resampling"; - - j--; - } + double a = d[0]*d[0] + d[1]*d[1] + d[2]*d[2]; + double b = 2 * (d[0] * (p[0] - lastV[0]) + d[1] * (p[1] - lastV[1]) + d[2] * (p[2] - lastV[2])); + double c = (p[0] - lastV[0])*(p[0] - lastV[0]) + (p[1] - lastV[1])*(p[1] - lastV[1]) + (p[2] - lastV[2])*(p[2] - lastV[2]) - pointDistance*pointDistance; -#pragma omp critical - { - vtkIdType id = vtkNewPoints->InsertNextPoint(newV.data_block()); - container->GetPointIds()->InsertNextId(id); - } - lastV = newV; - } - else if (j==vertices.size()-1 && new_dist>0.0001) - { -#pragma omp critical - { - vtkIdType id = vtkNewPoints->InsertNextPoint(vertices.at(j).data_block()); - container->GetPointIds()->InsertNextId(id); - } - } + double v1 =(-b + std::sqrt(b*b-4*a*c))/(2*a); + double v2 =(-b - std::sqrt(b*b-4*a*c))/(2*a); + + if (v1>0) + newV = p + d * v1; + else if (v2>0) + newV = p + d * v2; + else + MITK_INFO << "ERROR1 - linear resampling"; + + j--; } #pragma omp critical { - newFiberWeights->SetValue(vtkNewCells->GetNumberOfCells(), weight); - vtkNewCells->InsertNextCell(container); + vtkIdType id = vtkNewPoints->InsertNextPoint(newV.data_block()); + container->GetPointIds()->InsertNextId(id); + } + lastV = newV; + } + else if (j==vertices.size()-1 && new_dist>0.0001) + { +#pragma omp critical + { + vtkIdType id = vtkNewPoints->InsertNextPoint(vertices.at(j).data_block()); + container->GetPointIds()->InsertNextId(id); } + } } - if (vtkNewCells->GetNumberOfCells()>0) +#pragma omp critical { - SetFiberWeights(newFiberWeights); - m_FiberPolyData = vtkSmartPointer::New(); - m_FiberPolyData->SetPoints(vtkNewPoints); - m_FiberPolyData->SetLines(vtkNewCells); - this->SetFiberPolyData(m_FiberPolyData, true); + newFiberWeights->SetValue(vtkNewCells->GetNumberOfCells(), weight); + vtkNewCells->InsertNextCell(container); } + } + + if (vtkNewCells->GetNumberOfCells()>0) + { + SetFiberWeights(newFiberWeights); + m_FiberPolyData = vtkSmartPointer::New(); + m_FiberPolyData->SetPoints(vtkNewPoints); + m_FiberPolyData->SetLines(vtkNewCells); + this->SetFiberPolyData(m_FiberPolyData, true); + } } // reapply selected colorcoding in case PolyData structure has changed bool mitk::FiberBundle::Equals(mitk::FiberBundle* fib, double eps) { - if (fib==nullptr) - { - MITK_INFO << "Reference bundle is nullptr!"; - return false; - } + if (fib==nullptr) + { + MITK_INFO << "Reference bundle is nullptr!"; + return false; + } - if (m_NumFibers!=fib->GetNumFibers()) - { - MITK_INFO << "Unequal number of fibers!"; - MITK_INFO << m_NumFibers << " vs. " << fib->GetNumFibers(); - return false; - } + if (m_NumFibers!=fib->GetNumFibers()) + { + MITK_INFO << "Unequal number of fibers!"; + MITK_INFO << m_NumFibers << " vs. " << fib->GetNumFibers(); + return false; + } - for (int i=0; iGetCell(i); - int numPoints = cell->GetNumberOfPoints(); - vtkPoints* points = cell->GetPoints(); + for (int i=0; iGetCell(i); + int numPoints = cell->GetNumberOfPoints(); + vtkPoints* points = cell->GetPoints(); - vtkCell* cell2 = fib->GetFiberPolyData()->GetCell(i); - int numPoints2 = cell2->GetNumberOfPoints(); - vtkPoints* points2 = cell2->GetPoints(); + vtkCell* cell2 = fib->GetFiberPolyData()->GetCell(i); + int numPoints2 = cell2->GetNumberOfPoints(); + vtkPoints* points2 = cell2->GetPoints(); - if (numPoints2!=numPoints) - { - MITK_INFO << "Unequal number of points in fiber " << i << "!"; - MITK_INFO << numPoints2 << " vs. " << numPoints; - return false; - } + if (numPoints2!=numPoints) + { + MITK_INFO << "Unequal number of points in fiber " << i << "!"; + MITK_INFO << numPoints2 << " vs. " << numPoints; + return false; + } - for (int j=0; jGetPoint(j); - double* p2 = points2->GetPoint(j); - if (fabs(p1[0]-p2[0])>eps || fabs(p1[1]-p2[1])>eps || fabs(p1[2]-p2[2])>eps) - { - MITK_INFO << "Unequal points in fiber " << i << " at position " << j << "!"; - MITK_INFO << "p1: " << p1[0] << ", " << p1[1] << ", " << p1[2]; - MITK_INFO << "p2: " << p2[0] << ", " << p2[1] << ", " << p2[2]; - return false; - } - } + for (int j=0; jGetPoint(j); + double* p2 = points2->GetPoint(j); + if (fabs(p1[0]-p2[0])>eps || fabs(p1[1]-p2[1])>eps || fabs(p1[2]-p2[2])>eps) + { + MITK_INFO << "Unequal points in fiber " << i << " at position " << j << "!"; + MITK_INFO << "p1: " << p1[0] << ", " << p1[1] << ", " << p1[2]; + MITK_INFO << "p2: " << p2[0] << ", " << p2[1] << ", " << p2[2]; + return false; + } } + } - return true; + return true; } /* ESSENTIAL IMPLEMENTATION OF SUPERCLASS METHODS */ void mitk::FiberBundle::UpdateOutputInformation() { } void mitk::FiberBundle::SetRequestedRegionToLargestPossibleRegion() { } bool mitk::FiberBundle::RequestedRegionIsOutsideOfTheBufferedRegion() { - return false; + return false; } bool mitk::FiberBundle::VerifyRequestedRegion() { - return true; + return true; } void mitk::FiberBundle::SetRequestedRegion(const itk::DataObject* ) { } diff --git a/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundle/mitkFiberBundle.h b/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundle/mitkFiberBundle.h index 5d6ab36f98..bc65bf7c11 100644 --- a/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundle/mitkFiberBundle.h +++ b/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundle/mitkFiberBundle.h @@ -1,176 +1,177 @@ /*=================================================================== 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 _MITK_FiberBundle_H #define _MITK_FiberBundle_H //includes for MITK datastructure #include #include #include #include #include #include #include //includes storing fiberdata #include #include #include #include #include #include namespace mitk { /** * \brief Base Class for Fiber Bundles; */ class MITKFIBERTRACKING_EXPORT FiberBundle : public BaseData { public: typedef itk::Image ItkUcharImgType; // fiber colorcodings static const char* FIBER_ID_ARRAY; virtual void UpdateOutputInformation() override; virtual void SetRequestedRegionToLargestPossibleRegion() override; virtual bool RequestedRegionIsOutsideOfTheBufferedRegion() override; virtual bool VerifyRequestedRegion() override; virtual void SetRequestedRegion(const itk::DataObject*) override; mitkClassMacro( FiberBundle, BaseData ) itkFactorylessNewMacro(Self) itkCloneMacro(Self) mitkNewMacro1Param(Self, vtkSmartPointer) // custom constructor // colorcoding related methods + void ColorFibersByFiberWeights(); void ColorFibersByCurvature(bool minMaxNorm=true); void ColorFibersByScalarMap(mitk::Image::Pointer, bool opacity); template void ColorFibersByScalarMap(const mitk::PixelType pixelType, mitk::Image::Pointer, bool opacity); void ColorFibersByOrientation(); void SetFiberOpacity(vtkDoubleArray *FAValArray); void ResetFiberOpacity(); void SetFiberColors(vtkSmartPointer fiberColors); void SetFiberColors(float r, float g, float b, float alpha=255); vtkSmartPointer GetFiberColors() const { return m_FiberColors; } // fiber compression void Compress(float error = 0.0); // fiber resampling void ResampleSpline(float pointDistance=1); void ResampleSpline(float pointDistance, double tension, double continuity, double bias ); void ResampleLinear(double pointDistance=1); bool RemoveShortFibers(float lengthInMM); bool RemoveLongFibers(float lengthInMM); bool ApplyCurvatureThreshold(float minRadius, bool deleteFibers); void MirrorFibers(unsigned int axis); void RotateAroundAxis(double x, double y, double z); void TranslateFibers(double x, double y, double z); void ScaleFibers(double x, double y, double z, bool subtractCenter=true); void TransformFibers(double rx, double ry, double rz, double tx, double ty, double tz); void RemoveDir(vnl_vector_fixed dir, double threshold); itk::Point TransformPoint(vnl_vector_fixed< double, 3 > point, double rx, double ry, double rz, double tx, double ty, double tz); itk::Matrix< double, 3, 3 > TransformMatrix(itk::Matrix< double, 3, 3 > m, double rx, double ry, double rz); // add/subtract fibers FiberBundle::Pointer AddBundle(FiberBundle* fib); FiberBundle::Pointer SubtractBundle(FiberBundle* fib); // fiber subset extraction FiberBundle::Pointer ExtractFiberSubset(DataNode *roi, DataStorage* storage); std::vector ExtractFiberIdSubset(DataNode* roi, DataStorage* storage); FiberBundle::Pointer ExtractFiberSubset(ItkUcharImgType* mask, bool anyPoint, bool invert=false, bool bothEnds=true, float fraction=0.0); FiberBundle::Pointer RemoveFibersOutside(ItkUcharImgType* mask, bool invert=false); vtkSmartPointer GeneratePolyDataByIds( std::vector ); // TODO: make protected void GenerateFiberIds(); // TODO: make protected // get/set data vtkSmartPointer GetFiberWeights() const { return m_FiberWeights; } float GetFiberWeight(unsigned int fiber) const; void SetFiberWeights(float newWeight); void SetFiberWeight(unsigned int fiber, float weight); void SetFiberWeights(vtkSmartPointer weights); void SetFiberPolyData(vtkSmartPointer, bool updateGeometry = true); vtkSmartPointer GetFiberPolyData() const; itkGetConstMacro( NumFibers, int) //itkGetMacro( FiberSampling, int) itkGetConstMacro( MinFiberLength, float ) itkGetConstMacro( MaxFiberLength, float ) itkGetConstMacro( MeanFiberLength, float ) itkGetConstMacro( MedianFiberLength, float ) itkGetConstMacro( LengthStDev, float ) itkGetConstMacro( UpdateTime2D, itk::TimeStamp ) itkGetConstMacro( UpdateTime3D, itk::TimeStamp ) void RequestUpdate2D(){ m_UpdateTime2D.Modified(); } void RequestUpdate3D(){ m_UpdateTime3D.Modified(); } void RequestUpdate(){ m_UpdateTime2D.Modified(); m_UpdateTime3D.Modified(); } unsigned long GetNumberOfPoints() const; // copy fiber bundle mitk::FiberBundle::Pointer GetDeepCopy(); // compare fiber bundles bool Equals(FiberBundle* fib, double eps=0.01); itkSetMacro( ReferenceGeometry, mitk::BaseGeometry::Pointer ) itkGetConstMacro( ReferenceGeometry, mitk::BaseGeometry::Pointer ) protected: FiberBundle( vtkPolyData* fiberPolyData = nullptr ); virtual ~FiberBundle(); itk::Point GetItkPoint(double point[3]); // calculate geometry from fiber extent void UpdateFiberGeometry(); private: // actual fiber container vtkSmartPointer m_FiberPolyData; // contains fiber ids vtkSmartPointer m_FiberIdDataSet; int m_NumFibers; vtkSmartPointer m_FiberColors; vtkSmartPointer m_FiberWeights; std::vector< float > m_FiberLengths; float m_MinFiberLength; float m_MaxFiberLength; float m_MeanFiberLength; float m_MedianFiberLength; float m_LengthStDev; itk::TimeStamp m_UpdateTime2D; itk::TimeStamp m_UpdateTime3D; mitk::BaseGeometry::Pointer m_ReferenceGeometry; }; } // namespace mitk #endif /* _MITK_FiberBundle_H */ diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingView.cpp index 971cec5b78..92a2e09d0b 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingView.cpp @@ -1,1428 +1,1447 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ // Blueberry #include #include #include // Qmitk #include "QmitkFiberProcessingView.h" // Qt #include // MITK #include #include #include #include #include #include #include #include #include #include #include "usModuleRegistry.h" #include #include "mitkNodePredicateDataType.h" #include #include #include #include // ITK #include #include #include #include #include #include #include #define _USE_MATH_DEFINES #include const std::string QmitkFiberProcessingView::VIEW_ID = "org.mitk.views.fiberprocessing"; const std::string id_DataManager = "org.mitk.views.datamanager"; using namespace mitk; QmitkFiberProcessingView::QmitkFiberProcessingView() : QmitkAbstractView() , m_Controls( 0 ) , m_CircleCounter(0) , m_PolygonCounter(0) , m_UpsamplingFactor(1) { } // Destructor QmitkFiberProcessingView::~QmitkFiberProcessingView() { } void QmitkFiberProcessingView::CreateQtPartControl( QWidget *parent ) { // build up qt view, unless already done if ( !m_Controls ) { // create GUI widgets from the Qt Designer's .ui file m_Controls = new Ui::QmitkFiberProcessingViewControls; m_Controls->setupUi( parent ); connect( m_Controls->m_CircleButton, SIGNAL( clicked() ), this, SLOT( OnDrawCircle() ) ); connect( m_Controls->m_PolygonButton, SIGNAL( clicked() ), this, SLOT( OnDrawPolygon() ) ); connect(m_Controls->PFCompoANDButton, SIGNAL(clicked()), this, SLOT(GenerateAndComposite()) ); connect(m_Controls->PFCompoORButton, SIGNAL(clicked()), this, SLOT(GenerateOrComposite()) ); connect(m_Controls->PFCompoNOTButton, SIGNAL(clicked()), this, SLOT(GenerateNotComposite()) ); connect(m_Controls->m_GenerateRoiImage, SIGNAL(clicked()), this, SLOT(GenerateRoiImage()) ); connect(m_Controls->m_JoinBundles, SIGNAL(clicked()), this, SLOT(JoinBundles()) ); connect(m_Controls->m_SubstractBundles, SIGNAL(clicked()), this, SLOT(SubstractBundles()) ); connect(m_Controls->m_CopyBundle, SIGNAL(clicked()), this, SLOT(CopyBundles()) ); connect(m_Controls->m_ExtractFibersButton, SIGNAL(clicked()), this, SLOT(Extract())); connect(m_Controls->m_RemoveButton, SIGNAL(clicked()), this, SLOT(Remove())); connect(m_Controls->m_ModifyButton, SIGNAL(clicked()), this, SLOT(Modify())); connect(m_Controls->m_ExtractionMethodBox, SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui())); connect(m_Controls->m_RemovalMethodBox, SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui())); connect(m_Controls->m_ModificationMethodBox, SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui())); connect(m_Controls->m_ExtractionBoxMask, SIGNAL(currentIndexChanged(int)), this, SLOT(OnMaskExtractionChanged())); m_Controls->m_ColorMapBox->SetDataStorage(this->GetDataStorage()); mitk::TNodePredicateDataType::Pointer isMitkImage = mitk::TNodePredicateDataType::New(); mitk::NodePredicateDataType::Pointer isDwi = mitk::NodePredicateDataType::New("DiffusionImage"); mitk::NodePredicateDataType::Pointer isDti = mitk::NodePredicateDataType::New("TensorImage"); mitk::NodePredicateDataType::Pointer isQbi = mitk::NodePredicateDataType::New("QBallImage"); mitk::NodePredicateOr::Pointer isDiffusionImage = mitk::NodePredicateOr::New(isDwi, isDti); isDiffusionImage = mitk::NodePredicateOr::New(isDiffusionImage, isQbi); mitk::NodePredicateNot::Pointer noDiffusionImage = mitk::NodePredicateNot::New(isDiffusionImage); mitk::NodePredicateAnd::Pointer finalPredicate = mitk::NodePredicateAnd::New(isMitkImage, noDiffusionImage); m_Controls->m_ColorMapBox->SetPredicate(finalPredicate); m_Controls->label_17->setVisible(false); m_Controls->m_FiberExtractionFractionBox->setVisible(false); } UpdateGui(); } void QmitkFiberProcessingView::OnMaskExtractionChanged() { if (m_Controls->m_ExtractionBoxMask->currentIndex() == 2) { m_Controls->label_17->setVisible(true); m_Controls->m_FiberExtractionFractionBox->setVisible(true); m_Controls->m_BothEnds->setVisible(false); } else { m_Controls->label_17->setVisible(false); m_Controls->m_FiberExtractionFractionBox->setVisible(false); if (m_Controls->m_ExtractionBoxMask->currentIndex() != 3) m_Controls->m_BothEnds->setVisible(true); } } void QmitkFiberProcessingView::SetFocus() { m_Controls->toolBoxx->setFocus(); } void QmitkFiberProcessingView::Modify() { switch (m_Controls->m_ModificationMethodBox->currentIndex()) { case 0: { ResampleSelectedBundlesSpline(); break; } case 1: { ResampleSelectedBundlesLinear(); break; } case 2: { CompressSelectedBundles(); break; } case 3: { DoImageColorCoding(); break; } case 4: { MirrorFibers(); break; } case 5: { WeightFibers(); break; } case 6: { DoCurvatureColorCoding(); break; } + case 7: + { + DoWeightColorCoding(); + break; + } } } void QmitkFiberProcessingView::WeightFibers() { float weight = this->m_Controls->m_BundleWeightBox->value(); for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); fib->SetFiberWeights(weight); } } void QmitkFiberProcessingView::Remove() { switch (m_Controls->m_RemovalMethodBox->currentIndex()) { case 0: { RemoveDir(); break; } case 1: { PruneBundle(); break; } case 2: { ApplyCurvatureThreshold(); break; } case 3: { RemoveWithMask(false); break; } case 4: { RemoveWithMask(true); break; } } } void QmitkFiberProcessingView::Extract() { switch (m_Controls->m_ExtractionMethodBox->currentIndex()) { case 0: { ExtractWithPlanarFigure(); break; } case 1: { switch (m_Controls->m_ExtractionBoxMask->currentIndex()) { { case 0: ExtractWithMask(true, false); break; } { case 1: ExtractWithMask(true, true); break; } { case 2: ExtractWithMask(false, false); break; } { case 3: ExtractWithMask(false, true); break; } } break; } } } void QmitkFiberProcessingView::PruneBundle() { int minLength = this->m_Controls->m_PruneFibersMinBox->value(); int maxLength = this->m_Controls->m_PruneFibersMaxBox->value(); for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); if (!fib->RemoveShortFibers(minLength)) QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers."); else if (!fib->RemoveLongFibers(maxLength)) QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers."); } RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberProcessingView::ApplyCurvatureThreshold() { int angle = this->m_Controls->m_CurvSpinBox->value(); int dist = this->m_Controls->m_CurvDistanceSpinBox->value(); std::vector< DataNode::Pointer > nodes = m_SelectedFB; for (auto node : nodes) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); itk::FiberCurvatureFilter::Pointer filter = itk::FiberCurvatureFilter::New(); filter->SetInputFiberBundle(fib); filter->SetAngularDeviation(angle); filter->SetDistance(dist); filter->SetRemoveFibers(m_Controls->m_RemoveCurvedFibersBox->isChecked()); filter->Update(); mitk::FiberBundle::Pointer newFib = filter->GetOutputFiberBundle(); if (newFib->GetNumFibers()>0) { newFib->ColorFibersByOrientation(); node->SetData(newFib); } else QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers."); } RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberProcessingView::RemoveDir() { for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); vnl_vector_fixed dir; dir[0] = m_Controls->m_ExtractDirX->value(); dir[1] = m_Controls->m_ExtractDirY->value(); dir[2] = m_Controls->m_ExtractDirZ->value(); fib->RemoveDir(dir,cos((float)m_Controls->m_ExtractAngle->value()*M_PI/180)); } RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberProcessingView::RemoveWithMask(bool removeInside) { if (m_MaskImageNode.IsNull()) return; mitk::Image::Pointer mitkMask = dynamic_cast(m_MaskImageNode->GetData()); for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); QString name(node->GetName().c_str()); itkUCharImageType::Pointer mask = itkUCharImageType::New(); mitk::CastToItkImage(mitkMask, mask); mitk::FiberBundle::Pointer newFib = fib->RemoveFibersOutside(mask, removeInside); if (newFib->GetNumFibers()<=0) { QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers."); continue; } node->SetData(newFib); } RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberProcessingView::ExtractWithMask(bool onlyEnds, bool invert) { if (m_MaskImageNode.IsNull()) return; mitk::Image::Pointer mitkMask = dynamic_cast(m_MaskImageNode->GetData()); for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); QString name(node->GetName().c_str()); itkUCharImageType::Pointer mask = itkUCharImageType::New(); mitk::CastToItkImage(mitkMask, mask); mitk::FiberBundle::Pointer newFib = fib->ExtractFiberSubset(mask, !onlyEnds, invert, m_Controls->m_BothEnds->isChecked(), m_Controls->m_FiberExtractionFractionBox->value()); if (newFib->GetNumFibers()<=0) { QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers."); continue; } DataNode::Pointer newNode = DataNode::New(); newNode->SetData(newFib); if (invert) { name += "_not"; if (onlyEnds) name += "-ending-in-mask"; else name += "-passing-mask"; } else { if (onlyEnds) name += "_ending-in-mask"; else name += "_passing-mask"; } newNode->SetName(name.toStdString()); GetDataStorage()->Add(newNode); node->SetVisibility(false); } } void QmitkFiberProcessingView::GenerateRoiImage() { if (m_SelectedPF.empty()) return; mitk::BaseGeometry::Pointer geometry; if (!m_SelectedFB.empty()) { mitk::FiberBundle::Pointer fib = dynamic_cast(m_SelectedFB.front()->GetData()); geometry = fib->GetGeometry(); } else if (m_SelectedImage) geometry = m_SelectedImage->GetGeometry(); else return; itk::Vector spacing = geometry->GetSpacing(); spacing /= m_UpsamplingFactor; mitk::Point3D newOrigin = geometry->GetOrigin(); mitk::Geometry3D::BoundsArrayType bounds = geometry->GetBounds(); newOrigin[0] += bounds.GetElement(0); newOrigin[1] += bounds.GetElement(2); newOrigin[2] += bounds.GetElement(4); itk::Matrix direction; itk::ImageRegion<3> imageRegion; for (int i=0; i<3; i++) for (int j=0; j<3; j++) direction[j][i] = geometry->GetMatrixColumn(i)[j]/spacing[j]; imageRegion.SetSize(0, geometry->GetExtent(0)*m_UpsamplingFactor); imageRegion.SetSize(1, geometry->GetExtent(1)*m_UpsamplingFactor); imageRegion.SetSize(2, geometry->GetExtent(2)*m_UpsamplingFactor); m_PlanarFigureImage = itkUCharImageType::New(); m_PlanarFigureImage->SetSpacing( spacing ); // Set the image spacing m_PlanarFigureImage->SetOrigin( newOrigin ); // Set the image origin m_PlanarFigureImage->SetDirection( direction ); // Set the image direction m_PlanarFigureImage->SetRegions( imageRegion ); m_PlanarFigureImage->Allocate(); m_PlanarFigureImage->FillBuffer( 0 ); Image::Pointer tmpImage = Image::New(); tmpImage->InitializeByItk(m_PlanarFigureImage.GetPointer()); tmpImage->SetVolume(m_PlanarFigureImage->GetBufferPointer()); std::string name = m_SelectedPF.at(0)->GetName(); WritePfToImage(m_SelectedPF.at(0), tmpImage); for (unsigned int i=1; iGetName(); WritePfToImage(m_SelectedPF.at(i), tmpImage); } DataNode::Pointer node = DataNode::New(); tmpImage = Image::New(); tmpImage->InitializeByItk(m_PlanarFigureImage.GetPointer()); tmpImage->SetVolume(m_PlanarFigureImage->GetBufferPointer()); node->SetData(tmpImage); node->SetName(name); this->GetDataStorage()->Add(node); } void QmitkFiberProcessingView::WritePfToImage(mitk::DataNode::Pointer node, mitk::Image* image) { if (dynamic_cast(node->GetData())) { m_PlanarFigure = dynamic_cast(node->GetData()); AccessFixedDimensionByItk_2( image, InternalReorientImagePlane, 3, m_PlanarFigure->GetGeometry(), -1); AccessFixedDimensionByItk_2( m_InternalImage, InternalCalculateMaskFromPlanarFigure, 3, 2, node->GetName() ); } else if (dynamic_cast(node->GetData())) { DataStorage::SetOfObjects::ConstPointer children = GetDataStorage()->GetDerivations(node); for (unsigned int i=0; iSize(); i++) { WritePfToImage(children->at(i), image); } } } template < typename TPixel, unsigned int VImageDimension > void QmitkFiberProcessingView::InternalReorientImagePlane( const itk::Image< TPixel, VImageDimension > *image, mitk::BaseGeometry* planegeo3D, int additionalIndex ) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< float, VImageDimension > FloatImageType; typedef itk::ResampleImageFilter ResamplerType; typename ResamplerType::Pointer resampler = ResamplerType::New(); mitk::PlaneGeometry* planegeo = dynamic_cast(planegeo3D); float upsamp = m_UpsamplingFactor; float gausssigma = 0.5; // Spacing typename ResamplerType::SpacingType spacing = planegeo->GetSpacing(); spacing[0] = image->GetSpacing()[0] / upsamp; spacing[1] = image->GetSpacing()[1] / upsamp; spacing[2] = image->GetSpacing()[2]; resampler->SetOutputSpacing( spacing ); // Size typename ResamplerType::SizeType size; size[0] = planegeo->GetExtentInMM(0) / spacing[0]; size[1] = planegeo->GetExtentInMM(1) / spacing[1]; size[2] = 1; resampler->SetSize( size ); // Origin typename mitk::Point3D orig = planegeo->GetOrigin(); typename mitk::Point3D corrorig; planegeo3D->WorldToIndex(orig,corrorig); corrorig[0] += 0.5/upsamp; corrorig[1] += 0.5/upsamp; corrorig[2] += 0; planegeo3D->IndexToWorld(corrorig,corrorig); resampler->SetOutputOrigin(corrorig ); // Direction typename ResamplerType::DirectionType direction; typename mitk::AffineTransform3D::MatrixType matrix = planegeo->GetIndexToWorldTransform()->GetMatrix(); for(int c=0; cSetOutputDirection( direction ); // Gaussian interpolation if(gausssigma != 0) { double sigma[3]; for( unsigned int d = 0; d < 3; d++ ) sigma[d] = gausssigma * image->GetSpacing()[d]; double alpha = 2.0; typedef itk::GaussianInterpolateImageFunction GaussianInterpolatorType; typename GaussianInterpolatorType::Pointer interpolator = GaussianInterpolatorType::New(); interpolator->SetInputImage( image ); interpolator->SetParameters( sigma, alpha ); resampler->SetInterpolator( interpolator ); } else { typedef typename itk::LinearInterpolateImageFunction InterpolatorType; typename InterpolatorType::Pointer interpolator = InterpolatorType::New(); interpolator->SetInputImage( image ); resampler->SetInterpolator( interpolator ); } resampler->SetInput( image ); resampler->SetDefaultPixelValue(0); resampler->Update(); if(additionalIndex < 0) { this->m_InternalImage = mitk::Image::New(); this->m_InternalImage->InitializeByItk( resampler->GetOutput() ); this->m_InternalImage->SetVolume( resampler->GetOutput()->GetBufferPointer() ); } } template < typename TPixel, unsigned int VImageDimension > void QmitkFiberProcessingView::InternalCalculateMaskFromPlanarFigure( itk::Image< TPixel, VImageDimension > *image, unsigned int axis, std::string ) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::CastImageFilter< ImageType, itkUCharImageType > CastFilterType; // Generate mask image as new image with same header as input image and // initialize with "1". itkUCharImageType::Pointer newMaskImage = itkUCharImageType::New(); newMaskImage->SetSpacing( image->GetSpacing() ); // Set the image spacing newMaskImage->SetOrigin( image->GetOrigin() ); // Set the image origin newMaskImage->SetDirection( image->GetDirection() ); // Set the image direction newMaskImage->SetRegions( image->GetLargestPossibleRegion() ); newMaskImage->Allocate(); newMaskImage->FillBuffer( 1 ); // Generate VTK polygon from (closed) PlanarFigure polyline // (The polyline points are shifted by -0.5 in z-direction to make sure // that the extrusion filter, which afterwards elevates all points by +0.5 // in z-direction, creates a 3D object which is cut by the the plane z=0) const PlaneGeometry *planarFigurePlaneGeometry = m_PlanarFigure->GetPlaneGeometry(); const PlanarFigure::PolyLineType planarFigurePolyline = m_PlanarFigure->GetPolyLine( 0 ); const BaseGeometry *imageGeometry3D = m_InternalImage->GetGeometry( 0 ); vtkPolyData *polyline = vtkPolyData::New(); polyline->Allocate( 1, 1 ); // Determine x- and y-dimensions depending on principal axis int i0, i1; switch ( axis ) { case 0: i0 = 1; i1 = 2; break; case 1: i0 = 0; i1 = 2; break; case 2: default: i0 = 0; i1 = 1; break; } // Create VTK polydata object of polyline contour vtkPoints *points = vtkPoints::New(); PlanarFigure::PolyLineType::const_iterator it; unsigned int numberOfPoints = 0; for ( it = planarFigurePolyline.begin(); it != planarFigurePolyline.end(); ++it ) { Point3D point3D; // Convert 2D point back to the local index coordinates of the selected image Point2D point2D = *it; planarFigurePlaneGeometry->WorldToIndex(point2D, point2D); point2D[0] -= 0.5/m_UpsamplingFactor; point2D[1] -= 0.5/m_UpsamplingFactor; planarFigurePlaneGeometry->IndexToWorld(point2D, point2D); planarFigurePlaneGeometry->Map( point2D, point3D ); // Polygons (partially) outside of the image bounds can not be processed further due to a bug in vtkPolyDataToImageStencil if ( !imageGeometry3D->IsInside( point3D ) ) { float bounds[2] = {0,0}; bounds[0] = this->m_InternalImage->GetLargestPossibleRegion().GetSize().GetElement(i0); bounds[1] = this->m_InternalImage->GetLargestPossibleRegion().GetSize().GetElement(i1); imageGeometry3D->WorldToIndex( point3D, point3D ); if (point3D[i0]<0) point3D[i0] = 0.0; else if (point3D[i0]>bounds[0]) point3D[i0] = bounds[0]-0.001; if (point3D[i1]<0) point3D[i1] = 0.0; else if (point3D[i1]>bounds[1]) point3D[i1] = bounds[1]-0.001; points->InsertNextPoint( point3D[i0], point3D[i1], -0.5 ); numberOfPoints++; } else { imageGeometry3D->WorldToIndex( point3D, point3D ); // Add point to polyline array points->InsertNextPoint( point3D[i0], point3D[i1], -0.5 ); numberOfPoints++; } } polyline->SetPoints( points ); points->Delete(); vtkIdType *ptIds = new vtkIdType[numberOfPoints]; for ( vtkIdType i = 0; i < numberOfPoints; ++i ) ptIds[i] = i; polyline->InsertNextCell( VTK_POLY_LINE, numberOfPoints, ptIds ); // Extrude the generated contour polygon vtkLinearExtrusionFilter *extrudeFilter = vtkLinearExtrusionFilter::New(); extrudeFilter->SetInputData( polyline ); extrudeFilter->SetScaleFactor( 1 ); extrudeFilter->SetExtrusionTypeToNormalExtrusion(); extrudeFilter->SetVector( 0.0, 0.0, 1.0 ); // Make a stencil from the extruded polygon vtkPolyDataToImageStencil *polyDataToImageStencil = vtkPolyDataToImageStencil::New(); polyDataToImageStencil->SetInputConnection( extrudeFilter->GetOutputPort() ); // Export from ITK to VTK (to use a VTK filter) typedef itk::VTKImageImport< itkUCharImageType > ImageImportType; typedef itk::VTKImageExport< itkUCharImageType > ImageExportType; typename ImageExportType::Pointer itkExporter = ImageExportType::New(); itkExporter->SetInput( newMaskImage ); vtkImageImport *vtkImporter = vtkImageImport::New(); this->ConnectPipelines( itkExporter, vtkImporter ); vtkImporter->Update(); // Apply the generated image stencil to the input image vtkImageStencil *imageStencilFilter = vtkImageStencil::New(); imageStencilFilter->SetInputConnection( vtkImporter->GetOutputPort() ); imageStencilFilter->SetStencilConnection(polyDataToImageStencil->GetOutputPort() ); imageStencilFilter->ReverseStencilOff(); imageStencilFilter->SetBackgroundValue( 0 ); imageStencilFilter->Update(); // Export from VTK back to ITK vtkImageExport *vtkExporter = vtkImageExport::New(); vtkExporter->SetInputConnection( imageStencilFilter->GetOutputPort() ); vtkExporter->Update(); typename ImageImportType::Pointer itkImporter = ImageImportType::New(); this->ConnectPipelines( vtkExporter, itkImporter ); itkImporter->Update(); // calculate cropping bounding box m_InternalImageMask3D = itkImporter->GetOutput(); m_InternalImageMask3D->SetDirection(image->GetDirection()); itk::ImageRegionConstIterator itmask(m_InternalImageMask3D, m_InternalImageMask3D->GetLargestPossibleRegion()); itk::ImageRegionIterator itimage(image, image->GetLargestPossibleRegion()); itmask.GoToBegin(); itimage.GoToBegin(); typename ImageType::SizeType lowersize = {{itk::NumericTraits::max(),itk::NumericTraits::max(),itk::NumericTraits::max()}}; typename ImageType::SizeType uppersize = {{0,0,0}}; while( !itmask.IsAtEnd() ) { if(itmask.Get() == 0) itimage.Set(0); else { typename ImageType::IndexType index = itimage.GetIndex(); typename ImageType::SizeType signedindex; signedindex[0] = index[0]; signedindex[1] = index[1]; signedindex[2] = index[2]; lowersize[0] = signedindex[0] < lowersize[0] ? signedindex[0] : lowersize[0]; lowersize[1] = signedindex[1] < lowersize[1] ? signedindex[1] : lowersize[1]; lowersize[2] = signedindex[2] < lowersize[2] ? signedindex[2] : lowersize[2]; uppersize[0] = signedindex[0] > uppersize[0] ? signedindex[0] : uppersize[0]; uppersize[1] = signedindex[1] > uppersize[1] ? signedindex[1] : uppersize[1]; uppersize[2] = signedindex[2] > uppersize[2] ? signedindex[2] : uppersize[2]; } ++itmask; ++itimage; } typename ImageType::IndexType index; index[0] = lowersize[0]; index[1] = lowersize[1]; index[2] = lowersize[2]; typename ImageType::SizeType size; size[0] = uppersize[0] - lowersize[0] + 1; size[1] = uppersize[1] - lowersize[1] + 1; size[2] = uppersize[2] - lowersize[2] + 1; itk::ImageRegion<3> cropRegion = itk::ImageRegion<3>(index, size); // crop internal mask typedef itk::RegionOfInterestImageFilter< itkUCharImageType, itkUCharImageType > ROIMaskFilterType; typename ROIMaskFilterType::Pointer roi2 = ROIMaskFilterType::New(); roi2->SetRegionOfInterest(cropRegion); roi2->SetInput(m_InternalImageMask3D); roi2->Update(); m_InternalImageMask3D = roi2->GetOutput(); Image::Pointer tmpImage = Image::New(); tmpImage->InitializeByItk(m_InternalImageMask3D.GetPointer()); tmpImage->SetVolume(m_InternalImageMask3D->GetBufferPointer()); Image::Pointer tmpImage2 = Image::New(); tmpImage2->InitializeByItk(m_PlanarFigureImage.GetPointer()); const BaseGeometry *pfImageGeometry3D = tmpImage2->GetGeometry( 0 ); const BaseGeometry *intImageGeometry3D = tmpImage->GetGeometry( 0 ); typedef itk::ImageRegionIteratorWithIndex IteratorType; IteratorType imageIterator (m_InternalImageMask3D, m_InternalImageMask3D->GetRequestedRegion()); imageIterator.GoToBegin(); while ( !imageIterator.IsAtEnd() ) { unsigned char val = imageIterator.Value(); if (val>0) { itk::Index<3> index = imageIterator.GetIndex(); Point3D point; point[0] = index[0]; point[1] = index[1]; point[2] = index[2]; intImageGeometry3D->IndexToWorld(point, point); pfImageGeometry3D->WorldToIndex(point, point); point[i0] += 0.5; point[i1] += 0.5; index[0] = point[0]; index[1] = point[1]; index[2] = point[2]; if (pfImageGeometry3D->IsIndexInside(index)) m_PlanarFigureImage->SetPixel(index, 1); } ++imageIterator; } // Clean up VTK objects polyline->Delete(); extrudeFilter->Delete(); polyDataToImageStencil->Delete(); vtkImporter->Delete(); imageStencilFilter->Delete(); //vtkExporter->Delete(); // TODO: crashes when outcommented; memory leak?? delete[] ptIds; } void QmitkFiberProcessingView::UpdateGui() { m_Controls->m_FibLabel->setText("mandatory"); m_Controls->m_PfLabel->setText("needed for extraction"); m_Controls->m_InputData->setTitle("Please Select Input Data"); m_Controls->m_RemoveButton->setEnabled(false); m_Controls->m_PlanarFigureButtonsFrame->setEnabled(false); m_Controls->PFCompoANDButton->setEnabled(false); m_Controls->PFCompoORButton->setEnabled(false); m_Controls->PFCompoNOTButton->setEnabled(false); m_Controls->m_GenerateRoiImage->setEnabled(false); m_Controls->m_ExtractFibersButton->setEnabled(false); m_Controls->m_ModifyButton->setEnabled(false); m_Controls->m_CopyBundle->setEnabled(false); m_Controls->m_JoinBundles->setEnabled(false); m_Controls->m_SubstractBundles->setEnabled(false); // disable alle frames m_Controls->m_BundleWeightFrame->setVisible(false); m_Controls->m_ExtactionFramePF->setVisible(false); m_Controls->m_RemoveDirectionFrame->setVisible(false); m_Controls->m_RemoveLengthFrame->setVisible(false); m_Controls->m_RemoveCurvatureFrame->setVisible(false); m_Controls->m_SmoothFibersFrame->setVisible(false); m_Controls->m_CompressFibersFrame->setVisible(false); m_Controls->m_ColorFibersFrame->setVisible(false); m_Controls->m_MirrorFibersFrame->setVisible(false); m_Controls->m_MaskExtractionFrame->setVisible(false); bool pfSelected = !m_SelectedPF.empty(); bool fibSelected = !m_SelectedFB.empty(); bool multipleFibsSelected = (m_SelectedFB.size()>1); bool maskSelected = m_MaskImageNode.IsNotNull(); bool imageSelected = m_SelectedImage.IsNotNull(); // toggle visibility of elements according to selected method switch ( m_Controls->m_ExtractionMethodBox->currentIndex() ) { case 0: m_Controls->m_ExtactionFramePF->setVisible(true); break; case 1: m_Controls->m_MaskExtractionFrame->setVisible(true); break; } switch ( m_Controls->m_RemovalMethodBox->currentIndex() ) { case 0: m_Controls->m_RemoveDirectionFrame->setVisible(true); if ( fibSelected ) m_Controls->m_RemoveButton->setEnabled(true); break; case 1: m_Controls->m_RemoveLengthFrame->setVisible(true); if ( fibSelected ) m_Controls->m_RemoveButton->setEnabled(true); break; case 2: m_Controls->m_RemoveCurvatureFrame->setVisible(true); if ( fibSelected ) m_Controls->m_RemoveButton->setEnabled(true); break; case 3: break; case 4: break; } switch ( m_Controls->m_ModificationMethodBox->currentIndex() ) { case 0: m_Controls->m_SmoothFibersFrame->setVisible(true); break; case 1: m_Controls->m_SmoothFibersFrame->setVisible(true); break; case 2: m_Controls->m_CompressFibersFrame->setVisible(true); break; case 3: m_Controls->m_ColorFibersFrame->setVisible(true); break; case 4: m_Controls->m_MirrorFibersFrame->setVisible(true); if (m_SelectedSurfaces.size()>0) m_Controls->m_ModifyButton->setEnabled(true); break; case 5: m_Controls->m_BundleWeightFrame->setVisible(true); } // are fiber bundles selected? if ( fibSelected ) { m_Controls->m_CopyBundle->setEnabled(true); m_Controls->m_ModifyButton->setEnabled(true); m_Controls->m_PlanarFigureButtonsFrame->setEnabled(true); m_Controls->m_FibLabel->setText(QString(m_SelectedFB.at(0)->GetName().c_str())); // one bundle and one planar figure needed to extract fibers if (pfSelected && m_Controls->m_ExtractionMethodBox->currentIndex()==0) { m_Controls->m_InputData->setTitle("Input Data"); m_Controls->m_PfLabel->setText(QString(m_SelectedPF.at(0)->GetName().c_str())); m_Controls->m_ExtractFibersButton->setEnabled(true); } // more than two bundles needed to join/subtract if (multipleFibsSelected) { m_Controls->m_FibLabel->setText("multiple bundles selected"); m_Controls->m_JoinBundles->setEnabled(true); m_Controls->m_SubstractBundles->setEnabled(true); } if (maskSelected && m_Controls->m_ExtractionMethodBox->currentIndex()==1) { m_Controls->m_InputData->setTitle("Input Data"); m_Controls->m_PfLabel->setText(QString(m_MaskImageNode->GetName().c_str())); m_Controls->m_ExtractFibersButton->setEnabled(true); } if (maskSelected && (m_Controls->m_RemovalMethodBox->currentIndex()==3 || m_Controls->m_RemovalMethodBox->currentIndex()==4) ) { m_Controls->m_InputData->setTitle("Input Data"); m_Controls->m_PfLabel->setText(QString(m_MaskImageNode->GetName().c_str())); m_Controls->m_RemoveButton->setEnabled(true); } } // are planar figures selected? if (pfSelected) { if ( fibSelected || m_SelectedImage.IsNotNull()) m_Controls->m_GenerateRoiImage->setEnabled(true); if (m_SelectedPF.size() > 1) { m_Controls->PFCompoANDButton->setEnabled(true); m_Controls->PFCompoORButton->setEnabled(true); } else m_Controls->PFCompoNOTButton->setEnabled(true); } // is image selected if (imageSelected || maskSelected) { m_Controls->m_PlanarFigureButtonsFrame->setEnabled(true); } } void QmitkFiberProcessingView::NodeRemoved(const mitk::DataNode* node) { berry::IWorkbenchPart::Pointer nullPart; QList nodes; OnSelectionChanged(nullPart, nodes); } void QmitkFiberProcessingView::NodeAdded(const mitk::DataNode* node) { berry::IWorkbenchPart::Pointer nullPart; QList nodes; OnSelectionChanged(nullPart, nodes); } void QmitkFiberProcessingView::OnSelectionChanged(berry::IWorkbenchPart::Pointer /*part*/, const QList& nodes) { //reset existing Vectors containing FiberBundles and PlanarFigures from a previous selection m_SelectedFB.clear(); m_SelectedPF.clear(); m_SelectedSurfaces.clear(); m_SelectedImage = nullptr; m_MaskImageNode = nullptr; for (auto node: nodes) { if ( dynamic_cast(node->GetData()) ) m_SelectedFB.push_back(node); else if (dynamic_cast(node->GetData()) || dynamic_cast(node->GetData()) || dynamic_cast(node->GetData())) m_SelectedPF.push_back(node); else if (dynamic_cast(node->GetData())) { m_SelectedImage = dynamic_cast(node->GetData()); bool isBinary = false; node->GetPropertyValue("binary", isBinary); if (isBinary) m_MaskImageNode = node; } else if (dynamic_cast(node->GetData())) m_SelectedSurfaces.push_back(dynamic_cast(node->GetData())); } if (m_SelectedFB.empty() && m_SelectedSurfaces.empty()) { int maxLayer = 0; itk::VectorContainer::ConstPointer nodes = this->GetDataStorage()->GetAll(); for (unsigned int i=0; iSize(); i++) if (dynamic_cast(nodes->at(i)->GetData())) { mitk::DataStorage::SetOfObjects::ConstPointer sources = GetDataStorage()->GetSources(nodes->at(i)); if (sources->Size()>0) continue; int layer = 0; nodes->at(i)->GetPropertyValue("layer", layer); if (layer>=maxLayer) { maxLayer = layer; m_SelectedFB.clear(); m_SelectedFB.push_back(nodes->at(i)); } } } if (m_SelectedPF.empty()) { int maxLayer = 0; itk::VectorContainer::ConstPointer nodes = this->GetDataStorage()->GetAll(); for (unsigned int i=0; iSize(); i++) if (dynamic_cast(nodes->at(i)->GetData()) || dynamic_cast(nodes->at(i)->GetData()) || dynamic_cast(nodes->at(i)->GetData())) { mitk::DataStorage::SetOfObjects::ConstPointer sources = GetDataStorage()->GetSources(nodes->at(i)); if (sources->Size()>0) continue; int layer = 0; nodes->at(i)->GetPropertyValue("layer", layer); if (layer>=maxLayer) { maxLayer = layer; m_SelectedPF.clear(); m_SelectedPF.push_back(nodes->at(i)); } } } UpdateGui(); } void QmitkFiberProcessingView::OnDrawPolygon() { mitk::PlanarPolygon::Pointer figure = mitk::PlanarPolygon::New(); figure->ClosedOn(); this->AddFigureToDataStorage(figure, QString("Polygon%1").arg(++m_PolygonCounter)); } void QmitkFiberProcessingView::OnDrawCircle() { mitk::PlanarCircle::Pointer figure = mitk::PlanarCircle::New(); this->AddFigureToDataStorage(figure, QString("Circle%1").arg(++m_CircleCounter)); } void QmitkFiberProcessingView::AddFigureToDataStorage(mitk::PlanarFigure* figure, const QString& name, const char *, mitk::BaseProperty* ) { // initialize figure's geometry with empty geometry mitk::PlaneGeometry::Pointer emptygeometry = mitk::PlaneGeometry::New(); figure->SetPlaneGeometry( emptygeometry ); //set desired data to DataNode where Planarfigure is stored mitk::DataNode::Pointer newNode = mitk::DataNode::New(); newNode->SetName(name.toStdString()); newNode->SetData(figure); newNode->SetBoolProperty("planarfigure.3drendering", true); newNode->SetBoolProperty("planarfigure.3drendering.fill", true); mitk::PlanarFigureInteractor::Pointer figureInteractor = dynamic_cast(newNode->GetDataInteractor().GetPointer()); if(figureInteractor.IsNull()) { figureInteractor = mitk::PlanarFigureInteractor::New(); us::Module* planarFigureModule = us::ModuleRegistry::GetModule( "MitkPlanarFigure" ); figureInteractor->LoadStateMachine("PlanarFigureInteraction.xml", planarFigureModule ); figureInteractor->SetEventConfig( "PlanarFigureConfig.xml", planarFigureModule ); figureInteractor->SetDataNode(newNode); } // figure drawn on the topmost layer / image GetDataStorage()->Add(newNode ); for(unsigned int i = 0; i < m_SelectedPF.size(); i++) m_SelectedPF[i]->SetSelected(false); newNode->SetSelected(true); m_SelectedPF.clear(); m_SelectedPF.push_back(newNode); UpdateGui(); } void QmitkFiberProcessingView::ExtractWithPlanarFigure() { if ( m_SelectedFB.empty() || m_SelectedPF.empty() ){ QMessageBox::information( nullptr, "Warning", "No fibe bundle selected!"); return; } std::vector fiberBundles = m_SelectedFB; mitk::DataNode::Pointer planarFigure = m_SelectedPF.at(0); for (unsigned int i=0; i(fiberBundles.at(i)->GetData()); mitk::FiberBundle::Pointer extFB = fib->ExtractFiberSubset(planarFigure, GetDataStorage()); if (extFB->GetNumFibers()<=0) { QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers."); continue; } mitk::DataNode::Pointer node; node = mitk::DataNode::New(); node->SetData(extFB); QString name(fiberBundles.at(i)->GetName().c_str()); name += "*"; //name += planarFigure->GetName().c_str(); node->SetName(name.toStdString()); fiberBundles.at(i)->SetVisibility(false); GetDataStorage()->Add(node); } } void QmitkFiberProcessingView::GenerateAndComposite() { mitk::PlanarFigureComposite::Pointer PFCAnd = mitk::PlanarFigureComposite::New(); PFCAnd->setOperationType(mitk::PlanarFigureComposite::AND); mitk::DataNode::Pointer newPFCNode; newPFCNode = mitk::DataNode::New(); newPFCNode->SetName("AND"); newPFCNode->SetData(PFCAnd); AddCompositeToDatastorage(newPFCNode, m_SelectedPF); m_SelectedPF.clear(); m_SelectedPF.push_back(newPFCNode); UpdateGui(); } void QmitkFiberProcessingView::GenerateOrComposite() { mitk::PlanarFigureComposite::Pointer PFCOr = mitk::PlanarFigureComposite::New(); PFCOr->setOperationType(mitk::PlanarFigureComposite::OR); mitk::DataNode::Pointer newPFCNode; newPFCNode = mitk::DataNode::New(); newPFCNode->SetName("OR"); newPFCNode->SetData(PFCOr); AddCompositeToDatastorage(newPFCNode, m_SelectedPF); m_SelectedPF.clear(); m_SelectedPF.push_back(newPFCNode); UpdateGui(); } void QmitkFiberProcessingView::GenerateNotComposite() { mitk::PlanarFigureComposite::Pointer PFCNot = mitk::PlanarFigureComposite::New(); PFCNot->setOperationType(mitk::PlanarFigureComposite::NOT); mitk::DataNode::Pointer newPFCNode; newPFCNode = mitk::DataNode::New(); newPFCNode->SetName("NOT"); newPFCNode->SetData(PFCNot); AddCompositeToDatastorage(newPFCNode, m_SelectedPF); m_SelectedPF.clear(); m_SelectedPF.push_back(newPFCNode); UpdateGui(); } void QmitkFiberProcessingView::AddCompositeToDatastorage(mitk::DataNode::Pointer pfc, std::vector children, mitk::DataNode::Pointer parentNode ) { pfc->SetSelected(true); if (parentNode.IsNotNull()) GetDataStorage()->Add(pfc, parentNode); else GetDataStorage()->Add(pfc); for (auto child : children) { if (dynamic_cast(child->GetData())) { mitk::DataNode::Pointer newChild; newChild = mitk::DataNode::New(); newChild->SetData(dynamic_cast(child->GetData())); newChild->SetName( child->GetName() ); newChild->SetBoolProperty("planarfigure.3drendering", true); newChild->SetBoolProperty("planarfigure.3drendering.fill", true); GetDataStorage()->Add(newChild, pfc); GetDataStorage()->Remove(child); } else if (dynamic_cast(child->GetData())) { mitk::DataNode::Pointer newChild; newChild = mitk::DataNode::New(); newChild->SetData(dynamic_cast(child->GetData())); newChild->SetName( child->GetName() ); std::vector< mitk::DataNode::Pointer > grandChildVector; mitk::DataStorage::SetOfObjects::ConstPointer grandchildren = GetDataStorage()->GetDerivations(child); for( mitk::DataStorage::SetOfObjects::const_iterator it = grandchildren->begin(); it != grandchildren->end(); ++it ) grandChildVector.push_back(*it); AddCompositeToDatastorage(newChild, grandChildVector, pfc); GetDataStorage()->Remove(child); } } UpdateGui(); } void QmitkFiberProcessingView::CopyBundles() { if ( m_SelectedFB.empty() ){ QMessageBox::information( nullptr, "Warning", "Select at least one fiber bundle!"); MITK_WARN("QmitkFiberProcessingView") << "Select at least one fiber bundle!"; return; } for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); mitk::FiberBundle::Pointer newFib = fib->GetDeepCopy(); node->SetVisibility(false); QString name(""); name += QString(m_SelectedFB.at(0)->GetName().c_str()); name += "_copy"; mitk::DataNode::Pointer fbNode = mitk::DataNode::New(); fbNode->SetData(newFib); fbNode->SetName(name.toStdString()); fbNode->SetVisibility(true); GetDataStorage()->Add(fbNode); } UpdateGui(); } void QmitkFiberProcessingView::JoinBundles() { if ( m_SelectedFB.size()<2 ){ QMessageBox::information( nullptr, "Warning", "Select at least two fiber bundles!"); MITK_WARN("QmitkFiberProcessingView") << "Select at least two fiber bundles!"; return; } mitk::FiberBundle::Pointer newBundle = dynamic_cast(m_SelectedFB.at(0)->GetData()); m_SelectedFB.at(0)->SetVisibility(false); QString name(""); name += QString(m_SelectedFB.at(0)->GetName().c_str()); for (unsigned int i=1; iAddBundle(dynamic_cast(m_SelectedFB.at(i)->GetData())); name += "+"+QString(m_SelectedFB.at(i)->GetName().c_str()); m_SelectedFB.at(i)->SetVisibility(false); } mitk::DataNode::Pointer fbNode = mitk::DataNode::New(); fbNode->SetData(newBundle); fbNode->SetName(name.toStdString()); fbNode->SetVisibility(true); GetDataStorage()->Add(fbNode); UpdateGui(); } void QmitkFiberProcessingView::SubstractBundles() { if ( m_SelectedFB.size()<2 ){ QMessageBox::information( nullptr, "Warning", "Select at least two fiber bundles!"); MITK_WARN("QmitkFiberProcessingView") << "Select at least two fiber bundles!"; return; } mitk::FiberBundle::Pointer newBundle = dynamic_cast(m_SelectedFB.at(0)->GetData()); m_SelectedFB.at(0)->SetVisibility(false); QString name(""); name += QString(m_SelectedFB.at(0)->GetName().c_str()); for (unsigned int i=1; iSubtractBundle(dynamic_cast(m_SelectedFB.at(i)->GetData())); if (newBundle.IsNull()) break; name += "-"+QString(m_SelectedFB.at(i)->GetName().c_str()); m_SelectedFB.at(i)->SetVisibility(false); } if (newBundle.IsNull()) { QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers. Did you select the fiber bundles in the correct order? X-Y is not equal to Y-X!"); return; } mitk::DataNode::Pointer fbNode = mitk::DataNode::New(); fbNode->SetData(newBundle); fbNode->SetName(name.toStdString()); fbNode->SetVisibility(true); GetDataStorage()->Add(fbNode); UpdateGui(); } void QmitkFiberProcessingView::ResampleSelectedBundlesSpline() { double factor = this->m_Controls->m_SmoothFibersBox->value(); for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); fib->ResampleSpline(factor); } RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberProcessingView::ResampleSelectedBundlesLinear() { double factor = this->m_Controls->m_SmoothFibersBox->value(); for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); fib->ResampleLinear(factor); } RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberProcessingView::CompressSelectedBundles() { double factor = this->m_Controls->m_ErrorThresholdBox->value(); for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); fib->Compress(factor); fib->ColorFibersByOrientation(); } RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberProcessingView::DoImageColorCoding() { if (m_Controls->m_ColorMapBox->GetSelectedNode().IsNull()) { QMessageBox::information(nullptr, "Bundle coloring aborted:", "No image providing the scalar values for coloring the selected bundle available."); return; } for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); fib->ColorFibersByScalarMap(dynamic_cast(m_Controls->m_ColorMapBox->GetSelectedNode()->GetData()), m_Controls->m_FiberOpacityBox->isChecked()); } if (auto renderWindowPart = this->GetRenderWindowPart()) { renderWindowPart->RequestUpdate(); } } void QmitkFiberProcessingView::DoCurvatureColorCoding() { for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); fib->ColorFibersByCurvature(); } if (auto renderWindowPart = this->GetRenderWindowPart()) { renderWindowPart->RequestUpdate(); } } +void QmitkFiberProcessingView::DoWeightColorCoding() +{ + for (auto node : m_SelectedFB) + { + mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); + fib->ColorFibersByFiberWeights(); + } + + if (auto renderWindowPart = this->GetRenderWindowPart()) + { + renderWindowPart->RequestUpdate(); + } +} + void QmitkFiberProcessingView::MirrorFibers() { unsigned int axis = this->m_Controls->m_MirrorFibersBox->currentIndex(); for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); if (m_SelectedImage.IsNotNull()) fib->SetReferenceGeometry(m_SelectedImage->GetGeometry()); fib->MirrorFibers(axis); } for (auto surf : m_SelectedSurfaces) { vtkSmartPointer poly = surf->GetVtkPolyData(); vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); for (int i=0; iGetNumberOfPoints(); i++) { double* point = poly->GetPoint(i); point[axis] *= -1; vtkNewPoints->InsertNextPoint(point); } poly->SetPoints(vtkNewPoints); surf->CalculateBoundingBox(); } if (auto renderWindowPart = this->GetRenderWindowPart()) { renderWindowPart->RequestUpdate(); } } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingView.h b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingView.h index 970dd2c34f..7c83fc44ae 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingView.h +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingView.h @@ -1,190 +1,191 @@ /*=================================================================== 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 QmitkFiberProcessingView_h #define QmitkFiberProcessingView_h #include #include "ui_QmitkFiberProcessingViewControls.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include /*! \brief View to process fiber bundles. Supplies methods to extract fibers from the bundle, fiber resampling, mirroring, join and subtract bundles and much more. */ class QmitkFiberProcessingView : 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) Q_OBJECT public: typedef itk::Image< unsigned char, 3 > itkUCharImageType; static const std::string VIEW_ID; QmitkFiberProcessingView(); virtual ~QmitkFiberProcessingView(); virtual void CreateQtPartControl(QWidget *parent) override; /// /// Sets the focus to an internal widget. /// virtual void SetFocus() override; protected slots: void OnDrawCircle(); ///< add circle interactors etc. void OnDrawPolygon(); ///< add circle interactors etc. void GenerateAndComposite(); void GenerateOrComposite(); void GenerateNotComposite(); void CopyBundles(); ///< add copies of selected bundles to data storage void JoinBundles(); ///< merge selected fiber bundles void SubstractBundles(); ///< subtract bundle A from bundle B. Not commutative! Defined by order of selection. void GenerateRoiImage(); ///< generate binary image of selected planar figures. void Remove(); void Extract(); void Modify(); void UpdateGui(); ///< update button activity etc. dpending on current datamanager selection void OnMaskExtractionChanged(); virtual void AddFigureToDataStorage(mitk::PlanarFigure* figure, const QString& name, const char *propertyKey = nullptr, mitk::BaseProperty *property = nullptr ); protected: void MirrorFibers(); ///< mirror bundle on the specified plane void ResampleSelectedBundlesSpline(); ///< void ResampleSelectedBundlesLinear(); ///< void DoImageColorCoding(); ///< color fibers by selected scalar image + void DoWeightColorCoding(); ///< color fibers by their respective weights void DoCurvatureColorCoding(); ///< color fibers by curvature void CompressSelectedBundles(); ///< remove points below certain error threshold void WeightFibers(); void RemoveWithMask(bool removeInside); void RemoveDir(); void ApplyCurvatureThreshold(); ///< remove/split fibers with a too high curvature threshold void PruneBundle(); ///< remove too short/too long fibers void ExtractWithMask(bool onlyEnds, bool invert); void ExtractWithPlanarFigure(); /// \brief called by QmitkAbstractView when DataManager's selection has changed virtual void OnSelectionChanged(berry::IWorkbenchPart::Pointer part, const QList& nodes) override; Ui::QmitkFiberProcessingViewControls* m_Controls; /** Connection from VTK to ITK */ template void ConnectPipelines(VTK_Exporter* exporter, ITK_Importer importer) { importer->SetUpdateInformationCallback(exporter->GetUpdateInformationCallback()); importer->SetPipelineModifiedCallback(exporter->GetPipelineModifiedCallback()); importer->SetWholeExtentCallback(exporter->GetWholeExtentCallback()); importer->SetSpacingCallback(exporter->GetSpacingCallback()); importer->SetOriginCallback(exporter->GetOriginCallback()); importer->SetScalarTypeCallback(exporter->GetScalarTypeCallback()); importer->SetNumberOfComponentsCallback(exporter->GetNumberOfComponentsCallback()); importer->SetPropagateUpdateExtentCallback(exporter->GetPropagateUpdateExtentCallback()); importer->SetUpdateDataCallback(exporter->GetUpdateDataCallback()); importer->SetDataExtentCallback(exporter->GetDataExtentCallback()); importer->SetBufferPointerCallback(exporter->GetBufferPointerCallback()); importer->SetCallbackUserData(exporter->GetCallbackUserData()); } template void ConnectPipelines(ITK_Exporter exporter, VTK_Importer* importer) { importer->SetUpdateInformationCallback(exporter->GetUpdateInformationCallback()); importer->SetPipelineModifiedCallback(exporter->GetPipelineModifiedCallback()); importer->SetWholeExtentCallback(exporter->GetWholeExtentCallback()); importer->SetSpacingCallback(exporter->GetSpacingCallback()); importer->SetOriginCallback(exporter->GetOriginCallback()); importer->SetScalarTypeCallback(exporter->GetScalarTypeCallback()); importer->SetNumberOfComponentsCallback(exporter->GetNumberOfComponentsCallback()); importer->SetPropagateUpdateExtentCallback(exporter->GetPropagateUpdateExtentCallback()); importer->SetUpdateDataCallback(exporter->GetUpdateDataCallback()); importer->SetDataExtentCallback(exporter->GetDataExtentCallback()); importer->SetBufferPointerCallback(exporter->GetBufferPointerCallback()); importer->SetCallbackUserData(exporter->GetCallbackUserData()); } template < typename TPixel, unsigned int VImageDimension > void InternalCalculateMaskFromPlanarFigure( itk::Image< TPixel, VImageDimension > *image, unsigned int axis, std::string nodeName ); template < typename TPixel, unsigned int VImageDimension > void InternalReorientImagePlane( const itk::Image< TPixel, VImageDimension > *image, mitk::BaseGeometry* planegeo3D, int additionalIndex ); int m_CircleCounter; ///< used for data node naming int m_PolygonCounter; ///< used for data node naming std::vector m_SelectedFB; ///< selected fiber bundle nodes std::vector m_SelectedPF; ///< selected planar figure nodes std::vector m_SelectedSurfaces; mitk::Image::Pointer m_SelectedImage; mitk::Image::Pointer m_InternalImage; mitk::PlanarFigure::Pointer m_PlanarFigure; itkUCharImageType::Pointer m_InternalImageMask3D; itkUCharImageType::Pointer m_PlanarFigureImage; float m_UpsamplingFactor; ///< upsampling factor for all image generations mitk::DataNode::Pointer m_MaskImageNode; void AddCompositeToDatastorage(mitk::DataNode::Pointer pfc, std::vector children, mitk::DataNode::Pointer parentNode=nullptr); void debugPFComposition(mitk::PlanarFigureComposite::Pointer , int ); void WritePfToImage(mitk::DataNode::Pointer node, mitk::Image* image); mitk::DataNode::Pointer GenerateTractDensityImage(mitk::FiberBundle::Pointer fib, bool binary, bool absolute); mitk::DataNode::Pointer GenerateColorHeatmap(mitk::FiberBundle::Pointer fib); mitk::DataNode::Pointer GenerateFiberEndingsImage(mitk::FiberBundle::Pointer fib); mitk::DataNode::Pointer GenerateFiberEndingsPointSet(mitk::FiberBundle::Pointer fib); void NodeAdded( const mitk::DataNode* node ) override; void NodeRemoved(const mitk::DataNode* node) override; }; #endif // _QMITKFIBERTRACKINGVIEW_H_INCLUDED diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingViewControls.ui b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingViewControls.ui index c4e0f9faf5..ff22e91380 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingViewControls.ui +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingViewControls.ui @@ -1,1422 +1,1427 @@ QmitkFiberProcessingViewControls 0 0 386 540 Form - 1 + 0 0 0 - 355 - 363 + 354 + 347 Fiber Extraction Extract a fiber subset from the selected fiber bundle using manually placed planar figures as waypoints or binary regions of interest. false 0 0 200 16777215 11 Extract fibers passing through selected ROI or composite ROI. Select ROI and fiber bundle to execute. Extract Qt::Vertical 20 40 QFrame::NoFrame QFrame::Raised 9 9 9 9 0 false 0 0 16777215 16777215 11 Generate a binary image containing all selected ROIs. Select at least one ROI (planar figure) and a reference fiber bundle or image. Generate ROI Image 0 0 200 0 16777215 60 QFrame::NoFrame QFrame::Raised 0 0 0 0 false 60 16777215 Create OR composition with selected ROIs. OR Qt::Horizontal 40 20 false 60 16777215 Create NOT composition from selected ROI. NOT false 60 16777215 Create AND composition with selected ROIs. AND 0 0 200 0 16777215 60 QFrame::NoFrame QFrame::Raised 0 0 0 0 30 30 Draw circular ROI. Select reference fiber bundle to execute. :/QmitkDiffusionImaging/circle.png:/QmitkDiffusionImaging/circle.png 32 32 false true Qt::Horizontal 40 20 30 30 Draw polygonal ROI. Select reference fiber bundle to execute. :/QmitkDiffusionImaging/polygon.png:/QmitkDiffusionImaging/polygon.png 32 32 true true 0 0 Extract using planar figures Extract using binary ROI image QFrame::NoFrame QFrame::Raised 0 0 0 0 Extract fibers: Min. overlap: 0 0 Ending in mask Not ending in mask Passing mask Not passing mask Both ends true 1.000000000000000 0.100000000000000 0 0 - 355 - 380 + 354 + 342 Fiber Removal Remove fibers that satisfy certain criteria from the selected bundle. 0 0 Remove fibers in direction Remove fibers by length Remove fibers by curvature Remove fiber parts outside mask Remove fiber parts inside mask QFrame::NoFrame QFrame::Raised 0 0 0 0 Qt::Horizontal 40 20 Minimum fiber length in mm 0 999999999 20 Max. Length: Min. Length: Maximum fiber length in mm 0 999999999 300 QFrame::NoFrame QFrame::Raised 0 0 0 0 0 X: Y: Z: Angle: Angular deviation threshold in degree 1 90.000000000000000 1.000000000000000 25.000000000000000 QFrame::NoFrame QFrame::Raised 0 0 0 0 If unchecked, the fiber exceeding the threshold will be split in two instead of removed. Remove Fiber false QFrame::NoFrame QFrame::Raised 0 0 0 0 0 Max. Angular Deviation: Qt::Horizontal 40 20 Maximum angular deviation in degree 180.000000000000000 0.100000000000000 30.000000000000000 Distance: Distance in mm 1 999.000000000000000 1.000000000000000 10.000000000000000 false 0 0 200 16777215 11 Remove Qt::Vertical 20 40 0 0 - 280 - 349 + 368 + 309 Bundle Modification Modify the selected bundle with operations such as fiber resampling, FA coloring, etc. QFrame::NoFrame QFrame::Raised 0 0 0 0 0 Error threshold in mm: 999999999.000000000000000 0.100000000000000 0.100000000000000 QFrame::NoFrame QFrame::Raised 0 0 0 0 0 Sagittal Coronal Axial Select direction: QFrame::NoFrame QFrame::Raised 0 0 0 0 0 Scalar map: If checked, the image values are not only used to color the fibers but are also used as opaxity values. Values as opacity true 0 0 Resample fibers (spline) Resample fibers (linear) Compress fibers Color fibers by scalar map (e.g. FA) Mirror fibers - Weight Bundle + Weight bundle Color fibers by curvature + + + Color fibers by fiber weights + + QFrame::NoFrame QFrame::Raised 0 0 0 0 0 0.010000000000000 999999999.000000000000000 0.100000000000000 1.000000000000000 Point distance in mm: Qt::Vertical 20 40 false 0 0 200 16777215 11 Execute QFrame::NoFrame QFrame::Raised 0 0 0 0 0 Weight: 7 999999999.000000000000000 0.100000000000000 1.000000000000000 0 0 - 360 - 104 + 368 + 152 Bundle Operations Join, subtract or copy bundles. false 0 0 200 16777215 11 Returns all fibers contained in bundle X that are not contained in bundle Y (not commutative!). Select at least two fiber bundles to execute. Substract Qt::Vertical 20 40 false 0 0 200 16777215 11 Merge selected fiber bundles. Select at least two fiber bundles to execute. Join false 0 0 200 16777215 11 Merge selected fiber bundles. Select at least two fiber bundles to execute. Copy Please Select Input Data <html><head/><body><p><span style=" color:#ff0000;">mandatory</span></p></body></html> true <html><head/><body><p><span style=" color:#969696;">needed for extraction</span></p></body></html> true Input DTI Fiber Bundle: Binary seed ROI. If not specified, the whole image area is seeded. ROI: Qt::Vertical 20 40 QmitkDataStorageComboBox QComboBox
QmitkDataStorageComboBox.h
 diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberQuantificationView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberQuantificationView.cpp index 7b0bf6bbe0..fd48eb2f8c 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberQuantificationView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberQuantificationView.cpp @@ -1,475 +1,475 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ // Blueberry #include #include // Qmitk #include "QmitkFiberQuantificationView.h" // Qt #include // MITK #include #include #include #include #include #include #include #include #include #include // ITK #include #include #include #include #include #include #include #include #include #include const std::string QmitkFiberQuantificationView::VIEW_ID = "org.mitk.views.fiberquantification"; const std::string id_DataManager = "org.mitk.views.datamanager"; using namespace mitk; QmitkFiberQuantificationView::QmitkFiberQuantificationView() - : QmitkAbstractView() - , m_Controls( 0 ) - , m_UpsamplingFactor(5) + : QmitkAbstractView() + , m_Controls( 0 ) + , m_UpsamplingFactor(5) { } // Destructor QmitkFiberQuantificationView::~QmitkFiberQuantificationView() { } void QmitkFiberQuantificationView::CreateQtPartControl( QWidget *parent ) { - // build up qt view, unless already done - if ( !m_Controls ) - { - // create GUI widgets from the Qt Designer's .ui file - m_Controls = new Ui::QmitkFiberQuantificationViewControls; - m_Controls->setupUi( parent ); - - connect( m_Controls->m_ProcessFiberBundleButton, SIGNAL(clicked()), this, SLOT(ProcessSelectedBundles()) ); - connect( m_Controls->m_ExtractFiberPeaks, SIGNAL(clicked()), this, SLOT(CalculateFiberDirections()) ); - } + // build up qt view, unless already done + if ( !m_Controls ) + { + // create GUI widgets from the Qt Designer's .ui file + m_Controls = new Ui::QmitkFiberQuantificationViewControls; + m_Controls->setupUi( parent ); + + connect( m_Controls->m_ProcessFiberBundleButton, SIGNAL(clicked()), this, SLOT(ProcessSelectedBundles()) ); + connect( m_Controls->m_ExtractFiberPeaks, SIGNAL(clicked()), this, SLOT(CalculateFiberDirections()) ); + } } void QmitkFiberQuantificationView::SetFocus() { m_Controls->m_ProcessFiberBundleButton->setFocus(); } void QmitkFiberQuantificationView::CalculateFiberDirections() { - typedef itk::Image ItkUcharImgType; - typedef itk::Image< itk::Vector< float, 3>, 3 > ItkDirectionImage3DType; - typedef itk::VectorContainer< unsigned int, ItkDirectionImage3DType::Pointer > ItkDirectionImageContainerType; - - // load fiber bundle - mitk::FiberBundle::Pointer inputTractogram = dynamic_cast(m_SelectedFB.back()->GetData()); - - itk::TractsToVectorImageFilter::Pointer fOdfFilter = itk::TractsToVectorImageFilter::New(); + typedef itk::Image ItkUcharImgType; + typedef itk::Image< itk::Vector< float, 3>, 3 > ItkDirectionImage3DType; + typedef itk::VectorContainer< unsigned int, ItkDirectionImage3DType::Pointer > ItkDirectionImageContainerType; + + // load fiber bundle + mitk::FiberBundle::Pointer inputTractogram = dynamic_cast(m_SelectedFB.back()->GetData()); + + itk::TractsToVectorImageFilter::Pointer fOdfFilter = itk::TractsToVectorImageFilter::New(); + if (m_SelectedImage.IsNotNull()) + { + ItkUcharImgType::Pointer itkMaskImage = ItkUcharImgType::New(); + mitk::CastToItkImage(m_SelectedImage, itkMaskImage); + fOdfFilter->SetMaskImage(itkMaskImage); + } + + // extract directions from fiber bundle + fOdfFilter->SetFiberBundle(inputTractogram); + fOdfFilter->SetAngularThreshold(cos(m_Controls->m_AngularThreshold->value()*M_PI/180)); + fOdfFilter->SetNormalizeVectors(m_Controls->m_NormalizeDirectionsBox->isChecked()); + fOdfFilter->SetUseWorkingCopy(true); + fOdfFilter->SetCreateDirectionImages(m_Controls->m_DirectionImagesBox->isChecked()); + fOdfFilter->SetSizeThreshold(m_Controls->m_PeakThreshold->value()); + fOdfFilter->SetMaxNumDirections(m_Controls->m_MaxNumDirections->value()); + fOdfFilter->Update(); + + QString name = m_SelectedFB.back()->GetName().c_str(); + + if (m_Controls->m_VectorFieldBox->isChecked()) + { + float minSpacing = 1; if (m_SelectedImage.IsNotNull()) { - ItkUcharImgType::Pointer itkMaskImage = ItkUcharImgType::New(); - mitk::CastToItkImage(m_SelectedImage, itkMaskImage); - fOdfFilter->SetMaskImage(itkMaskImage); + mitk::Vector3D outImageSpacing = m_SelectedImage->GetGeometry()->GetSpacing(); + + if(outImageSpacing[0]SetFiberBundle(inputTractogram); - fOdfFilter->SetAngularThreshold(cos(m_Controls->m_AngularThreshold->value()*M_PI/180)); - fOdfFilter->SetNormalizeVectors(m_Controls->m_NormalizeDirectionsBox->isChecked()); - fOdfFilter->SetUseWorkingCopy(true); - fOdfFilter->SetCreateDirectionImages(m_Controls->m_DirectionImagesBox->isChecked()); - fOdfFilter->SetSizeThreshold(m_Controls->m_PeakThreshold->value()); - fOdfFilter->SetMaxNumDirections(m_Controls->m_MaxNumDirections->value()); - fOdfFilter->Update(); + mitk::FiberBundle::Pointer directions = fOdfFilter->GetOutputFiberBundle(); + mitk::DataNode::Pointer node = mitk::DataNode::New(); + node->SetData(directions); + node->SetName((name+"_VECTOR_FIELD").toStdString().c_str()); + node->SetProperty("Fiber2DSliceThickness", mitk::FloatProperty::New(minSpacing)); + node->SetProperty("Fiber2DfadeEFX", mitk::BoolProperty::New(false)); + node->SetProperty("color", mitk::ColorProperty::New(1.0f, 1.0f, 1.0f)); - QString name = m_SelectedFB.back()->GetName().c_str(); + GetDataStorage()->Add(node, m_SelectedFB.back()); + } - if (m_Controls->m_VectorFieldBox->isChecked()) - { - float minSpacing = 1; - if (m_SelectedImage.IsNotNull()) - { - mitk::Vector3D outImageSpacing = m_SelectedImage->GetGeometry()->GetSpacing(); - - if(outImageSpacing[0]GetOutputFiberBundle(); - mitk::DataNode::Pointer node = mitk::DataNode::New(); - node->SetData(directions); - node->SetName((name+"_VECTOR_FIELD").toStdString().c_str()); - node->SetProperty("Fiber2DSliceThickness", mitk::FloatProperty::New(minSpacing)); - node->SetProperty("Fiber2DfadeEFX", mitk::BoolProperty::New(false)); - node->SetProperty("color", mitk::ColorProperty::New(1.0f, 1.0f, 1.0f)); - - GetDataStorage()->Add(node, m_SelectedFB.back()); - } + if (m_Controls->m_NumDirectionsBox->isChecked()) + { + mitk::Image::Pointer mitkImage = mitk::Image::New(); + mitkImage->InitializeByItk( fOdfFilter->GetNumDirectionsImage().GetPointer() ); + mitkImage->SetVolume( fOdfFilter->GetNumDirectionsImage()->GetBufferPointer() ); - if (m_Controls->m_NumDirectionsBox->isChecked()) - { - mitk::Image::Pointer mitkImage = mitk::Image::New(); - mitkImage->InitializeByItk( fOdfFilter->GetNumDirectionsImage().GetPointer() ); - mitkImage->SetVolume( fOdfFilter->GetNumDirectionsImage()->GetBufferPointer() ); - - mitk::DataNode::Pointer node = mitk::DataNode::New(); - node->SetData(mitkImage); - node->SetName((name+"_NUM_DIRECTIONS").toStdString().c_str()); - GetDataStorage()->Add(node, m_SelectedFB.back()); - } + mitk::DataNode::Pointer node = mitk::DataNode::New(); + node->SetData(mitkImage); + node->SetName((name+"_NUM_DIRECTIONS").toStdString().c_str()); + GetDataStorage()->Add(node, m_SelectedFB.back()); + } - if (m_Controls->m_DirectionImagesBox->isChecked()) - { - itk::TractsToVectorImageFilter::ItkDirectionImageType::Pointer itkImg = fOdfFilter->GetDirectionImage(); + if (m_Controls->m_DirectionImagesBox->isChecked()) + { + itk::TractsToVectorImageFilter::ItkDirectionImageType::Pointer itkImg = fOdfFilter->GetDirectionImage(); - if (itkImg.IsNull()) - return; + if (itkImg.IsNull()) + return; - mitk::Image::Pointer mitkImage = mitk::Image::New(); - mitkImage->InitializeByItk( itkImg.GetPointer() ); - mitkImage->SetVolume( itkImg->GetBufferPointer() ); + mitk::Image::Pointer mitkImage = mitk::Image::New(); + mitkImage->InitializeByItk( itkImg.GetPointer() ); + mitkImage->SetVolume( itkImg->GetBufferPointer() ); - mitk::DataNode::Pointer node = mitk::DataNode::New(); - node->SetData(mitkImage); - node->SetName( (name+"_DIRECTIONS").toStdString().c_str()); - GetDataStorage()->Add(node, m_SelectedFB.back()); - } + mitk::DataNode::Pointer node = mitk::DataNode::New(); + node->SetData(mitkImage); + node->SetName( (name+"_DIRECTIONS").toStdString().c_str()); + GetDataStorage()->Add(node, m_SelectedFB.back()); + } } void QmitkFiberQuantificationView::UpdateGui() { - m_Controls->m_ProcessFiberBundleButton->setEnabled(!m_SelectedFB.empty()); - m_Controls->m_ExtractFiberPeaks->setEnabled(!m_SelectedFB.empty()); + m_Controls->m_ProcessFiberBundleButton->setEnabled(!m_SelectedFB.empty()); + m_Controls->m_ExtractFiberPeaks->setEnabled(!m_SelectedFB.empty()); } void QmitkFiberQuantificationView::OnSelectionChanged(berry::IWorkbenchPart::Pointer /*part*/, const QList& nodes) { - //reset existing Vectors containing FiberBundles and PlanarFigures from a previous selection - m_SelectedFB.clear(); - m_SelectedSurfaces.clear(); - m_SelectedImage = nullptr; - - for (mitk::DataNode::Pointer node: nodes) + //reset existing Vectors containing FiberBundles and PlanarFigures from a previous selection + m_SelectedFB.clear(); + m_SelectedSurfaces.clear(); + m_SelectedImage = nullptr; + + for (mitk::DataNode::Pointer node: nodes) + { + if ( dynamic_cast(node->GetData()) ) { - if ( dynamic_cast(node->GetData()) ) - { - m_SelectedFB.push_back(node); - } - else if (dynamic_cast(node->GetData())) - m_SelectedImage = dynamic_cast(node->GetData()); - else if (dynamic_cast(node->GetData())) - { - m_SelectedSurfaces.push_back(dynamic_cast(node->GetData())); - } + m_SelectedFB.push_back(node); } - UpdateGui(); - GenerateStats(); + else if (dynamic_cast(node->GetData())) + m_SelectedImage = dynamic_cast(node->GetData()); + else if (dynamic_cast(node->GetData())) + { + m_SelectedSurfaces.push_back(dynamic_cast(node->GetData())); + } + } + UpdateGui(); + GenerateStats(); } void QmitkFiberQuantificationView::GenerateStats() { - if ( m_SelectedFB.empty() ) - return; + if ( m_SelectedFB.empty() ) + return; - QString stats(""); + QString stats(""); - for( int i=0; i(node->GetData())) { - mitk::DataNode::Pointer node = m_SelectedFB[i]; - if (node.IsNotNull() && dynamic_cast(node->GetData())) - { - if (i>0) - stats += "\n-----------------------------\n"; - stats += QString(node->GetName().c_str()) + "\n"; - mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); - stats += "Number of fibers: "+ QString::number(fib->GetNumFibers()) + "\n"; - stats += "Number of points: "+ QString::number(fib->GetNumberOfPoints()) + "\n"; - stats += "Min. length: "+ QString::number(fib->GetMinFiberLength(),'f',1) + " mm\n"; - stats += "Max. length: "+ QString::number(fib->GetMaxFiberLength(),'f',1) + " mm\n"; - stats += "Mean length: "+ QString::number(fib->GetMeanFiberLength(),'f',1) + " mm\n"; - stats += "Median length: "+ QString::number(fib->GetMedianFiberLength(),'f',1) + " mm\n"; - stats += "Standard deviation: "+ QString::number(fib->GetLengthStDev(),'f',1) + " mm\n"; - - vtkSmartPointer weights = fib->GetFiberWeights(); - if (weights!=nullptr) - { - stats += "Detected fiber weights\n"; -// stats += "Detected fiber weights:\n"; -// float weight=-1; -// int c = 0; -// for (int i=0; iGetSize(); i++) -// if (!mitk::Equal(weights->GetValue(i),weight,0.00001)) -// { -// c++; -// if (weight>0) -// stats += QString::number(i) + ": "+ QString::number(weights->GetValue(i)) + "\n"; -// stats += QString::number(c) + ". Fibers " + QString::number(i+1) + "-"; -// weight = weights->GetValue(i); -// } -// stats += QString::number(weights->GetSize()) + ": "+ QString::number(weight) + "\n"; - } - else - stats += "No fiber weight array found.\n"; - } + if (i>0) + stats += "\n-----------------------------\n"; + stats += QString(node->GetName().c_str()) + "\n"; + mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); + stats += "Number of fibers: "+ QString::number(fib->GetNumFibers()) + "\n"; + stats += "Number of points: "+ QString::number(fib->GetNumberOfPoints()) + "\n"; + stats += "Min. length: "+ QString::number(fib->GetMinFiberLength(),'f',1) + " mm\n"; + stats += "Max. length: "+ QString::number(fib->GetMaxFiberLength(),'f',1) + " mm\n"; + stats += "Mean length: "+ QString::number(fib->GetMeanFiberLength(),'f',1) + " mm\n"; + stats += "Median length: "+ QString::number(fib->GetMedianFiberLength(),'f',1) + " mm\n"; + stats += "Standard deviation: "+ QString::number(fib->GetLengthStDev(),'f',1) + " mm\n"; + + vtkSmartPointer weights = fib->GetFiberWeights(); + if (weights!=nullptr) + { + float weight = -1; + int c = 0; + for (int i=0; iGetSize(); i++) + if (!mitk::Equal(weights->GetValue(i),weight,0.0001)) + { + weight = weights->GetValue(i); + c++; + if (c>1) + break; + } + if (c>1) + stats += "Detected fiber weights. Fibers are not weighted uniformly.\n"; + else + stats += "Fibers are weighted equally.\n"; + } + else + stats += "No fiber weight array found.\n"; } - this->m_Controls->m_StatsTextEdit->setText(stats); + } + this->m_Controls->m_StatsTextEdit->setText(stats); } void QmitkFiberQuantificationView::ProcessSelectedBundles() { - if ( m_SelectedFB.empty() ){ - QMessageBox::information( nullptr, "Warning", "No fibe bundle selected!"); - MITK_WARN("QmitkFiberQuantificationView") << "no fibe bundle selected"; - return; - } - - int generationMethod = m_Controls->m_GenerationBox->currentIndex(); - - for( int i=0; im_GenerationBox->currentIndex(); + + for( int i=0; i(node->GetData())) { - mitk::DataNode::Pointer node = m_SelectedFB[i]; - if (node.IsNotNull() && dynamic_cast(node->GetData())) - { - mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); - QString name(node->GetName().c_str()); - DataNode::Pointer newNode = nullptr; - switch(generationMethod){ - case 0: - newNode = GenerateTractDensityImage(fib, false, true); - name += "_TDI"; - break; - case 1: - newNode = GenerateTractDensityImage(fib, false, false); - name += "_TDI"; - break; - case 2: - newNode = GenerateTractDensityImage(fib, true, false); - name += "_envelope"; - break; - case 3: - newNode = GenerateColorHeatmap(fib); - break; - case 4: - newNode = GenerateFiberEndingsImage(fib); - name += "_fiber_endings"; - break; - case 5: - newNode = GenerateFiberEndingsPointSet(fib); - name += "_fiber_endings"; - break; - } - if (newNode.IsNotNull()) - { - newNode->SetName(name.toStdString()); - GetDataStorage()->Add(newNode); - } - } + mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); + QString name(node->GetName().c_str()); + DataNode::Pointer newNode = nullptr; + switch(generationMethod){ + case 0: + newNode = GenerateTractDensityImage(fib, false, true); + name += "_TDI"; + break; + case 1: + newNode = GenerateTractDensityImage(fib, false, false); + name += "_TDI"; + break; + case 2: + newNode = GenerateTractDensityImage(fib, true, false); + name += "_envelope"; + break; + case 3: + newNode = GenerateColorHeatmap(fib); + break; + case 4: + newNode = GenerateFiberEndingsImage(fib); + name += "_fiber_endings"; + break; + case 5: + newNode = GenerateFiberEndingsPointSet(fib); + name += "_fiber_endings"; + break; + } + if (newNode.IsNotNull()) + { + newNode->SetName(name.toStdString()); + GetDataStorage()->Add(newNode); + } } + } } // generate pointset displaying the fiber endings mitk::DataNode::Pointer QmitkFiberQuantificationView::GenerateFiberEndingsPointSet(mitk::FiberBundle::Pointer fib) { - mitk::PointSet::Pointer pointSet = mitk::PointSet::New(); - vtkSmartPointer fiberPolyData = fib->GetFiberPolyData(); - - int count = 0; - int numFibers = fib->GetNumFibers(); - for( int i=0; i fiberPolyData = fib->GetFiberPolyData(); + + int count = 0; + int numFibers = fib->GetNumFibers(); + for( int i=0; iGetCell(i); + int numPoints = cell->GetNumberOfPoints(); + vtkPoints* points = cell->GetPoints(); + + if (numPoints>0) + { + double* point = points->GetPoint(0); + itk::Point itkPoint; + itkPoint[0] = point[0]; + itkPoint[1] = point[1]; + itkPoint[2] = point[2]; + pointSet->InsertPoint(count, itkPoint); + count++; + } + if (numPoints>2) { - vtkCell* cell = fiberPolyData->GetCell(i); - int numPoints = cell->GetNumberOfPoints(); - vtkPoints* points = cell->GetPoints(); - - if (numPoints>0) - { - double* point = points->GetPoint(0); - itk::Point itkPoint; - itkPoint[0] = point[0]; - itkPoint[1] = point[1]; - itkPoint[2] = point[2]; - pointSet->InsertPoint(count, itkPoint); - count++; - } - if (numPoints>2) - { - double* point = points->GetPoint(numPoints-1); - itk::Point itkPoint; - itkPoint[0] = point[0]; - itkPoint[1] = point[1]; - itkPoint[2] = point[2]; - pointSet->InsertPoint(count, itkPoint); - count++; - } + double* point = points->GetPoint(numPoints-1); + itk::Point itkPoint; + itkPoint[0] = point[0]; + itkPoint[1] = point[1]; + itkPoint[2] = point[2]; + pointSet->InsertPoint(count, itkPoint); + count++; } + } - mitk::DataNode::Pointer node = mitk::DataNode::New(); - node->SetData( pointSet ); - return node; + mitk::DataNode::Pointer node = mitk::DataNode::New(); + node->SetData( pointSet ); + return node; } // generate image displaying the fiber endings mitk::DataNode::Pointer QmitkFiberQuantificationView::GenerateFiberEndingsImage(mitk::FiberBundle::Pointer fib) { - typedef unsigned int OutPixType; - typedef itk::Image OutImageType; + typedef unsigned int OutPixType; + typedef itk::Image OutImageType; + + typedef itk::TractsToFiberEndingsImageFilter< OutImageType > ImageGeneratorType; + ImageGeneratorType::Pointer generator = ImageGeneratorType::New(); + generator->SetFiberBundle(fib); + generator->SetUpsamplingFactor(m_Controls->m_UpsamplingSpinBox->value()); + if (m_SelectedImage.IsNotNull()) + { + OutImageType::Pointer itkImage = OutImageType::New(); + CastToItkImage(m_SelectedImage, itkImage); + generator->SetInputImage(itkImage); + generator->SetUseImageGeometry(true); + } + generator->Update(); + + // get output image + OutImageType::Pointer outImg = generator->GetOutput(); + mitk::Image::Pointer img = mitk::Image::New(); + img->InitializeByItk(outImg.GetPointer()); + img->SetVolume(outImg->GetBufferPointer()); + + // init data node + mitk::DataNode::Pointer node = mitk::DataNode::New(); + node->SetData(img); + return node; +} - typedef itk::TractsToFiberEndingsImageFilter< OutImageType > ImageGeneratorType; - ImageGeneratorType::Pointer generator = ImageGeneratorType::New(); +// generate rgba heatmap from fiber bundle +mitk::DataNode::Pointer QmitkFiberQuantificationView::GenerateColorHeatmap(mitk::FiberBundle::Pointer fib) +{ + typedef itk::RGBAPixel OutPixType; + typedef itk::Image OutImageType; + typedef itk::TractsToRgbaImageFilter< OutImageType > ImageGeneratorType; + ImageGeneratorType::Pointer generator = ImageGeneratorType::New(); + generator->SetFiberBundle(fib); + generator->SetUpsamplingFactor(m_Controls->m_UpsamplingSpinBox->value()); + if (m_SelectedImage.IsNotNull()) + { + itk::Image::Pointer itkImage = itk::Image::New(); + CastToItkImage(m_SelectedImage, itkImage); + generator->SetInputImage(itkImage); + generator->SetUseImageGeometry(true); + } + generator->Update(); + + // get output image + typedef itk::Image OutType; + OutType::Pointer outImg = generator->GetOutput(); + mitk::Image::Pointer img = mitk::Image::New(); + img->InitializeByItk(outImg.GetPointer()); + img->SetVolume(outImg->GetBufferPointer()); + + // init data node + mitk::DataNode::Pointer node = mitk::DataNode::New(); + node->SetData(img); + return node; +} + +// generate tract density image from fiber bundle +mitk::DataNode::Pointer QmitkFiberQuantificationView::GenerateTractDensityImage(mitk::FiberBundle::Pointer fib, bool binary, bool absolute) +{ + mitk::DataNode::Pointer node = mitk::DataNode::New(); + if (binary) + { + typedef unsigned char OutPixType; + typedef itk::Image OutImageType; + + itk::TractDensityImageFilter< OutImageType >::Pointer generator = itk::TractDensityImageFilter< OutImageType >::New(); generator->SetFiberBundle(fib); + generator->SetBinaryOutput(binary); + generator->SetOutputAbsoluteValues(absolute); generator->SetUpsamplingFactor(m_Controls->m_UpsamplingSpinBox->value()); if (m_SelectedImage.IsNotNull()) { - OutImageType::Pointer itkImage = OutImageType::New(); - CastToItkImage(m_SelectedImage, itkImage); - generator->SetInputImage(itkImage); - generator->SetUseImageGeometry(true); + OutImageType::Pointer itkImage = OutImageType::New(); + CastToItkImage(m_SelectedImage, itkImage); + generator->SetInputImage(itkImage); + generator->SetUseImageGeometry(true); + } generator->Update(); // get output image - OutImageType::Pointer outImg = generator->GetOutput(); + typedef itk::Image OutType; + OutType::Pointer outImg = generator->GetOutput(); mitk::Image::Pointer img = mitk::Image::New(); img->InitializeByItk(outImg.GetPointer()); img->SetVolume(outImg->GetBufferPointer()); // init data node - mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData(img); - return node; -} - -// generate rgba heatmap from fiber bundle -mitk::DataNode::Pointer QmitkFiberQuantificationView::GenerateColorHeatmap(mitk::FiberBundle::Pointer fib) -{ - typedef itk::RGBAPixel OutPixType; + } + else + { + typedef float OutPixType; typedef itk::Image OutImageType; - typedef itk::TractsToRgbaImageFilter< OutImageType > ImageGeneratorType; - ImageGeneratorType::Pointer generator = ImageGeneratorType::New(); + + itk::TractDensityImageFilter< OutImageType >::Pointer generator = itk::TractDensityImageFilter< OutImageType >::New(); generator->SetFiberBundle(fib); + generator->SetBinaryOutput(binary); + generator->SetOutputAbsoluteValues(absolute); generator->SetUpsamplingFactor(m_Controls->m_UpsamplingSpinBox->value()); if (m_SelectedImage.IsNotNull()) { - itk::Image::Pointer itkImage = itk::Image::New(); - CastToItkImage(m_SelectedImage, itkImage); - generator->SetInputImage(itkImage); - generator->SetUseImageGeometry(true); + OutImageType::Pointer itkImage = OutImageType::New(); + CastToItkImage(m_SelectedImage, itkImage); + generator->SetInputImage(itkImage); + generator->SetUseImageGeometry(true); + } + //generator->SetDoFiberResampling(false); generator->Update(); // get output image typedef itk::Image OutType; OutType::Pointer outImg = generator->GetOutput(); mitk::Image::Pointer img = mitk::Image::New(); img->InitializeByItk(outImg.GetPointer()); img->SetVolume(outImg->GetBufferPointer()); // init data node - mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData(img); - return node; -} - -// generate tract density image from fiber bundle -mitk::DataNode::Pointer QmitkFiberQuantificationView::GenerateTractDensityImage(mitk::FiberBundle::Pointer fib, bool binary, bool absolute) -{ - mitk::DataNode::Pointer node = mitk::DataNode::New(); - if (binary) - { - typedef unsigned char OutPixType; - typedef itk::Image OutImageType; - - itk::TractDensityImageFilter< OutImageType >::Pointer generator = itk::TractDensityImageFilter< OutImageType >::New(); - generator->SetFiberBundle(fib); - generator->SetBinaryOutput(binary); - generator->SetOutputAbsoluteValues(absolute); - generator->SetUpsamplingFactor(m_Controls->m_UpsamplingSpinBox->value()); - if (m_SelectedImage.IsNotNull()) - { - OutImageType::Pointer itkImage = OutImageType::New(); - CastToItkImage(m_SelectedImage, itkImage); - generator->SetInputImage(itkImage); - generator->SetUseImageGeometry(true); - - } - generator->Update(); - - // get output image - typedef itk::Image OutType; - OutType::Pointer outImg = generator->GetOutput(); - mitk::Image::Pointer img = mitk::Image::New(); - img->InitializeByItk(outImg.GetPointer()); - img->SetVolume(outImg->GetBufferPointer()); - - // init data node - node->SetData(img); - } - else - { - typedef float OutPixType; - typedef itk::Image OutImageType; - - itk::TractDensityImageFilter< OutImageType >::Pointer generator = itk::TractDensityImageFilter< OutImageType >::New(); - generator->SetFiberBundle(fib); - generator->SetBinaryOutput(binary); - generator->SetOutputAbsoluteValues(absolute); - generator->SetUpsamplingFactor(m_Controls->m_UpsamplingSpinBox->value()); - if (m_SelectedImage.IsNotNull()) - { - OutImageType::Pointer itkImage = OutImageType::New(); - CastToItkImage(m_SelectedImage, itkImage); - generator->SetInputImage(itkImage); - generator->SetUseImageGeometry(true); - - } - //generator->SetDoFiberResampling(false); - generator->Update(); - - // get output image - typedef itk::Image OutType; - OutType::Pointer outImg = generator->GetOutput(); - mitk::Image::Pointer img = mitk::Image::New(); - img->InitializeByItk(outImg.GetPointer()); - img->SetVolume(outImg->GetBufferPointer()); - - // init data node - node->SetData(img); - } - return node; + } + return node; }