diff --git a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberProcessingView.cpp b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberProcessingView.cpp index 60524b3de1..ebb6128025 100644 --- a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberProcessingView.cpp +++ b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberProcessingView.cpp @@ -1,1738 +1,1740 @@ /*========================================================================= 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 "QmitkFiberProcessingView.h" #include // Qt #include // MITK #include #include #include #include #include #include #include #include #include #include #include // ITK #include #include #include #include #include #include #include #include const std::string QmitkFiberProcessingView::VIEW_ID = "org.mitk.views.fiberprocessing"; const std::string id_DataManager = "org.mitk.views.datamanager"; using namespace mitk; QmitkFiberProcessingView::QmitkFiberProcessingView() : QmitkFunctionality() , m_Controls( 0 ) , m_MultiWidget( NULL ) , m_EllipseCounter(0) , m_PolygonCounter(0) , m_UpsamplingFactor(5) { } // Destructor QmitkFiberProcessingView::~QmitkFiberProcessingView() { } void QmitkFiberProcessingView::CreateQtPartControl( QWidget *parent ) { // build up qt view, unless already done if ( !m_Controls ) { // create GUI widgets from the Qt Designer's .ui file m_Controls = new Ui::QmitkFiberProcessingViewControls; m_Controls->setupUi( parent ); m_Controls->doExtractFibersButton->setDisabled(true); m_Controls->PFCompoANDButton->setDisabled(true); m_Controls->PFCompoORButton->setDisabled(true); m_Controls->PFCompoNOTButton->setDisabled(true); m_Controls->m_PlanarFigureButtonsFrame->setEnabled(false); m_Controls->m_RectangleButton->setVisible(false); connect( m_Controls->doExtractFibersButton, SIGNAL(clicked()), this, SLOT(DoFiberExtraction()) ); 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(GenerateAndComposite()) ); connect(m_Controls->PFCompoORButton, SIGNAL(clicked()), this, SLOT(GenerateOrComposite()) ); connect(m_Controls->PFCompoNOTButton, SIGNAL(clicked()), this, SLOT(GenerateNotComposite()) ); 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_Extract3dButton, SIGNAL(clicked()), this, SLOT(Extract3d())); connect( m_Controls->m_ProcessFiberBundleButton, SIGNAL(clicked()), this, SLOT(ProcessSelectedBundles()) ); connect( m_Controls->m_ResampleFibersButton, SIGNAL(clicked()), this, SLOT(ResampleSelectedBundles()) ); } } void QmitkFiberProcessingView::Extract3d() { std::vector nodes = this->GetDataManagerSelection(); if (nodes.empty()) return; mitk::FiberBundleX::Pointer fib = mitk::FiberBundleX::New(); mitk::Surface::Pointer roi = mitk::Surface::New(); bool fibB = false; bool roiB = false; for (int i=0; i(nodes.at(i)->GetData())) { fib = dynamic_cast(nodes.at(i)->GetData()); fibB = true; } else if (dynamic_cast(nodes.at(i)->GetData())) { roi = dynamic_cast(nodes.at(i)->GetData()); roiB = true; } } if (!fibB) return; if (!roiB) return; vtkSmartPointer polyRoi = roi->GetVtkPolyData(); vtkSmartPointer polyFib = fib->GetFiberPolyData(); vtkSmartPointer selectEnclosedPoints = vtkSmartPointer::New(); selectEnclosedPoints->SetInput(polyFib); selectEnclosedPoints->SetSurface(polyRoi); selectEnclosedPoints->Update(); vtkSmartPointer newPoly = vtkSmartPointer::New(); vtkSmartPointer newCellArray = vtkSmartPointer::New(); vtkSmartPointer newPoints = vtkSmartPointer::New(); vtkSmartPointer newPolyComplement = vtkSmartPointer::New(); vtkSmartPointer newCellArrayComplement = vtkSmartPointer::New(); vtkSmartPointer newPointsComplement = vtkSmartPointer::New(); vtkSmartPointer vLines = polyFib->GetLines(); vLines->InitTraversal(); int numberOfLines = vLines->GetNumberOfCells(); // each line for (int j=0; jGetNextCell ( numPoints, points ); bool isPassing = false; // each point of this line for (int k=0; kIsInside(points[k])) { isPassing = true; // fill new polydata vtkSmartPointer container = vtkSmartPointer::New(); for (int k=0; kGetPoint(points[k]); vtkIdType pointId = newPoints->InsertNextPoint(point); container->GetPointIds()->InsertNextId(pointId); } newCellArray->InsertNextCell(container); break; } } if (!isPassing) { vtkSmartPointer container = vtkSmartPointer::New(); for (int k=0; kGetPoint(points[k]); vtkIdType pointId = newPointsComplement->InsertNextPoint(point); container->GetPointIds()->InsertNextId(pointId); } newCellArrayComplement->InsertNextCell(container); } } newPoly->SetPoints(newPoints); newPoly->SetLines(newCellArray); mitk::FiberBundleX::Pointer fb = mitk::FiberBundleX::New(newPoly); DataNode::Pointer newNode = DataNode::New(); newNode->SetData(fb); newNode->SetName("passing surface"); GetDefaultDataStorage()->Add(newNode); newPolyComplement->SetPoints(newPointsComplement); newPolyComplement->SetLines(newCellArrayComplement); mitk::FiberBundleX::Pointer fbComplement = mitk::FiberBundleX::New(newPolyComplement); DataNode::Pointer newNodeComplement = DataNode::New(); newNodeComplement->SetData(fbComplement); newNodeComplement->SetName("not passing surface"); GetDefaultDataStorage()->Add(newNodeComplement); } void QmitkFiberProcessingView::GenerateRoiImage(){ if (m_SelectedImage.IsNull() || m_SelectedPF.empty()) return; mitk::Image* image = const_cast(m_SelectedImage.GetPointer()); UCharImageType::Pointer temp = UCharImageType::New(); mitk::CastToItkImage(m_SelectedImage, temp); m_PlanarFigureImage = UCharImageType::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; iInitializeByItk(m_PlanarFigureImage.GetPointer()); tmpImage->SetVolume(m_PlanarFigureImage->GetBufferPointer()); node->SetData(tmpImage); node->SetName("ROI Image"); this->GetDefaultDataStorage()->Add(node); } void QmitkFiberProcessingView::CompositeExtraction(mitk::DataNode::Pointer node, mitk::Image* image) { if (dynamic_cast(node.GetPointer()->GetData()) && !dynamic_cast(node.GetPointer()->GetData())) { m_PlanarFigure = dynamic_cast(node.GetPointer()->GetData()); AccessFixedDimensionByItk_2( image, InternalReorientImagePlane, 3, m_PlanarFigure->GetGeometry(), -1); // itk::Image< unsigned char, 3 >::Pointer outimage = itk::Image< unsigned char, 3 >::New(); // outimage->SetSpacing( m_PlanarFigure->GetGeometry()->GetSpacing()/m_UpsamplingFactor ); // Set the image spacing // mitk::Point3D origin = m_PlanarFigure->GetGeometry()->GetOrigin(); // mitk::Point3D indexOrigin; // m_PlanarFigure->GetGeometry()->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; // m_PlanarFigure->GetGeometry()->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] = m_PlanarFigure->GetGeometry()->GetMatrixColumn(i)[j]/m_PlanarFigure->GetGeometry()->GetSpacing().GetElement(i); // outimage->SetDirection( matrix ); // Set the image direction // itk::ImageRegion<3> upsampledRegion; // upsampledRegion.SetSize(0, m_PlanarFigure->GetGeometry()->GetParametricExtentInMM(0)/m_PlanarFigure->GetGeometry()->GetSpacing()[0]); // upsampledRegion.SetSize(1, m_PlanarFigure->GetGeometry()->GetParametricExtentInMM(1)/m_PlanarFigure->GetGeometry()->GetSpacing()[1]); // upsampledRegion.SetSize(2, 1); // typename itk::Image< unsigned char, 3 >::RegionType::SizeType upsampledSize = upsampledRegion.GetSize(); // for (unsigned int n = 0; n < 2; n++) // { // upsampledSize[n] = upsampledSize[n] * m_UpsamplingFactor; // } // upsampledRegion.SetSize( upsampledSize ); // outimage->SetRegions( upsampledRegion ); // outimage->Allocate(); // this->m_InternalImage = mitk::Image::New(); // this->m_InternalImage->InitializeByItk( outimage.GetPointer() ); // this->m_InternalImage->SetVolume( outimage->GetBufferPointer() ); AccessFixedDimensionByItk_2( m_InternalImage, InternalCalculateMaskFromPlanarFigure, 3, 2, node->GetName() ); } } template < typename TPixel, unsigned int VImageDimension > void QmitkFiberProcessingView::InternalReorientImagePlane( const itk::Image< TPixel, VImageDimension > *image, 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 QmitkFiberProcessingView::InternalCalculateMaskFromPlanarFigure( itk::Image< TPixel, VImageDimension > *image, unsigned int axis, std::string nodeName ) { MITK_INFO << "InternalCalculateMaskFromPlanarFigure() start"; typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::CastImageFilter< ImageType, UCharImageType > CastFilterType; // Generate mask image as new image with same header as input image and // initialize with "1". UCharImageType::Pointer newMaskImage = UCharImageType::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; if (point3D[i0]<0) point3D[i0] = 0.0; else if (point3D[i0]>bounds[0]) point3D[i0] = bounds[0]-0.001; if (point3D[i1]<0) point3D[i1] = 0.0; else if (point3D[i1]>bounds[1]) point3D[i1] = bounds[1]-0.001; points->InsertNextPoint( point3D[i0], point3D[i1], -0.5 ); numberOfPoints++; } else { imageGeometry3D->WorldToIndex( point3D, point3D ); // Add point to polyline array points->InsertNextPoint( point3D[i0], point3D[i1], -0.5 ); numberOfPoints++; } } polyline->SetPoints( points ); points->Delete(); vtkIdType *ptIds = new vtkIdType[numberOfPoints]; for ( vtkIdType i = 0; i < numberOfPoints; ++i ) { ptIds[i] = i; } polyline->InsertNextCell( VTK_POLY_LINE, numberOfPoints, ptIds ); // Extrude the generated contour polygon vtkLinearExtrusionFilter *extrudeFilter = vtkLinearExtrusionFilter::New(); extrudeFilter->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< UCharImageType > ImageImportType; typedef itk::VTKImageExport< UCharImageType > 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 mask typedef itk::RegionOfInterestImageFilter< UCharImageType, UCharImageType > ROIMaskFilterType; typename ROIMaskFilterType::Pointer roi2 = ROIMaskFilterType::New(); roi2->SetRegionOfInterest(cropRegion); roi2->SetInput(m_InternalImageMask3D); roi2->Update(); m_InternalImageMask3D = roi2->GetOutput(); Image::Pointer tmpImage = Image::New(); tmpImage->InitializeByItk(m_InternalImageMask3D.GetPointer()); tmpImage->SetVolume(m_InternalImageMask3D->GetBufferPointer()); Image::Pointer tmpImage2 = Image::New(); tmpImage2->InitializeByItk(m_PlanarFigureImage.GetPointer()); const 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); point[i0] += 0.5; point[i1] += 0.5; 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 QmitkFiberProcessingView::StdMultiWidgetAvailable (QmitkStdMultiWidget &stdMultiWidget) { m_MultiWidget = &stdMultiWidget; } void QmitkFiberProcessingView::StdMultiWidgetNotAvailable() { m_MultiWidget = NULL; } /* OnSelectionChanged is registered to SelectionService, therefore no need to implement SelectionService Listener explicitly */ void QmitkFiberProcessingView::UpdateGui() { // are fiber bundles selected? if ( m_SelectedFB.empty() ) { m_Controls->m_JoinBundles->setEnabled(false); m_Controls->m_SubstractBundles->setEnabled(false); m_Controls->m_ProcessFiberBundleButton->setEnabled(false); m_Controls->doExtractFibersButton->setEnabled(false); m_Controls->m_Extract3dButton->setEnabled(false); m_Controls->m_ResampleFibersButton->setEnabled(false); } else { m_Controls->m_PlanarFigureButtonsFrame->setEnabled(true); m_Controls->m_ProcessFiberBundleButton->setEnabled(true); m_Controls->m_ResampleFibersButton->setEnabled(true); if (m_Surfaces.size()>0) m_Controls->m_Extract3dButton->setEnabled(true); // one bundle and one planar figure needed to extract fibers if (!m_SelectedPF.empty()) m_Controls->doExtractFibersButton->setEnabled(true); // more than two bundles needed to join/subtract if (m_SelectedFB.size() > 1) { 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); } } // are planar figures selected? if ( m_SelectedPF.empty() ) { m_Controls->doExtractFibersButton->setEnabled(false); m_Controls->PFCompoANDButton->setEnabled(false); m_Controls->PFCompoORButton->setEnabled(false); m_Controls->PFCompoNOTButton->setEnabled(false); m_Controls->m_GenerateRoiImage->setEnabled(false); } else { if ( m_SelectedImage.IsNotNull() ) m_Controls->m_GenerateRoiImage->setEnabled(true); else m_Controls->m_GenerateRoiImage->setEnabled(false); if (m_SelectedPF.size() > 1) { m_Controls->PFCompoANDButton->setEnabled(true); m_Controls->PFCompoORButton->setEnabled(true); m_Controls->PFCompoNOTButton->setEnabled(false); } else { m_Controls->PFCompoANDButton->setEnabled(false); m_Controls->PFCompoORButton->setEnabled(false); m_Controls->PFCompoNOTButton->setEnabled(true); } } } void QmitkFiberProcessingView::OnSelectionChanged( std::vector nodes ) { if ( !this->IsVisible() ) return; //reset existing Vectors containing FiberBundles and PlanarFigures from a previous selection m_SelectedFB.clear(); m_SelectedPF.clear(); m_Surfaces.clear(); m_SelectedImage = NULL; for( std::vector::iterator it = nodes.begin(); it != nodes.end(); ++it ) { mitk::DataNode::Pointer node = *it; if ( dynamic_cast(node->GetData()) ) m_SelectedFB.push_back(node); else if (dynamic_cast(node->GetData())) m_SelectedPF.push_back(node); else if (dynamic_cast(node->GetData())) m_SelectedImage = dynamic_cast(node->GetData()); else if (dynamic_cast(node->GetData())) m_Surfaces.push_back(dynamic_cast(node->GetData())); } UpdateGui(); GenerateStats(); } void QmitkFiberProcessingView::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 QmitkFiberProcessingView::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 QmitkFiberProcessingView::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 QmitkFiberProcessingView::AddFigureToDataStorage(mitk::PlanarFigure* figure, const QString& name, const char *propertyKey, mitk::BaseProperty *property ) { // initialize figure's geometry with empty geometry mitk::PlaneGeometry::Pointer emptygeometry = mitk::PlaneGeometry::New(); figure->SetGeometry2D( emptygeometry ); //set desired data to DataNode where Planarfigure is stored mitk::DataNode::Pointer newNode = mitk::DataNode::New(); newNode->SetName(name.toStdString()); newNode->SetData(figure); newNode->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 QmitkFiberProcessingView::DoFiberExtraction() { if ( m_SelectedFB.empty() ){ QMessageBox::information( NULL, "Warning", "No fibe bundle selected!"); MITK_WARN("QmitkFiberProcessingView") << "no fibe bundle selected"; return; } for (int i=0; i(m_SelectedFB.at(i)->GetData()); mitk::PlanarFigure::Pointer roi = dynamic_cast (m_SelectedPF.at(0)->GetData()); mitk::FiberBundleX::Pointer extFB = fib->ExtractFiberSubset(roi); mitk::DataNode::Pointer node; node = mitk::DataNode::New(); node->SetData(extFB); - QString name(m_SelectedFB.at(0)->GetName().c_str()); - name += "_extracted"; + QString name(m_SelectedFB.at(i)->GetName().c_str()); + name += "_"; + name += m_SelectedPF.at(0)->GetName().c_str(); node->SetName(name.toStdString()); GetDataStorage()->Add(node); + m_SelectedFB.at(i)->SetVisibility(false); } } void QmitkFiberProcessingView::GenerateAndComposite() { mitk::PlanarFigureComposite::Pointer PFCAnd = mitk::PlanarFigureComposite::New(); mitk::PlaneGeometry* currentGeometry2D = dynamic_cast( const_cast(GetActiveStdMultiWidget()->GetRenderWindow1()->GetRenderer()->GetCurrentWorldGeometry2D())); PFCAnd->SetGeometry2D(currentGeometry2D); 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"); } AddCompositeToDatastorage(PFCAnd, NULL); } void QmitkFiberProcessingView::GenerateOrComposite() { mitk::PlanarFigureComposite::Pointer PFCOr = mitk::PlanarFigureComposite::New(); mitk::PlaneGeometry* currentGeometry2D = dynamic_cast( const_cast(GetActiveStdMultiWidget()->GetRenderWindow1()->GetRenderer()->GetCurrentWorldGeometry2D())); PFCOr->SetGeometry2D(currentGeometry2D); 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"); } AddCompositeToDatastorage(PFCOr, NULL); } void QmitkFiberProcessingView::GenerateNotComposite() { mitk::PlanarFigureComposite::Pointer PFCNot = mitk::PlanarFigureComposite::New(); mitk::PlaneGeometry* currentGeometry2D = dynamic_cast( const_cast(GetActiveStdMultiWidget()->GetRenderWindow1()->GetRenderer()->GetCurrentWorldGeometry2D())); PFCNot->SetGeometry2D(currentGeometry2D); 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"); } AddCompositeToDatastorage(PFCNot, NULL); } /* CLEANUP NEEDED */ void QmitkFiberProcessingView::AddCompositeToDatastorage(mitk::PlanarFigureComposite::Pointer pfcomp, mitk::DataNode::Pointer parentDataNode ) { mitk::DataNode::Pointer newPFCNode; newPFCNode = mitk::DataNode::New(); newPFCNode->SetName( pfcomp->getDisplayName() ); newPFCNode->SetData(pfcomp); newPFCNode->SetVisibility(true); switch (pfcomp->getOperationType()) { case 0: { if (!parentDataNode.IsNull()) { GetDataStorage()->Add(newPFCNode, parentDataNode); } else { 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); AddCompositeToDatastorage(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: { if (!parentDataNode.IsNull()) { MITK_INFO << "adding " << newPFCNode->GetName() << " to " << parentDataNode->GetName() ; GetDataStorage()->Add(newPFCNode, parentDataNode); } else { MITK_INFO << "adding " << newPFCNode->GetName(); GetDataStorage()->Add(newPFCNode); } 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); AddCompositeToDatastorage(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: { 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); AddCompositeToDatastorage(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 QmitkFiberProcessingView::JoinBundles() { if ( m_SelectedFB.size()<2 ){ QMessageBox::information( NULL, "Warning", "Select at least two fiber bundles!"); MITK_WARN("QmitkFiberProcessingView") << "Select at least two fiber bundles!"; return; } std::vector::const_iterator it = m_SelectedFB.begin(); mitk::FiberBundleX::Pointer newBundle = dynamic_cast((*it)->GetData()); QString name(""); name += QString((*it)->GetName().c_str()); ++it; for (it; it!=m_SelectedFB.end(); ++it) { newBundle = *newBundle+dynamic_cast((*it)->GetData()); name += "+"+QString((*it)->GetName().c_str()); } mitk::DataNode::Pointer fbNode = mitk::DataNode::New(); fbNode->SetData(newBundle); fbNode->SetName(name.toStdString()); fbNode->SetVisibility(true); GetDataStorage()->Add(fbNode); } void QmitkFiberProcessingView::SubstractBundles() { if ( m_SelectedFB.size()<2 ){ QMessageBox::information( NULL, "Warning", "Select at least two fiber bundles!"); MITK_WARN("QmitkFiberProcessingView") << "Select at least two fiber bundles!"; return; } std::vector::const_iterator it = m_SelectedFB.begin(); mitk::FiberBundleX::Pointer newBundle = dynamic_cast((*it)->GetData()); QString name(""); name += QString((*it)->GetName().c_str()); ++it; for (it; it!=m_SelectedFB.end(); ++it) { newBundle = *newBundle-dynamic_cast((*it)->GetData()); name += "-"+QString((*it)->GetName().c_str()); } mitk::DataNode::Pointer fbNode = mitk::DataNode::New(); fbNode->SetData(newBundle); fbNode->SetName(name.toStdString()); fbNode->SetVisibility(true); GetDataStorage()->Add(fbNode); } void QmitkFiberProcessingView::GenerateStats() { if ( m_SelectedFB.empty() ) return; QString stats(""); for( int i=0; i(node->GetData())) { if (i>0) stats += "\n-----------------------------\n"; stats += QString(node->GetName().c_str()) + "\n"; mitk::FiberBundleX::Pointer fib = dynamic_cast(node->GetData()); vtkSmartPointer fiberPolyData = fib->GetFiberPolyData(); vtkSmartPointer vLines = fiberPolyData->GetLines(); vLines->InitTraversal(); int numberOfLines = vLines->GetNumberOfCells(); stats += "Number of fibers: "+ QString::number(numberOfLines) + "\n"; float length = 0; std::vector lengths; for (int i=0; iGetNextCell ( numPoints, points ); float l=0; for (unsigned int j=0; j p1; itk::Point p2; fiberPolyData->GetPoint(points[j], p1.GetDataPointer()); fiberPolyData->GetPoint(points[j+1], p2.GetDataPointer()); float dist = p1.EuclideanDistanceTo(p2); length += dist; l += dist; } itk::Point p2; fiberPolyData->GetPoint(points[numPoints-1], p2.GetDataPointer()); lengths.push_back(l); } std::sort(lengths.begin(), lengths.end()); if (numberOfLines>0) length /= numberOfLines; float dev=0; int count = 0; vLines->InitTraversal(); for (int i=0; iGetNextCell ( numPoints, points ); float l=0; for (unsigned int j=0; j p1; itk::Point p2; fiberPolyData->GetPoint(points[j], p1.GetDataPointer()); fiberPolyData->GetPoint(points[j+1], p2.GetDataPointer()); float dist = p1.EuclideanDistanceTo(p2); l += dist; } dev += (length-l)*(length-l); count++; } if (numberOfLines>1) dev /= (numberOfLines-1); else dev = 0; stats += "Min. length: "+ QString::number(lengths.front(),'f',1) + " mm\n"; stats += "Max. length: "+ QString::number(lengths.back(),'f',1) + " mm\n"; stats += "Mean length: "+ QString::number(length,'f',1) + " mm\n"; stats += "Median length: "+ QString::number(lengths.at(lengths.size()/2),'f',1) + " mm\n"; stats += "Standard deviation: "+ QString::number(sqrt(dev),'f',1) + " mm\n"; } } this->m_Controls->m_StatsTextEdit->setText(stats); } void QmitkFiberProcessingView::ProcessSelectedBundles() { if ( m_SelectedFB.empty() ){ QMessageBox::information( NULL, "Warning", "No fibe bundle selected!"); MITK_WARN("QmitkFiberProcessingView") << "no fibe bundle selected"; return; } int generationMethod = m_Controls->m_GenerationBox->currentIndex(); for( int i=0; i(node->GetData())) { mitk::FiberBundleX::Pointer fib = dynamic_cast(node->GetData()); QString name(node->GetName().c_str()); DataNode::Pointer newNode = NULL; switch(generationMethod){ case 0: newNode = GenerateTractDensityImage(fib, false); name += "_TDI"; break; case 1: newNode = GenerateTractDensityImage(fib, true); name += "_envelope"; break; case 2: newNode = GenerateColorHeatmap(fib); break; case 3: newNode = GenerateFiberEndingsImage(fib); name += "_fiber_endings"; break; case 4: newNode = GenerateFiberEndingsPointSet(fib); name += "_fiber_endings"; break; } if (newNode.IsNotNull()) { newNode->SetName(name.toStdString()); GetDataStorage()->Add(newNode); } } } } // generate pointset displaying the fiber endings mitk::DataNode::Pointer QmitkFiberProcessingView::GenerateFiberEndingsPointSet(mitk::FiberBundleX::Pointer fib) { mitk::PointSet::Pointer pointSet = mitk::PointSet::New(); vtkSmartPointer fiberPolyData = fib->GetFiberPolyData(); vtkSmartPointer vLines = fiberPolyData->GetLines(); vLines->InitTraversal(); int count = 0; int numFibers = fib->GetNumFibers(); for( int i=0; iGetNextCell ( numPoints, points ); if (numPoints>0) { double* point = fiberPolyData->GetPoint(points[0]); itk::Point itkPoint; itkPoint[0] = point[0]; itkPoint[1] = point[1]; itkPoint[2] = point[2]; pointSet->InsertPoint(count, itkPoint); count++; } if (numPoints>2) { double* point = fiberPolyData->GetPoint(points[numPoints-1]); itk::Point itkPoint; itkPoint[0] = point[0]; itkPoint[1] = point[1]; itkPoint[2] = point[2]; pointSet->InsertPoint(count, itkPoint); count++; } } mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData( pointSet ); return node; } // generate image displaying the fiber endings mitk::DataNode::Pointer QmitkFiberProcessingView::GenerateFiberEndingsImage(mitk::FiberBundleX::Pointer fib) { typedef unsigned char OutPixType; typedef itk::Image OutImageType; typedef itk::TractsToFiberEndingsImageFilter< OutImageType > ImageGeneratorType; ImageGeneratorType::Pointer generator = ImageGeneratorType::New(); generator->SetFiberBundle(fib); generator->SetInvertImage(m_Controls->m_InvertCheckbox->isChecked()); generator->SetUpsamplingFactor(m_Controls->m_UpsamplingSpinBox->value()); if (m_SelectedImage.IsNotNull()) { OutImageType::Pointer itkImage = OutImageType::New(); CastToItkImage(m_SelectedImage, itkImage); generator->SetInputImage(itkImage); generator->SetUseImageGeometry(true); } generator->Update(); // get output image OutImageType::Pointer outImg = generator->GetOutput(); mitk::Image::Pointer img = mitk::Image::New(); img->InitializeByItk(outImg.GetPointer()); img->SetVolume(outImg->GetBufferPointer()); // init data node mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData(img); return node; } // generate rgba heatmap from fiber bundle mitk::DataNode::Pointer QmitkFiberProcessingView::GenerateColorHeatmap(mitk::FiberBundleX::Pointer fib) { typedef itk::RGBAPixel OutPixType; typedef itk::Image OutImageType; typedef itk::TractsToRgbaImageFilter< OutImageType > ImageGeneratorType; ImageGeneratorType::Pointer generator = ImageGeneratorType::New(); generator->SetFiberBundle(fib); generator->SetUpsamplingFactor(m_Controls->m_UpsamplingSpinBox->value()); if (m_SelectedImage.IsNotNull()) { itk::Image::Pointer itkImage = itk::Image::New(); CastToItkImage(m_SelectedImage, itkImage); generator->SetInputImage(itkImage); generator->SetUseImageGeometry(true); } generator->Update(); // get output image typedef itk::Image OutType; OutType::Pointer outImg = generator->GetOutput(); mitk::Image::Pointer img = mitk::Image::New(); img->InitializeByItk(outImg.GetPointer()); img->SetVolume(outImg->GetBufferPointer()); // init data node mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData(img); return node; } // generate tract density image from fiber bundle mitk::DataNode::Pointer QmitkFiberProcessingView::GenerateTractDensityImage(mitk::FiberBundleX::Pointer fib, bool binary) { typedef float OutPixType; typedef itk::Image OutImageType; itk::TractDensityImageFilter< OutImageType >::Pointer generator = itk::TractDensityImageFilter< OutImageType >::New(); generator->SetFiberBundle(fib); generator->SetBinaryOutput(binary); generator->SetInvertImage(m_Controls->m_InvertCheckbox->isChecked()); generator->SetUpsamplingFactor(m_Controls->m_UpsamplingSpinBox->value()); if (m_SelectedImage.IsNotNull()) { OutImageType::Pointer itkImage = OutImageType::New(); CastToItkImage(m_SelectedImage, itkImage); generator->SetInputImage(itkImage); generator->SetUseImageGeometry(true); } generator->Update(); // get output image typedef itk::Image OutType; OutType::Pointer outImg = generator->GetOutput(); mitk::Image::Pointer img = mitk::Image::New(); img->InitializeByItk(outImg.GetPointer()); img->SetVolume(outImg->GetBufferPointer()); // init data node mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData(img); return node; } void QmitkFiberProcessingView::ResampleSelectedBundles() { int factor = this->m_Controls->m_ResampleFibersSpinBox->value(); for (int i=0; i(m_SelectedFB.at(i)->GetData()); fib->DoFiberSmoothing(factor); } } diff --git a/Modules/DiffusionImaging/IODataStructures/FiberBundleX/mitkFiberBundleX.cpp b/Modules/DiffusionImaging/IODataStructures/FiberBundleX/mitkFiberBundleX.cpp index a65fbc222c..451c6aeb39 100644 --- a/Modules/DiffusionImaging/IODataStructures/FiberBundleX/mitkFiberBundleX.cpp +++ b/Modules/DiffusionImaging/IODataStructures/FiberBundleX/mitkFiberBundleX.cpp @@ -1,1264 +1,1246 @@ /*========================================================================= 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. =========================================================================*/ #include "mitkFiberBundleX.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include const char* mitk::FiberBundleX::COLORCODING_ORIENTATION_BASED = "Color_Orient"; //const char* mitk::FiberBundleX::COLORCODING_FA_AS_OPACITY = "Color_Orient_FA_Opacity"; const char* mitk::FiberBundleX::COLORCODING_FA_BASED = "FA_Values"; const char* mitk::FiberBundleX::COLORCODING_CUSTOM = "custom"; const char* mitk::FiberBundleX::FIBER_ID_ARRAY = "Fiber_IDs"; mitk::FiberBundleX::FiberBundleX( vtkPolyData* fiberPolyData ) - : m_currentColorCoding(NULL) - , m_NumFibers(0) + : m_currentColorCoding(NULL) + , m_NumFibers(0) { m_FiberPolyData = vtkSmartPointer::New(); if (fiberPolyData != NULL) { m_FiberPolyData->DeepCopy(fiberPolyData); this->DoColorCodingOrientationbased(); } if(m_FiberPolyData->GetPointData()->HasArray(COLORCODING_ORIENTATION_BASED)) MITK_DEBUG << "ok"; vtkUnsignedCharArray* tmpColors = (vtkUnsignedCharArray*) m_FiberPolyData->GetPointData()->GetArray(COLORCODING_ORIENTATION_BASED); if (tmpColors!=NULL) { int tmpColorss = tmpColors->GetNumberOfTuples(); int tmpColorc = tmpColors->GetNumberOfComponents(); } m_NumFibers = m_FiberPolyData->GetNumberOfLines(); this->UpdateFiberGeometry(); this->SetColorCoding(COLORCODING_ORIENTATION_BASED); this->DoGenerateFiberIds(); } mitk::FiberBundleX::~FiberBundleX() { } mitk::FiberBundleX::Pointer mitk::FiberBundleX::GetDeepCopy() { - mitk::FiberBundleX::Pointer newFib = mitk::FiberBundleX::New(m_FiberPolyData); + mitk::FiberBundleX::Pointer newFib = mitk::FiberBundleX::New(m_FiberPolyData); - if(m_FiberPolyData->GetPointData()->HasArray(COLORCODING_ORIENTATION_BASED)) - MITK_DEBUG << "ok"; + if(m_FiberPolyData->GetPointData()->HasArray(COLORCODING_ORIENTATION_BASED)) + MITK_DEBUG << "ok"; - vtkUnsignedCharArray* tmpColors = (vtkUnsignedCharArray*) m_FiberPolyData->GetPointData()->GetArray(COLORCODING_ORIENTATION_BASED); - int tmpColorss = tmpColors->GetNumberOfTuples(); - int tmpColorc = tmpColors->GetNumberOfComponents(); + vtkUnsignedCharArray* tmpColors = (vtkUnsignedCharArray*) m_FiberPolyData->GetPointData()->GetArray(COLORCODING_ORIENTATION_BASED); + int tmpColorss = tmpColors->GetNumberOfTuples(); + int tmpColorc = tmpColors->GetNumberOfComponents(); - newFib->SetColorCoding(m_currentColorCoding); - return newFib; + newFib->SetColorCoding(m_currentColorCoding); + return newFib; } vtkSmartPointer mitk::FiberBundleX::GeneratePolyDataByIds(std::vector fiberIds) { - MITK_DEBUG << "\n=====FINAL RESULT: fib_id ======\n"; - MITK_DEBUG << "Number of new Fibers: " << fiberIds.size(); - // iterate through the vectorcontainer hosting all desired fiber Ids - - vtkSmartPointer newFiberPolyData = vtkSmartPointer::New(); - vtkSmartPointer newLineSet = vtkSmartPointer::New(); - vtkSmartPointer newPointSet = vtkSmartPointer::New(); - - // if FA array available, initialize fa double array - // if color orient array is available init color array - vtkSmartPointer faValueArray; - vtkSmartPointer colorsT; - //colors and alpha value for each single point, RGBA = 4 components - unsigned char rgba[4] = {0,0,0,0}; - int componentSize = sizeof(rgba); - - if (m_FiberIdDataSet->GetPointData()->HasArray(COLORCODING_FA_BASED)){ - MITK_DEBUG << "FA VALUES AVAILABLE, init array for new fiberbundle"; - faValueArray = vtkSmartPointer::New(); - } + MITK_DEBUG << "\n=====FINAL RESULT: fib_id ======\n"; + MITK_DEBUG << "Number of new Fibers: " << fiberIds.size(); + // iterate through the vectorcontainer hosting all desired fiber Ids + + vtkSmartPointer newFiberPolyData = vtkSmartPointer::New(); + vtkSmartPointer newLineSet = vtkSmartPointer::New(); + vtkSmartPointer newPointSet = vtkSmartPointer::New(); + + // if FA array available, initialize fa double array + // if color orient array is available init color array + vtkSmartPointer faValueArray; + vtkSmartPointer colorsT; + //colors and alpha value for each single point, RGBA = 4 components + unsigned char rgba[4] = {0,0,0,0}; + int componentSize = sizeof(rgba); + + if (m_FiberIdDataSet->GetPointData()->HasArray(COLORCODING_FA_BASED)){ + MITK_DEBUG << "FA VALUES AVAILABLE, init array for new fiberbundle"; + faValueArray = vtkSmartPointer::New(); + } - if (m_FiberIdDataSet->GetPointData()->HasArray(COLORCODING_ORIENTATION_BASED)){ - MITK_DEBUG << "colorValues available, init array for new fiberbundle"; - colorsT = vtkUnsignedCharArray::New(); - colorsT->SetNumberOfComponents(componentSize); - colorsT->SetName(COLORCODING_ORIENTATION_BASED); - } + if (m_FiberIdDataSet->GetPointData()->HasArray(COLORCODING_ORIENTATION_BASED)){ + MITK_DEBUG << "colorValues available, init array for new fiberbundle"; + colorsT = vtkUnsignedCharArray::New(); + colorsT->SetNumberOfComponents(componentSize); + colorsT->SetName(COLORCODING_ORIENTATION_BASED); + } - std::vector::iterator finIt = fiberIds.begin(); - while ( finIt != fiberIds.end() ) - { - if (*finIt < 0){ - MITK_DEBUG << "!!!!!! ERROR !!!!!! ERROR !!!!!\n=======================\nERROR, fiberID can not be negative!!! check id Extraction!" << *finIt; - break; - } - vtkSmartPointer fiber = m_FiberIdDataSet->GetCell(*finIt);//->DeepCopy(fiber); - vtkSmartPointer fibPoints = fiber->GetPoints(); + std::vector::iterator finIt = fiberIds.begin(); + while ( finIt != fiberIds.end() ) + { + if (*finIt < 0 || *finIt>GetNumFibers()){ + MITK_INFO << "FiberID can not be negative or >NumFibers!!! check id Extraction!" << *finIt; + break; + } - vtkSmartPointer newFiber = vtkSmartPointer::New(); - newFiber->GetPointIds()->SetNumberOfIds( fibPoints->GetNumberOfPoints() ); + vtkSmartPointer fiber = m_FiberIdDataSet->GetCell(*finIt);//->DeepCopy(fiber); - for(int i=0; iGetNumberOfPoints(); i++) - { - // MITK_DEBUG << "id: " << fiber->GetPointId(i); - // MITK_DEBUG << fibPoints->GetPoint(i)[0] << " | " << fibPoints->GetPoint(i)[1] << " | " << fibPoints->GetPoint(i)[2]; - newFiber->GetPointIds()->SetId(i, newPointSet->GetNumberOfPoints()); - newPointSet->InsertNextPoint(fibPoints->GetPoint(i)[0], fibPoints->GetPoint(i)[1], fibPoints->GetPoint(i)[2]); + vtkSmartPointer fibPoints = fiber->GetPoints(); + vtkSmartPointer newFiber = vtkSmartPointer::New(); + newFiber->GetPointIds()->SetNumberOfIds( fibPoints->GetNumberOfPoints() ); - if (m_FiberIdDataSet->GetPointData()->HasArray(COLORCODING_FA_BASED)){ - // MITK_DEBUG << m_FiberIdDataSet->GetPointData()->GetArray(FA_VALUE_ARRAY)->GetTuple(fiber->GetPointId(i)); - } + for(int i=0; iGetNumberOfPoints(); i++) + { + // MITK_DEBUG << "id: " << fiber->GetPointId(i); + // MITK_DEBUG << fibPoints->GetPoint(i)[0] << " | " << fibPoints->GetPoint(i)[1] << " | " << fibPoints->GetPoint(i)[2]; + newFiber->GetPointIds()->SetId(i, newPointSet->GetNumberOfPoints()); + newPointSet->InsertNextPoint(fibPoints->GetPoint(i)[0], fibPoints->GetPoint(i)[1], fibPoints->GetPoint(i)[2]); - if (m_FiberIdDataSet->GetPointData()->HasArray(COLORCODING_ORIENTATION_BASED)){ - // MITK_DEBUG << "ColorValue: " << m_FiberIdDataSet->GetPointData()->GetArray(COLORCODING_ORIENTATION_BASED)->GetTuple(fiber->GetPointId(i))[0]; - } - } - newLineSet->InsertNextCell(newFiber); - ++finIt; + if (m_FiberIdDataSet->GetPointData()->HasArray(COLORCODING_FA_BASED)){ + // MITK_DEBUG << m_FiberIdDataSet->GetPointData()->GetArray(FA_VALUE_ARRAY)->GetTuple(fiber->GetPointId(i)); + } + + if (m_FiberIdDataSet->GetPointData()->HasArray(COLORCODING_ORIENTATION_BASED)){ + // MITK_DEBUG << "ColorValue: " << m_FiberIdDataSet->GetPointData()->GetArray(COLORCODING_ORIENTATION_BASED)->GetTuple(fiber->GetPointId(i))[0]; + } } - newFiberPolyData->SetPoints(newPointSet); - newFiberPolyData->SetLines(newLineSet); - MITK_DEBUG << "new fiberbundle polydata points: " << newFiberPolyData->GetNumberOfPoints(); - MITK_DEBUG << "new fiberbundle polydata lines: " << newFiberPolyData->GetNumberOfLines(); - MITK_DEBUG << "=====================\n"; + newLineSet->InsertNextCell(newFiber); + ++finIt; + } - // mitk::FiberBundleX::Pointer newFib = mitk::FiberBundleX::New(newFiberPolyData); - return newFiberPolyData; + newFiberPolyData->SetPoints(newPointSet); + newFiberPolyData->SetLines(newLineSet); + MITK_DEBUG << "new fiberbundle polydata points: " << newFiberPolyData->GetNumberOfPoints(); + MITK_DEBUG << "new fiberbundle polydata lines: " << newFiberPolyData->GetNumberOfLines(); + MITK_DEBUG << "=====================\n"; + + // mitk::FiberBundleX::Pointer newFib = mitk::FiberBundleX::New(newFiberPolyData); + return newFiberPolyData; } // merge two fiber bundles mitk::FiberBundleX::Pointer mitk::FiberBundleX::operator+(mitk::FiberBundleX* fib) { - vtkSmartPointer vNewPolyData = vtkSmartPointer::New(); - vtkSmartPointer vNewLines = vtkSmartPointer::New(); - vtkSmartPointer vNewPoints = vtkSmartPointer::New(); + vtkSmartPointer vNewPolyData = vtkSmartPointer::New(); + vtkSmartPointer vNewLines = vtkSmartPointer::New(); + vtkSmartPointer vNewPoints = vtkSmartPointer::New(); + + vtkSmartPointer vLines = m_FiberPolyData->GetLines(); + vLines->InitTraversal(); - vtkSmartPointer vLines = m_FiberPolyData->GetLines(); - vLines->InitTraversal(); + // add current fiber bundle + int numFibers = GetNumFibers(); + for( int i=0; iGetNextCell ( numPoints, points ); - // add current fiber bundle - int numFibers = GetNumFibers(); - for( int i=0; i container = vtkSmartPointer::New(); + for( int j=0; jGetNextCell ( numPoints, points ); - - vtkSmartPointer container = vtkSmartPointer::New(); - for( int j=0; jInsertNextPoint(m_FiberPolyData->GetPoint(points[j])); - container->GetPointIds()->InsertNextId(id); - } - vNewLines->InsertNextCell(container); + vtkIdType id = vNewPoints->InsertNextPoint(m_FiberPolyData->GetPoint(points[j])); + container->GetPointIds()->InsertNextId(id); } + vNewLines->InsertNextCell(container); + } - vLines = fib->m_FiberPolyData->GetLines(); - vLines->InitTraversal(); + vLines = fib->m_FiberPolyData->GetLines(); + vLines->InitTraversal(); - // add new fiber bundle - numFibers = fib->GetNumFibers(); - for( int i=0; iGetNextCell ( numPoints, points ); + // add new fiber bundle + numFibers = fib->GetNumFibers(); + for( int i=0; iGetNextCell ( numPoints, points ); - vtkSmartPointer container = vtkSmartPointer::New(); - for( int j=0; jInsertNextPoint(fib->m_FiberPolyData->GetPoint(points[j])); - container->GetPointIds()->InsertNextId(id); - } - vNewLines->InsertNextCell(container); + vtkSmartPointer container = vtkSmartPointer::New(); + for( int j=0; jInsertNextPoint(fib->m_FiberPolyData->GetPoint(points[j])); + container->GetPointIds()->InsertNextId(id); } + vNewLines->InsertNextCell(container); + } - // initialize polydata - vNewPolyData->SetPoints(vNewPoints); - vNewPolyData->SetLines(vNewLines); + // initialize polydata + vNewPolyData->SetPoints(vNewPoints); + vNewPolyData->SetLines(vNewLines); - // initialize fiber bundle - mitk::FiberBundleX::Pointer newFib = mitk::FiberBundleX::New(vNewPolyData); - return newFib; + // initialize fiber bundle + mitk::FiberBundleX::Pointer newFib = mitk::FiberBundleX::New(vNewPolyData); + return newFib; } // subtract two fiber bundles mitk::FiberBundleX::Pointer mitk::FiberBundleX::operator-(mitk::FiberBundleX* fib) { - vtkSmartPointer vNewPolyData = vtkSmartPointer::New(); - vtkSmartPointer vNewLines = vtkSmartPointer::New(); - vtkSmartPointer vNewPoints = vtkSmartPointer::New(); + vtkSmartPointer vNewPolyData = vtkSmartPointer::New(); + vtkSmartPointer vNewLines = vtkSmartPointer::New(); + vtkSmartPointer vNewPoints = vtkSmartPointer::New(); - vtkSmartPointer vLines = m_FiberPolyData->GetLines(); - vLines->InitTraversal(); + vtkSmartPointer vLines = m_FiberPolyData->GetLines(); + vLines->InitTraversal(); - // iterate over current fibers - int numFibers = GetNumFibers(); - for( int i=0; iGetNextCell ( numPoints, points ); + + vtkSmartPointer vLines2 = fib->m_FiberPolyData->GetLines(); + vLines2->InitTraversal(); + int numFibers2 = fib->GetNumFibers(); + bool contained = false; + for( int i2=0; i2GetNextCell ( numPoints, points ); - - vtkSmartPointer vLines2 = fib->m_FiberPolyData->GetLines(); - vLines2->InitTraversal(); - int numFibers2 = fib->GetNumFibers(); - bool contained = false; - for( int i2=0; i2GetNextCell ( numPoints2, points2 ); - - // check endpoints - itk::Point point_start = GetItkPoint(m_FiberPolyData->GetPoint(points[0])); - itk::Point point_end = GetItkPoint(m_FiberPolyData->GetPoint(points[numPoints-1])); - itk::Point point2_start = GetItkPoint(fib->m_FiberPolyData->GetPoint(points2[0])); - itk::Point point2_end = GetItkPoint(fib->m_FiberPolyData->GetPoint(points2[numPoints2-1])); - - if (point_start.SquaredEuclideanDistanceTo(point2_start)<=mitk::eps && point_end.SquaredEuclideanDistanceTo(point2_end)<=mitk::eps || - point_start.SquaredEuclideanDistanceTo(point2_end)<=mitk::eps && point_end.SquaredEuclideanDistanceTo(point2_start)<=mitk::eps) - { - // further checking ??? - contained = true; - } - } + vtkIdType numPoints2(0); + vtkIdType* points2(NULL); + vLines2->GetNextCell ( numPoints2, points2 ); + + // check endpoints + itk::Point point_start = GetItkPoint(m_FiberPolyData->GetPoint(points[0])); + itk::Point point_end = GetItkPoint(m_FiberPolyData->GetPoint(points[numPoints-1])); + itk::Point point2_start = GetItkPoint(fib->m_FiberPolyData->GetPoint(points2[0])); + itk::Point point2_end = GetItkPoint(fib->m_FiberPolyData->GetPoint(points2[numPoints2-1])); + + if (point_start.SquaredEuclideanDistanceTo(point2_start)<=mitk::eps && point_end.SquaredEuclideanDistanceTo(point2_end)<=mitk::eps || + point_start.SquaredEuclideanDistanceTo(point2_end)<=mitk::eps && point_end.SquaredEuclideanDistanceTo(point2_start)<=mitk::eps) + { + // further checking ??? + contained = true; + } + } - // add to result because fiber is not subtracted - if (!contained) - { - vtkSmartPointer container = vtkSmartPointer::New(); - for( int j=0; jInsertNextPoint(m_FiberPolyData->GetPoint(points[j])); - container->GetPointIds()->InsertNextId(id); - } - vNewLines->InsertNextCell(container); - } + // add to result because fiber is not subtracted + if (!contained) + { + vtkSmartPointer container = vtkSmartPointer::New(); + for( int j=0; jInsertNextPoint(m_FiberPolyData->GetPoint(points[j])); + container->GetPointIds()->InsertNextId(id); + } + vNewLines->InsertNextCell(container); } + } - // initialize polydata - vNewPolyData->SetPoints(vNewPoints); - vNewPolyData->SetLines(vNewLines); + // initialize polydata + vNewPolyData->SetPoints(vNewPoints); + vNewPolyData->SetLines(vNewLines); - // initialize fiber bundle - mitk::FiberBundleX::Pointer newFib = mitk::FiberBundleX::New(vNewPolyData); - return newFib; + // initialize fiber bundle + mitk::FiberBundleX::Pointer newFib = mitk::FiberBundleX::New(vNewPolyData); + return newFib; } itk::Point mitk::FiberBundleX::GetItkPoint(double point[3]) { - itk::Point itkPoint; - itkPoint[0] = point[0]; - itkPoint[1] = point[1]; - itkPoint[2] = point[2]; - return itkPoint; + itk::Point itkPoint; + itkPoint[0] = point[0]; + itkPoint[1] = point[1]; + itkPoint[2] = point[2]; + return itkPoint; } /* * set polydata (additional flag to recompute fiber geometry, default = true) */ void mitk::FiberBundleX::SetFiberPolyData(vtkSmartPointer fiberPD, bool updateGeometry) { - if (fiberPD == NULL) - this->m_FiberPolyData = vtkSmartPointer::New(); - else - { - m_FiberPolyData->DeepCopy(fiberPD); - DoColorCodingOrientationbased(); - } + if (fiberPD == NULL) + this->m_FiberPolyData = vtkSmartPointer::New(); + else + { + m_FiberPolyData->DeepCopy(fiberPD); + DoColorCodingOrientationbased(); + } - m_NumFibers = m_FiberPolyData->GetNumberOfLines(); + m_NumFibers = m_FiberPolyData->GetNumberOfLines(); - if (updateGeometry) - UpdateFiberGeometry(); - SetColorCoding(COLORCODING_ORIENTATION_BASED); - DoGenerateFiberIds(); + if (updateGeometry) + UpdateFiberGeometry(); + SetColorCoding(COLORCODING_ORIENTATION_BASED); + DoGenerateFiberIds(); } /* * return vtkPolyData */ vtkSmartPointer mitk::FiberBundleX::GetFiberPolyData() { - return m_FiberPolyData; + return m_FiberPolyData; } void mitk::FiberBundleX::DoColorCodingOrientationbased() { - //===== FOR WRITING A TEST ======================== - // colorT size == tupelComponents * tupelElements - // compare color results - // to cover this code 100% also polydata needed, where colorarray already exists - // + one fiber with exactly 1 point - // + one fiber with 0 points - //================================================= - - - /* make sure that processing colorcoding is only called when necessary */ - if ( m_FiberPolyData->GetPointData()->HasArray(COLORCODING_ORIENTATION_BASED) && - m_FiberPolyData->GetNumberOfPoints() == - m_FiberPolyData->GetPointData()->GetArray(COLORCODING_ORIENTATION_BASED)->GetNumberOfTuples() ) - { - // fiberstructure is already colorcoded - MITK_DEBUG << " NO NEED TO REGENERATE COLORCODING! " ; - this->ResetFiberColorOpacity(); - this->SetColorCoding(COLORCODING_ORIENTATION_BASED); - return; - } + //===== FOR WRITING A TEST ======================== + // colorT size == tupelComponents * tupelElements + // compare color results + // to cover this code 100% also polydata needed, where colorarray already exists + // + one fiber with exactly 1 point + // + one fiber with 0 points + //================================================= + + + /* make sure that processing colorcoding is only called when necessary */ + if ( m_FiberPolyData->GetPointData()->HasArray(COLORCODING_ORIENTATION_BASED) && + m_FiberPolyData->GetNumberOfPoints() == + m_FiberPolyData->GetPointData()->GetArray(COLORCODING_ORIENTATION_BASED)->GetNumberOfTuples() ) + { + // fiberstructure is already colorcoded + MITK_DEBUG << " NO NEED TO REGENERATE COLORCODING! " ; + this->ResetFiberColorOpacity(); + this->SetColorCoding(COLORCODING_ORIENTATION_BASED); + return; + } - /* Finally, execute color calculation */ - vtkPoints* extrPoints = NULL; - extrPoints = m_FiberPolyData->GetPoints(); - int numOfPoints = 0; - if (extrPoints!=NULL) - numOfPoints = extrPoints->GetNumberOfPoints(); + /* Finally, execute color calculation */ + vtkPoints* extrPoints = NULL; + extrPoints = m_FiberPolyData->GetPoints(); + int numOfPoints = 0; + if (extrPoints!=NULL) + numOfPoints = extrPoints->GetNumberOfPoints(); - //colors and alpha value for each single point, RGBA = 4 components - unsigned char rgba[4] = {0,0,0,0}; - // int componentSize = sizeof(rgba); - int componentSize = 4; + //colors and alpha value for each single point, RGBA = 4 components + unsigned char rgba[4] = {0,0,0,0}; + // int componentSize = sizeof(rgba); + int componentSize = 4; - vtkSmartPointer colorsT = vtkUnsignedCharArray::New(); - colorsT->Allocate(numOfPoints * componentSize); - colorsT->SetNumberOfComponents(componentSize); - colorsT->SetName(COLORCODING_ORIENTATION_BASED); + vtkSmartPointer colorsT = vtkUnsignedCharArray::New(); + colorsT->Allocate(numOfPoints * componentSize); + colorsT->SetNumberOfComponents(componentSize); + colorsT->SetName(COLORCODING_ORIENTATION_BASED); - /* checkpoint: does polydata contain any fibers */ - int numOfFibers = m_FiberPolyData->GetNumberOfLines(); - if (numOfFibers < 1) { - MITK_DEBUG << "\n ========= Number of Fibers is 0 and below ========= \n"; - return; - } + /* checkpoint: does polydata contain any fibers */ + int numOfFibers = m_FiberPolyData->GetNumberOfLines(); + if (numOfFibers < 1) { + MITK_DEBUG << "\n ========= Number of Fibers is 0 and below ========= \n"; + return; + } - /* extract single fibers of fiberBundle */ - vtkCellArray* fiberList = m_FiberPolyData->GetLines(); - fiberList->InitTraversal(); - for (int fi=0; fiGetLines(); + fiberList->InitTraversal(); + for (int fi=0; fiGetNextCell(pointsPerFiber, idList); + vtkIdType* idList; // contains the point id's of the line + vtkIdType pointsPerFiber; // number of points for current line + fiberList->GetNextCell(pointsPerFiber, idList); - // MITK_DEBUG << "Fib#: " << fi << " of " << numOfFibers << " pnts in fiber: " << pointsPerFiber ; + // MITK_DEBUG << "Fib#: " << fi << " of " << numOfFibers << " pnts in fiber: " << pointsPerFiber ; - /* single fiber checkpoints: is number of points valid */ - if (pointsPerFiber > 1) - { - /* operate on points of single fiber */ - for (int i=0; i 1) + { + /* operate on points of single fiber */ + for (int i=0; i 0) - { - /* The color value of the current point is influenced by the previous point and next point. */ - vnl_vector_fixed< double, 3 > currentPntvtk(extrPoints->GetPoint(idList[i])[0], extrPoints->GetPoint(idList[i])[1],extrPoints->GetPoint(idList[i])[2]); - vnl_vector_fixed< double, 3 > nextPntvtk(extrPoints->GetPoint(idList[i+1])[0], extrPoints->GetPoint(idList[i+1])[1], extrPoints->GetPoint(idList[i+1])[2]); - vnl_vector_fixed< double, 3 > prevPntvtk(extrPoints->GetPoint(idList[i-1])[0], extrPoints->GetPoint(idList[i-1])[1], extrPoints->GetPoint(idList[i-1])[2]); + if (i 0) + { + /* The color value of the current point is influenced by the previous point and next point. */ + vnl_vector_fixed< double, 3 > currentPntvtk(extrPoints->GetPoint(idList[i])[0], extrPoints->GetPoint(idList[i])[1],extrPoints->GetPoint(idList[i])[2]); + vnl_vector_fixed< double, 3 > nextPntvtk(extrPoints->GetPoint(idList[i+1])[0], extrPoints->GetPoint(idList[i+1])[1], extrPoints->GetPoint(idList[i+1])[2]); + vnl_vector_fixed< double, 3 > prevPntvtk(extrPoints->GetPoint(idList[i-1])[0], extrPoints->GetPoint(idList[i-1])[1], extrPoints->GetPoint(idList[i-1])[2]); - vnl_vector_fixed< double, 3 > diff1; - diff1 = currentPntvtk - nextPntvtk; + vnl_vector_fixed< double, 3 > diff1; + diff1 = currentPntvtk - nextPntvtk; - vnl_vector_fixed< double, 3 > diff2; - diff2 = currentPntvtk - prevPntvtk; + vnl_vector_fixed< double, 3 > diff2; + diff2 = currentPntvtk - prevPntvtk; - vnl_vector_fixed< double, 3 > diff; - diff = (diff1 - diff2) / 2.0; - diff.normalize(); + vnl_vector_fixed< double, 3 > diff; + diff = (diff1 - diff2) / 2.0; + diff.normalize(); - rgba[0] = (unsigned char) (255.0 * std::fabs(diff[0])); - rgba[1] = (unsigned char) (255.0 * std::fabs(diff[1])); - rgba[2] = (unsigned char) (255.0 * std::fabs(diff[2])); - rgba[3] = (unsigned char) (255.0); + rgba[0] = (unsigned char) (255.0 * std::fabs(diff[0])); + rgba[1] = (unsigned char) (255.0 * std::fabs(diff[1])); + rgba[2] = (unsigned char) (255.0 * std::fabs(diff[2])); + rgba[3] = (unsigned char) (255.0); - } else if (i==0) { - /* First point has no previous point, therefore only diff1 is taken */ + } else if (i==0) { + /* First point has no previous point, therefore only diff1 is taken */ - vnl_vector_fixed< double, 3 > currentPntvtk(extrPoints->GetPoint(idList[i])[0], extrPoints->GetPoint(idList[i])[1],extrPoints->GetPoint(idList[i])[2]); - vnl_vector_fixed< double, 3 > nextPntvtk(extrPoints->GetPoint(idList[i+1])[0], extrPoints->GetPoint(idList[i+1])[1], extrPoints->GetPoint(idList[i+1])[2]); + vnl_vector_fixed< double, 3 > currentPntvtk(extrPoints->GetPoint(idList[i])[0], extrPoints->GetPoint(idList[i])[1],extrPoints->GetPoint(idList[i])[2]); + vnl_vector_fixed< double, 3 > nextPntvtk(extrPoints->GetPoint(idList[i+1])[0], extrPoints->GetPoint(idList[i+1])[1], extrPoints->GetPoint(idList[i+1])[2]); - vnl_vector_fixed< double, 3 > diff1; - diff1 = currentPntvtk - nextPntvtk; - diff1.normalize(); + vnl_vector_fixed< double, 3 > diff1; + diff1 = currentPntvtk - nextPntvtk; + diff1.normalize(); - rgba[0] = (unsigned char) (255.0 * std::fabs(diff1[0])); - rgba[1] = (unsigned char) (255.0 * std::fabs(diff1[1])); - rgba[2] = (unsigned char) (255.0 * std::fabs(diff1[2])); - rgba[3] = (unsigned char) (255.0); + rgba[0] = (unsigned char) (255.0 * std::fabs(diff1[0])); + rgba[1] = (unsigned char) (255.0 * std::fabs(diff1[1])); + rgba[2] = (unsigned char) (255.0 * std::fabs(diff1[2])); + rgba[3] = (unsigned char) (255.0); - } else if (i==pointsPerFiber-1) { - /* Last point has no next point, therefore only diff2 is taken */ - vnl_vector_fixed< double, 3 > currentPntvtk(extrPoints->GetPoint(idList[i])[0], extrPoints->GetPoint(idList[i])[1],extrPoints->GetPoint(idList[i])[2]); - vnl_vector_fixed< double, 3 > prevPntvtk(extrPoints->GetPoint(idList[i-1])[0], extrPoints->GetPoint(idList[i-1])[1], extrPoints->GetPoint(idList[i-1])[2]); + } else if (i==pointsPerFiber-1) { + /* Last point has no next point, therefore only diff2 is taken */ + vnl_vector_fixed< double, 3 > currentPntvtk(extrPoints->GetPoint(idList[i])[0], extrPoints->GetPoint(idList[i])[1],extrPoints->GetPoint(idList[i])[2]); + vnl_vector_fixed< double, 3 > prevPntvtk(extrPoints->GetPoint(idList[i-1])[0], extrPoints->GetPoint(idList[i-1])[1], extrPoints->GetPoint(idList[i-1])[2]); - vnl_vector_fixed< double, 3 > diff2; - diff2 = currentPntvtk - prevPntvtk; - diff2.normalize(); + vnl_vector_fixed< double, 3 > diff2; + diff2 = currentPntvtk - prevPntvtk; + diff2.normalize(); - rgba[0] = (unsigned char) (255.0 * std::fabs(diff2[0])); - rgba[1] = (unsigned char) (255.0 * std::fabs(diff2[1])); - rgba[2] = (unsigned char) (255.0 * std::fabs(diff2[2])); - rgba[3] = (unsigned char) (255.0); + rgba[0] = (unsigned char) (255.0 * std::fabs(diff2[0])); + rgba[1] = (unsigned char) (255.0 * std::fabs(diff2[1])); + rgba[2] = (unsigned char) (255.0 * std::fabs(diff2[2])); + rgba[3] = (unsigned char) (255.0); - } + } - colorsT->InsertTupleValue(idList[i], rgba); + colorsT->InsertTupleValue(idList[i], rgba); - } //end for loop + } //end for loop - } else if (pointsPerFiber == 1) { - /* a single point does not define a fiber (use vertex mechanisms instead */ - continue; - // colorsT->InsertTupleValue(0, rgba); + } else if (pointsPerFiber == 1) { + /* a single point does not define a fiber (use vertex mechanisms instead */ + continue; + // colorsT->InsertTupleValue(0, rgba); - } else { - MITK_DEBUG << "Fiber with 0 points detected... please check your tractography algorithm!" ; - continue; + } else { + MITK_DEBUG << "Fiber with 0 points detected... please check your tractography algorithm!" ; + continue; - } + } - }//end for loop + }//end for loop - m_FiberPolyData->GetPointData()->AddArray(colorsT); + m_FiberPolyData->GetPointData()->AddArray(colorsT); - /*========================= + /*========================= - this is more relevant for renderer than for fiberbundleX datastructure - think about sourcing this to a explicit method which coordinates colorcoding */ - this->SetColorCoding(COLORCODING_ORIENTATION_BASED); - // =========================== + this->SetColorCoding(COLORCODING_ORIENTATION_BASED); + // =========================== - //mini test, shall be ported to MITK TESTINGS! - if (colorsT->GetSize() != numOfPoints*componentSize) - MITK_DEBUG << "ALLOCATION ERROR IN INITIATING COLOR ARRAY"; + //mini test, shall be ported to MITK TESTINGS! + if (colorsT->GetSize() != numOfPoints*componentSize) + MITK_DEBUG << "ALLOCATION ERROR IN INITIATING COLOR ARRAY"; } void mitk::FiberBundleX::DoColorCodingFAbased() { - if(m_FiberPolyData->GetPointData()->HasArray(COLORCODING_FA_BASED) != 1 ) - return; + if(m_FiberPolyData->GetPointData()->HasArray(COLORCODING_FA_BASED) != 1 ) + return; - this->SetColorCoding(COLORCODING_FA_BASED); - MITK_DEBUG << "FBX: done CC FA based"; - this->DoGenerateFiberIds(); + this->SetColorCoding(COLORCODING_FA_BASED); + MITK_DEBUG << "FBX: done CC FA based"; + this->DoGenerateFiberIds(); } void mitk::FiberBundleX::DoUseFAasColorOpacity() { - if(m_FiberPolyData->GetPointData()->HasArray(COLORCODING_FA_BASED) != 1 ) - return; + if(m_FiberPolyData->GetPointData()->HasArray(COLORCODING_FA_BASED) != 1 ) + return; - if(m_FiberPolyData->GetPointData()->HasArray(COLORCODING_ORIENTATION_BASED) != 1 ) - return; + if(m_FiberPolyData->GetPointData()->HasArray(COLORCODING_ORIENTATION_BASED) != 1 ) + return; - vtkDoubleArray* FAValArray = (vtkDoubleArray*) m_FiberPolyData->GetPointData()->GetArray(COLORCODING_FA_BASED); - vtkUnsignedCharArray* ColorArray = dynamic_cast (m_FiberPolyData->GetPointData()->GetArray(COLORCODING_ORIENTATION_BASED)); + vtkDoubleArray* FAValArray = (vtkDoubleArray*) m_FiberPolyData->GetPointData()->GetArray(COLORCODING_FA_BASED); + vtkUnsignedCharArray* ColorArray = dynamic_cast (m_FiberPolyData->GetPointData()->GetArray(COLORCODING_ORIENTATION_BASED)); - for(long i=0; iGetNumberOfTuples(); i++) { - double faValue = FAValArray->GetValue(i); - faValue = faValue * 255.0; - ColorArray->SetComponent(i,3, (unsigned char) faValue ); - } + for(long i=0; iGetNumberOfTuples(); i++) { + double faValue = FAValArray->GetValue(i); + faValue = faValue * 255.0; + ColorArray->SetComponent(i,3, (unsigned char) faValue ); + } - this->SetColorCoding(COLORCODING_ORIENTATION_BASED); - MITK_DEBUG << "FBX: done CC OPACITY"; - this->DoGenerateFiberIds(); + this->SetColorCoding(COLORCODING_ORIENTATION_BASED); + MITK_DEBUG << "FBX: done CC OPACITY"; + this->DoGenerateFiberIds(); } void mitk::FiberBundleX::ResetFiberColorOpacity() { - vtkUnsignedCharArray* ColorArray = dynamic_cast (m_FiberPolyData->GetPointData()->GetArray(COLORCODING_ORIENTATION_BASED)); - for(long i=0; iGetNumberOfTuples(); i++) { - ColorArray->SetComponent(i,3, 255.0 ); - } + vtkUnsignedCharArray* ColorArray = dynamic_cast (m_FiberPolyData->GetPointData()->GetArray(COLORCODING_ORIENTATION_BASED)); + for(long i=0; iGetNumberOfTuples(); i++) { + ColorArray->SetComponent(i,3, 255.0 ); + } } void mitk::FiberBundleX::SetFAMap(mitk::Image::Pointer FAimage) { - MITK_DEBUG << "SetFAMap"; - vtkSmartPointer faValues = vtkDoubleArray::New(); - faValues->SetName(COLORCODING_FA_BASED); - faValues->Allocate(m_FiberPolyData->GetNumberOfPoints()); - // MITK_DEBUG << faValues->GetNumberOfTuples(); - // MITK_DEBUG << faValues->GetSize(); - - faValues->SetNumberOfValues(m_FiberPolyData->GetNumberOfPoints()); - // MITK_DEBUG << faValues->GetNumberOfTuples(); - // MITK_DEBUG << faValues->GetSize(); - - vtkPoints* pointSet = m_FiberPolyData->GetPoints(); - for(long i=0; iGetNumberOfPoints(); ++i) - { - Point3D px; - px[0] = pointSet->GetPoint(i)[0]; - px[1] = pointSet->GetPoint(i)[1]; - px[2] = pointSet->GetPoint(i)[2]; - double faPixelValue = FAimage->GetPixelValueByWorldCoordinate(px) * 0.01; - // faValues->InsertNextTuple1(faPixelValue); - faValues->InsertValue(i, faPixelValue); - // MITK_DEBUG << faPixelValue; - // MITK_DEBUG << faValues->GetValue(i); + MITK_DEBUG << "SetFAMap"; + vtkSmartPointer faValues = vtkDoubleArray::New(); + faValues->SetName(COLORCODING_FA_BASED); + faValues->Allocate(m_FiberPolyData->GetNumberOfPoints()); + // MITK_DEBUG << faValues->GetNumberOfTuples(); + // MITK_DEBUG << faValues->GetSize(); + + faValues->SetNumberOfValues(m_FiberPolyData->GetNumberOfPoints()); + // MITK_DEBUG << faValues->GetNumberOfTuples(); + // MITK_DEBUG << faValues->GetSize(); + + vtkPoints* pointSet = m_FiberPolyData->GetPoints(); + for(long i=0; iGetNumberOfPoints(); ++i) + { + Point3D px; + px[0] = pointSet->GetPoint(i)[0]; + px[1] = pointSet->GetPoint(i)[1]; + px[2] = pointSet->GetPoint(i)[2]; + double faPixelValue = FAimage->GetPixelValueByWorldCoordinate(px) * 0.01; + // faValues->InsertNextTuple1(faPixelValue); + faValues->InsertValue(i, faPixelValue); + // MITK_DEBUG << faPixelValue; + // MITK_DEBUG << faValues->GetValue(i); - } + } - m_FiberPolyData->GetPointData()->AddArray(faValues); - this->DoGenerateFiberIds(); + m_FiberPolyData->GetPointData()->AddArray(faValues); + this->DoGenerateFiberIds(); - if(m_FiberPolyData->GetPointData()->HasArray(COLORCODING_FA_BASED)) - MITK_DEBUG << "FA VALUE ARRAY SET"; + if(m_FiberPolyData->GetPointData()->HasArray(COLORCODING_FA_BASED)) + MITK_DEBUG << "FA VALUE ARRAY SET"; - // vtkDoubleArray* valueArray = (vtkDoubleArray*) m_FiberPolyData->GetPointData()->GetArray(FA_VALUE_ARRAY); - // for(long i=0; iGetNumberOfPoints(); i++) - // { - // MITK_DEBUG << "value at pos "<< i << ": " << valueArray->GetValue(i); - // } + // vtkDoubleArray* valueArray = (vtkDoubleArray*) m_FiberPolyData->GetPointData()->GetArray(FA_VALUE_ARRAY); + // for(long i=0; iGetNumberOfPoints(); i++) + // { + // MITK_DEBUG << "value at pos "<< i << ": " << valueArray->GetValue(i); + // } } void mitk::FiberBundleX::DoGenerateFiberIds() { - if (m_FiberPolyData == NULL) - return; + if (m_FiberPolyData == NULL) + return; - vtkSmartPointer idFiberFilter = vtkSmartPointer::New(); - idFiberFilter->SetInput(m_FiberPolyData); - idFiberFilter->CellIdsOn(); - // idFiberFilter->PointIdsOn(); // point id's are not needed - idFiberFilter->SetIdsArrayName(FIBER_ID_ARRAY); - idFiberFilter->FieldDataOn(); - idFiberFilter->Update(); + vtkSmartPointer idFiberFilter = vtkSmartPointer::New(); + idFiberFilter->SetInput(m_FiberPolyData); + idFiberFilter->CellIdsOn(); + // idFiberFilter->PointIdsOn(); // point id's are not needed + idFiberFilter->SetIdsArrayName(FIBER_ID_ARRAY); + idFiberFilter->FieldDataOn(); + idFiberFilter->Update(); - m_FiberIdDataSet = idFiberFilter->GetOutput(); + m_FiberIdDataSet = idFiberFilter->GetOutput(); - MITK_DEBUG << "Generating Fiber Ids...[done] | " << m_FiberIdDataSet->GetNumberOfCells(); + MITK_DEBUG << "Generating Fiber Ids...[done] | " << m_FiberIdDataSet->GetNumberOfCells(); } mitk::FiberBundleX::Pointer mitk::FiberBundleX::ExtractFiberSubset(mitk::PlanarFigure::Pointer pf) { - return mitk::FiberBundleX::New(GeneratePolyDataByIds(ExtractFiberIdSubset(pf))); + std::vector tmp = ExtractFiberIdSubset(pf); + vtkSmartPointer pTmp = GeneratePolyDataByIds(tmp); + return mitk::FiberBundleX::New(pTmp); } std::vector mitk::FiberBundleX::ExtractFiberIdSubset(mitk::PlanarFigure::Pointer pf) { - MITK_DEBUG << "Extracting fibers!"; - // vector which is returned, contains all extracted FiberIds - std::vector FibersInROI; - - /* Handle type of planarfigure */ - // if incoming pf is a pfc - mitk::PlanarFigureComposite::Pointer pfcomp= dynamic_cast(pf.GetPointer()); - if (!pfcomp.IsNull()) { - // process requested boolean operation of PFC - switch (pfcomp->getOperationType()) { - case 0: - { - MITK_DEBUG << "AND PROCESSING"; - //AND - //temporarly store results of the child in this vector, we need that to accumulate the - std::vector childResults = this->ExtractFiberIdSubset(pfcomp->getChildAt(0)); - MITK_DEBUG << "first roi got fibers in ROI: " << childResults.size(); - MITK_DEBUG << "sorting..."; - std::sort(childResults.begin(), childResults.end()); - MITK_DEBUG << "sorting done"; - std::vector AND_Assamblage(childResults.size()); - //std::vector AND_Assamblage; - fill(AND_Assamblage.begin(), AND_Assamblage.end(), -1); - //AND_Assamblage.reserve(childResults.size()); //max size AND can reach anyway - - std::vector::iterator it; - for (int i=1; igetNumberOfChildren(); ++i) - { - std::vector tmpChild = this->ExtractFiberIdSubset(pfcomp->getChildAt(i)); - MITK_DEBUG << "ROI " << i << " has fibers in ROI: " << tmpChild.size(); - sort(tmpChild.begin(), tmpChild.end()); - - it = std::set_intersection(childResults.begin(), childResults.end(), - tmpChild.begin(), tmpChild.end(), - AND_Assamblage.begin() ); - } - - MITK_DEBUG << "resize Vector"; - long i=0; - while (i < AND_Assamblage.size() && AND_Assamblage[i] != -1){ //-1 represents a placeholder in the array - ++i; - } - AND_Assamblage.resize(i); - - MITK_DEBUG << "returning AND vector, size: " << AND_Assamblage.size(); - return AND_Assamblage; - // break; - - } - case 1: - { - //OR - std::vector OR_Assamblage = this->ExtractFiberIdSubset(pfcomp->getChildAt(0)); - std::vector::iterator it; - MITK_DEBUG << OR_Assamblage.size(); - - for (int i=1; igetNumberOfChildren(); ++i) { - it = OR_Assamblage.end(); - std::vector tmpChild = this->ExtractFiberIdSubset(pfcomp->getChildAt(i)); - OR_Assamblage.insert(it, tmpChild.begin(), tmpChild.end()); - MITK_DEBUG << "ROI " << i << " has fibers in ROI: " << tmpChild.size() << " OR Assamblage: " << OR_Assamblage.size(); - } - - sort(OR_Assamblage.begin(), OR_Assamblage.end()); - it = unique(OR_Assamblage.begin(), OR_Assamblage.end()); - OR_Assamblage.resize( it - OR_Assamblage.begin() ); - MITK_DEBUG << "returning OR vector, size: " << OR_Assamblage.size(); - - return OR_Assamblage; - } - case 2: - { - //NOT - //get IDs of all fibers - std::vector childResults; - childResults.reserve(this->GetNumFibers()); - vtkSmartPointer idSet = m_FiberIdDataSet->GetCellData()->GetArray(FIBER_ID_ARRAY); - MITK_DEBUG << "m_NumOfFib: " << this->GetNumFibers() << " cellIdNum: " << idSet->GetNumberOfTuples(); - for(long i=0; iGetNumFibers(); i++) - { - MITK_DEBUG << "i: " << i << " idset: " << idSet->GetTuple(i)[0]; - childResults.push_back(idSet->GetTuple(i)[0]); - } - - std::sort(childResults.begin(), childResults.end()); - std::vector NOT_Assamblage(childResults.size()); - //fill it with -1, otherwise 0 will be stored and 0 can also be an ID of fiber! - fill(NOT_Assamblage.begin(), NOT_Assamblage.end(), -1); - std::vector::iterator it; - - for (long i=0; igetNumberOfChildren(); ++i) - { - std::vector tmpChild = ExtractFiberIdSubset(pfcomp->getChildAt(i)); - sort(tmpChild.begin(), tmpChild.end()); - - it = std::set_difference(childResults.begin(), childResults.end(), - tmpChild.begin(), tmpChild.end(), - NOT_Assamblage.begin() ); - - } - - MITK_DEBUG << "resize Vector"; - long i=0; - while (NOT_Assamblage[i] != -1){ //-1 represents a placeholder in the array - ++i; - } - NOT_Assamblage.resize(i); - - return NOT_Assamblage; - } - default: - MITK_DEBUG << "we have an UNDEFINED composition... ERROR" ; - break; + MITK_DEBUG << "Extracting fibers!"; + // vector which is returned, contains all extracted FiberIds + std::vector FibersInROI; + + /* Handle type of planarfigure */ + // if incoming pf is a pfc + mitk::PlanarFigureComposite::Pointer pfcomp= dynamic_cast(pf.GetPointer()); + if (!pfcomp.IsNull()) { + // process requested boolean operation of PFC + switch (pfcomp->getOperationType()) { + case 0: + { + MITK_DEBUG << "AND PROCESSING"; + //AND + //temporarly store results of the child in this vector, we need that to accumulate the + std::vector childResults = this->ExtractFiberIdSubset(pfcomp->getChildAt(0)); + MITK_DEBUG << "first roi got fibers in ROI: " << childResults.size(); + MITK_DEBUG << "sorting..."; + std::sort(childResults.begin(), childResults.end()); + MITK_DEBUG << "sorting done"; + std::vector AND_Assamblage(childResults.size()); + //std::vector AND_Assamblage; + fill(AND_Assamblage.begin(), AND_Assamblage.end(), -1); + //AND_Assamblage.reserve(childResults.size()); //max size AND can reach anyway + + std::vector::iterator it; + for (int i=1; igetNumberOfChildren(); ++i) + { + std::vector tmpChild = this->ExtractFiberIdSubset(pfcomp->getChildAt(i)); + MITK_DEBUG << "ROI " << i << " has fibers in ROI: " << tmpChild.size(); + sort(tmpChild.begin(), tmpChild.end()); + + it = std::set_intersection(childResults.begin(), childResults.end(), + tmpChild.begin(), tmpChild.end(), + AND_Assamblage.begin() ); + } + + MITK_DEBUG << "resize Vector"; + long i=0; + while (i < AND_Assamblage.size() && AND_Assamblage[i] != -1){ //-1 represents a placeholder in the array + ++i; + } + AND_Assamblage.resize(i); + + MITK_DEBUG << "returning AND vector, size: " << AND_Assamblage.size(); + return AND_Assamblage; + // break; - } - } else + } + case 1: { + //OR + std::vector OR_Assamblage = this->ExtractFiberIdSubset(pfcomp->getChildAt(0)); + std::vector::iterator it; + MITK_DEBUG << OR_Assamblage.size(); + + for (int i=1; igetNumberOfChildren(); ++i) { + it = OR_Assamblage.end(); + std::vector tmpChild = this->ExtractFiberIdSubset(pfcomp->getChildAt(i)); + OR_Assamblage.insert(it, tmpChild.begin(), tmpChild.end()); + MITK_DEBUG << "ROI " << i << " has fibers in ROI: " << tmpChild.size() << " OR Assamblage: " << OR_Assamblage.size(); + } + + sort(OR_Assamblage.begin(), OR_Assamblage.end()); + it = unique(OR_Assamblage.begin(), OR_Assamblage.end()); + OR_Assamblage.resize( it - OR_Assamblage.begin() ); + MITK_DEBUG << "returning OR vector, size: " << OR_Assamblage.size(); + + return OR_Assamblage; + } + case 2: + { + //NOT + //get IDs of all fibers + std::vector childResults; + childResults.reserve(this->GetNumFibers()); + vtkSmartPointer idSet = m_FiberIdDataSet->GetCellData()->GetArray(FIBER_ID_ARRAY); + MITK_DEBUG << "m_NumOfFib: " << this->GetNumFibers() << " cellIdNum: " << idSet->GetNumberOfTuples(); + for(long i=0; iGetNumFibers(); i++) + { + MITK_DEBUG << "i: " << i << " idset: " << idSet->GetTuple(i)[0]; + childResults.push_back(idSet->GetTuple(i)[0]); + } + + std::sort(childResults.begin(), childResults.end()); + std::vector NOT_Assamblage(childResults.size()); + //fill it with -1, otherwise 0 will be stored and 0 can also be an ID of fiber! + fill(NOT_Assamblage.begin(), NOT_Assamblage.end(), -1); + std::vector::iterator it; + + for (long i=0; igetNumberOfChildren(); ++i) + { + std::vector tmpChild = ExtractFiberIdSubset(pfcomp->getChildAt(i)); + sort(tmpChild.begin(), tmpChild.end()); + + it = std::set_difference(childResults.begin(), childResults.end(), + tmpChild.begin(), tmpChild.end(), + NOT_Assamblage.begin() ); + + } + + MITK_DEBUG << "resize Vector"; + long i=0; + while (NOT_Assamblage[i] != -1){ //-1 represents a placeholder in the array + ++i; + } + NOT_Assamblage.resize(i); + + return NOT_Assamblage; + } + default: + MITK_DEBUG << "we have an UNDEFINED composition... ERROR" ; + break; - mitk::Geometry2D::ConstPointer pfgeometry = pf->GetGeometry2D(); - const mitk::PlaneGeometry* planeGeometry = dynamic_cast (pfgeometry.GetPointer()); - Vector3D planeNormal = planeGeometry->GetNormal(); - planeNormal.Normalize(); - Point3D planeOrigin = planeGeometry->GetOrigin(); - - MITK_DEBUG << "planeOrigin: " << planeOrigin[0] << " | " << planeOrigin[1] << " | " << planeOrigin[2] << endl; - MITK_DEBUG << "planeNormal: " << planeNormal[0] << " | " << planeNormal[1] << " | " << planeNormal[2] << endl; - - - /* init necessary vectors hosting pointIds and FiberIds */ - // contains all pointIds which are crossing the cutting plane - std::vector PointsOnPlane; - - // based on PointsOnPlane, all ROI relevant point IDs are stored here - std::vector PointsInROI; - - - /* Define cutting plane by ROI (PlanarFigure) */ - vtkSmartPointer plane = vtkSmartPointer::New(); - plane->SetOrigin(planeOrigin[0],planeOrigin[1],planeOrigin[2]); - plane->SetNormal(planeNormal[0],planeNormal[1],planeNormal[2]); - - //same plane but opposite normal direction. so point cloud will be reduced -> better performance - // vtkSmartPointer planeR = vtkSmartPointer::New(); - - //define new origin along the normal but close to the original one - // OriginNew = OriginOld + 1*Normal - // Vector3D extendedNormal; - // int multiplyFactor = 1; - // extendedNormal[0] = planeNormal[0] * multiplyFactor; - // extendedNormal[1] = planeNormal[1] * multiplyFactor; - // extendedNormal[2] = planeNormal[2] * multiplyFactor; - // Point3D RplaneOrigin = planeOrigin - extendedNormal; - // planeR->SetOrigin(RplaneOrigin[0],RplaneOrigin[1],RplaneOrigin[2]); - // planeR->SetNormal(-planeNormal[0],-planeNormal[1],-planeNormal[2]); - // MITK_DEBUG << "RPlaneOrigin: " << RplaneOrigin[0] << " | " << RplaneOrigin[1] - // << " | " << RplaneOrigin[2]; - - /* get all points/fibers cutting the plane */ - MITK_DEBUG << "start clipping"; - vtkSmartPointer clipper = vtkSmartPointer::New(); - clipper->SetInput(m_FiberIdDataSet); - clipper->SetClipFunction(plane); - clipper->GenerateClipScalarsOn(); - clipper->GenerateClippedOutputOn(); - vtkSmartPointer clipperout = clipper->GetClippedOutput(); - MITK_DEBUG << "end clipping"; - - /* for some reason clipperoutput is not initialized for futher processing + } + } + else + { + mitk::Geometry2D::ConstPointer pfgeometry = pf->GetGeometry2D(); + const mitk::PlaneGeometry* planeGeometry = dynamic_cast (pfgeometry.GetPointer()); + Vector3D planeNormal = planeGeometry->GetNormal(); + planeNormal.Normalize(); + Point3D planeOrigin = planeGeometry->GetOrigin(); + + MITK_DEBUG << "planeOrigin: " << planeOrigin[0] << " | " << planeOrigin[1] << " | " << planeOrigin[2] << endl; + MITK_DEBUG << "planeNormal: " << planeNormal[0] << " | " << planeNormal[1] << " | " << planeNormal[2] << endl; + + std::vector PointsOnPlane; // contains all pointIds which are crossing the cutting plane + std::vector PointsInROI; // based on PointsOnPlane, all ROI relevant point IDs are stored here + + /* Define cutting plane by ROI (PlanarFigure) */ + vtkSmartPointer plane = vtkSmartPointer::New(); + plane->SetOrigin(planeOrigin[0],planeOrigin[1],planeOrigin[2]); + plane->SetNormal(planeNormal[0],planeNormal[1],planeNormal[2]); + + //same plane but opposite normal direction. so point cloud will be reduced -> better performance + // vtkSmartPointer planeR = vtkSmartPointer::New(); + + //define new origin along the normal but close to the original one + // OriginNew = OriginOld + 1*Normal + // Vector3D extendedNormal; + // int multiplyFactor = 1; + // extendedNormal[0] = planeNormal[0] * multiplyFactor; + // extendedNormal[1] = planeNormal[1] * multiplyFactor; + // extendedNormal[2] = planeNormal[2] * multiplyFactor; + // Point3D RplaneOrigin = planeOrigin - extendedNormal; + // planeR->SetOrigin(RplaneOrigin[0],RplaneOrigin[1],RplaneOrigin[2]); + // planeR->SetNormal(-planeNormal[0],-planeNormal[1],-planeNormal[2]); + // MITK_DEBUG << "RPlaneOrigin: " << RplaneOrigin[0] << " | " << RplaneOrigin[1] + // << " | " << RplaneOrigin[2]; + + /* get all points/fibers cutting the plane */ + MITK_DEBUG << "start clipping"; + vtkSmartPointer clipper = vtkSmartPointer::New(); + clipper->SetInput(m_FiberIdDataSet); + clipper->SetClipFunction(plane); + clipper->GenerateClipScalarsOn(); + clipper->GenerateClippedOutputOn(); + vtkSmartPointer clipperout = clipper->GetClippedOutput(); + MITK_DEBUG << "end clipping"; + + /* for some reason clipperoutput is not initialized for futher processing * so far only writing out clipped polydata provides requested */ - // MITK_DEBUG << "writing clipper output"; - // vtkSmartPointer writerC = vtkSmartPointer::New(); - // writerC->SetInput(clipperout1); - // writerC->SetFileName("/vtkOutput/Clipping.vtk"); - // writerC->SetFileTypeToASCII(); - // writerC->Write(); - // MITK_DEBUG << "writing done"; - - MITK_DEBUG << "init and update clipperoutput"; - clipperout->GetPointData()->Initialize(); - clipperout->Update(); - MITK_DEBUG << "init and update clipperoutput completed"; - - // MITK_DEBUG << "start clippingRecursive"; - // vtkSmartPointer Rclipper = vtkSmartPointer::New(); - // Rclipper->SetInput(clipperout1); - // Rclipper->SetClipFunction(planeR); - // Rclipper->GenerateClipScalarsOn(); - // Rclipper->GenerateClippedOutputOn(); - // vtkSmartPointer clipperout = Rclipper->GetClippedOutput(); - // MITK_DEBUG << "end clipping recursive"; - - // MITK_DEBUG << "writing clipper output 2"; - // vtkSmartPointer writerC1 = vtkSmartPointer::New(); - // writerC1->SetInput(clipperout); - // writerC1->SetFileName("/vtkOutput/RClipping.vtk"); - // writerC1->SetFileTypeToASCII(); - // writerC1->Write(); - // MITK_DEBUG << "init and update clipperoutput"; - // clipperout->GetPointData()->Initialize(); - // clipperout->Update(); - // MITK_DEBUG << "init and update clipperoutput completed"; - - MITK_DEBUG << "STEP 1: find all points which have distance 0 to the given plane"; - /*======STEP 1====== + // MITK_DEBUG << "writing clipper output"; + // vtkSmartPointer writerC = vtkSmartPointer::New(); + // writerC->SetInput(clipperout1); + // writerC->SetFileName("/vtkOutput/Clipping.vtk"); + // writerC->SetFileTypeToASCII(); + // writerC->Write(); + // MITK_DEBUG << "writing done"; + + MITK_DEBUG << "init and update clipperoutput"; + clipperout->GetPointData()->Initialize(); + clipperout->Update(); + MITK_DEBUG << "init and update clipperoutput completed"; + + // MITK_DEBUG << "start clippingRecursive"; + // vtkSmartPointer Rclipper = vtkSmartPointer::New(); + // Rclipper->SetInput(clipperout1); + // Rclipper->SetClipFunction(planeR); + // Rclipper->GenerateClipScalarsOn(); + // Rclipper->GenerateClippedOutputOn(); + // vtkSmartPointer clipperout = Rclipper->GetClippedOutput(); + // MITK_DEBUG << "end clipping recursive"; + + // MITK_DEBUG << "writing clipper output 2"; + // vtkSmartPointer writerC1 = vtkSmartPointer::New(); + // writerC1->SetInput(clipperout); + // writerC1->SetFileName("/vtkOutput/RClipping.vtk"); + // writerC1->SetFileTypeToASCII(); + // writerC1->Write(); + // MITK_DEBUG << "init and update clipperoutput"; + // clipperout->GetPointData()->Initialize(); + // clipperout->Update(); + // MITK_DEBUG << "init and update clipperoutput completed"; + + MITK_DEBUG << "STEP 1: find all points which have distance 0 to the given plane"; + /*======STEP 1====== * extract all points, which are crossing the plane */ - // Scalar values describe the distance between each remaining point to the given plane. Values sorted by point index - vtkSmartPointer distanceList = clipperout->GetPointData()->GetScalars(); - vtkIdType sizeOfList = distanceList->GetNumberOfTuples(); - PointsOnPlane.reserve(sizeOfList); /* use reserve for high-performant push_back, no hidden copy procedures are processed then! + // Scalar values describe the distance between each remaining point to the given plane. Values sorted by point index + vtkSmartPointer distanceList = clipperout->GetPointData()->GetScalars(); + vtkIdType sizeOfList = distanceList->GetNumberOfTuples(); + PointsOnPlane.reserve(sizeOfList); /* use reserve for high-performant push_back, no hidden copy procedures are processed then! * size of list can be optimized by reducing allocation, but be aware of iterator and vector size*/ - for (int i=0; iGetTuple(i); - //std::cout << "distance of point " << i << " : " << distance[0] << std::endl; + for (int i=0; iGetTuple(i); - // check if point is on plane. - // 0.01 due to some approximation errors when calculating distance - if (distance[0] >= -0.01 && distance[0] <= 0.01) - { - //std::cout << "adding " << i << endl; - PointsOnPlane.push_back(i); //push back in combination with reserve is fastest way to fill vector with various values - - } - - } + // check if point is on plane. + // 0.01 due to some approximation errors when calculating distance + if (distance[0] >= -0.01 && distance[0] <= 0.01) + PointsOnPlane.push_back(i); + } - // DEBUG print out all interesting points, stop where array starts with value -1. after -1 no more interesting idx are set! - // std::vector::iterator rit = PointsOnPlane.begin(); - // while (rit != PointsOnPlane.end() ) { - // std::cout << "interesting point: " << *rit << " coord: " << clipperout->GetPoint(*rit)[0] << " | " << clipperout->GetPoint(*rit)[1] << " | " << clipperout->GetPoint(*rit)[2] << endl; - // rit++; - // } + // DEBUG print out all interesting points, stop where array starts with value -1. after -1 no more interesting idx are set! + // std::vector::iterator rit = PointsOnPlane.begin(); + // while (rit != PointsOnPlane.end() ) { + // std::cout << "interesting point: " << *rit << " coord: " << clipperout->GetPoint(*rit)[0] << " | " << clipperout->GetPoint(*rit)[1] << " | " << clipperout->GetPoint(*rit)[2] << endl; + // rit++; + // } - MITK_DEBUG << "Num Of points on plane: " << PointsOnPlane.size(); + MITK_DEBUG << "Num Of points on plane: " << PointsOnPlane.size(); - MITK_DEBUG << "Step 2: extract Interesting points with respect to given extraction planarFigure"; + MITK_DEBUG << "Step 2: extract Interesting points with respect to given extraction planarFigure"; - PointsInROI.reserve(PointsOnPlane.size()); - /*=======STEP 2===== + PointsInROI.reserve(PointsOnPlane.size()); + /*=======STEP 2===== * extract ROI relevant pointIds */ - mitk::PlanarCircle::Pointer circleName = mitk::PlanarCircle::New(); - mitk::PlanarPolygon::Pointer polyName = mitk::PlanarPolygon::New(); - if (pf->GetNameOfClass() == circleName->GetNameOfClass() ) - { - - //calculate circle radius - mitk::Point3D V1w = pf->GetWorldControlPoint(0); //centerPoint - mitk::Point3D V2w = pf->GetWorldControlPoint(1); //radiusPoint - - //calculate distance between those 2 and - double distPF; - distPF = sqrt((double) (V2w[0] - V1w[0]) * (V2w[0] - V1w[0]) + - (V2w[1] - V1w[1]) * (V2w[1] - V1w[1]) + - (V2w[2] - V1w[2]) * (V2w[2] - V1w[2])); - - //MITK_DEBUG << "Circle Radius: " << distPF; - - for (int i=0; iGetPoint(PointsOnPlane[i])[0] << " - " << V1w[0]; - // MITK_DEBUG << clipperout->GetPoint(PointsOnPlane[i])[1] << " - " << V1w[1]; - // MITK_DEBUG << clipperout->GetPoint(PointsOnPlane[i])[2] << " - " << V1w[2]; - - //distance between circle radius and given point - double XdistPnt = sqrt((double) (clipperout->GetPoint(PointsOnPlane[i])[0] - V1w[0]) * (clipperout->GetPoint(PointsOnPlane[i])[0] - V1w[0]) + - (clipperout->GetPoint(PointsOnPlane[i])[1] - V1w[1]) * (clipperout->GetPoint(PointsOnPlane[i])[1] - V1w[1]) + - (clipperout->GetPoint(PointsOnPlane[i])[2] - V1w[2]) * (clipperout->GetPoint(PointsOnPlane[i])[2] - V1w[2])) ; - - // MITK_DEBUG << "PntDistance to Radius: " << XdistPnt; - if( XdistPnt <= distPF) - { - // MITK_DEBUG << "point in Circle"; - PointsInROI.push_back(PointsOnPlane[i]); - } - - }//end for(i) - // MITK_DEBUG << "Points inside circle radius: " << PointsInROI.size(); - } + mitk::PlanarCircle::Pointer circleName = mitk::PlanarCircle::New(); + mitk::PlanarPolygon::Pointer polyName = mitk::PlanarPolygon::New(); + if (pf->GetNameOfClass() == circleName->GetNameOfClass() ) + { + //calculate circle radius + mitk::Point3D V1w = pf->GetWorldControlPoint(0); //centerPoint + mitk::Point3D V2w = pf->GetWorldControlPoint(1); //radiusPoint + + double distPF = V1w.EuclideanDistanceTo(V2w); +// distPF = sqrt((double) (V2w[0] - V1w[0]) * (V2w[0] - V1w[0]) + +// (V2w[1] - V1w[1]) * (V2w[1] - V1w[1]) + +// (V2w[2] - V1w[2]) * (V2w[2] - V1w[2])); + + for (int i=0; iGetPoint(PointsOnPlane[i])[0] - V1w[0]) * (clipperout->GetPoint(PointsOnPlane[i])[0] - V1w[0]) + + (clipperout->GetPoint(PointsOnPlane[i])[1] - V1w[1]) * (clipperout->GetPoint(PointsOnPlane[i])[1] - V1w[1]) + + (clipperout->GetPoint(PointsOnPlane[i])[2] - V1w[2]) * (clipperout->GetPoint(PointsOnPlane[i])[2] - V1w[2])) ; + + // MITK_DEBUG << "PntDistance to Radius: " << XdistPnt; + if( XdistPnt <= distPF) + PointsInROI.push_back(PointsOnPlane[i]); + } + // MITK_DEBUG << "Points inside circle radius: " << PointsInROI.size(); + } - MITK_DEBUG << "Step3: Identify fibers"; + MITK_DEBUG << "Step3: Identify fibers"; - /*======STEP 3======= + /*======STEP 3======= * identify fiberIds for points in ROI */ - // we need to access the fiberId Array, so make sure that this array is available - if (!clipperout->GetCellData()->HasArray(FIBER_ID_ARRAY)) - { - MITK_DEBUG << "ERROR: FiberID array does not exist, no correlation between points and fiberIds possible! Make sure calling GenerateFiberIds()"; - return FibersInROI; // FibersInRoi is empty then - } - - // prepare a structure where each point id is represented as an indexId. - // vector looks like: | pntId | fiberIdx | - std::vector< long > pointindexFiberMap; - - // walk through the whole subline section and create an vector sorted by point index - vtkCellArray *clipperlines = clipperout->GetLines(); - clipperlines->InitTraversal(); - long numOfLineCells = clipperlines->GetNumberOfCells(); - long numofClippedPoints = clipperout->GetNumberOfPoints(); - pointindexFiberMap.reserve(numofClippedPoints); + // we need to access the fiberId Array, so make sure that this array is available + if (!clipperout->GetCellData()->HasArray(FIBER_ID_ARRAY)) + { + MITK_DEBUG << "ERROR: FiberID array does not exist, no correlation between points and fiberIds possible! Make sure calling GenerateFiberIds()"; + return FibersInROI; // FibersInRoi is empty then + } + // prepare a structure where each point id is represented as an indexId. + // vector looks like: | pntId | fiberIdx | + std::vector< long > pointindexFiberMap; - //prepare resulting vector - FibersInROI.reserve(PointsInROI.size()); + // walk through the whole subline section and create an vector sorted by point index + vtkCellArray *clipperlines = clipperout->GetLines(); + clipperlines->InitTraversal(); + long numOfLineCells = clipperlines->GetNumberOfCells(); + long numofClippedPoints = clipperout->GetNumberOfPoints(); + pointindexFiberMap.reserve(numofClippedPoints); - MITK_DEBUG << "\n===== Pointindex based structure initialized ======\n"; - // go through resulting "sub"lines which are stored as cells, "i" corresponds to current line id. - for (int i=0, ic=0 ; iGetCell(ic, npts, pts); + MITK_DEBUG << "\n===== Pointindex based structure initialized ======\n"; - // go through point ids in hosting subline, "j" corresponds to current pointindex in current line i. eg. idx[0]=45; idx[1]=46 - for (long j=0; jGetCellData()->GetArray(FIBER_ID_ARRAY)->GetTuple(i)[0] << " to pointId: " << pts[j]; - pointindexFiberMap[ pts[j] ] = clipperout->GetCellData()->GetArray(FIBER_ID_ARRAY)->GetTuple(i)[0]; - // MITK_DEBUG << "in array: " << pointindexFiberMap[ pts[j] ]; - } + // go through resulting "sub"lines which are stored as cells, "i" corresponds to current line id. + for (int i=0, ic=0 ; iGetCell(ic, npts, pts); - MITK_DEBUG << "\n===== Pointindex based structure finalized ======\n"; - - // get all Points in ROI with according fiberID - for (long k = 0; k < PointsInROI.size(); k++) - { - //MITK_DEBUG << "point " << PointsInROI[k] << " belongs to fiber " << pointindexFiberMap[ PointsInROI[k] ]; - FibersInROI.push_back(pointindexFiberMap[ PointsInROI[k] ]); - } + // go through point ids in hosting subline, "j" corresponds to current pointindex in current line i. eg. idx[0]=45; idx[1]=46 + for (long j=0; jGetCellData()->GetArray(FIBER_ID_ARRAY)->GetTuple(i)[0] << " to pointId: " << pts[j]; + pointindexFiberMap[ pts[j] ] = clipperout->GetCellData()->GetArray(FIBER_ID_ARRAY)->GetTuple(i)[0]; + // MITK_DEBUG << "in array: " << pointindexFiberMap[ pts[j] ]; + } } - // detecting fiberId duplicates - MITK_DEBUG << "check for duplicates"; + MITK_DEBUG << "\n===== Pointindex based structure finalized ======\n"; - sort(FibersInROI.begin(), FibersInROI.end()); - bool hasDuplicats = false; - for(long i=0; i=0) + FibersInROI.push_back(pointindexFiberMap[ PointsInROI[k] ]); + else + MITK_INFO << "ERROR in ExtractFiberIdSubset; impossible fiber id detected"; } - if(hasDuplicats) - { - std::vector::iterator it; - it = unique (FibersInROI.begin(), FibersInROI.end()); - FibersInROI.resize( it - FibersInROI.begin() ); - } + } + + // detecting fiberId duplicates + MITK_DEBUG << "check for duplicates"; - return FibersInROI; + sort(FibersInROI.begin(), FibersInROI.end()); + bool hasDuplicats = false; + for(long i=0; i::iterator it; + it = unique (FibersInROI.begin(), FibersInROI.end()); + FibersInROI.resize( it - FibersInROI.begin() ); + } + + return FibersInROI; } void mitk::FiberBundleX::UpdateFiberGeometry() { - if (m_NumFibers<=0) - { - mitk::Geometry3D::Pointer geometry = mitk::Geometry3D::New(); - geometry->SetImageGeometry(true); - float b[] = {0, 1, 0, 1, 0, 1}; - geometry->SetFloatBounds(b); - SetGeometry(geometry); - return; - } - float min = itk::NumericTraits::min(); - float max = itk::NumericTraits::max(); - float b[] = {max, min, max, min, max, min}; + if (m_NumFibers<=0) + { + mitk::Geometry3D::Pointer geometry = mitk::Geometry3D::New(); + geometry->SetImageGeometry(true); + float b[] = {0, 1, 0, 1, 0, 1}; + geometry->SetFloatBounds(b); + SetGeometry(geometry); + return; + } + float min = itk::NumericTraits::min(); + float max = itk::NumericTraits::max(); + float b[] = {max, min, max, min, max, min}; - vtkCellArray* cells = m_FiberPolyData->GetLines(); - cells->InitTraversal(); - for (int i=0; iGetNumberOfCells(); i++) + vtkCellArray* cells = m_FiberPolyData->GetLines(); + cells->InitTraversal(); + for (int i=0; iGetNumberOfCells(); i++) + { + vtkCell* cell = m_FiberPolyData->GetCell(i); + int p = cell->GetNumberOfPoints(); + vtkPoints* points = cell->GetPoints(); + for (int j=0; jGetCell(i); - int p = cell->GetNumberOfPoints(); - vtkPoints* points = cell->GetPoints(); - for (int j=0; jGetPoint(j, p); - - if (p[0]b[1]) - b[1]=p[0]; - - if (p[1]b[3]) - b[3]=p[1]; - - if (p[2]b[5]) - b[5]=p[2]; - } + double p[3]; + points->GetPoint(j, p); + + if (p[0]b[1]) + b[1]=p[0]; + + if (p[1]b[3]) + b[3]=p[1]; + + if (p[2]b[5]) + b[5]=p[2]; } + } - // provide some border margin - for(int i=0; i<=4; i+=2) - b[i] -=10; - for(int i=1; i<=5; i+=2) - b[i] +=10; + // provide some border margin + for(int i=0; i<=4; i+=2) + b[i] -=10; + for(int i=1; i<=5; i+=2) + b[i] +=10; - mitk::Geometry3D::Pointer geometry = mitk::Geometry3D::New(); - geometry->SetImageGeometry(true); - geometry->SetFloatBounds(b); - this->SetGeometry(geometry); + mitk::Geometry3D::Pointer geometry = mitk::Geometry3D::New(); + geometry->SetImageGeometry(true); + geometry->SetFloatBounds(b); + this->SetGeometry(geometry); } QStringList mitk::FiberBundleX::GetAvailableColorCodings() { - QStringList availableColorCodings; - int numColors = m_FiberPolyData->GetPointData()->GetNumberOfArrays(); - for(int i=0; iGetPointData()->GetArrayName(i)); - } + QStringList availableColorCodings; + int numColors = m_FiberPolyData->GetPointData()->GetNumberOfArrays(); + for(int i=0; iGetPointData()->GetArrayName(i)); + } - //this controlstructure shall be implemented by the calling method - if (availableColorCodings.isEmpty()) - MITK_DEBUG << "no colorcodings available in fiberbundleX"; + //this controlstructure shall be implemented by the calling method + if (availableColorCodings.isEmpty()) + MITK_DEBUG << "no colorcodings available in fiberbundleX"; - // for(int i=0; im_currentColorCoding = (char*) COLORCODING_ORIENTATION_BASED; + if( strcmp (COLORCODING_ORIENTATION_BASED,requestedColorCoding) == 0 ) { + this->m_currentColorCoding = (char*) COLORCODING_ORIENTATION_BASED; - } else if( strcmp (COLORCODING_FA_BASED,requestedColorCoding) == 0 ) { - this->m_currentColorCoding = (char*) COLORCODING_FA_BASED; + } else if( strcmp (COLORCODING_FA_BASED,requestedColorCoding) == 0 ) { + this->m_currentColorCoding = (char*) COLORCODING_FA_BASED; - } else if( strcmp (COLORCODING_CUSTOM,requestedColorCoding) == 0 ) { - this->m_currentColorCoding = (char*) COLORCODING_CUSTOM; + } else if( strcmp (COLORCODING_CUSTOM,requestedColorCoding) == 0 ) { + this->m_currentColorCoding = (char*) COLORCODING_CUSTOM; - } else { - MITK_DEBUG << "FIBERBUNDLE X: UNKNOWN COLORCODING in FIBERBUNDLEX Datastructure"; - this->m_currentColorCoding = (char*) COLORCODING_CUSTOM; //will cause blank colorcoding of fibers - } + } else { + MITK_DEBUG << "FIBERBUNDLE X: UNKNOWN COLORCODING in FIBERBUNDLEX Datastructure"; + this->m_currentColorCoding = (char*) COLORCODING_CUSTOM; //will cause blank colorcoding of fibers + } } void mitk::FiberBundleX::DoFiberSmoothing(int pointsPerCm) { vtkSmartPointer vtkSmoothPoints = vtkPoints::New(); //in smoothpoints the interpolated points representing a fiber are stored. //in vtkcells all polylines are stored, actually all id's of them are stored vtkSmartPointer vtkSmoothCells = vtkCellArray::New(); //cellcontainer for smoothed lines vtkSmartPointer vLines = m_FiberPolyData->GetLines(); vLines->InitTraversal(); vtkIdType pointHelperCnt = 0; for (int i=0; iGetNextCell ( numPoints, pointIds ); vtkSmartPointer points = vtkSmartPointer::New(); float length = 0; itk::Point lastP; for (int j=0; jGetPoint(pointIds[j]); points->InsertNextPoint(p); if (j>0) length += sqrt(pow(p[0]-lastP[0], 2)+pow(p[1]-lastP[1], 2)+pow(p[2]-lastP[2], 2)); lastP[0] = p[0]; lastP[1] = p[1]; lastP[2] = p[2]; } length /=10; int sampling = pointsPerCm*length; /////PROCESS POLYLINE SMOOTHING///// vtkSmartPointer xSpline = vtkKochanekSpline::New(); vtkSmartPointer ySpline = vtkKochanekSpline::New(); vtkSmartPointer zSpline = vtkKochanekSpline::New(); vtkSmartPointer spline = vtkParametricSpline::New(); spline->SetXSpline(xSpline); spline->SetYSpline(ySpline); spline->SetZSpline(zSpline); spline->SetPoints(points); vtkSmartPointer functionSource = vtkParametricFunctionSource::New(); functionSource->SetParametricFunction(spline); functionSource->SetUResolution(sampling); functionSource->SetVResolution(sampling); functionSource->SetWResolution(sampling); functionSource->Update(); vtkPolyData* outputFunction = functionSource->GetOutput(); vtkPoints* tmpSmoothPnts = outputFunction->GetPoints(); //smoothPoints of current fiber vtkSmartPointer smoothLine = vtkPolyLine::New(); smoothLine->GetPointIds()->SetNumberOfIds(tmpSmoothPnts->GetNumberOfPoints()); for (int j=0; jGetNumberOfPoints(); j++) { smoothLine->GetPointIds()->SetId(j, j+pointHelperCnt); vtkSmoothPoints->InsertNextPoint(tmpSmoothPnts->GetPoint(j)); } vtkSmoothCells->InsertNextCell(smoothLine); pointHelperCnt += tmpSmoothPnts->GetNumberOfPoints(); } m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkSmoothPoints); m_FiberPolyData->SetLines(vtkSmoothCells); UpdateColorCoding(); UpdateFiberGeometry(); } void mitk::FiberBundleX::ResampleFibers() { - mitk::Geometry3D::Pointer geometry = GetGeometry(); - mitk::Vector3D spacing = geometry->GetSpacing(); - - float minSpacing = 1; - if(spacing[0]GetSpacing(); + + float minSpacing = 1; + if(spacing[0] newPoly = vtkSmartPointer::New(); - vtkSmartPointer newCellArray = vtkSmartPointer::New(); - vtkSmartPointer newPoints = vtkSmartPointer::New(); - - vtkSmartPointer vLines = m_FiberPolyData->GetLines(); - vLines->InitTraversal(); - int numberOfLines = m_NumFibers; + vtkSmartPointer newPoly = vtkSmartPointer::New(); + vtkSmartPointer newCellArray = vtkSmartPointer::New(); + vtkSmartPointer newPoints = vtkSmartPointer::New(); - for (int i=0; iGetNextCell ( numPoints, points ); + vtkSmartPointer vLines = m_FiberPolyData->GetLines(); + vLines->InitTraversal(); + int numberOfLines = m_NumFibers; - vtkSmartPointer container = vtkSmartPointer::New(); + for (int i=0; iGetNextCell ( numPoints, points ); - double* point = m_FiberPolyData->GetPoint(points[0]); - vtkIdType pointId = newPoints->InsertNextPoint(point); - container->GetPointIds()->InsertNextId(pointId); + vtkSmartPointer container = vtkSmartPointer::New(); - float dtau = 0; - int cur_p = 1; - itk::Vector dR; - float normdR = 0; + double* point = m_FiberPolyData->GetPoint(points[0]); + vtkIdType pointId = newPoints->InsertNextPoint(point); + container->GetPointIds()->InsertNextId(pointId); - for (;;) - { - while (dtau <= len && cur_p < numPoints) - { - itk::Vector v1; - point = m_FiberPolyData->GetPoint(points[cur_p-1]); - v1[0] = point[0]; - v1[1] = point[1]; - v1[2] = point[2]; - itk::Vector v2; - point = m_FiberPolyData->GetPoint(points[cur_p]); - v2[0] = point[0]; - v2[1] = point[1]; - v2[2] = point[2]; - - dR = v2 - v1; - normdR = std::sqrt(dR.GetSquaredNorm()); - dtau += normdR; - cur_p++; - } - - if (dtau >= len) - { - itk::Vector v1; - point = m_FiberPolyData->GetPoint(points[cur_p-1]); - v1[0] = point[0]; - v1[1] = point[1]; - v1[2] = point[2]; - - itk::Vector v2 = v1 - dR*( (dtau-len)/normdR ); - pointId = newPoints->InsertNextPoint(v2.GetDataPointer()); - container->GetPointIds()->InsertNextId(pointId); - } - else - { - point = m_FiberPolyData->GetPoint(points[numPoints-1]); - pointId = newPoints->InsertNextPoint(point); - container->GetPointIds()->InsertNextId(pointId); - break; - } - dtau = dtau-len; - } + float dtau = 0; + int cur_p = 1; + itk::Vector dR; + float normdR = 0; - newCellArray->InsertNextCell(container); + for (;;) + { + while (dtau <= len && cur_p < numPoints) + { + itk::Vector v1; + point = m_FiberPolyData->GetPoint(points[cur_p-1]); + v1[0] = point[0]; + v1[1] = point[1]; + v1[2] = point[2]; + itk::Vector v2; + point = m_FiberPolyData->GetPoint(points[cur_p]); + v2[0] = point[0]; + v2[1] = point[1]; + v2[2] = point[2]; + + dR = v2 - v1; + normdR = std::sqrt(dR.GetSquaredNorm()); + dtau += normdR; + cur_p++; + } + + if (dtau >= len) + { + itk::Vector v1; + point = m_FiberPolyData->GetPoint(points[cur_p-1]); + v1[0] = point[0]; + v1[1] = point[1]; + v1[2] = point[2]; + + itk::Vector v2 = v1 - dR*( (dtau-len)/normdR ); + pointId = newPoints->InsertNextPoint(v2.GetDataPointer()); + container->GetPointIds()->InsertNextId(pointId); + } + else + { + point = m_FiberPolyData->GetPoint(points[numPoints-1]); + pointId = newPoints->InsertNextPoint(point); + container->GetPointIds()->InsertNextId(pointId); + break; + } + dtau = dtau-len; } - newPoly->SetPoints(newPoints); - newPoly->SetLines(newCellArray); - m_FiberPolyData = newPoly; - UpdateFiberGeometry(); - UpdateColorCoding(); + newCellArray->InsertNextCell(container); + } + + newPoly->SetPoints(newPoints); + newPoly->SetLines(newCellArray); + m_FiberPolyData = newPoly; + UpdateFiberGeometry(); + UpdateColorCoding(); } // reapply selected colorcoding in case polydata structure has changed void mitk::FiberBundleX::UpdateColorCoding() { char* cc = GetCurrentColorCoding(); if( strcmp (COLORCODING_ORIENTATION_BASED,cc) == 0 ) - DoColorCodingOrientationbased(); + DoColorCodingOrientationbased(); else if( strcmp (COLORCODING_FA_BASED,cc) == 0 ) DoColorCodingFAbased(); } /* ESSENTIAL IMPLEMENTATION OF SUPERCLASS METHODS */ void mitk::FiberBundleX::UpdateOutputInformation() { } void mitk::FiberBundleX::SetRequestedRegionToLargestPossibleRegion() { } bool mitk::FiberBundleX::RequestedRegionIsOutsideOfTheBufferedRegion() { - return false; + return false; } bool mitk::FiberBundleX::VerifyRequestedRegion() { - return true; + return true; } void mitk::FiberBundleX::SetRequestedRegion( itk::DataObject *data ) { } diff --git a/Modules/DiffusionImaging/IODataStructures/FiberBundleX/mitkFiberBundleXReader.cpp b/Modules/DiffusionImaging/IODataStructures/FiberBundleX/mitkFiberBundleXReader.cpp index e01189ee76..7b4000ebd8 100644 --- a/Modules/DiffusionImaging/IODataStructures/FiberBundleX/mitkFiberBundleXReader.cpp +++ b/Modules/DiffusionImaging/IODataStructures/FiberBundleX/mitkFiberBundleXReader.cpp @@ -1,282 +1,282 @@ /*========================================================================= Program: Medical Imaging & Interaction Toolkit Language: C++ Date: $Date: 2009-07-14 19:11:20 +0200 (Tue, 14 Jul 2009) $ Version: $Revision: 18127 $ 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 "mitkFiberBundleXReader.h" #include #include #include #include #include #include #include #include #include #include namespace mitk { void FiberBundleXReader ::GenerateData() { if ( ( ! m_OutputCache ) ) { Superclass::SetNumberOfRequiredOutputs(0); this->GenerateOutputInformation(); } if (!m_OutputCache) { itkWarningMacro("Output cache is empty!"); } Superclass::SetNumberOfRequiredOutputs(1); Superclass::SetNthOutput(0, m_OutputCache.GetPointer()); } void FiberBundleXReader::GenerateOutputInformation() { try { const std::string& locale = "C"; const std::string& currLocale = setlocale( LC_ALL, NULL ); setlocale(LC_ALL, locale.c_str()); std::string ext = itksys::SystemTools::GetFilenameLastExtension(m_FileName); ext = itksys::SystemTools::LowerCase(ext); vtkSmartPointer chooser=vtkDataReader::New(); chooser->SetFileName(m_FileName.c_str() ); if( chooser->IsFilePolyData()) { MITK_INFO << "Reading vtk fiber bundle"; vtkSmartPointer reader = vtkPolyDataReader::New(); reader->SetFileName( m_FileName.c_str() ); reader->Update(); if ( reader->GetOutput() != NULL ) { vtkSmartPointer fiberPolyData = reader->GetOutput(); - vtkSmartPointer cleaner = vtkSmartPointer::New(); - cleaner->SetInput(fiberPolyData); - cleaner->Update(); - fiberPolyData = cleaner->GetOutput(); +// vtkSmartPointer cleaner = vtkSmartPointer::New(); +// cleaner->SetInput(fiberPolyData); +// cleaner->Update(); +// fiberPolyData = cleaner->GetOutput(); m_OutputCache = OutputType::New(fiberPolyData); } } else // try to read deprecated fiber bundle file format { MITK_INFO << "Reading xml fiber bundle"; vtkSmartPointer fiberPolyData = vtkSmartPointer::New(); vtkSmartPointer cellArray = vtkSmartPointer::New(); vtkSmartPointer points = vtkSmartPointer::New(); TiXmlDocument doc( m_FileName ); if(doc.LoadFile()) { TiXmlHandle hDoc(&doc); TiXmlElement* pElem; TiXmlHandle hRoot(0); pElem = hDoc.FirstChildElement().Element(); // save this for later hRoot = TiXmlHandle(pElem); pElem = hRoot.FirstChildElement("geometry").Element(); // read geometry mitk::Geometry3D::Pointer geometry = mitk::Geometry3D::New(); // read origin mitk::Point3D origin; double temp = 0; pElem->Attribute("origin_x", &temp); origin[0] = temp; pElem->Attribute("origin_y", &temp); origin[1] = temp; pElem->Attribute("origin_z", &temp); origin[2] = temp; geometry->SetOrigin(origin); // read spacing float spacing[3]; pElem->Attribute("spacing_x", &temp); spacing[0] = temp; pElem->Attribute("spacing_y", &temp); spacing[1] = temp; pElem->Attribute("spacing_z", &temp); spacing[2] = temp; geometry->SetSpacing(spacing); // read transform vtkMatrix4x4* m = vtkMatrix4x4::New(); pElem->Attribute("xx", &temp); m->SetElement(0,0,temp); pElem->Attribute("xy", &temp); m->SetElement(1,0,temp); pElem->Attribute("xz", &temp); m->SetElement(2,0,temp); pElem->Attribute("yx", &temp); m->SetElement(0,1,temp); pElem->Attribute("yy", &temp); m->SetElement(1,1,temp); pElem->Attribute("yz", &temp); m->SetElement(2,1,temp); pElem->Attribute("zx", &temp); m->SetElement(0,2,temp); pElem->Attribute("zy", &temp); m->SetElement(1,2,temp); pElem->Attribute("zz", &temp); m->SetElement(2,2,temp); m->SetElement(0,3,origin[0]); m->SetElement(1,3,origin[1]); m->SetElement(2,3,origin[2]); m->SetElement(3,3,1); geometry->SetIndexToWorldTransformByVtkMatrix(m); // read bounds float bounds[] = {0, 0, 0, 0, 0, 0}; pElem->Attribute("size_x", &temp); bounds[1] = temp; pElem->Attribute("size_y", &temp); bounds[3] = temp; pElem->Attribute("size_z", &temp); bounds[5] = temp; geometry->SetFloatBounds(bounds); geometry->SetImageGeometry(true); pElem = hRoot.FirstChildElement("fiber_bundle").FirstChild().Element(); for( pElem; pElem; pElem=pElem->NextSiblingElement()) { TiXmlElement* pElem2 = pElem->FirstChildElement(); vtkSmartPointer container = vtkSmartPointer::New(); for( pElem2; pElem2; pElem2=pElem2->NextSiblingElement()) { itk::Point point; pElem2->Attribute("pos_x", &temp); point[0] = temp; pElem2->Attribute("pos_y", &temp); point[1] = temp; pElem2->Attribute("pos_z", &temp); point[2] = temp; geometry->IndexToWorld(point, point); vtkIdType id = points->InsertNextPoint(point.GetDataPointer()); container->GetPointIds()->InsertNextId(id); } cellArray->InsertNextCell(container); } fiberPolyData->SetPoints(points); fiberPolyData->SetLines(cellArray); vtkSmartPointer cleaner = vtkSmartPointer::New(); cleaner->SetInput(fiberPolyData); cleaner->Update(); fiberPolyData = cleaner->GetOutput(); m_OutputCache = OutputType::New(fiberPolyData); } else { MITK_INFO << "could not open xml file"; throw "could not open xml file"; } } setlocale(LC_ALL, currLocale.c_str()); MITK_INFO << "Fiber bundle read"; } catch(...) { throw; } } void FiberBundleXReader::Update() { this->GenerateData(); } const char* FiberBundleXReader ::GetFileName() const { return m_FileName.c_str(); } void FiberBundleXReader ::SetFileName(const char* aFileName) { m_FileName = aFileName; } const char* FiberBundleXReader ::GetFilePrefix() const { return m_FilePrefix.c_str(); } void FiberBundleXReader ::SetFilePrefix(const char* aFilePrefix) { m_FilePrefix = aFilePrefix; } const char* FiberBundleXReader ::GetFilePattern() const { return m_FilePattern.c_str(); } void FiberBundleXReader ::SetFilePattern(const char* aFilePattern) { m_FilePattern = aFilePattern; } bool FiberBundleXReader ::CanReadFile(const std::string filename, const std::string /*filePrefix*/, const std::string /*filePattern*/) { // First check the extension if( filename == "" ) { return false; } std::string ext = itksys::SystemTools::GetFilenameLastExtension(filename); ext = itksys::SystemTools::LowerCase(ext); if (ext == ".fib") { return true; } return false; } } //namespace MITK diff --git a/Modules/DiffusionImaging/IODataStructures/FiberBundleX/mitkFiberBundleXWriter.cpp b/Modules/DiffusionImaging/IODataStructures/FiberBundleX/mitkFiberBundleXWriter.cpp index 9856581eaa..786531bd83 100644 --- a/Modules/DiffusionImaging/IODataStructures/FiberBundleX/mitkFiberBundleXWriter.cpp +++ b/Modules/DiffusionImaging/IODataStructures/FiberBundleX/mitkFiberBundleXWriter.cpp @@ -1,101 +1,101 @@ /*========================================================================= Program: Medical Imaging & Interaction Toolkit Language: C++ Date: $Date: 2008-12-10 18:05:13 +0100 (Mi, 10 Dez 2008) $ Version: $Revision: 15922 $ 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 "mitkFiberBundleXWriter.h" #include #include mitk::FiberBundleXWriter::FiberBundleXWriter() : m_FileName(""), m_FilePrefix(""), m_FilePattern(""), m_Success(false) { this->SetNumberOfRequiredInputs( 1 ); } mitk::FiberBundleXWriter::~FiberBundleXWriter() {} void mitk::FiberBundleXWriter::GenerateData() { try { const std::string& locale = "C"; const std::string& currLocale = setlocale( LC_ALL, NULL ); setlocale(LC_ALL, locale.c_str()); MITK_INFO << "Writing fiber bundle"; m_Success = false; InputType* input = this->GetInput(); if (input == NULL) { itkWarningMacro(<<"Sorry, input to FiberBundleXWriter is NULL!"); return; } if ( m_FileName == "" ) { itkWarningMacro( << "Sorry, filename has not been set!" ); return ; } vtkSmartPointer writer = vtkSmartPointer::New(); - vtkSmartPointer cleaner = vtkSmartPointer::New(); - cleaner->SetInput(input->GetFiberPolyData()); - cleaner->Update(); +// vtkSmartPointer cleaner = vtkSmartPointer::New(); +// cleaner->SetInput(input->GetFiberPolyData()); +// cleaner->Update(); - writer->SetInput(cleaner->GetOutput()); + writer->SetInput(input->GetFiberPolyData()); writer->SetFileName(m_FileName.c_str()); writer->SetFileTypeToASCII(); writer->Write(); setlocale(LC_ALL, currLocale.c_str()); m_Success = true; MITK_INFO << "Fiber bundle written"; } catch(...) { throw; } } void mitk::FiberBundleXWriter::SetInputFiberBundleX( InputType* diffVolumes ) { this->ProcessObject::SetNthInput( 0, diffVolumes ); } mitk::FiberBundleX* mitk::FiberBundleXWriter::GetInput() { if ( this->GetNumberOfInputs() < 1 ) { return NULL; } else { return dynamic_cast ( this->ProcessObject::GetInput( 0 ) ); } } std::vector mitk::FiberBundleXWriter::GetPossibleFileExtensions() { std::vector possibleFileExtensions; possibleFileExtensions.push_back(".fib"); possibleFileExtensions.push_back(".vtk"); return possibleFileExtensions; } diff --git a/Modules/DiffusionImaging/Tractography/GibbsTracking/BuildFibres.cpp b/Modules/DiffusionImaging/Tractography/GibbsTracking/BuildFibres.cpp index 7abc049411..0907ef1cbc 100644 --- a/Modules/DiffusionImaging/Tractography/GibbsTracking/BuildFibres.cpp +++ b/Modules/DiffusionImaging/Tractography/GibbsTracking/BuildFibres.cpp @@ -1,179 +1,179 @@ #ifndef _BUILDFIBRES #define _BUILDFIBRES //#include "matrix.h" #define _USE_MATH_DEFINES #include #include #include #include using namespace std; #define PI M_PI #include "ParticleGrid.cpp" #include #include #include #include #include #include #include #include class FiberBuilder { public: Particle *particles; int pcnt; int attrcnt; typedef vector< Particle* > ParticleContainerType; typedef vector< ParticleContainerType* > FiberContainerType; vtkSmartPointer m_VtkCellArray; vtkSmartPointer m_VtkPoints; typedef itk::Vector OdfVectorType; typedef itk::Image ItkQBallImgType; ItkQBallImgType::Pointer m_ItkQBallImage; float m_FiberLength; itk::Point m_LastPoint; FiberBuilder(float *points, int numPoints, double spacing[], ItkQBallImgType::Pointer image) { m_FiberLength = 0; m_ItkQBallImage = image; particles = (Particle*) malloc(sizeof(Particle)*numPoints); pcnt = numPoints; attrcnt = 10; for (int k = 0; k < numPoints; k++) { Particle *p = &(particles[k]); p->R = pVector(points[attrcnt*k]/spacing[0]-0.5, points[attrcnt*k+1]/spacing[1]-0.5,points[attrcnt*k+2]/spacing[2]-0.5); p->N = pVector(points[attrcnt*k+3],points[attrcnt*k+4],points[attrcnt*k+5]); p->cap = points[attrcnt*k+6]; p->len = points[attrcnt*k+7]; p->mID = (int) points[attrcnt*k+8]; p->pID = (int) points[attrcnt*k+9]; p->ID = k; p->label = 0; } m_VtkCellArray = vtkSmartPointer::New(); m_VtkPoints = vtkSmartPointer::New(); } ~FiberBuilder() { free(particles); } vtkSmartPointer iterate(int minFiberLength) { int cur_label = 1; int numFibers = 0; m_FiberLength = 0; for (int k = 0; k < pcnt;k++) { Particle *dp = &(particles[k]); if (dp->label == 0) { vtkSmartPointer container = vtkSmartPointer::New(); dp->label = cur_label; dp->numerator = 0; labelPredecessors(dp, container); labelSuccessors(dp, container); cur_label++; if(m_FiberLength >= minFiberLength) { m_VtkCellArray->InsertNextCell(container); numFibers++; } m_FiberLength = 0; } } vtkSmartPointer fiberPolyData = vtkSmartPointer::New(); fiberPolyData->SetPoints(m_VtkPoints); fiberPolyData->SetLines(m_VtkCellArray); - vtkSmartPointer cleaner = vtkSmartPointer::New(); - cleaner->SetInput(fiberPolyData); - cleaner->Update(); - fiberPolyData = cleaner->GetOutput(); +// vtkSmartPointer cleaner = vtkSmartPointer::New(); +// cleaner->SetInput(fiberPolyData); +// cleaner->Update(); +// fiberPolyData = cleaner->GetOutput(); return fiberPolyData; } void AddPoint(Particle *dp, vtkSmartPointer container) { if (dp->inserted) return; dp->inserted = true; itk::ContinuousIndex index; index[0] = dp->R[0]; index[1] = dp->R[1]; index[2] = dp->R[2]; itk::Point point; m_ItkQBallImage->TransformContinuousIndexToPhysicalPoint( index, point ); vtkIdType id = m_VtkPoints->InsertNextPoint(point.GetDataPointer()); container->GetPointIds()->InsertNextId(id); if(container->GetNumberOfPoints()>1) m_FiberLength += m_LastPoint.EuclideanDistanceTo(point); m_LastPoint = point; } void labelPredecessors(Particle *dp, vtkSmartPointer container) { if (dp->mID != -1 && dp->mID!=dp->ID) { if (dp->ID!=particles[dp->mID].pID) { if (dp->ID==particles[dp->mID].mID) { int tmp = particles[dp->mID].pID; particles[dp->mID].pID = particles[dp->mID].mID; particles[dp->mID].mID = tmp; } } if (particles[dp->mID].label == 0) { particles[dp->mID].label = dp->label; particles[dp->mID].numerator = dp->numerator-1; labelPredecessors(&(particles[dp->mID]), container); } } AddPoint(dp, container); } void labelSuccessors(Particle *dp, vtkSmartPointer container) { AddPoint(dp, container); if (dp->pID != -1 && dp->pID!=dp->ID) { if (dp->ID!=particles[dp->pID].mID) { if (dp->ID==particles[dp->pID].pID) { int tmp = particles[dp->pID].pID; particles[dp->pID].pID = particles[dp->pID].mID; particles[dp->pID].mID = tmp; } } if (particles[dp->pID].label == 0) { particles[dp->pID].label = dp->label; particles[dp->pID].numerator = dp->numerator+1; labelSuccessors(&(particles[dp->pID]), container); } } } }; #endif