diff --git a/CMakeLists.txt.user b/CMakeLists.txt.user new file mode 100644 index 0000000000..5e8fb941b0 --- /dev/null +++ b/CMakeLists.txt.user @@ -0,0 +1,173 @@ + + + + ProjectExplorer.Project.ActiveTarget + 0 + + + ProjectExplorer.Project.EditorSettings + + Default + + + + ProjectExplorer.Project.Target.0 + + Desktop + + CMakeProjectManager.DefaultCMakeTarget + 0 + 0 + 4 + + /local/neher/mitkBinaryCTK/MITK-superbuild/MITK-build + + + + + -j15 + + + false + Make + + CMakeProjectManager.MakeStep + + 1 + Build + + ProjectExplorer.BuildSteps.Build + + + + + clean + + + true + Make + + CMakeProjectManager.MakeStep + + 1 + Clean + + ProjectExplorer.BuildSteps.Clean + + 2 + false + + all + + CMakeProjectManager.CMakeBuildConfiguration + + 1 + + + 0 + Deploy + + ProjectExplorer.BuildSteps.Deploy + + 1 + No deployment + + ProjectExplorer.DefaultDeployConfiguration + + 1 + + 2 + solstice + + /local/neher/mitkBinaryCTK/MITK-superbuild/MITK-build/bin/solstice + false + + + /local/neher/mitkBinaryCTK/MITK-superbuild/MITK-build/bin + solstice + + CMakeProjectManager.CMakeRunConfiguration. + 3768 + true + false + + + 2 + solstice_ui + + /local/neher/mitkBinaryCTK/MITK-superbuild/MITK-build/bin/solstice_ui + false + + + /local/neher/mitkBinaryCTK/MITK-superbuild/MITK-build/bin + solstice_ui + + CMakeProjectManager.CMakeRunConfiguration. + 3768 + true + false + + + 2 + IGTTutorialStep1 + + /local/neher/mitkBinaryCTK/MITK-superbuild/MITK-build/bin/IGTTutorialStep1 + false + + + /local/neher/mitkBinaryCTK/MITK-superbuild/MITK-build/bin + IGTTutorialStep1 + + CMakeProjectManager.CMakeRunConfiguration. + 3768 + true + false + + + 2 + CoreApp + + /local/neher/mitkBinaryCTK/MITK-superbuild/MITK-build/bin/CoreApp + false + + + /local/neher/mitkBinaryCTK/MITK-superbuild/MITK-build/bin + CoreApp + + CMakeProjectManager.CMakeRunConfiguration. + 3768 + true + false + + + 2 + ExtApp + + /local/neher/mitkBinaryCTK/MITK-superbuild/MITK-build/bin/ExtApp + false + + + /local/neher/mitkBinaryCTK/MITK-superbuild/MITK-build/bin + ExtApp + + CMakeProjectManager.CMakeRunConfiguration. + 3768 + true + false + + 5 + + + + ProjectExplorer.Project.TargetCount + 1 + + + ProjectExplorer.Project.Updater.EnvironmentId + {9a4fa05b-ac40-4923-8a62-06d1bd92182a} + + + ProjectExplorer.Project.Updater.FileVersion + 8 + + diff --git a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/files.cmake b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/files.cmake index 11ea4a277a..30c992ecd8 100644 --- a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/files.cmake +++ b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/files.cmake @@ -1,76 +1,80 @@ SET(SRC_CPP_FILES QmitkODFDetailsWidget.cpp QmitkODFRenderWidget.cpp ) SET(INTERNAL_CPP_FILES mitkPluginActivator.cpp QmitkQBallReconstructionView.cpp QmitkPreprocessingView.cpp QmitkDiffusionDicomImportView.cpp QmitkDiffusionQuantificationView.cpp QmitkTensorReconstructionView.cpp QmitkDiffusionImagingPublicPerspective.cpp QmitkControlVisualizationPropertiesView.cpp QmitkODFDetailsView.cpp QmitkGlobalFiberTrackingView.cpp + QmitkFiberBundleOperationsView.cpp ) SET(UI_FILES src/internal/QmitkQBallReconstructionViewControls.ui src/internal/QmitkPreprocessingViewControls.ui src/internal/QmitkDiffusionDicomImportViewControls.ui src/internal/QmitkDiffusionQuantificationViewControls.ui src/internal/QmitkTensorReconstructionViewControls.ui src/internal/QmitkControlVisualizationPropertiesViewControls.ui src/internal/QmitkODFDetailsViewControls.ui src/internal/QmitkGlobalFiberTrackingViewControls.ui + src/internal/QmitkFiberBundleOperationsViewControls.ui ) SET(MOC_H_FILES src/internal/mitkPluginActivator.h src/internal/QmitkQBallReconstructionView.h src/internal/QmitkPreprocessingView.h src/internal/QmitkDiffusionDicomImportView.h src/internal/QmitkDiffusionImagingPublicPerspective.h src/internal/QmitkDiffusionQuantificationView.h src/internal/QmitkTensorReconstructionView.h src/internal/QmitkControlVisualizationPropertiesView.h src/internal/QmitkODFDetailsView.h src/QmitkODFRenderWidget.h src/QmitkODFDetailsWidget.h src/internal/QmitkGlobalFiberTrackingView.h + src/internal/QmitkFiberBundleOperationsView.h ) SET(CACHED_RESOURCE_FILES # list of resource files which can be used by the plug-in # system without loading the plug-ins shared library, # for example the icon used in the menu and tabs for the # plug-in views in the workbench plugin.xml resources/preprocessing.png resources/dwiimport.png resources/quantification.png resources/reconodf.png resources/recontensor.png resources/vizControls.png resources/OdfDetails.png resources/GlobalTracking.png + resources/FiberBundleOperations.png ) SET(QRC_FILES # uncomment the following line if you want to use Qt resources resources/QmitkDiffusionImaging.qrc ) SET(CPP_FILES ) foreach(file ${SRC_CPP_FILES}) SET(CPP_FILES ${CPP_FILES} src/${file}) endforeach(file ${SRC_CPP_FILES}) foreach(file ${INTERNAL_CPP_FILES}) SET(CPP_FILES ${CPP_FILES} src/internal/${file}) endforeach(file ${INTERNAL_CPP_FILES}) diff --git a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/plugin.xml b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/plugin.xml index fe30ab1d8a..4d54407c7c 100644 --- a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/plugin.xml +++ b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/plugin.xml @@ -1,68 +1,75 @@ + + + + diff --git a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/resources/FiberBundleOperations.png b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/resources/FiberBundleOperations.png new file mode 100644 index 0000000000..8187b57ac0 Binary files /dev/null and b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/resources/FiberBundleOperations.png differ diff --git a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/resources/QmitkDiffusionImaging.qrc b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/resources/QmitkDiffusionImaging.qrc index 791d03ae51..355986b909 100644 --- a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/resources/QmitkDiffusionImaging.qrc +++ b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/resources/QmitkDiffusionImaging.qrc @@ -1,23 +1,26 @@ qball.png tensor.png dwi.png dwiimport.png quantification.png reconodf.png recontensor.png texIntONIcon.png texIntOFFIcon.png vizControls.png Refresh_48.png QBallData24.png glyphsoff_C.png glyphsoff_S.png glyphsoff_T.png glyphson_C.png glyphson_S.png glyphson_T.png FiberBundle.png + rectangle.png + circle.png + polygon.png diff --git a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/resources/circle.png b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/resources/circle.png new file mode 100644 index 0000000000..9d88f3ab5d Binary files /dev/null and b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/resources/circle.png differ diff --git a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/resources/polygon.png b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/resources/polygon.png new file mode 100644 index 0000000000..f3c651d6fd Binary files /dev/null and b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/resources/polygon.png differ diff --git a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/resources/rectangle.png b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/resources/rectangle.png new file mode 100644 index 0000000000..33959262c0 Binary files /dev/null and b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/resources/rectangle.png differ diff --git a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberBundleOperationsView.cpp b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberBundleOperationsView.cpp new file mode 100644 index 0000000000..83e24bf4e7 --- /dev/null +++ b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberBundleOperationsView.cpp @@ -0,0 +1,1829 @@ +/*========================================================================= + + Program: Medical Imaging & Interaction Toolkit + Language: C++ + Date: $Date: 2010-03-31 16:40:27 +0200 (Mi, 31 Mrz 2010) $ + Version: $Revision: 21975 $ + + Copyright (c) German Cancer Research Center, Division of Medical and + Biological Informatics. All rights reserved. + See MITKCopyright.txt or http://www.mitk.org/copyright.html for details. + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notices for more information. + + =========================================================================*/ + + +// Blueberry +#include +#include + + + +// Qmitk +#include "QmitkFiberBundleOperationsView.h" +#include "QmitkStdMultiWidget.h" + +// Qt +#include + +//MITK +#include "mitkNodePredicateProperty.h" +//#include "mitkNodePredicateAND.h" +#include "mitkImageCast.h" + +#include "mitkPointSet.h" + + +#include "mitkPlanarCircle.h" +#include "mitkPlanarPolygon.h" +#include +#include "mitkPlanarFigureInteractor.h" +#include "mitkGlobalInteraction.h" +#include +#include +#include +#include +#include +#include +#include +#include + + +const std::string QmitkFiberBundleOperationsView::VIEW_ID = "org.mitk.views.fiberBundleOperations"; +const std::string id_DataManager = "org.mitk.views.datamanager"; +using namespace berry; +using namespace mitk; + + + +QmitkFiberBundleOperationsView::QmitkFiberBundleOperationsView() +: QmitkFunctionality() +, m_Controls( 0 ) +, m_MultiWidget( NULL ) +, m_EllipseCounter(0) +, m_PolygonCounter(0) +//, m_SelectedFBNodes( NULL ) +//, m_SelectedPFNodes(0) +, m_UpsamplingFactor(5) +{ + +} + +// Destructor +QmitkFiberBundleOperationsView::~QmitkFiberBundleOperationsView() +{ + +} + + +void QmitkFiberBundleOperationsView::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::QmitkFiberBundleOperationsViewControls; + m_Controls->setupUi( parent ); + m_Controls->doExtractFibersButton->setDisabled(true); + m_Controls->PFCompoANDButton->setDisabled(true); + m_Controls->PFCompoORButton->setDisabled(true); + m_Controls->PFCompoNOTButton->setDisabled(true); + m_Controls->PFCompoDELButton->setDisabled(true); + + connect( m_Controls->doExtractFibersButton, SIGNAL(clicked()), this, SLOT(DoFiberExtraction()) ); + //connect( m_Controls->comboBox_fiberAlgo, SIGNAL(selected()), this, SLOT(handleAlgoSelection() ); + connect( m_Controls->m_CircleButton, SIGNAL( clicked() ), this, SLOT( ActionDrawEllipseTriggered() ) ); + connect( m_Controls->m_PolygonButton, SIGNAL( clicked() ), this, SLOT( ActionDrawPolygonTriggered() ) ); + connect(m_Controls->PFCompoANDButton, SIGNAL(clicked()), this, SLOT(generatePFCompo_AND()) ); + connect(m_Controls->PFCompoORButton, SIGNAL(clicked()), this, SLOT(generatePFCompo_OR()) ); + connect(m_Controls->PFCompoNOTButton, SIGNAL(clicked()), this, SLOT(generatePFCompo_NOT()) ); + connect(m_Controls->PFCompoDELButton, SIGNAL(clicked()), this, SLOT(deletePFCompo()) ); + + connect(m_Controls->m_JoinBundles, SIGNAL(clicked()), this, SLOT(JoinBundles()) ); + connect(m_Controls->m_SubstractBundles, SIGNAL(clicked()), this, SLOT(SubstractBundles()) ); + connect(m_Controls->m_GenerateROIImage, SIGNAL(clicked()), this, SLOT(GenerateROIImage()) ); + + connect( m_Controls->m_GenerationStartButton, SIGNAL(clicked()), this, SLOT(GenerationStart()) ); + } +} + +void QmitkFiberBundleOperationsView::GenerateROIImage(){ + + if (m_Image.IsNull() || m_SelectedPF.empty()) + return; + + mitk::Image* image = const_cast(m_Image.GetPointer()); + + MaskImage3DType::Pointer temp = MaskImage3DType::New(); + mitk::CastToItkImage(m_Image, temp); + + m_PlanarFigureImage = MaskImage3DType::New(); + m_PlanarFigureImage->SetSpacing( temp->GetSpacing() ); // Set the image spacing + m_PlanarFigureImage->SetOrigin( temp->GetOrigin() ); // Set the image origin + m_PlanarFigureImage->SetDirection( temp->GetDirection() ); // Set the image direction + m_PlanarFigureImage->SetRegions( temp->GetLargestPossibleRegion() ); + m_PlanarFigureImage->Allocate(); + m_PlanarFigureImage->FillBuffer( 0 ); + + + for (int i=0; i(m_SelectedPF.at(i)->GetData()), image); + + DataNode::Pointer node = DataNode::New(); + Image::Pointer tmpImage = Image::New(); + tmpImage->InitializeByItk(m_PlanarFigureImage.GetPointer()); + tmpImage->SetVolume(m_PlanarFigureImage->GetBufferPointer()); + node->SetData(tmpImage); + node->SetName("planarFigureImage"); + this->GetDefaultDataStorage()->Add(node); +} + +void QmitkFiberBundleOperationsView::CompositeExtraction(mitk::PlanarFigure::Pointer planarFigure, mitk::Image* image) +{ + if (dynamic_cast(planarFigure.GetPointer())){ + mitk::PlanarFigureComposite::Pointer pfcomp = dynamic_cast(planarFigure.GetPointer()); + + for (int i=0; igetNumberOfChildren(); ++i) + { + CompositeExtraction(pfcomp->getChildAt(i), image); + } + } + else + { + m_PlanarFigure = planarFigure; + AccessFixedDimensionByItk_3( + image, + InternalReorientImagePlane, 3, + image->GetGeometry(), + m_PlanarFigure->GetGeometry(), -1); + + AccessFixedDimensionByItk_1( + m_InternalImage, + InternalCalculateMaskFromPlanarFigure, + 3, 2 ); + } +} + +template < typename TPixel, unsigned int VImageDimension > + void QmitkFiberBundleOperationsView::InternalReorientImagePlane( + const itk::Image< TPixel, VImageDimension > *image, + mitk::Geometry3D* imggeo, mitk::Geometry3D* planegeo3D, int additionalIndex ) +{ + + MITK_INFO << "InternalReorientImagePlane() start"; + + 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->GetParametricExtentInMM(0) / spacing[0]; + size[1] = planegeo->GetParametricExtentInMM(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::BSplineInterpolateImageFunction + // InterpolatorType; + typedef typename itk::LinearInterpolateImageFunction InterpolatorType; + + typename InterpolatorType::Pointer interpolator + = InterpolatorType::New(); + + interpolator->SetInputImage( image ); + + resampler->SetInterpolator( interpolator ); + + } + + // Other resampling options + resampler->SetInput( image ); + resampler->SetDefaultPixelValue(0); + + MITK_INFO << "Resampling requested image plane ... "; + resampler->Update(); + MITK_INFO << " ... done"; + + 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 QmitkFiberBundleOperationsView::InternalCalculateMaskFromPlanarFigure( + itk::Image< TPixel, VImageDimension > *image, unsigned int axis ) +{ + + MITK_INFO << "InternalCalculateMaskFromPlanarFigure() start"; + + typedef itk::Image< TPixel, VImageDimension > ImageType; + typedef itk::CastImageFilter< ImageType, MaskImage3DType > CastFilterType; + + // Generate mask image as new image with same header as input image and + // initialize with "1". + MaskImage3DType::Pointer newMaskImage = MaskImage3DType::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 Geometry2D *planarFigureGeometry2D = m_PlanarFigure->GetGeometry2D(); + const PlanarFigure::PolyLineType planarFigurePolyline = m_PlanarFigure->GetPolyLine( 0 ); + const Geometry3D *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; + std::vector indices; + + 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->Point; + planarFigureGeometry2D->WorldToIndex(point2D, point2D); + point2D[0] -= 0.5/m_UpsamplingFactor; + point2D[1] -= 0.5/m_UpsamplingFactor; + planarFigureGeometry2D->IndexToWorld(point2D, point2D); + planarFigureGeometry2D->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.5; + else if (point3D[i0]>bounds[0]) + point3D[i0] = bounds[0]-0.5; + + if (point3D[i1]<0) + point3D[i1] = 0.5; + else if (point3D[i1]>bounds[1]) + point3D[i1] = bounds[1]-0.5; + + points->InsertNextPoint( point3D[i0], point3D[i1], -0.5 ); + numberOfPoints++; + } + else + { + imageGeometry3D->WorldToIndex( point3D, point3D ); + point3D[i0] += 0.5; + point3D[i1] += 0.5; + + // 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->SetInput( 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->SetInput( extrudeFilter->GetOutput() ); + + + + // Export from ITK to VTK (to use a VTK filter) + typedef itk::VTKImageImport< MaskImage3DType > ImageImportType; + typedef itk::VTKImageExport< MaskImage3DType > 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->SetStencil( polyDataToImageStencil->GetOutput() ); + 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 = itmask.Begin(); + itimage = itimage.Begin(); + + typename ImageType::SizeType lowersize = {{9999999999,9999999999,9999999999}}; + 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 image +// typedef itk::RegionOfInterestImageFilter< ImageType, ImageType > ROIFilterType; +// typename ROIFilterType::Pointer roi = ROIFilterType::New(); +// roi->SetRegionOfInterest(cropRegion); +// roi->SetInput(image); +// roi->Update(); + +// m_InternalImage = mitk::Image::New(); +// m_InternalImage->InitializeByItk(roi->GetOutput()); +// m_InternalImage->SetVolume(roi->GetOutput()->GetBufferPointer()); + + // crop internal mask + typedef itk::RegionOfInterestImageFilter< MaskImage3DType, MaskImage3DType > ROIMaskFilterType; + typename ROIMaskFilterType::Pointer roi2 = ROIMaskFilterType::New(); + roi2->SetRegionOfInterest(cropRegion); + roi2->SetInput(m_InternalImageMask3D); + roi2->Update(); + m_InternalImageMask3D = roi2->GetOutput(); + + DataNode::Pointer node = DataNode::New(); + Image::Pointer tmpImage = Image::New(); + tmpImage->InitializeByItk(m_InternalImageMask3D.GetPointer()); + tmpImage->SetVolume(m_InternalImageMask3D->GetBufferPointer()); + node->SetData(tmpImage); + this->GetDefaultDataStorage()->Add(node); + + Image::Pointer tmpImage2 = Image::New(); + tmpImage2->InitializeByItk(m_PlanarFigureImage.GetPointer()); + const Geometry3D *pfImageGeometry3D = tmpImage2->GetGeometry( 0 ); + + const Geometry3D *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); + + index[0] = point[0]; + index[1] = point[1]; + index[2] = point[2]; + + 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 QmitkFiberBundleOperationsView::StdMultiWidgetAvailable (QmitkStdMultiWidget &stdMultiWidget) +{ + m_MultiWidget = &stdMultiWidget; +} + + +void QmitkFiberBundleOperationsView::StdMultiWidgetNotAvailable() +{ + m_MultiWidget = NULL; +} + +/* OnSelectionChanged is registered to SelectionService, therefore no need to + implement SelectionService Listener explicitly */ + +void QmitkFiberBundleOperationsView::OnSelectionChanged( std::vector nodes ) +{ + + if ( !this->IsVisible() ) + { + // do nothing if nobody wants to see me :-( + return; + } + + + //reset existing Vectors containing FiberBundles and PlanarFigures from a previous selection + m_SelectedFB.clear(); + m_SelectedPF.clear(); + m_Image = NULL; + + + //differ between 2 scenarios... + // 1) add PF to an existing Spaghetti + // 2) BOOLEAN OPERATORS HANDLING ... selection of multiple PF, and or multiple Spaghetti + + + //scenario 1 + + for( std::vector::iterator it = nodes.begin(); + it != nodes.end(); ++it ) + { + + mitk::DataNode::Pointer node = *it; + if ( dynamic_cast(node->GetData()) ) + { + // MITK_INFO << "onselectionchg(): " << node->GetData()->GetNameOfClass(); + m_SelectedFB.push_back(node); + mitk::FiberBundle::Pointer bundle = dynamic_cast(node->GetData()); + QString numFibers; + this->m_Controls->m_NumFibersLabel->setText(numFibers.setNum(bundle->GetNumTracts())); + + } else if (dynamic_cast(node->GetData())){ + + // MITK_INFO << "onselectionchg(): " << node->GetData()->GetNameOfClass(); + m_SelectedPF.push_back(node); + + } else if (dynamic_cast(node->GetData())){ + m_Image = dynamic_cast(node->GetData()); + } + + } + + if (m_SelectedPF.size() >= 1 && m_Image.IsNotNull()) + m_Controls->m_GenerateROIImage->setEnabled(true); + else + m_Controls->m_GenerateROIImage->setEnabled(false); + + //quick and dirty control structure for extraction ATM + if (m_SelectedPF.size() == 1) { + + + //############################################# + //### PLANAR FIGURE COMPOSIT SelectionLOGIC ### + //############################################# + + m_Controls->PFCompoANDButton->setDisabled(true); + m_Controls->PFCompoORButton->setDisabled(true); + m_Controls->PFCompoNOTButton->setEnabled(true); + m_Controls->PFCompoDELButton->setDisabled(true); + + //if planarFigureComposite selected, activate DEL button + if ( dynamic_cast(m_SelectedPF.at(0)->GetData()) ) + { + m_Controls->PFCompoDELButton->setEnabled(true); + } + + + } else if (m_SelectedPF.size() > 1) { + + //############################################# + //### PLANAR FIGURE COMPOSIT SelectionLOGIC ### + //############################################# + // if 2 PlanarFigure objects (PFCircle, PFPoly, ..., PFComposite) selected, + //then activate AND OR NOT PF Composite Buttons + + m_Controls->PFCompoANDButton->setEnabled(true); + m_Controls->PFCompoORButton->setEnabled(true); + m_Controls->PFCompoNOTButton->setDisabled(true); + m_Controls->PFCompoDELButton->setDisabled(true); + } + + + + if (m_SelectedFB.size() == 1 && m_SelectedPF.size() == 1) + { + + //############################################# + //##### EXTRACT FIBERBUNDLE SelectionLOGIC #### + //############################################# + + m_Controls->doExtractFibersButton->setEnabled(true); + + + + } else if (nodes.size() > 1 ) { + //scenario 2 + + + } else if (nodes.empty()) { + + //############################################# + //##### EXTRACT FIBERBUNDLE SelectionLOGIC #### + //############################################# + + m_Controls->doExtractFibersButton->setDisabled(true); + + + //############################################# + //### PLANAR FIGURE COMPOSIT SelectionLOGIC ### + //############################################# + + m_Controls->PFCompoANDButton->setDisabled(true); + m_Controls->PFCompoORButton->setDisabled(true); + m_Controls->PFCompoNOTButton->setDisabled(true); + + } + + if (m_SelectedFB.size() == 2) + { + + //############################################# + //##### JOIN/SUBSTRACT FIBERBUNDLES #### + //############################################# + + m_Controls->m_JoinBundles->setEnabled(true); + m_Controls->m_SubstractBundles->setEnabled(true); + + } + else + { + + m_Controls->m_JoinBundles->setEnabled(false); + m_Controls->m_SubstractBundles->setEnabled(false); + + } + +} + + +void QmitkFiberBundleOperationsView::ActionDrawPolygonTriggered() +{ + // bool checked = m_Controls->m_PolygonButton->isChecked(); + // if(!this->AssertDrawingIsPossible(checked)) + // return; + + mitk::PlanarPolygon::Pointer figure = mitk::PlanarPolygon::New(); + figure->ClosedOn(); + this->AddFigureToDataStorage(figure, QString("Polygon%1").arg(++m_PolygonCounter)); + + MITK_INFO << "PlanarPolygon created ..."; + + mitk::DataStorage::SetOfObjects::ConstPointer _NodeSet = this->GetDefaultDataStorage()->GetAll(); + mitk::DataNode* node = 0; + mitk::PlanarFigureInteractor::Pointer figureInteractor = 0; + mitk::PlanarFigure* figureP = 0; + + for(mitk::DataStorage::SetOfObjects::ConstIterator it=_NodeSet->Begin(); it!=_NodeSet->End() + ; it++) + { + node = const_cast(it->Value().GetPointer()); + figureP = dynamic_cast(node->GetData()); + + if(figureP) + { + figureInteractor = dynamic_cast(node->GetInteractor()); + + if(figureInteractor.IsNull()) + figureInteractor = mitk::PlanarFigureInteractor::New("PlanarFigureInteractor", node); + + mitk::GlobalInteraction::GetInstance()->AddInteractor(figureInteractor); + } + } + +} + + +void QmitkFiberBundleOperationsView::ActionDrawEllipseTriggered() +{ + //bool checked = m_Controls->m_CircleButton->isChecked(); + //if(!this->AssertDrawingIsPossible(checked)) + // return; + + mitk::PlanarCircle::Pointer figure = mitk::PlanarCircle::New(); + this->AddFigureToDataStorage(figure, QString("Circle%1").arg(++m_EllipseCounter)); + + this->GetDataStorage()->Modified(); + MITK_INFO << "PlanarCircle created ..."; + + //call + + mitk::DataStorage::SetOfObjects::ConstPointer _NodeSet = this->GetDefaultDataStorage()->GetAll(); + mitk::DataNode* node = 0; + mitk::PlanarFigureInteractor::Pointer figureInteractor = 0; + mitk::PlanarFigure* figureP = 0; + + for(mitk::DataStorage::SetOfObjects::ConstIterator it=_NodeSet->Begin(); it!=_NodeSet->End() + ; it++) + { + node = const_cast(it->Value().GetPointer()); + figureP = dynamic_cast(node->GetData()); + + if(figureP) + { + figureInteractor = dynamic_cast(node->GetInteractor()); + + if(figureInteractor.IsNull()) + figureInteractor = mitk::PlanarFigureInteractor::New("PlanarFigureInteractor", node); + + mitk::GlobalInteraction::GetInstance()->AddInteractor(figureInteractor); + } + } + + + + +} + +void QmitkFiberBundleOperationsView::Activated() +{ + + MITK_INFO << "FB OPerations ACTIVATED()"; + /* + + mitk::DataStorage::SetOfObjects::ConstPointer _NodeSet = this->GetDefaultDataStorage()->GetAll(); + mitk::DataNode* node = 0; + mitk::PlanarFigureInteractor::Pointer figureInteractor = 0; + mitk::PlanarFigure* figure = 0; + + for(mitk::DataStorage::SetOfObjects::ConstIterator it=_NodeSet->Begin(); it!=_NodeSet->End() + ; it++) + { + node = const_cast(it->Value().GetPointer()); + figure = dynamic_cast(node->GetData()); + + if(figure) + { + figureInteractor = dynamic_cast(node->GetInteractor()); + + if(figureInteractor.IsNull()) + figureInteractor = mitk::PlanarFigureInteractor::New("PlanarFigureInteractor", node); + + mitk::GlobalInteraction::GetInstance()->AddInteractor(figureInteractor); + } + } + + */ + + +} + +void QmitkFiberBundleOperationsView::AddFigureToDataStorage(mitk::PlanarFigure* figure, const QString& name, + const char *propertyKey, mitk::BaseProperty *property ) +{ + //set desired data to DataNode where Planarfigure is stored + mitk::DataNode::Pointer newNode = mitk::DataNode::New(); + newNode->SetName(name.toStdString()); + newNode->SetData(figure); + + newNode->AddProperty( "planarfigure.default.line.color", mitk::ColorProperty::New(1.0,0.0,0.0)); + newNode->AddProperty( "planarfigure.line.width", mitk::FloatProperty::New(2.0)); + newNode->AddProperty( "planarfigure.drawshadow", mitk::BoolProperty::New(true)); + + + + + + + newNode->AddProperty( "selected", mitk::BoolProperty::New(true) ); + newNode->AddProperty( "planarfigure.ishovering", mitk::BoolProperty::New(true) ); + newNode->AddProperty( "planarfigure.drawoutline", mitk::BoolProperty::New(true) ); + newNode->AddProperty( "planarfigure.drawquantities", mitk::BoolProperty::New(false) ); + newNode->AddProperty( "planarfigure.drawshadow", mitk::BoolProperty::New(true) ); + + newNode->AddProperty( "planarfigure.line.width", mitk::FloatProperty::New(3.0) ); + newNode->AddProperty( "planarfigure.shadow.widthmodifier", mitk::FloatProperty::New(2.0) ); + newNode->AddProperty( "planarfigure.outline.width", mitk::FloatProperty::New(2.0) ); + newNode->AddProperty( "planarfigure.helperline.width", mitk::FloatProperty::New(2.0) ); + +// PlanarFigureControlPointStyleProperty::Pointer styleProperty = +// dynamic_cast< PlanarFigureControlPointStyleProperty* >( node->GetProperty( "planarfigure.controlpointshape" ) ); +// if ( styleProperty.IsNotNull() ) +// { +// m_ControlPointShape = styleProperty->GetShape(); +// } + + newNode->AddProperty( "planarfigure.default.line.color", mitk::ColorProperty::New(1.0,1.0,1.0) ); + newNode->AddProperty( "planarfigure.default.line.opacity", mitk::FloatProperty::New(2.0) ); + newNode->AddProperty( "planarfigure.default.outline.color", mitk::ColorProperty::New(1.0,0.0,0.0) ); + newNode->AddProperty( "planarfigure.default.outline.opacity", mitk::FloatProperty::New(2.0) ); + newNode->AddProperty( "planarfigure.default.helperline.color", mitk::ColorProperty::New(1.0,0.0,0.0) ); + newNode->AddProperty( "planarfigure.default.helperline.opacity", mitk::FloatProperty::New(2.0) ); + newNode->AddProperty( "planarfigure.default.markerline.color", mitk::ColorProperty::New(0.0,0.0,0.0) ); + newNode->AddProperty( "planarfigure.default.markerline.opacity", mitk::FloatProperty::New(2.0) ); + newNode->AddProperty( "planarfigure.default.marker.color", mitk::ColorProperty::New(1.0,1.0,1.0) ); + newNode->AddProperty( "planarfigure.default.marker.opacity",mitk::FloatProperty::New(2.0) ); + + newNode->AddProperty( "planarfigure.hover.line.color", mitk::ColorProperty::New(1.0,0.0,0.0) ); + newNode->AddProperty( "planarfigure.hover.line.opacity", mitk::FloatProperty::New(2.0) ); + newNode->AddProperty( "planarfigure.hover.outline.color", mitk::ColorProperty::New(1.0,0.0,0.0) ); + newNode->AddProperty( "planarfigure.hover.outline.opacity", mitk::FloatProperty::New(2.0) ); + newNode->AddProperty( "planarfigure.hover.helperline.color", mitk::ColorProperty::New(1.0,0.0,0.0) ); + newNode->AddProperty( "planarfigure.hover.helperline.opacity", mitk::FloatProperty::New(2.0) ); + newNode->AddProperty( "planarfigure.hover.markerline.color", mitk::ColorProperty::New(1.0,0.0,0.0) ); + newNode->AddProperty( "planarfigure.hover.markerline.opacity", mitk::FloatProperty::New(2.0) ); + newNode->AddProperty( "planarfigure.hover.marker.color", mitk::ColorProperty::New(1.0,0.0,0.0) ); + newNode->AddProperty( "planarfigure.hover.marker.opacity", mitk::FloatProperty::New(2.0) ); + + newNode->AddProperty( "planarfigure.selected.line.color", mitk::ColorProperty::New(1.0,0.0,0.0) ); + newNode->AddProperty( "planarfigure.selected.line.opacity",mitk::FloatProperty::New(2.0) ); + newNode->AddProperty( "planarfigure.selected.outline.color", mitk::ColorProperty::New(1.0,0.0,0.0) ); + newNode->AddProperty( "planarfigure.selected.outline.opacity", mitk::FloatProperty::New(2.0)); + newNode->AddProperty( "planarfigure.selected.helperline.color", mitk::ColorProperty::New(1.0,0.0,0.0) ); + newNode->AddProperty( "planarfigure.selected.helperline.opacity",mitk::FloatProperty::New(2.0) ); + newNode->AddProperty( "planarfigure.selected.markerline.color", mitk::ColorProperty::New(1.0,0.0,0.0) ); + newNode->AddProperty( "planarfigure.selected.markerline.opacity", mitk::FloatProperty::New(2.0) ); + newNode->AddProperty( "planarfigure.selected.marker.color", mitk::ColorProperty::New(1.0,0.0,0.0) ); + newNode->AddProperty( "planarfigure.selected.marker.opacity",mitk::FloatProperty::New(2.0)); + + + + + + + + + + + + + + + + // Add custom property, if available + //if ( (propertyKey != NULL) && (property != NULL) ) + //{ + // newNode->AddProperty( propertyKey, property ); + //} + + //get current selected DataNode -which should be a FiberBundle- and set PlanarFigure as child + + //this->GetDataStorage()->GetNodes() + + + + // mitk::FiberBundle::Pointer selectedFBNode = m_SelectedFBNodes.at(0); + + // figure drawn on the topmost layer / image + this->GetDataStorage()->Add(newNode ); + + std::vector selectedNodes = GetDataManagerSelection(); + for(unsigned int i = 0; i < selectedNodes.size(); i++) + { + selectedNodes[i]->SetSelected(false); + } + //selectedNodes = m_SelectedPlanarFigureNodes->GetNodes(); + /*for(unsigned int i = 0; i < selectedNodes.size(); i++) + { + selectedNodes[i]->SetSelected(false); + } + */ + newNode->SetSelected(true); + + //Select(newNode); +} + +void QmitkFiberBundleOperationsView::DoFiberExtraction() +{ + mitk::FiberBundle::Pointer selFB = dynamic_cast(m_SelectedFB.at(0)->GetData()); + mitk::PlanarFigure::Pointer selPF = dynamic_cast (m_SelectedPF.at(0)->GetData()); + + + std::vector extFBset = selFB->extractFibersByPF(selPF); + + + //MITK_INFO << "returned vector in FBOperationsView: " << extFBset.size(); + // for(std::vector::iterator dispIt = extFBset.begin(); dispIt != extFBset.end(); dispIt++) + // { + // MITK_INFO << "vector DTI ID: " << *dispIt; + // + // } + + + mitk::FiberBundle::Pointer extFB = selFB->extractFibersById(extFBset); + MITK_INFO << " Number Of Tracts in sourceFiberBundle: " << selFB->GetNumTracts(); + MITK_INFO << " Number Of Tracts in extractedFiberBundle: " << extFB->GetNumTracts(); + + mitk::DataNode::Pointer fbNode; + fbNode = mitk::DataNode::New(); + fbNode->SetData(extFB); + fbNode->SetName("extGroupFinberBundle"); + fbNode->SetVisibility(true); + GetDataStorage()->Add(fbNode); + + + +} + +void QmitkFiberBundleOperationsView::generatePFCompo_AND() +{ + mitk::PlanarFigureComposite::Pointer PFCAnd = mitk::PlanarFigureComposite::New(); + PFCAnd->setOperationType(mitk::PFCOMPOSITION_AND_OPERATION); + + + for( std::vector::iterator it = m_SelectedPF.begin(); + it != m_SelectedPF.end(); ++it ) + { + + mitk::DataNode::Pointer nodePF = *it; + mitk::PlanarFigure::Pointer tmpPF = dynamic_cast( nodePF->GetData() ); + PFCAnd->addPlanarFigure( tmpPF ); + PFCAnd->addDataNode( nodePF ); + PFCAnd->setDisplayName("AND_COMPO"); + // MITK_INFO << "PFCAND(): added to AND PF" << nodePF->GetName(); + + } + + debugPFComposition(PFCAnd, 0); + + + this->addPFCompositionToDataStorage(PFCAnd, NULL /*parent*/); + +} + + +void QmitkFiberBundleOperationsView::debugPFComposition(mitk::PlanarFigureComposite::Pointer pfc, int itLevelStatus) +{ + int myLevel = itLevelStatus; + if (myLevel == 0) + { + MITK_INFO << "############################################## " ; + MITK_INFO << "######### DEBUG START ############## " ; + MITK_INFO << "############################################## " ; + } + MITK_INFO << "############################################## " ; + MITK_INFO << "Name: " << pfc->getDisplayName(); + MITK_INFO << "iterationLevel: " << myLevel; + MITK_INFO << "CompositionType: " << pfc->getOperationType(); + MITK_INFO << "Number of children: " << pfc->getNumberOfChildren(); + + + //iterate through pfcs children + for(int i=0; igetNumberOfChildren(); ++i) + { + + mitk::PlanarFigure::Pointer tmpPFchild = pfc->getChildAt(i); + mitk::DataNode::Pointer savedPFchildNode = pfc->getDataNodeAt(i); + + if (tmpPFchild == savedPFchildNode->GetData()) + { + MITK_INFO << "[OK] Pointers point to same Data..."; + + }else{ + MITK_INFO << "Pointers differ in equation"; + } + + MITK_INFO << "Level: " << myLevel << " ChildNr.: " << i ; + + mitk::PlanarFigureComposite::Pointer pfcompcastNode= dynamic_cast(savedPFchildNode->GetData()); + mitk::PlanarFigureComposite::Pointer pfcompcast= dynamic_cast(tmpPFchild.GetPointer()); + if( !pfcompcast.IsNull() ) + { // we have a composite as child + + if ( pfcompcastNode.IsNull() ) + { + MITK_INFO << "************** NODE DIFFER FROM PFC...ERROR! ***************"; + } else { + MITK_INFO << "[OK]...node contains right type "; + } + + + + itLevelStatus++; + MITK_INFO << "child is PFC...debug this PFC"; + debugPFComposition(pfcompcast, itLevelStatus); + + } else { + + + // we have a planarFigure as child + // figure out which type + mitk::PlanarCircle::Pointer circleName = mitk::PlanarCircle::New(); + mitk::PlanarRectangle::Pointer rectName = mitk::PlanarRectangle::New(); + mitk::PlanarPolygon::Pointer polyName = mitk::PlanarPolygon::New(); + + + if (tmpPFchild->GetNameOfClass() == circleName->GetNameOfClass() ) + { + MITK_INFO << "a circle child of " << pfc->getDisplayName() ; + + } else if (tmpPFchild->GetNameOfClass() == rectName->GetNameOfClass() ){ + + MITK_INFO << "a rectangle child of " << pfc->getDisplayName() ; + + } else if (tmpPFchild->GetNameOfClass() == polyName->GetNameOfClass() ) { + + MITK_INFO << "a polygon child of " << pfc->getDisplayName() ; + } + + MITK_INFO << "....................................................... " ; + + + + + } + + } //end for + if (myLevel == 0) + { + MITK_INFO << "############################################## " ; + MITK_INFO << "######### DEBUG END ############## " ; + MITK_INFO << "############################################## " ; + } + + + +} + +void QmitkFiberBundleOperationsView::generatePFCompo_OR() +{ + mitk::PlanarFigureComposite::Pointer PFCOr = mitk::PlanarFigureComposite::New(); + PFCOr->setOperationType(mitk::PFCOMPOSITION_OR_OPERATION); + + + for( std::vector::iterator it = m_SelectedPF.begin(); + it != m_SelectedPF.end(); ++it ) + { + + mitk::DataNode::Pointer nodePF = *it; + mitk::PlanarFigure::Pointer tmpPF = dynamic_cast( nodePF->GetData() ); + PFCOr->addPlanarFigure( tmpPF ); + PFCOr->addDataNode( nodePF ); + PFCOr->setDisplayName("OR_COMPO"); + // MITK_INFO << "PFCAND(): added to AND PF" << nodePF->GetName(); + + } + + debugPFComposition(PFCOr, 0); + + + this->addPFCompositionToDataStorage(PFCOr, NULL /*parent*/); + +} + +void QmitkFiberBundleOperationsView::generatePFCompo_NOT() +{ + mitk::PlanarFigureComposite::Pointer PFCNot = mitk::PlanarFigureComposite::New(); + PFCNot->setOperationType(mitk::PFCOMPOSITION_NOT_OPERATION); + + + for( std::vector::iterator it = m_SelectedPF.begin(); + it != m_SelectedPF.end(); ++it ) + { + + mitk::DataNode::Pointer nodePF = *it; + mitk::PlanarFigure::Pointer tmpPF = dynamic_cast( nodePF->GetData() ); + PFCNot->addPlanarFigure( tmpPF ); + PFCNot->addDataNode( nodePF ); + PFCNot->setDisplayName("NOT_COMPO"); + // MITK_INFO << "PFCAND(): added to AND PF" << nodePF->GetName(); + + } + + debugPFComposition(PFCNot, 0); + + + this->addPFCompositionToDataStorage(PFCNot, NULL /*parent*/); + +} + +void QmitkFiberBundleOperationsView::deletePFCompo() +{ + +} + +void QmitkFiberBundleOperationsView::addPFCompositionToDataStorage(mitk::PlanarFigureComposite::Pointer pfcomp, mitk::DataNode::Pointer parentDataNode ) +{ + + //a new planarFigureComposition arrived + //convert it into a dataNode + mitk::DataNode::Pointer newPFCNode; + newPFCNode = mitk::DataNode::New(); + newPFCNode->SetName( pfcomp->getDisplayName() ); + //MITK_INFO << "PFComp Name: " << pfcomp->getDisplayName() << " newPFCNodeName: " << newPFCNode->GetName(); + newPFCNode->SetData(pfcomp); + newPFCNode->SetVisibility(true); + + switch (pfcomp->getOperationType()) { + case 0: + { + // AND PLANARFIGURECOMPOSITE + // newPFCNode->SetName("AND_PFCombo"); + + if (!parentDataNode.IsNull()) { + MITK_INFO << "adding " << newPFCNode->GetName() << " to " << parentDataNode->GetName() ; + GetDataStorage()->Add(newPFCNode, parentDataNode); + + } else { + MITK_INFO << "adding " << newPFCNode->GetName(); + GetDataStorage()->Add(newPFCNode); + + } + + + //iterate through its childs + + for(int i=0; igetNumberOfChildren(); ++i) + { + mitk::PlanarFigure::Pointer tmpPFchild = pfcomp->getChildAt(i); + mitk::DataNode::Pointer savedPFchildNode = pfcomp->getDataNodeAt(i); + + mitk::PlanarFigureComposite::Pointer pfcompcast= dynamic_cast(tmpPFchild.GetPointer()); + if ( !pfcompcast.IsNull() ) + { // child is of type planar Figure composite + // make new node of the child, cuz later the child has to be removed of its old position in datamanager + // feed new dataNode with information of the savedDataNode, which is gonna be removed soon + mitk::DataNode::Pointer newChildPFCNode; + newChildPFCNode = mitk::DataNode::New(); + newChildPFCNode->SetData(tmpPFchild); + newChildPFCNode->SetName( savedPFchildNode->GetName() ); + pfcompcast->setDisplayName( savedPFchildNode->GetName() ); //name might be changed in DataManager by user + + //update inside vector the dataNodePointer + pfcomp->replaceDataNodeAt(i, newChildPFCNode); + + addPFCompositionToDataStorage(pfcompcast, newPFCNode); //the current PFCNode becomes the childs parent + + + // remove savedNode here, cuz otherwise its children will change their position in the dataNodeManager + // without having its parent anymore + //GetDataStorage()->Remove(savedPFchildNode); + + if ( GetDataStorage()->Exists(savedPFchildNode)) { + MITK_INFO << savedPFchildNode->GetName() << " exists in DS...trying to remove it"; + + }else{ + MITK_INFO << "[ERROR] does NOT exist, but can I read its Name? " << savedPFchildNode->GetName(); + + } + // remove old child position in dataStorage + GetDataStorage()->Remove(savedPFchildNode); + + + if ( GetDataStorage()->Exists(savedPFchildNode)) { + MITK_INFO << savedPFchildNode->GetName() << " still exists"; + } + + + } else { + + // child is not of type PlanarFigureComposite, so its one of the planarFigures + // create new dataNode containing the data of the old dataNode, but position in dataManager will be + // modified cuz we re setting a (new) parent. + mitk::DataNode::Pointer newPFchildNode = mitk::DataNode::New(); + newPFchildNode->SetName(savedPFchildNode->GetName() ); + newPFchildNode->SetData(tmpPFchild); + newPFchildNode->SetVisibility(true); + + // replace the dataNode in PFComp DataNodeVector + pfcomp->replaceDataNodeAt(i, newPFchildNode); + + + if ( GetDataStorage()->Exists(savedPFchildNode)) { + MITK_INFO << savedPFchildNode->GetName() << " exists in DS...trying to remove it"; + + }else{ + MITK_INFO << "[ERROR] does NOT exist, but can I read its Name? " << savedPFchildNode->GetName(); + + } + // remove old child position in dataStorage + GetDataStorage()->Remove(savedPFchildNode); + + + if ( GetDataStorage()->Exists(savedPFchildNode)) { + MITK_INFO << savedPFchildNode->GetName() << " still exists"; + } + + MITK_INFO << "adding " << newPFchildNode->GetName() << " to " << newPFCNode->GetName(); + //add new child to datamanager with its new position as child of newPFCNode parent + GetDataStorage()->Add(newPFchildNode, newPFCNode); + + } + + + + + } + GetDataStorage()->Modified(); + + + + + break; + } + case 1: + { + // AND PLANARFIGURECOMPOSITE + // newPFCNode->SetName("AND_PFCombo"); + + if (!parentDataNode.IsNull()) { + MITK_INFO << "adding " << newPFCNode->GetName() << " to " << parentDataNode->GetName() ; + GetDataStorage()->Add(newPFCNode, parentDataNode); + + } else { + MITK_INFO << "adding " << newPFCNode->GetName(); + GetDataStorage()->Add(newPFCNode); + + } + + + //iterate through its childs + + for(int i=0; igetNumberOfChildren(); ++i) + { + mitk::PlanarFigure::Pointer tmpPFchild = pfcomp->getChildAt(i); + mitk::DataNode::Pointer savedPFchildNode = pfcomp->getDataNodeAt(i); + + mitk::PlanarFigureComposite::Pointer pfcompcast= dynamic_cast(tmpPFchild.GetPointer()); + if ( !pfcompcast.IsNull() ) + { // child is of type planar Figure composite + // make new node of the child, cuz later the child has to be removed of its old position in datamanager + // feed new dataNode with information of the savedDataNode, which is gonna be removed soon + mitk::DataNode::Pointer newChildPFCNode; + newChildPFCNode = mitk::DataNode::New(); + newChildPFCNode->SetData(tmpPFchild); + newChildPFCNode->SetName( savedPFchildNode->GetName() ); + pfcompcast->setDisplayName( savedPFchildNode->GetName() ); //name might be changed in DataManager by user + + //update inside vector the dataNodePointer + pfcomp->replaceDataNodeAt(i, newChildPFCNode); + + addPFCompositionToDataStorage(pfcompcast, newPFCNode); //the current PFCNode becomes the childs parent + + + // remove savedNode here, cuz otherwise its children will change their position in the dataNodeManager + // without having its parent anymore + //GetDataStorage()->Remove(savedPFchildNode); + + if ( GetDataStorage()->Exists(savedPFchildNode)) { + MITK_INFO << savedPFchildNode->GetName() << " exists in DS...trying to remove it"; + + }else{ + MITK_INFO << "[ERROR] does NOT exist, but can I read its Name? " << savedPFchildNode->GetName(); + + } + // remove old child position in dataStorage + GetDataStorage()->Remove(savedPFchildNode); + + + if ( GetDataStorage()->Exists(savedPFchildNode)) { + MITK_INFO << savedPFchildNode->GetName() << " still exists"; + } + + + } else { + + // child is not of type PlanarFigureComposite, so its one of the planarFigures + // create new dataNode containing the data of the old dataNode, but position in dataManager will be + // modified cuz we re setting a (new) parent. + mitk::DataNode::Pointer newPFchildNode = mitk::DataNode::New(); + newPFchildNode->SetName(savedPFchildNode->GetName() ); + newPFchildNode->SetData(tmpPFchild); + newPFchildNode->SetVisibility(true); + + // replace the dataNode in PFComp DataNodeVector + pfcomp->replaceDataNodeAt(i, newPFchildNode); + + + if ( GetDataStorage()->Exists(savedPFchildNode)) { + MITK_INFO << savedPFchildNode->GetName() << " exists in DS...trying to remove it"; + + }else{ + MITK_INFO << "[ERROR] does NOT exist, but can I read its Name? " << savedPFchildNode->GetName(); + + } + // remove old child position in dataStorage + GetDataStorage()->Remove(savedPFchildNode); + + + if ( GetDataStorage()->Exists(savedPFchildNode)) { + MITK_INFO << savedPFchildNode->GetName() << " still exists"; + } + + MITK_INFO << "adding " << newPFchildNode->GetName() << " to " << newPFCNode->GetName(); + //add new child to datamanager with its new position as child of newPFCNode parent + GetDataStorage()->Add(newPFchildNode, newPFCNode); + + } + + + + + } + GetDataStorage()->Modified(); + + + + + break; + + } + case 2: + { + // AND PLANARFIGURECOMPOSITE + // newPFCNode->SetName("AND_PFCombo"); + + if (!parentDataNode.IsNull()) { + MITK_INFO << "adding " << newPFCNode->GetName() << " to " << parentDataNode->GetName() ; + GetDataStorage()->Add(newPFCNode, parentDataNode); + + } else { + MITK_INFO << "adding " << newPFCNode->GetName(); + GetDataStorage()->Add(newPFCNode); + + } + + + //iterate through its childs + + for(int i=0; igetNumberOfChildren(); ++i) + { + mitk::PlanarFigure::Pointer tmpPFchild = pfcomp->getChildAt(i); + mitk::DataNode::Pointer savedPFchildNode = pfcomp->getDataNodeAt(i); + + mitk::PlanarFigureComposite::Pointer pfcompcast= dynamic_cast(tmpPFchild.GetPointer()); + if ( !pfcompcast.IsNull() ) + { // child is of type planar Figure composite + // makeRemoveBundle new node of the child, cuz later the child has to be removed of its old position in datamanager + // feed new dataNode with information of the savedDataNode, which is gonna be removed soon + mitk::DataNode::Pointer newChildPFCNode; + newChildPFCNode = mitk::DataNode::New(); + newChildPFCNode->SetData(tmpPFchild); + newChildPFCNode->SetName( savedPFchildNode->GetName() ); + pfcompcast->setDisplayName( savedPFchildNode->GetName() ); //name might be changed in DataManager by user + + //update inside vector the dataNodePointer + pfcomp->replaceDataNodeAt(i, newChildPFCNode); + + addPFCompositionToDataStorage(pfcompcast, newPFCNode); //the current PFCNode becomes the childs parent + + + // remove savedNode here, cuz otherwise its children will change their position in the dataNodeManager + // without having its parent anymore + //GetDataStorage()->Remove(savedPFchildNode); + + if ( GetDataStorage()->Exists(savedPFchildNode)) { + MITK_INFO << savedPFchildNode->GetName() << " exists in DS...trying to remove it"; + + }else{ + MITK_INFO << "[ERROR] does NOT exist, but can I read its Name? " << savedPFchildNode->GetName(); + } + // remove old child position in dataStorage + GetDataStorage()->Remove(savedPFchildNode); + + + if ( GetDataStorage()->Exists(savedPFchildNode)) { + MITK_INFO << savedPFchildNode->GetName() << " still exists"; + } + + + } else { + + // child is not of type PlanarFigureComposite, so its one of the planarFigures + // create new dataNode containing the data of the old dataNode, but position in dataManager will be + // modified cuz we re setting a (new) parent. + mitk::DataNode::Pointer newPFchildNode = mitk::DataNode::New(); + newPFchildNode->SetName(savedPFchildNode->GetName() ); + newPFchildNode->SetData(tmpPFchild); + newPFchildNode->SetVisibility(true); + + // replace the dataNode in PFComp DataNodeVector + pfcomp->replaceDataNodeAt(i, newPFchildNode); + + + if ( GetDataStorage()->Exists(savedPFchildNode)) { + MITK_INFO << savedPFchildNode->GetName() << " exists in DS...trying to remove it"; + + }else{ + MITK_INFO << "[ERROR] does NOT exist, but can I read its Name? " << savedPFchildNode->GetName(); + + } + // remove old child position in dataStorage + GetDataStorage()->Remove(savedPFchildNode); + + + if ( GetDataStorage()->Exists(savedPFchildNode)) { + MITK_INFO << savedPFchildNode->GetName() << " still exists"; + } + + MITK_INFO << "adding " << newPFchildNode->GetName() << " to " << newPFCNode->GetName(); + //add new child to datamanager with its new position as child of newPFCNode parent + GetDataStorage()->Add(newPFchildNode, newPFCNode); + + } + + + + + } + GetDataStorage()->Modified(); + + + + + break; + } + default: + MITK_INFO << "we have an UNDEFINED composition... ERROR" ; + break; + } + + + + + + +} + + +void QmitkFiberBundleOperationsView::JoinBundles() +{ + mitk::FiberBundle::Pointer bundle1 = dynamic_cast(m_SelectedFB.at(0)->GetData()); + mitk::FiberBundle::Pointer bundle2 = dynamic_cast(m_SelectedFB.at(1)->GetData()); + + mitk::FiberBundle::Pointer newBundle = bundle1->JoinBundle(bundle2); + mitk::DataNode::Pointer fbNode = mitk::DataNode::New(); + fbNode->SetData(newBundle); + fbNode->SetName(m_SelectedFB.at(0)->GetName()+"+"+m_SelectedFB.at(1)->GetName()); + fbNode->SetVisibility(true); + GetDataStorage()->Add(fbNode); +} + +void QmitkFiberBundleOperationsView::SubstractBundles() +{ + mitk::FiberBundle::Pointer bundle1 = dynamic_cast(m_SelectedFB.at(0)->GetData()); + mitk::FiberBundle::Pointer bundle2 = dynamic_cast(m_SelectedFB.at(1)->GetData()); + + mitk::FiberBundle::Pointer newBundle = bundle1->SubstractBundle(bundle2); + mitk::DataNode::Pointer fbNode = mitk::DataNode::New(); + fbNode->SetData(newBundle); + fbNode->SetName(m_SelectedFB.at(0)->GetName()+"-"+m_SelectedFB.at(1)->GetName()); + fbNode->SetVisibility(true); + GetDataStorage()->Add(fbNode); +} + +void QmitkFiberBundleOperationsView::GenerationStart() +{ + int generationMethod = m_Controls->m_GenerationBox->currentIndex(); + + std::vector nodes = GetDataManagerSelection(); + if (nodes.empty()){ + QMessageBox::information( NULL, "Warning", "No data object selected!"); + MITK_WARN("QmitkFiberBundleOperationsView") << "no data object selected"; + return; + } + + for( std::vector::iterator it = nodes.begin(); it != nodes.end(); ++it ) + { + mitk::DataNode::Pointer node = *it; + + if (node.IsNotNull() && dynamic_cast(node->GetData())) + { + m_FiberBundle = dynamic_cast(node->GetData()); + m_FiberBundleNode = node; + switch(generationMethod){ + case 0: + GenerateGreyscaleHeatmap(true); + break; + case 1: + GenerateGreyscaleHeatmap(false); + break; + case 2: + GenerateColorHeatmap(); + break; + case 3: + GenerateFiberEndingsImage(); + break; + case 4: + GenerateFiberEndingsPointSet(); + break; + case 5: + DWIGenerationStart(); + break; + } + } + } +} + +// generate pointset displaying the fiber endings +void QmitkFiberBundleOperationsView::GenerateFiberEndingsPointSet() +{ + if(m_FiberBundle.IsNull()){ + QMessageBox::information( NULL, "Warning", "No fiber bundle selected!"); + MITK_WARN("QmitkGlobalFiberTrackingView") << "no fiber bundle selected"; + return; + } + + mitk::Geometry3D::Pointer geometry = m_FiberBundle->GetGeometry(); + + mitk::PointSet::Pointer pointSet = mitk::PointSet::New(); + + int numTracts = m_FiberBundle->GetNumTracts(); + int count = 0; + for( int i=0; iGetNumPoints(i); + itk::Point start = m_FiberBundle->GetPoint(i,0); + itk::Point world1; + geometry->IndexToWorld(start, world1); + pointSet->InsertPoint(count, world1); + count++; + // get fiber end point + if(numVertices>1) + { + itk::Point end = m_FiberBundle->GetPoint(i,numVertices-1); + itk::Point world; + geometry->IndexToWorld(end, world); + pointSet->InsertPoint(count, world); + count++; + } + } + + mitk::DataNode::Pointer pointSetNode = mitk::DataNode::New(); + pointSetNode->SetData( pointSet ); + QString name(m_FiberBundleNode->GetName().c_str()); + name += "_fiber_endings"; + pointSetNode->SetName(name.toStdString()); + pointSetNode->SetProperty( "opacity", mitk::FloatProperty::New( 1 ) ); + pointSetNode->SetProperty( "pointsize", mitk::FloatProperty::New( 0.3) ); + pointSetNode->SetColor( 1.0, 1.0, 1.0 ); + + GetDefaultDataStorage()->Add(pointSetNode); +} + +// generate image displaying the fiber endings +void QmitkFiberBundleOperationsView::GenerateFiberEndingsImage() +{ + if(m_FiberBundle.IsNull()){ + QMessageBox::information( NULL, "Warning", "No fiber bundle selected!"); + MITK_WARN("QmitkGlobalFiberTrackingView") << "no fiber bundle selected"; + return; + } + + typedef itk::RGBAPixel OutPixType; + + // run generator + typedef itk::TractsToFiberEndingsImageFilter ImageGeneratorType; + ImageGeneratorType::Pointer generator = ImageGeneratorType::New(); + generator->SetFiberBundle(m_FiberBundle); + + generator->SetUpsamplingFactor(m_Controls->m_UpsamplingSpinBox->value()); + generator->Update(); + + // get result + typedef itk::Image OutType; + OutType::Pointer outImg = generator->GetOutput(); + + mitk::Image::Pointer img = mitk::Image::New(); + img->InitializeByItk(outImg.GetPointer()); + img->SetVolume(outImg->GetBufferPointer()); + + // to datastorage + mitk::DataNode::Pointer node = mitk::DataNode::New(); + node->SetData(img); + QString name(m_FiberBundleNode->GetName().c_str()); + name += "_fiber_endings"; + node->SetName(name.toStdString()); + node->SetVisibility(true); + + mitk::LevelWindow opaclevwin; + opaclevwin.SetRangeMinMax(0,255); + opaclevwin.SetWindowBounds(0,0); + mitk::LevelWindowProperty::Pointer prop = mitk::LevelWindowProperty::New(opaclevwin); + node->AddProperty( "opaclevelwindow", prop ); + + GetDataStorage()->Add(node); +} + +// generate rgba heatmap from fiber bundle +void QmitkFiberBundleOperationsView::GenerateColorHeatmap() +{ + if(m_FiberBundle.IsNull() || m_FiberBundleNode.IsNull()) + { + QMessageBox::information( NULL, "Warning", "No fiber bundle selected!"); + MITK_WARN("QmitkGlobalFiberTrackingView") << "no fiber bundle selected"; + return; + } + + typedef itk::RGBAPixel OutPixType; + + // run generator + typedef itk::TractsToProbabilityImageFilter + ImageGeneratorType; + ImageGeneratorType::Pointer generator = ImageGeneratorType::New(); + //generator->SetInput(NULL); + generator->SetFiberBundle(m_FiberBundle); + + generator->SetUpsamplingFactor(m_Controls->m_UpsamplingSpinBox->value()); + generator->Update(); + + // get result + typedef itk::Image OutType; + OutType::Pointer outImg = generator->GetOutput(); + + mitk::Image::Pointer img = mitk::Image::New(); + img->InitializeByItk(outImg.GetPointer()); + img->SetVolume(outImg->GetBufferPointer()); + + // to datastorage + mitk::DataNode::Pointer node = mitk::DataNode::New(); + node->SetData(img); + QString name(m_FiberBundleNode->GetName().c_str()); + name += "_rgba_heatmap"; + node->SetName(name.toStdString()); + node->SetVisibility(true); + + mitk::LevelWindow opaclevwin; + opaclevwin.SetRangeMinMax(0,255); + opaclevwin.SetWindowBounds(0,0); + mitk::LevelWindowProperty::Pointer prop = + mitk::LevelWindowProperty::New(opaclevwin); + node->AddProperty( "opaclevelwindow", prop ); + + GetDataStorage()->Add(node); +} + +// generate greyscale heatmap from fiber bundle +void QmitkFiberBundleOperationsView::GenerateGreyscaleHeatmap(bool binary) +{ + if(m_FiberBundle.IsNull() || m_FiberBundleNode.IsNull()) + { + QMessageBox::information( NULL, "Warning", "No fiber bundle selected!"); + MITK_WARN("QmitkGlobalFiberTrackingView") << "no fiber bundle selected"; + return; + } + + typedef unsigned char OutPixType; + + // run generator + typedef itk::TractsToProbabilityImageFilter ImageGeneratorType; + ImageGeneratorType::Pointer generator = ImageGeneratorType::New(); + generator->SetFiberBundle(m_FiberBundle); + generator->SetInvertImage(m_Controls->m_InvertCheckbox->isChecked()); + generator->SetUpsamplingFactor(m_Controls->m_UpsamplingSpinBox->value()); + if (binary) + generator->SetBinaryEnvelope(true); + else + generator->SetBinaryEnvelope(false); + generator->Update(); + + // get result + typedef itk::Image OutType; + OutType::Pointer outImg = generator->GetOutput(); + + mitk::Image::Pointer img = mitk::Image::New(); + img->InitializeByItk(outImg.GetPointer()); + img->SetVolume(outImg->GetBufferPointer()); + + // to datastorage + mitk::DataNode::Pointer node = mitk::DataNode::New(); + node->SetData(img); + QString name(m_FiberBundleNode->GetName().c_str()); + if(binary) + name += "_enveloppe"; + else + name += "_heatmap"; + node->SetName(name.toStdString()); + node->SetVisibility(true); + + mitk::LevelWindow opaclevwin2; + opaclevwin2.SetRangeMinMax(0,255); + opaclevwin2.SetWindowBounds(0,0); + mitk::LevelWindowProperty::Pointer prop2 = + mitk::LevelWindowProperty::New(opaclevwin2); + node->AddProperty( "opaclevelwindow", prop2 ); + + GetDataStorage()->Add(node); +} + +// generate dwi from fiber bundle (experimental) +void QmitkFiberBundleOperationsView::DWIGenerationStart() +{ + // get fiber bundle and dwi image from data manager + typedef mitk::DiffusionImage DiffVolumesType; + std::vector nodes = GetDataManagerSelection(); + DiffVolumesType::Pointer originalDWI = NULL; + mitk::FiberBundle::Pointer fiberBundle = NULL; + for( std::vector::iterator it = nodes.begin(); it != nodes.end(); ++it ) + { + mitk::DataNode::Pointer node = *it; + + if (node.IsNotNull() && + dynamic_cast(node->GetData())) + { + originalDWI = dynamic_cast(node->GetData()); + continue; + } + if (node.IsNotNull() && + dynamic_cast(node->GetData())) + { + fiberBundle = dynamic_cast(node->GetData()); + } + } + if(fiberBundle.IsNull() || originalDWI.IsNull()){ + QMessageBox::information( NULL, "Warning", "Please load and select a dwi image and a fiber bundle."); + MITK_WARN("QmitkGlobalFiberTrackingView") << "please select a fiber bundle and a diffusion image"; + return; + } + + // CONSTRUCT CONTAINER WITH DIRECTIONS + typedef vnl_vector_fixed< double, 3 > GradientDirectionType; + typedef itk::VectorContainer< unsigned int, GradientDirectionType > GradientDirectionContainerType; + GradientDirectionContainerType::Pointer directions = originalDWI->GetDirections(); + + float bVal = originalDWI->GetB_Value(); + + typedef itk::VectorImage< short, 3 > DWIImageType; + DWIImageType::Pointer vectorImage = DWIImageType::New(); + itk::TractsToDWIImageFilter::Pointer filter = itk::TractsToDWIImageFilter::New(); + + filter->SetInput(originalDWI->GetVectorImage()); + filter->SetFiberBundle(fiberBundle); + filter->SetbD(m_Controls->m_UpsamplingSpinBox->value()); + filter->SetGradientDirections(directions); + filter->SetParticleWidth(0.2); + filter->GenerateData(); + vectorImage = filter->GetOutput(); + + DiffVolumesType::Pointer diffImage = DiffVolumesType::New(); + diffImage->SetDirections(directions); + diffImage->SetOriginalDirections(directions); + diffImage->SetVectorImage(vectorImage); + diffImage->SetB_Value(bVal); + diffImage->InitializeFromVectorImage(); + + mitk::DataNode::Pointer node = mitk::DataNode::New(); + node->SetData( diffImage ); + QString name(m_FiberBundleNode->GetName().c_str()); + name += "_dwi"; + GetDefaultDataStorage()->Add(node); +} + diff --git a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberBundleOperationsView.h b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberBundleOperationsView.h new file mode 100644 index 0000000000..c82de84190 --- /dev/null +++ b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberBundleOperationsView.h @@ -0,0 +1,191 @@ +/*========================================================================= + +Program: Medical Imaging & Interaction Toolkit +Language: C++ +Date: $Date: 2010-03-31 16:40:27 +0200 (Mi, 31 Mrz 2010) $ +Version: $Revision: 21975 $ + +Copyright (c) German Cancer Research Center, Division of Medical and +Biological Informatics. All rights reserved. +See MITKCopyright.txt or http://www.mitk.org/copyright.html for details. + +This software is distributed WITHOUT ANY WARRANTY; without even +the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR +PURPOSE. See the above copyright notices for more information. + +=========================================================================*/ + +#ifndef QmitkFiberBundleOperationsView_h +#define QmitkFiberBundleOperationsView_h + +#include +#include + +#include +#include "ui_QmitkFiberBundleOperationsViewControls.h" + +#include "mitkDataStorage.h" +#include "mitkDataStorageSelection.h" + +#include "mitkPlanarFigure.h" +#include "mitkFiberBundle.h" +#include "mitkPlanarFigureComposite.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + + +/*! +\brief QmitkFiberBundleView + +\warning This application module is not yet documented. Use "svn blame/praise/annotate" and ask the author to provide basic documentation. + +\sa QmitkFunctionality +\ingroup Functionalities +*/ +class QmitkFiberBundleOperationsView : public QmitkFunctionality +{ + + + // 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 > MaskImage3DType; + typedef itk::Image< float, 3 > FloatImageType; + + static const std::string VIEW_ID; + + QmitkFiberBundleOperationsView(); + virtual ~QmitkFiberBundleOperationsView(); + + virtual void CreateQtPartControl(QWidget *parent); + + virtual void StdMultiWidgetAvailable (QmitkStdMultiWidget &stdMultiWidget); + virtual void StdMultiWidgetNotAvailable(); + virtual void Activated(); + + protected slots: + + /// \brief Called when the user clicks the GUI button + void ActionDrawEllipseTriggered(); + void ActionDrawPolygonTriggered(); + void DoFiberExtraction(); + void generatePFCompo_AND(); + void generatePFCompo_OR(); + void generatePFCompo_NOT(); + void deletePFCompo(); + + void JoinBundles(); + void SubstractBundles(); + void GenerateROIImage(); + void GenerationStart(); + + virtual void AddFigureToDataStorage(mitk::PlanarFigure* figure, const QString& name, + const char *propertyKey = NULL, mitk::BaseProperty *property = NULL ); + + +protected: + + /// \brief called by QmitkFunctionality when DataManager's selection has changed + virtual void OnSelectionChanged( std::vector nodes ); + + Ui::QmitkFiberBundleOperationsViewControls* m_Controls; + + QmitkStdMultiWidget* m_MultiWidget; + + //void Select( mitk::DataNode::Pointer node, bool clearMaskOnFirstArgNULL=false, bool clearImageOnFirstArgNULL=false ); + + /** 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 ); + + template < typename TPixel, unsigned int VImageDimension > + void InternalReorientImagePlane( + const itk::Image< TPixel, VImageDimension > *image, mitk::Geometry3D* imggeo, mitk::Geometry3D* planegeo3D, int additionalIndex ); + +private: + + int m_EllipseCounter; + int m_PolygonCounter; + //contains the selected FiberBundles + std::vector m_SelectedFB; + + //contains the selected PlanarFigures + std::vector m_SelectedPF; + + mitk::Image::ConstPointer m_Image; + mitk::Image::Pointer m_InternalImage; + mitk::PlanarFigure::Pointer m_PlanarFigure; + float m_UpsamplingFactor; + MaskImage3DType::Pointer m_InternalImageMask3D; + MaskImage3DType::Pointer m_PlanarFigureImage; + mitk::FiberBundle::Pointer m_FiberBundle; + mitk::DataNode::Pointer m_FiberBundleNode; + + void addPFCompositionToDataStorage(mitk::PlanarFigureComposite::Pointer, mitk::DataNode::Pointer); + void debugPFComposition(mitk::PlanarFigureComposite::Pointer , int ); + void CompositeExtraction(mitk::PlanarFigure::Pointer planarFigure, mitk::Image* image); + void GenerateGreyscaleHeatmap(bool binary); + void GenerateColorHeatmap(); + void GenerateFiberEndingsImage(); + void GenerateFiberEndingsPointSet(); + void DWIGenerationStart(); +}; + + + +#endif // _QMITKFIBERTRACKINGVIEW_H_INCLUDED + diff --git a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberBundleOperationsViewControls.ui b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberBundleOperationsViewControls.ui new file mode 100644 index 0000000000..ef94229ae0 --- /dev/null +++ b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberBundleOperationsViewControls.ui @@ -0,0 +1,523 @@ + + + QmitkFiberBundleOperationsViewControls + + + + 0 + 0 + 665 + 587 + + + + Form + + + + 0 + + + + + + 0 + 0 + + + + + 200 + 50 + + + + + 300 + 60 + + + + QFrame::NoFrame + + + QFrame::Raised + + + + 12 + + + 0 + + + + + + 30 + 30 + + + + + + + + :/QmitkDiffusionImaging/circle.png:/QmitkDiffusionImaging/circle.png + + + + 32 + 32 + + + + false + + + true + + + + + + + + 30 + 30 + + + + + + + + :/QmitkDiffusionImaging/rectangle.png:/QmitkDiffusionImaging/rectangle.png + + + + 32 + 32 + + + + true + + + true + + + + + + + + 30 + 30 + + + + + + + + :/QmitkDiffusionImaging/polygon.png:/QmitkDiffusionImaging/polygon.png + + + + 32 + 32 + + + + true + + + true + + + + + + + Qt::Horizontal + + + + 40 + 20 + + + + + + + + + + + + 16777215 + 16777215 + + + + QFrame::NoFrame + + + QFrame::Raised + + + + 0 + + + 0 + + + + + false + + + + 60 + 16777215 + + + + AND + + + + + + + false + + + + 60 + 16777215 + + + + OR + + + + + + + false + + + + 60 + 16777215 + + + + NOT + + + + + + + false + + + + 60 + 16777215 + + + + DEL + + + + + + + + + + false + + + + 0 + 0 + + + + + 200 + 16777215 + + + + + 11 + + + + Extract Fibers from ROI + + + + + + + false + + + + 0 + 0 + + + + + 300 + 16777215 + + + + + 11 + + + + Generates One Binary Image Per Planar Figure And One Binary Image Containing All Planar Figures + + + Generate Binary ROI Images + + + + + + + QFrame::NoFrame + + + QFrame::Raised + + + + 0 + + + 0 + + + + + false + + + + 0 + 0 + + + + + 200 + 16777215 + + + + + 11 + + + + Join Bundles + + + + + + + false + + + + 0 + 0 + + + + + 200 + 16777215 + + + + + 11 + + + + Substract Bundles + + + + + + + + + + + 0 + 0 + + + + QFrame::StyledPanel + + + QFrame::Raised + + + + + + Num. Fibers: + + + + + + + 0 + + + + + + + + + + QFrame::StyledPanel + + + QFrame::Raised + + + + + + + 0 + 0 + + + + + Generate Fiber Bundle Enveloppe + + + + + Generate 2D Fiber Bundle Image (Greyscale) + + + + + Generate 2D Fiber Bundle Image (Color) + + + + + Generate 2D Fiber Endings Image + + + + + Generate Fiber Endings Pointset + + + + + Generate DWI from Fiber Bundle (experimental) + + + + + + + + 1 + + + 10 + + + 4 + + + + + + + true + + + + 0 + 0 + + + + + 200 + 16777215 + + + + + 11 + + + + Generate + + + + + + + Invert + + + + + + + + + + Qt::Vertical + + + + 20 + 40 + + + + + + + + + + + diff --git a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/mitkPluginActivator.cpp b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/mitkPluginActivator.cpp index 1715b75d19..115a4b9f1e 100644 --- a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/mitkPluginActivator.cpp +++ b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/mitkPluginActivator.cpp @@ -1,40 +1,42 @@ #include "mitkPluginActivator.h" #include #include "src/internal/QmitkDiffusionImagingPublicPerspective.h" #include "src/internal/QmitkQBallReconstructionView.h" #include "src/internal/QmitkPreprocessingView.h" #include "src/internal/QmitkDiffusionDicomImportView.h" #include "src/internal/QmitkDiffusionQuantificationView.h" #include "src/internal/QmitkTensorReconstructionView.h" #include "src/internal/QmitkControlVisualizationPropertiesView.h" #include "src/internal/QmitkODFDetailsView.h" #include "src/internal/QmitkGlobalFiberTrackingView.h" +#include "src/internal/QmitkFiberBundleOperationsView.h" namespace mitk { void PluginActivator::start(ctkPluginContext* context) { BERRY_REGISTER_EXTENSION_CLASS(QmitkDiffusionImagingPublicPerspective, context) BERRY_REGISTER_EXTENSION_CLASS(QmitkQBallReconstructionView, context) BERRY_REGISTER_EXTENSION_CLASS(QmitkPreprocessingView, context) BERRY_REGISTER_EXTENSION_CLASS(QmitkDiffusionDicomImport, context) BERRY_REGISTER_EXTENSION_CLASS(QmitkDiffusionQuantificationView, context) BERRY_REGISTER_EXTENSION_CLASS(QmitkTensorReconstructionView, context) BERRY_REGISTER_EXTENSION_CLASS(QmitkControlVisualizationPropertiesView, context) BERRY_REGISTER_EXTENSION_CLASS(QmitkODFDetailsView, context) BERRY_REGISTER_EXTENSION_CLASS(QmitkGlobalFiberTrackingView, context) + BERRY_REGISTER_EXTENSION_CLASS(QmitkFiberBundleOperationsView, context) } void PluginActivator::stop(ctkPluginContext* context) { Q_UNUSED(context) } } Q_EXPORT_PLUGIN2(org_mitk_gui_qt_diffusionimaging, mitk::PluginActivator) diff --git a/Modules/DiffusionImaging/Algorithms/itkGaussianInterpolateImageFunction.h b/Modules/DiffusionImaging/Algorithms/itkGaussianInterpolateImageFunction.h new file mode 100644 index 0000000000..80ba244874 --- /dev/null +++ b/Modules/DiffusionImaging/Algorithms/itkGaussianInterpolateImageFunction.h @@ -0,0 +1,180 @@ +/*========================================================================= + + Program: Insight Segmentation & Registration Toolkit + Module: $RCSfile: itkGaussianInterpolateImageFunction.h,v $ + Language: C++ + Date: $Date: $ + Version: $Revision: $ + + Copyright (c) Insight Software Consortium. All rights reserved. + See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details. + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notices for more information. + +=========================================================================*/ +#ifndef __itkGaussianInterpolateImageFunction_h +#define __itkGaussianInterpolateImageFunction_h + +#include "itkInterpolateImageFunction.h" + +#include "itkConceptChecking.h" +#include "itkFixedArray.h" +#include "vnl/vnl_erf.h" + +namespace itk +{ + +/** \class GaussianInterpolateImageFunction + * \brief Evaluates the Gaussian interpolation of an image. + * + * This class defines an N-dimensional Gaussian interpolation function using + * the vnl error function. The two parameters associated with this function + * are: + * 1. Sigma - a scalar array of size ImageDimension determining the width + * of the interpolation function. + * 2. Alpha - a scalar specifying the cutoff distance over which the function + * is calculated. + * + * \ingroup ImageFunctions ImageInterpolators + */ + +template +class ITK_EXPORT GaussianInterpolateImageFunction : + public InterpolateImageFunction +{ +public: + /** Standard class typedefs. */ + typedef GaussianInterpolateImageFunction Self; + typedef InterpolateImageFunction Superclass; + typedef SmartPointer Pointer; + typedef SmartPointer ConstPointer; + + /** Run-time type information (and related methods). */ + itkTypeMacro( GaussianInterpolateImageFunction, InterpolateImageFunction ); + + /** Method for creation through the object factory. */ + itkNewMacro( Self ); + + /** ImageDimension constant */ + itkStaticConstMacro( ImageDimension, unsigned int, + TInputImage::ImageDimension ); + + + /** OutputType typedef support. */ + typedef typename Superclass::OutputType OutputType; + + /** InputImageType typedef support. */ + typedef typename Superclass::InputImageType InputImageType; + + /** RealType typedef support. */ + typedef typename Superclass::RealType RealType; + + /** Index typedef support. */ + typedef typename Superclass::IndexType IndexType; + + /** ContinuousIndex typedef support. */ + typedef typename Superclass::ContinuousIndexType ContinuousIndexType; + + /** Array typedef support */ + typedef FixedArray ArrayType; + + /** + * Set input image + */ + virtual void SetInputImage( const TInputImage *image ) + { + Superclass::SetInputImage( image ); + this->ComputeBoundingBox(); + } + + /** + * Set/Get ivars + */ + virtual void SetSigma( const ArrayType s ) + { + itkDebugMacro( "setting Sigma to " << s ); + if( this->m_Sigma != s ) + { + this->m_Sigma = s; + this->ComputeBoundingBox(); + this->Modified(); + } + } + itkGetConstMacro( Sigma, ArrayType ); + void SetSigma( RealType *s ) + { + ArrayType sigma; + for( unsigned int d = 0; d < ImageDimension; d++ ) + { + sigma[d] = s[d]; + } + this->SetSigma( sigma ); + } + + virtual void SetAlpha( const RealType a ) + { + itkDebugMacro( "setting Alpha to " << a ); + if( this->m_Alpha != a ) + { + this->m_Alpha = a; + this->ComputeBoundingBox(); + this->Modified(); + } + } + itkGetConstMacro( Alpha, RealType ); + + void SetParameters( RealType *sigma, RealType alpha ) + { + this->SetSigma( sigma ); + this->SetAlpha( alpha ); + } + + /** + * Evaluate at the given index + */ + virtual OutputType EvaluateAtContinuousIndex( + const ContinuousIndexType & cindex ) const + { + return this->EvaluateAtContinuousIndex( cindex, NULL ); + } + + /** + * Evaluate function value and gradient at the given index + */ + virtual OutputType EvaluateAtContinuousIndex( + const ContinuousIndexType &, OutputType * ) const; + +protected: + GaussianInterpolateImageFunction(); + ~GaussianInterpolateImageFunction(){}; + void PrintSelf( std::ostream& os, Indent indent ) const; + +private: + GaussianInterpolateImageFunction( const Self& ); //purposely not implemented + void operator=( const Self& ); //purposely not implemented + + void ComputeBoundingBox(); + + void ComputeErrorFunctionArray( unsigned int dimension, RealType cindex, + vnl_vector &erfArray, vnl_vector &gerfArray, + bool evaluateGradient = false ) const; + + ArrayType m_Sigma; + RealType m_Alpha; + + ArrayType m_BoundingBoxStart; + ArrayType m_BoundingBoxEnd; + ArrayType m_ScalingFactor; + ArrayType m_CutoffDistance; +}; + +} // end namespace itk + +#ifndef ITK_MANUAL_INSTANTIATION +#include "itkGaussianInterpolateImageFunction.txx" +#endif + +#endif diff --git a/Modules/DiffusionImaging/Algorithms/itkGaussianInterpolateImageFunction.txx b/Modules/DiffusionImaging/Algorithms/itkGaussianInterpolateImageFunction.txx new file mode 100644 index 0000000000..5fe9c588b1 --- /dev/null +++ b/Modules/DiffusionImaging/Algorithms/itkGaussianInterpolateImageFunction.txx @@ -0,0 +1,228 @@ +/*========================================================================= + + Program: Insight Segmentation & Registration Toolkit + Module: $RCSfile: itkGaussianInterpolateImageFunction.txx,v $ + Language: C++ + Date: $Date: $ + Version: $Revision: $ + + Copyright (c) Insight Software Consortium. All rights reserved. + See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details. + + Portions of this code are covered under the VTK copyright. + See VTKCopyright.txt or http://www.kitware.com/VTKCopyright.htm for details. + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notices for more information. + +=========================================================================*/ +#ifndef __itkGaussianInterpolateImageFunction_txx +#define __itkGaussianInterpolateImageFunction_txx + +#include "itkGaussianInterpolateImageFunction.h" + +#include "itkImageRegionConstIteratorWithIndex.h" + +namespace itk +{ + +/** + * Constructor + */ +template +GaussianInterpolateImageFunction +::GaussianInterpolateImageFunction() +{ + this->m_Alpha = 1.0; + this->m_Sigma.Fill( 1.0 ); +} + +/** + * Standard "PrintSelf" method + */ +template +void +GaussianInterpolateImageFunction +::PrintSelf( std::ostream& os, Indent indent ) const +{ + Superclass::PrintSelf( os, indent ); + os << indent << "Alpha: " << this->m_Alpha << std::endl; + os << indent << "Sigma: " << this->m_Sigma << std::endl; +} + +template +void +GaussianInterpolateImageFunction +::ComputeBoundingBox() +{ + if( !this->GetInputImage() ) + { + return; + } + + for( unsigned int d = 0; d < ImageDimension; d++ ) + { + this->m_BoundingBoxStart[d] = -0.5; + this->m_BoundingBoxEnd[d] = static_cast( + this->GetInputImage()->GetBufferedRegion().GetSize()[d] ) - 0.5; + this->m_ScalingFactor[d] = 1.0 / ( vnl_math::sqrt2 * this->m_Sigma[d] / + this->GetInputImage()->GetSpacing()[d] ); + this->m_CutoffDistance[d] = this->m_Sigma[d] * this->m_Alpha / + this->GetInputImage()->GetSpacing()[d]; + } +} + +template +typename GaussianInterpolateImageFunction +::OutputType +GaussianInterpolateImageFunction +::EvaluateAtContinuousIndex( const ContinuousIndexType &cindex, + OutputType *grad ) const +{ + vnl_vector erfArray[ImageDimension]; + vnl_vector gerfArray[ImageDimension]; + + // Compute the ERF difference arrays + for( unsigned int d = 0; d < ImageDimension; d++ ) + { + bool evaluateGradient = false; + if( grad ) + { + evaluateGradient = true; + } + this->ComputeErrorFunctionArray( d, cindex[d], erfArray[d], + gerfArray[d], evaluateGradient ); + } + + RealType sum_me = 0.0; + RealType sum_m = 0.0; + ArrayType dsum_me; + ArrayType dsum_m; + ArrayType dw; + + dsum_m.Fill( 0.0 ); + dsum_me.Fill( 0.0 ); + + // Loop over the voxels in the region identified + ImageRegion region; + for( unsigned int d = 0; d < ImageDimension; d++ ) + { + int boundingBoxSize = static_cast( + this->m_BoundingBoxEnd[d] - this->m_BoundingBoxStart[d] + 0.5 ); + int begin = vnl_math_max( 0, static_cast( vcl_floor( cindex[d] - + this->m_BoundingBoxStart[d] - this->m_CutoffDistance[d] ) ) ); + int end = vnl_math_min( boundingBoxSize, static_cast( vcl_ceil( + cindex[d] - this->m_BoundingBoxStart[d] + this->m_CutoffDistance[d] ) ) ); + region.SetIndex( d, begin ); + region.SetSize( d, end - begin ); + } + + ImageRegionConstIteratorWithIndex It( + this->GetInputImage(), region ); + for( It.GoToBegin(); !It.IsAtEnd(); ++It ) + { + unsigned int j = It.GetIndex()[0]; + RealType w = erfArray[0][j]; + if( grad ) + { + dw[0] = gerfArray[0][j]; + for( unsigned int d = 1; d < ImageDimension; d++ ) + { + dw[d] = erfArray[0][j]; + } + } + for( unsigned int d = 1; d < ImageDimension; d++) + { + j = It.GetIndex()[d]; + w *= erfArray[d][j]; + if( grad ) + { + for( unsigned int q = 0; q < ImageDimension; q++ ) + { + if( d == q ) + { + dw[q] *= gerfArray[d][j]; + } + else + { + dw[q] *= erfArray[d][j]; + } + } + } + } + RealType V = It.Get(); + sum_me += V * w; + sum_m += w; + if( grad ) + { + for( unsigned int q = 0; q < ImageDimension; q++ ) + { + dsum_me[q] += V * dw[q]; + dsum_m[q] += dw[q]; + } + } + } + RealType rc = sum_me / sum_m; + + if( grad ) + { + for( unsigned int q = 0; q < ImageDimension; q++ ) + { + grad[q] = ( dsum_me[q] - rc * dsum_m[q] ) / sum_m; + grad[q] /= -vnl_math::sqrt2 * this->m_Sigma[q]; + } + } + + return rc; +} + +template +void +GaussianInterpolateImageFunction +::ComputeErrorFunctionArray( unsigned int dimension, RealType cindex, + vnl_vector &erfArray, vnl_vector &gerfArray, + bool evaluateGradient ) const +{ + // Determine the range of voxels along the line where to evaluate erf + int boundingBoxSize = static_cast( + this->m_BoundingBoxEnd[dimension] - this->m_BoundingBoxStart[dimension] + + 0.5 ); + int begin = vnl_math_max( 0, static_cast( vcl_floor( cindex - + this->m_BoundingBoxStart[dimension] - + this->m_CutoffDistance[dimension] ) ) ); + int end = vnl_math_min( boundingBoxSize, static_cast( vcl_ceil( cindex - + this->m_BoundingBoxStart[dimension] + + this->m_CutoffDistance[dimension] ) ) ); + + erfArray.set_size( boundingBoxSize ); + gerfArray.set_size( boundingBoxSize ); + + // Start at the first voxel + RealType t = ( this->m_BoundingBoxStart[dimension] - cindex + + static_cast( begin ) ) * this->m_ScalingFactor[dimension]; + RealType e_last = vnl_erf( t ); + RealType g_last = 0.0; + if( evaluateGradient ) + { + g_last = vnl_math::two_over_sqrtpi * vcl_exp( -vnl_math_sqr( t ) ); + } + + for( int i = begin; i < end; i++ ) + { + t += this->m_ScalingFactor[dimension]; + RealType e_now = vnl_erf( t ); + erfArray[i] = e_now - e_last; + if( evaluateGradient ) + { + RealType g_now = vnl_math::two_over_sqrtpi * vcl_exp( -vnl_math_sqr( t ) ); + gerfArray[i] = g_now - g_last; + g_last = g_now; + } + e_last = e_now; + } +} + +} // namespace itk + +#endif diff --git a/Modules/DiffusionImaging/Algorithms/itkTractsToDWIImageFilter.cpp b/Modules/DiffusionImaging/Algorithms/itkTractsToDWIImageFilter.cpp new file mode 100644 index 0000000000..b909578354 --- /dev/null +++ b/Modules/DiffusionImaging/Algorithms/itkTractsToDWIImageFilter.cpp @@ -0,0 +1,335 @@ +/*========================================================================= +Copyright (c) German Cancer Research Center, Division of Medical and +Biological Informatics. All rights reserved. +See MITKCopyright.txt or http://www.mitk.org/copyright.html for details. + +This software is distributed WITHOUT ANY WARRANTY; without even +the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR +PURPOSE. See the above copyright notices for more information. + +=========================================================================*/ +#include "itkTractsToDWIImageFilter.h" +#include "math.h" +#include "vtkActor.h" +#include + +struct LessDereference { + template + bool operator()(const T * lhs, const T * rhs) const { + return *lhs < *rhs; + } +}; + +namespace itk +{ + TractsToDWIImageFilter::TractsToDWIImageFilter() + { + this->SetNumberOfRequiredInputs(1); + } + + TractsToDWIImageFilter::~TractsToDWIImageFilter() + { + } + + void TractsToDWIImageFilter::GenerateParticleGrid() + { + MITK_INFO << "Generating particle grid from fiber bundle"; + m_ParticleGrid = mitk::ParticleGrid::New(); + + float* bounds = m_FiberBundle->GetBounds(); + int size[] = {(int)bounds[0],(int)bounds[1],(int)bounds[2]}; + m_ParticleGrid->SetSize(size); + + typedef itk::Point ContainerPointType; //no need to init, is no smartpointer + typedef itk::VectorContainer ContainerTractType; + typedef itk::VectorContainer< unsigned int, ContainerTractType::Pointer > ContainerType; //init via smartpointer + + ContainerType::Pointer tractContainer = m_FiberBundle->GetTractContainer(); + + for (int i=0; iSize(); i++) + { + ContainerTractType::Pointer tract = tractContainer->GetElement(i); + for (int j=0; jSize(); j++) + { + vnl_vector_fixed pos = tract->GetElement(j).GetVnlVector(); + vnl_vector_fixed next; + vnl_vector_fixed dir; + + if (jSize()-1) + next = tract->GetElement(j+1).GetVnlVector(); + else + next = tract->GetElement(j-1).GetVnlVector(); + + dir = next-pos; + dir.normalize(); + mitk::Particle* p = new mitk::Particle(); + p->SetPosition(pos); + p->SetDirection(dir); + m_ParticleGrid->AddParticle(p); + } + } + } + + void TractsToDWIImageFilter::GenerateData() + { + + DWIImageType::Pointer originalDwiImage = dynamic_cast(this->ProcessObject::GetInput(0)); + if (originalDwiImage.IsNull()) + { + MITK_INFO << "no dwi image to extract b0 specified"; + return; + } + short* originalDwiImageBuffer = (short*) originalDwiImage->GetBufferPointer(); + + if (this->m_FiberBundle.IsNotNull()) + this->GenerateParticleGrid(); + else + { + MITK_INFO << "no fiber bundle specified"; + return; + } + + MITK_INFO << "reconstructing dwi-image from particle grid"; + + float* bounds = m_FiberBundle->GetBounds(); + ImageRegion<3> region; + region.SetSize(0, bounds[0]); + region.SetSize(1, bounds[1]); + region.SetSize(2, bounds[2]); + + m_Size[0] = bounds[0]; + m_Size[1] = bounds[1]; + m_Size[2] = bounds[2]; + + mitk::Geometry3D::Pointer geom = this->m_FiberBundle->GetGeometry(); + + int numDirections = m_GradientDirections->Size(); + int notNullDirections = 0; + for (int i=0; iGetElement(i); + if (dir[0]!=0 || dir[1]!=0 || dir[2]!=0) + notNullDirections++; + } + MITK_INFO << "Gradient directions: " << notNullDirections; + MITK_INFO << "B0 images: " << numDirections-notNullDirections; + +// short *averageImageBuffer = new short[m_Size[0]*m_Size[1]*m_Size[2]]; +// for (int x=0; xGetElement(i); +// if (dir[0]!=0 || dir[1]!=0 || dir[2]!=0) +// val += originalDwiImageBuffer[index]; +// } +// averageImageBuffer[x + m_Size[0]*(y+m_Size[1]*z)] = (short)val/notNullDirections; +// } + + // allocate output image + m_OutImage = static_cast< DWIImageType* >(this->ProcessObject::GetOutput(0)); + m_OutImage->SetSpacing( geom->GetSpacing() ); // Set the image spacing + m_OutImage->SetOrigin( geom->GetOrigin() ); // Set the image origin + itk::Matrix matrix; + for (int i=0; i<3; i++) + for (int j=0; j<3; j++) + matrix[j][i] = geom->GetMatrixColumn(i)[j]; + m_OutImage->SetDirection( matrix ); // Set the image direction + m_OutImage->SetRegions( region ); // Set image region + m_OutImage->SetVectorLength( numDirections ); + m_OutImage->Allocate(); + + short* imageBuffer = (short*) m_OutImage->GetBufferPointer(); + + m_BesselApproxCoeff = this->ComputeFiberCorrelation(); + + int dist = 1; + if (this->m_ParticleWidth==0) + { + double voxsize[3]; + voxsize[0] = geom->GetSpacing().GetElement(0); + voxsize[1] = geom->GetSpacing().GetElement(1); + voxsize[2] = geom->GetSpacing().GetElement(2); + float vox_half; + if(voxsize[0]m_ParticleWidth*this->m_ParticleWidth; + + int bufferSize = m_Size[0]*m_Size[1]*m_Size[2]*numDirections; + int counter = 0; + short maxVal = 0; + + for (int x=0; x pos; + pos[0] = (float)x; + pos[1] = (float)y; + pos[2] = (float)z; + for(unsigned int i=0; i=notNullDirections){ + imageBuffer[index] = originalDwiImageBuffer[index]; + continue; + } + + vnl_vector_fixed temp = m_GradientDirections->GetElement(i); + + // gradient direction + vnl_vector_fixed n; + n[0] = temp[0]; + n[1] = temp[1]; + n[2] = temp[2]; + + for (int nx=-dist; nx<=dist; nx++) + if((x+nx)>=0 && (x+nx)=0 && (y+ny)=0 && (z+nz)GetParticleContainer(x+nx, y+ny, z+nz); + if(container.IsNotNull()) + for (int j=0; jSize(); j++) + { + mitk::Particle* particle = container->GetElement(j); + // Particle direction + vnl_vector_fixed ni = particle->GetDirection(); + // Particle position + vnl_vector_fixed xi = particle->GetPosition(); + + vnl_vector_fixed diff = pos-xi; + float angle = fabs(dot_product(n, ni)); + + val += std::exp(-diff.squared_magnitude()/width)*(std::exp(-m_bD*angle*angle)); + } + } + //MITK_INFO << val; averageImageBuffer[x + m_Size[0]*(y+m_Size[1]*z)] + imageBuffer[index] = (short)(val*100+1); + if (counter%((int)(bufferSize/10))==0) + MITK_INFO << 100*counter/bufferSize << "%"; + counter++; + if (imageBuffer[index] > maxVal) + maxVal = imageBuffer[index]; + } + } + //delete(averageImageBuffer); + for (int x=0; x=notNullDirections){ + imageBuffer[index] += maxVal; + continue; + } + } + } + } + + float TractsToDWIImageFilter::Bessel(float x) + { + float y = x*x; + float z = y*y; + float erg = m_BesselApproxCoeff[0]; + erg += y*m_BesselApproxCoeff[1]; + erg += z*m_BesselApproxCoeff[2]; + erg += z*y*m_BesselApproxCoeff[3]; + return erg; + } + + float* TractsToDWIImageFilter::ComputeFiberCorrelation(){ + + float bD = m_bD; + + vnl_matrix_fixed bDir = + *itk::PointShell >::DistributePointShell(); + + const int N = QBALL_ODFSIZE; + + vnl_matrix_fixed C = bDir.transpose()*bDir; + vnl_matrix_fixed Q = C; + for(int i=0; i P = Q*Q; + + std::vector pointer; + pointer.reserve(N*N); + double * start = C.data_block(); + double * end = start + N*N; + for (double * iter = start; iter != end; ++iter) + { + pointer.push_back(iter); + } + std::sort(pointer.begin(), pointer.end(), LessDereference()); + + vnl_vector_fixed alpha; + vnl_vector_fixed beta; + for (int i=0; im_Meanval_sq = (sum*sum)/N; + + vnl_vector_fixed alpha_0; + vnl_vector_fixed alpha_2; + vnl_vector_fixed alpha_4; + vnl_vector_fixed alpha_6; + for(int i=0; i T; + T.set_column(0,alpha_0); + T.set_column(1,alpha_2); + T.set_column(2,alpha_4); + T.set_column(3,alpha_6); + + vnl_vector_fixed coeff = vnl_matrix_inverse(T).pinverse()*beta; + + float* retval = new float[4]; + retval[0] = coeff(0); + retval[1] = coeff(1); + retval[2] = coeff(2); + retval[3] = coeff(3); + return retval; + } +} diff --git a/Modules/DiffusionImaging/Algorithms/itkTractsToDWIImageFilter.h b/Modules/DiffusionImaging/Algorithms/itkTractsToDWIImageFilter.h new file mode 100644 index 0000000000..c39faed42a --- /dev/null +++ b/Modules/DiffusionImaging/Algorithms/itkTractsToDWIImageFilter.h @@ -0,0 +1,105 @@ +/*========================================================================= + +Program: Medical Imaging & Interaction Toolkit +Module: $RCSfile$ +Language: C++ +Date: $Date$ +Version: $Revision$ + +Copyright (c) German Cancer Research Center, Division of Medical and +Biological Informatics. All rights reserved. +See MITKCopyright.txt or http://www.mitk.org/copyright.html for details. + +This software is distributed WITHOUT ANY WARRANTY; without even +the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR +PURPOSE. See the above copyright notices for more information. + +=========================================================================*/ +#ifndef __itkTractsToDWIImageFilter_h__ +#define __itkTractsToDWIImageFilter_h__ + +// MITK +#include +#include +#include + +// ITK +#include +#include +#include +#include + +typedef itk::VectorImage< short, 3 > DWIImageType; + +namespace itk +{ + class TractsToDWIImageFilter : + public ImageToImageFilter< DWIImageType, DWIImageType >{ + + public: + typedef TractsToDWIImageFilter Self; + typedef ImageToImageFilter< DWIImageType, DWIImageType > Superclass; + typedef SmartPointer< Self > Pointer; + typedef SmartPointer< const Self > ConstPointer; + + typedef itk::VectorContainer< int, mitk::Particle* > ParticleContainerType; + typedef mitk::ParticleGrid::Pointer ParticleGridType; + + typedef mitk::FiberBundle::Pointer FiberBundleType; + + typedef vnl_vector_fixed< double, 3 > GradientDirectionType; + typedef itk::VectorContainer< unsigned int, GradientDirectionType > GradientDirectionContainerType; + + itkNewMacro(Self); + itkTypeMacro( TractsToDWIImageFilter, ImageToImageFilter ); + + // odf type + typedef OrientationDistributionFunction OdfType; + + itkSetMacro(FiberBundle, FiberBundleType); + itkGetMacro(FiberBundle, FiberBundleType); + + itkSetMacro(ParticleGrid, ParticleGridType); + itkGetMacro(ParticleGrid, ParticleGridType); + + itkSetMacro(bD, int); + itkGetMacro(bD, int); + + itkSetMacro(ParticleWidth, float); + itkGetMacro(ParticleWidth, float); + + itkSetMacro(ParticleLength, float); + itkGetMacro(ParticleLength, float); + + itkSetMacro(ParticleWeight, float); + itkGetMacro(ParticleWeight, float); + + itkSetMacro(GradientDirections, GradientDirectionContainerType::Pointer); + itkGetMacro(GradientDirections, GradientDirectionContainerType::Pointer); + + void GenerateData(); + + protected: + + TractsToDWIImageFilter(); + virtual ~TractsToDWIImageFilter(); + + GradientDirectionContainerType::Pointer m_GradientDirections; + FiberBundleType m_FiberBundle; + ParticleGridType m_ParticleGrid; + DWIImageType::Pointer m_OutImage; + int m_Size[3]; + float m_Meanval_sq; + + float *m_BesselApproxCoeff; + float* ComputeFiberCorrelation(); + float Bessel(float x); + void GenerateParticleGrid(); + + int m_bD; + float m_ParticleWidth; + float m_ParticleLength; + float m_ParticleWeight; + }; +} +#endif diff --git a/Modules/DiffusionImaging/Algorithms/itkTractsToFiberEndingsImageFilter.cpp b/Modules/DiffusionImaging/Algorithms/itkTractsToFiberEndingsImageFilter.cpp new file mode 100644 index 0000000000..a7148a7ef7 --- /dev/null +++ b/Modules/DiffusionImaging/Algorithms/itkTractsToFiberEndingsImageFilter.cpp @@ -0,0 +1,253 @@ +#include "itkTractsToFiberEndingsImageFilter.h" + +#include "itkBSplineUpsampleImageFilter.h" + +#define __CEIL_UCHAR__(val) (val) = \ +( (val) < 0 ) ? ( 0 ) : ( ( (val)>255 ) ? ( 255 ) : ( (val) ) ); + +namespace itk{ + + template< class TInputImage, class TOutputPixelType > + TractsToFiberEndingsImageFilter< TInputImage, TOutputPixelType > + ::TractsToFiberEndingsImageFilter() + { + this->SetNumberOfRequiredInputs(0); + } + + template< class TInputImage, class TOutputPixelType > + TractsToFiberEndingsImageFilter< TInputImage, TOutputPixelType > + ::~TractsToFiberEndingsImageFilter() + { + } + + template< class TInputImage, class TOutputPixelType > + void TractsToFiberEndingsImageFilter< TInputImage, TOutputPixelType > + ::GenerateData() + { + MITK_INFO << "Generating 2D fiber endings image"; + bool isRgba = false; + if(&typeid(TOutputPixelType) == &typeid(itk::RGBAPixel)) + { + isRgba = true; + } + else if(&typeid(TOutputPixelType) != &typeid(unsigned char)) + { + MITK_INFO << "Only 'unsigned char' and 'itk::RGBAPixel supported as OutputPixelType"; + return; + } + + mitk::Geometry3D::Pointer geometry = m_FiberBundle->GetGeometry(); + + typename OutputImageType::Pointer outImage = + static_cast< OutputImageType * >(this->ProcessObject::GetOutput(0)); + + outImage->SetSpacing( geometry->GetSpacing()/m_UpsamplingFactor ); // Set the image spacing + + mitk::Point3D origin = geometry->GetOrigin(); + mitk::Point3D indexOrigin; + geometry->WorldToIndex(origin, indexOrigin); + indexOrigin[0] = indexOrigin[0] - .5 * (1.0-1.0/m_UpsamplingFactor); + indexOrigin[1] = indexOrigin[1] - .5 * (1.0-1.0/m_UpsamplingFactor); + indexOrigin[2] = indexOrigin[2] - .5 * (1.0-1.0/m_UpsamplingFactor); + mitk::Point3D newOrigin; + geometry->IndexToWorld(indexOrigin, newOrigin); + + outImage->SetOrigin( newOrigin ); // Set the image origin + itk::Matrix matrix; + for (int i=0; i<3; i++) + for (int j=0; j<3; j++) + matrix[j][i] = geometry->GetMatrixColumn(i)[j]; + outImage->SetDirection( matrix ); // Set the image direction + + float* bounds = m_FiberBundle->GetBounds(); + ImageRegion<3> upsampledRegion; + upsampledRegion.SetSize(0, bounds[0]); + upsampledRegion.SetSize(1, bounds[1]); + upsampledRegion.SetSize(2, bounds[2]); + + typename InputImageType::RegionType::SizeType upsampledSize = upsampledRegion.GetSize(); + for (unsigned int n = 0; n < 3; n++) + { + upsampledSize[n] = upsampledSize[n] * m_UpsamplingFactor; + } + upsampledRegion.SetSize( upsampledSize ); + outImage->SetRegions( upsampledRegion ); + + outImage->Allocate(); + // itk::RGBAPixel pix; + // pix.Set(0,0,0,0); + // outImage->FillBuffer(pix); + + int w = upsampledSize[0]; + int h = upsampledSize[1]; + int d = upsampledSize[2]; + + + unsigned char* accuout; + float* accu; + + accuout = reinterpret_cast(outImage->GetBufferPointer()); + + if(isRgba) + { + accu = new float[w*h*d*4]; + for (int i=0; iGetNumTracts(); + for( int i=0; i > vertices; + + // get fiber start point + int numVertices = m_FiberBundle->GetNumPoints(i); + itk::Point start = m_FiberBundle->GetPoint(i,0); + vertices.push_back(start); + + // get fiber end point + if(numVertices>1) + { + itk::Point end = m_FiberBundle->GetPoint(i,numVertices-1); + vertices.push_back(end); + } + + //////////////////// + // fill output image + + // for each vertex + for( int j=0; j vertex = vertices.at(j); + + // scaling coordinates (index coords scale with upsampling) + vertex[0] = (vertex[0]+0.5) * m_UpsamplingFactor; + vertex[1] = (vertex[1]+0.5) * m_UpsamplingFactor; + vertex[2] = (vertex[2]+0.5) * m_UpsamplingFactor; + + // int coordinates inside image? + int px = (int) (vertex[0]); + if (px < 0 || px >= w) + continue; + int py = (int) (vertex[1]); + if (py < 0 || py >= h) + continue; + int pz = (int) (vertex[2]); + if (pz < 0 || pz >= d) + continue; + + // float fraction of coordinates + float frac_x = vertex[0] - px; + float frac_y = vertex[1] - py; + float frac_z = vertex[2] - pz; + + float scale = 100 * pow((float)m_UpsamplingFactor,3); + + if(isRgba) + { + // add to r-channel in output image + accu[0+4*( px + w*(py + h*pz ))] += (1-frac_x)*(1-frac_y)*(1-frac_z) * scale; + accu[0+4*( px + w*(py+1+ h*pz ))] += (1-frac_x)*( frac_y)*(1-frac_z) * scale; + accu[0+4*( px + w*(py + h*pz+h))] += (1-frac_x)*(1-frac_y)*( frac_z) * scale; + accu[0+4*( px + w*(py+1+ h*pz+h))] += (1-frac_x)*( frac_y)*( frac_z) * scale; + accu[0+4*( px+1 + w*(py + h*pz ))] += ( frac_x)*(1-frac_y)*(1-frac_z) * scale; + accu[0+4*( px+1 + w*(py + h*pz+h))] += ( frac_x)*(1-frac_y)*( frac_z) * scale; + accu[0+4*( px+1 + w*(py+1+ h*pz ))] += ( frac_x)*( frac_y)*(1-frac_z) * scale; + accu[0+4*( px+1 + w*(py+1+ h*pz+h))] += ( frac_x)*( frac_y)*( frac_z) * scale; + + // add to a-channel in output image + accu[3+4*( px + w*(py + h*pz ))] = 1; + accu[3+4*( px + w*(py+1+ h*pz ))] = 1; + accu[3+4*( px + w*(py + h*pz+h))] = 1; + accu[3+4*( px + w*(py+1+ h*pz+h))] = 1; + accu[3+4*( px+1 + w*(py + h*pz ))] = 1; + accu[3+4*( px+1 + w*(py + h*pz+h))] = 1; + accu[3+4*( px+1 + w*(py+1+ h*pz ))] = 1; + accu[3+4*( px+1 + w*(py+1+ h*pz+h))] = 1; + } + else + { + accu[( px + w*(py + h*pz ))] += (1-frac_x)*(1-frac_y)*(1-frac_z) * scale; + accu[( px + w*(py+1+ h*pz ))] += (1-frac_x)*( frac_y)*(1-frac_z) * scale; + accu[( px + w*(py + h*pz+h))] += (1-frac_x)*(1-frac_y)*( frac_z) * scale; + accu[( px + w*(py+1+ h*pz+h))] += (1-frac_x)*( frac_y)*( frac_z) * scale; + accu[( px+1 + w*(py + h*pz ))] += ( frac_x)*(1-frac_y)*(1-frac_z) * scale; + accu[( px+1 + w*(py + h*pz+h))] += ( frac_x)*(1-frac_y)*( frac_z) * scale; + accu[( px+1 + w*(py+1+ h*pz ))] += ( frac_x)*( frac_y)*(1-frac_z) * scale; + accu[( px+1 + w*(py+1+ h*pz+h))] += ( frac_x)*( frac_y)*( frac_z) * scale; + } + + } + } + + float maxRgb = 0.000000001; + float maxInt = 0.000000001; + int numPix; + + if(isRgba) + { + numPix = w*h*d*4; + + // calc maxima + for(int i=0; i maxRgb) + { + maxRgb = accu[i]; + } + } + else + { + if(accu[i] > maxInt) + { + maxInt = accu[i]; + } + } + } + + // write output, normalized uchar 0..255 + for(int i=0; i maxInt) + { + maxInt = accu[i]; + } + } + + // write output, normalized uchar 0..255 + for(int i=0; i +class TractsToFiberEndingsImageFilter : + public ImageToImageFilter > +//huhu public ImageToImageFilter > +{ + +public: + typedef TractsToFiberEndingsImageFilter Self; + typedef ImageToImageFilter > Superclass; +//huhu typedef ImageToImageFilter > Superclass; + typedef SmartPointer< Self > Pointer; + typedef SmartPointer< const Self > ConstPointer; + + itkNewMacro(Self); + itkTypeMacro( TractsToFiberEndingsImageFilter, + ImageToImageFilter ); + + /** Types for the Output Image**/ + typedef TInputImage InputImageType; + + /** Types for the Output Image**/ + typedef itk::Image OutputImageType; + + /** Type for the directions **/ + typedef VectorContainer< unsigned int, vnl_vector_fixed< double, 3 > > + TractOrientationsType; + + /** Type for the directions container **/ + typedef VectorContainer< unsigned int, TractOrientationsType > + TractOrientationsContainerType; + + /** Upsampling factor **/ + itkSetMacro( UpsamplingFactor, unsigned int); + itkGetMacro( UpsamplingFactor, unsigned int); + + itkSetMacro(FiberBundle, mitk::FiberBundle::Pointer); + itkGetMacro(FiberBundle, mitk::FiberBundle::Pointer); + + void GenerateData(); + +protected: + + TractsToFiberEndingsImageFilter(); + virtual ~TractsToFiberEndingsImageFilter(); + + unsigned int m_UpsamplingFactor; + mitk::FiberBundle::Pointer m_FiberBundle; + +}; + +} + +#ifndef ITK_MANUAL_INSTANTIATION +#include "itkTractsToFiberEndingsImageFilter.cpp" +#endif + +#endif // __itkTractsToFiberEndingsImageFilter_h__ diff --git a/Modules/DiffusionImaging/IODataStructures/FiberBundle/mitkParticleGrid.cpp b/Modules/DiffusionImaging/IODataStructures/FiberBundle/mitkParticleGrid.cpp new file mode 100644 index 0000000000..5938650b65 --- /dev/null +++ b/Modules/DiffusionImaging/IODataStructures/FiberBundle/mitkParticleGrid.cpp @@ -0,0 +1,131 @@ +/*========================================================================= + Copyright (c) German Cancer Research Center, Division of Medical and + Biological Informatics. All rights reserved. + See MITKCopyright.txt or http://www.mitk.org/copyright.html for details. + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notices for more information. + + =========================================================================*/ + +#include "mitkParticleGrid.h" + +namespace mitk +{ + ParticleGrid::ParticleGrid(): + m_ParticleWeight(0.009) + { + m_ParticleGrid = ParticleGridType::New(); + m_ParticleList = ParticleContainerType::New(); + } + + ParticleGrid::~ParticleGrid() + { + for (int i=0; iSize(); i++) + { + delete(m_ParticleList->GetElement(i)); + } + } + + int ParticleGrid::GetNumParticles() + { + return this->m_ParticleList->Size(); + } + + void ParticleGrid::SetSize(int size[3]) + { + m_Size[0] = size[0]; + m_Size[1] = size[1]; + m_Size[2] = size[2]; + } + + void ParticleGrid::AddParticle(mitk::Particle* particle) + { + vnl_vector_fixed pos = particle->GetPosition(); + int x = pos[0]+0.5; + int y = pos[1]+0.5; + int z = pos[2]+0.5; + + if (!(x>=0 && x=0 && y=0 && zIndexExists(index)) + { + container = m_ParticleGrid->GetElement(index); + if(container.IsNull()) + { + container = ParticleContainerType::New(); + container->InsertElement(container->Size(), particle); + this->m_ParticleGrid->InsertElement(index, container); + } + else + { + container->InsertElement(container->Size(), particle); + } + } + else + { + container = ParticleContainerType::New(); + container->InsertElement(container->Size(), particle); + this->m_ParticleGrid->InsertElement(index, container); + } + + this->m_ParticleList->InsertElement(m_ParticleList->Size(),particle); + } + + Particle* ParticleGrid::GetParticle(int x, int y, int z, int pos) + { + int index = x + m_Size[0]*(y+m_Size[1]*z); + + if (!m_ParticleGrid->IndexExists(index)) + return NULL; + + ParticleContainerType::Pointer container = m_ParticleGrid->GetElement(index); + + if (!container->IndexExists(pos)) + return NULL; + + return container->GetElement(pos); + } + + ParticleGrid::ParticleContainerType::Pointer ParticleGrid::GetParticleContainer(int x, int y, int z) + { + int index = x + m_Size[0]*(y+m_Size[1]*z); + + return GetParticleContainer(index); + } + + ParticleGrid::ParticleContainerType::Pointer ParticleGrid::GetParticleContainer(int index) + { + if (!m_ParticleGrid->IndexExists(index)) + return NULL; + + return m_ParticleGrid->GetElement(index); + } + + /* NECESSARY IMPLEMENTATION OF SUPERCLASS METHODS */ + void ParticleGrid::UpdateOutputInformation() + { + + } + void ParticleGrid::SetRequestedRegionToLargestPossibleRegion() + { + + } + bool ParticleGrid::RequestedRegionIsOutsideOfTheBufferedRegion() + { + return false; + } + bool ParticleGrid::VerifyRequestedRegion() + { + return false; + } + void ParticleGrid::SetRequestedRegion( itk::DataObject *data ) + { + + } +} diff --git a/Modules/DiffusionImaging/IODataStructures/FiberBundle/mitkParticleGrid.h b/Modules/DiffusionImaging/IODataStructures/FiberBundle/mitkParticleGrid.h new file mode 100644 index 0000000000..f1db7d26eb --- /dev/null +++ b/Modules/DiffusionImaging/IODataStructures/FiberBundle/mitkParticleGrid.h @@ -0,0 +1,75 @@ + +/*========================================================================= + Program: Medical Imaging & Interaction Toolkit + Module: $RCSfile$ + Language: C++ + Date: $Date$ + Version: $Revision: 11989 $ + + Copyright (c) German Cancer Research Center, Division of Medical and + Biological Informatics. All rights reserved. + See MITKCopyright.txt or http://www.mitk.org/copyright.html for details. + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notices for more information. + + =========================================================================*/ + + +#ifndef _MITK_ParticleGrid_H +#define _MITK_ParticleGrid_H + +#include "mitkBaseData.h" +#include "mitkParticle.h" +#include "MitkDiffusionImagingExports.h" + +namespace mitk { + + /** + * \brief 3D grid conatiner for fiber Particles; */ + class MitkDiffusionImaging_EXPORT ParticleGrid : public BaseData + { + public: + + mitkClassMacro( ParticleGrid, BaseData ); + itkNewMacro( Self ); + + // grid types + typedef itk::VectorContainer< int, mitk::Particle* > ParticleContainerType; + typedef itk::VectorContainer< int, ParticleContainerType::Pointer > ParticleGridType; + + void AddParticle(mitk::Particle* particle); + mitk::Particle* GetParticle(int x, int y, int z, int pos); + ParticleContainerType::Pointer GetParticleContainer(int x, int y, int z); + ParticleContainerType::Pointer GetParticleContainer(int index); + + itkSetMacro(ParticleWeight, float); + itkGetMacro(ParticleWeight, float); + + itkGetMacro(ParticleList, ParticleContainerType::Pointer); + + void SetSize(int size[3]); + int GetNumParticles(); + + // virtual methods that need to be implemented + virtual void UpdateOutputInformation(); + virtual void SetRequestedRegionToLargestPossibleRegion(); + virtual bool RequestedRegionIsOutsideOfTheBufferedRegion(); + virtual bool VerifyRequestedRegion(); + virtual void SetRequestedRegion( itk::DataObject *data ); + + protected: + ParticleGrid(); + virtual ~ParticleGrid(); + + ParticleGridType::Pointer m_ParticleGrid; + ParticleContainerType::Pointer m_ParticleList; + int m_Size[3]; + + float m_ParticleWeight; + }; + +} // namespace mitk + +#endif /* _MITK_ParticleGrid_H */ diff --git a/Modules/DiffusionImaging/files.cmake b/Modules/DiffusionImaging/files.cmake index 59c190f5cc..20befd0ca1 100644 --- a/Modules/DiffusionImaging/files.cmake +++ b/Modules/DiffusionImaging/files.cmake @@ -1,104 +1,111 @@ SET(CPP_FILES # DicomImport DicomImport/mitkDicomDiffusionImageReader.cpp DicomImport/mitkGroupDiffusionHeadersFilter.cpp DicomImport/mitkDicomDiffusionImageHeaderReader.cpp DicomImport/mitkGEDicomDiffusionImageHeaderReader.cpp DicomImport/mitkPhilipsDicomDiffusionImageHeaderReader.cpp DicomImport/mitkSiemensDicomDiffusionImageHeaderReader.cpp DicomImport/mitkSiemensMosaicDicomDiffusionImageHeaderReader.cpp # DataStructures IODataStructures/mitkDiffusionImagingObjectFactory.cpp # DataStructures -> DWI IODataStructures/DiffusionWeightedImages/mitkDiffusionImageHeaderInformation.cpp IODataStructures/DiffusionWeightedImages/mitkDiffusionImageSource.cpp IODataStructures/DiffusionWeightedImages/mitkNrrdDiffusionImageReader.cpp IODataStructures/DiffusionWeightedImages/mitkNrrdDiffusionImageWriter.cpp IODataStructures/DiffusionWeightedImages/mitkNrrdDiffusionImageIOFactory.cpp IODataStructures/DiffusionWeightedImages/mitkNrrdDiffusionImageWriterFactory.cpp IODataStructures/DiffusionWeightedImages/mitkDiffusionImageSerializer.cpp # DataStructures -> QBall IODataStructures/QBallImages/mitkQBallImageSource.cpp IODataStructures/QBallImages/mitkNrrdQBallImageReader.cpp IODataStructures/QBallImages/mitkNrrdQBallImageWriter.cpp IODataStructures/QBallImages/mitkNrrdQBallImageIOFactory.cpp IODataStructures/QBallImages/mitkNrrdQBallImageWriterFactory.cpp IODataStructures/QBallImages/mitkQBallImage.cpp IODataStructures/QBallImages/mitkQBallImageSerializer.cpp # DataStructures -> Tensor IODataStructures/TensorImages/mitkTensorImageSource.cpp IODataStructures/TensorImages/mitkNrrdTensorImageReader.cpp IODataStructures/TensorImages/mitkNrrdTensorImageWriter.cpp IODataStructures/TensorImages/mitkNrrdTensorImageIOFactory.cpp IODataStructures/TensorImages/mitkNrrdTensorImageWriterFactory.cpp IODataStructures/TensorImages/mitkTensorImage.cpp IODataStructures/TensorImages/mitkTensorImageSerializer.cpp # DataStructures -> FiberBundle IODataStructures/FiberBundle/mitkFiberBundle.cpp IODataStructures/FiberBundle/mitkFiberBundleWriter.cpp IODataStructures/FiberBundle/mitkFiberBundleReader.cpp IODataStructures/FiberBundle/mitkFiberBundleIOFactory.cpp IODataStructures/FiberBundle/mitkFiberBundleWriterFactory.cpp IODataStructures/FiberBundle/mitkFiberBundleSerializer.cpp IODataStructures/FiberBundle/mitkParticle.cpp + IODataStructures/FiberBundle/mitkParticleGrid.cpp # DataStructures -> PlanarFigureComposite IODataStructures/PlanarFigureComposite/mitkPlanarFigureComposite.cpp # Rendering Rendering/vtkMaskedProgrammableGlyphFilter.cpp Rendering/mitkCompositeMapper.cpp Rendering/mitkVectorImageVtkGlyphMapper3D.cpp Rendering/vtkOdfSource.cxx Rendering/vtkThickPlane.cxx Rendering/mitkOdfNormalizationMethodProperty.cpp Rendering/mitkOdfScaleByProperty.cpp Rendering/mitkFiberBundleMapper3D.cpp # Interactions Interactions/mitkFiberBundleInteractor.cpp + + # Algorithms + Algorithms/itkTractsToDWIImageFilter.cpp ) SET(H_FILES # Rendering Rendering/mitkDiffusionImageMapper.h Rendering/mitkOdfVtkMapper2D.h # Reconstruction Reconstruction/itkDiffusionQballReconstructionImageFilter.h Reconstruction/mitkTeemDiffusionTensor3DReconstructionImageFilter.h Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.h Reconstruction/itkPointShell.h Reconstruction/itkOrientationDistributionFunction.h # IO Datastructures IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.h IODataStructures/FiberBundle/itkSlowPolyLineParametricPath.h # Tractography Tractography/itkGlobalTractographyFilter.h # Algorithms Algorithms/itkDiffusionQballGeneralizedFaImageFilter.h Algorithms/itkDiffusionQballPrepareVisualizationImageFilter.h Algorithms/itkTensorDerivedMeasurementsFilter.h Algorithms/itkBrainMaskExtractionImageFilter.h Algorithms/itkB0ImageExtractionImageFilter.h Algorithms/itkTensorImageToDiffusionImageFilter.h Algorithms/itkTensorToL2NormImageFilter.h Algorithms/itkTractsToProbabilityImageFilter.h + Algorithms/itkTractsToDWIImageFilter.h + Algorithms/itkTractsToFiberEndingsImageFilter.h + Algorithms/itkGaussianInterpolateImageFunction.h ) SET( TOOL_FILES ) IF(WIN32) ENDIF(WIN32) #MITK_MULTIPLEX_PICTYPE( Algorithms/mitkImageRegistrationMethod-TYPE.cpp )