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/itkFieldmapGeneratorFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFieldmapGeneratorFilter.cpp index f6f7516e67..15c27ca564 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFieldmapGeneratorFilter.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFieldmapGeneratorFilter.cpp @@ -1,79 +1,107 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Coindex[1]right (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 "itkFieldmapGeneratorFilter.h" #include #include #include namespace itk{ template< class OutputImageType > FieldmapGeneratorFilter< OutputImageType >::FieldmapGeneratorFilter() { m_Gradient.fill(0.0); m_Offset.fill(0.0); } template< class OutputImageType > FieldmapGeneratorFilter< OutputImageType >::~FieldmapGeneratorFilter() { } template< class OutputImageType > void FieldmapGeneratorFilter< OutputImageType >::BeforeThreadedGenerateData() { typename 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->Allocate(); outImage->FillBuffer(0); this->SetNthOutput(0, outImage); } template< class OutputImageType > void FieldmapGeneratorFilter< OutputImageType >::ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, ThreadIdType threadId) { typename OutputImageType::Pointer outImage = static_cast< OutputImageType * >(this->ProcessObject::GetOutput(0)); ImageRegionIterator< OutputImageType > oit(outImage, outputRegionForThread); + int szx = m_ImageRegion.GetSize(0); + int szy = m_ImageRegion.GetSize(1); + int szz = m_ImageRegion.GetSize(2); + + double max = szx*szx + szy*szy - szx*szy; + while( !oit.IsAtEnd() ) { double value = 0; IndexType idx = oit.GetIndex(); - for (int i=0; i<3; i++) - value += idx[i]*m_Gradient[i] + m_Offset[i]; +// for (int i=0; i<3; i++) +// value += std::exp(-(double)idx[i]/32.0)*m_Gradient[i] + m_Offset[i]; + + //value += idx[0]*idx[0] + idx[1]*idx[1] - idx[0]*idx[1]; for (int i=0; i vertex; outImage->TransformIndexToPhysicalPoint(idx, vertex); double dist = c.EuclideanDistanceTo(vertex); - value += m_Heights.at(i)*exp(-dist*dist/(2*m_Variances.at(i))); + + c[0] -= vertex[0]; + c[1] -= vertex[1]; + c[2] -= vertex[2]; + + double x = c[0]; + double y = c[1]; + double z = c[2]; + + dist = fabs(dist) - m_Variances.at(i); + if (dist>0) + value += (2*z*z-y*y-x*x)/pow(x*x+y*y+z*z,5/2); + + +// dist = fabs(dist) - m_Variances.at(i); +// if (dist<0) +// dist = 0; +// value += std::exp(-dist/m_Variances.at(i))*m_Heights.at(i); + + // value += m_Heights.at(i)*exp(-dist*dist/(2*m_Variances.at(i))); } - oit.Set(value); + + oit.Set(m_Gradient[0]*value/max); ++oit; } MITK_INFO << "Thread " << threadId << "finished processing"; } } 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..a432420dcb --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToDWIImageFilter.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToDWIImageFilter.cpp @@ -1,709 +1,699 @@ /*=================================================================== 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)); + MITK_INFO << outSize; 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..860206a99f --- /dev/null +++ b/Modules/DiffusionImaging/FiberTracking/MiniApps/CMakeLists.txt @@ -0,0 +1,47 @@ +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 + ) + + # 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..c527d8367d --- /dev/null +++ b/Modules/DiffusionImaging/FiberTracking/MiniApps/GibbsTracking.cpp @@ -0,0 +1,219 @@ +/*=================================================================== + +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 + +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[]) +{ + if (argc<6) + { + std::cout << std::endl; + std::cout << "Performes Gibbs Tracking" << std::endl; + std::cout << std::endl; + std::cout << "usage: " << argv[0] << " " << std::endl; + std::cout << std::endl; + return EXIT_FAILURE; + } + + try + { + RegisterDiffusionCoreObjectFactory(); + RegisterFiberTrackingObjectFactory(); + + std::cout << "input: " << argv[1] << std::endl; + std::cout << "param: " << argv[2] << std::endl; + std::cout << "output: " << argv[3] << std::endl; + + // 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( argv[1], s1, s2, false ); + + // try to cast to qball image + QString inImageName(argv[1]); + 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(argc==6) + { + QString convention(argv[5]); + 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(argc>=5) + { + typedef itk::Image MaskImgType; + std::cout << "mask: " << argv[4]<< std::endl; + infile = mitk::BaseDataIO::LoadBaseDataFromFile( argv[4], s1, s2, false ); + mitk::Image::Pointer mitkMaskImage = dynamic_cast(infile.at(0).GetPointer()); + MaskImgType::Pointer itk_mask = MaskImgType::New(); + mitk::CastToItkImage(mitkMaskImage, itk_mask); + gibbsTracker->SetMaskImage(itk_mask); + } + else{ + std::cout << "perform tracking without mask" << std::endl; + } + + gibbsTracker->SetDuplicateImage(false); + gibbsTracker->SetLoadParameterFile( argv[2] ); +// gibbsTracker->SetLutPath( "" ); + gibbsTracker->Update(); + + mitk::FiberBundleX::Pointer mitkFiberBundle = mitk::FiberBundleX::New(gibbsTracker->GetFiberBundle()); + + std::string outfilename = argv[3]; + MITK_INFO << "searching writer"; + 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..a27128840b --- /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[argc-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/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..5f6ac6b882 100644 --- a/Modules/DiffusionImaging/FiberTracking/files.cmake +++ b/Modules/DiffusionImaging/FiberTracking/files.cmake @@ -1,101 +1,100 @@ set(CPP_FILES # 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 # 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 # (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/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