diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/GibbsTracking/mitkParticleGrid.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/GibbsTracking/mitkParticleGrid.cpp index 4f2333d592..49a67a5f82 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/GibbsTracking/mitkParticleGrid.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/GibbsTracking/mitkParticleGrid.cpp @@ -1,418 +1,417 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkParticleGrid.h" #include #include using namespace mitk; ParticleGrid::ParticleGrid(ItkFloatImageType* image, float particleLength, int cellCapacity) { // initialize counters m_NumParticles = 0; m_NumConnections = 0; m_NumCellOverflows = 0; m_ParticleLength = particleLength; // define isotropic grid from voxel spacing and particle length float cellSize = 2*m_ParticleLength; m_GridSize[0] = image->GetLargestPossibleRegion().GetSize()[0]*image->GetSpacing()[0]/cellSize +1; m_GridSize[1] = image->GetLargestPossibleRegion().GetSize()[1]*image->GetSpacing()[1]/cellSize +1; m_GridSize[2] = image->GetLargestPossibleRegion().GetSize()[2]*image->GetSpacing()[2]/cellSize +1; m_GridScale[0] = 1/cellSize; m_GridScale[1] = 1/cellSize; m_GridScale[2] = 1/cellSize; m_CellCapacity = cellCapacity; // maximum number of particles per grid cell m_ContainerCapacity = 100000; // initial particle container capacity unsigned long numCells = m_GridSize[0]*m_GridSize[1]*m_GridSize[2]; // number of grid cells unsigned long gridSize = numCells*m_CellCapacity; if ( itk::NumericTraits::max() R) { if (m_NumParticles >= m_ContainerCapacity) { if (!ReallocateGrid()) return NULL; } int xint = int(R[0]*m_GridScale[0]); if (xint < 0) return NULL; if (xint >= m_GridSize[0]) return NULL; int yint = int(R[1]*m_GridScale[1]); if (yint < 0) return NULL; if (yint >= m_GridSize[1]) return NULL; int zint = int(R[2]*m_GridScale[2]); if (zint < 0) return NULL; if (zint >= m_GridSize[2]) return NULL; int idx = xint + m_GridSize[0]*(yint + m_GridSize[1]*zint); if (m_OccupationCount[idx] < m_CellCapacity) { Particle *p = &(m_Particles[m_NumParticles]); p->pos = R; p->mID = -1; p->pID = -1; m_NumParticles++; p->gridindex = m_CellCapacity*idx + m_OccupationCount[idx]; m_Grid[p->gridindex] = p; m_OccupationCount[idx]++; return p; } else { m_NumCellOverflows++; return NULL; } } bool ParticleGrid::TryUpdateGrid(int k) { Particle* p = &(m_Particles[k]); int xint = int(p->pos[0]*m_GridScale[0]); if (xint < 0) return false; if (xint >= m_GridSize[0]) return false; int yint = int(p->pos[1]*m_GridScale[1]); if (yint < 0) return false; if (yint >= m_GridSize[1]) return false; int zint = int(p->pos[2]*m_GridScale[2]); if (zint < 0) return false; if (zint >= m_GridSize[2]) return false; int idx = xint + m_GridSize[0]*(yint+ zint*m_GridSize[1]); int cellidx = p->gridindex/m_CellCapacity; if (idx != cellidx) // cell has changed { if (m_OccupationCount[idx] < m_CellCapacity) { // remove from old position in grid; int grdindex = p->gridindex; m_Grid[grdindex] = m_Grid[cellidx*m_CellCapacity + m_OccupationCount[cellidx]-1]; m_Grid[grdindex]->gridindex = grdindex; m_OccupationCount[cellidx]--; // insert at new position in grid p->gridindex = idx*m_CellCapacity + m_OccupationCount[idx]; m_Grid[p->gridindex] = p; m_OccupationCount[idx]++; return true; } else { m_NumCellOverflows++; return false; } } return true; } void ParticleGrid::RemoveParticle(int k) { Particle* p = &(m_Particles[k]); int gridIndex = p->gridindex; int cellIdx = gridIndex/m_CellCapacity; int idx = gridIndex%m_CellCapacity; // remove pending connections if (p->mID != -1) DestroyConnection(p,-1); if (p->pID != -1) DestroyConnection(p,+1); // remove from grid if (idx < m_OccupationCount[cellIdx]-1) { m_Grid[gridIndex] = m_Grid[cellIdx*m_CellCapacity+m_OccupationCount[cellIdx]-1]; m_Grid[cellIdx*m_CellCapacity+m_OccupationCount[cellIdx]-1] = NULL; m_Grid[gridIndex]->gridindex = gridIndex; } m_OccupationCount[cellIdx]--; // remove from container if (k < m_NumParticles-1) { Particle* last = &m_Particles[m_NumParticles-1]; // last particle // update connections of last particle because its index is changing if (last->mID!=-1) { if ( m_Particles[last->mID].mID == m_NumParticles-1 ) m_Particles[last->mID].mID = k; else if ( m_Particles[last->mID].pID == m_NumParticles-1 ) m_Particles[last->mID].pID = k; } if (last->pID!=-1) { if ( m_Particles[last->pID].mID == m_NumParticles-1 ) m_Particles[last->pID].mID = k; else if ( m_Particles[last->pID].pID == m_NumParticles-1 ) m_Particles[last->pID].pID = k; } m_Particles[k] = m_Particles[m_NumParticles-1]; // move very last particle to empty slot m_Particles[m_NumParticles-1].ID = m_NumParticles-1; // update ID of removed particle to match the index m_Particles[k].ID = k; // update ID of moved particle m_Grid[m_Particles[k].gridindex] = &m_Particles[k]; // update address of moved particle } m_NumParticles--; } void ParticleGrid::ComputeNeighbors(vnl_vector_fixed &R) { float xfrac = R[0]*m_GridScale[0]; float yfrac = R[1]*m_GridScale[1]; float zfrac = R[2]*m_GridScale[2]; int xint = int(xfrac); int yint = int(yfrac); int zint = int(zfrac); int dx = -1; if (xfrac-xint > 0.5) dx = 1; if (xint <= 0) { xint = 0; dx = 1; } if (xint >= m_GridSize[0]-1) { xint = m_GridSize[0]-1; dx = -1; } int dy = -1; if (yfrac-yint > 0.5) dy = 1; if (yint <= 0) {yint = 0; dy = 1; } if (yint >= m_GridSize[1]-1) {yint = m_GridSize[1]-1; dy = -1;} int dz = -1; if (zfrac-zint > 0.5) dz = 1; if (zint <= 0) {zint = 0; dz = 1; } if (zint >= m_GridSize[2]-1) {zint = m_GridSize[2]-1; dz = -1;} m_NeighbourTracker.cellidx[0] = xint + m_GridSize[0]*(yint+zint*m_GridSize[1]); m_NeighbourTracker.cellidx[1] = m_NeighbourTracker.cellidx[0] + dx; m_NeighbourTracker.cellidx[2] = m_NeighbourTracker.cellidx[1] + dy*m_GridSize[0]; m_NeighbourTracker.cellidx[3] = m_NeighbourTracker.cellidx[2] - dx; m_NeighbourTracker.cellidx[4] = m_NeighbourTracker.cellidx[0] + dz*m_GridSize[0]*m_GridSize[1]; m_NeighbourTracker.cellidx[5] = m_NeighbourTracker.cellidx[4] + dx; m_NeighbourTracker.cellidx[6] = m_NeighbourTracker.cellidx[5] + dy*m_GridSize[0]; m_NeighbourTracker.cellidx[7] = m_NeighbourTracker.cellidx[6] - dx; m_NeighbourTracker.cellidx_c[0] = m_CellCapacity*m_NeighbourTracker.cellidx[0]; m_NeighbourTracker.cellidx_c[1] = m_CellCapacity*m_NeighbourTracker.cellidx[1]; m_NeighbourTracker.cellidx_c[2] = m_CellCapacity*m_NeighbourTracker.cellidx[2]; m_NeighbourTracker.cellidx_c[3] = m_CellCapacity*m_NeighbourTracker.cellidx[3]; m_NeighbourTracker.cellidx_c[4] = m_CellCapacity*m_NeighbourTracker.cellidx[4]; m_NeighbourTracker.cellidx_c[5] = m_CellCapacity*m_NeighbourTracker.cellidx[5]; m_NeighbourTracker.cellidx_c[6] = m_CellCapacity*m_NeighbourTracker.cellidx[6]; m_NeighbourTracker.cellidx_c[7] = m_CellCapacity*m_NeighbourTracker.cellidx[7]; m_NeighbourTracker.cellcnt = 0; m_NeighbourTracker.pcnt = 0; } Particle* ParticleGrid::GetNextNeighbor() { if (m_NeighbourTracker.pcnt < m_OccupationCount[m_NeighbourTracker.cellidx[m_NeighbourTracker.cellcnt]]) { return m_Grid[m_NeighbourTracker.cellidx_c[m_NeighbourTracker.cellcnt] + (m_NeighbourTracker.pcnt++)]; } else { for(;;) { m_NeighbourTracker.cellcnt++; if (m_NeighbourTracker.cellcnt >= 8) return 0; if (m_OccupationCount[m_NeighbourTracker.cellidx[m_NeighbourTracker.cellcnt]] > 0) break; } m_NeighbourTracker.pcnt = 1; return m_Grid[m_NeighbourTracker.cellidx_c[m_NeighbourTracker.cellcnt]]; } } void ParticleGrid::CreateConnection(Particle *P1,int ep1, Particle *P2, int ep2) { if (ep1 == -1) P1->mID = P2->ID; else P1->pID = P2->ID; if (ep2 == -1) P2->mID = P1->ID; else P2->pID = P1->ID; m_NumConnections++; } void ParticleGrid::DestroyConnection(Particle *P1,int ep1, Particle *P2, int ep2) { if (ep1 == -1) P1->mID = -1; else P1->pID = -1; if (ep2 == -1) P2->mID = -1; else P2->pID = -1; m_NumConnections--; } void ParticleGrid::DestroyConnection(Particle *P1,int ep1) { Particle *P2 = 0; if (ep1 == 1) { P2 = &m_Particles[P1->pID]; P1->pID = -1; } else { P2 = &m_Particles[P1->mID]; P1->mID = -1; } if (P2->mID == P1->ID) P2->mID = -1; else P2->pID = -1; m_NumConnections--; } bool ParticleGrid::CheckConsistency() { for (int i=0; iID != i) { std::cout << "Particle ID error!" << std::endl; return false; } if (p->mID!=-1) { Particle* p2 = &m_Particles[p->mID]; if (p2->mID!=p->ID && p2->pID!=p->ID) { std::cout << "Connection inconsistent!" << std::endl; return false; } } if (p->pID!=-1) { Particle* p2 = &m_Particles[p->pID]; if (p2->mID!=p->ID && p2->pID!=p->ID) { std::cout << "Connection inconsistent!" << std::endl; return false; } } } return true; } diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkAddArtifactsToDwiImageFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkAddArtifactsToDwiImageFilter.cpp index 3404b7c1f2..8ad4dee7c7 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkAddArtifactsToDwiImageFilter.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkAddArtifactsToDwiImageFilter.cpp @@ -1,273 +1,266 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef __itkAddArtifactsToDwiImageFilter_txx #define __itkAddArtifactsToDwiImageFilter_txx #include #include #include #include "itkAddArtifactsToDwiImageFilter.h" #include #include #include #include #include #include #define _USE_MATH_DEFINES #include namespace itk { template< class TPixelType > AddArtifactsToDwiImageFilter< TPixelType > ::AddArtifactsToDwiImageFilter() : m_NoiseModel(NULL) - , m_RingingModel(NULL) , m_FrequencyMap(NULL) , m_kOffset(0) , m_tLine(1) , m_EddyGradientStrength(0.001) , m_SimulateEddyCurrents(false) , m_TE(100) + , m_Upsampling(1) { this->SetNumberOfRequiredInputs( 1 ); } template< class TPixelType > AddArtifactsToDwiImageFilter< TPixelType >::ComplexSliceType::Pointer AddArtifactsToDwiImageFilter< TPixelType >::RearrangeSlice(ComplexSliceType::Pointer slice) { ImageRegion<2> region = slice->GetLargestPossibleRegion(); typename ComplexSliceType::Pointer rearrangedSlice = ComplexSliceType::New(); rearrangedSlice->SetLargestPossibleRegion( region ); rearrangedSlice->SetBufferedRegion( region ); rearrangedSlice->SetRequestedRegion( region ); rearrangedSlice->Allocate(); int xHalf = region.GetSize(0)/2; int yHalf = region.GetSize(1)/2; for (int y=0; y pix = slice->GetPixel(idx); if( idx[0] < xHalf ) idx[0] = idx[0] + xHalf; else idx[0] = idx[0] - xHalf; if( idx[1] < yHalf ) idx[1] = idx[1] + yHalf; else idx[1] = idx[1] - yHalf; rearrangedSlice->SetPixel(idx, pix); } return rearrangedSlice; } template< class TPixelType > void AddArtifactsToDwiImageFilter< TPixelType > ::GenerateData() { typename DiffusionImageType::Pointer inputImage = static_cast< DiffusionImageType * >( this->ProcessObject::GetInput(0) ); itk::ImageRegion<3> inputRegion = inputImage->GetLargestPossibleRegion(); typename itk::ImageDuplicator::Pointer duplicator = itk::ImageDuplicator::New(); duplicator->SetInputImage( inputImage ); duplicator->Update(); typename DiffusionImageType::Pointer outputImage = duplicator->GetOutput(); // is input slize size even? int x=inputRegion.GetSize(0); int y=inputRegion.GetSize(1); if ( x%2 == 1 ) x += 1; if ( y%2 == 1 ) y += 1; // create slice object typename SliceType::Pointer slice = SliceType::New(); ImageRegion<2> sliceRegion; sliceRegion.SetSize(0, x); sliceRegion.SetSize(1, y); slice->SetLargestPossibleRegion( sliceRegion ); slice->SetBufferedRegion( sliceRegion ); slice->SetRequestedRegion( sliceRegion ); slice->Allocate(); slice->FillBuffer(0.0); ImageRegion<2> upsampledSliceRegion; - unsigned int upsampling = 1; - if (m_RingingModel!=NULL) + if (m_Upsampling>1.00001) { - upsampling = m_RingingModel->GetKspaceCropping(); - upsampledSliceRegion.SetSize(0, x*upsampling); - upsampledSliceRegion.SetSize(1, y*upsampling); + upsampledSliceRegion.SetSize(0, x*m_Upsampling); + upsampledSliceRegion.SetSize(1, y*m_Upsampling); } // frequency map slice typename SliceType::Pointer fMap = NULL; if (m_FrequencyMap.IsNotNull()) { fMap = SliceType::New(); fMap->SetLargestPossibleRegion( sliceRegion ); fMap->SetBufferedRegion( sliceRegion ); fMap->SetRequestedRegion( sliceRegion ); fMap->Allocate(); fMap->FillBuffer(0.0); } - if ( m_FrequencyMap.IsNotNull() || m_kOffset>0.0 || m_RingingModel!=NULL || m_SimulateEddyCurrents) + if ( m_FrequencyMap.IsNotNull() || m_kOffset>0.0 || m_Upsampling>1.00001 || m_SimulateEddyCurrents) { MatrixType transform = inputImage->GetDirection(); for (int i=0; i<3; i++) for (int j=0; j<3; j++) transform[i][j] *= inputImage->GetSpacing()[j]; MITK_INFO << "Adjusting complex signal"; MITK_INFO << "line readout time: " << m_tLine; MITK_INFO << "line offset: " << m_kOffset; if (m_FrequencyMap.IsNotNull()) MITK_INFO << "frequency map is set"; else MITK_INFO << "no frequency map set"; - if (m_RingingModel!=NULL) + if (m_Upsampling>1.00001) MITK_INFO << "Gibbs ringing enabled"; else MITK_INFO << "Gibbs ringing disabled"; if (m_SimulateEddyCurrents) MITK_INFO << "Simulating eddy currents"; boost::progress_display disp(inputImage->GetVectorLength()*inputRegion.GetSize(2)); for (int g=0; gGetVectorLength(); g++) for (int z=0; z compartmentSlices; // extract slice from channel g for (int y=0; yGetPixel(index3D)[g]; slice->SetPixel(index2D, pix2D); if (fMap.IsNotNull()) fMap->SetPixel(index2D, m_FrequencyMap->GetPixel(index3D)); } - if (upsampling>1) + if (m_Upsampling>1.00001) { itk::ResampleImageFilter::Pointer resampler = itk::ResampleImageFilter::New(); resampler->SetInput(slice); resampler->SetOutputParametersFromImage(slice); resampler->SetSize(upsampledSliceRegion.GetSize()); - resampler->SetOutputSpacing(slice->GetSpacing()/upsampling); + resampler->SetOutputSpacing(slice->GetSpacing()/m_Upsampling); resampler->Update(); typename SliceType::Pointer upslice = resampler->GetOutput(); compartmentSlices.push_back(upslice); } else compartmentSlices.push_back(slice); // fourier transform slice typename ComplexSliceType::Pointer fSlice; - typename itk::KspaceImageFilter< SliceType::PixelType >::Pointer idft = itk::KspaceImageFilter< SliceType::PixelType >::New(); + itk::Size<2> outSize; outSize.SetElement(0, x); outSize.SetElement(1, y); + typename itk::KspaceImageFilter< SliceType::PixelType >::Pointer idft = itk::KspaceImageFilter< SliceType::PixelType >::New(); idft->SetCompartmentImages(compartmentSlices); idft->SetkOffset(m_kOffset); idft->SettLine(m_tLine); idft->SetSimulateRelaxation(false); idft->SetFrequencyMap(fMap); idft->SetDiffusionGradientDirection(m_GradientList.at(g)); idft->SetSimulateEddyCurrents(m_SimulateEddyCurrents); idft->SetEddyGradientMagnitude(m_EddyGradientStrength); idft->SetTE(m_TE); idft->SetZ((double)z-(double)inputRegion.GetSize(2)/2.0); idft->SetDirectionMatrix(transform); + idft->SetOutSize(outSize); idft->Update(); - fSlice = idft->GetOutput(); - if (m_RingingModel!=NULL) - { - fSlice = RearrangeSlice(fSlice); - fSlice = m_RingingModel->AddArtifact(fSlice); - } - // inverse fourier transform slice typename SliceType::Pointer newSlice; typename itk::DftImageFilter< SliceType::PixelType >::Pointer dft = itk::DftImageFilter< SliceType::PixelType >::New(); dft->SetInput(fSlice); dft->Update(); newSlice = dft->GetOutput(); // put slice back into channel g for (int y=0; yGetPixel(index3D); typename SliceType::IndexType index2D; index2D[0]=x; index2D[1]=y; double signal = newSlice->GetPixel(index2D); if (signal>0) signal = floor(signal+0.5); else signal = ceil(signal-0.5); pix3D[g] = signal; outputImage->SetPixel(index3D, pix3D); } + + ++disp; } } if (m_NoiseModel!=NULL) { ImageRegionIterator it1 (outputImage, inputRegion); boost::progress_display disp2(inputRegion.GetNumberOfPixels()); while(!it1.IsAtEnd()) { ++disp2; typename DiffusionImageType::PixelType signal = it1.Get(); m_NoiseModel->AddNoise(signal); it1.Set(signal); ++it1; } } this->SetNthOutput(0, outputImage); } } #endif diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkAddArtifactsToDwiImageFilter.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkAddArtifactsToDwiImageFilter.h index 8beccd5648..7a467447ba 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkAddArtifactsToDwiImageFilter.h +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkAddArtifactsToDwiImageFilter.h @@ -1,99 +1,97 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef __itkAddArtifactsToDwiImageFilter_h_ #define __itkAddArtifactsToDwiImageFilter_h_ #include "FiberTrackingExports.h" #include #include #include #include #include -#include #include namespace itk{ /** * \brief Adds several artifacts to the input DWI. */ template< class TPixelType > class AddArtifactsToDwiImageFilter : public ImageToImageFilter< VectorImage< TPixelType, 3 >, VectorImage< TPixelType, 3 > > { public: typedef AddArtifactsToDwiImageFilter Self; typedef SmartPointer Pointer; typedef SmartPointer ConstPointer; typedef ImageToImageFilter< VectorImage< TPixelType, 3 >, VectorImage< TPixelType, 3 > > Superclass; /** Method for creation through the object factory. */ itkNewMacro(Self) /** Runtime information support. */ itkTypeMacro(AddArtifactsToDwiImageFilter, ImageToImageFilter) typedef typename Superclass::InputImageType DiffusionImageType; typedef mitk::DiffusionNoiseModel NoiseModelType; - typedef mitk::GibbsRingingArtifact GibbsRingingType; typedef itk::Image< double, 2 > SliceType; typedef typename itk::KspaceImageFilter< double >::OutputImageType ComplexSliceType; typedef itk::Image ItkDoubleImgType; typedef itk::Matrix MatrixType; - void SetRingingModel(GibbsRingingType* ringingModel){ m_RingingModel = ringingModel; } void SetNoiseModel(NoiseModelType* noiseModel){ m_NoiseModel = noiseModel; } itkSetMacro( FrequencyMap, ItkDoubleImgType::Pointer ) itkSetMacro( kOffset, double ) itkSetMacro( tLine, double ) itkSetMacro( SimulateEddyCurrents, bool ) itkSetMacro( EddyGradientStrength, double ) void SetGradientList(mitk::DiffusionSignalModel::GradientListType list) { m_GradientList=list; } itkSetMacro( TE, double ) + itkSetMacro( Upsampling, double ) protected: AddArtifactsToDwiImageFilter(); ~AddArtifactsToDwiImageFilter() {} typename ComplexSliceType::Pointer RearrangeSlice(typename ComplexSliceType::Pointer slice); void GenerateData(); - GibbsRingingType* m_RingingModel; NoiseModelType* m_NoiseModel; ItkDoubleImgType::Pointer m_FrequencyMap; double m_kOffset; double m_tLine; bool m_SimulateEddyCurrents; double m_EddyGradientStrength; mitk::DiffusionSignalModel::GradientListType m_GradientList; double m_TE; + double m_Upsampling; ///< causes ringing artifacts private: }; } #ifndef ITK_MANUAL_INSTANTIATION #include "itkAddArtifactsToDwiImageFilter.cpp" #endif #endif //__itkAddArtifactsToDwiImageFilter_h_ diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkEvaluateDirectionImagesFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkEvaluateDirectionImagesFilter.cpp new file mode 100755 index 0000000000..5fc4b37e71 --- /dev/null +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkEvaluateDirectionImagesFilter.cpp @@ -0,0 +1,377 @@ +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +All rights reserved. + +This software is distributed WITHOUT ANY WARRANTY; without +even the implied warranty of MERCHANTABILITY or FITNESS FOR +A PARTICULAR PURPOSE. + +See LICENSE.txt or http://www.mitk.org for details. + +===================================================================*/ + +#ifndef __itkEvaluateDirectionImagesFilter_cpp +#define __itkEvaluateDirectionImagesFilter_cpp + +#include "itkEvaluateDirectionImagesFilter.h" +#include +#include +#include +#include + +#define _USE_MATH_DEFINES +#include + +namespace itk { + +template< class PixelType > +EvaluateDirectionImagesFilter< PixelType > +::EvaluateDirectionImagesFilter(): + m_ImageSet(NULL), + m_ReferenceImageSet(NULL), + m_IgnoreMissingDirections(false), + m_Eps(0.0001) +{ + this->SetNumberOfOutputs(2); +} + +template< class PixelType > +void EvaluateDirectionImagesFilter< PixelType >::GenerateData() +{ + if (m_ImageSet.IsNull() || m_ReferenceImageSet.IsNull()) + return; + + DirectionImageContainerType::Pointer set1 = DirectionImageContainerType::New(); + DirectionImageContainerType::Pointer set2 = DirectionImageContainerType::New(); + + for (int i=0; iSize(); i++) + { + typename itk::ImageDuplicator< DirectionImageType >::Pointer duplicator = itk::ImageDuplicator< DirectionImageType >::New(); + duplicator->SetInputImage( m_ImageSet->GetElement(i) ); + duplicator->Update(); + set1->InsertElement(i, dynamic_cast(duplicator->GetOutput())); + } + for (int i=0; iSize(); i++) + { + typename itk::ImageDuplicator< DirectionImageType >::Pointer duplicator = itk::ImageDuplicator< DirectionImageType >::New(); + duplicator->SetInputImage( m_ReferenceImageSet->GetElement(i) ); + duplicator->Update(); + set2->InsertElement(i, dynamic_cast(duplicator->GetOutput())); + } + + m_ImageSet = set1; + m_ReferenceImageSet = set2; + + // angular error image + typename OutputImageType::Pointer outputImage = OutputImageType::New(); + outputImage->SetOrigin( m_ReferenceImageSet->GetElement(0)->GetOrigin() ); + outputImage->SetRegions( m_ReferenceImageSet->GetElement(0)->GetLargestPossibleRegion() ); + outputImage->SetSpacing( m_ReferenceImageSet->GetElement(0)->GetSpacing() ); + outputImage->SetDirection( m_ReferenceImageSet->GetElement(0)->GetDirection() ); + outputImage->Allocate(); + outputImage->FillBuffer(0.0); + this->SetNthOutput(0, outputImage); + + // length error image + outputImage = OutputImageType::New(); + outputImage->SetOrigin( m_ReferenceImageSet->GetElement(0)->GetOrigin() ); + outputImage->SetRegions( m_ReferenceImageSet->GetElement(0)->GetLargestPossibleRegion() ); + outputImage->SetSpacing( m_ReferenceImageSet->GetElement(0)->GetSpacing() ); + outputImage->SetDirection( m_ReferenceImageSet->GetElement(0)->GetDirection() ); + outputImage->Allocate(); + outputImage->FillBuffer(0.0); + this->SetNthOutput(1, outputImage); + + if (m_MaskImage.IsNull()) + { + m_MaskImage = UCharImageType::New(); + m_MaskImage->SetOrigin( outputImage->GetOrigin() ); + m_MaskImage->SetRegions( outputImage->GetLargestPossibleRegion() ); + m_MaskImage->SetSpacing( outputImage->GetSpacing() ); + m_MaskImage->SetDirection( outputImage->GetDirection() ); + m_MaskImage->Allocate(); + m_MaskImage->FillBuffer(1); + } + + m_MeanAngularError = 0.0; + m_MedianAngularError = 0; + m_MaxAngularError = 0.0; + m_MinAngularError = itk::NumericTraits::max(); + m_VarAngularError = 0.0; + m_AngularErrorVector.clear(); + + m_MeanLengthError = 0.0; + m_MedianLengthError = 0; + m_MaxLengthError = 0.0; + m_MinLengthError = itk::NumericTraits::max(); + m_VarLengthError = 0.0; + m_LengthErrorVector.clear(); + + if (m_ImageSet.IsNull() || m_ReferenceImageSet.IsNull()) + return; + + outputImage = static_cast< OutputImageType* >(this->ProcessObject::GetOutput(0)); + typename OutputImageType::Pointer outputImage2 = static_cast< OutputImageType* >(this->ProcessObject::GetOutput(1)); + + ImageRegionIterator< OutputImageType > oit(outputImage, outputImage->GetLargestPossibleRegion()); + ImageRegionIterator< OutputImageType > oit2(outputImage2, outputImage2->GetLargestPossibleRegion()); + ImageRegionIterator< UCharImageType > mit(m_MaskImage, m_MaskImage->GetLargestPossibleRegion()); + + int numImages = m_ImageSet->Size(); + int numReferenceImages = m_ReferenceImageSet->Size(); + + // fill missing directions with zeros + if (numImages>numReferenceImages) + { + DirectionType zeroDir; zeroDir.Fill(0.0); + for (int i=0; iSetOrigin( m_ReferenceImageSet->GetElement(0)->GetOrigin() ); + img->SetRegions( m_ReferenceImageSet->GetElement(0)->GetLargestPossibleRegion() ); + img->SetSpacing( m_ReferenceImageSet->GetElement(0)->GetSpacing() ); + img->SetDirection( m_ReferenceImageSet->GetElement(0)->GetDirection() ); + img->Allocate(); + img->FillBuffer(zeroDir); + m_ReferenceImageSet->InsertElement(m_ReferenceImageSet->Size(), img); + } + numReferenceImages = numImages; + } + else if (numReferenceImages>numImages) + { + DirectionType zeroDir; zeroDir.Fill(0.0); + for (int i=0; iSetOrigin( m_ReferenceImageSet->GetElement(0)->GetOrigin() ); + img->SetRegions( m_ReferenceImageSet->GetElement(0)->GetLargestPossibleRegion() ); + img->SetSpacing( m_ReferenceImageSet->GetElement(0)->GetSpacing() ); + img->SetDirection( m_ReferenceImageSet->GetElement(0)->GetDirection() ); + img->Allocate(); + img->FillBuffer(zeroDir); + m_ImageSet->InsertElement(m_ImageSet->Size(), img); + } + numImages = numReferenceImages; + } + int numDirections = numReferenceImages; + + // matrix containing the angular error between the directions + vnl_matrix< float > diffM; diffM.set_size(numDirections, numDirections); + + boost::progress_display disp(outputImage->GetLargestPossibleRegion().GetSize()[0]*outputImage->GetLargestPossibleRegion().GetSize()[1]*outputImage->GetLargestPossibleRegion().GetSize()[2]); + while( !oit.IsAtEnd() ) + { + ++disp; + if( mit.Get()!=1 ) + { + ++oit; + ++oit2; + ++mit; + continue; + } + typename OutputImageType::IndexType index = oit.GetIndex(); + + float maxAngularError = 1.0; + + diffM.fill(10); // initialize with invalid error value + + // get number of valid directions (length > 0) + int numRefDirs = 0; + int numTestDirs = 0; + for (int i=0; iGetElement(i)->GetPixel(index).GetVnlVector().magnitude() > m_Eps ) + numRefDirs++; + if (m_ImageSet->GetElement(i)->GetPixel(index).GetVnlVector().magnitude() > m_Eps ) + numTestDirs++; + } + + // i: index of reference direction + // j: index of test direction + for (int i=0; i refDir = m_ReferenceImageSet->GetElement(i)->GetPixel(index).GetVnlVector(); + + if (i testDir = m_ImageSet->GetElement(j)->GetPixel(index).GetVnlVector(); + if (j1.0) + diffM[i][j] = 1.0; + } + } + + float angularError = 0.0; + float lengthError = 0.0; + int counter = 0; + vnl_matrix< float > diffM_copy = diffM; + for (int k=0; k small error) + for (int i=0; ierror && diffM[i][j]<2) // found valid error entry + { + error = diffM[i][j]; + a = i; + b = j; + missingDir = false; + } + else if (diffM[i][j]<0 && error<0) // found missing direction + { + a = i; + b = j; + missingDir = true; + } + } + + if (a<0 || b<0 || m_IgnoreMissingDirections && missingDir) + continue; // no more directions found + + if (a>=numRefDirs && b>=numTestDirs) + { + MITK_INFO << "ERROR: missing test and reference direction. should not be possible. check code."; + continue; + } + + // remove processed directions from error matrix + diffM.set_row(a, 10.0); + diffM.set_column(b, 10.0); + + if (a>=numRefDirs) // missing reference direction (find next closest) + { + for (int i=0; ierror) + { + error = diffM_copy[i][b]; + a = i; + } + } + else if (b>=numTestDirs) // missing test direction (find next closest) + { + for (int i=0; ierror) + { + error = diffM_copy[a][i]; + b = i; + } + } + + float refLength = m_ReferenceImageSet->GetElement(a)->GetPixel(index).GetVnlVector().magnitude(); + float testLength = m_ImageSet->GetElement(b)->GetPixel(index).GetVnlVector().magnitude(); + + if (a>=numRefDirs || b>=numTestDirs || error<0) + error = 0; + + m_LengthErrorVector.push_back( fabs(refLength-testLength) ); + m_AngularErrorVector.push_back( acos(error)*180.0/M_PI ); + + m_MeanAngularError += m_AngularErrorVector.back(); + m_MeanLengthError += m_LengthErrorVector.back(); + + angularError += m_AngularErrorVector.back(); + lengthError += m_LengthErrorVector.back(); + counter++; + } + + if (counter>0) + { + lengthError /= counter; + angularError /= counter; + } + oit2.Set(lengthError); + oit.Set(angularError); + + ++oit; + ++oit2; + ++mit; + } + + std::sort( m_AngularErrorVector.begin(), m_AngularErrorVector.end() ); + m_MeanAngularError /= m_AngularErrorVector.size(); // mean + for (int i=0; im_MaxAngularError ) + m_MaxAngularError = m_AngularErrorVector.at(i); + if ( m_AngularErrorVector.at(i)1) + { + m_VarAngularError /= (m_AngularErrorVector.size()-1); // variance + + // median + if (m_AngularErrorVector.size()%2 == 0) + m_MedianAngularError = 0.5*( m_AngularErrorVector.at( m_AngularErrorVector.size()/2 ) + m_AngularErrorVector.at( m_AngularErrorVector.size()/2+1 ) ); + else + m_MedianAngularError = m_AngularErrorVector.at( (m_AngularErrorVector.size()+1)/2 ) ; + } + + std::sort( m_LengthErrorVector.begin(), m_LengthErrorVector.end() ); + m_MeanLengthError /= m_LengthErrorVector.size(); // mean + for (int i=0; im_MaxLengthError ) + m_MaxLengthError = m_LengthErrorVector.at(i); + if ( m_LengthErrorVector.at(i)1) + { + m_VarLengthError /= (m_LengthErrorVector.size()-1); // variance + + // median + if (m_LengthErrorVector.size()%2 == 0) + m_MedianLengthError = 0.5*( m_LengthErrorVector.at( m_LengthErrorVector.size()/2 ) + m_LengthErrorVector.at( m_LengthErrorVector.size()/2+1 ) ); + else + m_MedianLengthError = m_LengthErrorVector.at( (m_LengthErrorVector.size()+1)/2 ) ; + } +} + +} + +#endif // __itkEvaluateDirectionImagesFilter_cpp diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkEvaluateDirectionImagesFilter.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkEvaluateDirectionImagesFilter.h new file mode 100755 index 0000000000..964c17e5fc --- /dev/null +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkEvaluateDirectionImagesFilter.h @@ -0,0 +1,110 @@ +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +All rights reserved. + +This software is distributed WITHOUT ANY WARRANTY; without +even the implied warranty of MERCHANTABILITY or FITNESS FOR +A PARTICULAR PURPOSE. + +See LICENSE.txt or http://www.mitk.org for details. + +===================================================================*/ + +/*=================================================================== + +This file is based heavily on a corresponding ITK filter. + +===================================================================*/ +#ifndef __itkEvaluateDirectionImagesFilter_h_ +#define __itkEvaluateDirectionImagesFilter_h_ + +#include +#include +#include + +namespace itk{ +/** \class EvaluateDirectionImagesFilter + */ + +template< class PixelType > +class EvaluateDirectionImagesFilter : public ImageSource< Image< PixelType, 3 > > +{ + +public: + + typedef EvaluateDirectionImagesFilter Self; + typedef SmartPointer Pointer; + typedef SmartPointer ConstPointer; + typedef ImageSource< Image< PixelType, 3 > > Superclass; + typedef typename Superclass::OutputImageRegionType OutputImageRegionType; + typedef typename Superclass::OutputImageType OutputImageType; + + /** Method for creation through the object factory. */ + itkNewMacro(Self) + + /** Runtime information support. */ + itkTypeMacro(EvaluateDirectionImagesFilter, ImageToImageFilter) + + typedef Vector< float, 3 > DirectionType; + typedef Image< DirectionType, 3 > DirectionImageType; + typedef VectorContainer< int, DirectionImageType::Pointer > DirectionImageContainerType; + typedef Image< float, 3 > FloatImageType; + typedef Image< bool, 3 > BoolImageType; + typedef Image< unsigned char, 3 > UCharImageType; + + itkSetMacro( ImageSet , DirectionImageContainerType::Pointer) + itkSetMacro( ReferenceImageSet , DirectionImageContainerType::Pointer) + itkSetMacro( MaskImage , UCharImageType::Pointer) + itkSetMacro( IgnoreMissingDirections , bool) + + itkGetMacro( MeanAngularError, float) + itkGetMacro( MinAngularError, float) + itkGetMacro( MaxAngularError, float) + itkGetMacro( VarAngularError, float) + itkGetMacro( MedianAngularError, float) + + itkGetMacro( MeanLengthError, float) + itkGetMacro( MinLengthError, float) + itkGetMacro( MaxLengthError, float) + itkGetMacro( VarLengthError, float) + itkGetMacro( MedianLengthError, float) + +protected: + EvaluateDirectionImagesFilter(); + ~EvaluateDirectionImagesFilter() {} + + void GenerateData(); + + UCharImageType::Pointer m_MaskImage; + DirectionImageContainerType::Pointer m_ImageSet; + DirectionImageContainerType::Pointer m_ReferenceImageSet; + bool m_IgnoreMissingDirections; + double m_MeanAngularError; + double m_MedianAngularError; + double m_MaxAngularError; + double m_MinAngularError; + double m_VarAngularError; + std::vector< double > m_AngularErrorVector; + + double m_MeanLengthError; + double m_MedianLengthError; + double m_MaxLengthError; + double m_MinLengthError; + double m_VarLengthError; + std::vector< double > m_LengthErrorVector; + + double m_Eps; +}; + +} + +#ifndef ITK_MANUAL_INSTANTIATION +#include "itkEvaluateDirectionImagesFilter.cpp" +#endif + +#endif //__itkEvaluateDirectionImagesFilter_h_ + diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkKspaceImageFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkKspaceImageFilter.cpp index f61bfab819..80faef1148 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkKspaceImageFilter.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkKspaceImageFilter.cpp @@ -1,213 +1,236 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef __itkKspaceImageFilter_txx #define __itkKspaceImageFilter_txx #include #include #include #include "itkKspaceImageFilter.h" #include #include #include #include #define _USE_MATH_DEFINES #include namespace itk { template< class TPixelType > KspaceImageFilter< TPixelType > ::KspaceImageFilter() : m_tLine(1) , m_kOffset(0) , m_FrequencyMap(NULL) , m_SimulateRelaxation(true) , m_SimulateEddyCurrents(false) , m_Tau(70) , m_EddyGradientMagnitude(30) , m_IsBaseline(true) , m_SignalScale(1) { m_DiffusionGradientDirection.Fill(0.0); } template< class TPixelType > void KspaceImageFilter< TPixelType > ::BeforeThreadedGenerateData() { typename OutputImageType::Pointer outputImage = OutputImageType::New(); outputImage->SetSpacing( m_CompartmentImages.at(0)->GetSpacing() ); outputImage->SetOrigin( m_CompartmentImages.at(0)->GetOrigin() ); outputImage->SetDirection( m_CompartmentImages.at(0)->GetDirection() ); - outputImage->SetLargestPossibleRegion( m_CompartmentImages.at(0)->GetLargestPossibleRegion() ); - outputImage->SetBufferedRegion( m_CompartmentImages.at(0)->GetLargestPossibleRegion() ); - outputImage->SetRequestedRegion( m_CompartmentImages.at(0)->GetLargestPossibleRegion() ); + itk::ImageRegion<2> region; region.SetSize(0, m_OutSize[0]); region.SetSize(1, m_OutSize[1]); + outputImage->SetLargestPossibleRegion( region ); + outputImage->SetBufferedRegion( region ); + outputImage->SetRequestedRegion( region ); outputImage->Allocate(); - typename InputImageType::Pointer tempImage = InputImageType::New(); - tempImage->SetSpacing( m_CompartmentImages.at(0)->GetSpacing() ); - tempImage->SetOrigin( m_CompartmentImages.at(0)->GetOrigin() ); - tempImage->SetDirection( m_CompartmentImages.at(0)->GetDirection() ); - tempImage->SetLargestPossibleRegion( m_CompartmentImages.at(0)->GetLargestPossibleRegion() ); - tempImage->SetBufferedRegion( m_CompartmentImages.at(0)->GetLargestPossibleRegion() ); - tempImage->SetRequestedRegion( m_CompartmentImages.at(0)->GetLargestPossibleRegion() ); - tempImage->Allocate(); + m_TEMPIMAGE = InputImageType::New(); + m_TEMPIMAGE->SetSpacing( m_CompartmentImages.at(0)->GetSpacing() ); + m_TEMPIMAGE->SetOrigin( m_CompartmentImages.at(0)->GetOrigin() ); + m_TEMPIMAGE->SetDirection( m_CompartmentImages.at(0)->GetDirection() ); + m_TEMPIMAGE->SetLargestPossibleRegion( region ); + m_TEMPIMAGE->SetBufferedRegion( region ); + m_TEMPIMAGE->SetRequestedRegion( region ); + m_TEMPIMAGE->Allocate(); m_SimulateDistortions = true; if (m_FrequencyMap.IsNull()) { m_SimulateDistortions = false; m_FrequencyMap = InputImageType::New(); - m_FrequencyMap->SetSpacing( outputImage->GetSpacing() ); - m_FrequencyMap->SetOrigin( outputImage->GetOrigin() ); - m_FrequencyMap->SetDirection( outputImage->GetDirection() ); - m_FrequencyMap->SetLargestPossibleRegion( outputImage->GetLargestPossibleRegion() ); - m_FrequencyMap->SetBufferedRegion( outputImage->GetLargestPossibleRegion() ); - m_FrequencyMap->SetRequestedRegion( outputImage->GetLargestPossibleRegion() ); + m_FrequencyMap->SetSpacing( m_CompartmentImages.at(0)->GetSpacing() ); + m_FrequencyMap->SetOrigin( m_CompartmentImages.at(0)->GetOrigin() ); + m_FrequencyMap->SetDirection( m_CompartmentImages.at(0)->GetDirection() ); + m_FrequencyMap->SetLargestPossibleRegion( m_CompartmentImages.at(0)->GetLargestPossibleRegion() ); + m_FrequencyMap->SetBufferedRegion( m_CompartmentImages.at(0)->GetLargestPossibleRegion() ); + m_FrequencyMap->SetRequestedRegion( m_CompartmentImages.at(0)->GetLargestPossibleRegion() ); m_FrequencyMap->Allocate(); m_FrequencyMap->FillBuffer(0); } double gamma = 42576000; // Gyromagnetic ratio in Hz/T if (m_DiffusionGradientDirection.GetNorm()>0.001) { m_DiffusionGradientDirection.Normalize(); m_EddyGradientMagnitude /= 1000; // eddy gradient magnitude in T/m m_DiffusionGradientDirection = m_DiffusionGradientDirection * m_EddyGradientMagnitude * gamma; m_IsBaseline = false; } else m_EddyGradientMagnitude = gamma*m_EddyGradientMagnitude/1000; this->SetNthOutput(0, outputImage); - this->SetNthOutput(1, tempImage); } template< class TPixelType > void KspaceImageFilter< TPixelType > ::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, ThreadIdType threadId) { typename OutputImageType::Pointer outputImage = static_cast< OutputImageType * >(this->ProcessObject::GetOutput(0)); - typename InputImageType::Pointer tempImage = static_cast< InputImageType * >(this->ProcessObject::GetOutput(1)); ImageRegionIterator< OutputImageType > oit(outputImage, outputRegionForThread); typedef ImageRegionConstIterator< InputImageType > InputIteratorType; double szx = outputImage->GetLargestPossibleRegion().GetSize(0); double szy = outputImage->GetLargestPossibleRegion().GetSize(1); double numPix = szx*szy; double dt = m_tLine/szx; double fromMaxEcho = - m_tLine*szy/2; + double in_szx = m_CompartmentImages.at(0)->GetLargestPossibleRegion().GetSize(0); + double in_szy = m_CompartmentImages.at(0)->GetLargestPossibleRegion().GetSize(1); + int xOffset = in_szx-szx; + int yOffset = in_szy-szy; + while( !oit.IsAtEnd() ) { itk::Index< 2 > kIdx; kIdx[0] = oit.GetIndex()[0]; kIdx[1] = oit.GetIndex()[1]; double t = fromMaxEcho + ((double)kIdx[1]*szx+(double)kIdx[0])*dt; // dephasing time // calculate eddy current decay factors double eddyDecay = 0; if (m_SimulateEddyCurrents) eddyDecay = exp(-(m_TE/2 + t)/m_Tau) * t/1000; // calcualte signal relaxation factors std::vector< double > relaxFactor; if (m_SimulateRelaxation) { for (int i=0; iSetPixel(kIdx, t ); + temp_kx += m_kOffset; // add gradient delay induced offset vcl_complex s(0,0); InputIteratorType it(m_CompartmentImages.at(0), m_CompartmentImages.at(0)->GetLargestPossibleRegion() ); while( !it.IsAtEnd() ) { double x = it.GetIndex()[0]; double y = it.GetIndex()[1]; vcl_complex f(0, 0); // sum compartment signals and simulate relaxation for (int i=0; i( m_CompartmentImages.at(i)->GetPixel(it.GetIndex()) * relaxFactor.at(i) * m_SignalScale, 0); else f += std::complex( m_CompartmentImages.at(i)->GetPixel(it.GetIndex()) * m_SignalScale ); // simulate eddy currents and other distortions double omega_t = 0; if ( m_SimulateEddyCurrents ) { if (!m_IsBaseline) { itk::Vector< double, 3 > pos; pos[0] = x-szx/2; pos[1] = y-szy/2; pos[2] = m_Z; pos = m_DirectionMatrix*pos/1000; // vector from image center to current position (in meter) omega_t += (m_DiffusionGradientDirection[0]*pos[0]+m_DiffusionGradientDirection[1]*pos[1]+m_DiffusionGradientDirection[2]*pos[2])*eddyDecay; } else omega_t += m_EddyGradientMagnitude*eddyDecay; } if (m_SimulateDistortions) omega_t += m_FrequencyMap->GetPixel(it.GetIndex())*t/1000; + // add gibbs ringing offset (cropps k-space) + double tempOffsetX = 0; + if (temp_kx>=szx/2) + tempOffsetX = xOffset; + double temp_ky = kIdx[1]; + if (temp_ky>=szy/2) + temp_ky += yOffset; + // actual DFT term - s += f * exp( std::complex(0, 2 * M_PI * (temp_k*x/szx + (double)kIdx[1]*y/szy) + omega_t ) ); + s += f * exp( std::complex(0, 2 * M_PI * ((temp_kx+tempOffsetX)*x/in_szx + temp_ky*y/in_szy) + omega_t ) ); ++it; } s /= numPix; + + /* rearrange slice not needed + if( kIdx[0] < szx/2 ) + kIdx[0] = kIdx[0] + szx/2; + else + kIdx[0] = kIdx[0] - szx/2; + + if( kIdx[1] < szx/2 ) + kIdx[1] = kIdx[1] + szy/2; + else + kIdx[1] = kIdx[1] - szy/2;*/ + + m_TEMPIMAGE->SetPixel(kIdx, sqrt(s.real()*s.real()+s.imag()*s.imag())); outputImage->SetPixel(kIdx, s); ++oit; } // typedef itk::ImageFileWriter< InputImageType > WriterType; // typename WriterType::Pointer writer = WriterType::New(); -// writer->SetFileName("/local/dist.nrrd"); -// writer->SetInput(tempImage); +// writer->SetFileName("/local/kspace.nrrd"); +// writer->SetInput(m_TEMPIMAGE); // writer->Update(); } template< class TPixelType > void KspaceImageFilter< TPixelType > ::AfterThreadedGenerateData() { } } #endif diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkKspaceImageFilter.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkKspaceImageFilter.h index 9da9b12400..01bdb7d9bb 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkKspaceImageFilter.h +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkKspaceImageFilter.h @@ -1,115 +1,118 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ /*=================================================================== This file is based heavily on a corresponding ITK filter. ===================================================================*/ #ifndef __itkKspaceImageFilter_h_ #define __itkKspaceImageFilter_h_ #include "FiberTrackingExports.h" #include #include #include namespace itk{ /** * \brief Performes deterministic streamline tracking on the input tensor image. */ template< class TPixelType > class KspaceImageFilter : public ImageSource< Image< vcl_complex< TPixelType >, 2 > > { public: typedef KspaceImageFilter Self; typedef SmartPointer Pointer; typedef SmartPointer ConstPointer; typedef ImageSource< Image< vcl_complex< TPixelType >, 2 > > Superclass; /** Method for creation through the object factory. */ - itkNewMacro(Self); + itkNewMacro(Self) /** Runtime information support. */ itkTypeMacro(KspaceImageFilter, ImageToImageFilter) typedef typename itk::Image< double, 2 > InputImageType; typedef typename InputImageType::Pointer InputImagePointerType; typedef typename Superclass::OutputImageType OutputImageType; typedef typename Superclass::OutputImageRegionType OutputImageRegionType; typedef itk::Matrix MatrixType; itkSetMacro( FrequencyMap, typename InputImageType::Pointer ) itkSetMacro( tLine, double ) itkSetMacro( kOffset, double ) itkSetMacro( TE, double) itkSetMacro( Tinhom, double) itkSetMacro( Tau, double) itkSetMacro( SimulateRelaxation, bool ) itkSetMacro( SimulateEddyCurrents, bool ) itkSetMacro( Z, double ) itkSetMacro( DirectionMatrix, MatrixType ) itkSetMacro( SignalScale, double ) + itkSetMacro( OutSize, itk::Size<2> ) void SetT2( std::vector< double > t2Vector ) { m_T2=t2Vector; } void SetCompartmentImages( std::vector< InputImagePointerType > cImgs ) { m_CompartmentImages=cImgs; } void SetDiffusionGradientDirection(itk::Vector g) { m_DiffusionGradientDirection=g; } void SetEddyGradientMagnitude(double g_mag) { m_EddyGradientMagnitude=g_mag; } ///< in T/m protected: KspaceImageFilter(); ~KspaceImageFilter() {} void BeforeThreadedGenerateData(); void ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, ThreadIdType threadId); void AfterThreadedGenerateData(); bool m_SimulateRelaxation; bool m_SimulateDistortions; bool m_SimulateEddyCurrents; + typename InputImageType::Pointer m_TEMPIMAGE; typename InputImageType::Pointer m_FrequencyMap; double m_tLine; double m_kOffset; double m_Tinhom; double m_TE; std::vector< double > m_T2; std::vector< InputImagePointerType > m_CompartmentImages; itk::Vector m_DiffusionGradientDirection; double m_Tau; ///< eddy current decay constant double m_EddyGradientMagnitude; ///< in T/m double m_Z; MatrixType m_DirectionMatrix; bool m_IsBaseline; double m_SignalScale; + itk::Size<2> m_OutSize; private: }; } #ifndef ITK_MANUAL_INSTANTIATION #include "itkKspaceImageFilter.cpp" #endif #endif //__itkKspaceImageFilter_h_ diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToDWIImageFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToDWIImageFilter.cpp old mode 100644 new mode 100755 index 772cbcda29..9a3f838d71 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToDWIImageFilter.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToDWIImageFilter.cpp @@ -1,709 +1,698 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "itkTractsToDWIImageFilter.h" #include #include #include #include #include #include #include -#include #include #include #include #include #include #include #include #include #include #include #include namespace itk { TractsToDWIImageFilter::TractsToDWIImageFilter() : m_CircleDummy(false) , m_VolumeAccuracy(10) , m_Upsampling(1) , m_NumberOfRepetitions(1) , m_EnforcePureFiberVoxels(false) , m_InterpolationShrink(1000) , m_FiberRadius(0) , m_SignalScale(25) , m_kOffset(0) , m_tLine(1) , m_UseInterpolation(false) , m_SimulateRelaxation(true) , m_tInhom(50) , m_TE(100) , m_FrequencyMap(NULL) , m_EddyGradientStrength(0.001) , m_SimulateEddyCurrents(false) { m_Spacing.Fill(2.5); m_Origin.Fill(0.0); m_DirectionMatrix.SetIdentity(); m_ImageRegion.SetSize(0, 10); m_ImageRegion.SetSize(1, 10); m_ImageRegion.SetSize(2, 10); } TractsToDWIImageFilter::~TractsToDWIImageFilter() { } TractsToDWIImageFilter::DoubleDwiType::Pointer TractsToDWIImageFilter::DoKspaceStuff( std::vector< DoubleDwiType::Pointer >& images ) { // create slice object ImageRegion<2> sliceRegion; sliceRegion.SetSize(0, m_UpsampledImageRegion.GetSize()[0]); sliceRegion.SetSize(1, m_UpsampledImageRegion.GetSize()[1]); // frequency map slice SliceType::Pointer fMap = NULL; if (m_FrequencyMap.IsNotNull()) { fMap = SliceType::New(); ImageRegion<2> region; region.SetSize(0, m_UpsampledImageRegion.GetSize()[0]); region.SetSize(1, m_UpsampledImageRegion.GetSize()[1]); fMap->SetLargestPossibleRegion( region ); fMap->SetBufferedRegion( region ); fMap->SetRequestedRegion( region ); fMap->Allocate(); } DoubleDwiType::Pointer newImage = DoubleDwiType::New(); newImage->SetSpacing( m_Spacing ); newImage->SetOrigin( m_Origin ); newImage->SetDirection( m_DirectionMatrix ); newImage->SetLargestPossibleRegion( m_ImageRegion ); newImage->SetBufferedRegion( m_ImageRegion ); newImage->SetRequestedRegion( m_ImageRegion ); newImage->SetVectorLength( images.at(0)->GetVectorLength() ); newImage->Allocate(); MatrixType transform = m_DirectionMatrix; for (int i=0; i<3; i++) for (int j=0; j<3; j++) { if (j<2) transform[i][j] *= m_UpsampledSpacing[j]; else transform[i][j] *= m_Spacing[j]; } boost::progress_display disp(images.at(0)->GetVectorLength()*images.at(0)->GetLargestPossibleRegion().GetSize(2)); for (int g=0; gGetVectorLength(); g++) for (int z=0; zGetLargestPossibleRegion().GetSize(2); z++) { - ++disp; std::vector< SliceType::Pointer > compartmentSlices; std::vector< double > t2Vector; for (int i=0; i* signalModel; if (iSetLargestPossibleRegion( sliceRegion ); slice->SetBufferedRegion( sliceRegion ); slice->SetRequestedRegion( sliceRegion ); slice->Allocate(); slice->FillBuffer(0.0); // extract slice from channel g for (int y=0; yGetLargestPossibleRegion().GetSize(1); y++) for (int x=0; xGetLargestPossibleRegion().GetSize(0); x++) { SliceType::IndexType index2D; index2D[0]=x; index2D[1]=y; DoubleDwiType::IndexType index3D; index3D[0]=x; index3D[1]=y; index3D[2]=z; slice->SetPixel(index2D, images.at(i)->GetPixel(index3D)[g]); if (fMap.IsNotNull() && i==0) fMap->SetPixel(index2D, m_FrequencyMap->GetPixel(index3D)); } compartmentSlices.push_back(slice); t2Vector.push_back(signalModel->GetT2()); } // create k-sapce (inverse fourier transform slices) + itk::Size<2> outSize; outSize.SetElement(0, m_ImageRegion.GetSize(0)); outSize.SetElement(1, m_ImageRegion.GetSize(1)); itk::KspaceImageFilter< SliceType::PixelType >::Pointer idft = itk::KspaceImageFilter< SliceType::PixelType >::New(); idft->SetCompartmentImages(compartmentSlices); idft->SetT2(t2Vector); idft->SetkOffset(m_kOffset); idft->SettLine(m_tLine); idft->SetTE(m_TE); idft->SetTinhom(m_tInhom); idft->SetSimulateRelaxation(m_SimulateRelaxation); idft->SetSimulateEddyCurrents(m_SimulateEddyCurrents); idft->SetEddyGradientMagnitude(m_EddyGradientStrength); idft->SetZ((double)z-(double)images.at(0)->GetLargestPossibleRegion().GetSize(2)/2.0); idft->SetDirectionMatrix(transform); idft->SetDiffusionGradientDirection(m_FiberModels.at(0)->GetGradientDirection(g)); idft->SetFrequencyMap(fMap); idft->SetSignalScale(m_SignalScale); + idft->SetOutSize(outSize); idft->Update(); ComplexSliceType::Pointer fSlice; fSlice = idft->GetOutput(); - fSlice = RearrangeSlice(fSlice); - // add artifacts - for (int a=0; aAddArtifact(fSlice); + for (int i=0; iAddArtifact(fSlice); // fourier transform slice SliceType::Pointer newSlice; itk::DftImageFilter< SliceType::PixelType >::Pointer dft = itk::DftImageFilter< SliceType::PixelType >::New(); dft->SetInput(fSlice); dft->Update(); newSlice = dft->GetOutput(); // put slice back into channel g for (int y=0; yGetLargestPossibleRegion().GetSize(1); y++) for (int x=0; xGetLargestPossibleRegion().GetSize(0); x++) { DoubleDwiType::IndexType index3D; index3D[0]=x; index3D[1]=y; index3D[2]=z; SliceType::IndexType index2D; index2D[0]=x; index2D[1]=y; DoubleDwiType::PixelType pix3D = newImage->GetPixel(index3D); pix3D[g] = newSlice->GetPixel(index2D); newImage->SetPixel(index3D, pix3D); } + ++disp; } return newImage; } TractsToDWIImageFilter::ComplexSliceType::Pointer TractsToDWIImageFilter::RearrangeSlice(ComplexSliceType::Pointer slice) { ImageRegion<2> region = slice->GetLargestPossibleRegion(); ComplexSliceType::Pointer rearrangedSlice = ComplexSliceType::New(); rearrangedSlice->SetLargestPossibleRegion( region ); rearrangedSlice->SetBufferedRegion( region ); rearrangedSlice->SetRequestedRegion( region ); rearrangedSlice->Allocate(); int xHalf = region.GetSize(0)/2; int yHalf = region.GetSize(1)/2; for (int y=0; y pix = slice->GetPixel(idx); if( idx[0] < xHalf ) idx[0] = idx[0] + xHalf; else idx[0] = idx[0] - xHalf; if( idx[1] < yHalf ) idx[1] = idx[1] + yHalf; else idx[1] = idx[1] - yHalf; rearrangedSlice->SetPixel(idx, pix); } return rearrangedSlice; } void TractsToDWIImageFilter::GenerateData() { // check input data if (m_FiberBundle.IsNull()) itkExceptionMacro("Input fiber bundle is NULL!"); int numFibers = m_FiberBundle->GetNumFibers(); if (numFibers<=0) itkExceptionMacro("Input fiber bundle contains no fibers!"); if (m_FiberModels.empty()) itkExceptionMacro("No diffusion model for fiber compartments defined!"); if (m_NonFiberModels.empty()) itkExceptionMacro("No diffusion model for non-fiber compartments defined!"); int baselineIndex = m_FiberModels[0]->GetFirstBaselineIndex(); if (baselineIndex<0) itkExceptionMacro("No baseline index found!"); - // determine k-space undersampling - for (int i=0; i*>(m_KspaceArtifacts.at(i)) ) - m_Upsampling = dynamic_cast*>(m_KspaceArtifacts.at(i))->GetKspaceCropping(); + // check k-space undersampling if (m_Upsampling<1) m_Upsampling = 1; if (m_TissueMask.IsNotNull()) { // use input tissue mask m_Spacing = m_TissueMask->GetSpacing(); m_Origin = m_TissueMask->GetOrigin(); m_DirectionMatrix = m_TissueMask->GetDirection(); m_ImageRegion = m_TissueMask->GetLargestPossibleRegion(); - if (m_Upsampling>1) + if (m_Upsampling>1.00001) { + MITK_INFO << "Adding ringing artifacts. Image upsampling factor: " << m_Upsampling; ImageRegion<3> region = m_ImageRegion; region.SetSize(0, m_ImageRegion.GetSize(0)*m_Upsampling); region.SetSize(1, m_ImageRegion.GetSize(1)*m_Upsampling); itk::Vector spacing = m_Spacing; spacing[0] /= m_Upsampling; spacing[1] /= m_Upsampling; itk::RescaleIntensityImageFilter::Pointer rescaler = itk::RescaleIntensityImageFilter::New(); rescaler->SetInput(0,m_TissueMask); rescaler->SetOutputMaximum(100); rescaler->SetOutputMinimum(0); rescaler->Update(); itk::ResampleImageFilter::Pointer resampler = itk::ResampleImageFilter::New(); resampler->SetInput(rescaler->GetOutput()); resampler->SetOutputParametersFromImage(m_TissueMask); resampler->SetSize(region.GetSize()); resampler->SetOutputSpacing(spacing); resampler->Update(); m_TissueMask = resampler->GetOutput(); } MITK_INFO << "Using tissue mask"; } // initialize output dwi image OutputImageType::Pointer outImage = OutputImageType::New(); outImage->SetSpacing( m_Spacing ); outImage->SetOrigin( m_Origin ); outImage->SetDirection( m_DirectionMatrix ); outImage->SetLargestPossibleRegion( m_ImageRegion ); outImage->SetBufferedRegion( m_ImageRegion ); outImage->SetRequestedRegion( m_ImageRegion ); outImage->SetVectorLength( m_FiberModels[0]->GetNumGradients() ); outImage->Allocate(); OutputImageType::PixelType temp; temp.SetSize(m_FiberModels[0]->GetNumGradients()); temp.Fill(0.0); outImage->FillBuffer(temp); // is input slize size a power of two? int x=m_ImageRegion.GetSize(0); int y=m_ImageRegion.GetSize(1); if ( x%2 == 1 ) x += 1; if ( y%2 == 1 ) y += 1; // if not, adjust size and dimension (needed for FFT); zero-padding if (x!=m_ImageRegion.GetSize(0)) - { - MITK_INFO << "Adjusting image width: " << m_ImageRegion.GetSize(0) << " --> " << x << " --> " << x*m_Upsampling; m_ImageRegion.SetSize(0, x); - } if (y!=m_ImageRegion.GetSize(1)) - { - MITK_INFO << "Adjusting image height: " << m_ImageRegion.GetSize(1) << " --> " << y << " --> " << y*m_Upsampling; m_ImageRegion.SetSize(1, y); - } // apply undersampling to image parameters m_UpsampledSpacing = m_Spacing; m_UpsampledImageRegion = m_ImageRegion; m_UpsampledSpacing[0] /= m_Upsampling; m_UpsampledSpacing[1] /= m_Upsampling; m_UpsampledImageRegion.SetSize(0, m_ImageRegion.GetSize()[0]*m_Upsampling); m_UpsampledImageRegion.SetSize(1, m_ImageRegion.GetSize()[1]*m_Upsampling); // everything from here on is using the upsampled image parameters!!! if (m_TissueMask.IsNull()) { m_TissueMask = ItkUcharImgType::New(); m_TissueMask->SetSpacing( m_UpsampledSpacing ); m_TissueMask->SetOrigin( m_Origin ); m_TissueMask->SetDirection( m_DirectionMatrix ); m_TissueMask->SetLargestPossibleRegion( m_UpsampledImageRegion ); m_TissueMask->SetBufferedRegion( m_UpsampledImageRegion ); m_TissueMask->SetRequestedRegion( m_UpsampledImageRegion ); m_TissueMask->Allocate(); m_TissueMask->FillBuffer(1); } // resample frequency map if (m_FrequencyMap.IsNotNull()) { itk::ResampleImageFilter::Pointer resampler = itk::ResampleImageFilter::New(); resampler->SetInput(m_FrequencyMap); resampler->SetOutputParametersFromImage(m_FrequencyMap); resampler->SetSize(m_UpsampledImageRegion.GetSize()); resampler->SetOutputSpacing(m_UpsampledSpacing); resampler->Update(); m_FrequencyMap = resampler->GetOutput(); } // initialize volume fraction images m_VolumeFractions.clear(); for (int i=0; iSetSpacing( m_UpsampledSpacing ); tempimg->SetOrigin( m_Origin ); tempimg->SetDirection( m_DirectionMatrix ); tempimg->SetLargestPossibleRegion( m_UpsampledImageRegion ); tempimg->SetBufferedRegion( m_UpsampledImageRegion ); tempimg->SetRequestedRegion( m_UpsampledImageRegion ); tempimg->Allocate(); tempimg->FillBuffer(0); m_VolumeFractions.push_back(tempimg); } // resample fiber bundle for sufficient voxel coverage double segmentVolume = 0.0001; float minSpacing = 1; if(m_UpsampledSpacing[0]GetDeepCopy(); fiberBundle->ResampleFibers(minSpacing/m_VolumeAccuracy); double mmRadius = m_FiberRadius/1000; if (mmRadius>0) segmentVolume = M_PI*mmRadius*mmRadius*minSpacing/m_VolumeAccuracy; // generate double images to work with because we don't want to lose precision // we use a separate image for each compartment model std::vector< DoubleDwiType::Pointer > compartments; for (int i=0; iSetSpacing( m_UpsampledSpacing ); doubleDwi->SetOrigin( m_Origin ); doubleDwi->SetDirection( m_DirectionMatrix ); doubleDwi->SetLargestPossibleRegion( m_UpsampledImageRegion ); doubleDwi->SetBufferedRegion( m_UpsampledImageRegion ); doubleDwi->SetRequestedRegion( m_UpsampledImageRegion ); doubleDwi->SetVectorLength( m_FiberModels[0]->GetNumGradients() ); doubleDwi->Allocate(); DoubleDwiType::PixelType pix; pix.SetSize(m_FiberModels[0]->GetNumGradients()); pix.Fill(0.0); doubleDwi->FillBuffer(pix); compartments.push_back(doubleDwi); } double interpFact = 2*atan(-0.5*m_InterpolationShrink); double maxVolume = 0; - vtkSmartPointer fiberPolyData = fiberBundle->GetFiberPolyData(); - vtkSmartPointer vLines = fiberPolyData->GetLines(); - vLines->InitTraversal(); - MITK_INFO << "Generating signal of " << m_FiberModels.size() << " fiber compartments"; + vtkSmartPointer fiberPolyData = fiberBundle->GetFiberPolyData(); boost::progress_display disp(numFibers); for( int i=0; iGetNextCell ( numPoints, points ); + vtkCell* cell = fiberPolyData->GetCell(i); + int numPoints = cell->GetNumberOfPoints(); + vtkPoints* points = cell->GetPoints(); + if (numPoints<2) continue; for( int j=0; jGetPoint(points[j]); + double* temp = points->GetPoint(j); itk::Point vertex = GetItkPoint(temp); itk::Vector v = GetItkVector(temp); itk::Vector dir(3); if (jGetPoint(points[j+1]))-v; + dir = GetItkVector(points->GetPoint(j+1))-v; else - dir = v-GetItkVector(fiberPolyData->GetPoint(points[j-1])); + dir = v-GetItkVector(points->GetPoint(j-1)); itk::Index<3> idx; itk::ContinuousIndex contIndex; m_TissueMask->TransformPhysicalPointToIndex(vertex, idx); m_TissueMask->TransformPhysicalPointToContinuousIndex(vertex, contIndex); if (!m_UseInterpolation) // use nearest neighbour interpolation { if (!m_TissueMask->GetLargestPossibleRegion().IsInside(idx) || m_TissueMask->GetPixel(idx)<=0) continue; // generate signal for each fiber compartment for (int k=0; kSetFiberDirection(dir); DoubleDwiType::PixelType pix = doubleDwi->GetPixel(idx); pix += segmentVolume*m_FiberModels[k]->SimulateMeasurement(); doubleDwi->SetPixel(idx, pix ); if (pix[baselineIndex]>maxVolume) maxVolume = pix[baselineIndex]; } continue; } double frac_x = contIndex[0] - idx[0]; double frac_y = contIndex[1] - idx[1]; double frac_z = contIndex[2] - idx[2]; if (frac_x<0) { idx[0] -= 1; frac_x += 1; } if (frac_y<0) { idx[1] -= 1; frac_y += 1; } if (frac_z<0) { idx[2] -= 1; frac_z += 1; } frac_x = atan((0.5-frac_x)*m_InterpolationShrink)/interpFact + 0.5; frac_y = atan((0.5-frac_y)*m_InterpolationShrink)/interpFact + 0.5; frac_z = atan((0.5-frac_z)*m_InterpolationShrink)/interpFact + 0.5; // use trilinear interpolation itk::Index<3> newIdx; for (int x=0; x<2; x++) { frac_x = 1-frac_x; for (int y=0; y<2; y++) { frac_y = 1-frac_y; for (int z=0; z<2; z++) { frac_z = 1-frac_z; newIdx[0] = idx[0]+x; newIdx[1] = idx[1]+y; newIdx[2] = idx[2]+z; double frac = frac_x*frac_y*frac_z; // is position valid? if (!m_TissueMask->GetLargestPossibleRegion().IsInside(newIdx) || m_TissueMask->GetPixel(newIdx)<=0) continue; // generate signal for each fiber compartment for (int k=0; kSetFiberDirection(dir); DoubleDwiType::PixelType pix = doubleDwi->GetPixel(newIdx); pix += segmentVolume*frac*m_FiberModels[k]->SimulateMeasurement(); doubleDwi->SetPixel(newIdx, pix ); if (pix[baselineIndex]>maxVolume) maxVolume = pix[baselineIndex]; } } } } } + ++disp; } MITK_INFO << "Generating signal of " << m_NonFiberModels.size() << " non-fiber compartments"; ImageRegionIterator it3(m_TissueMask, m_TissueMask->GetLargestPossibleRegion()); boost::progress_display disp3(m_TissueMask->GetLargestPossibleRegion().GetNumberOfPixels()); double voxelVolume = m_UpsampledSpacing[0]*m_UpsampledSpacing[1]*m_UpsampledSpacing[2]; double fact = 1; if (m_FiberRadius<0.0001) fact = voxelVolume/maxVolume; while(!it3.IsAtEnd()) { - ++disp3; DoubleDwiType::IndexType index = it3.GetIndex(); if (it3.Get()>0) { // get fiber volume fraction DoubleDwiType::Pointer fiberDwi = compartments.at(0); DoubleDwiType::PixelType fiberPix = fiberDwi->GetPixel(index); // intra axonal compartment if (fact>1) // auto scale intra-axonal if no fiber radius is specified { fiberPix *= fact; fiberDwi->SetPixel(index, fiberPix); } double f = fiberPix[baselineIndex]; if (f>voxelVolume || f>0 && m_EnforcePureFiberVoxels) // more fiber than space in voxel? { fiberDwi->SetPixel(index, fiberPix*voxelVolume/f); for (int i=1; iSetPixel(index, pix); m_VolumeFractions.at(i)->SetPixel(index, 1); } } else { m_VolumeFractions.at(0)->SetPixel(index, f); double nonf = voxelVolume-f; // non-fiber volume double inter = 0; if (m_FiberModels.size()>1) inter = nonf * f/voxelVolume; // intra-axonal fraction of non fiber compartment scales linearly with f double other = nonf - inter; // rest of compartment double singleinter = inter/(m_FiberModels.size()-1); // adjust non-fiber and intra-axonal signal for (int i=1; iGetPixel(index); if (pix[baselineIndex]>0) pix /= pix[baselineIndex]; pix *= singleinter; doubleDwi->SetPixel(index, pix); m_VolumeFractions.at(i)->SetPixel(index, singleinter/voxelVolume); } for (int i=0; iGetPixel(index) + m_NonFiberModels[i]->SimulateMeasurement()*other*m_NonFiberModels[i]->GetWeight(); doubleDwi->SetPixel(index, pix); m_VolumeFractions.at(i+m_FiberModels.size())->SetPixel(index, other/voxelVolume*m_NonFiberModels[i]->GetWeight()); } } } ++it3; + ++disp3; } // do k-space stuff DoubleDwiType::Pointer doubleOutImage; - if (m_FrequencyMap.IsNotNull() || !m_KspaceArtifacts.empty() || m_kOffset>0 || m_SimulateRelaxation || m_SimulateEddyCurrents) + if (m_FrequencyMap.IsNotNull() || !m_KspaceArtifacts.empty() || m_kOffset>0 || m_SimulateRelaxation || m_SimulateEddyCurrents || m_Upsampling>1.00001) { MITK_INFO << "Adjusting complex signal"; doubleOutImage = DoKspaceStuff(compartments); m_SignalScale = 1; } else { MITK_INFO << "Summing compartments"; doubleOutImage = compartments.at(0); for (int i=1; i::Pointer adder = itk::AddImageFilter< DoubleDwiType, DoubleDwiType, DoubleDwiType>::New(); adder->SetInput1(doubleOutImage); adder->SetInput2(compartments.at(i)); adder->Update(); doubleOutImage = adder->GetOutput(); } } MITK_INFO << "Finalizing image"; unsigned int window = 0; unsigned int min = itk::NumericTraits::max(); ImageRegionIterator it4 (outImage, outImage->GetLargestPossibleRegion()); DoubleDwiType::PixelType signal; signal.SetSize(m_FiberModels[0]->GetNumGradients()); boost::progress_display disp4(outImage->GetLargestPossibleRegion().GetNumberOfPixels()); while(!it4.IsAtEnd()) { ++disp4; DWIImageType::IndexType index = it4.GetIndex(); signal = doubleOutImage->GetPixel(index)*m_SignalScale; if (m_NoiseModel->GetNoiseVariance() > 0) { DoubleDwiType::PixelType accu = signal; accu.Fill(0.0); for (int i=0; iAddNoise(temp); accu += temp; } signal = accu/m_NumberOfRepetitions; } for (int i=0; i0) signal[i] = floor(signal[i]+0.5); else signal[i] = ceil(signal[i]-0.5); if (!m_FiberModels.at(0)->IsBaselineIndex(i) && signal[i]>window) window = signal[i]; if (!m_FiberModels.at(0)->IsBaselineIndex(i) && signal[i]SetNthOutput(0, outImage); } itk::Point TractsToDWIImageFilter::GetItkPoint(double point[3]) { itk::Point itkPoint; itkPoint[0] = point[0]; itkPoint[1] = point[1]; itkPoint[2] = point[2]; return itkPoint; } itk::Vector TractsToDWIImageFilter::GetItkVector(double point[3]) { itk::Vector itkVector; itkVector[0] = point[0]; itkVector[1] = point[1]; itkVector[2] = point[2]; return itkVector; } vnl_vector_fixed TractsToDWIImageFilter::GetVnlVector(double point[3]) { vnl_vector_fixed vnlVector; vnlVector[0] = point[0]; vnlVector[1] = point[1]; vnlVector[2] = point[2]; return vnlVector; } vnl_vector_fixed TractsToDWIImageFilter::GetVnlVector(Vector& vector) { vnl_vector_fixed vnlVector; vnlVector[0] = vector[0]; vnlVector[1] = vector[1]; vnlVector[2] = vector[2]; return vnlVector; } } diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToDWIImageFilter.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToDWIImageFilter.h old mode 100644 new mode 100755 index d63ea9f277..c09dfbc671 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToDWIImageFilter.h +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToDWIImageFilter.h @@ -1,151 +1,151 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef __itkTractsToDWIImageFilter_h__ #define __itkTractsToDWIImageFilter_h__ // MITK #include #include #include #include -#include // ITK #include #include #include #include #include #include typedef itk::VectorImage< short, 3 > DWIImageType; namespace itk { /** * \brief Generates artificial diffusion weighted image volume from the input fiberbundle using a generic multicompartment model. */ class TractsToDWIImageFilter : public ImageSource< DWIImageType > { public: typedef TractsToDWIImageFilter Self; typedef ImageSource< DWIImageType > Superclass; typedef SmartPointer< Self > Pointer; typedef SmartPointer< const Self > ConstPointer; typedef itk::Image ItkDoubleImgType; typedef itk::Image ItkFloatImgType; typedef itk::Image ItkUcharImgType; typedef mitk::FiberBundleX::Pointer FiberBundleType; typedef itk::VectorImage< double, 3 > DoubleDwiType; typedef std::vector< mitk::KspaceArtifact* > KspaceArtifactList; typedef std::vector< mitk::DiffusionSignalModel* > DiffusionModelList; typedef itk::Matrix MatrixType; typedef mitk::DiffusionNoiseModel NoiseModelType; typedef itk::Image< double, 2 > SliceType; typedef itk::VnlForwardFFTImageFilter::OutputImageType ComplexSliceType; itkNewMacro(Self) itkTypeMacro( TractsToDWIImageFilter, ImageToImageFilter ) // input itkSetMacro( SignalScale, double ) itkSetMacro( FiberRadius, double ) itkSetMacro( InterpolationShrink, double ) ///< large values shrink (towards nearest neighbour interpolation), small values strech interpolation function (towards linear interpolation) itkSetMacro( VolumeAccuracy, unsigned int ) ///< determines fiber sampling density and thereby the accuracy of the fiber volume fraction itkSetMacro( FiberBundle, FiberBundleType ) ///< input fiber bundle itkSetMacro( Spacing, mitk::Vector3D ) ///< output image spacing itkSetMacro( Origin, mitk::Point3D ) ///< output image origin itkSetMacro( DirectionMatrix, MatrixType ) ///< output image rotation itkSetMacro( EnforcePureFiberVoxels, bool ) ///< treat all voxels containing at least one fiber as fiber-only (actually disable non-fiber compartments for this voxel). itkSetMacro( ImageRegion, ImageRegion<3> ) ///< output image size itkSetMacro( NumberOfRepetitions, unsigned int ) ///< number of acquisition repetitions to reduce noise (default is no additional repetition) itkSetMacro( TissueMask, ItkUcharImgType::Pointer ) ///< voxels outside of this binary mask contain only noise (are treated as air) void SetNoiseModel(NoiseModelType* noiseModel){ m_NoiseModel = noiseModel; } ///< generates the noise added to the image values void SetFiberModels(DiffusionModelList modelList){ m_FiberModels = modelList; } ///< generate signal of fiber compartments void SetNonFiberModels(DiffusionModelList modelList){ m_NonFiberModels = modelList; } ///< generate signal of non-fiber compartments void SetKspaceArtifacts(KspaceArtifactList artifactList){ m_KspaceArtifacts = artifactList; } mitk::LevelWindow GetLevelWindow(){ return m_LevelWindow; } itkSetMacro( FrequencyMap, ItkDoubleImgType::Pointer ) itkSetMacro( kOffset, double ) itkSetMacro( tLine, double ) itkSetMacro( tInhom, double ) itkSetMacro( TE, double ) itkSetMacro( UseInterpolation, bool ) itkSetMacro( SimulateEddyCurrents, bool ) itkSetMacro( SimulateRelaxation, bool ) itkSetMacro( EddyGradientStrength, double ) + itkSetMacro( Upsampling, double ) // output std::vector< ItkDoubleImgType::Pointer > GetVolumeFractions(){ return m_VolumeFractions; } void GenerateData(); protected: TractsToDWIImageFilter(); virtual ~TractsToDWIImageFilter(); itk::Point GetItkPoint(double point[3]); itk::Vector GetItkVector(double point[3]); vnl_vector_fixed GetVnlVector(double point[3]); vnl_vector_fixed GetVnlVector(Vector< float, 3 >& vector); /** Transform generated image compartment by compartment, channel by channel and slice by slice using FFT and add k-space artifacts. */ DoubleDwiType::Pointer DoKspaceStuff(std::vector< DoubleDwiType::Pointer >& images); /** Rearrange FFT output to shift low frequencies to the iamge center (correct itk). */ TractsToDWIImageFilter::ComplexSliceType::Pointer RearrangeSlice(ComplexSliceType::Pointer slice); itk::Vector m_Spacing; ///< output image spacing itk::Vector m_UpsampledSpacing; mitk::Point3D m_Origin; ///< output image origin MatrixType m_DirectionMatrix; ///< output image rotation ImageRegion<3> m_ImageRegion; ///< output image size ImageRegion<3> m_UpsampledImageRegion; ItkUcharImgType::Pointer m_TissueMask; ///< voxels outside of this binary mask contain only noise (are treated as air) ItkDoubleImgType::Pointer m_FrequencyMap; ///< map of the B0 inhomogeneities double m_kOffset; double m_tLine; double m_TE; double m_tInhom; FiberBundleType m_FiberBundle; ///< input fiber bundle DiffusionModelList m_FiberModels; ///< generate signal of fiber compartments DiffusionModelList m_NonFiberModels; ///< generate signal of non-fiber compartments KspaceArtifactList m_KspaceArtifacts; NoiseModelType* m_NoiseModel; ///< generates the noise added to the image values bool m_CircleDummy; unsigned int m_VolumeAccuracy; - unsigned int m_Upsampling; + double m_Upsampling; ///< causes ringing artifacts unsigned int m_NumberOfRepetitions; bool m_EnforcePureFiberVoxels; double m_InterpolationShrink; double m_FiberRadius; double m_SignalScale; mitk::LevelWindow m_LevelWindow; bool m_UseInterpolation; std::vector< ItkDoubleImgType::Pointer > m_VolumeFractions; ///< one double image for each compartment containing the corresponding volume fraction per voxel bool m_SimulateRelaxation; bool m_SimulateEddyCurrents; double m_EddyGradientStrength; }; } #include "itkTractsToDWIImageFilter.cpp" #endif diff --git a/Modules/DiffusionImaging/FiberTracking/CMakeLists.txt b/Modules/DiffusionImaging/FiberTracking/CMakeLists.txt index 833650d9e1..0b885f1b05 100644 --- a/Modules/DiffusionImaging/FiberTracking/CMakeLists.txt +++ b/Modules/DiffusionImaging/FiberTracking/CMakeLists.txt @@ -1,44 +1,45 @@ MITK_CHECK_MODULE(_missing_deps DiffusionCore MitkGraphAlgorithms) if(NOT _missing_deps) set(lut_url http://mitk.org/download/data/FibertrackingLUT.tar.gz) set(lut_tarball ${CMAKE_CURRENT_BINARY_DIR}/FibertrackingLUT.tar.gz) message("Downloading FiberTracking LUT ${lut_url}...") file(DOWNLOAD ${lut_url} ${lut_tarball} EXPECTED_MD5 38ecb6d4a826c9ebb0f4965eb9aeee44 TIMEOUT 20 STATUS status SHOW_PROGRESS ) list(GET status 0 status_code) list(GET status 1 status_msg) if(NOT status_code EQUAL 0) message(SEND_ERROR "${status_msg} (error code ${status_code})") else() message("done.") endif() file(MAKE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/Resources) message("Unpacking FiberTracking LUT tarball...") execute_process(COMMAND ${CMAKE_COMMAND} -E tar xzf ../FibertrackingLUT.tar.gz WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/Resources RESULT_VARIABLE result ERROR_VARIABLE err_msg) if(result) message(SEND_ERROR "Unpacking FibertrackingLUT.tar.gz failed: ${err_msg}") else() message("done.") endif() endif() MITK_CREATE_MODULE( FiberTracking SUBPROJECTS MITK-DTI INCLUDE_DIRS Algorithms Algorithms/GibbsTracking Algorithms/StochasticTracking IODataStructures IODataStructures/FiberBundleX IODataStructures/PlanarFigureComposite Interactions SignalModels Rendering ${CMAKE_CURRENT_BINARY_DIR} DEPENDS DiffusionCore MitkGraphAlgorithms ) if(MODULE_IS_ENABLED) add_subdirectory(Testing) + add_subdirectory(MiniApps) endif() diff --git a/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundleX/mitkFiberBundleX.cpp b/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundleX/mitkFiberBundleX.cpp old mode 100644 new mode 100755 index 88ad9c95f0..70bf2f97ec --- a/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundleX/mitkFiberBundleX.cpp +++ b/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundleX/mitkFiberBundleX.cpp @@ -1,1703 +1,1703 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #define _USE_MATH_DEFINES #include "mitkFiberBundleX.h" #include #include #include #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"; using namespace std; mitk::FiberBundleX::FiberBundleX( vtkPolyData* fiberPolyData ) : m_CurrentColorCoding(NULL) , m_NumFibers(0) , m_FiberSampling(0) { m_FiberPolyData = vtkSmartPointer::New(); if (fiberPolyData != NULL) { m_FiberPolyData = fiberPolyData; //m_FiberPolyData->DeepCopy(fiberPolyData); this->DoColorCodingOrientationBased(); } this->UpdateFiberGeometry(); this->SetColorCoding(COLORCODING_ORIENTATION_BASED); this->GenerateFiberIds(); } mitk::FiberBundleX::~FiberBundleX() { } mitk::FiberBundleX::Pointer mitk::FiberBundleX::GetDeepCopy() { mitk::FiberBundleX::Pointer newFib = mitk::FiberBundleX::New(m_FiberPolyData); 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(); } 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 || *finIt>GetNumFibers()){ MITK_INFO << "FiberID can not be negative or >NumFibers!!! check id Extraction!" << *finIt; break; } vtkSmartPointer fiber = m_FiberIdDataSet->GetCell(*finIt);//->DeepCopy(fiber); vtkSmartPointer fibPoints = fiber->GetPoints(); vtkSmartPointer newFiber = vtkSmartPointer::New(); newFiber->GetPointIds()->SetNumberOfIds( fibPoints->GetNumberOfPoints() ); 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_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]; } } newLineSet->InsertNextCell(newFiber); ++finIt; } 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::AddBundle(mitk::FiberBundleX* fib) { if (fib==NULL) { MITK_WARN << "trying to call AddBundle with NULL argument"; return NULL; } MITK_INFO << "Adding fibers"; vtkSmartPointer vNewPolyData = vtkSmartPointer::New(); vtkSmartPointer vNewLines = vtkSmartPointer::New(); vtkSmartPointer vNewPoints = vtkSmartPointer::New(); // add current fiber bundle for (int i=0; iGetNumberOfCells(); i++) { vtkCell* cell = m_FiberPolyData->GetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(j, p); vtkIdType id = vNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } vNewLines->InsertNextCell(container); } // add new fiber bundle for (int i=0; iGetFiberPolyData()->GetNumberOfCells(); i++) { vtkCell* cell = fib->GetFiberPolyData()->GetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(j, p); vtkIdType id = vNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } vNewLines->InsertNextCell(container); } // initialize polydata vNewPolyData->SetPoints(vNewPoints); vNewPolyData->SetLines(vNewLines); // initialize fiber bundle mitk::FiberBundleX::Pointer newFib = mitk::FiberBundleX::New(vNewPolyData); return newFib; } // subtract two fiber bundles mitk::FiberBundleX::Pointer mitk::FiberBundleX::SubtractBundle(mitk::FiberBundleX* fib) { MITK_INFO << "Subtracting fibers"; vtkSmartPointer vNewPolyData = vtkSmartPointer::New(); vtkSmartPointer vNewLines = vtkSmartPointer::New(); vtkSmartPointer vNewPoints = vtkSmartPointer::New(); // iterate over current fibers int numFibers = GetNumFibers(); boost::progress_display disp(numFibers); for( int i=0; iGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); if (points==NULL || numPoints<=0) continue; int numFibers2 = fib->GetNumFibers(); bool contained = false; for( int i2=0; i2GetFiberPolyData()->GetCell(i2); int numPoints2 = cell2->GetNumberOfPoints(); vtkPoints* points2 = cell2->GetPoints(); if (points2==NULL || numPoints2<=0) continue; // check endpoints itk::Point point_start = GetItkPoint(points->GetPoint(0)); itk::Point point_end = GetItkPoint(points->GetPoint(numPoints-1)); itk::Point point2_start = GetItkPoint(points2->GetPoint(0)); itk::Point point2_end = GetItkPoint(points2->GetPoint(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 ??? if (numPoints2==numPoints) contained = true; } } // add to result because fiber is not subtracted if (!contained) { vtkSmartPointer container = vtkSmartPointer::New(); for( int j=0; jInsertNextPoint(points->GetPoint(j)); container->GetPointIds()->InsertNextId(id); } vNewLines->InsertNextCell(container); } } if(vNewLines->GetNumberOfCells()==0) return NULL; // initialize polydata vNewPolyData->SetPoints(vNewPoints); vNewPolyData->SetLines(vNewLines); // 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; } /* * 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(); } m_NumFibers = m_FiberPolyData->GetNumberOfLines(); if (updateGeometry) UpdateFiberGeometry(); SetColorCoding(COLORCODING_ORIENTATION_BASED); GenerateFiberIds(); } /* * return vtkPolyData */ vtkSmartPointer mitk::FiberBundleX::GetFiberPolyData() { 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->ResetFiberOpacity(); 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(); //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 = vtkSmartPointer::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; } /* extract single fibers of fiberBundle */ vtkCellArray* fiberList = m_FiberPolyData->GetLines(); fiberList->InitTraversal(); for (int fi=0; fiGetNextCell(pointsPerFiber, idList); // 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 0) { /* The color value of the current point is influenced by the previous point and next point. */ vnl_vector_fixed< double, 3 > currentPntvtk(extrPoints->GetPoint(idList[i])[0], extrPoints->GetPoint(idList[i])[1],extrPoints->GetPoint(idList[i])[2]); vnl_vector_fixed< double, 3 > nextPntvtk(extrPoints->GetPoint(idList[i+1])[0], extrPoints->GetPoint(idList[i+1])[1], extrPoints->GetPoint(idList[i+1])[2]); vnl_vector_fixed< double, 3 > prevPntvtk(extrPoints->GetPoint(idList[i-1])[0], extrPoints->GetPoint(idList[i-1])[1], extrPoints->GetPoint(idList[i-1])[2]); vnl_vector_fixed< double, 3 > diff1; diff1 = currentPntvtk - nextPntvtk; vnl_vector_fixed< double, 3 > diff2; diff2 = currentPntvtk - prevPntvtk; vnl_vector_fixed< double, 3 > diff; diff = (diff1 - diff2) / 2.0; diff.normalize(); rgba[0] = (unsigned char) (255.0 * std::fabs(diff[0])); rgba[1] = (unsigned char) (255.0 * std::fabs(diff[1])); rgba[2] = (unsigned char) (255.0 * std::fabs(diff[2])); rgba[3] = (unsigned char) (255.0); } else if (i==0) { /* First point has no previous point, therefore only diff1 is taken */ vnl_vector_fixed< double, 3 > currentPntvtk(extrPoints->GetPoint(idList[i])[0], extrPoints->GetPoint(idList[i])[1],extrPoints->GetPoint(idList[i])[2]); vnl_vector_fixed< double, 3 > nextPntvtk(extrPoints->GetPoint(idList[i+1])[0], extrPoints->GetPoint(idList[i+1])[1], extrPoints->GetPoint(idList[i+1])[2]); vnl_vector_fixed< double, 3 > diff1; diff1 = currentPntvtk - nextPntvtk; diff1.normalize(); rgba[0] = (unsigned char) (255.0 * std::fabs(diff1[0])); rgba[1] = (unsigned char) (255.0 * std::fabs(diff1[1])); rgba[2] = (unsigned char) (255.0 * std::fabs(diff1[2])); rgba[3] = (unsigned char) (255.0); } else if (i==pointsPerFiber-1) { /* Last point has no next point, therefore only diff2 is taken */ vnl_vector_fixed< double, 3 > currentPntvtk(extrPoints->GetPoint(idList[i])[0], extrPoints->GetPoint(idList[i])[1],extrPoints->GetPoint(idList[i])[2]); vnl_vector_fixed< double, 3 > prevPntvtk(extrPoints->GetPoint(idList[i-1])[0], extrPoints->GetPoint(idList[i-1])[1], extrPoints->GetPoint(idList[i-1])[2]); vnl_vector_fixed< double, 3 > diff2; diff2 = currentPntvtk - prevPntvtk; diff2.normalize(); rgba[0] = (unsigned char) (255.0 * std::fabs(diff2[0])); rgba[1] = (unsigned char) (255.0 * std::fabs(diff2[1])); rgba[2] = (unsigned char) (255.0 * std::fabs(diff2[2])); rgba[3] = (unsigned char) (255.0); } colorsT->InsertTupleValue(idList[i], rgba); } //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 { MITK_DEBUG << "Fiber with 0 points detected... please check your tractography algorithm!" ; continue; } }//end for loop 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); // =========================== //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; this->SetColorCoding(COLORCODING_FA_BASED); MITK_DEBUG << "FBX: done CC FA based"; this->GenerateFiberIds(); } void mitk::FiberBundleX::DoUseFaFiberOpacity() { if(m_FiberPolyData->GetPointData()->HasArray(COLORCODING_FA_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)); 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->GenerateFiberIds(); } void mitk::FiberBundleX::ResetFiberOpacity() { vtkUnsignedCharArray* ColorArray = dynamic_cast (m_FiberPolyData->GetPointData()->GetArray(COLORCODING_ORIENTATION_BASED)); if (ColorArray==NULL) return; 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 = vtkSmartPointer::New(); faValues->SetName(COLORCODING_FA_BASED); faValues->Allocate(m_FiberPolyData->GetNumberOfPoints()); faValues->SetNumberOfValues(m_FiberPolyData->GetNumberOfPoints()); 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 = 1-FAimage->GetPixelValueByWorldCoordinate(px); faValues->InsertValue(i, faPixelValue); } m_FiberPolyData->GetPointData()->AddArray(faValues); this->GenerateFiberIds(); if(m_FiberPolyData->GetPointData()->HasArray(COLORCODING_FA_BASED)) MITK_DEBUG << "FA VALUE ARRAY SET"; } void mitk::FiberBundleX::GenerateFiberIds() { 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(); m_FiberIdDataSet = idFiberFilter->GetOutput(); MITK_DEBUG << "Generating Fiber Ids...[done] | " << m_FiberIdDataSet->GetNumberOfCells(); } mitk::FiberBundleX::Pointer mitk::FiberBundleX::ExtractFiberSubset(ItkUcharImgType* mask, bool anyPoint) { vtkSmartPointer polyData = m_FiberPolyData; if (anyPoint) { float minSpacing = 1; if(mask->GetSpacing()[0]GetSpacing()[1] && mask->GetSpacing()[0]GetSpacing()[2]) minSpacing = mask->GetSpacing()[0]; else if (mask->GetSpacing()[1] < mask->GetSpacing()[2]) minSpacing = mask->GetSpacing()[1]; else minSpacing = mask->GetSpacing()[2]; mitk::FiberBundleX::Pointer fibCopy = this->GetDeepCopy(); fibCopy->ResampleFibers(minSpacing/10); polyData = fibCopy->GetFiberPolyData(); } vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); MITK_INFO << "Extracting fibers"; boost::progress_display disp(m_NumFibers); for (int i=0; iGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkCell* cellOriginal = m_FiberPolyData->GetCell(i); - int numPointsOriginal = cell->GetNumberOfPoints(); - vtkPoints* pointsOriginal = cell->GetPoints(); + int numPointsOriginal = cellOriginal->GetNumberOfPoints(); + vtkPoints* pointsOriginal = cellOriginal->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); - if (numPoints>1) + if (numPoints>1 && numPointsOriginal) { if (anyPoint) { for (int j=0; jGetPoint(j); itk::Point itkP; itkP[0] = p[0]; itkP[1] = p[1]; itkP[2] = p[2]; itk::Index<3> idx; mask->TransformPhysicalPointToIndex(itkP, idx); - if ( mask->GetPixel(idx)>0 ) + if ( mask->GetPixel(idx)>0 && mask->GetLargestPossibleRegion().IsInside(idx) ) { for (int k=0; kGetPoint(k); vtkIdType id = vtkNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } break; } } } else { - double* start = points->GetPoint(0); + double* start = pointsOriginal->GetPoint(0); itk::Point itkStart; itkStart[0] = start[0]; itkStart[1] = start[1]; itkStart[2] = start[2]; itk::Index<3> idxStart; mask->TransformPhysicalPointToIndex(itkStart, idxStart); - double* end = points->GetPoint(numPoints-1); + double* end = pointsOriginal->GetPoint(numPointsOriginal-1); itk::Point itkEnd; itkEnd[0] = end[0]; itkEnd[1] = end[1]; itkEnd[2] = end[2]; itk::Index<3> idxEnd; mask->TransformPhysicalPointToIndex(itkEnd, idxEnd); - if ( mask->GetPixel(idxStart)>0 && mask->GetPixel(idxEnd)>0 ) + if ( mask->GetPixel(idxStart)>0 && mask->GetPixel(idxEnd)>0 && mask->GetLargestPossibleRegion().IsInside(idxStart) && mask->GetLargestPossibleRegion().IsInside(idxEnd) ) { - for (int j=0; jGetPoint(j); + double* p = pointsOriginal->GetPoint(j); vtkIdType id = vtkNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } } } } vtkNewCells->InsertNextCell(container); } if (vtkNewCells->GetNumberOfCells()<=0) return NULL; vtkSmartPointer newPolyData = vtkSmartPointer::New(); newPolyData->SetPoints(vtkNewPoints); newPolyData->SetLines(vtkNewCells); return mitk::FiberBundleX::New(newPolyData); } mitk::FiberBundleX::Pointer mitk::FiberBundleX::ExtractFiberSubset(mitk::PlanarFigure* pf) { if (pf==NULL) return NULL; std::vector tmp = ExtractFiberIdSubset(pf); if (tmp.size()<=0) return mitk::FiberBundleX::New(); vtkSmartPointer pTmp = GeneratePolyDataByIds(tmp); return mitk::FiberBundleX::New(pTmp); } std::vector mitk::FiberBundleX::ExtractFiberIdSubset(mitk::PlanarFigure* pf) { MITK_DEBUG << "Extracting fibers!"; // vector which is returned, contains all extracted FiberIds std::vector FibersInROI; if (pf==NULL) return FibersInROI; /* Handle type of planarfigure */ // if incoming pf is a pfc mitk::PlanarFigureComposite::Pointer pfcomp= dynamic_cast(pf); 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; } } 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]); /* 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"; 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! * size of list can be optimized by reducing allocation, but be aware of iterator and vector size*/ 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) PointsOnPlane.push_back(i); } MITK_DEBUG << "Num Of points on plane: " << PointsOnPlane.size(); MITK_DEBUG << "Step 2: extract Interesting points with respect to given extraction planarFigure"; 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 double distPF = V1w.EuclideanDistanceTo(V2w); 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])) ; if( XdistPnt <= distPF) PointsInROI.push_back(PointsOnPlane[i]); } } else if ( pf->GetNameOfClass() == polyName->GetNameOfClass() ) { //create vtkPolygon using controlpoints from planarFigure polygon vtkSmartPointer polygonVtk = vtkSmartPointer::New(); //get the control points from pf and insert them to vtkPolygon unsigned int nrCtrlPnts = pf->GetNumberOfControlPoints(); for (int i=0; iGetPoints()->InsertNextPoint((double)pf->GetWorldControlPoint(i)[0], (double)pf->GetWorldControlPoint(i)[1], (double)pf->GetWorldControlPoint(i)[2] ); } //prepare everything for using pointInPolygon function double n[3]; polygonVtk->ComputeNormal(polygonVtk->GetPoints()->GetNumberOfPoints(), static_cast(polygonVtk->GetPoints()->GetData()->GetVoidPointer(0)), n); double bounds[6]; polygonVtk->GetPoints()->GetBounds(bounds); for (int i=0; iGetPoint(PointsOnPlane[i])[0], clipperout->GetPoint(PointsOnPlane[i])[1], clipperout->GetPoint(PointsOnPlane[i])[2]}; int isInPolygon = polygonVtk->PointInPolygon(checkIn, polygonVtk->GetPoints()->GetNumberOfPoints() , static_cast(polygonVtk->GetPoints()->GetData()->GetVoidPointer(0)), bounds, n); if( isInPolygon ) PointsInROI.push_back(PointsOnPlane[i]); } } MITK_DEBUG << "Step3: Identify fibers"; // 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 } if (PointsInROI.size()<=0) return FibersInROI; // 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.resize(numofClippedPoints); //prepare resulting vector FibersInROI.reserve(PointsInROI.size()); 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); // 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] ]; } } 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] ]; if (pointindexFiberMap[ PointsInROI[k] ]<=GetNumFibers() && pointindexFiberMap[ PointsInROI[k] ]>=0) FibersInROI.push_back(pointindexFiberMap[ PointsInROI[k] ]); else MITK_INFO << "ERROR in ExtractFiberIdSubset; impossible fiber id detected"; } m_PointsRoi = PointsInROI; } // detecting fiberId duplicates MITK_DEBUG << "check for duplicates"; 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() { vtkSmartPointer cleaner = vtkSmartPointer::New(); cleaner->SetInput(m_FiberPolyData); cleaner->PointMergingOff(); cleaner->Update(); m_FiberPolyData = cleaner->GetOutput(); m_FiberLengths.clear(); m_MeanFiberLength = 0; m_MedianFiberLength = 0; m_LengthStDev = 0; m_NumFibers = m_FiberPolyData->GetNumberOfCells(); if (m_NumFibers<=0) // no fibers present; apply default geometry { m_MinFiberLength = 0; m_MaxFiberLength = 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::NonpositiveMin(); float max = itk::NumericTraits::max(); float b[] = {max, min, max, min, max, min}; for (int i=0; iGetNumberOfCells(); i++) { vtkCell* cell = m_FiberPolyData->GetCell(i); int p = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); float length = 0; for (int j=0; jGetPoint(j, p1); if (p1[0]b[1]) b[1]=p1[0]; if (p1[1]b[3]) b[3]=p1[1]; if (p1[2]b[5]) b[5]=p1[2]; // calculate statistics if (jGetPoint(j+1, p2); float dist = std::sqrt((p1[0]-p2[0])*(p1[0]-p2[0])+(p1[1]-p2[1])*(p1[1]-p2[1])+(p1[2]-p2[2])*(p1[2]-p2[2])); length += dist; } } m_FiberLengths.push_back(length); m_MeanFiberLength += length; if (i==0) { m_MinFiberLength = length; m_MaxFiberLength = length; } else { if (lengthm_MaxFiberLength) m_MaxFiberLength = length; } } m_MeanFiberLength /= m_NumFibers; std::vector< float > sortedLengths = m_FiberLengths; std::sort(sortedLengths.begin(), sortedLengths.end()); for (int i=0; i1) m_LengthStDev /= (m_NumFibers-1); else m_LengthStDev = 0; m_LengthStDev = std::sqrt(m_LengthStDev); m_MedianFiberLength = sortedLengths.at(m_NumFibers/2); // 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->SetFloatBounds(b); this->SetGeometry(geometry); } QStringList mitk::FiberBundleX::GetAvailableColorCodings() { 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"; return availableColorCodings; } char* mitk::FiberBundleX::GetCurrentColorCoding() { return m_CurrentColorCoding; } void mitk::FiberBundleX::SetColorCoding(const char* requestedColorCoding) { if (requestedColorCoding==NULL) return; MITK_DEBUG << "SetColorCoding:" << requestedColorCoding; 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_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 } } void mitk::FiberBundleX::RotateAroundAxis(double x, double y, double z) { MITK_INFO << "Rotating fibers"; x = x*M_PI/180; y = y*M_PI/180; z = z*M_PI/180; vnl_matrix_fixed< double, 3, 3 > rotX; rotX.set_identity(); rotX[1][1] = cos(x); rotX[2][2] = rotX[1][1]; rotX[1][2] = -sin(x); rotX[2][1] = -rotX[1][2]; vnl_matrix_fixed< double, 3, 3 > rotY; rotY.set_identity(); rotY[0][0] = cos(y); rotY[2][2] = rotY[0][0]; rotY[0][2] = sin(y); rotY[2][0] = -rotY[0][2]; vnl_matrix_fixed< double, 3, 3 > rotZ; rotZ.set_identity(); rotZ[0][0] = cos(z); rotZ[1][1] = rotZ[0][0]; rotZ[0][1] = -sin(z); rotZ[1][0] = -rotZ[0][1]; mitk::Geometry3D::Pointer geom = this->GetGeometry(); mitk::Point3D center = geom->GetCenter(); boost::progress_display disp(m_NumFibers); vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); for (int i=0; iGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(j); vnl_vector_fixed< double, 3 > dir; dir[0] = p[0]-center[0]; dir[1] = p[1]-center[1]; dir[2] = p[2]-center[2]; dir = rotZ*rotY*rotX*dir; dir[0] += center[0]; dir[1] += center[1]; dir[2] += center[2]; vtkIdType id = vtkNewPoints->InsertNextPoint(dir.data_block()); container->GetPointIds()->InsertNextId(id); } vtkNewCells->InsertNextCell(container); } m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); UpdateColorCoding(); UpdateFiberGeometry(); } void mitk::FiberBundleX::ScaleFibers(double x, double y, double z) { MITK_INFO << "Scaling fibers"; boost::progress_display disp(m_NumFibers); mitk::Geometry3D* geom = this->GetGeometry(); mitk::Point3D c = geom->GetCenter(); vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); for (int i=0; iGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(j); p[0] -= c[0]; p[1] -= c[1]; p[2] -= c[2]; p[0] *= x; p[1] *= y; p[2] *= z; p[0] += c[0]; p[1] += c[1]; p[2] += c[2]; vtkIdType id = vtkNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } vtkNewCells->InsertNextCell(container); } m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); UpdateColorCoding(); UpdateFiberGeometry(); } void mitk::FiberBundleX::TranslateFibers(double x, double y, double z) { MITK_INFO << "Translating fibers"; boost::progress_display disp(m_NumFibers); vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); for (int i=0; iGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(j); p[0] += x; p[1] += y; p[2] += z; vtkIdType id = vtkNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } vtkNewCells->InsertNextCell(container); } m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); UpdateColorCoding(); UpdateFiberGeometry(); } void mitk::FiberBundleX::MirrorFibers(unsigned int axis) { if (axis>2) return; MITK_INFO << "Mirroring fibers"; boost::progress_display disp(m_NumFibers); vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); for (int i=0; iGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(j); p[axis] = -p[axis]; vtkIdType id = vtkNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } vtkNewCells->InsertNextCell(container); } m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); UpdateColorCoding(); UpdateFiberGeometry(); } bool mitk::FiberBundleX::ApplyCurvatureThreshold(float minRadius, bool deleteFibers) { if (minRadius<0) return true; vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); MITK_INFO << "Applying curvature threshold"; boost::progress_display disp(m_FiberPolyData->GetNumberOfCells()); for (int i=0; iGetNumberOfCells(); i++) { ++disp ; vtkCell* cell = m_FiberPolyData->GetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); // calculate curvatures vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(j, p1); double p2[3]; points->GetPoint(j+1, p2); double p3[3]; points->GetPoint(j+2, p3); vnl_vector_fixed< float, 3 > v1, v2, v3; v1[0] = p2[0]-p1[0]; v1[1] = p2[1]-p1[1]; v1[2] = p2[2]-p1[2]; v2[0] = p3[0]-p2[0]; v2[1] = p3[1]-p2[1]; v2[2] = p3[2]-p2[2]; v3[0] = p1[0]-p3[0]; v3[1] = p1[1]-p3[1]; v3[2] = p1[2]-p3[2]; float a = v1.magnitude(); float b = v2.magnitude(); float c = v3.magnitude(); float r = a*b*c/std::sqrt((a+b+c)*(a+b-c)*(b+c-a)*(a-b+c)); // radius of triangle via Heron's formula (area of triangle) vtkIdType id = vtkNewPoints->InsertNextPoint(p1); container->GetPointIds()->InsertNextId(id); if (deleteFibers && rInsertNextCell(container); container = vtkSmartPointer::New(); } else if (j==numPoints-3) { id = vtkNewPoints->InsertNextPoint(p2); container->GetPointIds()->InsertNextId(id); id = vtkNewPoints->InsertNextPoint(p3); container->GetPointIds()->InsertNextId(id); vtkNewCells->InsertNextCell(container); } } } if (vtkNewCells->GetNumberOfCells()<=0) return false; m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); UpdateColorCoding(); UpdateFiberGeometry(); return true; } bool mitk::FiberBundleX::RemoveShortFibers(float lengthInMM) { if (lengthInMM<=0 || lengthInMMm_MaxFiberLength) // can't remove all fibers return false; vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); float min = m_MaxFiberLength; MITK_INFO << "Removing short fibers"; boost::progress_display disp(m_NumFibers); for (int i=0; iGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); if (m_FiberLengths.at(i)>=lengthInMM) { vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(j); vtkIdType id = vtkNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } vtkNewCells->InsertNextCell(container); if (m_FiberLengths.at(i)GetNumberOfCells()<=0) return false; m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); UpdateColorCoding(); UpdateFiberGeometry(); return true; } bool mitk::FiberBundleX::RemoveLongFibers(float lengthInMM) { if (lengthInMM<=0 || lengthInMM>m_MaxFiberLength) return true; if (lengthInMM vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); MITK_INFO << "Removing long fibers"; boost::progress_display disp(m_NumFibers); for (int i=0; iGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); if (m_FiberLengths.at(i)<=lengthInMM) { vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(j); vtkIdType id = vtkNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } vtkNewCells->InsertNextCell(container); } } if (vtkNewCells->GetNumberOfCells()<=0) return false; m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); UpdateColorCoding(); UpdateFiberGeometry(); return true; } void mitk::FiberBundleX::DoFiberSmoothing(int pointsPerCm, double tension, double continuity, double bias ) { vtkSmartPointer vtkSmoothPoints = vtkSmartPointer::New(); //in smoothpoints the interpolated points representing a fiber are stored. //in vtkcells all polylines are stored, actually all id's of them are stored vtkSmartPointer vtkSmoothCells = vtkSmartPointer::New(); //cellcontainer for smoothed lines vtkIdType pointHelperCnt = 0; MITK_INFO << "Resampling fibers"; boost::progress_display disp(m_NumFibers); for (int i=0; iGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkSmartPointer newPoints = vtkSmartPointer::New(); for (int j=0; jInsertNextPoint(points->GetPoint(j)); float length = m_FiberLengths.at(i); length /=10; int sampling = pointsPerCm*length; vtkSmartPointer xSpline = vtkSmartPointer::New(); vtkSmartPointer ySpline = vtkSmartPointer::New(); vtkSmartPointer zSpline = vtkSmartPointer::New(); xSpline->SetDefaultBias(bias); xSpline->SetDefaultTension(tension); xSpline->SetDefaultContinuity(continuity); ySpline->SetDefaultBias(bias); ySpline->SetDefaultTension(tension); ySpline->SetDefaultContinuity(continuity); zSpline->SetDefaultBias(bias); zSpline->SetDefaultTension(tension); zSpline->SetDefaultContinuity(continuity); vtkSmartPointer spline = vtkSmartPointer::New(); spline->SetXSpline(xSpline); spline->SetYSpline(ySpline); spline->SetZSpline(zSpline); spline->SetPoints(newPoints); vtkSmartPointer functionSource = vtkSmartPointer::New(); functionSource->SetParametricFunction(spline); functionSource->SetUResolution(sampling); functionSource->SetVResolution(sampling); functionSource->SetWResolution(sampling); functionSource->Update(); vtkPolyData* outputFunction = functionSource->GetOutput(); vtkPoints* tmpSmoothPnts = outputFunction->GetPoints(); //smoothPoints of current fiber vtkSmartPointer smoothLine = vtkSmartPointer::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(); m_FiberSampling = pointsPerCm; } void mitk::FiberBundleX::DoFiberSmoothing(int pointsPerCm) { DoFiberSmoothing(pointsPerCm, 0, 0, 0 ); } // Resample fiber to get equidistant points void mitk::FiberBundleX::ResampleFibers(float pointDistance) { if (pointDistance<=0.00001) return; vtkSmartPointer newPoly = vtkSmartPointer::New(); vtkSmartPointer newCellArray = vtkSmartPointer::New(); vtkSmartPointer newPoints = vtkSmartPointer::New(); int numberOfLines = m_NumFibers; MITK_INFO << "Resampling fibers"; boost::progress_display disp(m_NumFibers); for (int i=0; iGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); double* point = points->GetPoint(0); vtkIdType pointId = newPoints->InsertNextPoint(point); container->GetPointIds()->InsertNextId(pointId); float dtau = 0; int cur_p = 1; itk::Vector dR; float normdR = 0; for (;;) { while (dtau <= pointDistance && cur_p < numPoints) { itk::Vector v1; point = points->GetPoint(cur_p-1); v1[0] = point[0]; v1[1] = point[1]; v1[2] = point[2]; itk::Vector v2; point = points->GetPoint(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 >= pointDistance) { itk::Vector v1; point = points->GetPoint(cur_p-1); v1[0] = point[0]; v1[1] = point[1]; v1[2] = point[2]; itk::Vector v2 = v1 - dR*( (dtau-pointDistance)/normdR ); pointId = newPoints->InsertNextPoint(v2.GetDataPointer()); container->GetPointIds()->InsertNextId(pointId); } else { point = points->GetPoint(numPoints-1); pointId = newPoints->InsertNextPoint(point); container->GetPointIds()->InsertNextId(pointId); break; } dtau = dtau-pointDistance; } newCellArray->InsertNextCell(container); } newPoly->SetPoints(newPoints); newPoly->SetLines(newCellArray); m_FiberPolyData = newPoly; UpdateFiberGeometry(); UpdateColorCoding(); m_FiberSampling = 10/pointDistance; } // reapply selected colorcoding in case polydata structure has changed void mitk::FiberBundleX::UpdateColorCoding() { char* cc = GetCurrentColorCoding(); if( strcmp (COLORCODING_ORIENTATION_BASED,cc) == 0 ) DoColorCodingOrientationBased(); else if( strcmp (COLORCODING_FA_BASED,cc) == 0 ) DoColorCodingFaBased(); } // reapply selected colorcoding in case polydata structure has changed bool mitk::FiberBundleX::Equals(mitk::FiberBundleX* fib) { if (fib==NULL) return false; mitk::FiberBundleX::Pointer tempFib = this->SubtractBundle(fib); mitk::FiberBundleX::Pointer tempFib2 = fib->SubtractBundle(this); if (tempFib.IsNull() && tempFib2.IsNull()) return true; return false; } /* ESSENTIAL IMPLEMENTATION OF SUPERCLASS METHODS */ void mitk::FiberBundleX::UpdateOutputInformation() { } void mitk::FiberBundleX::SetRequestedRegionToLargestPossibleRegion() { } bool mitk::FiberBundleX::RequestedRegionIsOutsideOfTheBufferedRegion() { return false; } bool mitk::FiberBundleX::VerifyRequestedRegion() { return true; } void mitk::FiberBundleX::SetRequestedRegion(const itk::DataObject *data ) { } diff --git a/Modules/DiffusionImaging/FiberTracking/MiniApps/CMakeLists.txt b/Modules/DiffusionImaging/FiberTracking/MiniApps/CMakeLists.txt new file mode 100755 index 0000000000..8e7361bdee --- /dev/null +++ b/Modules/DiffusionImaging/FiberTracking/MiniApps/CMakeLists.txt @@ -0,0 +1,49 @@ +OPTION(BUILD_FiberTrackingMiniApps "Build commandline tools for fiber tracking" OFF) + +IF(BUILD_FiberTrackingMiniApps) + + # include necessary modules here MitkExt QmitkExt + MITK_CHECK_MODULE(_RESULT DiffusionCore FiberTracking ) + IF(_RESULT) + MESSAGE("Warning: FiberTrackingMiniApps is missing ${_RESULT}") + ELSE(_RESULT) + MITK_USE_MODULE( DiffusionCore FiberTracking ) + + # needed include directories + INCLUDE_DIRECTORIES( + ${CMAKE_CURRENT_SOURCE_DIR} + ${CMAKE_CURRENT_BINARY_DIR} + ${ALL_INCLUDE_DIRECTORIES}) + + PROJECT( mitkFiberTrackingMiniApps ) + + # fill in the standalone executables here + SET(FIBERTRACKINGMINIAPPS + mitkFiberTrackingMiniApps + ) + + # set additional files here + SET(FIBERTRACKING_ADDITIONAL_FILES + MiniAppManager.cpp + GibbsTracking.cpp + StreamlineTracking.cpp + TractometerAngularErrorTool.cpp + ) + + # create an executable foreach tool (only one at the moment) + FOREACH(tool ${FIBERTRACKINGMINIAPPS}) + ADD_EXECUTABLE( + ${tool} + ${tool}.cpp + ${FIBERTRACKING_ADDITIONAL_FILES} + ) + + TARGET_LINK_LIBRARIES( + ${tool} + ${ALL_LIBRARIES} ) + ENDFOREACH(tool) + ENDIF() + + MITK_INSTALL_TARGETS(EXECUTABLES mitkFiberTrackingMiniApps ) + +ENDIF(BUILD_FiberTrackingMiniApps) diff --git a/Modules/DiffusionImaging/FiberTracking/MiniApps/GibbsTracking.cpp b/Modules/DiffusionImaging/FiberTracking/MiniApps/GibbsTracking.cpp new file mode 100755 index 0000000000..f15d2bb328 --- /dev/null +++ b/Modules/DiffusionImaging/FiberTracking/MiniApps/GibbsTracking.cpp @@ -0,0 +1,217 @@ +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +All rights reserved. + +This software is distributed WITHOUT ANY WARRANTY; without +even the implied warranty of MERCHANTABILITY or FITNESS FOR +A PARTICULAR PURPOSE. + +See LICENSE.txt or http://www.mitk.org for details. + +===================================================================*/ + +#include "MiniAppManager.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "ctkCommandLineParser.h" + +template +typename itk::ShCoefficientImageImporter< float, shOrder >::QballImageType::Pointer TemplatedConvertShCoeffs(mitk::Image* mitkImg, int toolkit) +{ + typedef itk::ShCoefficientImageImporter< float, shOrder > FilterType; + typedef mitk::ImageToItk< itk::Image< float, 4 > > CasterType; + CasterType::Pointer caster = CasterType::New(); + caster->SetInput(mitkImg); + caster->Update(); + + typename FilterType::Pointer filter = FilterType::New(); + + switch (toolkit) + { + case 0: + filter->SetToolkit(FilterType::FSL); + break; + case 1: + filter->SetToolkit(FilterType::MRTRIX); + break; + default: + filter->SetToolkit(FilterType::FSL); + } + + filter->SetInputImage(caster->GetOutput()); + filter->GenerateData(); + return filter->GetQballImage(); +} + +int GibbsTracking(int argc, char* argv[]) +{ + ctkCommandLineParser parser; + parser.setArgumentPrefix("--", "-"); + parser.addArgument("input", "i", ctkCommandLineParser::String, "input image (tensor, Q-ball or FSL/MRTrix SH-coefficient image)", mitk::Any(), false); + parser.addArgument("parameters", "p", ctkCommandLineParser::String, "parameter file (.gtp)", mitk::Any(), false); + parser.addArgument("mask", "m", ctkCommandLineParser::String, "binary mask image"); + parser.addArgument("shConvention", "s", ctkCommandLineParser::String, "sh coefficient convention (FSL, MRtrix) (default = FSL)"); + parser.addArgument("outFile", "o", ctkCommandLineParser::String, "output fiber bundle (.fib)", mitk::Any(), false); + + map parsedArgs = parser.parseArguments(argc, argv); + if (parsedArgs.size()==0) + return EXIT_FAILURE; + + string inFileName = mitk::any_cast(parsedArgs["input"]); + string paramFileName = mitk::any_cast(parsedArgs["parameters"]); + string outFileName = mitk::any_cast(parsedArgs["outFile"]); + + try + { + RegisterDiffusionCoreObjectFactory(); + RegisterFiberTrackingObjectFactory(); + + // instantiate gibbs tracker + typedef itk::Vector OdfVectorType; + typedef itk::Image ItkQballImageType; + typedef itk::GibbsTrackingFilter GibbsTrackingFilterType; + GibbsTrackingFilterType::Pointer gibbsTracker = GibbsTrackingFilterType::New(); + + // load input image + const std::string s1="", s2=""; + std::vector infile = mitk::BaseDataIO::LoadBaseDataFromFile( inFileName, s1, s2, false ); + + // try to cast to qball image + QString inImageName(inFileName.c_str()); + if( inImageName.endsWith(".qbi") ) + { + MITK_INFO << "Loading qball image ..."; + mitk::QBallImage::Pointer mitkQballImage = dynamic_cast(infile.at(0).GetPointer()); + ItkQballImageType::Pointer itk_qbi = ItkQballImageType::New(); + mitk::CastToItkImage(mitkQballImage, itk_qbi); + gibbsTracker->SetQBallImage(itk_qbi.GetPointer()); + } + else if( inImageName.endsWith(".dti") ) + { + MITK_INFO << "Loading tensor image ..."; + typedef itk::Image< itk::DiffusionTensor3D, 3 > ItkTensorImage; + mitk::TensorImage::Pointer mitkTensorImage = dynamic_cast(infile.at(0).GetPointer()); + ItkTensorImage::Pointer itk_dti = ItkTensorImage::New(); + mitk::CastToItkImage(mitkTensorImage, itk_dti); + gibbsTracker->SetTensorImage(itk_dti); + } + else if ( inImageName.endsWith(".nii") ) + { + MITK_INFO << "Loading sh-coefficient image ..."; + mitk::Image::Pointer mitkImage = dynamic_cast(infile.at(0).GetPointer()); + + int nrCoeffs = mitkImage->GetLargestPossibleRegion().GetSize()[3]; + int c=3, d=2-2*nrCoeffs; + double D = c*c-4*d; + int shOrder; + if (D>0) + { + shOrder = (-c+sqrt(D))/2.0; + if (shOrder<0) + shOrder = (-c-sqrt(D))/2.0; + } + else if (D==0) + shOrder = -c/2.0; + + MITK_INFO << "using SH-order " << shOrder; + + int toolkitConvention = 0; + + if (parsedArgs.count("shConvention")) + { + QString convention(mitk::any_cast(parsedArgs["shConvention"]).c_str()); + if ( convention.compare("MRtrix", Qt::CaseInsensitive)==0 ) + { + toolkitConvention = 1; + MITK_INFO << "Using MRtrix style sh-coefficient convention"; + } + else + MITK_INFO << "Using FSL style sh-coefficient convention"; + } + else + MITK_INFO << "Using FSL style sh-coefficient convention"; + + switch (shOrder) + { + case 4: + gibbsTracker->SetQBallImage(TemplatedConvertShCoeffs<4>(mitkImage, toolkitConvention)); + break; + case 6: + gibbsTracker->SetQBallImage(TemplatedConvertShCoeffs<6>(mitkImage, toolkitConvention)); + break; + case 8: + gibbsTracker->SetQBallImage(TemplatedConvertShCoeffs<8>(mitkImage, toolkitConvention)); + break; + case 10: + gibbsTracker->SetQBallImage(TemplatedConvertShCoeffs<10>(mitkImage, toolkitConvention)); + break; + case 12: + gibbsTracker->SetQBallImage(TemplatedConvertShCoeffs<12>(mitkImage, toolkitConvention)); + break; + default: + MITK_INFO << "SH-order " << shOrder << " not supported"; + } + } + else + return EXIT_FAILURE; + + // global tracking + if (parsedArgs.count("mask")) + { + typedef itk::Image MaskImgType; + mitk::Image::Pointer mitkMaskImage = mitk::IOUtil::LoadImage(mitk::any_cast(parsedArgs["mask"])); + MaskImgType::Pointer itk_mask = MaskImgType::New(); + mitk::CastToItkImage(mitkMaskImage, itk_mask); + gibbsTracker->SetMaskImage(itk_mask); + } + + gibbsTracker->SetDuplicateImage(false); + gibbsTracker->SetLoadParameterFile( paramFileName ); +// gibbsTracker->SetLutPath( "" ); + gibbsTracker->Update(); + + mitk::FiberBundleX::Pointer mitkFiberBundle = mitk::FiberBundleX::New(gibbsTracker->GetFiberBundle()); + + mitk::CoreObjectFactory::FileWriterList fileWriters = mitk::CoreObjectFactory::GetInstance()->GetFileWriters(); + for (mitk::CoreObjectFactory::FileWriterList::iterator it = fileWriters.begin() ; it != fileWriters.end() ; ++it) + { + if ( (*it)->CanWriteBaseDataType(mitkFiberBundle.GetPointer()) ) { + (*it)->SetFileName( outFileName.c_str() ); + (*it)->DoWrite( mitkFiberBundle.GetPointer() ); + } + } + } + catch (itk::ExceptionObject e) + { + MITK_INFO << e; + return EXIT_FAILURE; + } + catch (std::exception e) + { + MITK_INFO << e.what(); + return EXIT_FAILURE; + } + catch (...) + { + MITK_INFO << "ERROR!?!"; + return EXIT_FAILURE; + } + MITK_INFO << "DONE"; + return EXIT_SUCCESS; +} +RegisterFiberTrackingMiniApp(GibbsTracking); diff --git a/Modules/DiffusionImaging/FiberTracking/MiniApps/MiniAppManager.cpp b/Modules/DiffusionImaging/FiberTracking/MiniApps/MiniAppManager.cpp new file mode 100755 index 0000000000..ff9b82b1ad --- /dev/null +++ b/Modules/DiffusionImaging/FiberTracking/MiniApps/MiniAppManager.cpp @@ -0,0 +1,91 @@ +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +All rights reserved. + +This software is distributed WITHOUT ANY WARRANTY; without +even the implied warranty of MERCHANTABILITY or FITNESS FOR +A PARTICULAR PURPOSE. + +See LICENSE.txt or http://www.mitk.org for details. + +===================================================================*/ + +#include "MiniAppManager.h" + +MiniAppManager* MiniAppManager::GetInstance() +{ + static MiniAppManager instance; + return &instance; +} + +// Attention: Name of the miniApp must be the last argument!!! +// it will be cut off from the rest of the arguments and then +// the app will be run +int MiniAppManager::RunMiniApp(int argc, char* argv[]) +{ + try + { + std::string nameOfMiniApp; + std::map< std::string, MiniAppFunction >::iterator it = m_Functions.begin(); + + if( argc < 2) + { + std::cout << "Please choose the mini app to execute: " << std::endl; + + for(int i=0; it != m_Functions.end(); ++i,++it) + { + std::cout << "(" << i << ")" << " " << it->first << std::endl; + } + std::cout << "Please select: "; + int choose; + std::cin >> choose; + + it = m_Functions.begin(); + std::advance(it, choose); + if( it != m_Functions.end() ) + nameOfMiniApp = it->first; + } + else + { + nameOfMiniApp = argv[1]; + //--argc; + } + + it = m_Functions.find(nameOfMiniApp); + if(it == m_Functions.end()) + { + std::ostringstream s; s << "MiniApp (" << nameOfMiniApp << ") not found!"; + throw std::invalid_argument(s.str().c_str()); + } + + MITK_INFO << "Start " << nameOfMiniApp << " .."; + MiniAppFunction func = it->second; + return func( argc, argv ); + } + + catch(std::exception& e) + { + MITK_ERROR << e.what(); + } + + catch(...) + { + MITK_ERROR << "Unknown error occurred"; + } + + return EXIT_FAILURE; +} + +///////////////////// +// MiniAppFunction // +///////////////////// +MiniAppManager::MiniAppFunction +MiniAppManager::AddFunction(const std::string& name, MiniAppFunction func) +{ + m_Functions.insert( std::pair(name, func) ); + return func; +} diff --git a/Modules/DiffusionImaging/FiberTracking/MiniApps/MiniAppManager.h b/Modules/DiffusionImaging/FiberTracking/MiniApps/MiniAppManager.h new file mode 100755 index 0000000000..989baad3cc --- /dev/null +++ b/Modules/DiffusionImaging/FiberTracking/MiniApps/MiniAppManager.h @@ -0,0 +1,38 @@ +#ifndef MiniAppManager_h +#define MiniAppManager_h + +#include + +struct MiniAppManager +{ + +public: + + typedef int (*MiniAppFunction)(int argc, char* argv[]); + +public: + + static MiniAppManager* GetInstance(); + + // Attention: Name of the miniApp must be the last argument!!! + // it will be cut off from the rest of the arguments and then + // the app will be run + int RunMiniApp(int argc, char* argv[]); + + // + // Add miniApp + // + MiniAppFunction AddFunction(const std::string& name, MiniAppFunction func); + +protected: + + std::map< std::string, MiniAppFunction > m_Functions; +}; + +// +// Register miniApps +// +#define RegisterFiberTrackingMiniApp(functionName) \ + static MiniAppManager::MiniAppFunction MiniApp##functionName = \ + MiniAppManager::GetInstance()->AddFunction(#functionName, &functionName) +#endif diff --git a/Modules/DiffusionImaging/FiberTracking/MiniApps/StreamlineTracking.cpp b/Modules/DiffusionImaging/FiberTracking/MiniApps/StreamlineTracking.cpp new file mode 100755 index 0000000000..f95d8dfe8d --- /dev/null +++ b/Modules/DiffusionImaging/FiberTracking/MiniApps/StreamlineTracking.cpp @@ -0,0 +1,174 @@ +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +All rights reserved. + +This software is distributed WITHOUT ANY WARRANTY; without +even the implied warranty of MERCHANTABILITY or FITNESS FOR +A PARTICULAR PURPOSE. + +See LICENSE.txt or http://www.mitk.org for details. + +===================================================================*/ + +#include "MiniAppManager.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "ctkCommandLineParser.h" + +int StreamlineTracking(int argc, char* argv[]) +{ + ctkCommandLineParser parser; + parser.setArgumentPrefix("--", "-"); + parser.addArgument("input", "i", ctkCommandLineParser::String, "input tensor image (.dti)", mitk::Any(), false); + parser.addArgument("seed", "s", ctkCommandLineParser::String, "binary seed image"); + parser.addArgument("mask", "m", ctkCommandLineParser::String, "binary mask image"); + parser.addArgument("minFA", "t", ctkCommandLineParser::Float, "minimum fractional anisotropy threshold", 0.15, true); + parser.addArgument("minCurv", "c", ctkCommandLineParser::Float, "minimum curvature radius in mm (default = 0.5*minimum-spacing)"); + parser.addArgument("stepSize", "s", ctkCommandLineParser::Float, "stepsize in mm (default = 0.1*minimum-spacing)"); + parser.addArgument("tendf", "f", ctkCommandLineParser::Float, "Weighting factor between first eigenvector (f=1 equals FACT tracking) and input vector dependent direction (f=0).", 1.0, true); + parser.addArgument("tendg", "g", ctkCommandLineParser::Float, "Weighting factor between input vector (g=0) and tensor deflection (g=1 equals TEND tracking)", 0.0, true); + parser.addArgument("numSeeds", "n", ctkCommandLineParser::Int, "Number of seeds per voxel.", 1, true); + parser.addArgument("minLength", "l", ctkCommandLineParser::Float, "minimum fiber length in mm", 20, true); + + parser.addArgument("interpolate", "a", ctkCommandLineParser::Bool, "Use linear interpolation", false, true); + parser.addArgument("outFile", "o", ctkCommandLineParser::String, "output fiber bundle (.fib)", mitk::Any(), false); + + map parsedArgs = parser.parseArguments(argc, argv); + if (parsedArgs.size()==0) + return EXIT_FAILURE; + + string dtiFileName = mitk::any_cast(parsedArgs["input"]); + string outFileName = mitk::any_cast(parsedArgs["outFile"]); + + float minFA = 0.15; + float minCurv = -1; + float stepSize = -1; + float tendf = 1; + float tendg = 0; + float minLength = 20; + int numSeeds = 1; + bool interpolate = false; + + if (parsedArgs.count("minCurv")) + minCurv = mitk::any_cast(parsedArgs["minCurv"]); + if (parsedArgs.count("minFA")) + minFA = mitk::any_cast(parsedArgs["minFA"]); + if (parsedArgs.count("stepSize")) + stepSize = mitk::any_cast(parsedArgs["stepSize"]); + if (parsedArgs.count("tendf")) + tendf = mitk::any_cast(parsedArgs["tendf"]); + if (parsedArgs.count("tendg")) + tendg = mitk::any_cast(parsedArgs["tendg"]); + if (parsedArgs.count("minLength")) + minLength = mitk::any_cast(parsedArgs["minLength"]); + if (parsedArgs.count("numSeeds")) + numSeeds = mitk::any_cast(parsedArgs["numSeeds"]); + + + if (parsedArgs.count("interpolate")) + interpolate = mitk::any_cast(parsedArgs["interpolate"]); + + + + try + { + RegisterDiffusionCoreObjectFactory(); + RegisterFiberTrackingObjectFactory(); + + // load input image + const std::string s1="", s2=""; + std::vector infile = mitk::BaseDataIO::LoadBaseDataFromFile( dtiFileName, s1, s2, false ); + + MITK_INFO << "Loading tensor image ..."; + typedef itk::Image< itk::DiffusionTensor3D, 3 > ItkTensorImage; + mitk::TensorImage::Pointer mitkTensorImage = dynamic_cast(infile.at(0).GetPointer()); + ItkTensorImage::Pointer itk_dti = ItkTensorImage::New(); + mitk::CastToItkImage(mitkTensorImage, itk_dti); + + MITK_INFO << "Loading seed image ..."; + typedef itk::Image< unsigned char, 3 > ItkUCharImageType; + mitk::Image::Pointer mitkSeedImage = NULL; + if (parsedArgs.count("seed")) + mitkSeedImage = mitk::IOUtil::LoadImage(mitk::any_cast(parsedArgs["seed"])); + + MITK_INFO << "Loading mask image ..."; + mitk::Image::Pointer mitkMaskImage = NULL; + if (parsedArgs.count("mask")) + mitkMaskImage = mitk::IOUtil::LoadImage(mitk::any_cast(parsedArgs["mask"])); + + // instantiate tracker + typedef itk::StreamlineTrackingFilter< float > FilterType; + FilterType::Pointer filter = FilterType::New(); + filter->SetInput(itk_dti); + filter->SetSeedsPerVoxel(numSeeds); + filter->SetFaThreshold(minFA); + filter->SetMinCurvatureRadius(minCurv); + filter->SetStepSize(stepSize); + filter->SetF(tendf); + filter->SetG(tendg); + filter->SetInterpolate(interpolate); + filter->SetMinTractLength(minLength); + + if (mitkSeedImage.IsNotNull()) + { + ItkUCharImageType::Pointer mask = ItkUCharImageType::New(); + mitk::CastToItkImage(mitkSeedImage, mask); + filter->SetSeedImage(mask); + } + + if (mitkMaskImage.IsNotNull()) + { + ItkUCharImageType::Pointer mask = ItkUCharImageType::New(); + mitk::CastToItkImage(mitkMaskImage, mask); + filter->SetMaskImage(mask); + } + + filter->Update(); + + vtkSmartPointer fiberBundle = filter->GetFiberPolyData(); + if ( fiberBundle->GetNumberOfLines()==0 ) + { + MITK_INFO << "No fibers reconstructed. Check parametrization."; + return EXIT_FAILURE; + } + mitk::FiberBundleX::Pointer fib = mitk::FiberBundleX::New(fiberBundle); + + mitk::CoreObjectFactory::FileWriterList fileWriters = mitk::CoreObjectFactory::GetInstance()->GetFileWriters(); + for (mitk::CoreObjectFactory::FileWriterList::iterator it = fileWriters.begin() ; it != fileWriters.end() ; ++it) + { + if ( (*it)->CanWriteBaseDataType(fib.GetPointer()) ) { + (*it)->SetFileName( outFileName.c_str() ); + (*it)->DoWrite( fib.GetPointer() ); + } + } + } + catch (itk::ExceptionObject e) + { + MITK_INFO << e; + return EXIT_FAILURE; + } + catch (std::exception e) + { + MITK_INFO << e.what(); + return EXIT_FAILURE; + } + catch (...) + { + MITK_INFO << "ERROR!?!"; + return EXIT_FAILURE; + } + MITK_INFO << "DONE"; + return EXIT_SUCCESS; +} +RegisterFiberTrackingMiniApp(StreamlineTracking); diff --git a/Modules/DiffusionImaging/FiberTracking/MiniApps/TractometerAngularErrorTool.cpp b/Modules/DiffusionImaging/FiberTracking/MiniApps/TractometerAngularErrorTool.cpp new file mode 100755 index 0000000000..00996996dd --- /dev/null +++ b/Modules/DiffusionImaging/FiberTracking/MiniApps/TractometerAngularErrorTool.cpp @@ -0,0 +1,247 @@ +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +All rights reserved. + +This software is distributed WITHOUT ANY WARRANTY; without +even the implied warranty of MERCHANTABILITY or FITNESS FOR +A PARTICULAR PURPOSE. + +See LICENSE.txt or http://www.mitk.org for details. + +===================================================================*/ + +#include "MiniAppManager.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include "ctkCommandLineParser.h" +#include "ctkCommandLineParser.cpp" +#include +#include +#include +#include +#include +#include +#include + +#define _USE_MATH_DEFINES +#include + +int TractometerAngularErrorTool(int argc, char* argv[]) +{ + ctkCommandLineParser parser; + parser.setArgumentPrefix("--", "-"); + parser.addArgument("input", "i", ctkCommandLineParser::String, "input tractogram (.fib, vtk ascii file format)", mitk::Any(), false); + parser.addArgument("reference", "r", ctkCommandLineParser::StringList, "reference direction images", mitk::Any(), false); + parser.addArgument("out", "o", ctkCommandLineParser::String, "output root", mitk::Any(), false); + parser.addArgument("mask", "m", ctkCommandLineParser::String, "mask image"); + parser.addArgument("athresh", "a", ctkCommandLineParser::Float, "angular threshold in degrees. closer fiber directions are regarded as one direction and clustered together.", 25, true); + parser.addArgument("verbose", "v", ctkCommandLineParser::Bool, "output optional and intermediate calculation results"); + + map parsedArgs = parser.parseArguments(argc, argv); + if (parsedArgs.size()==0) + return EXIT_FAILURE; + + ctkCommandLineParser::StringContainerType referenceImages = mitk::any_cast(parsedArgs["reference"]); + + string fibFile = mitk::any_cast(parsedArgs["input"]); + + string maskImage(""); + if (parsedArgs.count("mask")) + maskImage = mitk::any_cast(parsedArgs["mask"]); + + float angularThreshold = 25; + if (parsedArgs.count("athresh")) + angularThreshold = mitk::any_cast(parsedArgs["athresh"]); + + string outRoot = mitk::any_cast(parsedArgs["out"]); + + bool verbose = false; + if (parsedArgs.count("verbose")) + verbose = mitk::any_cast(parsedArgs["verbose"]); + + try + { + RegisterDiffusionCoreObjectFactory(); + RegisterFiberTrackingObjectFactory(); + + typedef itk::Image ItkUcharImgType; + typedef itk::Image< itk::Vector< float, 3>, 3 > ItkDirectionImage3DType; + typedef itk::VectorContainer< int, ItkDirectionImage3DType::Pointer > ItkDirectionImageContainerType; + typedef itk::EvaluateDirectionImagesFilter< float > EvaluationFilterType; + + // load fiber bundle + mitk::FiberBundleX::Pointer inputTractogram = dynamic_cast(mitk::IOUtil::LoadDataNode(fibFile)->GetData()); + + // load reference directions + ItkDirectionImageContainerType::Pointer referenceImageContainer = ItkDirectionImageContainerType::New(); + for (int i=0; i(mitk::IOUtil::LoadDataNode(referenceImages.at(i))->GetData()); + typedef mitk::ImageToItk< ItkDirectionImage3DType > CasterType; + CasterType::Pointer caster = CasterType::New(); + caster->SetInput(img); + caster->Update(); + ItkDirectionImage3DType::Pointer itkImg = caster->GetOutput(); + referenceImageContainer->InsertElement(referenceImageContainer->Size(),itkImg); + } + catch(...){ MITK_INFO << "could not load: " << referenceImages.at(i); } + } + + // load/create mask image + ItkUcharImgType::Pointer itkMaskImage = ItkUcharImgType::New(); + if (maskImage.compare("")==0) + { + ItkDirectionImage3DType::Pointer dirImg = referenceImageContainer->GetElement(0); + itkMaskImage->SetSpacing( dirImg->GetSpacing() ); + itkMaskImage->SetOrigin( dirImg->GetOrigin() ); + itkMaskImage->SetDirection( dirImg->GetDirection() ); + itkMaskImage->SetLargestPossibleRegion( dirImg->GetLargestPossibleRegion() ); + itkMaskImage->SetBufferedRegion( dirImg->GetLargestPossibleRegion() ); + itkMaskImage->SetRequestedRegion( dirImg->GetLargestPossibleRegion() ); + itkMaskImage->Allocate(); + itkMaskImage->FillBuffer(1); + } + else + { + mitk::Image::Pointer mitkMaskImage = dynamic_cast(mitk::IOUtil::LoadDataNode(maskImage)->GetData()); + CastToItkImage(mitkMaskImage, itkMaskImage); + } + + + // extract directions from fiber bundle + itk::TractsToVectorImageFilter::Pointer fOdfFilter = itk::TractsToVectorImageFilter::New(); + fOdfFilter->SetFiberBundle(inputTractogram); + fOdfFilter->SetMaskImage(itkMaskImage); + fOdfFilter->SetAngularThreshold(cos(angularThreshold*M_PI/180)); + fOdfFilter->SetNormalizeVectors(true); + fOdfFilter->SetUseWorkingCopy(false); + fOdfFilter->Update(); + ItkDirectionImageContainerType::Pointer directionImageContainer = fOdfFilter->GetDirectionImageContainer(); + + if (verbose) + { + // write vector field + mitk::FiberBundleX::Pointer directions = fOdfFilter->GetOutputFiberBundle(); + mitk::CoreObjectFactory::FileWriterList fileWriters = mitk::CoreObjectFactory::GetInstance()->GetFileWriters(); + for (mitk::CoreObjectFactory::FileWriterList::iterator it = fileWriters.begin() ; it != fileWriters.end() ; ++it) + { + if ( (*it)->CanWriteBaseDataType(directions.GetPointer()) ) { + QString outfilename(outRoot.c_str()); + outfilename += "_VECTOR_FIELD.fib"; + (*it)->SetFileName( outfilename.toStdString().c_str() ); + (*it)->DoWrite( directions.GetPointer() ); + } + } + + // write direction images + for (int i=0; iSize(); i++) + { + itk::TractsToVectorImageFilter::ItkDirectionImageType::Pointer itkImg = directionImageContainer->GetElement(i); + typedef itk::ImageFileWriter< itk::TractsToVectorImageFilter::ItkDirectionImageType > WriterType; + WriterType::Pointer writer = WriterType::New(); + QString outfilename(outRoot.c_str()); + outfilename += "_DIRECTION_"; + outfilename += QString::number(i); + outfilename += ".nrrd"; + MITK_INFO << "writing " << outfilename.toStdString(); + writer->SetFileName(outfilename.toStdString().c_str()); + writer->SetInput(itkImg); + writer->Update(); + } + + // write num direction image + { + ItkUcharImgType::Pointer numDirImage = fOdfFilter->GetNumDirectionsImage(); + typedef itk::ImageFileWriter< ItkUcharImgType > WriterType; + WriterType::Pointer writer = WriterType::New(); + QString outfilename(outRoot.c_str()); + outfilename += "_NUM_DIRECTIONS.nrrd"; + MITK_INFO << "writing " << outfilename.toStdString(); + writer->SetFileName(outfilename.toStdString().c_str()); + writer->SetInput(numDirImage); + writer->Update(); + } + } + + // evaluate directions + EvaluationFilterType::Pointer evaluationFilter = EvaluationFilterType::New(); + evaluationFilter->SetImageSet(directionImageContainer); + evaluationFilter->SetReferenceImageSet(referenceImageContainer); + evaluationFilter->SetMaskImage(itkMaskImage); + evaluationFilter->Update(); + + if (verbose) + { + EvaluationFilterType::OutputImageType::Pointer angularErrorImage = evaluationFilter->GetOutput(0); + typedef itk::ImageFileWriter< EvaluationFilterType::OutputImageType > WriterType; + WriterType::Pointer writer = WriterType::New(); + QString outfilename(outRoot.c_str()); + outfilename += "_ERROR_IMAGE.nrrd"; + writer->SetFileName(outfilename.toStdString().c_str()); + writer->SetInput(angularErrorImage); + writer->Update(); + } + + QString logFile(outRoot.c_str()); logFile += "_ANGULAR_ERROR.csv"; + QFile file(logFile); + file.open(QIODevice::WriteOnly | QIODevice::Text); + QTextStream out(&file); + QString sens = QString("Mean:"); + sens += ","; + sens += QString::number(evaluationFilter->GetMeanAngularError()); + sens += ";"; + + sens += QString("Median:"); + sens += ","; + sens += QString::number(evaluationFilter->GetMedianAngularError()); + sens += ";"; + + sens += QString("Maximum:"); + sens += ","; + sens += QString::number(evaluationFilter->GetMaxAngularError()); + sens += ";"; + + sens += QString("Minimum:"); + sens += ","; + sens += QString::number(evaluationFilter->GetMinAngularError()); + sens += ";"; + + sens += QString("STDEV:"); + sens += ","; + sens += QString::number(std::sqrt(evaluationFilter->GetVarAngularError())); + sens += ";"; + + out << sens; + + MITK_INFO << "DONE"; + } + catch (itk::ExceptionObject e) + { + MITK_INFO << e; + return EXIT_FAILURE; + } + catch (std::exception e) + { + MITK_INFO << e.what(); + return EXIT_FAILURE; + } + catch (...) + { + MITK_INFO << "ERROR!?!"; + return EXIT_FAILURE; + } + return EXIT_SUCCESS; +} +RegisterFiberTrackingMiniApp(TractometerAngularErrorTool); diff --git a/Modules/DiffusionImaging/FiberTracking/MiniApps/ctkCommandLineParser.cpp b/Modules/DiffusionImaging/FiberTracking/MiniApps/ctkCommandLineParser.cpp new file mode 100755 index 0000000000..72da974c5b --- /dev/null +++ b/Modules/DiffusionImaging/FiberTracking/MiniApps/ctkCommandLineParser.cpp @@ -0,0 +1,727 @@ +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +All rights reserved. + +This software is distributed WITHOUT ANY WARRANTY; without +even the implied warranty of MERCHANTABILITY or FITNESS FOR +A PARTICULAR PURPOSE. + +See LICENSE.txt or http://www.mitk.org for details. + +===================================================================*/ +/*========================================================================= + + Library: CTK + + Copyright (c) Kitware Inc. + + Licensed under the Apache License, Version 2.0 (the "License"); + you may not use this file except in compliance with the License. + You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0.txt + + Unless required by applicable law or agreed to in writing, software + distributed under the License is distributed on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + See the License for the specific language governing permissions and + limitations under the License. + +=========================================================================*/ + +// STL includes +#include + + +// CTK includes +#include "ctkCommandLineParser.h" + +using namespace std; + +namespace +{ +// -------------------------------------------------------------------------- +class CommandLineParserArgumentDescription +{ +public: + + + CommandLineParserArgumentDescription( + const string& longArg, const string& longArgPrefix, + const string& shortArg, const string& shortArgPrefix, + ctkCommandLineParser::Type type, const string& argHelp, + const mitk::Any& defaultValue, bool ignoreRest, + bool deprecated, bool optional) + : LongArg(longArg), LongArgPrefix(longArgPrefix), + ShortArg(shortArg), ShortArgPrefix(shortArgPrefix), + ArgHelp(argHelp), IgnoreRest(ignoreRest), NumberOfParametersToProcess(0), + Deprecated(deprecated), DefaultValue(defaultValue), Value(type), ValueType(type), Optional(optional) + { + Value = defaultValue; + + switch (type) + { + case ctkCommandLineParser::String: + { + NumberOfParametersToProcess = 1; + } + break; + case ctkCommandLineParser::Bool: + { + NumberOfParametersToProcess = 0; + } + break; + case ctkCommandLineParser::StringList: + { + NumberOfParametersToProcess = -1; + } + break; + case ctkCommandLineParser::Int: + { + NumberOfParametersToProcess = 1; + } + break; + case ctkCommandLineParser::Float: + { + NumberOfParametersToProcess = 1; + } + break; + default: + MITK_INFO << "Type not supported: " << static_cast(type); + } + + } + + ~CommandLineParserArgumentDescription(){} + + bool addParameter(const string& value); + + string helpText(); + + string LongArg; + string LongArgPrefix; + string ShortArg; + string ShortArgPrefix; + string ArgHelp; + bool IgnoreRest; + int NumberOfParametersToProcess; + bool Deprecated; + bool Optional; + + mitk::Any DefaultValue; + mitk::Any Value; + ctkCommandLineParser::Type ValueType; +}; + +// -------------------------------------------------------------------------- +bool CommandLineParserArgumentDescription::addParameter(const string &value) +{ + switch (ValueType) + { + case ctkCommandLineParser::String: + { + Value = value; + } + break; + case ctkCommandLineParser::Bool: + { + if (value.compare("true")==0) + Value = true; + else + Value = false; + } + break; + case ctkCommandLineParser::StringList: + { + try + { + ctkCommandLineParser::StringContainerType list = mitk::any_cast(Value); + list.push_back(value); + Value = list; + } + catch(...) + { + ctkCommandLineParser::StringContainerType list; + list.push_back(value); + Value = list; + } + } + break; + case ctkCommandLineParser::Int: + { + stringstream ss(value); + int i; + ss >> i; + Value = i; + } + break; + case ctkCommandLineParser::Float: + { + stringstream ss(value); + float f; + ss >> f; + Value = f; + } + break; + default: + return false; + } + + return true; +} + +// -------------------------------------------------------------------------- +string CommandLineParserArgumentDescription::helpText() +{ + string text; + + string shortAndLongArg; + if (!this->ShortArg.empty()) + { + shortAndLongArg = " "; + shortAndLongArg += this->ShortArgPrefix; + shortAndLongArg += this->ShortArg; + } + + if (!this->LongArg.empty()) + { + if (this->ShortArg.empty()) + shortAndLongArg.append(" "); + else + shortAndLongArg.append(", "); + + shortAndLongArg += this->LongArgPrefix; + shortAndLongArg += this->LongArg; + } + + text = text + shortAndLongArg + ", " + this->ArgHelp; + + if (this->Optional) + text += " (optional)"; + + if (!this->DefaultValue.Empty()) + { + text = text + ", (default: " + this->DefaultValue.ToString() + ")"; + } + text += "\n"; + return text; +} + +} + +// -------------------------------------------------------------------------- +// ctkCommandLineParser::ctkInternal class + +// -------------------------------------------------------------------------- +class ctkCommandLineParser::ctkInternal +{ +public: + ctkInternal() + : Debug(false), FieldWidth(0), StrictMode(false) + {} + + ~ctkInternal() { } + + CommandLineParserArgumentDescription* argumentDescription(const string& argument); + + vector ArgumentDescriptionList; + map ArgNameToArgumentDescriptionMap; + map > GroupToArgumentDescriptionListMap; + + StringContainerType UnparsedArguments; + StringContainerType ProcessedArguments; + string ErrorString; + bool Debug; + int FieldWidth; + string LongPrefix; + string ShortPrefix; + string CurrentGroup; + string DisableQSettingsLongArg; + string DisableQSettingsShortArg; + bool StrictMode; +}; + +// -------------------------------------------------------------------------- +// ctkCommandLineParser::ctkInternal methods + +// -------------------------------------------------------------------------- +CommandLineParserArgumentDescription* +ctkCommandLineParser::ctkInternal::argumentDescription(const string& argument) +{ + string unprefixedArg = argument; + + if (!LongPrefix.empty() && argument.compare(0, LongPrefix.size(), LongPrefix)==0) + { + // Case when (ShortPrefix + UnPrefixedArgument) matches LongPrefix + if (argument == LongPrefix && !ShortPrefix.empty() && argument.compare(0, ShortPrefix.size(), ShortPrefix)==0) + { + unprefixedArg = argument.substr(ShortPrefix.size(),argument.size()); + } + else + { + unprefixedArg = argument.substr(LongPrefix.size(),argument.size()); + } + } + else if (!ShortPrefix.empty() && argument.compare(0, ShortPrefix.size(), ShortPrefix)==0) + { + unprefixedArg = argument.substr(ShortPrefix.size(),argument.size()); + } + else if (!LongPrefix.empty() && !ShortPrefix.empty()) + { + return 0; + } + + if (ArgNameToArgumentDescriptionMap.count(unprefixedArg)) + { + return this->ArgNameToArgumentDescriptionMap[unprefixedArg]; + } + return 0; +} + +// -------------------------------------------------------------------------- +// ctkCommandLineParser methods + +// -------------------------------------------------------------------------- +ctkCommandLineParser::ctkCommandLineParser() +{ + this->Internal = new ctkInternal(); +} + +// -------------------------------------------------------------------------- +ctkCommandLineParser::~ctkCommandLineParser() +{ + delete this->Internal; +} + +// -------------------------------------------------------------------------- +map ctkCommandLineParser::parseArguments(const StringContainerType& arguments, + bool* ok) +{ + // Reset + this->Internal->UnparsedArguments.clear(); + this->Internal->ProcessedArguments.clear(); + this->Internal->ErrorString.clear(); + // foreach (CommandLineParserArgumentDescription* desc, this->Internal->ArgumentDescriptionList) + for (int i=0; iArgumentDescriptionList.size(); i++) + { + CommandLineParserArgumentDescription* desc = Internal->ArgumentDescriptionList.at(i); + desc->Value = mitk::Any(desc->ValueType); + if (!desc->DefaultValue.Empty()) + { + desc->Value = desc->DefaultValue; + } + } + + bool error = false; + bool ignoreRest = false; + CommandLineParserArgumentDescription * currentArgDesc = 0; + vector parsedArgDescriptions; + for(int i = 1; i < arguments.size(); ++i) + { + string argument = arguments.at(i); + if (this->Internal->Debug) { MITK_DEBUG << "Processing" << argument; } + + // should argument be ignored ? + if (ignoreRest) + { + if (this->Internal->Debug) + { + MITK_DEBUG << " Skipping: IgnoreRest flag was been set"; + } + this->Internal->UnparsedArguments.push_back(argument); + continue; + } + + // Skip if the argument does not start with the defined prefix + if (!(argument.compare(0, Internal->LongPrefix.size(), Internal->LongPrefix)==0 + || argument.compare(0, Internal->ShortPrefix.size(), Internal->ShortPrefix)==0)) + { + if (this->Internal->StrictMode) + { + this->Internal->ErrorString = "Unknown argument "; + this->Internal->ErrorString += argument; + error = true; + break; + } + if (this->Internal->Debug) + { + MITK_DEBUG << " Skipping: It does not start with the defined prefix"; + } + this->Internal->UnparsedArguments.push_back(argument); + continue; + } + + // Skip if argument has already been parsed ... + bool alreadyProcessed = false; + for (int i=0; iProcessedArguments.size(); i++) + if (argument.compare(Internal->ProcessedArguments.at(i))==0) + { + alreadyProcessed = true; + break; + } + + if (alreadyProcessed) + { + if (this->Internal->StrictMode) + { + this->Internal->ErrorString = "Argument "; + this->Internal->ErrorString += argument; + this->Internal->ErrorString += " already processed !"; + error = true; + break; + } + if (this->Internal->Debug) + { + MITK_DEBUG << " Skipping: Already processed !"; + } + continue; + } + + // Retrieve corresponding argument description + currentArgDesc = this->Internal->argumentDescription(argument); + + // Is there a corresponding argument description ? + if (currentArgDesc) + { + // If the argument is deprecated, print the help text but continue processing + if (currentArgDesc->Deprecated) + { + MITK_WARN << "Deprecated argument " << argument << ": " << currentArgDesc->ArgHelp; + } + else + { + parsedArgDescriptions.push_back(currentArgDesc); + } + + this->Internal->ProcessedArguments.push_back(currentArgDesc->ShortArg); + this->Internal->ProcessedArguments.push_back(currentArgDesc->LongArg); + int numberOfParametersToProcess = currentArgDesc->NumberOfParametersToProcess; + ignoreRest = currentArgDesc->IgnoreRest; + if (this->Internal->Debug && ignoreRest) + { + MITK_DEBUG << " IgnoreRest flag is True"; + } + + // Is the number of parameters associated with the argument being processed known ? + if (numberOfParametersToProcess == 0) + { + currentArgDesc->addParameter("true"); + } + else if (numberOfParametersToProcess > 0) + { + string missingParameterError = + "Argument %1 has %2 value(s) associated whereas exacly %3 are expected."; + for(int j=1; j <= numberOfParametersToProcess; ++j) + { + if (i + j >= arguments.size()) + { +// this->Internal->ErrorString = +// missingParameterError.arg(argument).arg(j-1).arg(numberOfParametersToProcess); +// if (this->Internal->Debug) { MITK_DEBUG << this->Internal->ErrorString; } + if (ok) { *ok = false; } + return map(); + } + string parameter = arguments.at(i + j); + if (this->Internal->Debug) + { + MITK_DEBUG << " Processing parameter" << j << ", value:" << parameter; + } + if (this->argumentAdded(parameter)) + { +// this->Internal->ErrorString = +// missingParameterError.arg(argument).arg(j-1).arg(numberOfParametersToProcess); +// if (this->Internal->Debug) { MITK_DEBUG << this->Internal->ErrorString; } + if (ok) { *ok = false; } + return map(); + } + if (!currentArgDesc->addParameter(parameter)) + { +// this->Internal->ErrorString = string( +// "Value(s) associated with argument %1 are incorrect. %2"). +// arg(argument).arg(currentArgDesc->ExactMatchFailedMessage); +// if (this->Internal->Debug) { MITK_DEBUG << this->Internal->ErrorString; } + if (ok) { *ok = false; } + return map(); + } + } + // Update main loop increment + i = i + numberOfParametersToProcess; + } + else if (numberOfParametersToProcess == -1) + { + if (this->Internal->Debug) + { + MITK_DEBUG << " Proccessing StringList ..."; + } + int j = 1; + while(j + i < arguments.size()) + { + if (this->argumentAdded(arguments.at(j + i))) + { + if (this->Internal->Debug) + { + MITK_DEBUG << " No more parameter for" << argument; + } + break; + } + string parameter = arguments.at(j + i); + + if (parameter.compare(0, Internal->LongPrefix.size(), Internal->LongPrefix)==0 + || parameter.compare(0, Internal->ShortPrefix.size(), Internal->ShortPrefix)==0) + { + j--; + break; + } + + if (this->Internal->Debug) + { + MITK_DEBUG << " Processing parameter" << j << ", value:" << parameter; + } + if (!currentArgDesc->addParameter(parameter)) + { +// this->Internal->ErrorString = string( +// "Value(s) associated with argument %1 are incorrect. %2"). +// arg(argument).arg(currentArgDesc->ExactMatchFailedMessage); +// if (this->Internal->Debug) { MITK_DEBUG << this->Internal->ErrorString; } + if (ok) { *ok = false; } + return map(); + } + j++; + } + // Update main loop increment + i = i + j; + } + } + else + { + if (this->Internal->StrictMode) + { + this->Internal->ErrorString = "Unknown argument "; + this->Internal->ErrorString += argument; + error = true; + break; + } + if (this->Internal->Debug) + { + MITK_DEBUG << " Skipping: Unknown argument"; + } + this->Internal->UnparsedArguments.push_back(argument); + } + } + + if (ok) + { + *ok = !error; + } + + map parsedArguments; + + int obligatoryArgs = 0; + vector::iterator it; + for(it = Internal->ArgumentDescriptionList.begin(); it != Internal->ArgumentDescriptionList.end();++it) + { + CommandLineParserArgumentDescription* desc = *it; + + if(!desc->Optional) + obligatoryArgs++; + } + + int parsedObligatoryArgs = 0; + for(it = parsedArgDescriptions.begin(); it != parsedArgDescriptions.end();++it) + { + CommandLineParserArgumentDescription* desc = *it; + + string key; + if (!desc->LongArg.empty()) + { + key = desc->LongArg; + } + else + { + key = desc->ShortArg; + } + + if(!desc->Optional) + parsedObligatoryArgs++; + + std::pair elem; elem.first = key; elem.second = desc->Value; + parsedArguments.insert(elem); + } + + if (obligatoryArgs>parsedObligatoryArgs) + { + parsedArguments.clear(); + cout << helpText(); + } + + return parsedArguments; +} + +// ------------------------------------------------------------------------- +map ctkCommandLineParser::parseArguments(int argc, char** argv, bool* ok) +{ + StringContainerType arguments; + + // Create a StringContainerType of arguments + for(int i = 0; i < argc; ++i) + arguments.push_back(argv[i]); + + return this->parseArguments(arguments, ok); +} + +// ------------------------------------------------------------------------- +string ctkCommandLineParser::errorString() const +{ + return this->Internal->ErrorString; +} + +// ------------------------------------------------------------------------- +const ctkCommandLineParser::StringContainerType& ctkCommandLineParser::unparsedArguments() const +{ + return this->Internal->UnparsedArguments; +} + +// -------------------------------------------------------------------------- +void ctkCommandLineParser::addArgument(const string& longarg, const string& shortarg, + Type type, const string& argHelp, + const mitk::Any& defaultValue, bool optional, bool ignoreRest, + bool deprecated) +{ + if (longarg.empty() && shortarg.empty()) { return; } + + /* Make sure it's not already added */ + bool added = this->Internal->ArgNameToArgumentDescriptionMap.count(longarg); + if (added) { return; } + + added = this->Internal->ArgNameToArgumentDescriptionMap.count(shortarg); + if (added) { return; } + + CommandLineParserArgumentDescription* argDesc = + new CommandLineParserArgumentDescription(longarg, this->Internal->LongPrefix, + shortarg, this->Internal->ShortPrefix, type, + argHelp, defaultValue, ignoreRest, deprecated, optional); + + int argWidth = 0; + if (!longarg.empty()) + { + this->Internal->ArgNameToArgumentDescriptionMap[longarg] = argDesc; + argWidth += longarg.size() + this->Internal->LongPrefix.size(); + } + if (!shortarg.empty()) + { + this->Internal->ArgNameToArgumentDescriptionMap[shortarg] = argDesc; + argWidth += shortarg.size() + this->Internal->ShortPrefix.size() + 2; + } + argWidth += 5; + + // Set the field width for the arguments + if (argWidth > this->Internal->FieldWidth) + { + this->Internal->FieldWidth = argWidth; + } + + this->Internal->ArgumentDescriptionList.push_back(argDesc); + this->Internal->GroupToArgumentDescriptionListMap[this->Internal->CurrentGroup].push_back(argDesc); +} + +// -------------------------------------------------------------------------- +void ctkCommandLineParser::addDeprecatedArgument( + const string& longarg, const string& shortarg, const string& argHelp) +{ + addArgument(longarg, shortarg, StringList, argHelp, mitk::Any(), false, true); +} + +// -------------------------------------------------------------------------- +int ctkCommandLineParser::fieldWidth() const +{ + return this->Internal->FieldWidth; +} + +// -------------------------------------------------------------------------- +void ctkCommandLineParser::beginGroup(const string& description) +{ + this->Internal->CurrentGroup = description; +} + +// -------------------------------------------------------------------------- +void ctkCommandLineParser::endGroup() +{ + this->Internal->CurrentGroup.clear(); +} + +// -------------------------------------------------------------------------- +string ctkCommandLineParser::helpText() const +{ + string text; + vector deprecatedArgs; + + // Loop over grouped argument descriptions + map >::iterator it; + for(it = Internal->GroupToArgumentDescriptionListMap.begin(); it != Internal->GroupToArgumentDescriptionListMap.end();++it) + { + if (!(*it).first.empty()) + { + text = text + "\n" + (*it).first + "\n"; + } + + vector::iterator it2; + for(it2 = (*it).second.begin(); it2 != (*it).second.end(); ++it2) + { + CommandLineParserArgumentDescription* argDesc = *it2; + if (argDesc->Deprecated) + { + deprecatedArgs.push_back(argDesc); + } + else + { + text += argDesc->helpText(); + } + } + } + + if (!deprecatedArgs.empty()) + { + text += "\nDeprecated arguments:\n"; + vector::iterator it2; + for(it2 = deprecatedArgs.begin(); it2 != deprecatedArgs.end(); ++it2) + { + CommandLineParserArgumentDescription* argDesc = *it2; + text += argDesc->helpText(); + } + } + + return text; +} + +// -------------------------------------------------------------------------- +bool ctkCommandLineParser::argumentAdded(const string& argument) const +{ + return this->Internal->ArgNameToArgumentDescriptionMap.count(argument); +} + +// -------------------------------------------------------------------------- +bool ctkCommandLineParser::argumentParsed(const string& argument) const +{ + for (int i=0; iProcessedArguments.size(); i++) + if (argument.compare(Internal->ProcessedArguments.at(i))==0) + return true; + return false; +} + +// -------------------------------------------------------------------------- +void ctkCommandLineParser::setArgumentPrefix(const string& longPrefix, const string& shortPrefix) +{ + this->Internal->LongPrefix = longPrefix; + this->Internal->ShortPrefix = shortPrefix; +} + +// -------------------------------------------------------------------------- +void ctkCommandLineParser::setStrictModeEnabled(bool strictMode) +{ + this->Internal->StrictMode = strictMode; +} + diff --git a/Modules/DiffusionImaging/FiberTracking/MiniApps/ctkCommandLineParser.h b/Modules/DiffusionImaging/FiberTracking/MiniApps/ctkCommandLineParser.h new file mode 100755 index 0000000000..e7eb221379 --- /dev/null +++ b/Modules/DiffusionImaging/FiberTracking/MiniApps/ctkCommandLineParser.h @@ -0,0 +1,422 @@ +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +All rights reserved. + +This software is distributed WITHOUT ANY WARRANTY; without +even the implied warranty of MERCHANTABILITY or FITNESS FOR +A PARTICULAR PURPOSE. + +See LICENSE.txt or http://www.mitk.org for details. + +===================================================================*/ +/*========================================================================= + + Library: CTK + + Copyright (c) Kitware Inc. + + Licensed under the Apache License, Version 2.0 (the "License"); + you may not use this file except in compliance with the License. + You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0.txt + + Unless required by applicable law or agreed to in writing, software + distributed under the License is distributed on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + See the License for the specific language governing permissions and + limitations under the License. + +=========================================================================*/ + +#ifndef __ctkCommandLineParser_h +#define __ctkCommandLineParser_h + +#include +#include + + +/** + * \ingroup Core + * + * The CTK command line parser. + * + * Use this class to add information about the command line arguments + * your program understands and to easily parse them from a given list + * of strings. + * + * This parser provides the following features: + * + *
    + *
  • Add arguments by supplying a long name and/or a short name. + * Arguments are validated using a regular expression. They can have + * a default value and a help string.
    • + *
  • Deprecated arguments.
    • + *
  • Custom regular expressions for argument validation.
    • + *
  • Set different argument name prefixes for native platform look and feel.
    • + *
  • QSettings support. Default values for arguments can be read from + * a QSettings object.
    • + *
  • Create a help text for the command line arguments with support for + * grouping arguments.
    • + *
+ * + * Here is an example how to use this class inside a main function: + * + * \code + * #include + * #include + * #include + * + * int main(int argc, char** argv) + * { + * QCoreApplication app(argc, argv); + * // This is used by QSettings + * QCoreApplication::setOrganizationName("MyOrg"); + * QCoreApplication::setApplicationName("MyApp"); + * + * ctkCommandLineParser parser; + * // Use Unix-style argument names + * parser.setArgumentPrefix("--", "-"); + * // Enable QSettings support + * parser.enableSettings("disable-settings"); + * + * // Add command line argument names + * parser.addArgument("disable-settings", "", mitk::Any::Bool, "Do not use QSettings"); + * parser.addArgument("help", "h", mitk::Any::Bool, "Show this help text"); + * parser.addArgument("search-paths", "s", mitk::Any::StringList, "A list of paths to search"); + * + * // Parse the command line arguments + * bool ok = false; + * map parsedArgs = parser.parseArguments(QCoreApplication::arguments(), &ok); + * if (!ok) + * { + * QTextStream(stderr, QIODevice::WriteOnly) << "Error parsing arguments: " + * << parser.errorString() << "\n"; + * return EXIT_FAILURE; + * } + * + * // Show a help message + * if (parsedArgs.contains("help") || parsedArgs.contains("h")) + * { + * QTextStream(stdout, QIODevice::WriteOnly) << parser.helpText(); + * return EXIT_SUCCESS; + * } + * + * // Do something + * + * return EXIT_SUCCESS; + * } + * \endcode + */ + +using namespace std; + +class ctkCommandLineParser +{ + +public: + + enum Type { + String = 0, + Bool = 1, + StringList = 2, + Int = 3, + Float = 4 + }; + + typedef std::vector< std::string > StringContainerType; + + /** + * Constructs a parser instance. + * + * If QSettings support is enabled by a call to enableSettings() + * a default constructed QSettings instance will be used when parsing + * the command line arguments. Make sure to call QCoreApplication::setOrganizationName() + * and QCoreApplication::setApplicationName() before using default + * constructed QSettings objects. + * + * @param newParent The QObject parent. + */ + ctkCommandLineParser(); + + ~ctkCommandLineParser(); + + /** + * Parse a given list of command line arguments. + * + * This method parses a list of string elements considering the known arguments + * added by calls to addArgument(). If any one of the argument + * values does not match the corresponding regular expression, + * ok is set to false and an empty map object is returned. + * + * The keys in the returned map object correspond to the long argument string, + * if it is not empty. Otherwise, the short argument string is used as key. The + * mitk::Any values can safely be converted to the type specified in the + * addArgument() method call. + * + * @param arguments A StringContainerType containing command line arguments. Usually + * given by QCoreApplication::arguments(). + * @param ok A pointer to a boolean variable. Will be set to true + * if all regular expressions matched, false otherwise. + * @return A map object mapping the long argument (if empty, the short one) + * to a mitk::Any containing the value. + */ + + map parseArguments(const StringContainerType &arguments, bool* ok = 0); + + /** + * Convenient method allowing to parse a given list of command line arguments. + * @see parseArguments(const StringContainerType &, bool*) + */ + map parseArguments(int argc, char** argv, bool* ok = 0); + + /** + * Returns a detailed error description if a call to parseArguments() + * failed. + * + * @return The error description, empty if no error occured. + * @see parseArguments(const StringContainerType&, bool*) + */ + string errorString() const; + + /** + * This method returns all unparsed arguments, i.e. all arguments + * for which no long or short name has been registered via a call + * to addArgument(). + * + * @see addArgument() + * + * @return A list containing unparsed arguments. + */ + const StringContainerType& unparsedArguments() const; + + /** + * Checks if the given argument has been added via a call + * to addArgument(). + * + * @see addArgument() + * + * @param argument The argument to be checked. + * @return true if the argument was added, false + * otherwise. + */ + bool argumentAdded(const string& argument) const; + + /** + * Checks if the given argument has been parsed successfully by a previous + * call to parseArguments(). + * + * @param argument The argument to be checked. + * @return true if the argument was parsed, false + * otherwise. + */ + bool argumentParsed(const string& argument) const; + + /** + * Adds a command line argument. An argument can have a long name + * (like --long-argument-name), a short name (like -l), or both. The type + * of the argument can be specified by using the type parameter. + * The following types are supported: + * + * + * + * + * + * + * + * + *
Type# of parametersDefault regular exprExample
mitk::Any::String1.*--test-string StringParameter
mitk::Any::Bool0does not apply--enable-something
mitk::Any::StringList-1.*--test-list string1 string2
mitk::Any::Int1-?[0-9]+--test-int -5
+ * + * The regular expressions are used to validate the parameters of command line + * arguments. You can restrict the valid set of parameters by calling + * setExactMatchRegularExpression() for your argument. + * + * Optionally, a help string and a default value can be provided for the argument. If + * the mitk::Any type of the default value does not match type, an + * exception is thrown. Arguments with default values are always returned by + * parseArguments(). + * + * You can also declare an argument deprecated, by setting deprecated + * to true. Alternatively you can add a deprecated argument by calling + * addDeprecatedArgument(). + * + * If the long or short argument has already been added, or if both are empty strings, + * the method call has no effect. + * + * @param longarg The long argument name. + * @param shortarg The short argument name. + * @param type The argument type (see the list above for supported types). + * @param argHelp A help string describing the argument. + * @param defaultValue A default value for the argument. + * @param ignoreRest All arguments after the current one will be ignored. + * @param deprecated Declares the argument deprecated. + * + * @see setExactMatchRegularExpression() + * @see addDeprecatedArgument() + * @throws std::logic_error If the mitk::Any type of defaultValue + * does not match type, a std::logic_error is thrown. + */ + void addArgument(const string& longarg, const string& shortarg, + Type type, const string& argHelp = string(), + const mitk::Any& defaultValue = mitk::Any(), bool optional=true, + bool ignoreRest = false, bool deprecated = false); + + /** + * Adds a deprecated command line argument. If a deprecated argument is provided + * on the command line, argHelp is displayed in the console and + * processing continues with the next argument. + * + * Deprecated arguments are grouped separately at the end of the help text + * returned by helpText(). + * + * @param longarg The long argument name. + * @param shortarg The short argument name. + * @param argHelp A help string describing alternatives to the deprecated argument. + */ + void addDeprecatedArgument(const string& longarg, const string& shortarg, + const string& argHelp); + + /** + * Sets a custom regular expression for validating argument parameters. The method + * errorString() can be used the get the last error description. + * + * @param argument The previously added long or short argument name. + * @param expression A regular expression which the arugment parameters must match. + * @param exactMatchFailedMessage An error message explaining why the parameter did + * not match. + * + * @return true if the argument was found and the regular expression was set, + * false otherwise. + * + * @see errorString() + */ + bool setExactMatchRegularExpression(const string& argument, const string& expression, + const string& exactMatchFailedMessage); + + /** + * The field width for the argument names without the help text. + * + * @return The argument names field width in the help text. + */ + int fieldWidth() const; + + /** + * Creates a help text containing properly formatted argument names and help strings + * provided by calls to addArgument(). The arguments can be grouped by + * using beginGroup() and endGroup(). + * + * @param charPad The padding character. + * @return The formatted help text. + */ + string helpText() const; + + /** + * Sets the argument prefix for long and short argument names. This can be used + * to create native command line arguments without changing the calls to + * addArgument(). For example on Unix-based systems, long argument + * names start with "--" and short names with "-", while on Windows argument names + * always start with "/". + * + * Note that all methods in ctkCommandLineParser which take an argument name + * expect the name as it was supplied to addArgument. + * + * Example usage: + * + * \code + * ctkCommandLineParser parser; + * parser.setArgumentPrefix("--", "-"); + * parser.addArgument("long-argument", "l", mitk::Any::String); + * StringContainerType args; + * args << "program name" << "--long-argument Hi"; + * parser.parseArguments(args); + * \endcode + * + * @param longPrefix The prefix for long argument names. + * @param shortPrefix The prefix for short argument names. + */ + void setArgumentPrefix(const string& longPrefix, const string& shortPrefix); + + /** + * Begins a new group for documenting arguments. All newly added arguments via + * addArgument() will be put in the new group. You can close the + * current group by calling endGroup() or be opening a new group. + * + * Note that groups cannot be nested and all arguments which do not belong to + * a group will be listed at the top of the text created by helpText(). + * + * @param description The description of the group + */ + void beginGroup(const string& description); + + /** + * Ends the current group. + * + * @see beginGroup(const string&) + */ + void endGroup(); + + /** + * Enables QSettings support in ctkCommandLineParser. If an argument name is found + * in the QSettings instance with a valid mitk::Any, the value is considered as + * a default value and overwrites default values registered with + * addArgument(). User supplied values on the command line overwrite + * values in the QSettings instance, except for arguments with multiple parameters + * which are merged with QSettings values. Call mergeSettings(false) + * to disable merging. + * + * See ctkCommandLineParser(QSettings*) for information about how to + * supply a QSettings instance. + * + * Additionally, a long and short argument name can be specified which will disable + * QSettings support if supplied on the command line. The argument name must be + * registered as a regular argument via addArgument(). + * + * @param disableLongArg Long argument name. + * @param disableShortArg Short argument name. + * + * @see ctkCommandLineParser(QSettings*) + */ + void enableSettings(const string& disableLongArg = "", + const string& disableShortArg = ""); + + /** + * Controlls the merging behavior of user values and QSettings values. + * + * If merging is on (the default), user supplied values for an argument + * which can take more than one parameter are merged with values stored + * in the QSettings instance. If merging is off, the user values overwrite + * the QSettings values. + * + * @param merge true enables QSettings merging, false + * disables it. + */ + void mergeSettings(bool merge); + + /** + * Can be used to check if QSettings support has been enabled by a call to + * enableSettings(). + * + * @return true if QSettings support is enabled, false + * otherwise. + */ + bool settingsEnabled() const; + + + /** + * Can be used to teach the parser to stop parsing the arguments and return False when + * an unknown argument is encountered. By default StrictMode is disabled. + * + * @see parseArguments(const StringContainerType &, bool*) + */ + void setStrictModeEnabled(bool strictMode); + +private: + class ctkInternal; + ctkInternal * Internal; +}; + +#endif diff --git a/Modules/DiffusionImaging/FiberTracking/MiniApps/mitkFiberTrackingMiniApps.cpp b/Modules/DiffusionImaging/FiberTracking/MiniApps/mitkFiberTrackingMiniApps.cpp new file mode 100755 index 0000000000..c7d6ac601a --- /dev/null +++ b/Modules/DiffusionImaging/FiberTracking/MiniApps/mitkFiberTrackingMiniApps.cpp @@ -0,0 +1,22 @@ +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +All rights reserved. + +This software is distributed WITHOUT ANY WARRANTY; without +even the implied warranty of MERCHANTABILITY or FITNESS FOR +A PARTICULAR PURPOSE. + +See LICENSE.txt or http://www.mitk.org for details. + +===================================================================*/ + +#include + +int main(int argc, char* argv[]) +{ + return MiniAppManager::GetInstance()->RunMiniApp(argc, argv); +} diff --git a/Modules/DiffusionImaging/FiberTracking/SignalModels/mitkGibbsRingingArtifact.cpp b/Modules/DiffusionImaging/FiberTracking/SignalModels/mitkGibbsRingingArtifact.cpp deleted file mode 100644 index fc094304c2..0000000000 --- a/Modules/DiffusionImaging/FiberTracking/SignalModels/mitkGibbsRingingArtifact.cpp +++ /dev/null @@ -1,62 +0,0 @@ -/*=================================================================== - -The Medical Imaging Interaction Toolkit (MITK) - -Copyright (c) German Cancer Research Center, -Division of Medical and Biological Informatics. -All rights reserved. - -This software is distributed WITHOUT ANY WARRANTY; without -even the implied warranty of MERCHANTABILITY or FITNESS FOR -A PARTICULAR PURPOSE. - -See LICENSE.txt or http://www.mitk.org for details. - -===================================================================*/ - -#include - -template< class ScalarType > -GibbsRingingArtifact< ScalarType >::GibbsRingingArtifact() - : m_KspaceCropping(0.1) -{ - -} - -template< class ScalarType > -GibbsRingingArtifact< ScalarType >::~GibbsRingingArtifact() -{ - -} - -template< class ScalarType > -typename GibbsRingingArtifact< ScalarType >::ComplexSliceType::Pointer GibbsRingingArtifact< ScalarType >::AddArtifact(typename ComplexSliceType::Pointer slice) -{ - itk::ImageRegion<2> region = slice->GetLargestPossibleRegion(); - - int x = region.GetSize()[0]/m_KspaceCropping; - int y = region.GetSize()[1]/m_KspaceCropping; - - typename ComplexSliceType::Pointer newSlice = ComplexSliceType::New(); - itk::ImageRegion<2> newRegion; - newRegion.SetSize(0, x); - newRegion.SetSize(1, y); - - newSlice->SetLargestPossibleRegion( newRegion ); - newSlice->SetBufferedRegion( newRegion ); - newSlice->SetRequestedRegion( newRegion ); - newSlice->Allocate(); - - itk::ImageRegionIterator it(newSlice, newRegion); - while(!it.IsAtEnd()) - { - typename ComplexSliceType::IndexType idx; - idx[0] = region.GetSize()[0]/2-x/2+it.GetIndex()[0]; - idx[1] = region.GetSize()[1]/2-y/2+it.GetIndex()[1]; - - it.Set(slice->GetPixel(idx)); - - ++it; - } - return newSlice; -} diff --git a/Modules/DiffusionImaging/FiberTracking/SignalModels/mitkGibbsRingingArtifact.h b/Modules/DiffusionImaging/FiberTracking/SignalModels/mitkGibbsRingingArtifact.h deleted file mode 100644 index 9115346e66..0000000000 --- a/Modules/DiffusionImaging/FiberTracking/SignalModels/mitkGibbsRingingArtifact.h +++ /dev/null @@ -1,58 +0,0 @@ -/*=================================================================== - -The Medical Imaging Interaction Toolkit (MITK) - -Copyright (c) German Cancer Research Center, -Division of Medical and Biological Informatics. -All rights reserved. - -This software is distributed WITHOUT ANY WARRANTY; without -even the implied warranty of MERCHANTABILITY or FITNESS FOR -A PARTICULAR PURPOSE. - -See LICENSE.txt or http://www.mitk.org for details. - -===================================================================*/ - -#ifndef _MITK_GibbsRingingArtifact_H -#define _MITK_GibbsRingingArtifact_H - -#include -#include -#include -#include - -namespace mitk { - -/** - * \brief Class to add gibbs ringing artifact to input complex slice - * - */ -template< class ScalarType > -class GibbsRingingArtifact : public KspaceArtifact< ScalarType > -{ -public: - - GibbsRingingArtifact(); - ~GibbsRingingArtifact(); - - typedef typename KspaceArtifact< ScalarType >::ComplexSliceType ComplexSliceType; - - /** Adds Gibbs ringing to the input slice (by zero padding). **/ - typename ComplexSliceType::Pointer AddArtifact(typename ComplexSliceType::Pointer slice); - - void SetKspaceCropping(int factor){ m_KspaceCropping=factor; } - int GetKspaceCropping(){ return m_KspaceCropping; } - -protected: - - int m_KspaceCropping; - -}; - -#include "mitkGibbsRingingArtifact.cpp" - -} - -#endif - diff --git a/Modules/DiffusionImaging/FiberTracking/files.cmake b/Modules/DiffusionImaging/FiberTracking/files.cmake index 1792d5ee6d..912c88e510 100644 --- a/Modules/DiffusionImaging/FiberTracking/files.cmake +++ b/Modules/DiffusionImaging/FiberTracking/files.cmake @@ -1,101 +1,105 @@ set(CPP_FILES -# DataStructures -> FiberBundleX + MiniApps/ctkCommandLineParser.h + + # DataStructures -> FiberBundleX IODataStructures/FiberBundleX/mitkFiberBundleX.cpp IODataStructures/FiberBundleX/mitkFiberBundleXWriter.cpp IODataStructures/FiberBundleX/mitkFiberBundleXReader.cpp IODataStructures/FiberBundleX/mitkFiberBundleXIOFactory.cpp IODataStructures/FiberBundleX/mitkFiberBundleXWriterFactory.cpp IODataStructures/FiberBundleX/mitkFiberBundleXSerializer.cpp IODataStructures/FiberBundleX/mitkFiberBundleXThreadMonitor.cpp # DataStructures -> PlanarFigureComposite IODataStructures/PlanarFigureComposite/mitkPlanarFigureComposite.cpp # DataStructures IODataStructures/mitkFiberTrackingObjectFactory.cpp # Rendering Rendering/mitkFiberBundleXMapper2D.cpp Rendering/mitkFiberBundleXMapper3D.cpp Rendering/mitkFiberBundleXThreadMonitorMapper3D.cpp #Rendering/mitkPlanarFigureMapper3D.cpp # Interactions Interactions/mitkFiberBundleInteractor.cpp # Tractography Algorithms/GibbsTracking/mitkParticleGrid.cpp Algorithms/GibbsTracking/mitkMetropolisHastingsSampler.cpp Algorithms/GibbsTracking/mitkEnergyComputer.cpp Algorithms/GibbsTracking/mitkGibbsEnergyComputer.cpp Algorithms/GibbsTracking/mitkFiberBuilder.cpp Algorithms/GibbsTracking/mitkSphereInterpolator.cpp ) set(H_FILES + MiniApps/ctkCommandLineParser.h + # Rendering Rendering/mitkFiberBundleXMapper3D.h Rendering/mitkFiberBundleXMapper2D.h Rendering/mitkFiberBundleXThreadMonitorMapper3D.h #Rendering/mitkPlanarFigureMapper3D.h # DataStructures -> FiberBundleX IODataStructures/FiberBundleX/mitkFiberBundleX.h IODataStructures/FiberBundleX/mitkFiberBundleXWriter.h IODataStructures/FiberBundleX/mitkFiberBundleXReader.h IODataStructures/FiberBundleX/mitkFiberBundleXIOFactory.h IODataStructures/FiberBundleX/mitkFiberBundleXWriterFactory.h IODataStructures/FiberBundleX/mitkFiberBundleXSerializer.h IODataStructures/FiberBundleX/mitkFiberBundleXThreadMonitor.h IODataStructures/mitkFiberTrackingObjectFactory.h # Algorithms Algorithms/itkTractDensityImageFilter.h Algorithms/itkTractsToFiberEndingsImageFilter.h Algorithms/itkTractsToRgbaImageFilter.h Algorithms/itkElectrostaticRepulsionDiffusionGradientReductionFilter.h Algorithms/itkFibersFromPlanarFiguresFilter.h Algorithms/itkTractsToDWIImageFilter.h Algorithms/itkTractsToVectorImageFilter.h Algorithms/itkKspaceImageFilter.h Algorithms/itkDftImageFilter.h Algorithms/itkAddArtifactsToDwiImageFilter.h Algorithms/itkFieldmapGeneratorFilter.h + Algorithms/itkEvaluateDirectionImagesFilter.h # (old) Tractography Algorithms/itkGibbsTrackingFilter.h Algorithms/itkStochasticTractographyFilter.h Algorithms/itkStreamlineTrackingFilter.h Algorithms/GibbsTracking/mitkParticle.h Algorithms/GibbsTracking/mitkParticleGrid.h Algorithms/GibbsTracking/mitkMetropolisHastingsSampler.h Algorithms/GibbsTracking/mitkSimpSamp.h Algorithms/GibbsTracking/mitkEnergyComputer.h Algorithms/GibbsTracking/mitkGibbsEnergyComputer.h Algorithms/GibbsTracking/mitkSphereInterpolator.h Algorithms/GibbsTracking/mitkFiberBuilder.h # Signal Models SignalModels/mitkDiffusionSignalModel.h SignalModels/mitkTensorModel.h SignalModels/mitkBallModel.h SignalModels/mitkDotModel.h SignalModels/mitkAstroStickModel.h SignalModels/mitkStickModel.h SignalModels/mitkDiffusionNoiseModel.h SignalModels/mitkRicianNoiseModel.h SignalModels/mitkKspaceArtifact.h - SignalModels/mitkGibbsRingingArtifact.h ) set(RESOURCE_FILES # Binary directory resources FiberTrackingLUTBaryCoords.bin FiberTrackingLUTIndices.bin # Shaders Shaders/mitkShaderFiberClipping.xml ) diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/QmitkZeppelinModelParametersWidgetControls.ui b/Plugins/org.mitk.gui.qt.diffusionimaging/src/QmitkZeppelinModelParametersWidgetControls.ui index 418098a551..ceca723203 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/QmitkZeppelinModelParametersWidgetControls.ui +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/QmitkZeppelinModelParametersWidgetControls.ui @@ -1,124 +1,124 @@ QmitkZeppelinModelParametersWidgetControls 0 0 410 116 0 0 Form 0 Fractional anisotropy of resulting tensor. <html><head/><body><p>-</p></body></html> Diffusivity along second and third eigenvector. 5 1.000000000000000 0.000100000000000 0.000250000000000 T2 relaxation time of this compartment (in milliseconds). 10000 - 100 + 200 Diffusivity along largest eigenvector. 5 1.000000000000000 0.000100000000000 0.001000000000000 <html><head/><body><p><span style=" font-style:italic;">d</span><span style=" vertical-align:sub;">⟂</span>:</p></body></html> <html><head/><body><p><span style=" font-style:italic;">d</span><span style=" vertical-align:sub;">||</span>:</p></body></html> <html><head/><body><p><span style=" font-style:italic;">T2</span> - Relaxation:</p></body></html> <html><head/><body><p><span style=" font-style:italic;">FA</span>:</p></body></html> m_T2box m_D1box m_D2box diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberfoxView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberfoxView.cpp old mode 100644 new mode 100755 index 1ced21f898..ff0bd67487 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberfoxView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberfoxView.cpp @@ -1,1794 +1,1789 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ //misc #define _USE_MATH_DEFINES #include // Blueberry #include #include // Qmitk #include "QmitkFiberfoxView.h" // MITK #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include -#include #include #include #include #include #include #include #include #define _USE_MATH_DEFINES #include const std::string QmitkFiberfoxView::VIEW_ID = "org.mitk.views.fiberfoxview"; QmitkFiberfoxView::QmitkFiberfoxView() : QmitkAbstractView() , m_Controls( 0 ) , m_SelectedImage( NULL ) { } // Destructor QmitkFiberfoxView::~QmitkFiberfoxView() { } void QmitkFiberfoxView::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::QmitkFiberfoxViewControls; m_Controls->setupUi( parent ); m_Controls->m_StickWidget1->setVisible(true); m_Controls->m_StickWidget2->setVisible(false); m_Controls->m_ZeppelinWidget1->setVisible(false); m_Controls->m_ZeppelinWidget2->setVisible(false); m_Controls->m_TensorWidget1->setVisible(false); m_Controls->m_TensorWidget2->setVisible(false); m_Controls->m_BallWidget1->setVisible(true); m_Controls->m_BallWidget2->setVisible(false); m_Controls->m_AstrosticksWidget1->setVisible(false); m_Controls->m_AstrosticksWidget2->setVisible(false); m_Controls->m_DotWidget1->setVisible(false); m_Controls->m_DotWidget2->setVisible(false); m_Controls->m_Comp4FractionFrame->setVisible(false); m_Controls->m_DiffusionPropsMessage->setVisible(false); m_Controls->m_GeometryMessage->setVisible(false); m_Controls->m_AdvancedSignalOptionsFrame->setVisible(false); m_Controls->m_AdvancedFiberOptionsFrame->setVisible(false); m_Controls->m_VarianceBox->setVisible(false); m_Controls->m_GibbsRingingFrame->setVisible(false); m_Controls->m_NoiseFrame->setVisible(false); m_Controls->m_GhostFrame->setVisible(false); m_Controls->m_DistortionsFrame->setVisible(false); m_Controls->m_EddyFrame->setVisible(false); m_Controls->m_FrequencyMapBox->SetDataStorage(this->GetDataStorage()); mitk::TNodePredicateDataType::Pointer isMitkImage = mitk::TNodePredicateDataType::New(); mitk::NodePredicateDataType::Pointer isDwi = mitk::NodePredicateDataType::New("DiffusionImage"); mitk::NodePredicateDataType::Pointer isDti = mitk::NodePredicateDataType::New("TensorImage"); mitk::NodePredicateDataType::Pointer isQbi = mitk::NodePredicateDataType::New("QBallImage"); mitk::NodePredicateOr::Pointer isDiffusionImage = mitk::NodePredicateOr::New(isDwi, isDti); isDiffusionImage = mitk::NodePredicateOr::New(isDiffusionImage, isQbi); mitk::NodePredicateNot::Pointer noDiffusionImage = mitk::NodePredicateNot::New(isDiffusionImage); mitk::NodePredicateAnd::Pointer finalPredicate = mitk::NodePredicateAnd::New(isMitkImage, noDiffusionImage); m_Controls->m_FrequencyMapBox->SetPredicate(finalPredicate); connect((QObject*) m_Controls->m_GenerateImageButton, SIGNAL(clicked()), (QObject*) this, SLOT(GenerateImage())); connect((QObject*) m_Controls->m_GenerateFibersButton, SIGNAL(clicked()), (QObject*) this, SLOT(GenerateFibers())); connect((QObject*) m_Controls->m_CircleButton, SIGNAL(clicked()), (QObject*) this, SLOT(OnDrawROI())); connect((QObject*) m_Controls->m_FlipButton, SIGNAL(clicked()), (QObject*) this, SLOT(OnFlipButton())); connect((QObject*) m_Controls->m_JoinBundlesButton, SIGNAL(clicked()), (QObject*) this, SLOT(JoinBundles())); connect((QObject*) m_Controls->m_VarianceBox, SIGNAL(valueChanged(double)), (QObject*) this, SLOT(OnVarianceChanged(double))); connect((QObject*) m_Controls->m_DistributionBox, SIGNAL(currentIndexChanged(int)), (QObject*) this, SLOT(OnDistributionChanged(int))); connect((QObject*) m_Controls->m_FiberDensityBox, SIGNAL(valueChanged(int)), (QObject*) this, SLOT(OnFiberDensityChanged(int))); connect((QObject*) m_Controls->m_FiberSamplingBox, SIGNAL(valueChanged(int)), (QObject*) this, SLOT(OnFiberSamplingChanged(int))); connect((QObject*) m_Controls->m_TensionBox, SIGNAL(valueChanged(double)), (QObject*) this, SLOT(OnTensionChanged(double))); connect((QObject*) m_Controls->m_ContinuityBox, SIGNAL(valueChanged(double)), (QObject*) this, SLOT(OnContinuityChanged(double))); connect((QObject*) m_Controls->m_BiasBox, SIGNAL(valueChanged(double)), (QObject*) this, SLOT(OnBiasChanged(double))); connect((QObject*) m_Controls->m_AddGibbsRinging, SIGNAL(stateChanged(int)), (QObject*) this, SLOT(OnAddGibbsRinging(int))); connect((QObject*) m_Controls->m_AddNoise, SIGNAL(stateChanged(int)), (QObject*) this, SLOT(OnAddNoise(int))); connect((QObject*) m_Controls->m_AddGhosts, SIGNAL(stateChanged(int)), (QObject*) this, SLOT(OnAddGhosts(int))); connect((QObject*) m_Controls->m_AddDistortions, SIGNAL(stateChanged(int)), (QObject*) this, SLOT(OnAddDistortions(int))); connect((QObject*) m_Controls->m_AddEddy, SIGNAL(stateChanged(int)), (QObject*) this, SLOT(OnAddEddy(int))); connect((QObject*) m_Controls->m_ConstantRadiusBox, SIGNAL(stateChanged(int)), (QObject*) this, SLOT(OnConstantRadius(int))); connect((QObject*) m_Controls->m_CopyBundlesButton, SIGNAL(clicked()), (QObject*) this, SLOT(CopyBundles())); connect((QObject*) m_Controls->m_TransformBundlesButton, SIGNAL(clicked()), (QObject*) this, SLOT(ApplyTransform())); connect((QObject*) m_Controls->m_AlignOnGrid, SIGNAL(clicked()), (QObject*) this, SLOT(AlignOnGrid())); connect((QObject*) m_Controls->m_Compartment1Box, SIGNAL(currentIndexChanged(int)), (QObject*) this, SLOT(Comp1ModelFrameVisibility(int))); connect((QObject*) m_Controls->m_Compartment2Box, SIGNAL(currentIndexChanged(int)), (QObject*) this, SLOT(Comp2ModelFrameVisibility(int))); connect((QObject*) m_Controls->m_Compartment3Box, SIGNAL(currentIndexChanged(int)), (QObject*) this, SLOT(Comp3ModelFrameVisibility(int))); connect((QObject*) m_Controls->m_Compartment4Box, SIGNAL(currentIndexChanged(int)), (QObject*) this, SLOT(Comp4ModelFrameVisibility(int))); connect((QObject*) m_Controls->m_AdvancedOptionsBox, SIGNAL( stateChanged(int)), (QObject*) this, SLOT(ShowAdvancedOptions(int))); connect((QObject*) m_Controls->m_AdvancedOptionsBox_2, SIGNAL( stateChanged(int)), (QObject*) this, SLOT(ShowAdvancedOptions(int))); } } void QmitkFiberfoxView::ShowAdvancedOptions(int state) { if (state) { m_Controls->m_AdvancedFiberOptionsFrame->setVisible(true); m_Controls->m_AdvancedSignalOptionsFrame->setVisible(true); m_Controls->m_AdvancedOptionsBox->setChecked(true); m_Controls->m_AdvancedOptionsBox_2->setChecked(true); } else { m_Controls->m_AdvancedFiberOptionsFrame->setVisible(false); m_Controls->m_AdvancedSignalOptionsFrame->setVisible(false); m_Controls->m_AdvancedOptionsBox->setChecked(false); m_Controls->m_AdvancedOptionsBox_2->setChecked(false); } } void QmitkFiberfoxView::Comp1ModelFrameVisibility(int index) { m_Controls->m_StickWidget1->setVisible(false); m_Controls->m_ZeppelinWidget1->setVisible(false); m_Controls->m_TensorWidget1->setVisible(false); switch (index) { case 0: m_Controls->m_StickWidget1->setVisible(true); break; case 1: m_Controls->m_ZeppelinWidget1->setVisible(true); break; case 2: m_Controls->m_TensorWidget1->setVisible(true); break; } } void QmitkFiberfoxView::Comp2ModelFrameVisibility(int index) { m_Controls->m_StickWidget2->setVisible(false); m_Controls->m_ZeppelinWidget2->setVisible(false); m_Controls->m_TensorWidget2->setVisible(false); switch (index) { case 0: break; case 1: m_Controls->m_StickWidget2->setVisible(true); break; case 2: m_Controls->m_ZeppelinWidget2->setVisible(true); break; case 3: m_Controls->m_TensorWidget2->setVisible(true); break; } } void QmitkFiberfoxView::Comp3ModelFrameVisibility(int index) { m_Controls->m_BallWidget1->setVisible(false); m_Controls->m_AstrosticksWidget1->setVisible(false); m_Controls->m_DotWidget1->setVisible(false); switch (index) { case 0: m_Controls->m_BallWidget1->setVisible(true); break; case 1: m_Controls->m_AstrosticksWidget1->setVisible(true); break; case 2: m_Controls->m_DotWidget1->setVisible(true); break; } } void QmitkFiberfoxView::Comp4ModelFrameVisibility(int index) { m_Controls->m_BallWidget2->setVisible(false); m_Controls->m_AstrosticksWidget2->setVisible(false); m_Controls->m_DotWidget2->setVisible(false); m_Controls->m_Comp4FractionFrame->setVisible(false); switch (index) { case 0: break; case 1: m_Controls->m_BallWidget2->setVisible(true); m_Controls->m_Comp4FractionFrame->setVisible(true); break; case 2: m_Controls->m_AstrosticksWidget2->setVisible(true); m_Controls->m_Comp4FractionFrame->setVisible(true); break; case 3: m_Controls->m_DotWidget2->setVisible(true); m_Controls->m_Comp4FractionFrame->setVisible(true); break; } } void QmitkFiberfoxView::OnConstantRadius(int value) { if (value>0 && m_Controls->m_RealTimeFibers->isChecked()) GenerateFibers(); } void QmitkFiberfoxView::OnAddEddy(int value) { if (value>0) m_Controls->m_EddyFrame->setVisible(true); else m_Controls->m_EddyFrame->setVisible(false); } void QmitkFiberfoxView::OnAddDistortions(int value) { if (value>0) m_Controls->m_DistortionsFrame->setVisible(true); else m_Controls->m_DistortionsFrame->setVisible(false); } void QmitkFiberfoxView::OnAddGhosts(int value) { if (value>0) m_Controls->m_GhostFrame->setVisible(true); else m_Controls->m_GhostFrame->setVisible(false); } void QmitkFiberfoxView::OnAddNoise(int value) { if (value>0) m_Controls->m_NoiseFrame->setVisible(true); else m_Controls->m_NoiseFrame->setVisible(false); } void QmitkFiberfoxView::OnAddGibbsRinging(int value) { if (value>0) m_Controls->m_GibbsRingingFrame->setVisible(true); else m_Controls->m_GibbsRingingFrame->setVisible(false); } void QmitkFiberfoxView::OnDistributionChanged(int value) { if (value==1) m_Controls->m_VarianceBox->setVisible(true); else m_Controls->m_VarianceBox->setVisible(false); if (m_Controls->m_RealTimeFibers->isChecked()) GenerateFibers(); } void QmitkFiberfoxView::OnVarianceChanged(double value) { if (m_Controls->m_RealTimeFibers->isChecked()) GenerateFibers(); } void QmitkFiberfoxView::OnFiberDensityChanged(int value) { if (m_Controls->m_RealTimeFibers->isChecked()) GenerateFibers(); } void QmitkFiberfoxView::OnFiberSamplingChanged(int value) { if (m_Controls->m_RealTimeFibers->isChecked()) GenerateFibers(); } void QmitkFiberfoxView::OnTensionChanged(double value) { if (m_Controls->m_RealTimeFibers->isChecked()) GenerateFibers(); } void QmitkFiberfoxView::OnContinuityChanged(double value) { if (m_Controls->m_RealTimeFibers->isChecked()) GenerateFibers(); } void QmitkFiberfoxView::OnBiasChanged(double value) { if (m_Controls->m_RealTimeFibers->isChecked()) GenerateFibers(); } void QmitkFiberfoxView::AlignOnGrid() { for (int i=0; i(m_SelectedFiducials.at(i)->GetData()); mitk::Point3D wc0 = pe->GetWorldControlPoint(0); mitk::DataStorage::SetOfObjects::ConstPointer parentFibs = GetDataStorage()->GetSources(m_SelectedFiducials.at(i)); for( mitk::DataStorage::SetOfObjects::const_iterator it = parentFibs->begin(); it != parentFibs->end(); ++it ) { mitk::DataNode::Pointer pFibNode = *it; if ( pFibNode.IsNotNull() && dynamic_cast(pFibNode->GetData()) ) { mitk::DataStorage::SetOfObjects::ConstPointer parentImgs = GetDataStorage()->GetSources(pFibNode); for( mitk::DataStorage::SetOfObjects::const_iterator it2 = parentImgs->begin(); it2 != parentImgs->end(); ++it2 ) { mitk::DataNode::Pointer pImgNode = *it2; if ( pImgNode.IsNotNull() && dynamic_cast(pImgNode->GetData()) ) { mitk::Image::Pointer img = dynamic_cast(pImgNode->GetData()); mitk::Geometry3D::Pointer geom = img->GetGeometry(); itk::Index<3> idx; geom->WorldToIndex(wc0, idx); mitk::Point3D cIdx; cIdx[0]=idx[0]; cIdx[1]=idx[1]; cIdx[2]=idx[2]; mitk::Point3D world; geom->IndexToWorld(cIdx,world); mitk::Vector3D trans = world - wc0; pe->GetGeometry()->Translate(trans); break; } } break; } } } for( int i=0; iGetSources(fibNode); for( mitk::DataStorage::SetOfObjects::const_iterator it = sources->begin(); it != sources->end(); ++it ) { mitk::DataNode::Pointer imgNode = *it; if ( imgNode.IsNotNull() && dynamic_cast(imgNode->GetData()) ) { mitk::DataStorage::SetOfObjects::ConstPointer derivations = GetDataStorage()->GetDerivations(fibNode); for( mitk::DataStorage::SetOfObjects::const_iterator it2 = derivations->begin(); it2 != derivations->end(); ++it2 ) { mitk::DataNode::Pointer fiducialNode = *it2; if ( fiducialNode.IsNotNull() && dynamic_cast(fiducialNode->GetData()) ) { mitk::PlanarEllipse::Pointer pe = dynamic_cast(fiducialNode->GetData()); mitk::Point3D wc0 = pe->GetWorldControlPoint(0); mitk::Image::Pointer img = dynamic_cast(imgNode->GetData()); mitk::Geometry3D::Pointer geom = img->GetGeometry(); itk::Index<3> idx; geom->WorldToIndex(wc0, idx); mitk::Point3D cIdx; cIdx[0]=idx[0]; cIdx[1]=idx[1]; cIdx[2]=idx[2]; mitk::Point3D world; geom->IndexToWorld(cIdx,world); mitk::Vector3D trans = world - wc0; pe->GetGeometry()->Translate(trans); } } break; } } } for( int i=0; i(m_SelectedImages.at(i)->GetData()); mitk::DataStorage::SetOfObjects::ConstPointer derivations = GetDataStorage()->GetDerivations(m_SelectedImages.at(i)); for( mitk::DataStorage::SetOfObjects::const_iterator it = derivations->begin(); it != derivations->end(); ++it ) { mitk::DataNode::Pointer fibNode = *it; if ( fibNode.IsNotNull() && dynamic_cast(fibNode->GetData()) ) { mitk::DataStorage::SetOfObjects::ConstPointer derivations2 = GetDataStorage()->GetDerivations(fibNode); for( mitk::DataStorage::SetOfObjects::const_iterator it2 = derivations2->begin(); it2 != derivations2->end(); ++it2 ) { mitk::DataNode::Pointer fiducialNode = *it2; if ( fiducialNode.IsNotNull() && dynamic_cast(fiducialNode->GetData()) ) { mitk::PlanarEllipse::Pointer pe = dynamic_cast(fiducialNode->GetData()); mitk::Point3D wc0 = pe->GetWorldControlPoint(0); mitk::Geometry3D::Pointer geom = img->GetGeometry(); itk::Index<3> idx; geom->WorldToIndex(wc0, idx); mitk::Point3D cIdx; cIdx[0]=idx[0]; cIdx[1]=idx[1]; cIdx[2]=idx[2]; mitk::Point3D world; geom->IndexToWorld(cIdx,world); mitk::Vector3D trans = world - wc0; pe->GetGeometry()->Translate(trans); } } } } } mitk::RenderingManager::GetInstance()->RequestUpdateAll(); if (m_Controls->m_RealTimeFibers->isChecked()) GenerateFibers(); } void QmitkFiberfoxView::OnFlipButton() { if (m_SelectedFiducial.IsNull()) return; std::map::iterator it = m_DataNodeToPlanarFigureData.find(m_SelectedFiducial.GetPointer()); if( it != m_DataNodeToPlanarFigureData.end() ) { QmitkPlanarFigureData& data = it->second; data.m_Flipped += 1; data.m_Flipped %= 2; } if (m_Controls->m_RealTimeFibers->isChecked()) GenerateFibers(); } QmitkFiberfoxView::GradientListType QmitkFiberfoxView::GenerateHalfShell(int NPoints) { NPoints *= 2; GradientListType pointshell; int numB0 = NPoints/20; if (numB0==0) numB0=1; GradientType g; g.Fill(0.0); for (int i=0; i theta; theta.set_size(NPoints); vnl_vector phi; phi.set_size(NPoints); double C = sqrt(4*M_PI); phi(0) = 0.0; phi(NPoints-1) = 0.0; for(int i=0; i0 && i std::vector > QmitkFiberfoxView::MakeGradientList() { std::vector > retval; vnl_matrix_fixed* U = itk::PointShell >::DistributePointShell(); // Add 0 vector for B0 int numB0 = ndirs/10; if (numB0==0) numB0=1; itk::Vector v; v.Fill(0.0); for (int i=0; i v; v[0] = U->get(0,i); v[1] = U->get(1,i); v[2] = U->get(2,i); retval.push_back(v); } return retval; } void QmitkFiberfoxView::OnAddBundle() { if (m_SelectedImage.IsNull()) return; mitk::DataStorage::SetOfObjects::ConstPointer children = GetDataStorage()->GetDerivations(m_SelectedImage); mitk::FiberBundleX::Pointer bundle = mitk::FiberBundleX::New(); mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData( bundle ); QString name = QString("Bundle_%1").arg(children->size()); node->SetName(name.toStdString()); m_SelectedBundles.push_back(node); UpdateGui(); GetDataStorage()->Add(node, m_SelectedImage); } void QmitkFiberfoxView::OnDrawROI() { if (m_SelectedBundles.empty()) OnAddBundle(); if (m_SelectedBundles.empty()) return; mitk::DataStorage::SetOfObjects::ConstPointer children = GetDataStorage()->GetDerivations(m_SelectedBundles.at(0)); mitk::PlanarEllipse::Pointer figure = mitk::PlanarEllipse::New(); mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData( figure ); QList nodes = this->GetDataManagerSelection(); for( int i=0; iSetSelected(false); m_SelectedFiducial = node; QString name = QString("Fiducial_%1").arg(children->size()); node->SetName(name.toStdString()); node->SetSelected(true); GetDataStorage()->Add(node, m_SelectedBundles.at(0)); this->DisableCrosshairNavigation(); mitk::PlanarFigureInteractor::Pointer figureInteractor = dynamic_cast(node->GetInteractor()); if(figureInteractor.IsNull()) figureInteractor = mitk::PlanarFigureInteractor::New("PlanarFigureInteractor", node); mitk::GlobalInteraction::GetInstance()->AddInteractor(figureInteractor); UpdateGui(); } bool CompareLayer(mitk::DataNode::Pointer i,mitk::DataNode::Pointer j) { int li = -1; i->GetPropertyValue("layer", li); int lj = -1; j->GetPropertyValue("layer", lj); return liGetSources(m_SelectedFiducial); for( mitk::DataStorage::SetOfObjects::const_iterator it = parents->begin(); it != parents->end(); ++it ) if(dynamic_cast((*it)->GetData())) m_SelectedBundles.push_back(*it); if (m_SelectedBundles.empty()) return; } vector< vector< mitk::PlanarEllipse::Pointer > > fiducials; vector< vector< unsigned int > > fliplist; for (int i=0; iGetDerivations(m_SelectedBundles.at(i)); std::vector< mitk::DataNode::Pointer > childVector; for( mitk::DataStorage::SetOfObjects::const_iterator it = children->begin(); it != children->end(); ++it ) childVector.push_back(*it); sort(childVector.begin(), childVector.end(), CompareLayer); vector< mitk::PlanarEllipse::Pointer > fib; vector< unsigned int > flip; float radius = 1; int count = 0; for( std::vector< mitk::DataNode::Pointer >::const_iterator it = childVector.begin(); it != childVector.end(); ++it ) { mitk::DataNode::Pointer node = *it; if ( node.IsNotNull() && dynamic_cast(node->GetData()) ) { mitk::PlanarEllipse* ellipse = dynamic_cast(node->GetData()); if (m_Controls->m_ConstantRadiusBox->isChecked()) { ellipse->SetTreatAsCircle(true); mitk::Point2D c = ellipse->GetControlPoint(0); mitk::Point2D p = ellipse->GetControlPoint(1); mitk::Vector2D v = p-c; if (count==0) { radius = v.GetVnlVector().magnitude(); ellipse->SetControlPoint(1, p); } else { v.Normalize(); v *= radius; ellipse->SetControlPoint(1, c+v); } } fib.push_back(ellipse); std::map::iterator it = m_DataNodeToPlanarFigureData.find(node.GetPointer()); if( it != m_DataNodeToPlanarFigureData.end() ) { QmitkPlanarFigureData& data = it->second; flip.push_back(data.m_Flipped); } else flip.push_back(0); } count++; } if (fib.size()>1) { fiducials.push_back(fib); fliplist.push_back(flip); } else if (fib.size()>0) m_SelectedBundles.at(i)->SetData( mitk::FiberBundleX::New() ); mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } itk::FibersFromPlanarFiguresFilter::Pointer filter = itk::FibersFromPlanarFiguresFilter::New(); filter->SetFiducials(fiducials); filter->SetFlipList(fliplist); switch(m_Controls->m_DistributionBox->currentIndex()){ case 0: filter->SetFiberDistribution(itk::FibersFromPlanarFiguresFilter::DISTRIBUTE_UNIFORM); break; case 1: filter->SetFiberDistribution(itk::FibersFromPlanarFiguresFilter::DISTRIBUTE_GAUSSIAN); filter->SetVariance(m_Controls->m_VarianceBox->value()); break; } filter->SetDensity(m_Controls->m_FiberDensityBox->value()); filter->SetTension(m_Controls->m_TensionBox->value()); filter->SetContinuity(m_Controls->m_ContinuityBox->value()); filter->SetBias(m_Controls->m_BiasBox->value()); filter->SetFiberSampling(m_Controls->m_FiberSamplingBox->value()); filter->Update(); vector< mitk::FiberBundleX::Pointer > fiberBundles = filter->GetFiberBundles(); for (int i=0; iSetData( fiberBundles.at(i) ); if (fiberBundles.at(i)->GetNumFibers()>50000) m_SelectedBundles.at(i)->SetVisibility(false); } mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberfoxView::GenerateImage() { itk::ImageRegion<3> imageRegion; imageRegion.SetSize(0, m_Controls->m_SizeX->value()); imageRegion.SetSize(1, m_Controls->m_SizeY->value()); imageRegion.SetSize(2, m_Controls->m_SizeZ->value()); mitk::Vector3D spacing; spacing[0] = m_Controls->m_SpacingX->value(); spacing[1] = m_Controls->m_SpacingY->value(); spacing[2] = m_Controls->m_SpacingZ->value(); mitk::Point3D origin; origin[0] = spacing[0]/2; origin[1] = spacing[1]/2; origin[2] = spacing[2]/2; itk::Matrix directionMatrix; directionMatrix.SetIdentity(); if (m_SelectedImage.IsNotNull()) { mitk::Image* img = dynamic_cast(m_SelectedImage->GetData()); itk::Image< float, 3 >::Pointer itkImg = itk::Image< float, 3 >::New(); CastToItkImage< itk::Image< float, 3 > >(img, itkImg); imageRegion = itkImg->GetLargestPossibleRegion(); spacing = itkImg->GetSpacing(); origin = itkImg->GetOrigin(); directionMatrix = itkImg->GetDirection(); } DiffusionSignalModel::GradientListType gradientList; double bVal = 1000; if (m_SelectedDWI.IsNull()) { gradientList = GenerateHalfShell(m_Controls->m_NumGradientsBox->value());; bVal = m_Controls->m_BvalueBox->value(); } else { mitk::DiffusionImage::Pointer dwi = dynamic_cast*>(m_SelectedDWI->GetData()); imageRegion = dwi->GetVectorImage()->GetLargestPossibleRegion(); spacing = dwi->GetVectorImage()->GetSpacing(); origin = dwi->GetVectorImage()->GetOrigin(); directionMatrix = dwi->GetVectorImage()->GetDirection(); bVal = dwi->GetB_Value(); mitk::DiffusionImage::GradientDirectionContainerType::Pointer dirs = dwi->GetDirections(); for (int i=0; iSize(); i++) { DiffusionSignalModel::GradientType g; g[0] = dirs->at(i)[0]; g[1] = dirs->at(i)[1]; g[2] = dirs->at(i)[2]; gradientList.push_back(g); } } if (m_SelectedBundles.empty()) { if (m_SelectedDWI.IsNotNull()) // add artifacts to existing diffusion weighted image { for (int i=0; i*>(m_SelectedImages.at(i)->GetData())) continue; mitk::DiffusionImage::Pointer diffImg = dynamic_cast*>(m_SelectedImages.at(i)->GetData()); double noiseVariance = 0; if (m_Controls->m_AddNoise->isChecked()) { noiseVariance = m_Controls->m_NoiseLevel->value(); artifactModelString += "_NOISE"; artifactModelString += QString::number(noiseVariance); } mitk::RicianNoiseModel noiseModel; noiseModel.SetNoiseVariance(noiseVariance); if ( this->m_Controls->m_TEbox->value() < imageRegion.GetSize(1)*m_Controls->m_LineReadoutTimeBox->value() ) { this->m_Controls->m_TEbox->setValue( imageRegion.GetSize(1)*m_Controls->m_LineReadoutTimeBox->value() ); QMessageBox::information( NULL, "Warning", "Echo time is too short! Time not sufficient to read slice. Automaticall adjusted to "+QString::number(this->m_Controls->m_TEbox->value())+" ms"); } double lineReadoutTime = m_Controls->m_LineReadoutTimeBox->value(); // add N/2 ghosting double kOffset = 0; if (m_Controls->m_AddGhosts->isChecked()) { artifactModelString += "_GHOST"; kOffset = m_Controls->m_kOffsetBox->value(); } if ( this->m_Controls->m_TEbox->value() < imageRegion.GetSize(1)*m_Controls->m_LineReadoutTimeBox->value() ) { this->m_Controls->m_TEbox->setValue( imageRegion.GetSize(1)*m_Controls->m_LineReadoutTimeBox->value() ); QMessageBox::information( NULL, "Warning", "Echo time is too short! Time not sufficient to read slice. Automaticall adjusted to "+QString::number(this->m_Controls->m_TEbox->value())+" ms"); } itk::AddArtifactsToDwiImageFilter< short >::Pointer filter = itk::AddArtifactsToDwiImageFilter< short >::New(); filter->SetInput(diffImg->GetVectorImage()); filter->SettLine(lineReadoutTime); filter->SetkOffset(kOffset); filter->SetNoiseModel(&noiseModel); filter->SetGradientList(gradientList); filter->SetTE(this->m_Controls->m_TEbox->value()); if (m_Controls->m_AddEddy->isChecked()) { filter->SetSimulateEddyCurrents(true); filter->SetEddyGradientStrength(m_Controls->m_EddyGradientStrength->value()); artifactModelString += "_EDDY"; } - mitk::GibbsRingingArtifact gibbsModel; if (m_Controls->m_AddGibbsRinging->isChecked()) { - resultNode->AddProperty("Fiberfox.k-Space-Undersampling", IntProperty::New(m_Controls->m_KspaceUndersamplingBox->currentText().toInt())); - gibbsModel.SetKspaceCropping((double)m_Controls->m_KspaceUndersamplingBox->currentText().toInt()); - filter->SetRingingModel(&gibbsModel); + resultNode->AddProperty("Fiberfox.kRinging-Upsampling", DoubleProperty::New(m_Controls->m_ImageUpsamplingBox->value())); + filter->SetUpsampling(m_Controls->m_ImageUpsamplingBox->value()); } if (m_Controls->m_AddDistortions->isChecked() && m_Controls->m_FrequencyMapBox->GetSelectedNode().IsNotNull()) { mitk::DataNode::Pointer fMapNode = m_Controls->m_FrequencyMapBox->GetSelectedNode(); mitk::Image* img = dynamic_cast(fMapNode->GetData()); ItkDoubleImgType::Pointer itkImg = ItkDoubleImgType::New(); CastToItkImage< ItkDoubleImgType >(img, itkImg); if (imageRegion.GetSize(0)==itkImg->GetLargestPossibleRegion().GetSize(0) && imageRegion.GetSize(1)==itkImg->GetLargestPossibleRegion().GetSize(1) && imageRegion.GetSize(2)==itkImg->GetLargestPossibleRegion().GetSize(2)) { filter->SetFrequencyMap(itkImg); artifactModelString += "_DISTORTED"; } } filter->Update(); resultNode->AddProperty("Fiberfox.Line-Offset", DoubleProperty::New(kOffset)); resultNode->AddProperty("Fiberfox.Noise-Variance", DoubleProperty::New(noiseVariance)); resultNode->AddProperty("Fiberfox.Tinhom", IntProperty::New(m_Controls->m_T2starBox->value())); resultNode->AddProperty("binary", BoolProperty::New(false)); mitk::DiffusionImage::Pointer image = mitk::DiffusionImage::New(); image->SetVectorImage( filter->GetOutput() ); image->SetB_Value(diffImg->GetB_Value()); image->SetDirections(diffImg->GetDirections()); image->InitializeFromVectorImage(); resultNode->SetData( image ); resultNode->SetName(m_SelectedImages.at(i)->GetName()+artifactModelString.toStdString()); GetDataStorage()->Add(resultNode); } return; } mitk::Image::Pointer image = mitk::ImageGenerator::GenerateGradientImage( m_Controls->m_SizeX->value(), m_Controls->m_SizeY->value(), m_Controls->m_SizeZ->value(), m_Controls->m_SpacingX->value(), m_Controls->m_SpacingY->value(), m_Controls->m_SpacingZ->value()); mitk::Geometry3D* geom = image->GetGeometry(); geom->SetOrigin(origin); mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData( image ); node->SetName("Dummy"); unsigned int window = m_Controls->m_SizeX->value()*m_Controls->m_SizeY->value()*m_Controls->m_SizeZ->value(); unsigned int level = window/2; mitk::LevelWindow lw; lw.SetLevelWindow(level, window); node->SetProperty( "levelwindow", mitk::LevelWindowProperty::New( lw ) ); GetDataStorage()->Add(node); m_SelectedImage = node; mitk::BaseData::Pointer basedata = node->GetData(); if (basedata.IsNotNull()) { mitk::RenderingManager::GetInstance()->InitializeViews( basedata->GetTimeSlicedGeometry(), mitk::RenderingManager::REQUEST_UPDATE_ALL, true ); mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } UpdateGui(); return; } for (int i=0; im_Compartment4Box->currentIndex()>0) { comp4Weight = m_Controls->m_Comp4FractionBox->value(); comp3Weight -= comp4Weight; } mitk::StickModel stickModel1; mitk::StickModel stickModel2; mitk::TensorModel zeppelinModel1; mitk::TensorModel zeppelinModel2; mitk::TensorModel tensorModel1; mitk::TensorModel tensorModel2; mitk::BallModel ballModel1; mitk::BallModel ballModel2; mitk::AstroStickModel astrosticksModel1; mitk::AstroStickModel astrosticksModel2; mitk::DotModel dotModel1; mitk::DotModel dotModel2; // compartment 1 switch (m_Controls->m_Compartment1Box->currentIndex()) { case 0: MITK_INFO << "Using stick model"; stickModel1.SetGradientList(gradientList); stickModel1.SetDiffusivity(m_Controls->m_StickWidget1->GetD()); stickModel1.SetT2(m_Controls->m_StickWidget1->GetT2()); fiberModelList.push_back(&stickModel1); signalModelString += "Stick"; resultNode->AddProperty("Fiberfox.Compartment1.Description", StringProperty::New("Intra-axonal compartment") ); resultNode->AddProperty("Fiberfox.Compartment1.Model", StringProperty::New("Stick") ); resultNode->AddProperty("Fiberfox.Compartment1.D", DoubleProperty::New(m_Controls->m_StickWidget1->GetD()) ); resultNode->AddProperty("Fiberfox.Compartment1.T2", DoubleProperty::New(stickModel1.GetT2()) ); break; case 1: MITK_INFO << "Using zeppelin model"; zeppelinModel1.SetGradientList(gradientList); zeppelinModel1.SetBvalue(bVal); zeppelinModel1.SetDiffusivity1(m_Controls->m_ZeppelinWidget1->GetD1()); zeppelinModel1.SetDiffusivity2(m_Controls->m_ZeppelinWidget1->GetD2()); zeppelinModel1.SetDiffusivity3(m_Controls->m_ZeppelinWidget1->GetD2()); zeppelinModel1.SetT2(m_Controls->m_ZeppelinWidget1->GetT2()); fiberModelList.push_back(&zeppelinModel1); signalModelString += "Zeppelin"; resultNode->AddProperty("Fiberfox.Compartment1.Description", StringProperty::New("Intra-axonal compartment") ); resultNode->AddProperty("Fiberfox.Compartment1.Model", StringProperty::New("Zeppelin") ); resultNode->AddProperty("Fiberfox.Compartment1.D1", DoubleProperty::New(m_Controls->m_ZeppelinWidget1->GetD1()) ); resultNode->AddProperty("Fiberfox.Compartment1.D2", DoubleProperty::New(m_Controls->m_ZeppelinWidget1->GetD2()) ); resultNode->AddProperty("Fiberfox.Compartment1.T2", DoubleProperty::New(zeppelinModel1.GetT2()) ); break; case 2: MITK_INFO << "Using tensor model"; tensorModel1.SetGradientList(gradientList); tensorModel1.SetBvalue(bVal); tensorModel1.SetDiffusivity1(m_Controls->m_TensorWidget1->GetD1()); tensorModel1.SetDiffusivity2(m_Controls->m_TensorWidget1->GetD2()); tensorModel1.SetDiffusivity3(m_Controls->m_TensorWidget1->GetD3()); tensorModel1.SetT2(m_Controls->m_TensorWidget1->GetT2()); fiberModelList.push_back(&tensorModel1); signalModelString += "Tensor"; resultNode->AddProperty("Fiberfox.Compartment1.Description", StringProperty::New("Intra-axonal compartment") ); resultNode->AddProperty("Fiberfox.Compartment1.Model", StringProperty::New("Tensor") ); resultNode->AddProperty("Fiberfox.Compartment1.D1", DoubleProperty::New(m_Controls->m_TensorWidget1->GetD1()) ); resultNode->AddProperty("Fiberfox.Compartment1.D2", DoubleProperty::New(m_Controls->m_TensorWidget1->GetD2()) ); resultNode->AddProperty("Fiberfox.Compartment1.D3", DoubleProperty::New(m_Controls->m_TensorWidget1->GetD3()) ); resultNode->AddProperty("Fiberfox.Compartment1.T2", DoubleProperty::New(zeppelinModel1.GetT2()) ); break; } // compartment 2 switch (m_Controls->m_Compartment2Box->currentIndex()) { case 0: break; case 1: stickModel2.SetGradientList(gradientList); stickModel2.SetDiffusivity(m_Controls->m_StickWidget2->GetD()); stickModel2.SetT2(m_Controls->m_StickWidget2->GetT2()); fiberModelList.push_back(&stickModel2); signalModelString += "Stick"; resultNode->AddProperty("Fiberfox.Compartment2.Description", StringProperty::New("Inter-axonal compartment") ); resultNode->AddProperty("Fiberfox.Compartment2.Model", StringProperty::New("Stick") ); resultNode->AddProperty("Fiberfox.Compartment2.D", DoubleProperty::New(m_Controls->m_StickWidget2->GetD()) ); resultNode->AddProperty("Fiberfox.Compartment2.T2", DoubleProperty::New(stickModel2.GetT2()) ); break; case 2: zeppelinModel2.SetGradientList(gradientList); zeppelinModel2.SetBvalue(bVal); zeppelinModel2.SetDiffusivity1(m_Controls->m_ZeppelinWidget2->GetD1()); zeppelinModel2.SetDiffusivity2(m_Controls->m_ZeppelinWidget2->GetD2()); zeppelinModel2.SetDiffusivity3(m_Controls->m_ZeppelinWidget2->GetD2()); zeppelinModel2.SetT2(m_Controls->m_ZeppelinWidget2->GetT2()); fiberModelList.push_back(&zeppelinModel2); signalModelString += "Zeppelin"; resultNode->AddProperty("Fiberfox.Compartment2.Description", StringProperty::New("Inter-axonal compartment") ); resultNode->AddProperty("Fiberfox.Compartment2.Model", StringProperty::New("Zeppelin") ); resultNode->AddProperty("Fiberfox.Compartment2.D1", DoubleProperty::New(m_Controls->m_ZeppelinWidget2->GetD1()) ); resultNode->AddProperty("Fiberfox.Compartment2.D2", DoubleProperty::New(m_Controls->m_ZeppelinWidget2->GetD2()) ); resultNode->AddProperty("Fiberfox.Compartment2.T2", DoubleProperty::New(zeppelinModel2.GetT2()) ); break; case 3: tensorModel2.SetGradientList(gradientList); tensorModel2.SetBvalue(bVal); tensorModel2.SetDiffusivity1(m_Controls->m_TensorWidget2->GetD1()); tensorModel2.SetDiffusivity2(m_Controls->m_TensorWidget2->GetD2()); tensorModel2.SetDiffusivity3(m_Controls->m_TensorWidget2->GetD3()); tensorModel2.SetT2(m_Controls->m_TensorWidget2->GetT2()); fiberModelList.push_back(&tensorModel2); signalModelString += "Tensor"; resultNode->AddProperty("Fiberfox.Compartment2.Description", StringProperty::New("Inter-axonal compartment") ); resultNode->AddProperty("Fiberfox.Compartment2.Model", StringProperty::New("Tensor") ); resultNode->AddProperty("Fiberfox.Compartment2.D1", DoubleProperty::New(m_Controls->m_TensorWidget2->GetD1()) ); resultNode->AddProperty("Fiberfox.Compartment2.D2", DoubleProperty::New(m_Controls->m_TensorWidget2->GetD2()) ); resultNode->AddProperty("Fiberfox.Compartment2.D3", DoubleProperty::New(m_Controls->m_TensorWidget2->GetD3()) ); resultNode->AddProperty("Fiberfox.Compartment2.T2", DoubleProperty::New(zeppelinModel2.GetT2()) ); break; } // compartment 3 switch (m_Controls->m_Compartment3Box->currentIndex()) { case 0: ballModel1.SetGradientList(gradientList); ballModel1.SetBvalue(bVal); ballModel1.SetDiffusivity(m_Controls->m_BallWidget1->GetD()); ballModel1.SetT2(m_Controls->m_BallWidget1->GetT2()); ballModel1.SetWeight(comp3Weight); nonFiberModelList.push_back(&ballModel1); signalModelString += "Ball"; resultNode->AddProperty("Fiberfox.Compartment3.Description", StringProperty::New("Extra-axonal compartment 1") ); resultNode->AddProperty("Fiberfox.Compartment3.Model", StringProperty::New("Ball") ); resultNode->AddProperty("Fiberfox.Compartment3.D", DoubleProperty::New(m_Controls->m_BallWidget1->GetD()) ); resultNode->AddProperty("Fiberfox.Compartment3.T2", DoubleProperty::New(ballModel1.GetT2()) ); break; case 1: astrosticksModel1.SetGradientList(gradientList); astrosticksModel1.SetBvalue(bVal); astrosticksModel1.SetDiffusivity(m_Controls->m_AstrosticksWidget1->GetD()); astrosticksModel1.SetT2(m_Controls->m_AstrosticksWidget1->GetT2()); astrosticksModel1.SetRandomizeSticks(m_Controls->m_AstrosticksWidget1->GetRandomizeSticks()); astrosticksModel1.SetWeight(comp3Weight); nonFiberModelList.push_back(&astrosticksModel1); signalModelString += "Astrosticks"; resultNode->AddProperty("Fiberfox.Compartment3.Description", StringProperty::New("Extra-axonal compartment 1") ); resultNode->AddProperty("Fiberfox.Compartment3.Model", StringProperty::New("Astrosticks") ); resultNode->AddProperty("Fiberfox.Compartment3.D", DoubleProperty::New(m_Controls->m_AstrosticksWidget1->GetD()) ); resultNode->AddProperty("Fiberfox.Compartment3.T2", DoubleProperty::New(astrosticksModel1.GetT2()) ); resultNode->AddProperty("Fiberfox.Compartment3.RandomSticks", BoolProperty::New(m_Controls->m_AstrosticksWidget1->GetRandomizeSticks()) ); break; case 2: dotModel1.SetGradientList(gradientList); dotModel1.SetT2(m_Controls->m_DotWidget1->GetT2()); dotModel1.SetWeight(comp3Weight); nonFiberModelList.push_back(&dotModel1); signalModelString += "Dot"; resultNode->AddProperty("Fiberfox.Compartment3.Description", StringProperty::New("Extra-axonal compartment 1") ); resultNode->AddProperty("Fiberfox.Compartment3.Model", StringProperty::New("Dot") ); resultNode->AddProperty("Fiberfox.Compartment3.T2", DoubleProperty::New(dotModel1.GetT2()) ); break; } // compartment 4 switch (m_Controls->m_Compartment4Box->currentIndex()) { case 0: break; case 1: ballModel2.SetGradientList(gradientList); ballModel2.SetBvalue(bVal); ballModel2.SetDiffusivity(m_Controls->m_BallWidget2->GetD()); ballModel2.SetT2(m_Controls->m_BallWidget2->GetT2()); ballModel2.SetWeight(comp4Weight); nonFiberModelList.push_back(&ballModel2); signalModelString += "Ball"; resultNode->AddProperty("Fiberfox.Compartment4.Description", StringProperty::New("Extra-axonal compartment 2") ); resultNode->AddProperty("Fiberfox.Compartment4.Model", StringProperty::New("Ball") ); resultNode->AddProperty("Fiberfox.Compartment4.D", DoubleProperty::New(m_Controls->m_BallWidget2->GetD()) ); resultNode->AddProperty("Fiberfox.Compartment4.T2", DoubleProperty::New(ballModel2.GetT2()) ); break; case 2: astrosticksModel2.SetGradientList(gradientList); astrosticksModel2.SetBvalue(bVal); astrosticksModel2.SetDiffusivity(m_Controls->m_AstrosticksWidget2->GetD()); astrosticksModel2.SetT2(m_Controls->m_AstrosticksWidget2->GetT2()); astrosticksModel2.SetRandomizeSticks(m_Controls->m_AstrosticksWidget2->GetRandomizeSticks()); astrosticksModel2.SetWeight(comp4Weight); nonFiberModelList.push_back(&astrosticksModel2); signalModelString += "Astrosticks"; resultNode->AddProperty("Fiberfox.Compartment4.Description", StringProperty::New("Extra-axonal compartment 2") ); resultNode->AddProperty("Fiberfox.Compartment4.Model", StringProperty::New("Astrosticks") ); resultNode->AddProperty("Fiberfox.Compartment4.D", DoubleProperty::New(m_Controls->m_AstrosticksWidget2->GetD()) ); resultNode->AddProperty("Fiberfox.Compartment4.T2", DoubleProperty::New(astrosticksModel2.GetT2()) ); resultNode->AddProperty("Fiberfox.Compartment4.RandomSticks", BoolProperty::New(m_Controls->m_AstrosticksWidget2->GetRandomizeSticks()) ); break; case 3: dotModel2.SetGradientList(gradientList); dotModel2.SetT2(m_Controls->m_DotWidget2->GetT2()); dotModel2.SetWeight(comp4Weight); nonFiberModelList.push_back(&dotModel2); signalModelString += "Dot"; resultNode->AddProperty("Fiberfox.Compartment4.Description", StringProperty::New("Extra-axonal compartment 2") ); resultNode->AddProperty("Fiberfox.Compartment4.Model", StringProperty::New("Dot") ); resultNode->AddProperty("Fiberfox.Compartment4.T2", DoubleProperty::New(dotModel2.GetT2()) ); break; } itk::TractsToDWIImageFilter::KspaceArtifactList artifactList; // artifact models QString artifactModelString(""); double noiseVariance = 0; if (m_Controls->m_AddNoise->isChecked()) { noiseVariance = m_Controls->m_NoiseLevel->value(); artifactModelString += "_NOISE"; artifactModelString += QString::number(noiseVariance); resultNode->AddProperty("Fiberfox.Noise-Variance", DoubleProperty::New(noiseVariance)); } mitk::RicianNoiseModel noiseModel; noiseModel.SetNoiseVariance(noiseVariance); - mitk::GibbsRingingArtifact gibbsModel; if (m_Controls->m_AddGibbsRinging->isChecked()) { artifactModelString += "_RINGING"; - resultNode->AddProperty("Fiberfox.k-Space-Undersampling", IntProperty::New(m_Controls->m_KspaceUndersamplingBox->currentText().toInt())); - gibbsModel.SetKspaceCropping((double)m_Controls->m_KspaceUndersamplingBox->currentText().toInt()); - artifactList.push_back(&gibbsModel); + resultNode->AddProperty("Fiberfox.Ringing-Upsampling", DoubleProperty::New(m_Controls->m_ImageUpsamplingBox->value())); + tractsToDwiFilter->SetUpsampling(m_Controls->m_ImageUpsamplingBox->value()); } if ( this->m_Controls->m_TEbox->value() < imageRegion.GetSize(1)*m_Controls->m_LineReadoutTimeBox->value() ) { this->m_Controls->m_TEbox->setValue( imageRegion.GetSize(1)*m_Controls->m_LineReadoutTimeBox->value() ); QMessageBox::information( NULL, "Warning", "Echo time is too short! Time not sufficient to read slice. Automaticall adjusted to "+QString::number(this->m_Controls->m_TEbox->value())+" ms"); } double lineReadoutTime = m_Controls->m_LineReadoutTimeBox->value(); // adjusting line readout time to the adapted image size needed for the DFT int y = imageRegion.GetSize(1); if ( y%2 == 1 ) y += 1; if ( y>imageRegion.GetSize(1) ) lineReadoutTime *= (double)imageRegion.GetSize(1)/y; // add N/2 ghosting double kOffset = 0; if (m_Controls->m_AddGhosts->isChecked()) { artifactModelString += "_GHOST"; kOffset = m_Controls->m_kOffsetBox->value(); resultNode->AddProperty("Fiberfox.Line-Offset", DoubleProperty::New(kOffset)); } // add distortions if (m_Controls->m_AddDistortions->isChecked() && m_Controls->m_FrequencyMapBox->GetSelectedNode().IsNotNull()) { mitk::DataNode::Pointer fMapNode = m_Controls->m_FrequencyMapBox->GetSelectedNode(); mitk::Image* img = dynamic_cast(fMapNode->GetData()); ItkDoubleImgType::Pointer itkImg = ItkDoubleImgType::New(); CastToItkImage< ItkDoubleImgType >(img, itkImg); if (imageRegion.GetSize(0)==itkImg->GetLargestPossibleRegion().GetSize(0) && imageRegion.GetSize(1)==itkImg->GetLargestPossibleRegion().GetSize(1) && imageRegion.GetSize(2)==itkImg->GetLargestPossibleRegion().GetSize(2)) { tractsToDwiFilter->SetFrequencyMap(itkImg); artifactModelString += "_DISTORTED"; } } if (m_Controls->m_AddEddy->isChecked()) { tractsToDwiFilter->SetSimulateEddyCurrents(true); tractsToDwiFilter->SetEddyGradientStrength(m_Controls->m_EddyGradientStrength->value()); artifactModelString += "_EDDY"; } mitk::FiberBundleX::Pointer fiberBundle = dynamic_cast(m_SelectedBundles.at(i)->GetData()); if (fiberBundle->GetNumFibers()<=0) continue; if (m_Controls->m_RelaxationBox->isChecked()) artifactModelString += "_RELAX"; tractsToDwiFilter->SetSimulateRelaxation(m_Controls->m_RelaxationBox->isChecked()); tractsToDwiFilter->SetImageRegion(imageRegion); tractsToDwiFilter->SetSpacing(spacing); tractsToDwiFilter->SetOrigin(origin); tractsToDwiFilter->SetDirectionMatrix(directionMatrix); tractsToDwiFilter->SetFiberBundle(fiberBundle); tractsToDwiFilter->SetFiberModels(fiberModelList); tractsToDwiFilter->SetNonFiberModels(nonFiberModelList); tractsToDwiFilter->SetNoiseModel(&noiseModel); tractsToDwiFilter->SetKspaceArtifacts(artifactList); tractsToDwiFilter->SetkOffset(kOffset); tractsToDwiFilter->SettLine(m_Controls->m_LineReadoutTimeBox->value()); tractsToDwiFilter->SettInhom(this->m_Controls->m_T2starBox->value()); tractsToDwiFilter->SetTE(this->m_Controls->m_TEbox->value()); tractsToDwiFilter->SetNumberOfRepetitions(m_Controls->m_RepetitionsBox->value()); tractsToDwiFilter->SetEnforcePureFiberVoxels(m_Controls->m_EnforcePureFiberVoxelsBox->isChecked()); tractsToDwiFilter->SetInterpolationShrink(m_Controls->m_InterpolationShrink->value()); tractsToDwiFilter->SetFiberRadius(m_Controls->m_FiberRadius->value()); tractsToDwiFilter->SetSignalScale(m_Controls->m_SignalScaleBox->value()); if (m_Controls->m_InterpolationShrink->value()<1000) tractsToDwiFilter->SetUseInterpolation(true); if (m_TissueMask.IsNotNull()) { ItkUcharImgType::Pointer mask = ItkUcharImgType::New(); mitk::CastToItkImage(m_TissueMask, mask); tractsToDwiFilter->SetTissueMask(mask); } tractsToDwiFilter->Update(); mitk::DiffusionImage::Pointer image = mitk::DiffusionImage::New(); image->SetVectorImage( tractsToDwiFilter->GetOutput() ); image->SetB_Value(bVal); image->SetDirections(gradientList); image->InitializeFromVectorImage(); resultNode->SetData( image ); resultNode->SetName(m_SelectedBundles.at(i)->GetName() +"_D"+QString::number(imageRegion.GetSize(0)).toStdString() +"-"+QString::number(imageRegion.GetSize(1)).toStdString() +"-"+QString::number(imageRegion.GetSize(2)).toStdString() +"_S"+QString::number(spacing[0]).toStdString() +"-"+QString::number(spacing[1]).toStdString() +"-"+QString::number(spacing[2]).toStdString() +"_b"+QString::number(bVal).toStdString() +"_"+signalModelString.toStdString() +artifactModelString.toStdString()); GetDataStorage()->Add(resultNode, m_SelectedBundles.at(i)); resultNode->AddProperty("Fiberfox.InterpolationShrink", IntProperty::New(m_Controls->m_InterpolationShrink->value())); resultNode->AddProperty("Fiberfox.SignalScale", IntProperty::New(m_Controls->m_SignalScaleBox->value())); resultNode->AddProperty("Fiberfox.FiberRadius", IntProperty::New(m_Controls->m_FiberRadius->value())); resultNode->AddProperty("Fiberfox.Tinhom", IntProperty::New(m_Controls->m_T2starBox->value())); resultNode->AddProperty("Fiberfox.Repetitions", IntProperty::New(m_Controls->m_RepetitionsBox->value())); resultNode->AddProperty("Fiberfox.b-value", DoubleProperty::New(bVal)); resultNode->AddProperty("Fiberfox.Model", StringProperty::New(signalModelString.toStdString())); resultNode->AddProperty("Fiberfox.PureFiberVoxels", BoolProperty::New(m_Controls->m_EnforcePureFiberVoxelsBox->isChecked())); resultNode->AddProperty("binary", BoolProperty::New(false)); resultNode->SetProperty( "levelwindow", mitk::LevelWindowProperty::New(tractsToDwiFilter->GetLevelWindow()) ); if (m_Controls->m_VolumeFractionsBox->isChecked()) { std::vector< itk::TractsToDWIImageFilter::ItkDoubleImgType::Pointer > volumeFractions = tractsToDwiFilter->GetVolumeFractions(); for (int k=0; kInitializeByItk(volumeFractions.at(k).GetPointer()); image->SetVolume(volumeFractions.at(k)->GetBufferPointer()); mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData( image ); node->SetName(m_SelectedBundles.at(i)->GetName()+"_CompartmentVolume-"+QString::number(k).toStdString()); GetDataStorage()->Add(node, m_SelectedBundles.at(i)); } } mitk::BaseData::Pointer basedata = resultNode->GetData(); if (basedata.IsNotNull()) { mitk::RenderingManager::GetInstance()->InitializeViews( basedata->GetTimeSlicedGeometry(), mitk::RenderingManager::REQUEST_UPDATE_ALL, true ); mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } } } void QmitkFiberfoxView::ApplyTransform() { vector< mitk::DataNode::Pointer > selectedBundles; for( int i=0; iGetDerivations(m_SelectedImages.at(i)); for( mitk::DataStorage::SetOfObjects::const_iterator it = derivations->begin(); it != derivations->end(); ++it ) { mitk::DataNode::Pointer fibNode = *it; if ( fibNode.IsNotNull() && dynamic_cast(fibNode->GetData()) ) selectedBundles.push_back(fibNode); } } if (selectedBundles.empty()) selectedBundles = m_SelectedBundles2; if (!selectedBundles.empty()) { std::vector::const_iterator it = selectedBundles.begin(); for (it; it!=selectedBundles.end(); ++it) { mitk::FiberBundleX::Pointer fib = dynamic_cast((*it)->GetData()); fib->RotateAroundAxis(m_Controls->m_XrotBox->value(), m_Controls->m_YrotBox->value(), m_Controls->m_ZrotBox->value()); fib->TranslateFibers(m_Controls->m_XtransBox->value(), m_Controls->m_YtransBox->value(), m_Controls->m_ZtransBox->value()); fib->ScaleFibers(m_Controls->m_XscaleBox->value(), m_Controls->m_YscaleBox->value(), m_Controls->m_ZscaleBox->value()); // handle child fiducials if (m_Controls->m_IncludeFiducials->isChecked()) { mitk::DataStorage::SetOfObjects::ConstPointer derivations = GetDataStorage()->GetDerivations(*it); for( mitk::DataStorage::SetOfObjects::const_iterator it2 = derivations->begin(); it2 != derivations->end(); ++it2 ) { mitk::DataNode::Pointer fiducialNode = *it2; if ( fiducialNode.IsNotNull() && dynamic_cast(fiducialNode->GetData()) ) { mitk::PlanarEllipse* pe = dynamic_cast(fiducialNode->GetData()); mitk::Geometry3D* geom = pe->GetGeometry(); // translate mitk::Vector3D world; world[0] = m_Controls->m_XtransBox->value(); world[1] = m_Controls->m_YtransBox->value(); world[2] = m_Controls->m_ZtransBox->value(); geom->Translate(world); // calculate rotation matrix double x = m_Controls->m_XrotBox->value()*M_PI/180; double y = m_Controls->m_YrotBox->value()*M_PI/180; double z = m_Controls->m_ZrotBox->value()*M_PI/180; itk::Matrix< float, 3, 3 > rotX; rotX.SetIdentity(); rotX[1][1] = cos(x); rotX[2][2] = rotX[1][1]; rotX[1][2] = -sin(x); rotX[2][1] = -rotX[1][2]; itk::Matrix< float, 3, 3 > rotY; rotY.SetIdentity(); rotY[0][0] = cos(y); rotY[2][2] = rotY[0][0]; rotY[0][2] = sin(y); rotY[2][0] = -rotY[0][2]; itk::Matrix< float, 3, 3 > rotZ; rotZ.SetIdentity(); rotZ[0][0] = cos(z); rotZ[1][1] = rotZ[0][0]; rotZ[0][1] = -sin(z); rotZ[1][0] = -rotZ[0][1]; itk::Matrix< float, 3, 3 > rot = rotZ*rotY*rotX; // transform control point coordinate into geometry translation geom->SetOrigin(pe->GetWorldControlPoint(0)); mitk::Point2D cp; cp.Fill(0.0); pe->SetControlPoint(0, cp); // rotate fiducial geom->GetIndexToWorldTransform()->SetMatrix(rot*geom->GetIndexToWorldTransform()->GetMatrix()); // implicit translation mitk::Vector3D trans; trans[0] = geom->GetOrigin()[0]-fib->GetGeometry()->GetCenter()[0]; trans[1] = geom->GetOrigin()[1]-fib->GetGeometry()->GetCenter()[1]; trans[2] = geom->GetOrigin()[2]-fib->GetGeometry()->GetCenter()[2]; mitk::Vector3D newWc = rot*trans; newWc = newWc-trans; geom->Translate(newWc); } } } } } else { for (int i=0; i(m_SelectedFiducials.at(i)->GetData()); mitk::Geometry3D* geom = pe->GetGeometry(); // translate mitk::Vector3D world; world[0] = m_Controls->m_XtransBox->value(); world[1] = m_Controls->m_YtransBox->value(); world[2] = m_Controls->m_ZtransBox->value(); geom->Translate(world); // calculate rotation matrix double x = m_Controls->m_XrotBox->value()*M_PI/180; double y = m_Controls->m_YrotBox->value()*M_PI/180; double z = m_Controls->m_ZrotBox->value()*M_PI/180; itk::Matrix< float, 3, 3 > rotX; rotX.SetIdentity(); rotX[1][1] = cos(x); rotX[2][2] = rotX[1][1]; rotX[1][2] = -sin(x); rotX[2][1] = -rotX[1][2]; itk::Matrix< float, 3, 3 > rotY; rotY.SetIdentity(); rotY[0][0] = cos(y); rotY[2][2] = rotY[0][0]; rotY[0][2] = sin(y); rotY[2][0] = -rotY[0][2]; itk::Matrix< float, 3, 3 > rotZ; rotZ.SetIdentity(); rotZ[0][0] = cos(z); rotZ[1][1] = rotZ[0][0]; rotZ[0][1] = -sin(z); rotZ[1][0] = -rotZ[0][1]; itk::Matrix< float, 3, 3 > rot = rotZ*rotY*rotX; // transform control point coordinate into geometry translation geom->SetOrigin(pe->GetWorldControlPoint(0)); mitk::Point2D cp; cp.Fill(0.0); pe->SetControlPoint(0, cp); // rotate fiducial geom->GetIndexToWorldTransform()->SetMatrix(rot*geom->GetIndexToWorldTransform()->GetMatrix()); } if (m_Controls->m_RealTimeFibers->isChecked()) GenerateFibers(); } mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberfoxView::CopyBundles() { if ( m_SelectedBundles.size()<1 ){ QMessageBox::information( NULL, "Warning", "Select at least one fiber bundle!"); MITK_WARN("QmitkFiberProcessingView") << "Select at least one fiber bundle!"; return; } std::vector::const_iterator it = m_SelectedBundles.begin(); for (it; it!=m_SelectedBundles.end(); ++it) { // find parent image mitk::DataNode::Pointer parentNode; mitk::DataStorage::SetOfObjects::ConstPointer parentImgs = GetDataStorage()->GetSources(*it); for( mitk::DataStorage::SetOfObjects::const_iterator it2 = parentImgs->begin(); it2 != parentImgs->end(); ++it2 ) { mitk::DataNode::Pointer pImgNode = *it2; if ( pImgNode.IsNotNull() && dynamic_cast(pImgNode->GetData()) ) { parentNode = pImgNode; break; } } mitk::FiberBundleX::Pointer fib = dynamic_cast((*it)->GetData()); mitk::FiberBundleX::Pointer newBundle = fib->GetDeepCopy(); QString name((*it)->GetName().c_str()); name += "_copy"; mitk::DataNode::Pointer fbNode = mitk::DataNode::New(); fbNode->SetData(newBundle); fbNode->SetName(name.toStdString()); fbNode->SetVisibility(true); if (parentNode.IsNotNull()) GetDataStorage()->Add(fbNode, parentNode); else GetDataStorage()->Add(fbNode); // copy child fiducials if (m_Controls->m_IncludeFiducials->isChecked()) { mitk::DataStorage::SetOfObjects::ConstPointer derivations = GetDataStorage()->GetDerivations(*it); for( mitk::DataStorage::SetOfObjects::const_iterator it2 = derivations->begin(); it2 != derivations->end(); ++it2 ) { mitk::DataNode::Pointer fiducialNode = *it2; if ( fiducialNode.IsNotNull() && dynamic_cast(fiducialNode->GetData()) ) { mitk::PlanarEllipse::Pointer pe = mitk::PlanarEllipse::New(); pe->DeepCopy(dynamic_cast(fiducialNode->GetData())); mitk::DataNode::Pointer newNode = mitk::DataNode::New(); newNode->SetData(pe); newNode->SetName(fiducialNode->GetName()); GetDataStorage()->Add(newNode, fbNode); } } } } mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberfoxView::JoinBundles() { if ( m_SelectedBundles.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_SelectedBundles.begin(); mitk::FiberBundleX::Pointer newBundle = dynamic_cast((*it)->GetData()); QString name(""); name += QString((*it)->GetName().c_str()); ++it; for (it; it!=m_SelectedBundles.end(); ++it) { newBundle = newBundle->AddBundle(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); mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberfoxView::UpdateGui() { m_Controls->m_FiberBundleLabel->setText("mandatory"); m_Controls->m_GeometryFrame->setEnabled(true); m_Controls->m_GeometryMessage->setVisible(false); m_Controls->m_DiffusionPropsMessage->setVisible(false); m_Controls->m_FiberGenMessage->setVisible(true); m_Controls->m_TransformBundlesButton->setEnabled(false); m_Controls->m_CopyBundlesButton->setEnabled(false); m_Controls->m_GenerateFibersButton->setEnabled(false); m_Controls->m_FlipButton->setEnabled(false); m_Controls->m_CircleButton->setEnabled(false); m_Controls->m_BvalueBox->setEnabled(true); m_Controls->m_NumGradientsBox->setEnabled(true); m_Controls->m_JoinBundlesButton->setEnabled(false); m_Controls->m_AlignOnGrid->setEnabled(false); if (m_SelectedFiducial.IsNotNull()) { m_Controls->m_TransformBundlesButton->setEnabled(true); m_Controls->m_FlipButton->setEnabled(true); m_Controls->m_AlignOnGrid->setEnabled(true); } if (m_SelectedImage.IsNotNull() || !m_SelectedBundles.empty()) { m_Controls->m_TransformBundlesButton->setEnabled(true); m_Controls->m_CircleButton->setEnabled(true); m_Controls->m_FiberGenMessage->setVisible(false); m_Controls->m_AlignOnGrid->setEnabled(true); } if (m_TissueMask.IsNotNull() || m_SelectedImage.IsNotNull()) { m_Controls->m_GeometryMessage->setVisible(true); m_Controls->m_GeometryFrame->setEnabled(false); } if (m_SelectedDWI.IsNotNull()) { m_Controls->m_DiffusionPropsMessage->setVisible(true); m_Controls->m_BvalueBox->setEnabled(false); m_Controls->m_NumGradientsBox->setEnabled(false); m_Controls->m_GeometryMessage->setVisible(true); m_Controls->m_GeometryFrame->setEnabled(false); } if (!m_SelectedBundles.empty()) { m_Controls->m_CopyBundlesButton->setEnabled(true); m_Controls->m_GenerateFibersButton->setEnabled(true); m_Controls->m_FiberBundleLabel->setText(m_SelectedBundles.at(0)->GetName().c_str()); if (m_SelectedBundles.size()>1) m_Controls->m_JoinBundlesButton->setEnabled(true); } } void QmitkFiberfoxView::OnSelectionChanged( berry::IWorkbenchPart::Pointer, const QList& nodes ) { m_SelectedBundles2.clear(); m_SelectedImages.clear(); m_SelectedFiducials.clear(); m_SelectedFiducial = NULL; m_TissueMask = NULL; m_SelectedBundles.clear(); m_SelectedImage = NULL; m_SelectedDWI = NULL; m_Controls->m_TissueMaskLabel->setText("optional"); // iterate all selected objects, adjust warning visibility for( int i=0; i*>(node->GetData()) ) { m_SelectedDWI = node; m_SelectedImage = node; m_SelectedImages.push_back(node); } else if( node.IsNotNull() && dynamic_cast(node->GetData()) ) { m_SelectedImages.push_back(node); m_SelectedImage = node; bool isBinary = false; node->GetPropertyValue("binary", isBinary); if (isBinary) { m_TissueMask = dynamic_cast(node->GetData()); m_Controls->m_TissueMaskLabel->setText(node->GetName().c_str()); } } else if ( node.IsNotNull() && dynamic_cast(node->GetData()) ) { m_SelectedBundles2.push_back(node); if (m_Controls->m_RealTimeFibers->isChecked()) { m_SelectedBundles.push_back(node); mitk::FiberBundleX::Pointer newFib = dynamic_cast(node->GetData()); if (newFib->GetNumFibers()!=m_Controls->m_FiberDensityBox->value()) GenerateFibers(); } else m_SelectedBundles.push_back(node); } else if ( node.IsNotNull() && dynamic_cast(node->GetData()) ) { m_SelectedFiducials.push_back(node); m_SelectedFiducial = node; m_SelectedBundles.clear(); mitk::DataStorage::SetOfObjects::ConstPointer parents = GetDataStorage()->GetSources(node); for( mitk::DataStorage::SetOfObjects::const_iterator it = parents->begin(); it != parents->end(); ++it ) { mitk::DataNode::Pointer pNode = *it; if ( pNode.IsNotNull() && dynamic_cast(pNode->GetData()) ) m_SelectedBundles.push_back(pNode); } } } UpdateGui(); } void QmitkFiberfoxView::EnableCrosshairNavigation() { MITK_DEBUG << "EnableCrosshairNavigation"; // enable the crosshair navigation if (mitk::ILinkedRenderWindowPart* linkedRenderWindow = dynamic_cast(this->GetRenderWindowPart())) { MITK_DEBUG << "enabling linked navigation"; linkedRenderWindow->EnableLinkedNavigation(true); // linkedRenderWindow->EnableSlicingPlanes(true); } if (m_Controls->m_RealTimeFibers->isChecked()) GenerateFibers(); } void QmitkFiberfoxView::DisableCrosshairNavigation() { MITK_DEBUG << "DisableCrosshairNavigation"; // disable the crosshair navigation during the drawing if (mitk::ILinkedRenderWindowPart* linkedRenderWindow = dynamic_cast(this->GetRenderWindowPart())) { MITK_DEBUG << "disabling linked navigation"; linkedRenderWindow->EnableLinkedNavigation(false); // linkedRenderWindow->EnableSlicingPlanes(false); } } void QmitkFiberfoxView::NodeRemoved(const mitk::DataNode* node) { mitk::DataNode* nonConstNode = const_cast(node); std::map::iterator it = m_DataNodeToPlanarFigureData.find(nonConstNode); if( it != m_DataNodeToPlanarFigureData.end() ) { QmitkPlanarFigureData& data = it->second; // remove observers data.m_Figure->RemoveObserver( data.m_EndPlacementObserverTag ); data.m_Figure->RemoveObserver( data.m_SelectObserverTag ); data.m_Figure->RemoveObserver( data.m_StartInteractionObserverTag ); data.m_Figure->RemoveObserver( data.m_EndInteractionObserverTag ); m_DataNodeToPlanarFigureData.erase( it ); } } void QmitkFiberfoxView::NodeAdded( const mitk::DataNode* node ) { // add observer for selection in renderwindow mitk::PlanarFigure* figure = dynamic_cast(node->GetData()); bool isPositionMarker (false); node->GetBoolProperty("isContourMarker", isPositionMarker); if( figure && !isPositionMarker ) { MITK_DEBUG << "figure added. will add interactor if needed."; mitk::PlanarFigureInteractor::Pointer figureInteractor = dynamic_cast(node->GetInteractor()); mitk::DataNode* nonConstNode = const_cast( node ); if(figureInteractor.IsNull()) { figureInteractor = mitk::PlanarFigureInteractor::New("PlanarFigureInteractor", nonConstNode); } else { // just to be sure that the interactor is not added twice mitk::GlobalInteraction::GetInstance()->RemoveInteractor(figureInteractor); } MITK_DEBUG << "adding interactor to globalinteraction"; mitk::GlobalInteraction::GetInstance()->AddInteractor(figureInteractor); MITK_DEBUG << "will now add observers for planarfigure"; QmitkPlanarFigureData data; data.m_Figure = figure; // // add observer for event when figure has been placed typedef itk::SimpleMemberCommand< QmitkFiberfoxView > SimpleCommandType; // SimpleCommandType::Pointer initializationCommand = SimpleCommandType::New(); // initializationCommand->SetCallbackFunction( this, &QmitkFiberfoxView::PlanarFigureInitialized ); // data.m_EndPlacementObserverTag = figure->AddObserver( mitk::EndPlacementPlanarFigureEvent(), initializationCommand ); // add observer for event when figure is picked (selected) typedef itk::MemberCommand< QmitkFiberfoxView > MemberCommandType; MemberCommandType::Pointer selectCommand = MemberCommandType::New(); selectCommand->SetCallbackFunction( this, &QmitkFiberfoxView::PlanarFigureSelected ); data.m_SelectObserverTag = figure->AddObserver( mitk::SelectPlanarFigureEvent(), selectCommand ); // add observer for event when interaction with figure starts SimpleCommandType::Pointer startInteractionCommand = SimpleCommandType::New(); startInteractionCommand->SetCallbackFunction( this, &QmitkFiberfoxView::DisableCrosshairNavigation); data.m_StartInteractionObserverTag = figure->AddObserver( mitk::StartInteractionPlanarFigureEvent(), startInteractionCommand ); // add observer for event when interaction with figure starts SimpleCommandType::Pointer endInteractionCommand = SimpleCommandType::New(); endInteractionCommand->SetCallbackFunction( this, &QmitkFiberfoxView::EnableCrosshairNavigation); data.m_EndInteractionObserverTag = figure->AddObserver( mitk::EndInteractionPlanarFigureEvent(), endInteractionCommand ); m_DataNodeToPlanarFigureData[nonConstNode] = data; } } void QmitkFiberfoxView::PlanarFigureSelected( itk::Object* object, const itk::EventObject& ) { mitk::TNodePredicateDataType::Pointer isPf = mitk::TNodePredicateDataType::New(); mitk::DataStorage::SetOfObjects::ConstPointer allPfs = this->GetDataStorage()->GetSubset( isPf ); for ( mitk::DataStorage::SetOfObjects::const_iterator it = allPfs->begin(); it!=allPfs->end(); ++it) { mitk::DataNode* node = *it; if( node->GetData() == object ) { node->SetSelected(true); m_SelectedFiducial = node; } else node->SetSelected(false); } UpdateGui(); this->RequestRenderWindowUpdate(); } void QmitkFiberfoxView::SetFocus() { m_Controls->m_CircleButton->setFocus(); } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberfoxView.h b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberfoxView.h old mode 100644 new mode 100755 diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberfoxViewControls.ui b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberfoxViewControls.ui old mode 100644 new mode 100755 index 9b378012bd..63972bd20f --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberfoxViewControls.ui +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberfoxViewControls.ui @@ -1,2227 +1,2204 @@ QmitkFiberfoxViewControls 0 0 493 1565 Form 0 Fiber Definition Qt::Vertical 20 40 color: rgb(255, 0, 0); Please select an image or an existing fiber bundle to draw the fiber fiducials. If you can't provide a suitable image, generate one using the "Signal Generation" tab. Qt::AutoText Qt::AlignJustify|Qt::AlignVCenter true Fiducial Options All fiducials are treated as circles with the same radius as the first fiducial. Use Constant Fiducial Radius false false Align selected fiducials with voxel grid. Shifts selected fiducials to nearest voxel center. Align With Grid Operations false Copy Bundles false Transform Selection QFrame::NoFrame QFrame::Raised 0 Y false Rotation angle (in degree) around x-axis. -360.000000000000000 360.000000000000000 0.100000000000000 Axis: false Rotation angle (in degree) around y-axis. -360.000000000000000 360.000000000000000 0.100000000000000 Translation: false Translation (in mm) in direction of the z-axis. -1000.000000000000000 1000.000000000000000 0.100000000000000 Translation (in mm) in direction of the y-axis. -1000.000000000000000 1000.000000000000000 0.100000000000000 X false Rotation: false Z false Rotation angle (in degree) around z-axis. -360.000000000000000 360.000000000000000 0.100000000000000 Translation (in mm) in direction of the x-axis. -1000.000000000000000 1000.000000000000000 0.100000000000000 Scaling: false Scaling factor for selected fiber bundle along the x-axis. 0.010000000000000 10.000000000000000 0.010000000000000 1.000000000000000 Scaling factor for selected fiber bundle along the y-axis. 0.010000000000000 10.000000000000000 0.010000000000000 1.000000000000000 Scaling factor for selected fiber bundle along the z-axis. 0.010000000000000 10.000000000000000 0.010000000000000 1.000000000000000 false Join Bundles If checked, the fiducials belonging to the modified bundle are also modified. Include Fiducials true Fiber Options QFrame::NoFrame QFrame::Raised 0 QFrame::NoFrame QFrame::Raised 0 Tension: false Fiber Sampling: false 3 -1.000000000000000 1.000000000000000 0.100000000000000 0.000000000000000 Fiber sampling points (per cm) 1 100 1 10 3 -1.000000000000000 1.000000000000000 0.100000000000000 0.000000000000000 Bias: false Continuity: false 3 -1.000000000000000 1.000000000000000 0.100000000000000 0.000000000000000 QFrame::NoFrame QFrame::Raised 0 6 #Fibers: false Specify number of fibers to generate for the selected bundle. 1 1000000 100 100 false Generate Fibers QFrame::NoFrame QFrame::Raised 0 Select fiber distribution inside of the fiducials. Uniform Gaussian Fiber Distribution: false Variance of the gaussian 3 0.001000000000000 10.000000000000000 0.010000000000000 0.100000000000000 QFrame::NoFrame QFrame::Raised 0 Disable to only generate fibers if "Generate Fibers" button is pressed. Real Time Fibers true Disable to only generate fibers if "Generate Fibers" button is pressed. Advanced Options false QFrame::NoFrame QFrame::Raised 0 false 30 30 Draw elliptical fiducial. :/QmitkDiffusionImaging/circle.png:/QmitkDiffusionImaging/circle.png 32 32 false true false 30 30 Flip fiber waypoints of selcted fiducial around one axis. :/QmitkDiffusionImaging/refresh.xpm:/QmitkDiffusionImaging/refresh.xpm 32 32 false true Qt::Horizontal 40 20 Signal Generation Intra-axonal Compartment Select signal model for intra-axonal compartment. Stick Model Zeppelin Model Tensor Model Data Fiber Bundle: false <html><head/><body><p><span style=" color:#ff0000;">mandatory</span></p></body></html> true Tissue Mask: false <html><head/><body><p><span style=" color:#969696;">optional</span></p></body></html> true Extra-axonal Compartments Select signal model for extra-axonal compartment. Ball Model Astrosticks Model Dot Model Select signal model for extra-axonal compartment. -- Ball Model Astrosticks Model Dot Model Qt::Horizontal QFrame::NoFrame QFrame::Raised 0 Weighting factor between the two extra-axonal compartments. 1.000000000000000 0.100000000000000 0.300000000000000 Compartment Fraction: true Start DWI generation from selected fiebr bundle. If no fiber bundle is selected, a grayscale image containing a simple gradient is generated. Generate Image Image Settings QFrame::NoFrame QFrame::Raised 0 6 <html><head/><body><p><span style=" font-style:italic;">TE</span>, <span style=" font-style:italic;">T</span><span style=" font-style:italic; vertical-align:sub;">inhom</span> and <span style=" font-style:italic;">T2</span> will have no effect if unchecked.</p></body></html> Simulate Signal Relaxation true Repetitions: T2* relaxation time (in milliseconds). 100.000000000000000 0.100000000000000 1.000000000000000 Fiber Radius: Fiber radius used to calculate volume fractions (in µm). Set to 0 for automatic radius estimation. 0 1000 0 TE in milliseconds 1 10000 1 100 Interpolation Shrink: Line Readout Time: false <html><head/><body><p>Echo Time <span style=" font-style:italic;">TE</span>: </p></body></html> false Disable partial volume. Treat voxel content as fiber-only if at least one fiber is present. Disable Partial Volume Effects false Output one image per compartment containing the corresponding volume fractions per voxel. Output Volume Fractions false <html><head/><body><p><span style=" font-style:italic;">T</span><span style=" font-style:italic; vertical-align:sub;">inhom</span> Relaxation: </p></body></html> false Number of signal averages. Increase to reduce noise. 1 100 1 1 Relaxation time due to magnetic field inhomogeneities (T2', in milliseconds). 1 10000 1 50 TE in milliseconds 1 10000 1 100 <html><head/><body><p>Large values shrink (towards nearest neighbour interpolation), small values strech interpolation function (towards linear interpolation). 1000 equals nearest neighbour interpolation.</p></body></html> 1 1000 1000 Signal Scale: color: rgb(255, 0, 0); Using geometry of selected image! color: rgb(255, 0, 0); Using gradients of selected DWI! QFrame::NoFrame QFrame::Raised 0 3 0.100000000000000 50.000000000000000 0.100000000000000 2.000000000000000 Image Spacing: 3 0.100000000000000 50.000000000000000 0.100000000000000 2.000000000000000 3 0.100000000000000 50.000000000000000 0.100000000000000 2.000000000000000 Image Dimensions: Fiber sampling factor which determines the accuracy of the calculated fiber and non-fiber volume fractions. 1 1000 1 11 Fiber sampling factor which determines the accuracy of the calculated fiber and non-fiber volume fractions. 1 1000 1 11 Fiber sampling factor which determines the accuracy of the calculated fiber and non-fiber volume fractions. 1 1000 1 3 QFrame::NoFrame QFrame::Raised 0 6 Gradient Directions: Number of gradient directions distributed over the half sphere. 0 10000 1 30 b-Value: false b-value in mm/s² 0 10000 100 1000 Advanced Options Qt::Vertical 20 40 Noise and other Artifacts + + + + Add N/2 Ghosts + + + false + + + true QFrame::NoFrame QFrame::Raised 6 0 K-Space Line Offset: false A larger offset increases the inensity of the ghost image. 3 1.000000000000000 0.010000000000000 0.100000000000000 QFrame::NoFrame QFrame::Raised 0 Variance: Variance of Rician noise model. 4 0.000000000000000 100000.000000000000000 0.001000000000000 25.000000000000000 + + Add ringing artifacts occuring at strong edges in the image. + Add Gibbs Ringing false Add Rician Noise false - + true QFrame::NoFrame QFrame::Raised 6 0 - k-Space Undersampling: + Upsampling: false - + - Image is upsampled using this factor, afterwards fourier transformed, cropped to the original size and then inverse fourier transformed. + Larger values increase the ringing range. Beware of performance issues! - - 0 + + 2 + + + 1.000000000000000 + + + 10.000000000000000 + + + 0.100000000000000 + + + 2.000000000000000 - - - 2 - - - - - 4 - - - - - 8 - - - - - 16 - - - - - 32 - - - - - 64 - - - - - 128 - - - - - 256 - - Add Distortions false true QFrame::NoFrame QFrame::Raised QFormLayout::AllNonFixedFieldsGrow 6 0 Gradient Strength: false A larger offset increases the inensity of the ghost image. 5 1000.000000000000000 0.001000000000000 - 0.015000000000000 + 0.001000000000000 true QFrame::NoFrame QFrame::Raised 6 0 Frequency Map: false Select image specifying the frequency inhomogeneities (in Hz). - - - - Add N/2 Ghosts - - - false - - - + + !!!EXPERIMENTAL!!! + Add Eddy Current Effects false Inter-axonal Compartment Select signal model for intra-axonal compartment. -- Stick Model Zeppelin Model Tensor Model QmitkDataStorageComboBox QComboBox
QmitkDataStorageComboBox.h
 QmitkTensorModelParametersWidget QWidget
QmitkTensorModelParametersWidget.h
1 QmitkStickModelParametersWidget QWidget
QmitkStickModelParametersWidget.h
1 QmitkZeppelinModelParametersWidget QWidget
QmitkZeppelinModelParametersWidget.h
1 QmitkBallModelParametersWidget QWidget
QmitkBallModelParametersWidget.h
1 QmitkAstrosticksModelParametersWidget QWidget
QmitkAstrosticksModelParametersWidget.h
1 QmitkDotModelParametersWidget QWidget
QmitkDotModelParametersWidget.h
1 m_CircleButton m_FlipButton m_RealTimeFibers m_AdvancedOptionsBox m_DistributionBox m_VarianceBox m_FiberDensityBox m_FiberSamplingBox m_TensionBox m_ContinuityBox m_BiasBox m_GenerateFibersButton m_ConstantRadiusBox m_AlignOnGrid m_XrotBox m_YrotBox m_ZrotBox m_XtransBox m_YtransBox m_ZtransBox m_XscaleBox m_YscaleBox m_ZscaleBox m_TransformBundlesButton m_CopyBundlesButton m_JoinBundlesButton m_IncludeFiducials m_GenerateImageButton m_SizeX m_SizeY m_SizeZ m_SpacingX m_SpacingY m_SpacingZ m_NumGradientsBox m_BvalueBox m_AdvancedOptionsBox_2 m_RepetitionsBox m_SignalScaleBox m_TEbox m_LineReadoutTimeBox m_T2starBox m_FiberRadius m_InterpolationShrink m_EnforcePureFiberVoxelsBox m_Compartment1Box m_Compartment2Box m_Compartment3Box m_Compartment4Box m_Comp4FractionBox m_AddNoise m_NoiseLevel m_AddGhosts m_kOffsetBox m_AddDistortions m_FrequencyMapBox m_AddGibbsRinging - m_KspaceUndersamplingBox tabWidget