diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToDWIImageFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToDWIImageFilter.cpp index 56b0cbf3db..a2f8596050 100755 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToDWIImageFilter.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToDWIImageFilter.cpp @@ -1,781 +1,849 @@ /*=================================================================== 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 #include #include #include #include #include +#include +#include namespace itk { template< class PixelType > TractsToDWIImageFilter< PixelType >::TractsToDWIImageFilter() : m_CircleDummy(false) , m_VolumeAccuracy(10) , m_AddGibbsRinging(false) , 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_Spikes(0) , m_Wrap(1.0) , m_SpikeAmplitude(1) , m_AddMotionArtifact(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); m_MaxTranslation.Fill(0.0); m_MaxRotation.Fill(0.0); m_RandGen = itk::Statistics::MersenneTwisterRandomVariateGenerator::New(); m_RandGen->SetSeed(); } template< class PixelType > TractsToDWIImageFilter< PixelType >::~TractsToDWIImageFilter() { } template< class PixelType > TractsToDWIImageFilter< PixelType >::DoubleDwiType::Pointer TractsToDWIImageFilter< PixelType >::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]); Vector< double, 2 > sliceSpacing; sliceSpacing[0] = m_UpsampledSpacing[0]; sliceSpacing[1] = m_UpsampledSpacing[1]; // frequency map slice SliceType::Pointer fMapSlice = NULL; if (m_FrequencyMap.IsNotNull()) { fMapSlice = SliceType::New(); ImageRegion<2> region; region.SetSize(0, m_UpsampledImageRegion.GetSize()[0]); region.SetSize(1, m_UpsampledImageRegion.GetSize()[1]); fMapSlice->SetLargestPossibleRegion( region ); fMapSlice->SetBufferedRegion( region ); fMapSlice->SetRequestedRegion( region ); fMapSlice->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]; } std::vector< int > spikeVolume; for (int i=0; iGetVectorLength()); std::sort (spikeVolume.begin(), spikeVolume.end()); std::reverse (spikeVolume.begin(), spikeVolume.end()); boost::progress_display disp(images.at(0)->GetVectorLength()*images.at(0)->GetLargestPossibleRegion().GetSize(2)); for (unsigned int g=0; gGetVectorLength(); g++) { std::vector< int > spikeSlice; while (!spikeVolume.empty() && spikeVolume.back()==g) { spikeSlice.push_back(rand()%images.at(0)->GetLargestPossibleRegion().GetSize(2)); spikeVolume.pop_back(); } std::sort (spikeSlice.begin(), spikeSlice.end()); std::reverse (spikeSlice.begin(), spikeSlice.end()); for (unsigned int z=0; zGetLargestPossibleRegion().GetSize(2); z++) { std::vector< SliceType::Pointer > compartmentSlices; std::vector< double > t2Vector; for (unsigned int i=0; i* signalModel; if (iSetLargestPossibleRegion( sliceRegion ); slice->SetBufferedRegion( sliceRegion ); slice->SetRequestedRegion( sliceRegion ); slice->SetSpacing(sliceSpacing); slice->Allocate(); slice->FillBuffer(0.0); // extract slice from channel g for (unsigned int y=0; yGetLargestPossibleRegion().GetSize(1); y++) for (unsigned 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 (fMapSlice.IsNotNull() && i==0) fMapSlice->SetPixel(index2D, m_FrequencyMap->GetPixel(index3D)); } compartmentSlices.push_back(slice); t2Vector.push_back(signalModel->GetT2()); } // create k-sapce (inverse fourier transform slices) itk::Size<2> outSize; outSize.SetElement(0, m_ImageRegion.GetSize(0)); outSize.SetElement(1, m_ImageRegion.GetSize(1)); itk::KspaceImageFilter< SliceType::PixelType >::Pointer idft = itk::KspaceImageFilter< SliceType::PixelType >::New(); idft->SetCompartmentImages(compartmentSlices); idft->SetT2(t2Vector); idft->SetkOffset(m_kOffset); idft->SettLine(m_tLine); idft->SetTE(m_TE); idft->SetTinhom(m_tInhom); idft->SetSimulateRelaxation(m_SimulateRelaxation); idft->SetSimulateEddyCurrents(m_SimulateEddyCurrents); idft->SetEddyGradientMagnitude(m_EddyGradientStrength); idft->SetZ((double)z-(double)images.at(0)->GetLargestPossibleRegion().GetSize(2)/2.0); idft->SetDirectionMatrix(transform); idft->SetDiffusionGradientDirection(m_FiberModels.at(0)->GetGradientDirection(g)); idft->SetFrequencyMap(fMapSlice); idft->SetSignalScale(m_SignalScale); idft->SetOutSize(outSize); int numSpikes = 0; while (!spikeSlice.empty() && spikeSlice.back()==z) { numSpikes++; spikeSlice.pop_back(); } idft->SetSpikes(numSpikes); idft->SetSpikeAmplitude(m_SpikeAmplitude); idft->Update(); ComplexSliceType::Pointer fSlice; fSlice = idft->GetOutput(); for (unsigned 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 (unsigned int y=0; yGetLargestPossibleRegion().GetSize(1); y++) for (unsigned 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; } template< class PixelType > void TractsToDWIImageFilter< PixelType >::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_EnforcePureFiberVoxels) while (m_FiberModels.size()>1) m_FiberModels.pop_back(); 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!"); // initialize output dwi image ImageRegion<3> croppedRegion = m_ImageRegion; croppedRegion.SetSize(1, croppedRegion.GetSize(1)*m_Wrap); itk::Point shiftedOrigin = m_Origin; shiftedOrigin[1] += (m_ImageRegion.GetSize(1)-croppedRegion.GetSize(1))*m_Spacing[1]/2; typename OutputImageType::Pointer outImage = OutputImageType::New(); outImage->SetSpacing( m_Spacing ); outImage->SetOrigin( shiftedOrigin ); outImage->SetDirection( m_DirectionMatrix ); outImage->SetLargestPossibleRegion( croppedRegion ); outImage->SetBufferedRegion( croppedRegion ); outImage->SetRequestedRegion( croppedRegion ); outImage->SetVectorLength( m_FiberModels[0]->GetNumGradients() ); outImage->Allocate(); typename OutputImageType::PixelType temp; temp.SetSize(m_FiberModels[0]->GetNumGradients()); temp.Fill(0.0); outImage->FillBuffer(temp); // ADJUST GEOMETRY FOR FURTHER PROCESSING // is input slize size a power of two? unsigned int x=m_ImageRegion.GetSize(0); unsigned int y=m_ImageRegion.GetSize(1); if ( x%2 == 1 ) m_ImageRegion.SetSize(0, x+1); if ( y%2 == 1 ) m_ImageRegion.SetSize(1, y+1); // apply in-plane upsampling double upsampling = 1; if (m_AddGibbsRinging) { MITK_INFO << "Adding ringing artifacts."; upsampling = 2; } m_UpsampledSpacing = m_Spacing; m_UpsampledSpacing[0] /= upsampling; m_UpsampledSpacing[1] /= upsampling; m_UpsampledImageRegion = m_ImageRegion; m_UpsampledImageRegion.SetSize(0, m_ImageRegion.GetSize()[0]*upsampling); m_UpsampledImageRegion.SetSize(1, m_ImageRegion.GetSize()[1]*upsampling); + m_UpsampledOrigin = m_Origin; + m_UpsampledOrigin[0] -= m_Spacing[0]/2; m_UpsampledOrigin[0] += m_UpsampledSpacing[0]/2; + m_UpsampledOrigin[1] -= m_Spacing[1]/2; m_UpsampledOrigin[1] += m_UpsampledSpacing[1]/2; + m_UpsampledOrigin[2] -= m_Spacing[2]/2; m_UpsampledOrigin[2] += m_UpsampledSpacing[2]/2; // generate double images to store the individual compartment signals std::vector< DoubleDwiType::Pointer > compartments; for (unsigned int i=0; iSetSpacing( m_UpsampledSpacing ); - doubleDwi->SetOrigin( m_Origin ); + doubleDwi->SetOrigin( m_UpsampledOrigin ); 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); } // initialize volume fraction images m_VolumeFractions.clear(); for (unsigned int i=0; iSetSpacing( m_UpsampledSpacing ); - doubleImg->SetOrigin( m_Origin ); + doubleImg->SetOrigin( m_UpsampledOrigin ); doubleImg->SetDirection( m_DirectionMatrix ); doubleImg->SetLargestPossibleRegion( m_UpsampledImageRegion ); doubleImg->SetBufferedRegion( m_UpsampledImageRegion ); doubleImg->SetRequestedRegion( m_UpsampledImageRegion ); doubleImg->Allocate(); doubleImg->FillBuffer(0); m_VolumeFractions.push_back(doubleImg); } // resample mask image and frequency map to fit upsampled geometry if (m_AddGibbsRinging) { if (m_TissueMask.IsNotNull()) { // rescale mask image (otherwise there are problems with the resampling) itk::RescaleIntensityImageFilter::Pointer rescaler = itk::RescaleIntensityImageFilter::New(); rescaler->SetInput(0,m_TissueMask); rescaler->SetOutputMaximum(100); rescaler->SetOutputMinimum(0); rescaler->Update(); // resample mask image itk::ResampleImageFilter::Pointer resampler = itk::ResampleImageFilter::New(); resampler->SetInput(rescaler->GetOutput()); resampler->SetOutputParametersFromImage(m_TissueMask); resampler->SetSize(m_UpsampledImageRegion.GetSize()); resampler->SetOutputSpacing(m_UpsampledSpacing); + resampler->SetOutputOrigin(m_UpsampledOrigin); resampler->Update(); m_TissueMask = resampler->GetOutput(); } // 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->SetOutputOrigin(m_UpsampledOrigin); resampler->Update(); m_FrequencyMap = resampler->GetOutput(); } } + // no input tissue mask is set -> create default + bool maskImageSet = true; if (m_TissueMask.IsNull()) { + MITK_INFO << "No tissue mask set"; m_TissueMask = ItkUcharImgType::New(); m_TissueMask->SetSpacing( m_UpsampledSpacing ); - m_TissueMask->SetOrigin( m_Origin ); + m_TissueMask->SetOrigin( m_UpsampledOrigin ); 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); + maskImageSet = false; } + else + MITK_INFO << "Using tissue mask"; m_ImageRegion = croppedRegion; x=m_ImageRegion.GetSize(0); y=m_ImageRegion.GetSize(1); if ( x%2 == 1 ) m_ImageRegion.SetSize(0, x+1); if ( y%2 == 1 ) m_ImageRegion.SetSize(1, y+1); // 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; double interpFact = 2*atan(-0.5*m_InterpolationShrink); double maxVolume = 0; double voxelVolume = m_UpsampledSpacing[0]*m_UpsampledSpacing[1]*m_UpsampledSpacing[2]; MITK_INFO << "Generating signal of " << m_FiberModels.size() << " fiber compartments"; boost::progress_display disp(numFibers*m_FiberModels.at(0)->GetNumGradients()); if (m_AddMotionArtifact) { MITK_INFO << "Adding motion artifacts"; MITK_INFO << "Maximum rotation: " << m_MaxRotation; MITK_INFO << "Maxmimum translation: " << m_MaxTranslation; } maxVolume = 0; ofstream logFile; logFile.open("fiberfox_motion.log"); - logFile << "0,0,0,0\n"; + logFile << "0 rotation: 0,0,0; translation: 0,0,0\n"; // astrosticks überarbeiten!!!!!!!!!!!!!!!!!!!!! for (unsigned int i=0; i* >(m_NonFiberModels.at(i))) { mitk::AstroStickModel* model = dynamic_cast< mitk::AstroStickModel* >(m_NonFiberModels.at(i)); model->SetSeed(8111984); } + // get transform for motion artifacts FiberBundleType fiberBundleTransformed = fiberBundle; + VectorType rotation = m_MaxRotation/m_FiberModels.at(0)->GetNumGradients(); + VectorType translation = m_MaxTranslation/m_FiberModels.at(0)->GetNumGradients(); + + ItkUcharImgType::Pointer tempTissueMask = ItkUcharImgType::New(); + itk::ImageDuplicator::Pointer duplicator = itk::ImageDuplicator::New(); + duplicator->SetInputImage(m_TissueMask); + duplicator->Update(); + tempTissueMask = duplicator->GetOutput(); + + ImageRegion<3> upsampledImageRegion = m_UpsampledImageRegion; + itk::Vector upsampledSpacing = m_UpsampledSpacing; + upsampledSpacing[0] /= 4; + upsampledSpacing[1] /= 4; + upsampledSpacing[2] /= 4; + upsampledImageRegion.SetSize(0, m_UpsampledImageRegion.GetSize()[0]*4); + upsampledImageRegion.SetSize(1, m_UpsampledImageRegion.GetSize()[1]*4); + upsampledImageRegion.SetSize(2, m_UpsampledImageRegion.GetSize()[2]*4); + itk::Point upsampledOrigin = m_UpsampledOrigin; + upsampledOrigin[0] -= m_UpsampledSpacing[0]/2; upsampledOrigin[0] += upsampledSpacing[0]/2; + upsampledOrigin[1] -= m_UpsampledSpacing[1]/2; upsampledOrigin[1] += upsampledSpacing[1]/2; + upsampledOrigin[2] -= m_UpsampledSpacing[2]/2; upsampledOrigin[2] += upsampledSpacing[2]/2; + ItkUcharImgType::Pointer upsampledTissueMask = ItkUcharImgType::New(); + itk::ResampleImageFilter::Pointer upsampler = itk::ResampleImageFilter::New(); + upsampler->SetInput(m_TissueMask); + upsampler->SetOutputParametersFromImage(m_TissueMask); + upsampler->SetSize(upsampledImageRegion.GetSize()); + upsampler->SetOutputSpacing(upsampledSpacing); + upsampler->SetOutputOrigin(upsampledOrigin); + itk::NearestNeighborInterpolateImageFunction::Pointer nn_interpolator + = itk::NearestNeighborInterpolateImageFunction::New(); + upsampler->SetInterpolator(nn_interpolator); + upsampler->Update(); + upsampledTissueMask = upsampler->GetOutput(); + for (int g=0; gGetNumGradients(); g++) { - vtkSmartPointer fiberPolyData = fiberBundleTransformed->GetFiberPolyData(); + vtkPolyData* fiberPolyData = fiberBundleTransformed->GetFiberPolyData(); ItkDoubleImgType::Pointer intraAxonalVolume = ItkDoubleImgType::New(); intraAxonalVolume->SetSpacing( m_UpsampledSpacing ); - intraAxonalVolume->SetOrigin( m_Origin ); + intraAxonalVolume->SetOrigin( m_UpsampledOrigin ); intraAxonalVolume->SetDirection( m_DirectionMatrix ); intraAxonalVolume->SetLargestPossibleRegion( m_UpsampledImageRegion ); intraAxonalVolume->SetBufferedRegion( m_UpsampledImageRegion ); intraAxonalVolume->SetRequestedRegion( m_UpsampledImageRegion ); intraAxonalVolume->Allocate(); intraAxonalVolume->FillBuffer(0); + // generate fiber signal for( int i=0; iGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); if (numPoints<2) continue; for( int j=0; jGetPoint(j); itk::Point vertex = GetItkPoint(temp); itk::Vector v = GetItkVector(temp); itk::Vector dir(3); if (jGetPoint(j+1))-v; else dir = v-GetItkVector(points->GetPoint(j-1)); if (dir.GetSquaredNorm()<0.0001) continue; itk::Index<3> idx; itk::ContinuousIndex contIndex; - m_TissueMask->TransformPhysicalPointToIndex(vertex, idx); - m_TissueMask->TransformPhysicalPointToContinuousIndex(vertex, contIndex); + tempTissueMask->TransformPhysicalPointToIndex(vertex, idx); + tempTissueMask->TransformPhysicalPointToContinuousIndex(vertex, contIndex); if (!m_UseInterpolation) // use nearest neighbour interpolation { - if (!m_TissueMask->GetLargestPossibleRegion().IsInside(idx) || m_TissueMask->GetPixel(idx)<=0) + if (!tempTissueMask->GetLargestPossibleRegion().IsInside(idx) || tempTissueMask->GetPixel(idx)<=0) continue; // generate signal for each fiber compartment for (unsigned int k=0; kSetFiberDirection(dir); DoubleDwiType::PixelType pix = doubleDwi->GetPixel(idx); pix[g] += segmentVolume*m_FiberModels[k]->SimulateMeasurement(g); doubleDwi->SetPixel(idx, pix ); double vol = intraAxonalVolume->GetPixel(idx) + segmentVolume; intraAxonalVolume->SetPixel(idx, vol ); if (g==0 && vol>maxVolume) maxVolume = vol; } continue; } - double frac_x = contIndex[0] - idx[0]; - double frac_y = contIndex[1] - idx[1]; - double frac_z = contIndex[2] - idx[2]; + 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) + if (!tempTissueMask->GetLargestPossibleRegion().IsInside(newIdx) || tempTissueMask->GetPixel(newIdx)<=0) continue; // generate signal for each fiber compartment for (unsigned int k=0; kSetFiberDirection(dir); DoubleDwiType::PixelType pix = doubleDwi->GetPixel(newIdx); pix[g] += segmentVolume*frac*m_FiberModels[k]->SimulateMeasurement(g); doubleDwi->SetPixel(newIdx, pix ); double vol = intraAxonalVolume->GetPixel(idx) + segmentVolume; intraAxonalVolume->SetPixel(idx, vol ); if (g==0 && vol>maxVolume) maxVolume = vol; } } } } } ++disp; } - ImageRegionIterator it3(m_TissueMask, m_TissueMask->GetLargestPossibleRegion()); - + // generate non-fiber signal + ImageRegionIterator it3(tempTissueMask, tempTissueMask->GetLargestPossibleRegion()); double fact = 1; if (m_FiberRadius<0.0001) fact = voxelVolume/maxVolume; - while(!it3.IsAtEnd()) { if (it3.Get()>0) { DoubleDwiType::IndexType index = it3.GetIndex(); // 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[g] *= fact; fiberDwi->SetPixel(index, fiberPix); } double f = intraAxonalVolume->GetPixel(index)*fact; if (f>voxelVolume || (f>0.0 && m_EnforcePureFiberVoxels) ) // more fiber than space in voxel? { fiberPix[g] *= voxelVolume/f; fiberDwi->SetPixel(index, fiberPix); m_VolumeFractions.at(0)->SetPixel(index, 1); } else { m_VolumeFractions.at(0)->SetPixel(index, f/voxelVolume); double nonf = voxelVolume-f; // non-fiber volume double inter = 0; if (m_FiberModels.size()>1) inter = nonf * f/voxelVolume; // inter-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 (unsigned int i=1; iGetPixel(index); if (f>0) pix[g] /= f; pix[g] *= singleinter; doubleDwi->SetPixel(index, pix); m_VolumeFractions.at(i)->SetPixel(index, singleinter/voxelVolume); } for (unsigned int i=0; iGetPixel(index); pix[g] += m_NonFiberModels[i]->SimulateMeasurement(g)*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; } + // move fibers if (m_AddMotionArtifact) { if (m_RandomMotion) { fiberBundleTransformed = fiberBundle->GetDeepCopy(); - - double rotX = m_RandGen->GetVariateWithClosedRange(m_MaxRotation[0]*2)-m_MaxRotation[0]; - double rotY = m_RandGen->GetVariateWithClosedRange(m_MaxRotation[1]*2)-m_MaxRotation[1]; - double rotZ = m_RandGen->GetVariateWithClosedRange(m_MaxRotation[2]*2)-m_MaxRotation[2]; - fiberBundleTransformed->TranslateFibers(rotX, rotY, rotZ); - logFile << g+1 << " rotation:" << rotX << "," << rotY << "," << rotZ << "\n"; - - double transX = m_RandGen->GetVariateWithClosedRange(m_MaxTranslation[0]*2)-m_MaxTranslation[0]; - double transY = m_RandGen->GetVariateWithClosedRange(m_MaxTranslation[1]*2)-m_MaxTranslation[1]; - double transZ = m_RandGen->GetVariateWithClosedRange(m_MaxTranslation[2]*2)-m_MaxTranslation[2]; - fiberBundleTransformed->TranslateFibers(transX, transY, transZ); - logFile << g+1 << " translation:" << transX << "," << transY << "," << transZ << "\n"; + rotation[0] = m_RandGen->GetVariateWithClosedRange(m_MaxRotation[0]*2)-m_MaxRotation[0]; + rotation[1] = m_RandGen->GetVariateWithClosedRange(m_MaxRotation[1]*2)-m_MaxRotation[1]; + rotation[2] = m_RandGen->GetVariateWithClosedRange(m_MaxRotation[2]*2)-m_MaxRotation[2]; + translation[0] = m_RandGen->GetVariateWithClosedRange(m_MaxTranslation[0]*2)-m_MaxTranslation[0]; + translation[1] = m_RandGen->GetVariateWithClosedRange(m_MaxTranslation[1]*2)-m_MaxTranslation[1]; + translation[2] = m_RandGen->GetVariateWithClosedRange(m_MaxTranslation[2]*2)-m_MaxTranslation[2]; } - else + + // rotate mask image + if (maskImageSet) { - logFile << g+1 << " rotation:" << m_MaxRotation[0]/m_FiberModels.at(0)->GetNumGradients() << ", " << m_MaxRotation[1]/m_FiberModels.at(0)->GetNumGradients() << ", " << m_MaxRotation[2]/m_FiberModels.at(0)->GetNumGradients(); - fiberBundleTransformed->RotateAroundAxis(m_MaxRotation[0]/m_FiberModels.at(0)->GetNumGradients(), m_MaxRotation[1]/m_FiberModels.at(0)->GetNumGradients(), m_MaxRotation[2]/m_FiberModels.at(0)->GetNumGradients()); + ImageRegionIterator maskIt(upsampledTissueMask, upsampledTissueMask->GetLargestPossibleRegion()); + tempTissueMask->FillBuffer(0); - logFile << g+1 << " translation:" << m_MaxTranslation[0]/m_FiberModels.at(0)->GetNumGradients() << ", " << m_MaxTranslation[1]/m_FiberModels.at(0)->GetNumGradients() << ", " << m_MaxTranslation[2]/m_FiberModels.at(0)->GetNumGradients(); - fiberBundleTransformed->TranslateFibers(m_MaxTranslation[0]/m_FiberModels.at(0)->GetNumGradients(),m_MaxTranslation[1]/m_FiberModels.at(0)->GetNumGradients(),m_MaxTranslation[2]/m_FiberModels.at(0)->GetNumGradients()); + while(!maskIt.IsAtEnd()) + { + if (maskIt.Get()<=0) + { + ++maskIt; + continue; + } + + DoubleDwiType::IndexType index = maskIt.GetIndex(); + itk::Point point; + upsampledTissueMask->TransformIndexToPhysicalPoint(index, point); + if (m_RandomMotion) + point = fiberBundle->TransformPoint(point.GetVnlVector(), rotation[0],rotation[1],rotation[2],translation[0],translation[1],translation[2]); + else + point = fiberBundle->TransformPoint(point.GetVnlVector(), rotation[0]*(g+1),rotation[1]*(g+1),rotation[2]*(g+1),translation[0]*(g+1),translation[1]*(g+1),translation[2]*(g+1)); + + tempTissueMask->TransformPhysicalPointToIndex(point, index); + if (tempTissueMask->GetLargestPossibleRegion().IsInside(index)) + tempTissueMask->SetPixel(index,100); + ++maskIt; + } } + + // rotate fibers + logFile << g+1 << " rotation:" << rotation[0] << "," << rotation[1] << "," << rotation[2] << ";"; + logFile << " translation:" << translation[0] << "," << translation[1] << "," << translation[2] << "\n"; + fiberBundleTransformed->TransformFibers(rotation[0],rotation[1],rotation[2],translation[0],translation[1],translation[2]); } } logFile.close(); // do k-space stuff DoubleDwiType::Pointer doubleOutImage; if (m_Spikes>0 || m_FrequencyMap.IsNotNull() || !m_KspaceArtifacts.empty() || m_kOffset>0 || m_SimulateRelaxation || m_SimulateEddyCurrents || m_AddGibbsRinging || m_Wrap<1.0) { MITK_INFO << "Adjusting complex signal"; doubleOutImage = DoKspaceStuff(compartments); m_SignalScale = 1; } else { MITK_INFO << "Summing compartments"; doubleOutImage = compartments.at(0); for (unsigned 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; typename OutputImageType::IndexType index = it4.GetIndex(); signal = doubleOutImage->GetPixel(index)*m_SignalScale; if (m_NoiseModel->GetNoiseVariance() > 0) { DoubleDwiType::PixelType accu = signal; accu.Fill(0.0); for (unsigned int i=0; iAddNoise(temp); accu += temp; } signal = accu/m_NumberOfRepetitions; } for (unsigned 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); } template< class PixelType > itk::Point TractsToDWIImageFilter< PixelType >::GetItkPoint(double point[3]) { itk::Point itkPoint; itkPoint[0] = point[0]; itkPoint[1] = point[1]; itkPoint[2] = point[2]; return itkPoint; } template< class PixelType > itk::Vector TractsToDWIImageFilter< PixelType >::GetItkVector(double point[3]) { itk::Vector itkVector; itkVector[0] = point[0]; itkVector[1] = point[1]; itkVector[2] = point[2]; return itkVector; } template< class PixelType > vnl_vector_fixed TractsToDWIImageFilter< PixelType >::GetVnlVector(double point[3]) { vnl_vector_fixed vnlVector; vnlVector[0] = point[0]; vnlVector[1] = point[1]; vnlVector[2] = point[2]; return vnlVector; } template< class PixelType > vnl_vector_fixed TractsToDWIImageFilter< PixelType >::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 index cd4b8e3703..9cb14a49be 100755 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToDWIImageFilter.h +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToDWIImageFilter.h @@ -1,167 +1,168 @@ /*=================================================================== 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 // ITK #include #include #include #include #include #include namespace itk { /** * \brief Generates artificial diffusion weighted image volume from the input fiberbundle using a generic multicompartment model. */ template< class PixelType > class TractsToDWIImageFilter : public ImageSource< itk::VectorImage< PixelType, 3 > > { public: typedef TractsToDWIImageFilter Self; typedef ImageSource< itk::VectorImage< PixelType, 3 > > Superclass; typedef SmartPointer< Self > Pointer; typedef SmartPointer< const Self > ConstPointer; typedef typename Superclass::OutputImageType OutputImageType; 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; typedef itk::Vector< double,3> VectorType; typedef itk::Point< double,3> PointType; 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, VectorType ) ///< output image spacing itkSetMacro( Origin, PointType ) ///< 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( AddGibbsRinging, bool ) itkSetMacro( Spikes, int ) itkSetMacro( SpikeAmplitude, double ) itkSetMacro( Wrap, double ) itkSetMacro( MaxTranslation, VectorType ) itkSetMacro( MaxRotation, VectorType ) itkSetMacro( AddMotionArtifact, bool ) itkSetMacro( RandomMotion, bool ) // 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); itk::Vector m_Spacing; ///< output image spacing itk::Vector m_UpsampledSpacing; itk::Point m_Origin; ///< output image origin + itk::Point m_UpsampledOrigin; 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; bool m_AddGibbsRinging; ///< 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; int m_Spikes; double m_SpikeAmplitude; double m_Wrap; VectorType m_MaxTranslation; VectorType m_MaxRotation; bool m_AddMotionArtifact; bool m_RandomMotion; itk::Statistics::MersenneTwisterRandomVariateGenerator::Pointer m_RandGen; }; } #include "itkTractsToDWIImageFilter.cpp" #endif diff --git a/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundleX/mitkFiberBundleX.cpp b/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundleX/mitkFiberBundleX.cpp index 99602169e8..15e7b36cab 100755 --- a/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundleX/mitkFiberBundleX.cpp +++ b/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundleX/mitkFiberBundleX.cpp @@ -1,1791 +1,1926 @@ /*=================================================================== 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 "mitkImagePixelReadAccessor.h" #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) { - mitkPixelTypeMultiplex1( SetFAMap, FAimage->GetPixelType(), FAimage ); + mitkPixelTypeMultiplex1( SetFAMap, FAimage->GetPixelType(), FAimage ); } template void mitk::FiberBundleX::SetFAMap(const mitk::PixelType pixelType, 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()); - - mitk::ImagePixelReadAccessor readFAimage (FAimage, FAimage->GetVolumeData(0)); - - 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-readFAimage.GetPixelByWorldCoordinates(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"; + MITK_DEBUG << "SetFAMap"; + vtkSmartPointer faValues = vtkSmartPointer::New(); + faValues->SetName(COLORCODING_FA_BASED); + faValues->Allocate(m_FiberPolyData->GetNumberOfPoints()); + faValues->SetNumberOfValues(m_FiberPolyData->GetNumberOfPoints()); + + mitk::ImagePixelReadAccessor readFAimage (FAimage, FAimage->GetVolumeData(0)); + + 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-readFAimage.GetPixelByWorldCoordinates(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 = cellOriginal->GetNumberOfPoints(); vtkPoints* pointsOriginal = cellOriginal->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); 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 && mask->GetLargestPossibleRegion().IsInside(idx) ) { for (int k=0; kGetPoint(k); vtkIdType id = vtkNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } break; } } } else { double* start = pointsOriginal->GetPoint(0); itk::Point itkStart; itkStart[0] = start[0]; itkStart[1] = start[1]; itkStart[2] = start[2]; itk::Index<3> idxStart; mask->TransformPhysicalPointToIndex(itkStart, idxStart); double* end = pointsOriginal->GetPoint(numPointsOriginal-1); itk::Point itkEnd; itkEnd[0] = end[0]; itkEnd[1] = end[1]; itkEnd[2] = end[2]; itk::Index<3> idxEnd; mask->TransformPhysicalPointToIndex(itkEnd, idxEnd); if ( mask->GetPixel(idxStart)>0 && mask->GetPixel(idxEnd)>0 && mask->GetLargestPossibleRegion().IsInside(idxStart) && mask->GetLargestPossibleRegion().IsInside(idxEnd) ) { for (int j=0; jGetPoint(j); vtkIdType id = vtkNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } } } } 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::RemoveFibersOutside(ItkUcharImgType* mask, bool invert) { float minSpacing = 1; if(mask->GetSpacing()[0]GetSpacing()[1] && mask->GetSpacing()[0]GetSpacing()[2]) minSpacing = mask->GetSpacing()[0]; else if (mask->GetSpacing()[1] < mask->GetSpacing()[2]) minSpacing = mask->GetSpacing()[1]; else minSpacing = mask->GetSpacing()[2]; mitk::FiberBundleX::Pointer fibCopy = this->GetDeepCopy(); fibCopy->ResampleFibers(minSpacing/10); vtkSmartPointer polyData =fibCopy->GetFiberPolyData(); vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); MITK_INFO << "Cutting fibers"; boost::progress_display disp(m_NumFibers); for (int i=0; iGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); if (numPoints>1) { int newNumPoints = 0; for (int j=0; jGetPoint(j); itk::Point itkP; itkP[0] = p[0]; itkP[1] = p[1]; itkP[2] = p[2]; itk::Index<3> idx; mask->TransformPhysicalPointToIndex(itkP, idx); if ( mask->GetPixel(idx)>0 && mask->GetLargestPossibleRegion().IsInside(idx) && !invert ) { vtkIdType id = vtkNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); newNumPoints++; } else if ( (mask->GetPixel(idx)<=0 || !mask->GetLargestPossibleRegion().IsInside(idx)) && invert ) { vtkIdType id = vtkNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); newNumPoints++; } else if (newNumPoints>0) { vtkNewCells->InsertNextCell(container); newNumPoints = 0; container = vtkSmartPointer::New(); } } if (newNumPoints>0) 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])) ; + (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); } std::vector mitk::FiberBundleX::GetAvailableColorCodings() { std::vector 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.empty()) 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 } } +itk::Matrix< double, 3, 3 > mitk::FiberBundleX::TransformMatrix(itk::Matrix< double, 3, 3 > m, double rx, double ry, double rz) +{ + rx = rx*M_PI/180; + ry = ry*M_PI/180; + rz = rz*M_PI/180; + + itk::Matrix< double, 3, 3 > rotX; rotX.SetIdentity(); + rotX[1][1] = cos(rx); + rotX[2][2] = rotX[1][1]; + rotX[1][2] = -sin(rx); + rotX[2][1] = -rotX[1][2]; + + itk::Matrix< double, 3, 3 > rotY; rotY.SetIdentity(); + rotY[0][0] = cos(ry); + rotY[2][2] = rotY[0][0]; + rotY[0][2] = sin(ry); + rotY[2][0] = -rotY[0][2]; + + itk::Matrix< double, 3, 3 > rotZ; rotZ.SetIdentity(); + rotZ[0][0] = cos(rz); + rotZ[1][1] = rotZ[0][0]; + rotZ[0][1] = -sin(rz); + rotZ[1][0] = -rotZ[0][1]; + + itk::Matrix< double, 3, 3 > rot = rotZ*rotY*rotX; + + m = rot*m; + + return m; +} + +itk::Point mitk::FiberBundleX::TransformPoint(vnl_vector_fixed< double, 3 > point, double rx, double ry, double rz, double tx, double ty, double tz) +{ + rx = rx*M_PI/180; + ry = ry*M_PI/180; + rz = rz*M_PI/180; + + vnl_matrix_fixed< double, 3, 3 > rotX; rotX.set_identity(); + rotX[1][1] = cos(rx); + rotX[2][2] = rotX[1][1]; + rotX[1][2] = -sin(rx); + rotX[2][1] = -rotX[1][2]; + + vnl_matrix_fixed< double, 3, 3 > rotY; rotY.set_identity(); + rotY[0][0] = cos(ry); + rotY[2][2] = rotY[0][0]; + rotY[0][2] = sin(ry); + rotY[2][0] = -rotY[0][2]; + + vnl_matrix_fixed< double, 3, 3 > rotZ; rotZ.set_identity(); + rotZ[0][0] = cos(rz); + rotZ[1][1] = rotZ[0][0]; + rotZ[0][1] = -sin(rz); + rotZ[1][0] = -rotZ[0][1]; + + vnl_matrix_fixed< double, 3, 3 > rot = rotZ*rotY*rotX; + + mitk::Geometry3D::Pointer geom = this->GetGeometry(); + mitk::Point3D center = geom->GetCenter(); + + point[0] -= center[0]; + point[1] -= center[1]; + point[2] -= center[2]; + point = rot*point; + point[0] += center[0]+tx; + point[1] += center[1]+ty; + point[2] += center[2]+tz; + itk::Point out; out[0] = point[0]; out[1] = point[1]; out[2] = point[2]; + return out; +} + +void mitk::FiberBundleX::TransformFibers(double rx, double ry, double rz, double tx, double ty, double tz) +{ + rx = rx*M_PI/180; + ry = ry*M_PI/180; + rz = rz*M_PI/180; + + vnl_matrix_fixed< double, 3, 3 > rotX; rotX.set_identity(); + rotX[1][1] = cos(rx); + rotX[2][2] = rotX[1][1]; + rotX[1][2] = -sin(rx); + rotX[2][1] = -rotX[1][2]; + + vnl_matrix_fixed< double, 3, 3 > rotY; rotY.set_identity(); + rotY[0][0] = cos(ry); + rotY[2][2] = rotY[0][0]; + rotY[0][2] = sin(ry); + rotY[2][0] = -rotY[0][2]; + + vnl_matrix_fixed< double, 3, 3 > rotZ; rotZ.set_identity(); + rotZ[0][0] = cos(rz); + rotZ[1][1] = rotZ[0][0]; + rotZ[0][1] = -sin(rz); + rotZ[1][0] = -rotZ[0][1]; + + vnl_matrix_fixed< double, 3, 3 > rot = rotZ*rotY*rotX; + + mitk::Geometry3D::Pointer geom = this->GetGeometry(); + mitk::Point3D center = geom->GetCenter(); + + vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); + vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); + + for (int i=0; iGetCell(i); + int numPoints = cell->GetNumberOfPoints(); + vtkPoints* points = cell->GetPoints(); + + vtkSmartPointer container = vtkSmartPointer::New(); + for (int j=0; jGetPoint(j); + vnl_vector_fixed< double, 3 > dir; + dir[0] = p[0]-center[0]; + dir[1] = p[1]-center[1]; + dir[2] = p[2]-center[2]; + dir = rot*dir; + dir[0] += center[0]+tx; + dir[1] += center[1]+ty; + dir[2] += center[2]+tz; + vtkIdType id = vtkNewPoints->InsertNextPoint(dir.data_block()); + container->GetPointIds()->InsertNextId(id); + } + vtkNewCells->InsertNextCell(container); + } + + m_FiberPolyData = vtkSmartPointer::New(); + m_FiberPolyData->SetPoints(vtkNewPoints); + m_FiberPolyData->SetLines(vtkNewCells); + UpdateColorCoding(); + UpdateFiberGeometry(); +} + void mitk::FiberBundleX::RotateAroundAxis(double x, double y, double z) { x = x*M_PI/180; y = y*M_PI/180; z = z*M_PI/180; vnl_matrix_fixed< double, 3, 3 > rotX; rotX.set_identity(); rotX[1][1] = cos(x); rotX[2][2] = rotX[1][1]; rotX[1][2] = -sin(x); rotX[2][1] = -rotX[1][2]; vnl_matrix_fixed< double, 3, 3 > rotY; rotY.set_identity(); rotY[0][0] = cos(y); rotY[2][2] = rotY[0][0]; rotY[0][2] = sin(y); rotY[2][0] = -rotY[0][2]; vnl_matrix_fixed< double, 3, 3 > rotZ; rotZ.set_identity(); rotZ[0][0] = cos(z); rotZ[1][1] = rotZ[0][0]; rotZ[0][1] = -sin(z); rotZ[1][0] = -rotZ[0][1]; mitk::Geometry3D::Pointer geom = this->GetGeometry(); mitk::Point3D center = geom->GetCenter(); vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); for (int i=0; iGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; 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) { 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) { MITK_INFO << "Removing short fibers"; if (lengthInMM<=0 || lengthInMMm_MaxFiberLength) // can't remove all fibers { MITK_WARN << "Process aborted. No fibers would be left!"; return false; } vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); float min = m_MaxFiberLength; boost::progress_display disp(m_NumFibers); for (int i=0; iGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); 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(float pointDistance, double tension, double continuity, double bias ) { if (pointDistance<=0) return; 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 << "Smoothing 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); int sampling = std::ceil(length/pointDistance); vtkSmartPointer xSpline = vtkSmartPointer::New(); vtkSmartPointer ySpline = vtkSmartPointer::New(); vtkSmartPointer zSpline = vtkSmartPointer::New(); xSpline->SetDefaultBias(bias); xSpline->SetDefaultTension(tension); xSpline->SetDefaultContinuity(continuity); ySpline->SetDefaultBias(bias); ySpline->SetDefaultTension(tension); ySpline->SetDefaultContinuity(continuity); zSpline->SetDefaultBias(bias); zSpline->SetDefaultTension(tension); zSpline->SetDefaultContinuity(continuity); vtkSmartPointer 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 = 10/pointDistance; } void mitk::FiberBundleX::DoFiberSmoothing(float pointDistance) { DoFiberSmoothing(pointDistance, 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/IODataStructures/FiberBundleX/mitkFiberBundleX.h b/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundleX/mitkFiberBundleX.h index 3b4e3479eb..31361727a3 100644 --- a/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundleX/mitkFiberBundleX.h +++ b/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundleX/mitkFiberBundleX.h @@ -1,160 +1,164 @@ /*=================================================================== 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_FiberBundleX_H #define _MITK_FiberBundleX_H //includes for MITK datastructure #include #include "FiberTrackingExports.h" #include //includes storing fiberdata -#include //may be replaced by class precompile argument -#include // may be replaced by class -#include // my be replaced by class +#include +#include +#include #include +#include //#include #include namespace mitk { /** * \brief Base Class for Fiber Bundles; */ class FiberTracking_EXPORT FiberBundleX : public BaseData { public: typedef itk::Image ItkUcharImgType; // fiber colorcodings static const char* COLORCODING_ORIENTATION_BASED; static const char* COLORCODING_FA_BASED; static const char* COLORCODING_CUSTOM; static const char* FIBER_ID_ARRAY; virtual void UpdateOutputInformation(); virtual void SetRequestedRegionToLargestPossibleRegion(); virtual bool RequestedRegionIsOutsideOfTheBufferedRegion(); virtual bool VerifyRequestedRegion(); virtual void SetRequestedRegion(const itk::DataObject *data ); mitkClassMacro( FiberBundleX, BaseData ) itkNewMacro( Self ) mitkNewMacro1Param(Self, vtkSmartPointer) // custom constructor // colorcoding related methods void SetColorCoding(const char*); void SetFAMap(mitk::Image::Pointer); template void SetFAMap(const mitk::PixelType pixelType, mitk::Image::Pointer); void DoColorCodingOrientationBased(); void DoColorCodingFaBased(); void DoUseFaFiberOpacity(); void ResetFiberOpacity(); // fiber smoothing/resampling void ResampleFibers(float pointDistance = 1); void DoFiberSmoothing(float pointDistance); void DoFiberSmoothing(float pointDistance, double tension, double continuity, double bias ); bool RemoveShortFibers(float lengthInMM); bool RemoveLongFibers(float lengthInMM); bool ApplyCurvatureThreshold(float minRadius, bool deleteFibers); void MirrorFibers(unsigned int axis); void RotateAroundAxis(double x, double y, double z); void TranslateFibers(double x, double y, double z); void ScaleFibers(double x, double y, double z); + void TransformFibers(double rx, double ry, double rz, double tx, double ty, double tz); + itk::Point TransformPoint(vnl_vector_fixed< double, 3 > point, double rx, double ry, double rz, double tx, double ty, double tz); + itk::Matrix< double, 3, 3 > TransformMatrix(itk::Matrix< double, 3, 3 > m, double rx, double ry, double rz); // add/subtract fibers FiberBundleX::Pointer AddBundle(FiberBundleX* fib); FiberBundleX::Pointer SubtractBundle(FiberBundleX* fib); // fiber subset extraction FiberBundleX::Pointer ExtractFiberSubset(PlanarFigure *pf); std::vector ExtractFiberIdSubset(PlanarFigure* pf); FiberBundleX::Pointer ExtractFiberSubset(ItkUcharImgType* mask, bool anyPoint); FiberBundleX::Pointer RemoveFibersOutside(ItkUcharImgType* mask, bool invert=false); vtkSmartPointer GeneratePolyDataByIds( std::vector ); // TODO: make protected void GenerateFiberIds(); // TODO: make protected // get/set data void SetFiberPolyData(vtkSmartPointer, bool updateGeometry = true); vtkSmartPointer GetFiberPolyData(); std::vector< std::string > GetAvailableColorCodings(); char* GetCurrentColorCoding(); itkGetMacro( NumFibers, int) itkGetMacro( FiberSampling, int) itkGetMacro( MinFiberLength, float ) itkGetMacro( MaxFiberLength, float ) itkGetMacro( MeanFiberLength, float ) itkGetMacro( MedianFiberLength, float ) itkGetMacro( LengthStDev, float ) std::vector GetPointsRoi() { return m_PointsRoi; } // copy fiber bundle mitk::FiberBundleX::Pointer GetDeepCopy(); // compare fiber bundles bool Equals(FiberBundleX* fib); protected: FiberBundleX( vtkPolyData* fiberPolyData = NULL ); virtual ~FiberBundleX(); itk::Point GetItkPoint(double point[3]); // calculate geometry from fiber extent void UpdateFiberGeometry(); // calculate colorcoding values according to m_CurrentColorCoding void UpdateColorCoding(); private: // actual fiber container vtkSmartPointer m_FiberPolyData; // contains fiber ids vtkSmartPointer m_FiberIdDataSet; char* m_CurrentColorCoding; int m_NumFibers; std::vector< float > m_FiberLengths; float m_MinFiberLength; float m_MaxFiberLength; float m_MeanFiberLength; float m_MedianFiberLength; float m_LengthStDev; int m_FiberSampling; std::vector m_PointsRoi; // this global variable needs to be refactored }; } // namespace mitk #endif /* _MITK_FiberBundleX_H */ diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberfoxView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberfoxView.cpp index d13af15074..491b134f8a 100755 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberfoxView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberfoxView.cpp @@ -1,2054 +1,2055 @@ /*=================================================================== 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 "usModuleRegistry.h" #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_NoiseFrame->setVisible(false); m_Controls->m_GhostFrame->setVisible(false); m_Controls->m_DistortionsFrame->setVisible(false); m_Controls->m_EddyFrame->setVisible(false); m_Controls->m_SpikeFrame->setVisible(false); m_Controls->m_AliasingFrame->setVisible(false); m_Controls->m_MotionArtifactFrame->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(double)), (QObject*) this, SLOT(OnFiberSamplingChanged(double))); 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_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_AddSpikes, SIGNAL(stateChanged(int)), (QObject*) this, SLOT(OnAddSpikes(int))); connect((QObject*) m_Controls->m_AddAliasing, SIGNAL(stateChanged(int)), (QObject*) this, SLOT(OnAddAliasing(int))); connect((QObject*) m_Controls->m_AddMotion, SIGNAL(stateChanged(int)), (QObject*) this, SLOT(OnAddMotion(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))); connect((QObject*) m_Controls->m_SaveParametersButton, SIGNAL(clicked()), (QObject*) this, SLOT(SaveParameters())); connect((QObject*) m_Controls->m_LoadParametersButton, SIGNAL(clicked()), (QObject*) this, SLOT(LoadParameters())); } } void QmitkFiberfoxView::UpdateImageParameters() { m_ImageGenParameters.artifactList.clear(); m_ImageGenParameters.nonFiberModelList.clear(); m_ImageGenParameters.fiberModelList.clear(); m_ImageGenParameters.signalModelString = ""; m_ImageGenParameters.artifactModelString = ""; m_ImageGenParameters.resultNode = mitk::DataNode::New(); m_ImageGenParameters.frequencyMap = NULL; m_ImageGenParameters.gradientDirections.clear(); m_ImageGenParameters.spikes = 0; m_ImageGenParameters.spikeAmplitude = 1; m_ImageGenParameters.wrap = 1; if (m_SelectedDWI.IsNotNull()) // use parameters of selected DWI { mitk::DiffusionImage::Pointer dwi = dynamic_cast*>(m_SelectedDWI->GetData()); m_ImageGenParameters.imageRegion = dwi->GetVectorImage()->GetLargestPossibleRegion(); m_ImageGenParameters.imageSpacing = dwi->GetVectorImage()->GetSpacing(); m_ImageGenParameters.imageOrigin = dwi->GetVectorImage()->GetOrigin(); m_ImageGenParameters.imageDirection = dwi->GetVectorImage()->GetDirection(); m_ImageGenParameters.b_value = dwi->GetB_Value(); mitk::DiffusionImage::GradientDirectionContainerType::Pointer dirs = dwi->GetDirections(); m_ImageGenParameters.numGradients = 0; 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]; m_ImageGenParameters.gradientDirections.push_back(g); if (dirs->at(i).magnitude()>0.0001) m_ImageGenParameters.numGradients++; } } else if (m_SelectedImage.IsNotNull()) // use geometry of selected image { mitk::Image::Pointer img = dynamic_cast(m_SelectedImage->GetData()); itk::Image< float, 3 >::Pointer itkImg = itk::Image< float, 3 >::New(); CastToItkImage< itk::Image< float, 3 > >(img, itkImg); m_ImageGenParameters.imageRegion = itkImg->GetLargestPossibleRegion(); m_ImageGenParameters.imageSpacing = itkImg->GetSpacing(); m_ImageGenParameters.imageOrigin = itkImg->GetOrigin(); m_ImageGenParameters.imageDirection = itkImg->GetDirection(); m_ImageGenParameters.numGradients = m_Controls->m_NumGradientsBox->value(); m_ImageGenParameters.gradientDirections = GenerateHalfShell(m_Controls->m_NumGradientsBox->value()); m_ImageGenParameters.b_value = m_Controls->m_BvalueBox->value(); } else // use GUI parameters { m_ImageGenParameters.imageRegion.SetSize(0, m_Controls->m_SizeX->value()); m_ImageGenParameters.imageRegion.SetSize(1, m_Controls->m_SizeY->value()); m_ImageGenParameters.imageRegion.SetSize(2, m_Controls->m_SizeZ->value()); m_ImageGenParameters.imageSpacing[0] = m_Controls->m_SpacingX->value(); m_ImageGenParameters.imageSpacing[1] = m_Controls->m_SpacingY->value(); m_ImageGenParameters.imageSpacing[2] = m_Controls->m_SpacingZ->value(); m_ImageGenParameters.imageOrigin[0] = m_ImageGenParameters.imageSpacing[0]/2; m_ImageGenParameters.imageOrigin[1] = m_ImageGenParameters.imageSpacing[1]/2; m_ImageGenParameters.imageOrigin[2] = m_ImageGenParameters.imageSpacing[2]/2; m_ImageGenParameters.imageDirection.SetIdentity(); m_ImageGenParameters.numGradients = m_Controls->m_NumGradientsBox->value(); m_ImageGenParameters.gradientDirections = GenerateHalfShell(m_Controls->m_NumGradientsBox->value());; m_ImageGenParameters.b_value = m_Controls->m_BvalueBox->value(); } // signal relaxation m_ImageGenParameters.doSimulateRelaxation = m_Controls->m_RelaxationBox->isChecked(); if (m_ImageGenParameters.doSimulateRelaxation && m_SelectedBundles.size()>0 ) m_ImageGenParameters.artifactModelString += "_RELAX"; // N/2 ghosts if (m_Controls->m_AddGhosts->isChecked()) { m_ImageGenParameters.artifactModelString += "_GHOST"; m_ImageGenParameters.kspaceLineOffset = m_Controls->m_kOffsetBox->value(); } else m_ImageGenParameters.kspaceLineOffset = 0; // Aliasing if (m_Controls->m_AddAliasing->isChecked()) { m_ImageGenParameters.artifactModelString += "_ALIASING"; m_ImageGenParameters.wrap = 1/m_Controls->m_WrapBox->value(); } // Motion m_ImageGenParameters.doAddMotion = m_Controls->m_AddMotion->isChecked(); m_ImageGenParameters.randomMotion = m_Controls->m_RandomMotion->isChecked(); m_ImageGenParameters.translation[0] = m_Controls->m_MaxTranslationBoxX->value(); m_ImageGenParameters.translation[1] = m_Controls->m_MaxTranslationBoxY->value(); m_ImageGenParameters.translation[2] = m_Controls->m_MaxTranslationBoxZ->value(); m_ImageGenParameters.rotation[0] = m_Controls->m_MaxRotationBoxX->value(); m_ImageGenParameters.rotation[1] = m_Controls->m_MaxRotationBoxY->value(); m_ImageGenParameters.rotation[2] = m_Controls->m_MaxRotationBoxZ->value(); if ( m_Controls->m_AddMotion->isChecked() && m_SelectedBundles.size()>0 ) m_ImageGenParameters.artifactModelString += "_MOTION"; // other imaging parameters m_ImageGenParameters.tLine = m_Controls->m_LineReadoutTimeBox->value(); m_ImageGenParameters.tInhom = m_Controls->m_T2starBox->value(); m_ImageGenParameters.tEcho = m_Controls->m_TEbox->value(); m_ImageGenParameters.repetitions = m_Controls->m_RepetitionsBox->value(); m_ImageGenParameters.doDisablePartialVolume = m_Controls->m_EnforcePureFiberVoxelsBox->isChecked(); m_ImageGenParameters.interpolationShrink = m_Controls->m_InterpolationShrink->value(); m_ImageGenParameters.axonRadius = m_Controls->m_FiberRadius->value(); m_ImageGenParameters.signalScale = m_Controls->m_SignalScaleBox->value(); if (m_Controls->m_AddSpikes->isChecked()) { m_ImageGenParameters.spikes = m_Controls->m_SpikeNumBox->value(); m_ImageGenParameters.spikeAmplitude = m_Controls->m_SpikeScaleBox->value(); m_ImageGenParameters.artifactModelString += "_SPIKES"; } // adjust echo time if needed if ( m_ImageGenParameters.tEcho < m_ImageGenParameters.imageRegion.GetSize(1)*m_ImageGenParameters.tLine ) { this->m_Controls->m_TEbox->setValue( m_ImageGenParameters.imageRegion.GetSize(1)*m_ImageGenParameters.tLine ); m_ImageGenParameters.tEcho = m_Controls->m_TEbox->value(); QMessageBox::information( NULL, "Warning", "Echo time is too short! Time not sufficient to read slice. Automaticall adjusted to "+QString::number(m_ImageGenParameters.tEcho)+" ms"); } // rician noise if (m_Controls->m_AddNoise->isChecked()) { double noiseVariance = m_Controls->m_NoiseLevel->value(); m_ImageGenParameters.ricianNoiseModel.SetNoiseVariance(noiseVariance); m_ImageGenParameters.artifactModelString += "_NOISE"; m_ImageGenParameters.artifactModelString += QString::number(noiseVariance); m_ImageGenParameters.resultNode->AddProperty("Fiberfox.Noise-Variance", DoubleProperty::New(noiseVariance)); } else m_ImageGenParameters.ricianNoiseModel.SetNoiseVariance(0); // gibbs ringing m_ImageGenParameters.addGibbsRinging = m_Controls->m_AddGibbsRinging->isChecked(); if (m_Controls->m_AddGibbsRinging->isChecked()) m_ImageGenParameters.artifactModelString += "_RINGING"; // adjusting line readout time to the adapted image size needed for the DFT int y = m_ImageGenParameters.imageRegion.GetSize(1); if ( y%2 == 1 ) y += 1; if ( y>m_ImageGenParameters.imageRegion.GetSize(1) ) m_ImageGenParameters.tLine *= (double)m_ImageGenParameters.imageRegion.GetSize(1)/y; // 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 (m_ImageGenParameters.imageRegion.GetSize(0)==itkImg->GetLargestPossibleRegion().GetSize(0) && m_ImageGenParameters.imageRegion.GetSize(1)==itkImg->GetLargestPossibleRegion().GetSize(1) && m_ImageGenParameters.imageRegion.GetSize(2)==itkImg->GetLargestPossibleRegion().GetSize(2)) { m_ImageGenParameters.frequencyMap = itkImg; m_ImageGenParameters.artifactModelString += "_DISTORTED"; } } m_ImageGenParameters.doSimulateEddyCurrents = m_Controls->m_AddEddy->isChecked(); m_ImageGenParameters.eddyStrength = 0; if (m_Controls->m_AddEddy->isChecked()) { m_ImageGenParameters.eddyStrength = m_Controls->m_EddyGradientStrength->value(); m_ImageGenParameters.artifactModelString += "_EDDY"; } // signal models m_ImageGenParameters.comp3Weight = 1; m_ImageGenParameters.comp4Weight = 0; if (m_Controls->m_Compartment4Box->currentIndex()>0) { m_ImageGenParameters.comp4Weight = m_Controls->m_Comp4FractionBox->value(); m_ImageGenParameters.comp3Weight -= m_ImageGenParameters.comp4Weight; } // compartment 1 switch (m_Controls->m_Compartment1Box->currentIndex()) { case 0: m_StickModel1.SetGradientList(m_ImageGenParameters.gradientDirections); m_StickModel1.SetBvalue(m_ImageGenParameters.b_value); m_StickModel1.SetDiffusivity(m_Controls->m_StickWidget1->GetD()); m_StickModel1.SetT2(m_Controls->m_StickWidget1->GetT2()); m_ImageGenParameters.fiberModelList.push_back(&m_StickModel1); m_ImageGenParameters.signalModelString += "Stick"; m_ImageGenParameters.resultNode->AddProperty("Fiberfox.Compartment1.Description", StringProperty::New("Intra-axonal compartment") ); m_ImageGenParameters.resultNode->AddProperty("Fiberfox.Compartment1.Model", StringProperty::New("Stick") ); m_ImageGenParameters.resultNode->AddProperty("Fiberfox.Compartment1.D", DoubleProperty::New(m_Controls->m_StickWidget1->GetD()) ); m_ImageGenParameters.resultNode->AddProperty("Fiberfox.Compartment1.T2", DoubleProperty::New(m_StickModel1.GetT2()) ); break; case 1: m_ZeppelinModel1.SetGradientList(m_ImageGenParameters.gradientDirections); m_ZeppelinModel1.SetBvalue(m_ImageGenParameters.b_value); m_ZeppelinModel1.SetDiffusivity1(m_Controls->m_ZeppelinWidget1->GetD1()); m_ZeppelinModel1.SetDiffusivity2(m_Controls->m_ZeppelinWidget1->GetD2()); m_ZeppelinModel1.SetDiffusivity3(m_Controls->m_ZeppelinWidget1->GetD2()); m_ZeppelinModel1.SetT2(m_Controls->m_ZeppelinWidget1->GetT2()); m_ImageGenParameters.fiberModelList.push_back(&m_ZeppelinModel1); m_ImageGenParameters.signalModelString += "Zeppelin"; m_ImageGenParameters.resultNode->AddProperty("Fiberfox.Compartment1.Description", StringProperty::New("Intra-axonal compartment") ); m_ImageGenParameters.resultNode->AddProperty("Fiberfox.Compartment1.Model", StringProperty::New("Zeppelin") ); m_ImageGenParameters.resultNode->AddProperty("Fiberfox.Compartment1.D1", DoubleProperty::New(m_Controls->m_ZeppelinWidget1->GetD1()) ); m_ImageGenParameters.resultNode->AddProperty("Fiberfox.Compartment1.D2", DoubleProperty::New(m_Controls->m_ZeppelinWidget1->GetD2()) ); m_ImageGenParameters.resultNode->AddProperty("Fiberfox.Compartment1.T2", DoubleProperty::New(m_ZeppelinModel1.GetT2()) ); break; case 2: m_TensorModel1.SetGradientList(m_ImageGenParameters.gradientDirections); m_TensorModel1.SetBvalue(m_ImageGenParameters.b_value); m_TensorModel1.SetDiffusivity1(m_Controls->m_TensorWidget1->GetD1()); m_TensorModel1.SetDiffusivity2(m_Controls->m_TensorWidget1->GetD2()); m_TensorModel1.SetDiffusivity3(m_Controls->m_TensorWidget1->GetD3()); m_TensorModel1.SetT2(m_Controls->m_TensorWidget1->GetT2()); m_ImageGenParameters.fiberModelList.push_back(&m_TensorModel1); m_ImageGenParameters.signalModelString += "Tensor"; m_ImageGenParameters.resultNode->AddProperty("Fiberfox.Compartment1.Description", StringProperty::New("Intra-axonal compartment") ); m_ImageGenParameters.resultNode->AddProperty("Fiberfox.Compartment1.Model", StringProperty::New("Tensor") ); m_ImageGenParameters.resultNode->AddProperty("Fiberfox.Compartment1.D1", DoubleProperty::New(m_Controls->m_TensorWidget1->GetD1()) ); m_ImageGenParameters.resultNode->AddProperty("Fiberfox.Compartment1.D2", DoubleProperty::New(m_Controls->m_TensorWidget1->GetD2()) ); m_ImageGenParameters.resultNode->AddProperty("Fiberfox.Compartment1.D3", DoubleProperty::New(m_Controls->m_TensorWidget1->GetD3()) ); m_ImageGenParameters.resultNode->AddProperty("Fiberfox.Compartment1.T2", DoubleProperty::New(m_ZeppelinModel1.GetT2()) ); break; } // compartment 2 switch (m_Controls->m_Compartment2Box->currentIndex()) { case 0: break; case 1: m_StickModel2.SetGradientList(m_ImageGenParameters.gradientDirections); m_StickModel2.SetBvalue(m_ImageGenParameters.b_value); m_StickModel2.SetDiffusivity(m_Controls->m_StickWidget2->GetD()); m_StickModel2.SetT2(m_Controls->m_StickWidget2->GetT2()); m_ImageGenParameters.fiberModelList.push_back(&m_StickModel2); m_ImageGenParameters.signalModelString += "Stick"; m_ImageGenParameters.resultNode->AddProperty("Fiberfox.Compartment2.Description", StringProperty::New("Inter-axonal compartment") ); m_ImageGenParameters.resultNode->AddProperty("Fiberfox.Compartment2.Model", StringProperty::New("Stick") ); m_ImageGenParameters.resultNode->AddProperty("Fiberfox.Compartment2.D", DoubleProperty::New(m_Controls->m_StickWidget2->GetD()) ); m_ImageGenParameters.resultNode->AddProperty("Fiberfox.Compartment2.T2", DoubleProperty::New(m_StickModel2.GetT2()) ); break; case 2: m_ZeppelinModel2.SetGradientList(m_ImageGenParameters.gradientDirections); m_ZeppelinModel2.SetBvalue(m_ImageGenParameters.b_value); m_ZeppelinModel2.SetDiffusivity1(m_Controls->m_ZeppelinWidget2->GetD1()); m_ZeppelinModel2.SetDiffusivity2(m_Controls->m_ZeppelinWidget2->GetD2()); m_ZeppelinModel2.SetDiffusivity3(m_Controls->m_ZeppelinWidget2->GetD2()); m_ZeppelinModel2.SetT2(m_Controls->m_ZeppelinWidget2->GetT2()); m_ImageGenParameters.fiberModelList.push_back(&m_ZeppelinModel2); m_ImageGenParameters.signalModelString += "Zeppelin"; m_ImageGenParameters.resultNode->AddProperty("Fiberfox.Compartment2.Description", StringProperty::New("Inter-axonal compartment") ); m_ImageGenParameters.resultNode->AddProperty("Fiberfox.Compartment2.Model", StringProperty::New("Zeppelin") ); m_ImageGenParameters.resultNode->AddProperty("Fiberfox.Compartment2.D1", DoubleProperty::New(m_Controls->m_ZeppelinWidget2->GetD1()) ); m_ImageGenParameters.resultNode->AddProperty("Fiberfox.Compartment2.D2", DoubleProperty::New(m_Controls->m_ZeppelinWidget2->GetD2()) ); m_ImageGenParameters.resultNode->AddProperty("Fiberfox.Compartment2.T2", DoubleProperty::New(m_ZeppelinModel2.GetT2()) ); break; case 3: m_TensorModel2.SetGradientList(m_ImageGenParameters.gradientDirections); m_TensorModel2.SetBvalue(m_ImageGenParameters.b_value); m_TensorModel2.SetDiffusivity1(m_Controls->m_TensorWidget2->GetD1()); m_TensorModel2.SetDiffusivity2(m_Controls->m_TensorWidget2->GetD2()); m_TensorModel2.SetDiffusivity3(m_Controls->m_TensorWidget2->GetD3()); m_TensorModel2.SetT2(m_Controls->m_TensorWidget2->GetT2()); m_ImageGenParameters.fiberModelList.push_back(&m_TensorModel2); m_ImageGenParameters.signalModelString += "Tensor"; m_ImageGenParameters.resultNode->AddProperty("Fiberfox.Compartment2.Description", StringProperty::New("Inter-axonal compartment") ); m_ImageGenParameters.resultNode->AddProperty("Fiberfox.Compartment2.Model", StringProperty::New("Tensor") ); m_ImageGenParameters.resultNode->AddProperty("Fiberfox.Compartment2.D1", DoubleProperty::New(m_Controls->m_TensorWidget2->GetD1()) ); m_ImageGenParameters.resultNode->AddProperty("Fiberfox.Compartment2.D2", DoubleProperty::New(m_Controls->m_TensorWidget2->GetD2()) ); m_ImageGenParameters.resultNode->AddProperty("Fiberfox.Compartment2.D3", DoubleProperty::New(m_Controls->m_TensorWidget2->GetD3()) ); m_ImageGenParameters.resultNode->AddProperty("Fiberfox.Compartment2.T2", DoubleProperty::New(m_ZeppelinModel2.GetT2()) ); break; } // compartment 3 switch (m_Controls->m_Compartment3Box->currentIndex()) { case 0: m_BallModel1.SetGradientList(m_ImageGenParameters.gradientDirections); m_BallModel1.SetBvalue(m_ImageGenParameters.b_value); m_BallModel1.SetDiffusivity(m_Controls->m_BallWidget1->GetD()); m_BallModel1.SetT2(m_Controls->m_BallWidget1->GetT2()); m_BallModel1.SetWeight(m_ImageGenParameters.comp3Weight); m_ImageGenParameters.nonFiberModelList.push_back(&m_BallModel1); m_ImageGenParameters.signalModelString += "Ball"; m_ImageGenParameters.resultNode->AddProperty("Fiberfox.Compartment3.Description", StringProperty::New("Extra-axonal compartment 1") ); m_ImageGenParameters.resultNode->AddProperty("Fiberfox.Compartment3.Model", StringProperty::New("Ball") ); m_ImageGenParameters.resultNode->AddProperty("Fiberfox.Compartment3.D", DoubleProperty::New(m_Controls->m_BallWidget1->GetD()) ); m_ImageGenParameters.resultNode->AddProperty("Fiberfox.Compartment3.T2", DoubleProperty::New(m_BallModel1.GetT2()) ); break; case 1: m_AstrosticksModel1.SetGradientList(m_ImageGenParameters.gradientDirections); m_AstrosticksModel1.SetBvalue(m_ImageGenParameters.b_value); m_AstrosticksModel1.SetDiffusivity(m_Controls->m_AstrosticksWidget1->GetD()); m_AstrosticksModel1.SetT2(m_Controls->m_AstrosticksWidget1->GetT2()); m_AstrosticksModel1.SetRandomizeSticks(m_Controls->m_AstrosticksWidget1->GetRandomizeSticks()); m_AstrosticksModel1.SetWeight(m_ImageGenParameters.comp3Weight); m_ImageGenParameters.nonFiberModelList.push_back(&m_AstrosticksModel1); m_ImageGenParameters.signalModelString += "Astrosticks"; m_ImageGenParameters.resultNode->AddProperty("Fiberfox.Compartment3.Description", StringProperty::New("Extra-axonal compartment 1") ); m_ImageGenParameters.resultNode->AddProperty("Fiberfox.Compartment3.Model", StringProperty::New("Astrosticks") ); m_ImageGenParameters.resultNode->AddProperty("Fiberfox.Compartment3.D", DoubleProperty::New(m_Controls->m_AstrosticksWidget1->GetD()) ); m_ImageGenParameters.resultNode->AddProperty("Fiberfox.Compartment3.T2", DoubleProperty::New(m_AstrosticksModel1.GetT2()) ); m_ImageGenParameters.resultNode->AddProperty("Fiberfox.Compartment3.RandomSticks", BoolProperty::New(m_Controls->m_AstrosticksWidget1->GetRandomizeSticks()) ); break; case 2: m_DotModel1.SetGradientList(m_ImageGenParameters.gradientDirections); m_DotModel1.SetT2(m_Controls->m_DotWidget1->GetT2()); m_DotModel1.SetWeight(m_ImageGenParameters.comp3Weight); m_ImageGenParameters.nonFiberModelList.push_back(&m_DotModel1); m_ImageGenParameters.signalModelString += "Dot"; m_ImageGenParameters.resultNode->AddProperty("Fiberfox.Compartment3.Description", StringProperty::New("Extra-axonal compartment 1") ); m_ImageGenParameters.resultNode->AddProperty("Fiberfox.Compartment3.Model", StringProperty::New("Dot") ); m_ImageGenParameters.resultNode->AddProperty("Fiberfox.Compartment3.T2", DoubleProperty::New(m_DotModel1.GetT2()) ); break; } // compartment 4 switch (m_Controls->m_Compartment4Box->currentIndex()) { case 0: break; case 1: m_BallModel2.SetGradientList(m_ImageGenParameters.gradientDirections); m_BallModel2.SetBvalue(m_ImageGenParameters.b_value); m_BallModel2.SetDiffusivity(m_Controls->m_BallWidget2->GetD()); m_BallModel2.SetT2(m_Controls->m_BallWidget2->GetT2()); m_BallModel2.SetWeight(m_ImageGenParameters.comp4Weight); m_ImageGenParameters.nonFiberModelList.push_back(&m_BallModel2); m_ImageGenParameters.signalModelString += "Ball"; m_ImageGenParameters.resultNode->AddProperty("Fiberfox.Compartment4.Description", StringProperty::New("Extra-axonal compartment 2") ); m_ImageGenParameters.resultNode->AddProperty("Fiberfox.Compartment4.Model", StringProperty::New("Ball") ); m_ImageGenParameters.resultNode->AddProperty("Fiberfox.Compartment4.D", DoubleProperty::New(m_Controls->m_BallWidget2->GetD()) ); m_ImageGenParameters.resultNode->AddProperty("Fiberfox.Compartment4.T2", DoubleProperty::New(m_BallModel2.GetT2()) ); break; case 2: m_AstrosticksModel2.SetGradientList(m_ImageGenParameters.gradientDirections); m_AstrosticksModel2.SetBvalue(m_ImageGenParameters.b_value); m_AstrosticksModel2.SetDiffusivity(m_Controls->m_AstrosticksWidget2->GetD()); m_AstrosticksModel2.SetT2(m_Controls->m_AstrosticksWidget2->GetT2()); m_AstrosticksModel2.SetRandomizeSticks(m_Controls->m_AstrosticksWidget2->GetRandomizeSticks()); m_AstrosticksModel2.SetWeight(m_ImageGenParameters.comp4Weight); m_ImageGenParameters.nonFiberModelList.push_back(&m_AstrosticksModel2); m_ImageGenParameters.signalModelString += "Astrosticks"; m_ImageGenParameters.resultNode->AddProperty("Fiberfox.Compartment4.Description", StringProperty::New("Extra-axonal compartment 2") ); m_ImageGenParameters.resultNode->AddProperty("Fiberfox.Compartment4.Model", StringProperty::New("Astrosticks") ); m_ImageGenParameters.resultNode->AddProperty("Fiberfox.Compartment4.D", DoubleProperty::New(m_Controls->m_AstrosticksWidget2->GetD()) ); m_ImageGenParameters.resultNode->AddProperty("Fiberfox.Compartment4.T2", DoubleProperty::New(m_AstrosticksModel2.GetT2()) ); m_ImageGenParameters.resultNode->AddProperty("Fiberfox.Compartment4.RandomSticks", BoolProperty::New(m_Controls->m_AstrosticksWidget2->GetRandomizeSticks()) ); break; case 3: m_DotModel2.SetGradientList(m_ImageGenParameters.gradientDirections); m_DotModel2.SetT2(m_Controls->m_DotWidget2->GetT2()); m_DotModel2.SetWeight(m_ImageGenParameters.comp4Weight); m_ImageGenParameters.nonFiberModelList.push_back(&m_DotModel2); m_ImageGenParameters.signalModelString += "Dot"; m_ImageGenParameters.resultNode->AddProperty("Fiberfox.Compartment4.Description", StringProperty::New("Extra-axonal compartment 2") ); m_ImageGenParameters.resultNode->AddProperty("Fiberfox.Compartment4.Model", StringProperty::New("Dot") ); m_ImageGenParameters.resultNode->AddProperty("Fiberfox.Compartment4.T2", DoubleProperty::New(m_DotModel2.GetT2()) ); break; } m_ImageGenParameters.resultNode->AddProperty("Fiberfox.InterpolationShrink", IntProperty::New(m_ImageGenParameters.interpolationShrink)); m_ImageGenParameters.resultNode->AddProperty("Fiberfox.SignalScale", IntProperty::New(m_ImageGenParameters.signalScale)); m_ImageGenParameters.resultNode->AddProperty("Fiberfox.FiberRadius", IntProperty::New(m_ImageGenParameters.axonRadius)); m_ImageGenParameters.resultNode->AddProperty("Fiberfox.Tinhom", IntProperty::New(m_ImageGenParameters.tInhom)); m_ImageGenParameters.resultNode->AddProperty("Fiberfox.Repetitions", IntProperty::New(m_ImageGenParameters.repetitions)); m_ImageGenParameters.resultNode->AddProperty("Fiberfox.b-value", DoubleProperty::New(m_ImageGenParameters.b_value)); m_ImageGenParameters.resultNode->AddProperty("Fiberfox.Model", StringProperty::New(m_ImageGenParameters.signalModelString.toStdString())); m_ImageGenParameters.resultNode->AddProperty("Fiberfox.PureFiberVoxels", BoolProperty::New(m_ImageGenParameters.doDisablePartialVolume)); m_ImageGenParameters.resultNode->AddProperty("binary", BoolProperty::New(false)); } void QmitkFiberfoxView::SaveParameters() { UpdateImageParameters(); QString filename = QFileDialog::getSaveFileName( 0, tr("Save Parameters"), QDir::currentPath()+"/param.ffp", tr("Fiberfox Parameters (*.ffp)") ); if(filename.isEmpty() || filename.isNull()) return; if(!filename.endsWith(".ffp")) filename += ".ffp"; boost::property_tree::ptree parameters; // fiber generation parameters parameters.put("fiberfox.fibers.realtime", m_Controls->m_RealTimeFibers->isChecked()); parameters.put("fiberfox.fibers.showadvanced", m_Controls->m_AdvancedOptionsBox->isChecked()); parameters.put("fiberfox.fibers.distribution", m_Controls->m_DistributionBox->currentIndex()); parameters.put("fiberfox.fibers.variance", m_Controls->m_VarianceBox->value()); parameters.put("fiberfox.fibers.density", m_Controls->m_FiberDensityBox->value()); parameters.put("fiberfox.fibers.spline.sampling", m_Controls->m_FiberSamplingBox->value()); parameters.put("fiberfox.fibers.spline.tension", m_Controls->m_TensionBox->value()); parameters.put("fiberfox.fibers.spline.continuity", m_Controls->m_ContinuityBox->value()); parameters.put("fiberfox.fibers.spline.bias", m_Controls->m_BiasBox->value()); parameters.put("fiberfox.fibers.constantradius", m_Controls->m_ConstantRadiusBox->isChecked()); parameters.put("fiberfox.fibers.rotation.x", m_Controls->m_XrotBox->value()); parameters.put("fiberfox.fibers.rotation.y", m_Controls->m_YrotBox->value()); parameters.put("fiberfox.fibers.rotation.z", m_Controls->m_ZrotBox->value()); parameters.put("fiberfox.fibers.translation.x", m_Controls->m_XtransBox->value()); parameters.put("fiberfox.fibers.translation.y", m_Controls->m_YtransBox->value()); parameters.put("fiberfox.fibers.translation.z", m_Controls->m_ZtransBox->value()); parameters.put("fiberfox.fibers.scale.x", m_Controls->m_XscaleBox->value()); parameters.put("fiberfox.fibers.scale.y", m_Controls->m_YscaleBox->value()); parameters.put("fiberfox.fibers.scale.z", m_Controls->m_ZscaleBox->value()); parameters.put("fiberfox.fibers.includeFiducials", m_Controls->m_IncludeFiducials->isChecked()); parameters.put("fiberfox.fibers.includeFiducials", m_Controls->m_IncludeFiducials->isChecked()); // image generation parameters parameters.put("fiberfox.image.basic.size.x", m_ImageGenParameters.imageRegion.GetSize(0)); parameters.put("fiberfox.image.basic.size.y", m_ImageGenParameters.imageRegion.GetSize(1)); parameters.put("fiberfox.image.basic.size.z", m_ImageGenParameters.imageRegion.GetSize(2)); parameters.put("fiberfox.image.basic.spacing.x", m_ImageGenParameters.imageSpacing[0]); parameters.put("fiberfox.image.basic.spacing.y", m_ImageGenParameters.imageSpacing[1]); parameters.put("fiberfox.image.basic.spacing.z", m_ImageGenParameters.imageSpacing[2]); parameters.put("fiberfox.image.basic.numgradients", m_ImageGenParameters.numGradients); parameters.put("fiberfox.image.basic.bvalue", m_ImageGenParameters.b_value); parameters.put("fiberfox.image.showadvanced", m_Controls->m_AdvancedOptionsBox_2->isChecked()); parameters.put("fiberfox.image.repetitions", m_ImageGenParameters.repetitions); parameters.put("fiberfox.image.signalScale", m_ImageGenParameters.signalScale); parameters.put("fiberfox.image.tEcho", m_ImageGenParameters.tEcho); parameters.put("fiberfox.image.tLine", m_Controls->m_LineReadoutTimeBox->value()); parameters.put("fiberfox.image.tInhom", m_ImageGenParameters.tInhom); parameters.put("fiberfox.image.axonRadius", m_ImageGenParameters.axonRadius); parameters.put("fiberfox.image.interpolationShrink", m_ImageGenParameters.interpolationShrink); parameters.put("fiberfox.image.doSimulateRelaxation", m_ImageGenParameters.doSimulateRelaxation); parameters.put("fiberfox.image.doDisablePartialVolume", m_ImageGenParameters.doDisablePartialVolume); parameters.put("fiberfox.image.outputvolumefractions", m_Controls->m_VolumeFractionsBox->isChecked()); parameters.put("fiberfox.image.artifacts.addnoise", m_Controls->m_AddNoise->isChecked()); parameters.put("fiberfox.image.artifacts.noisevariance", m_Controls->m_NoiseLevel->value()); parameters.put("fiberfox.image.artifacts.addghost", m_Controls->m_AddGhosts->isChecked()); parameters.put("fiberfox.image.artifacts.kspaceLineOffset", m_Controls->m_kOffsetBox->value()); parameters.put("fiberfox.image.artifacts.distortions", m_Controls->m_AddDistortions->isChecked()); parameters.put("fiberfox.image.artifacts.addeddy", m_Controls->m_AddEddy->isChecked()); parameters.put("fiberfox.image.artifacts.eddyStrength", m_Controls->m_EddyGradientStrength->value()); parameters.put("fiberfox.image.artifacts.addringing", m_Controls->m_AddGibbsRinging->isChecked()); parameters.put("fiberfox.image.artifacts.addspikes", m_Controls->m_AddSpikes->isChecked()); parameters.put("fiberfox.image.artifacts.spikesnum", m_Controls->m_SpikeNumBox->value()); parameters.put("fiberfox.image.artifacts.spikesscale", m_Controls->m_SpikeScaleBox->value()); parameters.put("fiberfox.image.artifacts.addaliasing", m_Controls->m_AddAliasing->isChecked()); parameters.put("fiberfox.image.artifacts.aliasingfactor", m_Controls->m_WrapBox->value()); parameters.put("fiberfox.image.compartment1.index", m_Controls->m_Compartment1Box->currentIndex()); parameters.put("fiberfox.image.compartment2.index", m_Controls->m_Compartment2Box->currentIndex()); parameters.put("fiberfox.image.compartment3.index", m_Controls->m_Compartment3Box->currentIndex()); parameters.put("fiberfox.image.compartment4.index", m_Controls->m_Compartment4Box->currentIndex()); parameters.put("fiberfox.image.compartment1.stick.d", m_Controls->m_StickWidget1->GetD()); parameters.put("fiberfox.image.compartment1.stick.t2", m_Controls->m_StickWidget1->GetT2()); parameters.put("fiberfox.image.compartment1.zeppelin.d1", m_Controls->m_ZeppelinWidget1->GetD1()); parameters.put("fiberfox.image.compartment1.zeppelin.d2", m_Controls->m_ZeppelinWidget1->GetD2()); parameters.put("fiberfox.image.compartment1.zeppelin.t2", m_Controls->m_ZeppelinWidget1->GetT2()); parameters.put("fiberfox.image.compartment1.tensor.d1", m_Controls->m_TensorWidget1->GetD1()); parameters.put("fiberfox.image.compartment1.tensor.d2", m_Controls->m_TensorWidget1->GetD2()); parameters.put("fiberfox.image.compartment1.tensor.d3", m_Controls->m_TensorWidget1->GetD3()); parameters.put("fiberfox.image.compartment1.tensor.t2", m_Controls->m_TensorWidget1->GetT2()); parameters.put("fiberfox.image.compartment2.stick.d", m_Controls->m_StickWidget2->GetD()); parameters.put("fiberfox.image.compartment2.stick.t2", m_Controls->m_StickWidget2->GetT2()); parameters.put("fiberfox.image.compartment2.zeppelin.d1", m_Controls->m_ZeppelinWidget2->GetD1()); parameters.put("fiberfox.image.compartment2.zeppelin.d2", m_Controls->m_ZeppelinWidget2->GetD2()); parameters.put("fiberfox.image.compartment2.zeppelin.t2", m_Controls->m_ZeppelinWidget2->GetT2()); parameters.put("fiberfox.image.compartment2.tensor.d1", m_Controls->m_TensorWidget2->GetD1()); parameters.put("fiberfox.image.compartment2.tensor.d2", m_Controls->m_TensorWidget2->GetD2()); parameters.put("fiberfox.image.compartment2.tensor.d3", m_Controls->m_TensorWidget2->GetD3()); parameters.put("fiberfox.image.compartment2.tensor.t2", m_Controls->m_TensorWidget2->GetT2()); parameters.put("fiberfox.image.compartment3.ball.d", m_Controls->m_BallWidget1->GetD()); parameters.put("fiberfox.image.compartment3.ball.t2", m_Controls->m_BallWidget1->GetT2()); parameters.put("fiberfox.image.compartment3.astrosticks.d", m_Controls->m_AstrosticksWidget1->GetD()); parameters.put("fiberfox.image.compartment3.astrosticks.t2", m_Controls->m_AstrosticksWidget1->GetT2()); parameters.put("fiberfox.image.compartment3.astrosticks.randomize", m_Controls->m_AstrosticksWidget1->GetRandomizeSticks()); parameters.put("fiberfox.image.compartment3.dot.t2", m_Controls->m_DotWidget1->GetT2()); parameters.put("fiberfox.image.compartment4.ball.d", m_Controls->m_BallWidget2->GetD()); parameters.put("fiberfox.image.compartment4.ball.t2", m_Controls->m_BallWidget2->GetT2()); parameters.put("fiberfox.image.compartment4.astrosticks.d", m_Controls->m_AstrosticksWidget2->GetD()); parameters.put("fiberfox.image.compartment4.astrosticks.t2", m_Controls->m_AstrosticksWidget2->GetT2()); parameters.put("fiberfox.image.compartment4.astrosticks.randomize", m_Controls->m_AstrosticksWidget2->GetRandomizeSticks()); parameters.put("fiberfox.image.compartment4.dot.t2", m_Controls->m_DotWidget2->GetT2()); parameters.put("fiberfox.image.compartment4.weight", m_Controls->m_Comp4FractionBox->value()); boost::property_tree::xml_parser::write_xml(filename.toStdString(), parameters); } void QmitkFiberfoxView::LoadParameters() { QString filename = QFileDialog::getOpenFileName(0, tr("Load Parameters"), QDir::currentPath(), tr("Fiberfox Parameters (*.ffp)") ); if(filename.isEmpty() || filename.isNull()) return; boost::property_tree::ptree parameters; boost::property_tree::xml_parser::read_xml(filename.toStdString(), parameters); BOOST_FOREACH( boost::property_tree::ptree::value_type const& v1, parameters.get_child("fiberfox") ) { if( v1.first == "fibers" ) { m_Controls->m_RealTimeFibers->setChecked(v1.second.get("realtime")); m_Controls->m_AdvancedOptionsBox->setChecked(v1.second.get("showadvanced")); m_Controls->m_DistributionBox->setCurrentIndex(v1.second.get("distribution")); m_Controls->m_VarianceBox->setValue(v1.second.get("variance")); m_Controls->m_FiberDensityBox->setValue(v1.second.get("density")); m_Controls->m_IncludeFiducials->setChecked(v1.second.get("includeFiducials")); m_Controls->m_ConstantRadiusBox->setChecked(v1.second.get("constantradius")); BOOST_FOREACH( boost::property_tree::ptree::value_type const& v2, v1.second ) { if( v2.first == "spline" ) { m_Controls->m_FiberSamplingBox->setValue(v2.second.get("sampling")); m_Controls->m_TensionBox->setValue(v2.second.get("tension")); m_Controls->m_ContinuityBox->setValue(v2.second.get("continuity")); m_Controls->m_BiasBox->setValue(v2.second.get("bias")); } if( v2.first == "rotation" ) { m_Controls->m_XrotBox->setValue(v2.second.get("x")); m_Controls->m_YrotBox->setValue(v2.second.get("y")); m_Controls->m_ZrotBox->setValue(v2.second.get("z")); } if( v2.first == "translation" ) { m_Controls->m_XtransBox->setValue(v2.second.get("x")); m_Controls->m_YtransBox->setValue(v2.second.get("y")); m_Controls->m_ZtransBox->setValue(v2.second.get("z")); } if( v2.first == "scale" ) { m_Controls->m_XscaleBox->setValue(v2.second.get("x")); m_Controls->m_YscaleBox->setValue(v2.second.get("y")); m_Controls->m_ZscaleBox->setValue(v2.second.get("z")); } } } if( v1.first == "image" ) { m_Controls->m_SizeX->setValue(v1.second.get("basic.size.x")); m_Controls->m_SizeY->setValue(v1.second.get("basic.size.y")); m_Controls->m_SizeZ->setValue(v1.second.get("basic.size.z")); m_Controls->m_SpacingX->setValue(v1.second.get("basic.spacing.x")); m_Controls->m_SpacingY->setValue(v1.second.get("basic.spacing.y")); m_Controls->m_SpacingZ->setValue(v1.second.get("basic.spacing.z")); m_Controls->m_NumGradientsBox->setValue(v1.second.get("basic.numgradients")); m_Controls->m_BvalueBox->setValue(v1.second.get("basic.bvalue")); m_Controls->m_AdvancedOptionsBox_2->setChecked(v1.second.get("showadvanced")); m_Controls->m_RepetitionsBox->setValue(v1.second.get("repetitions")); m_Controls->m_SignalScaleBox->setValue(v1.second.get("signalScale")); m_Controls->m_TEbox->setValue(v1.second.get("tEcho")); m_Controls->m_LineReadoutTimeBox->setValue(v1.second.get("tLine")); m_Controls->m_T2starBox->setValue(v1.second.get("tInhom")); m_Controls->m_FiberRadius->setValue(v1.second.get("axonRadius")); m_Controls->m_InterpolationShrink->setValue(v1.second.get("interpolationShrink")); m_Controls->m_RelaxationBox->setChecked(v1.second.get("doSimulateRelaxation")); m_Controls->m_EnforcePureFiberVoxelsBox->setChecked(v1.second.get("doDisablePartialVolume")); m_Controls->m_VolumeFractionsBox->setChecked(v1.second.get("outputvolumefractions")); m_Controls->m_AddNoise->setChecked(v1.second.get("artifacts.addnoise")); m_Controls->m_NoiseLevel->setValue(v1.second.get("artifacts.noisevariance")); m_Controls->m_AddGhosts->setChecked(v1.second.get("artifacts.addghost")); m_Controls->m_kOffsetBox->setValue(v1.second.get("artifacts.kspaceLineOffset")); m_Controls->m_AddAliasing->setChecked(v1.second.get("artifacts.addaliasing")); m_Controls->m_WrapBox->setValue(v1.second.get("artifacts.aliasingfactor")); m_Controls->m_AddDistortions->setChecked(v1.second.get("artifacts.distortions")); m_Controls->m_AddSpikes->setChecked(v1.second.get("artifacts.addspikes")); m_Controls->m_SpikeNumBox->setValue(v1.second.get("artifacts.spikesnum")); m_Controls->m_SpikeScaleBox->setValue(v1.second.get("artifacts.spikesscale")); m_Controls->m_AddEddy->setChecked(v1.second.get("artifacts.addeddy")); m_Controls->m_EddyGradientStrength->setValue(v1.second.get("artifacts.eddyStrength")); m_Controls->m_AddGibbsRinging->setChecked(v1.second.get("artifacts.addringing")); m_Controls->m_Compartment1Box->setCurrentIndex(v1.second.get("compartment1.index")); m_Controls->m_Compartment2Box->setCurrentIndex(v1.second.get("compartment2.index")); m_Controls->m_Compartment3Box->setCurrentIndex(v1.second.get("compartment3.index")); m_Controls->m_Compartment4Box->setCurrentIndex(v1.second.get("compartment4.index")); m_Controls->m_StickWidget1->SetD(v1.second.get("compartment1.stick.d")); m_Controls->m_StickWidget1->SetT2(v1.second.get("compartment1.stick.t2")); m_Controls->m_ZeppelinWidget1->SetD1(v1.second.get("compartment1.zeppelin.d1")); m_Controls->m_ZeppelinWidget1->SetD2(v1.second.get("compartment1.zeppelin.d2")); m_Controls->m_ZeppelinWidget1->SetT2(v1.second.get("compartment1.zeppelin.t2")); m_Controls->m_TensorWidget1->SetD1(v1.second.get("compartment1.tensor.d1")); m_Controls->m_TensorWidget1->SetD2(v1.second.get("compartment1.tensor.d2")); m_Controls->m_TensorWidget1->SetD3(v1.second.get("compartment1.tensor.d3")); m_Controls->m_TensorWidget1->SetT2(v1.second.get("compartment1.tensor.t2")); m_Controls->m_StickWidget2->SetD(v1.second.get("compartment2.stick.d")); m_Controls->m_StickWidget2->SetT2(v1.second.get("compartment2.stick.t2")); m_Controls->m_ZeppelinWidget2->SetD1(v1.second.get("compartment2.zeppelin.d1")); m_Controls->m_ZeppelinWidget2->SetD2(v1.second.get("compartment2.zeppelin.d2")); m_Controls->m_ZeppelinWidget2->SetT2(v1.second.get("compartment2.zeppelin.t2")); m_Controls->m_TensorWidget2->SetD1(v1.second.get("compartment2.tensor.d1")); m_Controls->m_TensorWidget2->SetD2(v1.second.get("compartment2.tensor.d2")); m_Controls->m_TensorWidget2->SetD3(v1.second.get("compartment2.tensor.d3")); m_Controls->m_TensorWidget2->SetT2(v1.second.get("compartment2.tensor.t2")); m_Controls->m_BallWidget1->SetD(v1.second.get("compartment3.ball.d")); m_Controls->m_BallWidget1->SetT2(v1.second.get("compartment3.ball.t2")); m_Controls->m_AstrosticksWidget1->SetD(v1.second.get("compartment3.astrosticks.d")); m_Controls->m_AstrosticksWidget1->SetT2(v1.second.get("compartment3.astrosticks.t2")); m_Controls->m_AstrosticksWidget1->SetRandomizeSticks(v1.second.get("compartment3.astrosticks.randomize")); m_Controls->m_DotWidget1->SetT2(v1.second.get("compartment3.dot.t2")); m_Controls->m_BallWidget2->SetD(v1.second.get("compartment4.ball.d")); m_Controls->m_BallWidget2->SetT2(v1.second.get("compartment4.ball.t2")); m_Controls->m_AstrosticksWidget2->SetD(v1.second.get("compartment4.astrosticks.d")); m_Controls->m_AstrosticksWidget2->SetT2(v1.second.get("compartment4.astrosticks.t2")); m_Controls->m_AstrosticksWidget2->SetRandomizeSticks(v1.second.get("compartment4.astrosticks.randomize")); m_Controls->m_DotWidget2->SetT2(v1.second.get("compartment4.dot.t2")); m_Controls->m_Comp4FractionBox->setValue(v1.second.get("compartment4.weight")); } } UpdateImageParameters(); } 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::OnAddMotion(int value) { if (value>0) m_Controls->m_MotionArtifactFrame->setVisible(true); else m_Controls->m_MotionArtifactFrame->setVisible(false); } void QmitkFiberfoxView::OnAddAliasing(int value) { if (value>0) m_Controls->m_AliasingFrame->setVisible(true); else m_Controls->m_AliasingFrame->setVisible(false); } void QmitkFiberfoxView::OnAddSpikes(int value) { if (value>0) m_Controls->m_SpikeFrame->setVisible(true); else m_Controls->m_SpikeFrame->setVisible(false); } 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::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(double 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); this->DisableCrosshairNavigation(); mitk::PlanarFigureInteractor::Pointer figureInteractor = dynamic_cast(node->GetDataInteractor().GetPointer()); if(figureInteractor.IsNull()) { figureInteractor = mitk::PlanarFigureInteractor::New(); us::Module* planarFigureModule = us::ModuleRegistry::GetModule( "PlanarFigure" ); figureInteractor->LoadStateMachine("PlanarFigureInteraction.xml", planarFigureModule ); figureInteractor->SetEventConfig( "PlanarFigureConfig.xml", planarFigureModule ); figureInteractor->SetDataNode( node ); } UpdateGui(); GetDataStorage()->Add(node, m_SelectedBundles.at(0)); } 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 (unsigned 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() { UpdateImageParameters(); if (m_SelectedBundles.empty()) { if (m_SelectedDWI.IsNotNull()) // add artifacts to existing diffusion weighted image { for (unsigned int i=0; i*>(m_SelectedImages.at(i)->GetData())) continue; m_SelectedDWI = m_SelectedImages.at(i); UpdateImageParameters(); mitk::DiffusionImage::Pointer diffImg = dynamic_cast*>(m_SelectedImages.at(i)->GetData()); mitk::RicianNoiseModel noiseModel; noiseModel.SetNoiseVariance(m_ImageGenParameters.ricianNoiseModel.GetNoiseVariance()); itk::AddArtifactsToDwiImageFilter< short >::Pointer filter = itk::AddArtifactsToDwiImageFilter< short >::New(); filter->SetInput(diffImg->GetVectorImage()); filter->SettLine(m_ImageGenParameters.tLine); filter->SetkOffset(m_ImageGenParameters.kspaceLineOffset); filter->SetNoiseModel(&noiseModel); filter->SetGradientList(m_ImageGenParameters.gradientDirections); filter->SetTE(m_ImageGenParameters.tEcho); filter->SetSimulateEddyCurrents(m_ImageGenParameters.doSimulateEddyCurrents); filter->SetEddyGradientStrength(m_ImageGenParameters.eddyStrength); filter->SetAddGibbsRinging(m_ImageGenParameters.addGibbsRinging); filter->SetFrequencyMap(m_ImageGenParameters.frequencyMap); filter->SetSpikeAmplitude(m_ImageGenParameters.spikeAmplitude); filter->SetSpikes(m_ImageGenParameters.spikes); filter->SetWrap(m_ImageGenParameters.wrap); filter->Update(); mitk::DiffusionImage::Pointer image = mitk::DiffusionImage::New(); image->SetVectorImage( filter->GetOutput() ); image->SetB_Value(diffImg->GetB_Value()); image->SetDirections(diffImg->GetDirections()); image->InitializeFromVectorImage(); m_ImageGenParameters.resultNode->SetData( image ); m_ImageGenParameters.resultNode->SetName(m_SelectedImages.at(i)->GetName()+m_ImageGenParameters.artifactModelString.toStdString()); GetDataStorage()->Add(m_ImageGenParameters.resultNode); } m_SelectedDWI = m_SelectedImages.front(); 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(m_ImageGenParameters.imageOrigin); 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; i(m_SelectedBundles.at(i)->GetData()); if (fiberBundle->GetNumFibers()<=0) continue; itk::TractsToDWIImageFilter< short >::Pointer tractsToDwiFilter = itk::TractsToDWIImageFilter< short >::New(); tractsToDwiFilter->SetSimulateEddyCurrents(m_ImageGenParameters.doSimulateEddyCurrents); tractsToDwiFilter->SetEddyGradientStrength(m_ImageGenParameters.eddyStrength); tractsToDwiFilter->SetAddGibbsRinging(m_ImageGenParameters.addGibbsRinging); tractsToDwiFilter->SetSimulateRelaxation(m_ImageGenParameters.doSimulateRelaxation); tractsToDwiFilter->SetImageRegion(m_ImageGenParameters.imageRegion); tractsToDwiFilter->SetSpacing(m_ImageGenParameters.imageSpacing); tractsToDwiFilter->SetOrigin(m_ImageGenParameters.imageOrigin); tractsToDwiFilter->SetDirectionMatrix(m_ImageGenParameters.imageDirection); tractsToDwiFilter->SetFiberBundle(fiberBundle); tractsToDwiFilter->SetFiberModels(m_ImageGenParameters.fiberModelList); tractsToDwiFilter->SetNonFiberModels(m_ImageGenParameters.nonFiberModelList); tractsToDwiFilter->SetNoiseModel(&m_ImageGenParameters.ricianNoiseModel); tractsToDwiFilter->SetKspaceArtifacts(m_ImageGenParameters.artifactList); tractsToDwiFilter->SetkOffset(m_ImageGenParameters.kspaceLineOffset); tractsToDwiFilter->SettLine(m_ImageGenParameters.tLine); tractsToDwiFilter->SettInhom(m_ImageGenParameters.tInhom); tractsToDwiFilter->SetTE(m_ImageGenParameters.tEcho); tractsToDwiFilter->SetNumberOfRepetitions(m_ImageGenParameters.repetitions); tractsToDwiFilter->SetEnforcePureFiberVoxels(m_ImageGenParameters.doDisablePartialVolume); tractsToDwiFilter->SetInterpolationShrink(m_ImageGenParameters.interpolationShrink); tractsToDwiFilter->SetFiberRadius(m_ImageGenParameters.axonRadius); tractsToDwiFilter->SetSignalScale(m_ImageGenParameters.signalScale); if (m_ImageGenParameters.interpolationShrink>0) tractsToDwiFilter->SetUseInterpolation(true); tractsToDwiFilter->SetTissueMask(m_ImageGenParameters.maskImage); tractsToDwiFilter->SetFrequencyMap(m_ImageGenParameters.frequencyMap); tractsToDwiFilter->SetSpikeAmplitude(m_ImageGenParameters.spikeAmplitude); tractsToDwiFilter->SetSpikes(m_ImageGenParameters.spikes); tractsToDwiFilter->SetWrap(m_ImageGenParameters.wrap); tractsToDwiFilter->SetAddMotionArtifact(m_ImageGenParameters.doAddMotion); tractsToDwiFilter->SetMaxTranslation(m_ImageGenParameters.translation); tractsToDwiFilter->SetMaxRotation(m_ImageGenParameters.rotation); + tractsToDwiFilter->SetRandomMotion(m_ImageGenParameters.randomMotion); tractsToDwiFilter->Update(); mitk::DiffusionImage::Pointer image = mitk::DiffusionImage::New(); image->SetVectorImage( tractsToDwiFilter->GetOutput() ); image->SetB_Value(m_ImageGenParameters.b_value); image->SetDirections(m_ImageGenParameters.gradientDirections); image->InitializeFromVectorImage(); m_ImageGenParameters.resultNode->SetData( image ); m_ImageGenParameters.resultNode->SetName(m_SelectedBundles.at(i)->GetName() +"_D"+QString::number(m_ImageGenParameters.imageRegion.GetSize(0)).toStdString() +"-"+QString::number(m_ImageGenParameters.imageRegion.GetSize(1)).toStdString() +"-"+QString::number(m_ImageGenParameters.imageRegion.GetSize(2)).toStdString() +"_S"+QString::number(m_ImageGenParameters.imageSpacing[0]).toStdString() +"-"+QString::number(m_ImageGenParameters.imageSpacing[1]).toStdString() +"-"+QString::number(m_ImageGenParameters.imageSpacing[2]).toStdString() +"_b"+QString::number(m_ImageGenParameters.b_value).toStdString() +"_"+m_ImageGenParameters.signalModelString.toStdString() +m_ImageGenParameters.artifactModelString.toStdString()); GetDataStorage()->Add(m_ImageGenParameters.resultNode, m_SelectedBundles.at(i)); m_ImageGenParameters.resultNode->SetProperty( "levelwindow", mitk::LevelWindowProperty::New(tractsToDwiFilter->GetLevelWindow()) ); if (m_Controls->m_VolumeFractionsBox->isChecked()) { std::vector< itk::TractsToDWIImageFilter< short >::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 = m_ImageGenParameters.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_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_ImageGenParameters.maskImage = 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; mitk::Image::Pointer image = dynamic_cast(node->GetData()); bool isbinary = false; node->GetPropertyValue("binary", isbinary); if (isbinary) { mitk::CastToItkImage(image, m_ImageGenParameters.maskImage); 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 (dynamic_cast(node->GetData())) { m_SelectedBundles.clear(); m_SelectedBundles2.clear(); } else if (dynamic_cast(node->GetData())) m_SelectedImages.clear(); 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->GetDataInteractor().GetPointer()); mitk::DataNode* nonConstNode = const_cast( node ); if(figureInteractor.IsNull()) { figureInteractor = mitk::PlanarFigureInteractor::New(); us::Module* planarFigureModule = us::ModuleRegistry::GetModule( "PlanarFigure" ); figureInteractor->LoadStateMachine("PlanarFigureInteraction.xml", planarFigureModule ); figureInteractor->SetEventConfig( "PlanarFigureConfig.xml", planarFigureModule ); figureInteractor->SetDataNode( nonConstNode ); } 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(); }