diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFitFibersToImageFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFitFibersToImageFilter.cpp index ba6dae923c..30bcdab106 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFitFibersToImageFilter.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFitFibersToImageFilter.cpp @@ -1,762 +1,760 @@ #include "itkFitFibersToImageFilter.h" #include namespace itk{ FitFibersToImageFilter::FitFibersToImageFilter() : m_PeakImage(nullptr) , m_MaskImage(nullptr) , m_FitIndividualFibers(true) , m_GradientTolerance(1e-5) , m_Lambda(1.0) , m_MaxIterations(20) , m_FiberSampling(10) , m_Coverage(0) , m_Overshoot(0) , m_RMSE(0.0) , m_FilterOutliers(true) , m_MeanWeight(1.0) , m_MedianWeight(1.0) , m_MinWeight(1.0) , m_MaxWeight(1.0) , m_Verbose(true) , m_DeepCopy(true) , m_ResampleFibers(true) , m_NumUnknowns(0) , m_NumResiduals(0) , m_NumCoveredDirections(0) , m_SignalModel(nullptr) , sz_x(0) , sz_y(0) , sz_z(0) , m_MeanTractDensity(0) , m_MeanSignal(0) , fiber_count(0) , m_Regularization(VnlCostFunction::REGU::Local_MSE) { this->SetNumberOfRequiredOutputs(3); } FitFibersToImageFilter::~FitFibersToImageFilter() { } void FitFibersToImageFilter::CreateDiffSystem() { sz_x = m_DiffImage->GetLargestPossibleRegion().GetSize(0); sz_y = m_DiffImage->GetLargestPossibleRegion().GetSize(1); sz_z = m_DiffImage->GetLargestPossibleRegion().GetSize(2); dim_four_size = m_DiffImage->GetVectorLength(); int num_voxels = sz_x*sz_y*sz_z; float minSpacing = 1; if(m_DiffImage->GetSpacing()[0]GetSpacing()[1] && m_DiffImage->GetSpacing()[0]GetSpacing()[2]) minSpacing = m_DiffImage->GetSpacing()[0]; else if (m_DiffImage->GetSpacing()[1] < m_DiffImage->GetSpacing()[2]) minSpacing = m_DiffImage->GetSpacing()[1]; else minSpacing = m_DiffImage->GetSpacing()[2]; if (m_ResampleFibers) for (unsigned int bundle=0; bundleGetDeepCopy(); m_Tractograms.at(bundle)->ResampleLinear(minSpacing/m_FiberSampling); std::cout.rdbuf (old); } m_NumResiduals = num_voxels * dim_four_size; MITK_INFO << "Num. unknowns: " << m_NumUnknowns; MITK_INFO << "Num. residuals: " << m_NumResiduals; MITK_INFO << "Creating system ..."; A.set_size(m_NumResiduals, m_NumUnknowns); b.set_size(m_NumResiduals); b.fill(0.0); m_MeanTractDensity = 0; m_MeanSignal = 0; m_NumCoveredDirections = 0; fiber_count = 0; vnl_vector voxel_indicator; voxel_indicator.set_size(sz_x*sz_y*sz_z); voxel_indicator.fill(0); for (unsigned int bundle=0; bundle polydata = m_Tractograms.at(bundle)->GetFiberPolyData(); for (int i=0; iGetNumFibers(); ++i) { vtkCell* cell = polydata->GetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); if (numPoints<2) MITK_INFO << "FIBER WITH ONLY ONE POINT ENCOUNTERED!"; for (int j=0; jGetPoint(j); PointType3 p; p[0]=p1[0]; p[1]=p1[1]; p[2]=p1[2]; itk::Index<3> idx3; m_DiffImage->TransformPhysicalPointToIndex(p, idx3); if (!m_DiffImage->GetLargestPossibleRegion().IsInside(idx3) || (m_MaskImage.IsNotNull() && m_MaskImage->GetPixel(idx3)==0)) continue; double* p2 = points->GetPoint(j+1); mitk::DiffusionSignalModel<>::GradientType fiber_dir; fiber_dir[0] = p[0]-p2[0]; fiber_dir[1] = p[1]-p2[1]; fiber_dir[2] = p[2]-p2[2]; fiber_dir.Normalize(); int x = idx3[0]; int y = idx3[1]; int z = idx3[2]; - - m_SignalModel->SetFiberDirection(fiber_dir); - mitk::DiffusionSignalModel<>::PixelType simulated_pixel = m_SignalModel->SimulateMeasurement(); + mitk::DiffusionSignalModel<>::PixelType simulated_pixel = m_SignalModel->SimulateMeasurement(fiber_dir); VectorImgType::PixelType measured_pixel = m_DiffImage->GetPixel(idx3); double simulated_mean = 0; double measured_mean = 0; int num_nonzero_g = 0; for (int g=0; gGetGradientDirection(g).GetNorm()GetLargestPossibleRegion().GetSize(0); sz_y = m_PeakImage->GetLargestPossibleRegion().GetSize(1); sz_z = m_PeakImage->GetLargestPossibleRegion().GetSize(2); dim_four_size = m_PeakImage->GetLargestPossibleRegion().GetSize(3)/3 + 1; // +1 for zero - peak int num_voxels = sz_x*sz_y*sz_z; float minSpacing = 1; if(m_PeakImage->GetSpacing()[0]GetSpacing()[1] && m_PeakImage->GetSpacing()[0]GetSpacing()[2]) minSpacing = m_PeakImage->GetSpacing()[0]; else if (m_PeakImage->GetSpacing()[1] < m_PeakImage->GetSpacing()[2]) minSpacing = m_PeakImage->GetSpacing()[1]; else minSpacing = m_PeakImage->GetSpacing()[2]; if (m_ResampleFibers) for (unsigned int bundle=0; bundleGetDeepCopy(); m_Tractograms.at(bundle)->ResampleLinear(minSpacing/m_FiberSampling); std::cout.rdbuf (old); } m_NumResiduals = num_voxels * dim_four_size; MITK_INFO << "Num. unknowns: " << m_NumUnknowns; MITK_INFO << "Num. residuals: " << m_NumResiduals; MITK_INFO << "Creating system ..."; A.set_size(m_NumResiduals, m_NumUnknowns); b.set_size(m_NumResiduals); b.fill(0.0); m_MeanTractDensity = 0; m_MeanSignal = 0; m_NumCoveredDirections = 0; fiber_count = 0; for (unsigned int bundle=0; bundle polydata = m_Tractograms.at(bundle)->GetFiberPolyData(); for (int i=0; iGetNumFibers(); ++i) { vtkCell* cell = polydata->GetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); if (numPoints<2) MITK_INFO << "FIBER WITH ONLY ONE POINT ENCOUNTERED!"; for (int j=0; jGetPoint(j); PointType4 p; p[0]=p1[0]; p[1]=p1[1]; p[2]=p1[2]; p[3]=0; itk::Index<4> idx4; m_PeakImage->TransformPhysicalPointToIndex(p, idx4); itk::Index<3> idx3; idx3[0] = idx4[0]; idx3[1] = idx4[1]; idx3[2] = idx4[2]; if (!m_PeakImage->GetLargestPossibleRegion().IsInside(idx4) || (m_MaskImage.IsNotNull() && m_MaskImage->GetPixel(idx3)==0)) continue; double* p2 = points->GetPoint(j+1); vnl_vector_fixed fiber_dir; fiber_dir[0] = p[0]-p2[0]; fiber_dir[1] = p[1]-p2[1]; fiber_dir[2] = p[2]-p2[2]; fiber_dir.normalize(); double w = 1; int peak_id = dim_four_size-1; vnl_vector_fixed odf_peak = GetClosestPeak(idx4, m_PeakImage, fiber_dir, peak_id, w); float peak_mag = odf_peak.magnitude(); int x = idx4[0]; int y = idx4[1]; int z = idx4[2]; unsigned int linear_index = x + sz_x*y + sz_x*sz_y*z + sz_x*sz_y*sz_z*peak_id; if (b[linear_index] == 0 && peak_idGetNumFibers(); } else m_FilterOutliers = false; if (m_NumUnknowns<1) { MITK_INFO << "No fibers in tractogram."; return; } fiber_count = 0; sz_x = 0; sz_y = 0; sz_z = 0; m_MeanTractDensity = 0; m_MeanSignal = 0; if (m_PeakImage.IsNotNull()) CreatePeakSystem(); else if (m_DiffImage.IsNotNull()) CreateDiffSystem(); else mitkThrow() << "No input image set!"; double init_lambda = fiber_count; // initialization for lambda estimation itk::TimeProbe clock; clock.Start(); cost = VnlCostFunction(m_NumUnknowns); cost.SetProblem(A, b, init_lambda, m_Regularization); m_Weights.set_size(m_NumUnknowns); m_Weights.fill( 0.0 ); vnl_lbfgsb minimizer(cost); vnl_vector l; l.set_size(m_NumUnknowns); l.fill(0); vnl_vector bound_selection; bound_selection.set_size(m_NumUnknowns); bound_selection.fill(1); minimizer.set_bound_selection(bound_selection); minimizer.set_lower_bound(l); minimizer.set_projected_gradient_tolerance(m_GradientTolerance); MITK_INFO << "Regularization type: " << m_Regularization; if (m_Regularization!=VnlCostFunction::REGU::NONE) // REMOVE FOR NEW FIT AND SET cost.m_Lambda = m_Lambda { MITK_INFO << "Estimating regularization"; minimizer.set_trace(false); minimizer.set_max_function_evals(2); minimizer.minimize(m_Weights); vnl_vector dx; dx.set_size(m_NumUnknowns); dx.fill(0.0); cost.calc_regularization_gradient(m_Weights, dx); double r = dx.magnitude()/m_Weights.magnitude(); // wtf??? cost.m_Lambda *= m_Lambda*55.0/r; MITK_INFO << r << " - " << m_Lambda*55.0/r; if (cost.m_Lambda>10e7) { MITK_INFO << "Regularization estimation failed. Using default value."; cost.m_Lambda = fiber_count; } } MITK_INFO << "Using regularization factor of " << cost.m_Lambda << " (λ: " << m_Lambda << ")"; MITK_INFO << "Fitting fibers"; minimizer.set_trace(m_Verbose); minimizer.set_max_function_evals(m_MaxIterations); minimizer.minimize(m_Weights); std::vector< double > weights; if (m_FilterOutliers) { for (auto w : m_Weights) weights.push_back(w); std::sort(weights.begin(), weights.end()); MITK_INFO << "Setting upper weight bound to " << weights.at(m_NumUnknowns*0.99); vnl_vector u; u.set_size(m_NumUnknowns); u.fill(weights.at(m_NumUnknowns*0.99)); minimizer.set_upper_bound(u); bound_selection.fill(2); minimizer.set_bound_selection(bound_selection); minimizer.minimize(m_Weights); weights.clear(); } for (auto w : m_Weights) weights.push_back(w); std::sort(weights.begin(), weights.end()); m_MeanWeight = m_Weights.mean(); m_MedianWeight = weights.at(m_NumUnknowns*0.5); m_MinWeight = weights.at(0); m_MaxWeight = weights.at(m_NumUnknowns-1); MITK_INFO << "*************************"; MITK_INFO << "Weight statistics"; MITK_INFO << "Sum: " << m_Weights.sum(); MITK_INFO << "Mean: " << m_MeanWeight; MITK_INFO << "1% quantile: " << weights.at(m_NumUnknowns*0.01); MITK_INFO << "5% quantile: " << weights.at(m_NumUnknowns*0.05); MITK_INFO << "25% quantile: " << weights.at(m_NumUnknowns*0.25); MITK_INFO << "Median: " << m_MedianWeight; MITK_INFO << "75% quantile: " << weights.at(m_NumUnknowns*0.75); MITK_INFO << "95% quantile: " << weights.at(m_NumUnknowns*0.95); MITK_INFO << "99% quantile: " << weights.at(m_NumUnknowns*0.99); MITK_INFO << "Min: " << m_MinWeight; MITK_INFO << "Max: " << m_MaxWeight; MITK_INFO << "*************************"; MITK_INFO << "NumEvals: " << minimizer.get_num_evaluations(); MITK_INFO << "NumIterations: " << minimizer.get_num_iterations(); MITK_INFO << "Residual cost: " << minimizer.get_end_error(); m_RMSE = cost.S->get_rms_error(m_Weights); MITK_INFO << "Final RMS: " << m_RMSE; clock.Stop(); int h = clock.GetTotal()/3600; int m = ((int)clock.GetTotal()%3600)/60; int s = (int)clock.GetTotal()%60; MITK_INFO << "Optimization took " << h << "h, " << m << "m and " << s << "s"; MITK_INFO << "Weighting fibers"; m_RmsDiffPerBundle.set_size(m_Tractograms.size()); std::streambuf *old = cout.rdbuf(); // <-- save std::stringstream ss; std::cout.rdbuf (ss.rdbuf()); if (m_FitIndividualFibers) { unsigned int fiber_count = 0; for (unsigned int bundle=0; bundle temp_weights; temp_weights.set_size(m_Weights.size()); temp_weights.copy_in(m_Weights.data_block()); for (int i=0; iGetNumFibers(); i++) { m_Tractograms.at(bundle)->SetFiberWeight(i, m_Weights[fiber_count]); temp_weights[fiber_count] = 0; ++fiber_count; } double d_rms = cost.S->get_rms_error(temp_weights) - m_RMSE; m_RmsDiffPerBundle[bundle] = d_rms; m_Tractograms.at(bundle)->Compress(0.1); m_Tractograms.at(bundle)->ColorFibersByFiberWeights(false, true); } } else { for (unsigned int i=0; i temp_weights; temp_weights.set_size(m_Weights.size()); temp_weights.copy_in(m_Weights.data_block()); temp_weights[i] = 0; double d_rms = cost.S->get_rms_error(temp_weights) - m_RMSE; m_RmsDiffPerBundle[i] = d_rms; m_Tractograms.at(i)->SetFiberWeights(m_Weights[i]); m_Tractograms.at(i)->Compress(0.1); m_Tractograms.at(i)->ColorFibersByFiberWeights(false, true); } } std::cout.rdbuf (old); // transform back A *= m_MeanSignal/100.0; b *= m_MeanSignal/100.0; MITK_INFO << "Generating output images ..."; if (m_PeakImage.IsNotNull()) GenerateOutputPeakImages(); else if (m_DiffImage.IsNotNull()) GenerateOutputDiffImages(); m_Coverage = m_Coverage/m_MeanSignal; m_Overshoot = m_Overshoot/m_MeanSignal; MITK_INFO << std::fixed << "Coverage: " << setprecision(2) << 100.0*m_Coverage << "%"; MITK_INFO << std::fixed << "Overshoot: " << setprecision(2) << 100.0*m_Overshoot << "%"; } void FitFibersToImageFilter::GenerateOutputDiffImages() { VectorImgType::PixelType pix; pix.SetSize(m_DiffImage->GetVectorLength()); pix.Fill(0); itk::ImageDuplicator< VectorImgType >::Pointer duplicator = itk::ImageDuplicator< VectorImgType >::New(); duplicator->SetInputImage(m_DiffImage); duplicator->Update(); m_UnderexplainedImageDiff = duplicator->GetOutput(); m_UnderexplainedImageDiff->FillBuffer(pix); duplicator->SetInputImage(m_UnderexplainedImageDiff); duplicator->Update(); m_OverexplainedImageDiff = duplicator->GetOutput(); m_OverexplainedImageDiff->FillBuffer(pix); duplicator->SetInputImage(m_OverexplainedImageDiff); duplicator->Update(); m_ResidualImageDiff = duplicator->GetOutput(); m_ResidualImageDiff->FillBuffer(pix); duplicator->SetInputImage(m_ResidualImageDiff); duplicator->Update(); m_FittedImageDiff = duplicator->GetOutput(); m_FittedImageDiff->FillBuffer(pix); vnl_vector fitted_b; fitted_b.set_size(b.size()); cost.S->multiply(m_Weights, fitted_b); itk::ImageRegionIterator it1 = itk::ImageRegionIterator(m_DiffImage, m_DiffImage->GetLargestPossibleRegion()); itk::ImageRegionIterator it2 = itk::ImageRegionIterator(m_FittedImageDiff, m_FittedImageDiff->GetLargestPossibleRegion()); itk::ImageRegionIterator it3 = itk::ImageRegionIterator(m_ResidualImageDiff, m_ResidualImageDiff->GetLargestPossibleRegion()); itk::ImageRegionIterator it4 = itk::ImageRegionIterator(m_UnderexplainedImageDiff, m_UnderexplainedImageDiff->GetLargestPossibleRegion()); itk::ImageRegionIterator it5 = itk::ImageRegionIterator(m_OverexplainedImageDiff, m_OverexplainedImageDiff->GetLargestPossibleRegion()); m_MeanSignal = 0; m_Coverage = 0; m_Overshoot = 0; while( !it2.IsAtEnd() ) { itk::Index<3> idx3 = it2.GetIndex(); VectorImgType::PixelType original_pix =it1.Get(); VectorImgType::PixelType fitted_pix =it2.Get(); VectorImgType::PixelType residual_pix =it3.Get(); VectorImgType::PixelType underexplained_pix =it4.Get(); VectorImgType::PixelType overexplained_pix =it5.Get(); int num_nonzero_g = 0; double original_mean = 0; for (int g=0; gGetGradientDirection(g).GetNorm()>=mitk::eps ) { original_mean += original_pix[g]; ++num_nonzero_g; } } original_mean /= num_nonzero_g; for (int g=0; g=0) { underexplained_pix[g] = residual_pix[g]; m_Coverage += fitted_b[linear_index] + original_mean; } m_MeanSignal += b[linear_index] + original_mean; } it2.Set(fitted_pix); it3.Set(residual_pix); it4.Set(underexplained_pix); it5.Set(overexplained_pix); ++it1; ++it2; ++it3; ++it4; ++it5; } } VnlCostFunction::REGU FitFibersToImageFilter::GetRegularization() const { return m_Regularization; } void FitFibersToImageFilter::SetRegularization(const VnlCostFunction::REGU &Regularization) { m_Regularization = Regularization; } void FitFibersToImageFilter::GenerateOutputPeakImages() { itk::ImageDuplicator< PeakImgType >::Pointer duplicator = itk::ImageDuplicator< PeakImgType >::New(); duplicator->SetInputImage(m_PeakImage); duplicator->Update(); m_UnderexplainedImage = duplicator->GetOutput(); m_UnderexplainedImage->FillBuffer(0.0); duplicator->SetInputImage(m_UnderexplainedImage); duplicator->Update(); m_OverexplainedImage = duplicator->GetOutput(); m_OverexplainedImage->FillBuffer(0.0); duplicator->SetInputImage(m_OverexplainedImage); duplicator->Update(); m_ResidualImage = duplicator->GetOutput(); m_ResidualImage->FillBuffer(0.0); duplicator->SetInputImage(m_ResidualImage); duplicator->Update(); m_FittedImage = duplicator->GetOutput(); m_FittedImage->FillBuffer(0.0); vnl_vector fitted_b; fitted_b.set_size(b.size()); cost.S->multiply(m_Weights, fitted_b); for (unsigned int r=0; r idx4; unsigned int linear_index = r; idx4[0] = linear_index % sz_x; linear_index /= sz_x; idx4[1] = linear_index % sz_y; linear_index /= sz_y; idx4[2] = linear_index % sz_z; linear_index /= sz_z; int peak_id = linear_index % dim_four_size; if (peak_id peak_dir; idx4[3] = peak_id*3; peak_dir[0] = m_PeakImage->GetPixel(idx4); idx4[3] += 1; peak_dir[1] = m_PeakImage->GetPixel(idx4); idx4[3] += 1; peak_dir[2] = m_PeakImage->GetPixel(idx4); peak_dir.normalize(); peak_dir *= fitted_b[r]; idx4[3] = peak_id*3; m_FittedImage->SetPixel(idx4, peak_dir[0]); idx4[3] += 1; m_FittedImage->SetPixel(idx4, peak_dir[1]); idx4[3] += 1; m_FittedImage->SetPixel(idx4, peak_dir[2]); } } m_MeanSignal = 0; m_Coverage = 0; m_Overshoot = 0; itk::Index<4> idx4; for (idx4[0]=0; idx4[0] idx3; idx3[0] = idx4[0]; idx3[1] = idx4[1]; idx3[2] = idx4[2]; if (m_MaskImage.IsNotNull() && m_MaskImage->GetPixel(idx3)==0) continue; vnl_vector_fixed peak_dir; vnl_vector_fixed fitted_dir; vnl_vector_fixed overshoot_dir; for (idx4[3]=0; idx4[3]<(itk::IndexValueType)m_PeakImage->GetLargestPossibleRegion().GetSize(3); ++idx4[3]) { peak_dir[idx4[3]%3] = m_PeakImage->GetPixel(idx4); fitted_dir[idx4[3]%3] = m_FittedImage->GetPixel(idx4); m_ResidualImage->SetPixel(idx4, m_PeakImage->GetPixel(idx4) - m_FittedImage->GetPixel(idx4)); if (idx4[3]%3==2) { m_MeanSignal += peak_dir.magnitude(); itk::Index<4> tidx= idx4; if (peak_dir.magnitude()>fitted_dir.magnitude()) { m_Coverage += fitted_dir.magnitude(); m_UnderexplainedImage->SetPixel(tidx, peak_dir[2]-fitted_dir[2]); tidx[3]--; m_UnderexplainedImage->SetPixel(tidx, peak_dir[1]-fitted_dir[1]); tidx[3]--; m_UnderexplainedImage->SetPixel(tidx, peak_dir[0]-fitted_dir[0]); } else { overshoot_dir[0] = fitted_dir[0]-peak_dir[0]; overshoot_dir[1] = fitted_dir[1]-peak_dir[1]; overshoot_dir[2] = fitted_dir[2]-peak_dir[2]; m_Coverage += peak_dir.magnitude(); m_Overshoot += overshoot_dir.magnitude(); m_OverexplainedImage->SetPixel(tidx, overshoot_dir[2]); tidx[3]--; m_OverexplainedImage->SetPixel(tidx, overshoot_dir[1]); tidx[3]--; m_OverexplainedImage->SetPixel(tidx, overshoot_dir[0]); } } } } } vnl_vector_fixed FitFibersToImageFilter::GetClosestPeak(itk::Index<4> idx, PeakImgType::Pointer peak_image , vnl_vector_fixed fiber_dir, int& id, double& w ) { int m_NumDirs = peak_image->GetLargestPossibleRegion().GetSize()[3]/3; vnl_vector_fixed out_dir; out_dir.fill(0); float angle = 0.9; for (int i=0; i dir; idx[3] = i*3; dir[0] = peak_image->GetPixel(idx); idx[3] += 1; dir[1] = peak_image->GetPixel(idx); idx[3] += 1; dir[2] = peak_image->GetPixel(idx); float mag = dir.magnitude(); if (magangle) { angle = fabs(a); w = angle; if (a<0) out_dir = -dir; else out_dir = dir; out_dir *= mag; id = i; } } return out_dir; } std::vector FitFibersToImageFilter::GetTractograms() const { return m_Tractograms; } void FitFibersToImageFilter::SetTractograms(const std::vector &tractograms) { m_Tractograms = tractograms; } void FitFibersToImageFilter::SetSignalModel(mitk::DiffusionSignalModel<> *SignalModel) { m_SignalModel = SignalModel; } } diff --git a/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkAstroStickModel.cpp b/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkAstroStickModel.cpp index 24de236142..3a818c7cea 100644 --- a/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkAstroStickModel.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkAstroStickModel.cpp @@ -1,123 +1,123 @@ /*=================================================================== 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 #include #include using namespace mitk; template< class ScalarType > AstroStickModel< ScalarType >::AstroStickModel() : m_BValue(1000) , m_Diffusivity(0.001) , m_NumSticks(42) , m_RandomizeSticks(false) { this->m_RandGen = itk::Statistics::MersenneTwisterRandomVariateGenerator::New(); this->m_RandGen->SetSeed(); vnl_matrix_fixed* sticks = itk::PointShell<42, vnl_matrix_fixed >::DistributePointShell(); for (unsigned int i=0; iget(0,i); stick[1] = sticks->get(1,i); stick[2] = sticks->get(2,i); stick.Normalize(); m_Sticks.push_back(stick); } } template< class ScalarType > AstroStickModel< ScalarType >::~AstroStickModel() { } template< class ScalarType > -ScalarType AstroStickModel< ScalarType >::SimulateMeasurement(unsigned int dir) +ScalarType AstroStickModel< ScalarType >::SimulateMeasurement(unsigned int dir, GradientType& ) { ScalarType signal = 0; if (dir>=this->m_GradientList.size()) return signal; ScalarType b = -m_BValue*m_Diffusivity; if (m_RandomizeSticks) // random number of sticks m_NumSticks = 30 + this->m_RandGen->GetIntegerVariate()%31; GradientType g = this->m_GradientList[dir]; if (g.GetNorm()>0.0001) // is weighted direction { for (unsigned int j=0; j typename AstroStickModel< ScalarType >::GradientType AstroStickModel< ScalarType >::GetRandomDirection() { GradientType vec; vec[0] = this->m_RandGen->GetNormalVariate(); vec[1] = this->m_RandGen->GetNormalVariate(); vec[2] = this->m_RandGen->GetNormalVariate(); vec.Normalize(); return vec; } template< class ScalarType > -typename AstroStickModel< ScalarType >::PixelType AstroStickModel< ScalarType >::SimulateMeasurement() +typename AstroStickModel< ScalarType >::PixelType AstroStickModel< ScalarType >::SimulateMeasurement(GradientType& ) { PixelType signal; signal.SetSize(this->m_GradientList.size()); ScalarType b = -m_BValue*m_Diffusivity; if (m_RandomizeSticks) m_NumSticks = 30 + this->m_RandGen->GetIntegerVariate()%31; for( unsigned int i=0; im_GradientList.size(); i++) { GradientType g = this->m_GradientList[i]; if (g.GetNorm()>0.0001) { for (unsigned int j=0; j #include namespace mitk { /** * \brief Generates the diffusion signal using a collection of idealised cylinder with zero radius: e^(-bd(ng)²) * */ template< class ScalarType = double > class AstroStickModel : public DiffusionSignalModel< ScalarType > { public: AstroStickModel(); template< class OtherType >AstroStickModel(AstroStickModel* model) { this->m_CompartmentId = model->m_CompartmentId; this->m_T2 = model->GetT2(); this->m_FiberDirection = model->GetFiberDirection(); this->m_GradientList = model->GetGradientList(); this->m_VolumeFractionImage = model->GetVolumeFractionImage(); this->m_RandGen = model->GetRandomGenerator(); this->m_BValue = model->GetBvalue(); this->m_Diffusivity = model->GetDiffusivity(); this->m_Sticks = model->GetSticks(); this->m_NumSticks = model->GetNumSticks(); this->m_RandomizeSticks = model->GetRandomizeSticks(); } ~AstroStickModel(); typedef typename DiffusionSignalModel< ScalarType >::PixelType PixelType; typedef typename DiffusionSignalModel< ScalarType >::GradientType GradientType; typedef typename DiffusionSignalModel< ScalarType >::GradientListType GradientListType; /** Actual signal generation **/ - PixelType SimulateMeasurement(); - ScalarType SimulateMeasurement(unsigned int dir); - - - void SetFiberDirection(GradientType fiberDirection){ this->m_FiberDirection = fiberDirection; } + PixelType SimulateMeasurement(GradientType& fiberDirection); + ScalarType SimulateMeasurement(unsigned int dir, GradientType& fiberDirection); void SetRandomizeSticks(bool randomize=true){ m_RandomizeSticks=randomize; } ///< Random stick configuration in each voxel bool GetRandomizeSticks() { return m_RandomizeSticks; } void SetBvalue(double bValue) { m_BValue = bValue; } ///< b-value used to generate the artificial signal double GetBvalue() { return m_BValue; } void SetDiffusivity(double diffusivity) { m_Diffusivity = diffusivity; } ///< Scalar diffusion constant double GetDiffusivity() { return m_Diffusivity; } void SetNumSticks(unsigned int order) { vnl_matrix sticks; switch (order) { case 1: m_NumSticks = 12; sticks = itk::PointShell<12, vnl_matrix_fixed >::DistributePointShell()->as_matrix(); break; case 2: m_NumSticks = 42; sticks = itk::PointShell<42, vnl_matrix_fixed >::DistributePointShell()->as_matrix(); break; case 3: m_NumSticks = 92; sticks = itk::PointShell<92, vnl_matrix_fixed >::DistributePointShell()->as_matrix(); break; case 4: m_NumSticks = 162; sticks = itk::PointShell<162, vnl_matrix_fixed >::DistributePointShell()->as_matrix(); break; case 5: m_NumSticks = 252; sticks = itk::PointShell<252, vnl_matrix_fixed >::DistributePointShell()->as_matrix(); break; default: m_NumSticks = 42; sticks = itk::PointShell<42, vnl_matrix_fixed >::DistributePointShell()->as_matrix(); break; } for (int i=0; i #include #include using namespace mitk; template< class ScalarType > BallModel< ScalarType >::BallModel() : m_Diffusivity(0.001) , m_BValue(1000) { } template< class ScalarType > BallModel< ScalarType >::~BallModel() { } template< class ScalarType > -ScalarType BallModel< ScalarType >::SimulateMeasurement(unsigned int dir) +ScalarType BallModel< ScalarType >::SimulateMeasurement(unsigned int dir, GradientType& ) { ScalarType signal = 0; if (dir>=this->m_GradientList.size()) return signal; GradientType g = this->m_GradientList[dir]; ScalarType bVal = g.GetNorm(); bVal *= bVal; if (bVal>0.0001) signal = std::exp( -m_BValue * bVal * m_Diffusivity ); else signal = 1; return signal; } template< class ScalarType > -typename BallModel< ScalarType >::PixelType BallModel< ScalarType >::SimulateMeasurement() +typename BallModel< ScalarType >::PixelType BallModel< ScalarType >::SimulateMeasurement(GradientType& ) { PixelType signal; signal.SetSize(this->m_GradientList.size()); for( unsigned int i=0; im_GradientList.size(); i++) { GradientType g = this->m_GradientList[i]; ScalarType bVal = g.GetNorm(); bVal *= bVal; if (bVal>0.0001) signal[i] = std::exp( -m_BValue * bVal * m_Diffusivity ); else signal[i] = 1; } return signal; } diff --git a/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkBallModel.h b/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkBallModel.h index 3629525dd3..b59865454e 100644 --- a/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkBallModel.h +++ b/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkBallModel.h @@ -1,76 +1,74 @@ /*=================================================================== 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_BallModel_H #define _MITK_BallModel_H #include namespace mitk { /** * \brief Generates direction independent diffusion measurement employing a scalar diffusion constant d: e^(-bd) * */ template< class ScalarType = double > class BallModel : public DiffusionSignalModel< ScalarType > { public: BallModel(); template< class OtherType >BallModel(BallModel* model) { this->m_CompartmentId = model->m_CompartmentId; this->m_T2 = model->GetT2(); this->m_FiberDirection = model->GetFiberDirection(); this->m_GradientList = model->GetGradientList(); this->m_VolumeFractionImage = model->GetVolumeFractionImage(); this->m_RandGen = model->GetRandomGenerator(); this->m_BValue = model->GetBvalue(); this->m_Diffusivity = model->GetDiffusivity(); } ~BallModel(); typedef typename DiffusionSignalModel< ScalarType >::PixelType PixelType; typedef typename DiffusionSignalModel< ScalarType >::GradientType GradientType; typedef typename DiffusionSignalModel< ScalarType >::GradientListType GradientListType; /** Actual signal generation **/ - PixelType SimulateMeasurement(); - ScalarType SimulateMeasurement(unsigned int dir); + PixelType SimulateMeasurement(GradientType& fiberDirection); + ScalarType SimulateMeasurement(unsigned int dir, GradientType& fiberDirection); void SetDiffusivity(double D) { m_Diffusivity = D; } double GetDiffusivity() { return m_Diffusivity; } void SetBvalue(double bValue) { m_BValue = bValue; } ///< b-value used to generate the artificial signal double GetBvalue() { return m_BValue; } - void SetFiberDirection(GradientType fiberDirection){ this->m_FiberDirection = fiberDirection; } - protected: double m_Diffusivity; ///< Scalar diffusion constant double m_BValue; ///< b-value used to generate the artificial signal }; } #include "mitkBallModel.cpp" #endif diff --git a/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkDiffusionSignalModel.h b/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkDiffusionSignalModel.h index 7ad90a30dd..f9cb445093 100644 --- a/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkDiffusionSignalModel.h +++ b/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkDiffusionSignalModel.h @@ -1,113 +1,108 @@ /*=================================================================== 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_DiffusionSignalModel_H #define _MITK_DiffusionSignalModel_H #include #include #include #include #include #include #include namespace mitk { /** * \brief Abstract class for diffusion signal models * */ template< class ScalarType = double > class DiffusionSignalModel { public: DiffusionSignalModel() : m_T2(100) , m_T1(0) {} ~DiffusionSignalModel(){} typedef itk::Image ItkDoubleImgType; typedef itk::VariableLengthVector< ScalarType > PixelType; typedef itk::Vector GradientType; typedef std::vector GradientListType; typedef itk::Statistics::MersenneTwisterRandomVariateGenerator ItkRandGenType; typedef mitk::DiffusionPropertyHelper DPH; /** Realizes actual signal generation. Has to be implemented in subclass. **/ - virtual PixelType SimulateMeasurement() = 0; - virtual ScalarType SimulateMeasurement(unsigned int dir) = 0; - - - virtual void SetFiberDirection(GradientType fiberDirection) = 0; - GradientType GetFiberDirection(){ return m_FiberDirection; } + virtual PixelType SimulateMeasurement(GradientType& fiberDirection) = 0; + virtual ScalarType SimulateMeasurement(unsigned int dir, GradientType& fiberDirection) = 0; void SetGradientList(DPH::GradientDirectionsContainerType* gradients) { m_GradientList.clear(); for ( unsigned int i=0; iSize(); ++i ) { DPH::GradientDirectionType g_vnl = gradients->GetElement(i); GradientType g_itk; g_itk[0] = g_vnl[0]; g_itk[1] = g_vnl[1]; g_itk[2] = g_vnl[2]; m_GradientList.push_back(g_itk); } } void SetGradientList(GradientListType gradientList) { this->m_GradientList = gradientList; } GradientListType GetGradientList(){ return m_GradientList; } GradientType GetGradientDirection(int i) { return m_GradientList.at(i); } void SetT2(double T2) { m_T2 = T2; } double GetT2() { return m_T2; } void SetT1(double T1) { m_T1 = T1; } double GetT1() { return m_T1; } void SetVolumeFractionImage(ItkDoubleImgType::Pointer img){ m_VolumeFractionImage = img; } ItkDoubleImgType::Pointer GetVolumeFractionImage(){ return m_VolumeFractionImage; } void SetRandomGenerator(ItkRandGenType::Pointer randgen){ m_RandGen = randgen; } ItkRandGenType::Pointer GetRandomGenerator(){ return m_RandGen; } void SetSeed(int s) ///< set seed for random generator { if (m_RandGen.IsNull()) m_RandGen = itk::Statistics::MersenneTwisterRandomVariateGenerator::New(); m_RandGen->SetSeed(s); } unsigned int m_CompartmentId; ///< GUI flag. Which compartment is this model assigned to? protected: - GradientType m_FiberDirection; ///< Needed to generate anisotropc signal to determin direction of anisotropy GradientListType m_GradientList; ///< Diffusion gradient direction container double m_T2; ///< Tissue specific transversal relaxation time double m_T1; ///< Tissue specific longitudinal relaxation time ItkDoubleImgType::Pointer m_VolumeFractionImage; ///< Tissue specific volume fraction for each voxel (only relevant for non fiber compartments) ItkRandGenType::Pointer m_RandGen; ///< Random number generator }; } #endif diff --git a/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkDotModel.cpp b/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkDotModel.cpp index 7d91147d48..1eb8d4af29 100644 --- a/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkDotModel.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkDotModel.cpp @@ -1,47 +1,47 @@ /*=================================================================== 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 #include #include using namespace mitk; template< class ScalarType > DotModel< ScalarType >::DotModel() { } template< class ScalarType > DotModel< ScalarType >::~DotModel() { } template< class ScalarType > -ScalarType DotModel< ScalarType >::SimulateMeasurement(unsigned int /*dir*/) +ScalarType DotModel< ScalarType >::SimulateMeasurement(unsigned int /*dir*/, GradientType& ) { return 1; } template< class ScalarType > -typename DotModel< ScalarType >::PixelType DotModel< ScalarType >::SimulateMeasurement() +typename DotModel< ScalarType >::PixelType DotModel< ScalarType >::SimulateMeasurement(GradientType& ) { PixelType signal; signal.SetSize(this->m_GradientList.size()); signal.Fill(1); return signal; } diff --git a/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkDotModel.h b/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkDotModel.h index 1410833182..90272bb3c1 100644 --- a/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkDotModel.h +++ b/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkDotModel.h @@ -1,65 +1,63 @@ /*=================================================================== 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_DotModel_H #define _MITK_DotModel_H #include namespace mitk { /** * \brief Generates constant direction independent signal. * */ template< class ScalarType = double > class DotModel : public DiffusionSignalModel< ScalarType > { public: DotModel(); template< class OtherType >DotModel(DotModel* model) { this->m_CompartmentId = model->m_CompartmentId; this->m_T2 = model->GetT2(); this->m_FiberDirection = model->GetFiberDirection(); this->m_GradientList = model->GetGradientList(); this->m_VolumeFractionImage = model->GetVolumeFractionImage(); this->m_RandGen = model->GetRandomGenerator(); } ~DotModel(); typedef typename DiffusionSignalModel< ScalarType >::PixelType PixelType; typedef typename DiffusionSignalModel< ScalarType >::GradientType GradientType; typedef typename DiffusionSignalModel< ScalarType >::GradientListType GradientListType; /** Actual signal generation **/ - PixelType SimulateMeasurement(); - ScalarType SimulateMeasurement(unsigned int dir); - - void SetFiberDirection(GradientType fiberDirection){ this->m_FiberDirection = fiberDirection; } + PixelType SimulateMeasurement(GradientType& fiberDirection); + ScalarType SimulateMeasurement(unsigned int dir, GradientType& fiberDirection); protected: }; } #include "mitkDotModel.cpp" #endif diff --git a/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkRawShModel.cpp b/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkRawShModel.cpp index fc0b8c239d..a32bb2b460 100644 --- a/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkRawShModel.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkRawShModel.cpp @@ -1,390 +1,394 @@ /*=================================================================== 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 #include #include #include #include #include #include #include #include #include #include #include using namespace mitk; using namespace boost::math; template< class ScalarType > RawShModel< ScalarType >::RawShModel() - : m_ShOrder(0) - , m_ModelIndex(-1) - , m_MaxNumKernels(1000) + : m_ShOrder(0) + , m_ModelIndex(-1) + , m_MaxNumKernels(1000) { - this->m_RandGen = itk::Statistics::MersenneTwisterRandomVariateGenerator::New(); - this->m_RandGen->SetSeed(); - m_AdcRange.first = 0; - m_AdcRange.second = 0.004; - m_FaRange.first = 0; - m_FaRange.second = 1; + this->m_RandGen = itk::Statistics::MersenneTwisterRandomVariateGenerator::New(); + this->m_RandGen->SetSeed(); + m_AdcRange.first = 0; + m_AdcRange.second = 0.004; + m_FaRange.first = 0; + m_FaRange.second = 1; } template< class ScalarType > RawShModel< ScalarType >::~RawShModel() { } template< class ScalarType > void RawShModel< ScalarType >::Clear() { - m_ShCoefficients.clear(); - m_PrototypeMaxDirection.clear(); - m_B0Signal.clear(); + m_ShCoefficients.clear(); + m_PrototypeMaxDirection.clear(); + m_B0Signal.clear(); } template< class ScalarType > void RawShModel< ScalarType >::RandomModel() { - m_ModelIndex = this->m_RandGen->GetIntegerVariate(m_B0Signal.size()-1); + m_ModelIndex = this->m_RandGen->GetIntegerVariate(m_B0Signal.size()-1); } template< class ScalarType > unsigned int RawShModel< ScalarType >::GetNumberOfKernels() { - return m_B0Signal.size(); + return m_B0Signal.size(); } template< class ScalarType > bool RawShModel< ScalarType >::SampleKernels(Image::Pointer diffImg, ItkUcharImageType::Pointer maskImage, TensorImageType::Pointer tensorImage, QballFilterType::CoefficientImageType::Pointer itkFeatureImage, ItkDoubleImageType::Pointer adcImage) { - if (diffImg.IsNull()) - return false; + if (diffImg.IsNull()) + return false; - typedef itk::VectorImage ITKDiffusionImageType; - ITKDiffusionImageType::Pointer itkVectorImagePointer = ITKDiffusionImageType::New(); + typedef itk::VectorImage ITKDiffusionImageType; + ITKDiffusionImageType::Pointer itkVectorImagePointer = ITKDiffusionImageType::New(); mitk::CastToItkImage(diffImg, itkVectorImagePointer); - const int shOrder = 2; - - if (tensorImage.IsNull()) + const int shOrder = 2; + + if (tensorImage.IsNull()) + { + typedef itk::DiffusionTensor3DReconstructionImageFilter< short, short, double > TensorReconstructionImageFilterType; + + TensorReconstructionImageFilterType::Pointer filter = TensorReconstructionImageFilterType::New(); + filter->SetBValue(static_cast( diffImg->GetProperty(mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str()).GetPointer() )->GetValue()); + filter->SetGradientImage( static_cast( diffImg->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() )->GetGradientDirectionsContainer(), itkVectorImagePointer ); + filter->Update(); + tensorImage = filter->GetOutput(); + } + + const int NumCoeffs = (shOrder*shOrder + shOrder + 2)/2 + shOrder; + if (itkFeatureImage.IsNull()) + { + QballFilterType::Pointer qballfilter = QballFilterType::New(); + qballfilter->SetBValue(static_cast( diffImg->GetProperty(mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str()).GetPointer() )->GetValue()); + qballfilter->SetGradientImage( static_cast( diffImg->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() )->GetGradientDirectionsContainer(), itkVectorImagePointer ); + qballfilter->SetLambda(0.006); + qballfilter->SetNormalizationMethod(QballFilterType::QBAR_RAW_SIGNAL); + qballfilter->Update(); + itkFeatureImage = qballfilter->GetCoefficientImage(); + } + + if (adcImage.IsNull()) + { + itk::AdcImageFilter< short, double >::Pointer adcFilter = itk::AdcImageFilter< short, double >::New(); + adcFilter->SetInput(itkVectorImagePointer); + adcFilter->SetGradientDirections(static_cast( diffImg->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() )->GetGradientDirectionsContainer()); + adcFilter->SetB_value(static_cast( diffImg->GetProperty(mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str()).GetPointer() )->GetValue()); + adcFilter->Update(); + adcImage = adcFilter->GetOutput(); + } + + int b0Index = 0; + for (unsigned int i=0; i( diffImg->GetProperty(mitk::DiffusionPropertyHelper::ORIGINALGRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() )->GetGradientDirectionsContainer()->Size(); i++) + if ( static_cast( diffImg->GetProperty(mitk::DiffusionPropertyHelper::ORIGINALGRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() )->GetGradientDirectionsContainer()->GetElement(i).magnitude()<0.001 ) { - typedef itk::DiffusionTensor3DReconstructionImageFilter< short, short, double > TensorReconstructionImageFilterType; - - TensorReconstructionImageFilterType::Pointer filter = TensorReconstructionImageFilterType::New(); - filter->SetBValue(static_cast( diffImg->GetProperty(mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str()).GetPointer() )->GetValue()); - filter->SetGradientImage( static_cast( diffImg->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() )->GetGradientDirectionsContainer(), itkVectorImagePointer ); - filter->Update(); - tensorImage = filter->GetOutput(); + b0Index = i; + break; } - const int NumCoeffs = (shOrder*shOrder + shOrder + 2)/2 + shOrder; - if (itkFeatureImage.IsNull()) + double max = 0; + { + itk::ImageRegionIterator< itk::VectorImage< short, 3 > > it(itkVectorImagePointer, itkVectorImagePointer->GetLargestPossibleRegion()); + while(!it.IsAtEnd()) { - QballFilterType::Pointer qballfilter = QballFilterType::New(); - qballfilter->SetBValue(static_cast( diffImg->GetProperty(mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str()).GetPointer() )->GetValue()); - qballfilter->SetGradientImage( static_cast( diffImg->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() )->GetGradientDirectionsContainer(), itkVectorImagePointer ); - qballfilter->SetLambda(0.006); - qballfilter->SetNormalizationMethod(QballFilterType::QBAR_RAW_SIGNAL); - qballfilter->Update(); - itkFeatureImage = qballfilter->GetCoefficientImage(); + if (maskImage.IsNotNull() && maskImage->GetPixel(it.GetIndex())<=0) + { + ++it; + continue; + } + if (it.Get()[b0Index]>max) + max = it.Get()[b0Index]; + ++it; } - - if (adcImage.IsNull()) + } + + MITK_INFO << "Sampling signal kernels."; + unsigned int count = 0; + itk::ImageRegionIterator< itk::Image< itk::DiffusionTensor3D< double >, 3 > > it(tensorImage, tensorImage->GetLargestPossibleRegion()); + while(!it.IsAtEnd()) + { + bool skipPixel = false; + if (maskImage.IsNotNull() && maskImage->GetPixel(it.GetIndex())<=0) { - itk::AdcImageFilter< short, double >::Pointer adcFilter = itk::AdcImageFilter< short, double >::New(); - adcFilter->SetInput(itkVectorImagePointer); - adcFilter->SetGradientDirections(static_cast( diffImg->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() )->GetGradientDirectionsContainer()); - adcFilter->SetB_value(static_cast( diffImg->GetProperty(mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str()).GetPointer() )->GetValue()); - adcFilter->Update(); - adcImage = adcFilter->GetOutput(); + ++it; + continue; } - - int b0Index = 0; - for (unsigned int i=0; i( diffImg->GetProperty(mitk::DiffusionPropertyHelper::ORIGINALGRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() )->GetGradientDirectionsContainer()->Size(); i++) - if ( static_cast( diffImg->GetProperty(mitk::DiffusionPropertyHelper::ORIGINALGRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() )->GetGradientDirectionsContainer()->GetElement(i).magnitude()<0.001 ) - { - b0Index = i; - break; - } - - double max = 0; + for (unsigned int i=0; i( diffImg->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() )->GetGradientDirectionsContainer()->Size(); i++) { - itk::ImageRegionIterator< itk::VectorImage< short, 3 > > it(itkVectorImagePointer, itkVectorImagePointer->GetLargestPossibleRegion()); - while(!it.IsAtEnd()) - { - if (maskImage.IsNotNull() && maskImage->GetPixel(it.GetIndex())<=0) - { - ++it; - continue; - } - if (it.Get()[b0Index]>max) - max = it.Get()[b0Index]; - ++it; - } + if (itkVectorImagePointer->GetPixel(it.GetIndex())[i]!=itkVectorImagePointer->GetPixel(it.GetIndex())[i] || itkVectorImagePointer->GetPixel(it.GetIndex())[i]<=0 || itkVectorImagePointer->GetPixel(it.GetIndex())[i]>itkVectorImagePointer->GetPixel(it.GetIndex())[b0Index]) + { + skipPixel = true; + break; + } } - - MITK_INFO << "Sampling signal kernels."; - unsigned int count = 0; - itk::ImageRegionIterator< itk::Image< itk::DiffusionTensor3D< double >, 3 > > it(tensorImage, tensorImage->GetLargestPossibleRegion()); - while(!it.IsAtEnd()) + if (skipPixel) { - bool skipPixel = false; - if (maskImage.IsNotNull() && maskImage->GetPixel(it.GetIndex())<=0) - { - ++it; - continue; - } - for (unsigned int i=0; i( diffImg->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() )->GetGradientDirectionsContainer()->Size(); i++) - { - if (itkVectorImagePointer->GetPixel(it.GetIndex())[i]!=itkVectorImagePointer->GetPixel(it.GetIndex())[i] || itkVectorImagePointer->GetPixel(it.GetIndex())[i]<=0 || itkVectorImagePointer->GetPixel(it.GetIndex())[i]>itkVectorImagePointer->GetPixel(it.GetIndex())[b0Index]) - { - skipPixel = true; - break; - } - } - if (skipPixel) - { - ++it; - continue; - } + ++it; + continue; + } - typedef itk::DiffusionTensor3D TensorType; - TensorType::EigenValuesArrayType eigenvalues; - TensorType::EigenVectorsMatrixType eigenvectors; - TensorType tensor = it.Get(); - double FA = tensor.GetFractionalAnisotropy(); - double ADC = adcImage->GetPixel(it.GetIndex()); - QballFilterType::CoefficientImageType::PixelType itkv = itkFeatureImage->GetPixel(it.GetIndex()); - vnl_vector_fixed< double, NumCoeffs > coeffs; - for (unsigned int c=0; cGetMaxNumKernels()>this->GetNumberOfKernels() && FA>m_FaRange.first && FAm_AdcRange.first && ADCSetShCoefficients( coeffs, (double)itkVectorImagePointer->GetPixel(it.GetIndex())[b0Index]/max )) - { - tensor.ComputeEigenAnalysis(eigenvalues, eigenvectors); - itk::Vector dir; - dir[0] = eigenvectors(2, 0); - dir[1] = eigenvectors(2, 1); - dir[2] = eigenvectors(2, 2); - m_PrototypeMaxDirection.push_back(dir); - count++; - MITK_INFO << "KERNEL: " << it.GetIndex() << " (" << count << ")"; - } - } - ++it; + typedef itk::DiffusionTensor3D TensorType; + TensorType::EigenValuesArrayType eigenvalues; + TensorType::EigenVectorsMatrixType eigenvectors; + TensorType tensor = it.Get(); + double FA = tensor.GetFractionalAnisotropy(); + double ADC = adcImage->GetPixel(it.GetIndex()); + QballFilterType::CoefficientImageType::PixelType itkv = itkFeatureImage->GetPixel(it.GetIndex()); + vnl_vector_fixed< double, NumCoeffs > coeffs; + for (unsigned int c=0; cGetMaxNumKernels()>this->GetNumberOfKernels() && FA>m_FaRange.first && FAm_AdcRange.first && ADCSetShCoefficients( coeffs, (double)itkVectorImagePointer->GetPixel(it.GetIndex())[b0Index]/max )) + { + tensor.ComputeEigenAnalysis(eigenvalues, eigenvectors); + itk::Vector dir; + dir[0] = eigenvectors(2, 0); + dir[1] = eigenvectors(2, 1); + dir[2] = eigenvectors(2, 2); + m_PrototypeMaxDirection.push_back(dir); + count++; + MITK_INFO << "KERNEL: " << it.GetIndex() << " (" << count << ")"; + } } + ++it; + } - if (m_ShCoefficients.size()>0) - return true; - return false; + if (m_ShCoefficients.size()>0) + return true; + return false; } // convert cartesian to spherical coordinates template< class ScalarType > -void RawShModel< ScalarType >::Cart2Sph( GradientListType gradients ) +vnl_matrix RawShModel::Cart2Sph(GradientListType &gradients ) { - m_SphCoords.set_size(gradients.size(), 2); - m_SphCoords.fill(0.0); - - for (unsigned int i=0; i sphCoords; + sphCoords.set_size(gradients.size(), 2); + sphCoords.fill(0.0); + + for (unsigned int i=0; i 0.0001 ) { - GradientType g = gradients[i]; - if( g.GetNorm() > 0.0001 ) - { - gradients[i].Normalize(); - m_SphCoords(i,0) = acos(gradients[i][2]); // theta - m_SphCoords(i,1) = atan2(gradients[i][1], gradients[i][0]); // phi - } + gradients[i].Normalize(); + sphCoords(i,0) = acos(gradients[i][2]); // theta + sphCoords(i,1) = atan2(gradients[i][1], gradients[i][0]); // phi } + } + + return sphCoords; } template< class ScalarType > -void RawShModel< ScalarType >::SetFiberDirection(GradientType fiberDirection) +vnl_matrix RawShModel< ScalarType >::SetFiberDirection(GradientType& fiberDirection) { - this->m_FiberDirection = fiberDirection; - this->m_FiberDirection.Normalize(); - - RandomModel(); - GradientType axis = itk::CrossProduct(this->m_FiberDirection, m_PrototypeMaxDirection.at(m_ModelIndex)); - axis.Normalize(); - - vnl_quaternion rotation(axis.GetVnlVector(), acos(dot_product(this->m_FiberDirection.GetVnlVector(), m_PrototypeMaxDirection.at(m_ModelIndex).GetVnlVector()))); - rotation.normalize(); - - GradientListType gradients; - for (unsigned int i=0; im_GradientList.size(); i++) + RandomModel(); + GradientType axis = itk::CrossProduct(fiberDirection, m_PrototypeMaxDirection.at(m_ModelIndex)); + axis.Normalize(); + + vnl_quaternion rotation(axis.GetVnlVector(), acos(dot_product(fiberDirection.GetVnlVector(), m_PrototypeMaxDirection.at(m_ModelIndex).GetVnlVector()))); + rotation.normalize(); + + GradientListType gradients; + for (unsigned int i=0; im_GradientList.size(); i++) + { + GradientType dir = this->m_GradientList.at(i); + if( dir.GetNorm() > 0.0001 ) { - GradientType dir = this->m_GradientList.at(i); - if( dir.GetNorm() > 0.0001 ) - { - dir.Normalize(); - vnl_vector_fixed< double, 3 > vnlDir = rotation.rotate(dir.GetVnlVector()); - dir[0] = vnlDir[0]; dir[1] = vnlDir[1]; dir[2] = vnlDir[2]; - dir.Normalize(); - } - gradients.push_back(dir); + dir.Normalize(); + vnl_vector_fixed< double, 3 > vnlDir = rotation.rotate(dir.GetVnlVector()); + dir[0] = vnlDir[0]; dir[1] = vnlDir[1]; dir[2] = vnlDir[2]; + dir.Normalize(); } + gradients.push_back(dir); + } - Cart2Sph( gradients ); + return Cart2Sph( gradients ); } template< class ScalarType > bool RawShModel< ScalarType >::SetShCoefficients(vnl_vector< double > shCoefficients, double b0 ) { - m_ShOrder = 2; - while ( (m_ShOrder*m_ShOrder + m_ShOrder + 2)/2 + m_ShOrder <= shCoefficients.size() ) - m_ShOrder += 2; - m_ShOrder -= 2; - - m_ModelIndex = m_B0Signal.size(); - m_B0Signal.push_back(b0); - m_ShCoefficients.push_back(shCoefficients); - - // itk::OrientationDistributionFunction odf; - // GradientListType gradients; - // for (unsigned int i=0; im_GradientList ); - // PixelType signal = SimulateMeasurement(); - - // int minDirIdx = 0; - // double min = itk::NumericTraits::max(); - // for (unsigned int i=0; ib0 || signal[i]<0) - // { - // MITK_INFO << "Corrupted signal value detected. Kernel rejected."; - // m_B0Signal.pop_back(); - // m_ShCoefficients.pop_back(); - // return false; - // } - // if (signal[i]m_GradientList.at(minDirIdx); - // maxDir.Normalize(); - // m_PrototypeMaxDirection.push_back(maxDir); - - Cart2Sph( this->m_GradientList ); - m_ModelIndex = -1; - - return true; + m_ShOrder = 2; + while ( (m_ShOrder*m_ShOrder + m_ShOrder + 2)/2 + m_ShOrder <= shCoefficients.size() ) + m_ShOrder += 2; + m_ShOrder -= 2; + + m_ModelIndex = m_B0Signal.size(); + m_B0Signal.push_back(b0); + m_ShCoefficients.push_back(shCoefficients); + + // itk::OrientationDistributionFunction odf; + // GradientListType gradients; + // for (unsigned int i=0; im_GradientList ); + // PixelType signal = SimulateMeasurement(); + + // int minDirIdx = 0; + // double min = itk::NumericTraits::max(); + // for (unsigned int i=0; ib0 || signal[i]<0) + // { + // MITK_INFO << "Corrupted signal value detected. Kernel rejected."; + // m_B0Signal.pop_back(); + // m_ShCoefficients.pop_back(); + // return false; + // } + // if (signal[i]m_GradientList.at(minDirIdx); + // maxDir.Normalize(); + // m_PrototypeMaxDirection.push_back(maxDir); + + Cart2Sph( this->m_GradientList ); + m_ModelIndex = -1; + + return true; } template< class ScalarType > -ScalarType RawShModel< ScalarType >::SimulateMeasurement(unsigned int dir) +ScalarType RawShModel< ScalarType >::SimulateMeasurement(unsigned int dir, GradientType& fiberDirection) { - ScalarType signal = 0; + ScalarType signal = 0; - if (m_ModelIndex==-1) - RandomModel(); + if (m_ModelIndex==-1) + RandomModel(); - if (dir>=this->m_GradientList.size()) - return signal; + vnl_matrix sphCoords = this->SetFiberDirection(fiberDirection); - int j, m; double mag, plm; + if (dir>=this->m_GradientList.size()) + return signal; - if (this->m_GradientList[dir].GetNorm()>0.001) - { - j=0; - for (int l=0; l<=static_cast(m_ShOrder); l=l+2) - for (m=-l; m<=l; m++) - { - plm = legendre_p(l,abs(m),cos(m_SphCoords(dir,0))); - mag = sqrt((double)(2*l+1)/(4.0*M_PI)*factorial(l-abs(m))/factorial(l+abs(m)))*plm; - double basis; - - if (m<0) - basis = sqrt(2.0)*mag*cos(fabs((double)m)*m_SphCoords(dir,1)); - else if (m==0) - basis = mag; - else - basis = pow(-1.0, m)*sqrt(2.0)*mag*sin(m*m_SphCoords(dir,1)); - - signal += m_ShCoefficients.at(m_ModelIndex)[j]*basis; - j++; - } - } - else - signal = m_B0Signal.at(m_ModelIndex); + int j, m; double mag, plm; + + if (this->m_GradientList[dir].GetNorm()>0.001) + { + j=0; + for (int l=0; l<=static_cast(m_ShOrder); l=l+2) + for (m=-l; m<=l; m++) + { + plm = legendre_p(l,abs(m),cos(sphCoords(dir,0))); + mag = sqrt((double)(2*l+1)/(4.0*M_PI)*factorial(l-abs(m))/factorial(l+abs(m)))*plm; + double basis; + + if (m<0) + basis = sqrt(2.0)*mag*cos(fabs((double)m)*sphCoords(dir,1)); + else if (m==0) + basis = mag; + else + basis = pow(-1.0, m)*sqrt(2.0)*mag*sin(m*sphCoords(dir,1)); - m_ModelIndex = -1; + signal += m_ShCoefficients.at(m_ModelIndex)[j]*basis; + j++; + } + } + else + signal = m_B0Signal.at(m_ModelIndex); - return signal; + m_ModelIndex = -1; + + return signal; } template< class ScalarType > -typename RawShModel< ScalarType >::PixelType RawShModel< ScalarType >::SimulateMeasurement() +typename RawShModel< ScalarType >::PixelType RawShModel< ScalarType >::SimulateMeasurement(GradientType& fiberDirection) { - if (m_ModelIndex==-1) - RandomModel(); + if (m_ModelIndex==-1) + RandomModel(); + + PixelType signal; + signal.SetSize(this->m_GradientList.size()); - PixelType signal; - signal.SetSize(this->m_GradientList.size()); + int M = this->m_GradientList.size(); + int j, m; double mag, plm; - int M = this->m_GradientList.size(); - int j, m; double mag, plm; + vnl_matrix< double > shBasis; + shBasis.set_size(M, m_ShCoefficients.at(m_ModelIndex).size()); + shBasis.fill(0.0); - vnl_matrix< double > shBasis; - shBasis.set_size(M, m_ShCoefficients.at(m_ModelIndex).size()); - shBasis.fill(0.0); + vnl_matrix sphCoords = this->SetFiberDirection(fiberDirection); - for (int p=0; pm_GradientList[p].GetNorm()>0.001) { - if (this->m_GradientList[p].GetNorm()>0.001) + j=0; + for (int l=0; l<=static_cast(m_ShOrder); l=l+2) + for (m=-l; m<=l; m++) { - j=0; - for (int l=0; l<=static_cast(m_ShOrder); l=l+2) - for (m=-l; m<=l; m++) - { - plm = legendre_p(l,abs(m),cos(m_SphCoords(p,0))); - mag = sqrt((double)(2*l+1)/(4.0*M_PI)*factorial(l-abs(m))/factorial(l+abs(m)))*plm; - - if (m<0) - shBasis(p,j) = sqrt(2.0)*mag*cos(fabs((double)m)*m_SphCoords(p,1)); - else if (m==0) - shBasis(p,j) = mag; - else - shBasis(p,j) = pow(-1.0, m)*sqrt(2.0)*mag*sin(m*m_SphCoords(p,1)); - - j++; - } - - - double val = 0; - for (unsigned int cidx=0; cidx(l,abs(m),cos(sphCoords(p,0))); + mag = sqrt((double)(2*l+1)/(4.0*M_PI)*factorial(l-abs(m))/factorial(l+abs(m)))*plm; + + if (m<0) + shBasis(p,j) = sqrt(2.0)*mag*cos(fabs((double)m)*sphCoords(p,1)); + else if (m==0) + shBasis(p,j) = mag; + else + shBasis(p,j) = pow(-1.0, m)*sqrt(2.0)*mag*sin(m*sphCoords(p,1)); + + j++; } - else - signal[p] = m_B0Signal.at(m_ModelIndex); + + + double val = 0; + for (unsigned int cidx=0; cidx #include #include namespace mitk { /** * \brief The spherical harmonic representation of a prototype diffusion weighted MR signal is used to obtain the direction dependent signal. * */ template< class ScalarType = double > class RawShModel : public DiffusionSignalModel< ScalarType > { public: RawShModel(); template< class OtherType >RawShModel(RawShModel* model) { this->m_CompartmentId = model->m_CompartmentId; this->m_T2 = model->GetT2(); this->m_FiberDirection = model->GetFiberDirection(); this->m_GradientList = model->GetGradientList(); this->m_VolumeFractionImage = model->GetVolumeFractionImage(); this->m_RandGen = model->GetRandomGenerator(); this->m_AdcRange = model->GetAdcRange(); this->m_FaRange = model->GetFaRange(); this->m_ShCoefficients = model->GetShCoefficients(); this->m_B0Signal = model->GetB0Signals(); this->m_SphCoords = model->GetSphericalCoordinates(); this->m_ShOrder = model->GetShOrder(); this->m_ModelIndex = model->GetModelIndex(); this->m_MaxNumKernels = model->GetMaxNumKernels(); } ~RawShModel(); typedef itk::Image< double, 3 > ItkDoubleImageType; typedef itk::Image< unsigned char, 3 > ItkUcharImageType; typedef itk::Image< itk::DiffusionTensor3D< double >, 3 > TensorImageType; typedef itk::AnalyticalDiffusionQballReconstructionImageFilter QballFilterType; typedef typename DiffusionSignalModel< ScalarType >::PixelType PixelType; typedef typename DiffusionSignalModel< ScalarType >::GradientType GradientType; typedef typename DiffusionSignalModel< ScalarType >::GradientListType GradientListType; typedef itk::Matrix< double, 3, 3 > MatrixType; /** Actual signal generation **/ - PixelType SimulateMeasurement(); - ScalarType SimulateMeasurement(unsigned int dir); + PixelType SimulateMeasurement(GradientType& fiberDirection); + ScalarType SimulateMeasurement(unsigned int dir, GradientType& fiberDirection); bool SetShCoefficients(vnl_vector< double > shCoefficients, double b0); - void SetFiberDirection(GradientType fiberDirection); + vnl_matrix SetFiberDirection(GradientType& fiberDirection); void SetFaRange(double min, double max){ m_FaRange.first = min; m_FaRange.second = max; } void SetAdcRange(double min, double max){ m_AdcRange.first = min; m_AdcRange.second = max; } void SetMaxNumKernels(unsigned int max){ m_MaxNumKernels = max; } unsigned int GetNumberOfKernels(); std::pair< double, double > GetFaRange(){ return m_FaRange; } std::pair< double, double > GetAdcRange(){ return m_AdcRange; } unsigned int GetMaxNumKernels(){ return m_MaxNumKernels; } void Clear(); std::vector< vnl_vector< double > > GetShCoefficients(){ return m_ShCoefficients; } std::vector< double > GetB0Signals(){ return m_B0Signal; } - vnl_matrix GetSphericalCoordinates(){ return m_SphCoords; } unsigned int GetShOrder(){ return m_ShOrder; } int GetModelIndex(){ return m_ModelIndex; } double GetBaselineSignal(int index){ return m_B0Signal.at(index); } vnl_vector< double > GetCoefficients(int listIndex){ return m_ShCoefficients.at(listIndex); } bool SampleKernels(Image::Pointer diffImg, ItkUcharImageType::Pointer maskImage, TensorImageType::Pointer tensorImage=nullptr, QballFilterType::CoefficientImageType::Pointer coeffImage=nullptr, ItkDoubleImageType::Pointer adcImage=nullptr); protected: - void Cart2Sph( GradientListType gradients ); + vnl_matrix Cart2Sph( GradientListType& gradients ); void RandomModel(); std::vector< vnl_vector< double > > m_ShCoefficients; std::vector< double > m_B0Signal; std::vector< GradientType > m_PrototypeMaxDirection; - vnl_matrix m_SphCoords; std::pair< double, double > m_AdcRange; std::pair< double, double > m_FaRange; unsigned int m_ShOrder; int m_ModelIndex; unsigned int m_MaxNumKernels; }; } #include "mitkRawShModel.cpp" #endif diff --git a/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkStickModel.cpp b/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkStickModel.cpp index e2e3ff6d2d..d41a052a9f 100644 --- a/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkStickModel.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkStickModel.cpp @@ -1,78 +1,75 @@ /*=================================================================== 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 #include #include using namespace mitk; template< class ScalarType > StickModel< ScalarType >::StickModel() - : m_Diffusivity(0.001) - , m_BValue(1000) + : m_Diffusivity(0.001) + , m_BValue(1000) { } template< class ScalarType > StickModel< ScalarType >::~StickModel() { } template< class ScalarType > -ScalarType StickModel< ScalarType >::SimulateMeasurement(unsigned int dir) +ScalarType StickModel< ScalarType >::SimulateMeasurement(unsigned int dir, GradientType& fiberDirection) { - ScalarType signal = 0; + ScalarType signal = 0; - if (dir>=this->m_GradientList.size()) - return signal; - - this->m_FiberDirection.Normalize(); + if (dir>=this->m_GradientList.size()) + return signal; - GradientType g = this->m_GradientList[dir]; - if (g.GetNorm()>0.0001) - { - ScalarType dot = this->m_FiberDirection*g; - signal = std::exp( -m_BValue*m_Diffusivity*dot*dot ); // skip * bVal becaus bVal is already encoded in the dot product (norm of g encodes b-value relative to baseline b-value m_BValue) - } - else - signal = 1; + GradientType g = this->m_GradientList[dir]; + if (g.GetNorm()>0.0001) + { + ScalarType dot = fiberDirection*g; + signal = std::exp( -m_BValue*m_Diffusivity*dot*dot ); // skip * bVal becaus bVal is already encoded in the dot product (norm of g encodes b-value relative to baseline b-value m_BValue) + } + else + signal = 1; - return signal; + return signal; } template< class ScalarType > -typename StickModel< ScalarType >::PixelType StickModel< ScalarType >::SimulateMeasurement() +typename StickModel< ScalarType >::PixelType StickModel< ScalarType >::SimulateMeasurement(GradientType &fiberDirection) { - this->m_FiberDirection.Normalize(); - PixelType signal; - signal.SetSize(this->m_GradientList.size()); + PixelType signal; + signal.SetSize(this->m_GradientList.size()); - for( unsigned int i=0; im_GradientList.size(); i++) + for( unsigned int i=0; im_GradientList.size(); i++) + { + GradientType g = this->m_GradientList[i]; + if (g.GetNorm()>0.0001) { - GradientType g = this->m_GradientList[i]; - if (g.GetNorm()>0.0001) - { - ScalarType dot = this->m_FiberDirection*g; - signal[i] = std::exp( -m_BValue*m_Diffusivity*dot*dot ); // skip * bVal becaus bVal is already encoded in the dot product (norm of g encodes b-value relative to baseline b-value m_BValue) - } - else - signal[i] = 1; + ScalarType dot = fiberDirection*g; + signal[i] = std::exp( -m_BValue*m_Diffusivity*dot*dot ); // skip * bVal becaus bVal is already encoded in the dot product (norm of g encodes b-value relative to baseline b-value m_BValue) } + else + signal[i] = 1; + } - return signal; + return signal; } diff --git a/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkStickModel.h b/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkStickModel.h index bc61237f84..5d327bb28b 100644 --- a/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkStickModel.h +++ b/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkStickModel.h @@ -1,75 +1,72 @@ /*=================================================================== 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_StickModel_H #define _MITK_StickModel_H #include namespace mitk { /** * \brief Generates the diffusion signal using an idealised cylinder with zero radius: e^(-bd(ng)²) * */ template< class ScalarType = double > class StickModel : public DiffusionSignalModel< ScalarType > { public: StickModel(); template< class OtherType >StickModel(StickModel* model) { this->m_CompartmentId = model->m_CompartmentId; this->m_T2 = model->GetT2(); - this->m_FiberDirection = model->GetFiberDirection(); this->m_GradientList = model->GetGradientList(); this->m_VolumeFractionImage = model->GetVolumeFractionImage(); this->m_RandGen = model->GetRandomGenerator(); this->m_BValue = model->GetBvalue(); this->m_Diffusivity = model->GetDiffusivity(); } ~StickModel(); typedef typename DiffusionSignalModel< ScalarType >::PixelType PixelType; typedef typename DiffusionSignalModel< ScalarType >::GradientType GradientType; typedef typename DiffusionSignalModel< ScalarType >::GradientListType GradientListType; /** Actual signal generation **/ - PixelType SimulateMeasurement(); - ScalarType SimulateMeasurement(unsigned int dir); + PixelType SimulateMeasurement(GradientType& fiberDirection); + ScalarType SimulateMeasurement(unsigned int dir, GradientType& fiberDirection); void SetBvalue(double bValue) { m_BValue = bValue; } ///< b-value used to generate the artificial signal double GetBvalue() { return m_BValue; } void SetDiffusivity(double diffusivity) { m_Diffusivity = diffusivity; } ///< Scalar diffusion constant double GetDiffusivity() { return m_Diffusivity; } - void SetFiberDirection(GradientType fiberDirection){ this->m_FiberDirection = fiberDirection; } - protected: double m_Diffusivity; ///< Scalar diffusion constant double m_BValue; ///< b-value used to generate the artificial signal }; } #include "mitkStickModel.cpp" #endif diff --git a/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkTensorModel.cpp b/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkTensorModel.cpp index b86e5e6886..89d5a49871 100644 --- a/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkTensorModel.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkTensorModel.cpp @@ -1,127 +1,125 @@ /*=================================================================== 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 #include #include using namespace mitk; template< class ScalarType > TensorModel< ScalarType >::TensorModel() : m_BValue(1000) { m_KernelDirection[0]=1; m_KernelDirection[1]=0; m_KernelDirection[2]=0; m_KernelTensorMatrix.fill(0.0); m_KernelTensorMatrix[0][0] = 0.002; m_KernelTensorMatrix[1][1] = 0.0005; m_KernelTensorMatrix[2][2] = 0.0005; } template< class ScalarType > TensorModel< ScalarType >::~TensorModel() { } template< class ScalarType > -ScalarType TensorModel< ScalarType >::SimulateMeasurement(unsigned int dir) +ScalarType TensorModel< ScalarType >::SimulateMeasurement(unsigned int dir, GradientType &fiberDirection) { ScalarType signal = 0; if (dir>=this->m_GradientList.size()) return signal; ItkTensorType tensor; tensor.Fill(0.0); - this->m_FiberDirection.Normalize(); - vnl_vector_fixed axis = itk::CrossProduct(m_KernelDirection, this->m_FiberDirection).GetVnlVector(); axis.normalize(); - vnl_quaternion rotation(axis, acos(m_KernelDirection*this->m_FiberDirection)); + vnl_vector_fixed axis = itk::CrossProduct(m_KernelDirection, fiberDirection).GetVnlVector(); axis.normalize(); + vnl_quaternion rotation(axis, acos(m_KernelDirection*fiberDirection)); rotation.normalize(); vnl_matrix_fixed matrix = rotation.rotation_matrix_transpose(); vnl_matrix_fixed tensorMatrix = matrix.transpose()*m_KernelTensorMatrix*matrix; tensor[0] = tensorMatrix[0][0]; tensor[1] = tensorMatrix[0][1]; tensor[2] = tensorMatrix[0][2]; tensor[3] = tensorMatrix[1][1]; tensor[4] = tensorMatrix[1][2]; tensor[5] = tensorMatrix[2][2]; GradientType g = this->m_GradientList[dir]; if (g.GetNorm()>0.0001) { // calc g^T * D * g itk::DiffusionTensor3D< ScalarType > S; S[0] = g[0]*g[0]; S[1] = g[1]*g[0]; S[2] = g[2]*g[0]; S[3] = g[1]*g[1]; S[4] = g[2]*g[1]; S[5] = g[2]*g[2]; ScalarType D_scalar = tensor[0]*S[0] + tensor[1]*S[1] + tensor[2]*S[2] + tensor[1]*S[1] + tensor[3]*S[3] + tensor[4]*S[4] + tensor[2]*S[2] + tensor[4]*S[4] + tensor[5]*S[5]; // check for corrupted tensor and generate signal if (D_scalar>=0) signal = std::exp ( -m_BValue * D_scalar ); // skip * bVal becaus bVal is already encoded in g^T*g (norm of g encodes b-value relative to baseline b-value m_BValue) } else signal = 1; return signal; } template< class ScalarType > -typename TensorModel< ScalarType >::PixelType TensorModel< ScalarType >::SimulateMeasurement() +typename TensorModel< ScalarType >::PixelType TensorModel< ScalarType >::SimulateMeasurement(GradientType& fiberDirection) { PixelType signal; signal.SetSize(this->m_GradientList.size()); signal.Fill(0.0); ItkTensorType tensor; tensor.Fill(0.0); - this->m_FiberDirection.Normalize(); - vnl_vector_fixed axis = itk::CrossProduct(m_KernelDirection, this->m_FiberDirection).GetVnlVector(); axis.normalize(); - vnl_quaternion rotation(axis, acos(m_KernelDirection*this->m_FiberDirection)); + vnl_vector_fixed axis = itk::CrossProduct(m_KernelDirection, fiberDirection).GetVnlVector(); axis.normalize(); + vnl_quaternion rotation(axis, acos(m_KernelDirection*fiberDirection)); rotation.normalize(); vnl_matrix_fixed matrix = rotation.rotation_matrix_transpose(); vnl_matrix_fixed tensorMatrix = matrix.transpose()*m_KernelTensorMatrix*matrix; tensor[0] = tensorMatrix[0][0]; tensor[1] = tensorMatrix[0][1]; tensor[2] = tensorMatrix[0][2]; tensor[3] = tensorMatrix[1][1]; tensor[4] = tensorMatrix[1][2]; tensor[5] = tensorMatrix[2][2]; for( unsigned int i=0; im_GradientList.size(); i++) { GradientType g = this->m_GradientList[i]; if (g.GetNorm()>0.0001) { // calc g^T * D * g itk::DiffusionTensor3D< ScalarType > S; S[0] = g[0]*g[0]; S[1] = g[1]*g[0]; S[2] = g[2]*g[0]; S[3] = g[1]*g[1]; S[4] = g[2]*g[1]; S[5] = g[2]*g[2]; ScalarType D_scalar = tensor[0]*S[0] + tensor[1]*S[1] + tensor[2]*S[2] + tensor[1]*S[1] + tensor[3]*S[3] + tensor[4]*S[4] + tensor[2]*S[2] + tensor[4]*S[4] + tensor[5]*S[5]; // check for corrupted tensor and generate signal if (D_scalar>=0) signal[i] = std::exp ( -m_BValue * D_scalar ); // skip * bVal becaus bVal is already encoded in g^T*g (norm of g encodes b-value relative to baseline b-value m_BValue) } else signal[i] = 1; } return signal; } diff --git a/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkTensorModel.h b/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkTensorModel.h index 4d9e73a953..69fec6e328 100644 --- a/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkTensorModel.h +++ b/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkTensorModel.h @@ -1,87 +1,86 @@ /*=================================================================== 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_TensorModel_H #define _MITK_TensorModel_H #include #include namespace mitk { /** * \brief Generates diffusion measurement employing a second rank tensor model: e^(-bg^TDg) * */ template< class ScalarType = double > class TensorModel : public DiffusionSignalModel< ScalarType > { public: TensorModel(); template< class OtherType >TensorModel(TensorModel* model) { this->m_CompartmentId = model->m_CompartmentId; this->m_T2 = model->GetT2(); this->m_FiberDirection = model->GetFiberDirection(); this->m_GradientList = model->GetGradientList(); this->m_VolumeFractionImage = model->GetVolumeFractionImage(); this->m_RandGen = model->GetRandomGenerator(); this->m_BValue = model->GetBvalue(); this->m_KernelDirection = model->GetKernelDirection(); this->m_KernelTensorMatrix = model->GetKernelTensorMatrix(); } ~TensorModel(); typedef typename DiffusionSignalModel< ScalarType >::PixelType PixelType; typedef itk::DiffusionTensor3D< ScalarType > ItkTensorType; typedef typename DiffusionSignalModel< ScalarType >::GradientType GradientType; typedef typename DiffusionSignalModel< ScalarType >::GradientListType GradientListType; /** Actual signal generation **/ - PixelType SimulateMeasurement(); - ScalarType SimulateMeasurement(unsigned int dir); + PixelType SimulateMeasurement(GradientType& fiberDirection); + ScalarType SimulateMeasurement(unsigned int dir, GradientType& fiberDirection); void SetBvalue(double bValue) { m_BValue = bValue; } ///< b-value used to generate the artificial signal double GetBvalue() { return m_BValue; } void SetDiffusivity1(double d1){ m_KernelTensorMatrix[0][0] = d1; } void SetDiffusivity2(double d2){ m_KernelTensorMatrix[1][1] = d2; } void SetDiffusivity3(double d3){ m_KernelTensorMatrix[2][2] = d3; } double GetDiffusivity1() { return m_KernelTensorMatrix[0][0]; } double GetDiffusivity2() { return m_KernelTensorMatrix[1][1]; } double GetDiffusivity3() { return m_KernelTensorMatrix[2][2]; } - void SetFiberDirection(GradientType fiberDirection){ this->m_FiberDirection = fiberDirection; } GradientType GetKernelDirection(){ return m_KernelDirection; } vnl_matrix_fixed GetKernelTensorMatrix(){ return m_KernelTensorMatrix; } protected: /** Calculates tensor matrix from FA and ADC **/ void UpdateKernelTensor(); GradientType m_KernelDirection; ///< Direction of the kernel tensors principal eigenvector vnl_matrix_fixed m_KernelTensorMatrix; ///< 3x3 matrix containing the kernel tensor values double m_BValue; ///< b-value used to generate the artificial signal }; } #include "mitkTensorModel.cpp" #endif diff --git a/Modules/DiffusionImaging/FiberTracking/Fiberfox/itkTractsToDWIImageFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Fiberfox/itkTractsToDWIImageFilter.cpp index 57db32ab0e..cf7c84852e 100755 --- a/Modules/DiffusionImaging/FiberTracking/Fiberfox/itkTractsToDWIImageFilter.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Fiberfox/itkTractsToDWIImageFilter.cpp @@ -1,1674 +1,1681 @@ /*=================================================================== 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 #include #include #include #include #include #include #include #include #include #include #include #include namespace itk { template< class PixelType > TractsToDWIImageFilter< PixelType >::TractsToDWIImageFilter() : m_FiberBundle(nullptr) , m_StatusText("") , m_UseConstantRandSeed(false) , m_RandGen(itk::Statistics::MersenneTwisterRandomVariateGenerator::New()) { m_RandGen->SetSeed(); m_DoubleInterpolator = itk::LinearInterpolateImageFunction< ItkDoubleImgType, float >::New(); + m_NullDir.Fill(0); } template< class PixelType > TractsToDWIImageFilter< PixelType >::~TractsToDWIImageFilter() { } template< class PixelType > TractsToDWIImageFilter< PixelType >::DoubleDwiType::Pointer TractsToDWIImageFilter< PixelType >:: - SimulateKspaceAcquisition( std::vector< DoubleDwiType::Pointer >& images ) + SimulateKspaceAcquisition( std::vector< DoubleDwiType::Pointer >& compartment_images ) { unsigned int numFiberCompartments = m_Parameters.m_FiberModelList.size(); // create slice object ImageRegion<2> sliceRegion; sliceRegion.SetSize(0, m_WorkingImageRegion.GetSize()[0]); sliceRegion.SetSize(1, m_WorkingImageRegion.GetSize()[1]); Vector< double, 2 > sliceSpacing; sliceSpacing[0] = m_WorkingSpacing[0]; sliceSpacing[1] = m_WorkingSpacing[1]; - DoubleDwiType::PixelType nullPix; nullPix.SetSize(images.at(0)->GetVectorLength()); nullPix.Fill(0.0); + DoubleDwiType::PixelType nullPix; nullPix.SetSize(compartment_images.at(0)->GetVectorLength()); nullPix.Fill(0.0); auto magnitudeDwiImage = DoubleDwiType::New(); magnitudeDwiImage->SetSpacing( m_Parameters.m_SignalGen.m_ImageSpacing ); magnitudeDwiImage->SetOrigin( m_Parameters.m_SignalGen.m_ImageOrigin ); magnitudeDwiImage->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection ); magnitudeDwiImage->SetLargestPossibleRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); magnitudeDwiImage->SetBufferedRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); magnitudeDwiImage->SetRequestedRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); - magnitudeDwiImage->SetVectorLength( images.at(0)->GetVectorLength() ); + magnitudeDwiImage->SetVectorLength( compartment_images.at(0)->GetVectorLength() ); magnitudeDwiImage->Allocate(); magnitudeDwiImage->FillBuffer(nullPix); m_PhaseImage = DoubleDwiType::New(); m_PhaseImage->SetSpacing( m_Parameters.m_SignalGen.m_ImageSpacing ); m_PhaseImage->SetOrigin( m_Parameters.m_SignalGen.m_ImageOrigin ); m_PhaseImage->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection ); m_PhaseImage->SetLargestPossibleRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); m_PhaseImage->SetBufferedRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); m_PhaseImage->SetRequestedRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); - m_PhaseImage->SetVectorLength( images.at(0)->GetVectorLength() ); + m_PhaseImage->SetVectorLength( compartment_images.at(0)->GetVectorLength() ); m_PhaseImage->Allocate(); m_PhaseImage->FillBuffer(nullPix); m_KspaceImage = DoubleDwiType::New(); m_KspaceImage->SetSpacing( m_Parameters.m_SignalGen.m_ImageSpacing ); m_KspaceImage->SetOrigin( m_Parameters.m_SignalGen.m_ImageOrigin ); m_KspaceImage->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection ); m_KspaceImage->SetLargestPossibleRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); m_KspaceImage->SetBufferedRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); m_KspaceImage->SetRequestedRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); m_KspaceImage->SetVectorLength( m_Parameters.m_SignalGen.m_NumberOfCoils ); m_KspaceImage->Allocate(); m_KspaceImage->FillBuffer(nullPix); std::vector< unsigned int > spikeVolume; for (unsigned int i=0; iGetIntegerVariate()%(images.at(0)->GetVectorLength())); + spikeVolume.push_back(m_RandGen->GetIntegerVariate()%(compartment_images.at(0)->GetVectorLength())); std::sort (spikeVolume.begin(), spikeVolume.end()); std::reverse (spikeVolume.begin(), spikeVolume.end()); // calculate coil positions double a = m_Parameters.m_SignalGen.m_ImageRegion.GetSize(0)*m_Parameters.m_SignalGen.m_ImageSpacing[0]; double b = m_Parameters.m_SignalGen.m_ImageRegion.GetSize(1)*m_Parameters.m_SignalGen.m_ImageSpacing[1]; double c = m_Parameters.m_SignalGen.m_ImageRegion.GetSize(2)*m_Parameters.m_SignalGen.m_ImageSpacing[2]; double diagonal = sqrt(a*a+b*b)/1000; // image diagonal in m m_CoilPointset = mitk::PointSet::New(); std::vector< itk::Vector > coilPositions; itk::Vector pos; pos.Fill(0.0); pos[1] = -diagonal/2; itk::Vector center; center[0] = a/2-m_Parameters.m_SignalGen.m_ImageSpacing[0]/2; center[1] = b/2-m_Parameters.m_SignalGen.m_ImageSpacing[2]/2; center[2] = c/2-m_Parameters.m_SignalGen.m_ImageSpacing[1]/2; for (int c=0; cInsertPoint(c, pos*1000 + m_Parameters.m_SignalGen.m_ImageOrigin.GetVectorFromOrigin() + center ); double rz = 360.0/m_Parameters.m_SignalGen.m_NumberOfCoils * M_PI/180; 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]; pos.SetVnlVector(rotZ*pos.GetVnlVector()); } PrintToLog("0% 10 20 30 40 50 60 70 80 90 100%", false, true, false); PrintToLog("|----|----|----|----|----|----|----|----|----|----|\n*", false, false, false); unsigned long lastTick = 0; - boost::progress_display disp(images.at(0)->GetVectorLength()*images.at(0)->GetLargestPossibleRegion().GetSize(2)); - for (unsigned int g=0; gGetVectorLength(); g++) + boost::progress_display disp(compartment_images.at(0)->GetVectorLength()*compartment_images.at(0)->GetLargestPossibleRegion().GetSize(2)); + +#pragma omp parallel for + for (unsigned int g=0; gGetVectorLength(); g++) { + if (this->GetAbortGenerateData()) + continue; + std::vector< unsigned int > spikeSlice; +#pragma omp critical while (!spikeVolume.empty() && spikeVolume.back()==g) { - spikeSlice.push_back(m_RandGen->GetIntegerVariate()%images.at(0)->GetLargestPossibleRegion().GetSize(2)); + spikeSlice.push_back(m_RandGen->GetIntegerVariate()%compartment_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++) + for (unsigned int z=0; zGetLargestPossibleRegion().GetSize(2); z++) { - std::vector< SliceType::Pointer > compartmentSlices; + std::vector< Double2DImageType::Pointer > compartment_slices; std::vector< double > t2Vector; std::vector< double > t1Vector; - 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++) + 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; + Double2DImageType::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]); + slice->SetPixel(index2D, compartment_images.at(i)->GetPixel(index3D)[g]); } - compartmentSlices.push_back(slice); + compartment_slices.push_back(slice); t2Vector.push_back(signalModel->GetT2()); t1Vector.push_back(signalModel->GetT1()); } int numSpikes = 0; while (!spikeSlice.empty() && spikeSlice.back()==z) { numSpikes++; spikeSlice.pop_back(); } int spikeCoil = m_RandGen->GetIntegerVariate()%m_Parameters.m_SignalGen.m_NumberOfCoils; if (this->GetAbortGenerateData()) - return nullptr; + continue; -#pragma omp parallel for for (int c=0; c::New(); - idft->SetCompartmentImages(compartmentSlices); + auto idft = itk::KspaceImageFilter< Double2DImageType::PixelType >::New(); + idft->SetCompartmentImages(compartment_slices); idft->SetT2(t2Vector); idft->SetT1(t1Vector); idft->SetUseConstantRandSeed(m_UseConstantRandSeed); idft->SetParameters(&m_Parameters); - idft->SetZ((double)z-(double)( images.at(0)->GetLargestPossibleRegion().GetSize(2) - -images.at(0)->GetLargestPossibleRegion().GetSize(2)%2 ) / 2.0); + idft->SetZ((double)z-(double)( compartment_images.at(0)->GetLargestPossibleRegion().GetSize(2) + -compartment_images.at(0)->GetLargestPossibleRegion().GetSize(2)%2 ) / 2.0); idft->SetZidx(z); idft->SetCoilPosition(coilPositions.at(c)); idft->SetFiberBundle(m_FiberBundleWorkingCopy); idft->SetTranslation(m_Translations.at(g)); idft->SetRotation(m_Rotations.at(g)); idft->SetDiffusionGradientDirection(m_Parameters.m_SignalGen.GetGradientDirection(g)); if (c==spikeCoil) idft->SetSpikesPerSlice(numSpikes); idft->Update(); #pragma omp critical if (c==spikeCoil && numSpikes>0) { m_SpikeLog += "Volume " + boost::lexical_cast(g) + " Coil " + boost::lexical_cast(c) + "\n"; m_SpikeLog += idft->GetSpikeLog(); } - ComplexSliceType::Pointer fSlice; + Complex2DImageType::Pointer fSlice; fSlice = idft->GetOutput(); // fourier transform slice - ComplexSliceType::Pointer newSlice; - auto dft = itk::DftImageFilter< SliceType::PixelType >::New(); + Complex2DImageType::Pointer newSlice; + auto dft = itk::DftImageFilter< Double2DImageType::PixelType >::New(); dft->SetInput(fSlice); dft->SetParameters(m_Parameters); 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; - ComplexSliceType::IndexType index2D; index2D[0]=x; index2D[1]=y; + Complex2DImageType::IndexType index2D; index2D[0]=x; index2D[1]=y; - ComplexSliceType::PixelType cPix = newSlice->GetPixel(index2D); + Complex2DImageType::PixelType cPix = newSlice->GetPixel(index2D); double magn = sqrt(cPix.real()*cPix.real()+cPix.imag()*cPix.imag()); double phase = 0; if (cPix.real()!=0) phase = atan( cPix.imag()/cPix.real() ); DoubleDwiType::PixelType dwiPix = magnitudeDwiImage->GetPixel(index3D); DoubleDwiType::PixelType phasePix = m_PhaseImage->GetPixel(index3D); if (m_Parameters.m_SignalGen.m_NumberOfCoils>1) { dwiPix[g] += magn*magn; phasePix[g] += phase*phase; } else { dwiPix[g] = magn; phasePix[g] = phase; } -#pragma omp critical +//#pragma omp critical { magnitudeDwiImage->SetPixel(index3D, dwiPix); m_PhaseImage->SetPixel(index3D, phasePix); // k-space image if (g==0) { DoubleDwiType::PixelType kspacePix = m_KspaceImage->GetPixel(index3D); kspacePix[c] = idft->GetKSpaceImage()->GetPixel(index2D); m_KspaceImage->SetPixel(index3D, kspacePix); } } } } if (m_Parameters.m_SignalGen.m_NumberOfCoils>1) { -#ifdef WIN32 -#pragma omp parallel for -#else -#pragma omp parallel for collapse(2) -#endif for (int y=0; y(magnitudeDwiImage->GetLargestPossibleRegion().GetSize(1)); y++) for (int x=0; x(magnitudeDwiImage->GetLargestPossibleRegion().GetSize(0)); x++) { DoubleDwiType::IndexType index3D; index3D[0]=x; index3D[1]=y; index3D[2]=z; DoubleDwiType::PixelType magPix = magnitudeDwiImage->GetPixel(index3D); magPix[g] = sqrt(magPix[g]/m_Parameters.m_SignalGen.m_NumberOfCoils); DoubleDwiType::PixelType phasePix = m_PhaseImage->GetPixel(index3D); phasePix[g] = sqrt(phasePix[g]/m_Parameters.m_SignalGen.m_NumberOfCoils); -#pragma omp critical +//#pragma omp critical { magnitudeDwiImage->SetPixel(index3D, magPix); m_PhaseImage->SetPixel(index3D, phasePix); } } } ++disp; unsigned long newTick = 50*disp.count()/disp.expected_count(); for (unsigned long tick = 0; tick<(newTick-lastTick); tick++) PrintToLog("*", false, false, false); lastTick = newTick; } } + + PrintToLog("\n", false); return magnitudeDwiImage; } template< class PixelType > TractsToDWIImageFilter< PixelType >::ItkDoubleImgType::Pointer TractsToDWIImageFilter< PixelType >:: NormalizeInsideMask(ItkDoubleImgType::Pointer image) { double max = itk::NumericTraits< double >::min(); double min = itk::NumericTraits< double >::max(); itk::ImageRegionIterator< ItkDoubleImgType > it(image, image->GetLargestPossibleRegion()); while(!it.IsAtEnd()) { if (m_Parameters.m_SignalGen.m_MaskImage.IsNotNull() && m_Parameters.m_SignalGen.m_MaskImage->GetPixel(it.GetIndex())<=0) { it.Set(0.0); ++it; continue; } - // if (it.Get()>900) - // it.Set(900); if (it.Get()>max) max = it.Get(); if (it.Get()::New(); scaler->SetInput(image); scaler->SetShift(-min); scaler->SetScale(1.0/(max-min)); scaler->Update(); return scaler->GetOutput(); } template< class PixelType > void TractsToDWIImageFilter< PixelType >::CheckVolumeFractionImages() { m_UseRelativeNonFiberVolumeFractions = false; // check for fiber volume fraction maps unsigned int fibVolImages = 0; for (std::size_t i=0; iGetVolumeFractionImage().IsNotNull()) { PrintToLog("Using volume fraction map for fiber compartment " + boost::lexical_cast(i+1)); fibVolImages++; } } // check for non-fiber volume fraction maps unsigned int nonfibVolImages = 0; for (std::size_t i=0; iGetVolumeFractionImage().IsNotNull()) { PrintToLog("Using volume fraction map for non-fiber compartment " + boost::lexical_cast(i+1)); nonfibVolImages++; } } // not all fiber compartments are using volume fraction maps // --> non-fiber volume fractions are assumed to be relative to the // non-fiber volume and not absolute voxel-volume fractions. // this means if two non-fiber compartments are used but only one of them // has an associated volume fraction map, the repesctive other volume fraction map // can be determined as inverse (1-val) of the present volume fraction map- if ( fibVolImages::New(); inverter->SetMaximum(1.0); if ( m_Parameters.m_NonFiberModelList[0]->GetVolumeFractionImage().IsNull() && m_Parameters.m_NonFiberModelList[1]->GetVolumeFractionImage().IsNotNull() ) { // m_Parameters.m_NonFiberModelList[1]->SetVolumeFractionImage( // NormalizeInsideMask( m_Parameters.m_NonFiberModelList[1]->GetVolumeFractionImage() ) ); inverter->SetInput( m_Parameters.m_NonFiberModelList[1]->GetVolumeFractionImage() ); inverter->Update(); m_Parameters.m_NonFiberModelList[0]->SetVolumeFractionImage(inverter->GetOutput()); } else if ( m_Parameters.m_NonFiberModelList[1]->GetVolumeFractionImage().IsNull() && m_Parameters.m_NonFiberModelList[0]->GetVolumeFractionImage().IsNotNull() ) { // m_Parameters.m_NonFiberModelList[0]->SetVolumeFractionImage( // NormalizeInsideMask( m_Parameters.m_NonFiberModelList[0]->GetVolumeFractionImage() ) ); inverter->SetInput( m_Parameters.m_NonFiberModelList[0]->GetVolumeFractionImage() ); inverter->Update(); m_Parameters.m_NonFiberModelList[1]->SetVolumeFractionImage(inverter->GetOutput()); } else { itkExceptionMacro("Something went wrong in automatically calculating the missing non-fiber volume fraction image!" " Did you use two non fiber compartments but only one volume fraction image?" " Then it should work and this error is really strange."); } m_UseRelativeNonFiberVolumeFractions = true; nonfibVolImages++; } // Up to two fiber compartments are allowed without volume fraction maps since the volume fractions can then be determined automatically if (m_Parameters.m_FiberModelList.size()>2 && fibVolImages!=m_Parameters.m_FiberModelList.size()) itkExceptionMacro("More than two fiber compartment selected but no corresponding volume fraction maps set!"); // One non-fiber compartment is allowed without volume fraction map since the volume fraction can then be determined automatically if (m_Parameters.m_NonFiberModelList.size()>1 && nonfibVolImages!=m_Parameters.m_NonFiberModelList.size()) itkExceptionMacro("More than one non-fiber compartment selected but no volume fraction maps set!"); if (fibVolImages0) { PrintToLog("Not all fiber compartments are using an associated volume fraction image.\n" "Assuming non-fiber volume fraction images to contain values relative to the" " remaining non-fiber volume, not absolute values."); m_UseRelativeNonFiberVolumeFractions = true; // itk::ImageFileWriter::Pointer wr = itk::ImageFileWriter::New(); // wr->SetInput(m_Parameters.m_NonFiberModelList[1]->GetVolumeFractionImage()); // wr->SetFileName("/local/volumefraction.nrrd"); // wr->Update(); } // initialize the images that store the output volume fraction of each compartment m_VolumeFractions.clear(); for (std::size_t i=0; iSetSpacing( m_WorkingSpacing ); doubleImg->SetOrigin( m_WorkingOrigin ); doubleImg->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection ); doubleImg->SetLargestPossibleRegion( m_WorkingImageRegion ); doubleImg->SetBufferedRegion( m_WorkingImageRegion ); doubleImg->SetRequestedRegion( m_WorkingImageRegion ); doubleImg->Allocate(); doubleImg->FillBuffer(0); m_VolumeFractions.push_back(doubleImg); } } template< class PixelType > void TractsToDWIImageFilter< PixelType >::InitializeData() { m_Rotations.clear(); m_Translations.clear(); m_MotionLog = ""; m_SpikeLog = ""; // initialize output dwi image m_Parameters.m_SignalGen.m_CroppedRegion = m_Parameters.m_SignalGen.m_ImageRegion; m_Parameters.m_SignalGen.m_CroppedRegion.SetSize( 1, m_Parameters.m_SignalGen.m_CroppedRegion.GetSize(1) *m_Parameters.m_SignalGen.m_CroppingFactor); itk::Point shiftedOrigin = m_Parameters.m_SignalGen.m_ImageOrigin; shiftedOrigin[1] += (m_Parameters.m_SignalGen.m_ImageRegion.GetSize(1) -m_Parameters.m_SignalGen.m_CroppedRegion.GetSize(1))*m_Parameters.m_SignalGen.m_ImageSpacing[1]/2; m_OutputImage = OutputImageType::New(); m_OutputImage->SetSpacing( m_Parameters.m_SignalGen.m_ImageSpacing ); m_OutputImage->SetOrigin( shiftedOrigin ); m_OutputImage->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection ); m_OutputImage->SetLargestPossibleRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); m_OutputImage->SetBufferedRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); m_OutputImage->SetRequestedRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); m_OutputImage->SetVectorLength( m_Parameters.m_SignalGen.GetNumVolumes() ); m_OutputImage->Allocate(); typename OutputImageType::PixelType temp; temp.SetSize(m_Parameters.m_SignalGen.GetNumVolumes()); temp.Fill(0.0); m_OutputImage->FillBuffer(temp); - MITK_INFO << "Output image parameters:"; - MITK_INFO << "Size:" << m_Parameters.m_SignalGen.m_CroppedRegion; - MITK_INFO << "Spacing:" << m_Parameters.m_SignalGen.m_ImageSpacing; - MITK_INFO << "Origin:" << shiftedOrigin; - MITK_INFO << "Matrix:" << m_Parameters.m_SignalGen.m_ImageDirection; - // Apply in-plane upsampling for Gibbs ringing artifact double upsampling = 1; if (m_Parameters.m_SignalGen.m_DoAddGibbsRinging) upsampling = 2; m_WorkingSpacing = m_Parameters.m_SignalGen.m_ImageSpacing; m_WorkingSpacing[0] /= upsampling; m_WorkingSpacing[1] /= upsampling; m_WorkingImageRegion = m_Parameters.m_SignalGen.m_ImageRegion; m_WorkingImageRegion.SetSize(0, m_Parameters.m_SignalGen.m_ImageRegion.GetSize()[0]*upsampling); m_WorkingImageRegion.SetSize(1, m_Parameters.m_SignalGen.m_ImageRegion.GetSize()[1]*upsampling); m_WorkingOrigin = m_Parameters.m_SignalGen.m_ImageOrigin; m_WorkingOrigin[0] -= m_Parameters.m_SignalGen.m_ImageSpacing[0]/2; m_WorkingOrigin[0] += m_WorkingSpacing[0]/2; m_WorkingOrigin[1] -= m_Parameters.m_SignalGen.m_ImageSpacing[1]/2; m_WorkingOrigin[1] += m_WorkingSpacing[1]/2; m_WorkingOrigin[2] -= m_Parameters.m_SignalGen.m_ImageSpacing[2]/2; m_WorkingOrigin[2] += m_WorkingSpacing[2]/2; m_VoxelVolume = m_WorkingSpacing[0]*m_WorkingSpacing[1]*m_WorkingSpacing[2]; // generate double images to store the individual compartment signals m_CompartmentImages.clear(); int numFiberCompartments = m_Parameters.m_FiberModelList.size(); int numNonFiberCompartments = m_Parameters.m_NonFiberModelList.size(); for (int i=0; iSetSpacing( m_WorkingSpacing ); doubleDwi->SetOrigin( m_WorkingOrigin ); doubleDwi->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection ); doubleDwi->SetLargestPossibleRegion( m_WorkingImageRegion ); doubleDwi->SetBufferedRegion( m_WorkingImageRegion ); doubleDwi->SetRequestedRegion( m_WorkingImageRegion ); doubleDwi->SetVectorLength( m_Parameters.m_SignalGen.GetNumVolumes() ); doubleDwi->Allocate(); DoubleDwiType::PixelType pix; pix.SetSize(m_Parameters.m_SignalGen.GetNumVolumes()); pix.Fill(0.0); doubleDwi->FillBuffer(pix); m_CompartmentImages.push_back(doubleDwi); } if (m_FiberBundle.IsNull() && m_InputImage.IsNotNull()) { m_CompartmentImages.clear(); m_Parameters.m_SignalGen.m_DoAddMotion = false; m_Parameters.m_SignalGen.m_DoSimulateRelaxation = false; PrintToLog("Simulating acquisition for input diffusion-weighted image.", false); auto caster = itk::CastImageFilter< OutputImageType, DoubleDwiType >::New(); caster->SetInput(m_InputImage); caster->Update(); if (m_Parameters.m_SignalGen.m_DoAddGibbsRinging) { PrintToLog("Upsampling input diffusion-weighted image for Gibbs ringing simulation.", false); auto resampler = itk::ResampleDwiImageFilter< double >::New(); resampler->SetInput(caster->GetOutput()); itk::Vector< double, 3 > samplingFactor; samplingFactor[0] = upsampling; samplingFactor[1] = upsampling; samplingFactor[2] = 1; resampler->SetSamplingFactor(samplingFactor); resampler->SetInterpolation(itk::ResampleDwiImageFilter< double >::Interpolate_WindowedSinc); resampler->Update(); m_CompartmentImages.push_back(resampler->GetOutput()); } else m_CompartmentImages.push_back(caster->GetOutput()); for (unsigned int g=0; g::New(); rescaler->SetInput(0,m_Parameters.m_SignalGen.m_MaskImage); rescaler->SetOutputMaximum(100); rescaler->SetOutputMinimum(0); rescaler->Update(); // resample mask image auto resampler = itk::ResampleImageFilter::New(); resampler->SetInput(rescaler->GetOutput()); resampler->SetOutputParametersFromImage(m_Parameters.m_SignalGen.m_MaskImage); resampler->SetSize(m_WorkingImageRegion.GetSize()); resampler->SetOutputSpacing(m_WorkingSpacing); resampler->SetOutputOrigin(m_WorkingOrigin); auto nn_interpolator = itk::NearestNeighborInterpolateImageFunction::New(); resampler->SetInterpolator(nn_interpolator); resampler->Update(); m_Parameters.m_SignalGen.m_MaskImage = resampler->GetOutput(); } // resample frequency map if (m_Parameters.m_SignalGen.m_FrequencyMap.IsNotNull()) { auto resampler = itk::ResampleImageFilter::New(); resampler->SetInput(m_Parameters.m_SignalGen.m_FrequencyMap); resampler->SetOutputParametersFromImage(m_Parameters.m_SignalGen.m_FrequencyMap); resampler->SetSize(m_WorkingImageRegion.GetSize()); resampler->SetOutputSpacing(m_WorkingSpacing); resampler->SetOutputOrigin(m_WorkingOrigin); auto nn_interpolator = itk::NearestNeighborInterpolateImageFunction::New(); resampler->SetInterpolator(nn_interpolator); resampler->Update(); m_Parameters.m_SignalGen.m_FrequencyMap = resampler->GetOutput(); } } m_MaskImageSet = true; if (m_Parameters.m_SignalGen.m_MaskImage.IsNull()) { // no input tissue mask is set -> create default PrintToLog("No tissue mask set", false); m_Parameters.m_SignalGen.m_MaskImage = ItkUcharImgType::New(); m_Parameters.m_SignalGen.m_MaskImage->SetSpacing( m_WorkingSpacing ); m_Parameters.m_SignalGen.m_MaskImage->SetOrigin( m_WorkingOrigin ); m_Parameters.m_SignalGen.m_MaskImage->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection ); m_Parameters.m_SignalGen.m_MaskImage->SetLargestPossibleRegion( m_WorkingImageRegion ); m_Parameters.m_SignalGen.m_MaskImage->SetBufferedRegion( m_WorkingImageRegion ); m_Parameters.m_SignalGen.m_MaskImage->SetRequestedRegion( m_WorkingImageRegion ); m_Parameters.m_SignalGen.m_MaskImage->Allocate(); m_Parameters.m_SignalGen.m_MaskImage->FillBuffer(100); m_MaskImageSet = false; } else { if (m_Parameters.m_SignalGen.m_MaskImage->GetLargestPossibleRegion()!=m_WorkingImageRegion) { itkExceptionMacro("Mask image and specified DWI geometry are not matching!"); } PrintToLog("Using tissue mask", false); } if (m_Parameters.m_SignalGen.m_DoAddMotion) { if (m_Parameters.m_SignalGen.m_DoRandomizeMotion) { PrintToLog("Random motion artifacts:", false); PrintToLog("Maximum rotation: +/-" + boost::lexical_cast(m_Parameters.m_SignalGen.m_Rotation) + "°", false); PrintToLog("Maximum translation: +/-" + boost::lexical_cast(m_Parameters.m_SignalGen.m_Translation) + "mm", false); } else { PrintToLog("Linear motion artifacts:", false); PrintToLog("Maximum rotation: " + boost::lexical_cast(m_Parameters.m_SignalGen.m_Rotation) + "°", false); PrintToLog("Maximum translation: " + boost::lexical_cast(m_Parameters.m_SignalGen.m_Translation) + "mm", false); } } if ( m_Parameters.m_SignalGen.m_MotionVolumes.empty() ) { // no motion in first volume m_Parameters.m_SignalGen.m_MotionVolumes.push_back(false); // motion in all other volumes while ( m_Parameters.m_SignalGen.m_MotionVolumes.size() < m_Parameters.m_SignalGen.GetNumVolumes() ) { m_Parameters.m_SignalGen.m_MotionVolumes.push_back(true); } } // we need to know for every volume if there is motion. if this information is missing, then set corresponding fal to false while ( m_Parameters.m_SignalGen.m_MotionVolumes.size()::New(); duplicator->SetInputImage(m_Parameters.m_SignalGen.m_MaskImage); duplicator->Update(); m_TransformedMaskImage = duplicator->GetOutput(); // second upsampling needed for motion artifacts ImageRegion<3> upsampledImageRegion = m_WorkingImageRegion; DoubleVectorType upsampledSpacing = m_WorkingSpacing; upsampledSpacing[0] /= 4; upsampledSpacing[1] /= 4; upsampledSpacing[2] /= 4; upsampledImageRegion.SetSize(0, m_WorkingImageRegion.GetSize()[0]*4); upsampledImageRegion.SetSize(1, m_WorkingImageRegion.GetSize()[1]*4); upsampledImageRegion.SetSize(2, m_WorkingImageRegion.GetSize()[2]*4); itk::Point upsampledOrigin = m_WorkingOrigin; upsampledOrigin[0] -= m_WorkingSpacing[0]/2; upsampledOrigin[0] += upsampledSpacing[0]/2; upsampledOrigin[1] -= m_WorkingSpacing[1]/2; upsampledOrigin[1] += upsampledSpacing[1]/2; upsampledOrigin[2] -= m_WorkingSpacing[2]/2; upsampledOrigin[2] += upsampledSpacing[2]/2; m_UpsampledMaskImage = ItkUcharImgType::New(); auto upsampler = itk::ResampleImageFilter::New(); upsampler->SetInput(m_Parameters.m_SignalGen.m_MaskImage); upsampler->SetOutputParametersFromImage(m_Parameters.m_SignalGen.m_MaskImage); upsampler->SetSize(upsampledImageRegion.GetSize()); upsampler->SetOutputSpacing(upsampledSpacing); upsampler->SetOutputOrigin(upsampledOrigin); auto nn_interpolator = itk::NearestNeighborInterpolateImageFunction::New(); upsampler->SetInterpolator(nn_interpolator); upsampler->Update(); m_UpsampledMaskImage = upsampler->GetOutput(); } template< class PixelType > void TractsToDWIImageFilter< PixelType >::InitializeFiberData() { // resample fiber bundle for sufficient voxel coverage PrintToLog("Resampling fibers ..."); m_SegmentVolume = 0.0001; float minSpacing = 1; if( m_WorkingSpacing[0]GetDeepCopy(); double volumeAccuracy = 10; m_FiberBundleWorkingCopy->ResampleLinear(minSpacing/volumeAccuracy); m_mmRadius = m_Parameters.m_SignalGen.m_AxonRadius/1000; auto caster = itk::CastImageFilter< itk::Image, itk::Image >::New(); caster->SetInput(m_TransformedMaskImage); caster->Update(); auto density_calculator = itk::TractDensityImageFilter< itk::Image >::New(); density_calculator->SetFiberBundle(m_FiberBundleWorkingCopy); density_calculator->SetInputImage(caster->GetOutput()); density_calculator->SetBinaryOutput(false); density_calculator->SetUseImageGeometry(true); density_calculator->SetDoFiberResampling(false); density_calculator->SetOutputAbsoluteValues(true); density_calculator->SetWorkOnFiberCopy(false); density_calculator->Update(); float max_density = density_calculator->GetMaxDensity(); if (m_mmRadius>0) { m_SegmentVolume = M_PI*m_mmRadius*m_mmRadius*minSpacing/volumeAccuracy; std::stringstream stream; stream << std::fixed << setprecision(2) << max_density * m_SegmentVolume; std::string s = stream.str(); PrintToLog("\nMax. fiber volume: " + s + "mm².", false, true, true); } else { std::stringstream stream; stream << std::fixed << setprecision(2) << max_density * m_SegmentVolume; std::string s = stream.str(); PrintToLog("\nMax. fiber volume: " + s + "mm² (before rescaling to voxel volume).", false, true, true); } float voxel_volume = m_WorkingSpacing[0]*m_WorkingSpacing[1]*m_WorkingSpacing[2]; float new_seg_vol = voxel_volume/max_density; float new_fib_radius = 1000*std::sqrt(new_seg_vol*volumeAccuracy/(minSpacing*M_PI)); std::stringstream stream; stream << std::fixed << setprecision(2) << new_fib_radius; std::string s = stream.str(); PrintToLog("\nA full fiber voxel corresponds to a fiber radius of ~" + s + "µm, given the current fiber configuration.", false, true, true); // a second fiber bundle is needed to store the transformed version of the m_FiberBundleWorkingCopy m_FiberBundleTransformed = m_FiberBundleWorkingCopy; } template< class PixelType > bool TractsToDWIImageFilter< PixelType >::PrepareLogFile() { assert( ! m_Logfile.is_open() ); std::string filePath; std::string fileName; // Get directory name: if (m_Parameters.m_Misc.m_OutputPath.size() > 0) { filePath = m_Parameters.m_Misc.m_OutputPath; if( *(--(filePath.cend())) != '/') { filePath.push_back('/'); } } else { filePath = mitk::IOUtil::GetTempPath() + '/'; } // check if directory exists, else use /tmp/: if( itksys::SystemTools::FileIsDirectory( filePath ) ) { while( *(--(filePath.cend())) == '/') { filePath.pop_back(); } filePath = filePath + '/'; } else { filePath = mitk::IOUtil::GetTempPath() + '/'; } // Get file name: if( ! m_Parameters.m_Misc.m_ResultNode->GetName().empty() ) { fileName = m_Parameters.m_Misc.m_ResultNode->GetName(); } else { fileName = ""; } if( ! m_Parameters.m_Misc.m_OutputPrefix.empty() ) { fileName = m_Parameters.m_Misc.m_OutputPrefix + fileName; } else { fileName = "fiberfox"; } // check if file already exists and DO NOT overwrite existing files: std::string NameTest = fileName; int c = 0; while( itksys::SystemTools::FileExists( filePath + '/' + fileName + ".log" ) && c <= std::numeric_limits::max() ) { fileName = NameTest + "_" + boost::lexical_cast(c); ++c; } try { m_Logfile.open( ( filePath + '/' + fileName + ".log" ).c_str() ); } catch (const std::ios_base::failure &fail) { MITK_ERROR << "itkTractsToDWIImageFilter.cpp: Exception " << fail.what() << " while trying to open file" << filePath << '/' << fileName << ".log"; return false; } if ( m_Logfile.is_open() ) { PrintToLog( "Logfile: " + filePath + '/' + fileName + ".log", false ); return true; } else { m_StatusText += "Logfile could not be opened!\n"; MITK_ERROR << "itkTractsToDWIImageFilter.cpp: Logfile could not be opened!"; return false; } } template< class PixelType > void TractsToDWIImageFilter< PixelType >::GenerateData() { // prepare logfile if ( ! PrepareLogFile() ) { this->SetAbortGenerateData( true ); return; } m_TimeProbe.Start(); // check input data if (m_FiberBundle.IsNull() && m_InputImage.IsNull()) itkExceptionMacro("Input fiber bundle and input diffusion-weighted image is nullptr!"); if (m_Parameters.m_FiberModelList.empty() && m_InputImage.IsNull()) itkExceptionMacro("No diffusion model for fiber compartments defined and input diffusion-weighted" " image is nullptr! At least one fiber compartment is necessary to simulate diffusion."); if (m_Parameters.m_NonFiberModelList.empty() && m_InputImage.IsNull()) itkExceptionMacro("No diffusion model for non-fiber compartments defined and input diffusion-weighted" " image is nullptr! At least one non-fiber compartment is necessary to simulate diffusion."); int baselineIndex = m_Parameters.m_SignalGen.GetFirstBaselineIndex(); if (baselineIndex<0) { itkExceptionMacro("No baseline index found!"); } if (!m_Parameters.m_SignalGen.m_SimulateKspaceAcquisition) // No upsampling of input image needed if no k-space simulation is performed { m_Parameters.m_SignalGen.m_DoAddGibbsRinging = false; } if (m_UseConstantRandSeed) // always generate the same random numbers? { m_RandGen->SetSeed(0); } else { m_RandGen->SetSeed(); } InitializeData(); if ( m_FiberBundle.IsNotNull() ) // if no fiber bundle is found, we directly proceed to the k-space acquisition simulation { CheckVolumeFractionImages(); InitializeFiberData(); int numFiberCompartments = m_Parameters.m_FiberModelList.size(); int numNonFiberCompartments = m_Parameters.m_NonFiberModelList.size(); double maxVolume = 0; unsigned long lastTick = 0; int signalModelSeed = m_RandGen->GetIntegerVariate(); PrintToLog("\n", false, false); PrintToLog("Generating " + boost::lexical_cast(numFiberCompartments+numNonFiberCompartments) + "-compartment diffusion-weighted signal."); std::vector< int > bVals = m_Parameters.m_SignalGen.GetBvalues(); PrintToLog("b-values: ", false, false, true); for (auto v : bVals) PrintToLog(boost::lexical_cast(v) + " ", false, false, true); PrintToLog("\n", false, false, true); PrintToLog("\n", false, false, true); int numFibers = m_FiberBundleWorkingCopy->GetNumFibers(); boost::progress_display disp(numFibers*m_Parameters.m_SignalGen.GetNumVolumes()); PrintToLog("0% 10 20 30 40 50 60 70 80 90 100%", false, true, false); PrintToLog("|----|----|----|----|----|----|----|----|----|----|\n*", false, false, false); for (unsigned int g=0; gSetSeed(signalModelSeed); for (std::size_t i=0; iSetSeed(signalModelSeed); // storing voxel-wise intra-axonal volume in mm³ auto intraAxonalVolumeImage = ItkDoubleImgType::New(); intraAxonalVolumeImage->SetSpacing( m_WorkingSpacing ); intraAxonalVolumeImage->SetOrigin( m_WorkingOrigin ); intraAxonalVolumeImage->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection ); intraAxonalVolumeImage->SetLargestPossibleRegion( m_WorkingImageRegion ); intraAxonalVolumeImage->SetBufferedRegion( m_WorkingImageRegion ); intraAxonalVolumeImage->SetRequestedRegion( m_WorkingImageRegion ); intraAxonalVolumeImage->Allocate(); intraAxonalVolumeImage->FillBuffer(0); maxVolume = 0; + if (this->GetAbortGenerateData()) + continue; + vtkPolyData* fiberPolyData = m_FiberBundleTransformed->GetFiberPolyData(); // generate fiber signal (if there are any fiber models present) if (!m_Parameters.m_FiberModelList.empty()) + { +#pragma omp parallel for for( int i=0; iGetAbortGenerateData()) + continue; + float fiberWeight = m_FiberBundleTransformed->GetFiberWeight(i); - vtkCell* cell = fiberPolyData->GetCell(i); - int numPoints = cell->GetNumberOfPoints(); - vtkPoints* points = cell->GetPoints(); + + int numPoints = -1; + std::vector< itk::Vector > points_copy; +#pragma omp critical + { + vtkCell* cell = fiberPolyData->GetCell(i); + numPoints = cell->GetNumberOfPoints(); + vtkPoints* points = cell->GetPoints(); + for (int j=0; jGetPoint(j))); + } if (numPoints<2) continue; for( int j=0; jGetAbortGenerateData()) { - PrintToLog("\n", false, false); - PrintToLog("Simulation aborted"); - return; + j=numPoints; + continue; } - double* temp = points->GetPoint(j); - itk::Point vertex = GetItkPoint(temp); - itk::Vector v = GetItkVector(temp); + itk::Point vertex = points_copy.at(j); + itk::Vector v = points_copy.at(j); itk::Vector dir(3); - if (jGetPoint(j+1))-v; } - else { dir = v-GetItkVector(points->GetPoint(j-1)); } + if (j idx; itk::ContinuousIndex contIndex; m_TransformedMaskImage->TransformPhysicalPointToIndex(vertex, idx); m_TransformedMaskImage->TransformPhysicalPointToContinuousIndex(vertex, contIndex); if (!m_TransformedMaskImage->GetLargestPossibleRegion().IsInside(idx) || m_TransformedMaskImage->GetPixel(idx)<=0) - { continue; - } + dir.Normalize(); // generate signal for each fiber compartment for (int k=0; kSetFiberDirection(dir); DoubleDwiType::PixelType pix = m_CompartmentImages.at(k)->GetPixel(idx); - pix[g] += fiberWeight*m_SegmentVolume*m_Parameters.m_FiberModelList[k]->SimulateMeasurement(g); + pix[g] += fiberWeight*m_SegmentVolume*m_Parameters.m_FiberModelList[k]->SimulateMeasurement(g, dir); m_CompartmentImages.at(k)->SetPixel(idx, pix); } // update fiber volume image double vol = intraAxonalVolumeImage->GetPixel(idx) + m_SegmentVolume*fiberWeight; intraAxonalVolumeImage->SetPixel(idx, vol); // we assume that the first volume is always unweighted! if (vol>maxVolume) { maxVolume = vol; } } // progress report ++disp; unsigned long newTick = 50*disp.count()/disp.expected_count(); for (unsigned int tick = 0; tick<(newTick-lastTick); tick++) { PrintToLog("*", false, false, false); } lastTick = newTick; } + } // generate non-fiber signal ImageRegionIterator it3(m_TransformedMaskImage, m_TransformedMaskImage->GetLargestPossibleRegion()); double fact = 1; // density correction factor in mm³ if (m_Parameters.m_SignalGen.m_AxonRadius<0.0001 || maxVolume>m_VoxelVolume) // the fullest voxel is always completely full fact = m_VoxelVolume/maxVolume; while(!it3.IsAtEnd()) { if (it3.Get()>0) { DoubleDwiType::IndexType index = it3.GetIndex(); itk::Point point; m_TransformedMaskImage->TransformIndexToPhysicalPoint(index, point); if ( m_Parameters.m_SignalGen.m_DoAddMotion && g>=0 && m_Parameters.m_SignalGen.m_MotionVolumes[g] ) { if (m_Parameters.m_SignalGen.m_DoRandomizeMotion) { point = m_FiberBundleWorkingCopy->TransformPoint( point.GetVnlVector(), -m_Rotation[0], -m_Rotation[1], -m_Rotation[2], -m_Translation[0], -m_Translation[1], -m_Translation[2] ); } else { point = m_FiberBundleWorkingCopy->TransformPoint( point.GetVnlVector(), -m_Rotation[0]*m_MotionCounter, -m_Rotation[1]*m_MotionCounter, -m_Rotation[2]*m_MotionCounter, -m_Translation[0]*m_MotionCounter, -m_Translation[1]*m_MotionCounter, -m_Translation[2]*m_MotionCounter ); } } double iAxVolume = intraAxonalVolumeImage->GetPixel(index); // if volume fraction image is set use it, otherwise use scaling factor to obtain one full fiber voxel double fact2 = fact; if ( m_Parameters.m_FiberModelList[0]->GetVolumeFractionImage()!=nullptr && iAxVolume>0.0001 ) { m_DoubleInterpolator->SetInputImage(m_Parameters.m_FiberModelList[0]->GetVolumeFractionImage()); double val = mitk::imv::GetImageValue(point, true, m_DoubleInterpolator); if (val<0) mitkThrow() << "Volume fraction image (index 1) contains negative values (intra-axonal compartment)!"; fact2 = m_VoxelVolume*val/iAxVolume; } // adjust intra-axonal image value for (int i=0; iGetPixel(index); pix[g] *= fact2; m_CompartmentImages.at(i)->SetPixel(index, pix); } // simulate other compartments SimulateExtraAxonalSignal(index, iAxVolume*fact2, g); } ++it3; } } PrintToLog("\n", false); - if (this->GetAbortGenerateData()) - { - PrintToLog("\n", false, false); - PrintToLog("Simulation aborted"); - return; - } + } + + if (this->GetAbortGenerateData()) + { + PrintToLog("\n", false, false); + PrintToLog("Simulation aborted"); + return; } DoubleDwiType::Pointer doubleOutImage; double signalScale = m_Parameters.m_SignalGen.m_SignalScale; if ( m_Parameters.m_SignalGen.m_SimulateKspaceAcquisition ) // do k-space stuff { PrintToLog("\n", false, false); PrintToLog("Simulating k-space acquisition using " +boost::lexical_cast(m_Parameters.m_SignalGen.m_NumberOfCoils) +" coil(s)"); switch (m_Parameters.m_SignalGen.m_AcquisitionType) { case SignalGenerationParameters::SingleShotEpi: { PrintToLog("Acquisition type: single shot EPI", false); break; } case SignalGenerationParameters::SpinEcho: { PrintToLog("Acquisition type: classic spin echo with cartesian k-space trajectory", false); break; } default: { PrintToLog("Acquisition type: single shot EPI", false); break; } } if (m_Parameters.m_SignalGen.m_NoiseVariance>0) PrintToLog("Simulating complex Gaussian noise", false); if (m_Parameters.m_SignalGen.m_DoSimulateRelaxation) PrintToLog("Simulating signal relaxation", false); if (m_Parameters.m_SignalGen.m_FrequencyMap.IsNotNull()) PrintToLog("Simulating distortions", false); if (m_Parameters.m_SignalGen.m_DoAddGibbsRinging) PrintToLog("Simulating ringing artifacts", false); if (m_Parameters.m_SignalGen.m_EddyStrength>0) PrintToLog("Simulating eddy currents", false); if (m_Parameters.m_SignalGen.m_Spikes>0) PrintToLog("Simulating spikes", false); if (m_Parameters.m_SignalGen.m_CroppingFactor<1.0) PrintToLog("Simulating aliasing artifacts", false); if (m_Parameters.m_SignalGen.m_KspaceLineOffset>0) PrintToLog("Simulating ghosts", false); doubleOutImage = SimulateKspaceAcquisition(m_CompartmentImages); signalScale = 1; // already scaled in SimulateKspaceAcquisition() } else // don't do k-space stuff, just sum compartments { PrintToLog("Summing compartments"); doubleOutImage = m_CompartmentImages.at(0); for (unsigned int i=1; i::New(); adder->SetInput1(doubleOutImage); adder->SetInput2(m_CompartmentImages.at(i)); adder->Update(); doubleOutImage = adder->GetOutput(); } } if (this->GetAbortGenerateData()) { PrintToLog("\n", false, false); PrintToLog("Simulation aborted"); return; } PrintToLog("Finalizing image"); if (signalScale>1) PrintToLog(" Scaling signal", false); if (m_Parameters.m_NoiseModel) PrintToLog(" Adding noise", false); unsigned int window = 0; unsigned int min = itk::NumericTraits::max(); ImageRegionIterator it4 (m_OutputImage, m_OutputImage->GetLargestPossibleRegion()); DoubleDwiType::PixelType signal; signal.SetSize(m_Parameters.m_SignalGen.GetNumVolumes()); boost::progress_display disp2(m_OutputImage->GetLargestPossibleRegion().GetNumberOfPixels()); PrintToLog("0% 10 20 30 40 50 60 70 80 90 100%", false, true, false); PrintToLog("|----|----|----|----|----|----|----|----|----|----|\n*", false, false, false); int lastTick = 0; while(!it4.IsAtEnd()) { if (this->GetAbortGenerateData()) { PrintToLog("\n", false, false); PrintToLog("Simulation aborted"); return; } ++disp2; unsigned long newTick = 50*disp2.count()/disp2.expected_count(); for (unsigned long tick = 0; tick<(newTick-lastTick); tick++) PrintToLog("*", false, false, false); lastTick = newTick; typename OutputImageType::IndexType index = it4.GetIndex(); signal = doubleOutImage->GetPixel(index)*signalScale; if (m_Parameters.m_NoiseModel) m_Parameters.m_NoiseModel->AddNoise(signal); for (unsigned int i=0; i0) signal[i] = floor(signal[i]+0.5); else signal[i] = ceil(signal[i]-0.5); if ( (!m_Parameters.m_SignalGen.IsBaselineIndex(i) || signal.Size()==1) && signal[i]>window) window = signal[i]; if ( (!m_Parameters.m_SignalGen.IsBaselineIndex(i) || signal.Size()==1) && signal[i]SetNthOutput(0, m_OutputImage); PrintToLog("\n", false); PrintToLog("Finished simulation"); m_TimeProbe.Stop(); if (m_Parameters.m_SignalGen.m_DoAddMotion) { PrintToLog("\nHead motion log:", false); PrintToLog(m_MotionLog, false, false); } if (m_Parameters.m_SignalGen.m_Spikes>0) { PrintToLog("\nSpike log:", false); PrintToLog(m_SpikeLog, false, false); } if (m_Logfile.is_open()) m_Logfile.close(); } template< class PixelType > void TractsToDWIImageFilter< PixelType >::PrintToLog(std::string m, bool addTime, bool linebreak, bool stdOut) { // timestamp if (addTime) { m_Logfile << this->GetTime() << " > "; m_StatusText += this->GetTime() + " > "; if (stdOut) std::cout << this->GetTime() << " > "; } // message if (m_Logfile.is_open()) m_Logfile << m; m_StatusText += m; if (stdOut) std::cout << m; // new line if (linebreak) { if (m_Logfile.is_open()) m_Logfile << "\n"; m_StatusText += "\n"; if (stdOut) std::cout << "\n"; } } template< class PixelType > void TractsToDWIImageFilter< PixelType >::SimulateMotion(int g) { // is motion artifact enabled? // is the current volume g affected by motion? if ( m_Parameters.m_SignalGen.m_DoAddMotion && m_Parameters.m_SignalGen.m_MotionVolumes[g] && g(m_Parameters.m_SignalGen.GetNumVolumes()) ) { if ( m_Parameters.m_SignalGen.m_DoRandomizeMotion ) { // either undo last transform or work on fresh copy of untransformed fibers m_FiberBundleTransformed = m_FiberBundleWorkingCopy->GetDeepCopy(); m_Rotation[0] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_SignalGen.m_Rotation[0]*2) -m_Parameters.m_SignalGen.m_Rotation[0]; m_Rotation[1] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_SignalGen.m_Rotation[1]*2) -m_Parameters.m_SignalGen.m_Rotation[1]; m_Rotation[2] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_SignalGen.m_Rotation[2]*2) -m_Parameters.m_SignalGen.m_Rotation[2]; m_Translation[0] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_SignalGen.m_Translation[0]*2) -m_Parameters.m_SignalGen.m_Translation[0]; m_Translation[1] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_SignalGen.m_Translation[1]*2) -m_Parameters.m_SignalGen.m_Translation[1]; m_Translation[2] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_SignalGen.m_Translation[2]*2) -m_Parameters.m_SignalGen.m_Translation[2]; } else { m_Rotation = m_Parameters.m_SignalGen.m_Rotation / m_NumMotionVolumes; m_Translation = m_Parameters.m_SignalGen.m_Translation / m_NumMotionVolumes; m_MotionCounter++; } // move mask image if (m_MaskImageSet) { ImageRegionIterator maskIt(m_UpsampledMaskImage, m_UpsampledMaskImage->GetLargestPossibleRegion()); m_TransformedMaskImage->FillBuffer(0); while(!maskIt.IsAtEnd()) { if (maskIt.Get()<=0) { ++maskIt; continue; } DoubleDwiType::IndexType index = maskIt.GetIndex(); itk::Point point; m_UpsampledMaskImage->TransformIndexToPhysicalPoint(index, point); if (m_Parameters.m_SignalGen.m_DoRandomizeMotion) { point = m_FiberBundleWorkingCopy->TransformPoint(point.GetVnlVector(), m_Rotation[0],m_Rotation[1],m_Rotation[2], m_Translation[0],m_Translation[1],m_Translation[2]); } else { point = m_FiberBundleWorkingCopy->TransformPoint(point.GetVnlVector(), m_Rotation[0]*m_MotionCounter,m_Rotation[1]*m_MotionCounter,m_Rotation[2]*m_MotionCounter, m_Translation[0]*m_MotionCounter,m_Translation[1]*m_MotionCounter,m_Translation[2]*m_MotionCounter); } m_TransformedMaskImage->TransformPhysicalPointToIndex(point, index); if (m_TransformedMaskImage->GetLargestPossibleRegion().IsInside(index)) { m_TransformedMaskImage->SetPixel(index,100); } ++maskIt; } } if (m_Parameters.m_SignalGen.m_DoRandomizeMotion) { m_Rotations.push_back(m_Rotation); m_Translations.push_back(m_Translation); m_MotionLog += boost::lexical_cast(g) + " rotation: " + boost::lexical_cast(m_Rotation[0]) + "," + boost::lexical_cast(m_Rotation[1]) + "," + boost::lexical_cast(m_Rotation[2]) + ";"; m_MotionLog += " translation: " + boost::lexical_cast(m_Translation[0]) + "," + boost::lexical_cast(m_Translation[1]) + "," + boost::lexical_cast(m_Translation[2]) + "\n"; } else { m_Rotations.push_back(m_Rotation*m_MotionCounter); m_Translations.push_back(m_Translation*m_MotionCounter); m_MotionLog += boost::lexical_cast(g) + " rotation: " + boost::lexical_cast(m_Rotation[0]*m_MotionCounter) + "," + boost::lexical_cast(m_Rotation[1]*m_MotionCounter) + "," + boost::lexical_cast(m_Rotation[2]*m_MotionCounter) + ";"; m_MotionLog += " translation: " + boost::lexical_cast(m_Translation[0]*m_MotionCounter) + "," + boost::lexical_cast(m_Translation[1]*m_MotionCounter) + "," + boost::lexical_cast(m_Translation[2]*m_MotionCounter) + "\n"; } m_FiberBundleTransformed->TransformFibers(m_Rotation[0],m_Rotation[1],m_Rotation[2],m_Translation[0],m_Translation[1],m_Translation[2]); } else { m_Rotation.Fill(0.0); m_Translation.Fill(0.0); m_Rotations.push_back(m_Rotation); m_Translations.push_back(m_Translation); m_MotionLog += boost::lexical_cast(g) + " rotation: " + boost::lexical_cast(m_Rotation[0]) + "," + boost::lexical_cast(m_Rotation[1]) + "," + boost::lexical_cast(m_Rotation[2]) + ";"; m_MotionLog += " translation: " + boost::lexical_cast(m_Translation[0]) + "," + boost::lexical_cast(m_Translation[1]) + "," + boost::lexical_cast(m_Translation[2]) + "\n"; } } template< class PixelType > void TractsToDWIImageFilter< PixelType >:: SimulateExtraAxonalSignal(ItkUcharImgType::IndexType index, double intraAxonalVolume, int g) { int numFiberCompartments = m_Parameters.m_FiberModelList.size(); int numNonFiberCompartments = m_Parameters.m_NonFiberModelList.size(); if (intraAxonalVolume>0.0001 && m_Parameters.m_SignalGen.m_DoDisablePartialVolume) // only fiber in voxel { DoubleDwiType::PixelType pix = m_CompartmentImages.at(0)->GetPixel(index); if (g>=0) pix[g] *= m_VoxelVolume/intraAxonalVolume; else pix *= m_VoxelVolume/intraAxonalVolume; m_CompartmentImages.at(0)->SetPixel(index, pix); if (g==0) m_VolumeFractions.at(0)->SetPixel(index, 1); for (int i=1; iGetPixel(index); if (g>=0) pix[g] = 0.0; else pix.Fill(0.0); m_CompartmentImages.at(i)->SetPixel(index, pix); } } else { if (g==0) { m_VolumeFractions.at(0)->SetPixel(index, intraAxonalVolume/m_VoxelVolume); } // get non-transformed point (remove headmotion tranformation) // this point can then be transformed to each of the original images, regardless of their geometry itk::Point point; m_TransformedMaskImage->TransformIndexToPhysicalPoint(index, point); if ( m_Parameters.m_SignalGen.m_DoAddMotion && g>=0 && m_Parameters.m_SignalGen.m_MotionVolumes[g] ) { if (m_Parameters.m_SignalGen.m_DoRandomizeMotion) { point = m_FiberBundleWorkingCopy->TransformPoint(point.GetVnlVector(), -m_Rotation[0],-m_Rotation[1],-m_Rotation[2], -m_Translation[0],-m_Translation[1],-m_Translation[2]); } else { point = m_FiberBundleWorkingCopy->TransformPoint(point.GetVnlVector(), -m_Rotation[0]*m_MotionCounter,-m_Rotation[1]*m_MotionCounter,-m_Rotation[2]*m_MotionCounter, -m_Translation[0]*m_MotionCounter,-m_Translation[1]*m_MotionCounter,-m_Translation[2]*m_MotionCounter); } } if (m_Parameters.m_SignalGen.m_DoDisablePartialVolume) { int maxVolumeIndex = 0; double maxWeight = 0; for (int i=0; i1) { m_DoubleInterpolator->SetInputImage(m_Parameters.m_NonFiberModelList[i]->GetVolumeFractionImage()); double val = mitk::imv::GetImageValue(point, true, m_DoubleInterpolator); if (val<0) mitkThrow() << "Volume fraction image (index " << i << ") contains values less than zero!"; else weight = val; } if (weight>maxWeight) { maxWeight = weight; maxVolumeIndex = i; } } DoubleDwiType::Pointer doubleDwi = m_CompartmentImages.at(maxVolumeIndex+numFiberCompartments); DoubleDwiType::PixelType pix = doubleDwi->GetPixel(index); if (g>=0) - pix[g] += m_Parameters.m_NonFiberModelList[maxVolumeIndex]->SimulateMeasurement(g)*m_VoxelVolume; + pix[g] += m_Parameters.m_NonFiberModelList[maxVolumeIndex]->SimulateMeasurement(g, m_NullDir)*m_VoxelVolume; else - pix += m_Parameters.m_NonFiberModelList[maxVolumeIndex]->SimulateMeasurement()*m_VoxelVolume; + pix += m_Parameters.m_NonFiberModelList[maxVolumeIndex]->SimulateMeasurement(m_NullDir)*m_VoxelVolume; doubleDwi->SetPixel(index, pix); if (g==0) m_VolumeFractions.at(maxVolumeIndex+numFiberCompartments)->SetPixel(index, 1); } else { double extraAxonalVolume = m_VoxelVolume-intraAxonalVolume; // non-fiber volume if (extraAxonalVolume<0) { if (extraAxonalVolume<-0.001) MITK_ERROR << "Corrupted intra-axonal signal voxel detected. Fiber volume larger voxel volume! " << m_VoxelVolume << "<" << intraAxonalVolume; extraAxonalVolume = 0; } double interAxonalVolume = 0; if (numFiberCompartments>1) interAxonalVolume = extraAxonalVolume * intraAxonalVolume/m_VoxelVolume; // inter-axonal fraction of non fiber compartment double other = extraAxonalVolume - interAxonalVolume; // rest of compartment if (other<0) { if (other<-0.001) MITK_ERROR << "Corrupted signal voxel detected. Fiber volume larger voxel volume!"; other = 0; interAxonalVolume = extraAxonalVolume; } double compartmentSum = intraAxonalVolume; // adjust non-fiber and intra-axonal signal for (int i=1; iGetPixel(index); if (intraAxonalVolume>0) // remove scaling by intra-axonal volume from inter-axonal compartment { if (g>=0) pix[g] /= intraAxonalVolume; else pix /= intraAxonalVolume; } else { if (g>=0) pix[g] = 0; else pix *= 0; } if (m_Parameters.m_FiberModelList[i]->GetVolumeFractionImage()!=nullptr) { m_DoubleInterpolator->SetInputImage(m_Parameters.m_FiberModelList[i]->GetVolumeFractionImage()); double val = mitk::imv::GetImageValue(point, true, m_DoubleInterpolator); if (val<0) mitkThrow() << "Volume fraction image (index " << i+1 << ") contains negative values!"; else weight = val*m_VoxelVolume; } compartmentSum += weight; if (g>=0) pix[g] *= weight; else pix *= weight; m_CompartmentImages.at(i)->SetPixel(index, pix); if (g==0) m_VolumeFractions.at(i)->SetPixel(index, weight/m_VoxelVolume); } for (int i=0; iGetPixel(index); if (m_Parameters.m_NonFiberModelList[i]->GetVolumeFractionImage()!=nullptr) { m_DoubleInterpolator->SetInputImage(m_Parameters.m_NonFiberModelList[i]->GetVolumeFractionImage()); double val = mitk::imv::GetImageValue(point, true, m_DoubleInterpolator); if (val<0) mitkThrow() << "Volume fraction image (index " << numFiberCompartments+i+1 << ") contains negative values (non-fiber compartment)!"; else weight = val*m_VoxelVolume; if (m_UseRelativeNonFiberVolumeFractions) weight *= other/m_VoxelVolume; } compartmentSum += weight; if (g>=0) - pix[g] += m_Parameters.m_NonFiberModelList[i]->SimulateMeasurement(g)*weight; + pix[g] += m_Parameters.m_NonFiberModelList[i]->SimulateMeasurement(g, m_NullDir)*weight; else - pix += m_Parameters.m_NonFiberModelList[i]->SimulateMeasurement()*weight; + pix += m_Parameters.m_NonFiberModelList[i]->SimulateMeasurement(m_NullDir)*weight; m_CompartmentImages.at(i+numFiberCompartments)->SetPixel(index, pix); if (g==0) m_VolumeFractions.at(i+numFiberCompartments)->SetPixel(index, weight/m_VoxelVolume); } if (compartmentSum/m_VoxelVolume>1.05) MITK_ERROR << "Compartments do not sum to 1 in voxel " << index << " (" << compartmentSum/m_VoxelVolume << ")"; } } } 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; } template< class PixelType > double TractsToDWIImageFilter< PixelType >::RoundToNearest(double num) { return (num > 0.0) ? floor(num + 0.5) : ceil(num - 0.5); } template< class PixelType > std::string TractsToDWIImageFilter< PixelType >::GetTime() { m_TimeProbe.Stop(); unsigned long total = RoundToNearest(m_TimeProbe.GetTotal()); unsigned long hours = total/3600; unsigned long minutes = (total%3600)/60; unsigned long seconds = total%60; std::string out = ""; out.append(boost::lexical_cast(hours)); out.append(":"); out.append(boost::lexical_cast(minutes)); out.append(":"); out.append(boost::lexical_cast(seconds)); m_TimeProbe.Start(); return out; } } diff --git a/Modules/DiffusionImaging/FiberTracking/Fiberfox/itkTractsToDWIImageFilter.h b/Modules/DiffusionImaging/FiberTracking/Fiberfox/itkTractsToDWIImageFilter.h index 3e631ebf57..89d8d3dc26 100755 --- a/Modules/DiffusionImaging/FiberTracking/Fiberfox/itkTractsToDWIImageFilter.h +++ b/Modules/DiffusionImaging/FiberTracking/Fiberfox/itkTractsToDWIImageFilter.h @@ -1,164 +1,165 @@ /*=================================================================== 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__ #include #include #include #include #include #include #include #include #include #include namespace itk { /** * \brief Generates artificial diffusion weighted image volume from the input fiberbundle using a generic multicompartment model. * See "Fiberfox: Facilitating the creation of realistic white matter software phantoms" (DOI: 10.1002/mrm.25045) for details. */ 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 ItkDoubleImgType4D; typedef itk::Image ItkDoubleImgType; typedef itk::Image ItkFloatImgType; typedef itk::Image ItkUcharImgType; typedef mitk::FiberBundle::Pointer FiberBundleType; typedef itk::VectorImage< double, 3 > DoubleDwiType; typedef itk::Matrix MatrixType; - typedef itk::Image< double, 2 > SliceType; - typedef itk::VnlForwardFFTImageFilter::OutputImageType ComplexSliceType; + typedef itk::Image< double, 2 > Double2DImageType; + typedef itk::VnlForwardFFTImageFilter::OutputImageType Complex2DImageType; typedef itk::VectorImage< vcl_complex< double >, 3 > ComplexDwiType; typedef itk::Vector< double,3> DoubleVectorType; itkFactorylessNewMacro(Self) itkCloneMacro(Self) itkTypeMacro( TractsToDWIImageFilter, ImageSource ) /** Input */ itkSetMacro( FiberBundle, FiberBundleType ) ///< Input fiber bundle itkSetMacro( InputImage, typename OutputImageType::Pointer ) ///< Input diffusion-weighted image. If no fiber bundle is set, then the acquisition is simulated for this image without a new diffusion simulation. itkSetMacro( UseConstantRandSeed, bool ) ///< Seed for random generator. void SetParameters( FiberfoxParameters param ) ///< Simulation parameters. { m_Parameters = param; } /** Output */ FiberfoxParameters GetParameters(){ return m_Parameters; } std::vector< ItkDoubleImgType::Pointer > GetVolumeFractions() ///< one double image for each compartment containing the corresponding volume fraction per voxel { return m_VolumeFractions; } mitk::LevelWindow GetLevelWindow() ///< Level window is determined from the output image { return m_LevelWindow; } itkGetMacro( StatusText, std::string ) itkGetMacro( PhaseImage, DoubleDwiType::Pointer ) itkGetMacro( KspaceImage, DoubleDwiType::Pointer ) itkGetMacro( CoilPointset, mitk::PointSet::Pointer ) 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); double RoundToNearest(double num); std::string GetTime(); bool PrepareLogFile(); /** Prepares the log file and returns true if successful or false if failed. */ void PrintToLog(std::string m, bool addTime=true, bool linebreak=true, bool stdOut=true); /** Transform generated image compartment by compartment, channel by channel and slice by slice using DFT and add k-space artifacts/effects. */ DoubleDwiType::Pointer SimulateKspaceAcquisition(std::vector< DoubleDwiType::Pointer >& images); /** Generate signal of non-fiber compartments. */ void SimulateExtraAxonalSignal(ItkUcharImgType::IndexType index, double intraAxonalVolume, int g=-1); /** Move fibers to simulate headmotion */ void SimulateMotion(int g=-1); void CheckVolumeFractionImages(); ItkDoubleImgType::Pointer NormalizeInsideMask(ItkDoubleImgType::Pointer image); void InitializeData(); void InitializeFiberData(); // input mitk::FiberfoxParameters m_Parameters; FiberBundleType m_FiberBundle; typename OutputImageType::Pointer m_InputImage; // output typename OutputImageType::Pointer m_OutputImage; typename DoubleDwiType::Pointer m_PhaseImage; typename DoubleDwiType::Pointer m_KspaceImage; mitk::LevelWindow m_LevelWindow; std::vector< ItkDoubleImgType::Pointer > m_VolumeFractions; std::string m_StatusText; // MISC itk::TimeProbe m_TimeProbe; bool m_UseConstantRandSeed; bool m_MaskImageSet; ofstream m_Logfile; std::string m_MotionLog; std::string m_SpikeLog; // signal generation FiberBundleType m_FiberBundleWorkingCopy; ///< we work on an upsampled version of the input bundle FiberBundleType m_FiberBundleTransformed; ///< transformed bundle simulating headmotion itk::Vector m_WorkingSpacing; itk::Point m_WorkingOrigin; ImageRegion<3> m_WorkingImageRegion; double m_VoxelVolume; std::vector< DoubleDwiType::Pointer > m_CompartmentImages; ItkUcharImgType::Pointer m_TransformedMaskImage; ///< copy of mask image (changes for each motion step) ItkUcharImgType::Pointer m_UpsampledMaskImage; ///< helper image for motion simulation DoubleVectorType m_Rotation; DoubleVectorType m_Translation; std::vector< DoubleVectorType > m_Rotations; /// m_Translations; ///::Pointer m_DoubleInterpolator; + itk::Vector m_NullDir; }; } #ifndef ITK_MANUAL_INSTANTIATION #include "itkTractsToDWIImageFilter.cpp" #endif #endif diff --git a/Modules/DiffusionImaging/FiberTracking/Testing/mitkFiberfoxSignalGenerationTest.cpp b/Modules/DiffusionImaging/FiberTracking/Testing/mitkFiberfoxSignalGenerationTest.cpp index 2a6365c5ab..1c94105b9a 100644 --- a/Modules/DiffusionImaging/FiberTracking/Testing/mitkFiberfoxSignalGenerationTest.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Testing/mitkFiberfoxSignalGenerationTest.cpp @@ -1,259 +1,266 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include - +#include #include "mitkTestFixture.h" class mitkFiberfoxSignalGenerationTestSuite : public mitk::TestFixture { - CPPUNIT_TEST_SUITE(mitkFiberfoxSignalGenerationTestSuite); - //MITK_TEST(Test0); - //MITK_TEST(Test1); - MITK_TEST(Test2); - MITK_TEST(Test3); - MITK_TEST(Test4); - MITK_TEST(Test5); - //MITK_TEST(Test6); - MITK_TEST(Test7); - MITK_TEST(Test8); - CPPUNIT_TEST_SUITE_END(); + CPPUNIT_TEST_SUITE(mitkFiberfoxSignalGenerationTestSuite); + MITK_TEST(Test0); + MITK_TEST(Test1); + MITK_TEST(Test2); + MITK_TEST(Test3); + MITK_TEST(Test4); + MITK_TEST(Test5); + MITK_TEST(Test6); + MITK_TEST(Test7); + MITK_TEST(Test8); + CPPUNIT_TEST_SUITE_END(); - typedef itk::VectorImage< short, 3> ItkDwiType; + typedef itk::VectorImage< short, 3> ItkDwiType; private: public: - /** Members used inside the different (sub-)tests. All members are initialized via setUp().*/ - FiberBundle::Pointer m_FiberBundle; - std::vector< FiberfoxParameters > m_Parameters; - std::vector< mitk::Image::Pointer > m_RefImages; - - void setUp() override - { - m_FiberBundle = dynamic_cast(mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/Fiberfox/Signalgen.fib"))[0].GetPointer()); - - { - FiberfoxParameters parameters; - parameters.LoadParameters(GetTestDataFilePath("DiffusionImaging/Fiberfox/params/param1.ffp")); - m_Parameters.push_back(parameters); - m_RefImages.push_back(dynamic_cast(mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/Fiberfox/params/param1.dwi"))[0].GetPointer())); - } - - { - FiberfoxParameters parameters; - parameters.LoadParameters(GetTestDataFilePath("DiffusionImaging/Fiberfox/params/param2.ffp")); - m_Parameters.push_back(parameters); - m_RefImages.push_back(dynamic_cast(mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/Fiberfox/params/param2.dwi"))[0].GetPointer())); - } + /** Members used inside the different (sub-)tests. All members are initialized via setUp().*/ + FiberBundle::Pointer m_FiberBundle; + std::vector< FiberfoxParameters > m_Parameters; + std::vector< mitk::Image::Pointer > m_RefImages; - { - FiberfoxParameters parameters; - parameters.LoadParameters(GetTestDataFilePath("DiffusionImaging/Fiberfox/params/param3.ffp")); - m_Parameters.push_back(parameters); - m_RefImages.push_back(dynamic_cast(mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/Fiberfox/params/param3.dwi"))[0].GetPointer())); - } + void setUp() override + { + omp_set_num_threads(1); - { - FiberfoxParameters parameters; - parameters.LoadParameters(GetTestDataFilePath("DiffusionImaging/Fiberfox/params/param4.ffp")); - m_Parameters.push_back(parameters); - m_RefImages.push_back(dynamic_cast(mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/Fiberfox/params/param4.dwi"))[0].GetPointer())); - } + m_FiberBundle = dynamic_cast(mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/Fiberfox/Signalgen.fib"))[0].GetPointer()); - { - FiberfoxParameters parameters; - parameters.LoadParameters(GetTestDataFilePath("DiffusionImaging/Fiberfox/params/param5.ffp")); - m_Parameters.push_back(parameters); - m_RefImages.push_back(dynamic_cast(mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/Fiberfox/params/param5.dwi"))[0].GetPointer())); - } - - { - FiberfoxParameters parameters; - parameters.LoadParameters(GetTestDataFilePath("DiffusionImaging/Fiberfox/params/param6.ffp")); - m_Parameters.push_back(parameters); - m_RefImages.push_back(dynamic_cast(mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/Fiberfox/params/param6.dwi"))[0].GetPointer())); - } - - { - FiberfoxParameters parameters; - parameters.LoadParameters(GetTestDataFilePath("DiffusionImaging/Fiberfox/params/param7.ffp")); - m_Parameters.push_back(parameters); - m_RefImages.push_back(dynamic_cast(mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/Fiberfox/params/param7.dwi"))[0].GetPointer())); - } - - { - FiberfoxParameters parameters; - parameters.LoadParameters(GetTestDataFilePath("DiffusionImaging/Fiberfox/params/param8.ffp")); - m_Parameters.push_back(parameters); - m_RefImages.push_back(dynamic_cast(mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/Fiberfox/params/param8.dwi"))[0].GetPointer())); - } - - { - FiberfoxParameters parameters; - parameters.LoadParameters(GetTestDataFilePath("DiffusionImaging/Fiberfox/params/param9.ffp")); - m_Parameters.push_back(parameters); - m_RefImages.push_back(dynamic_cast(mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/Fiberfox/params/param9.dwi"))[0].GetPointer())); - } - } - - void tearDown() override { - + FiberfoxParameters parameters; + parameters.LoadParameters(GetTestDataFilePath("DiffusionImaging/Fiberfox/params/param1.ffp")); + m_Parameters.push_back(parameters); + m_RefImages.push_back(dynamic_cast(mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/Fiberfox/params/param1.dwi"))[0].GetPointer())); } - bool CompareDwi(itk::VectorImage< short, 3 >* dwi1, itk::VectorImage< short, 3 >* dwi2) { - bool out = true; - typedef itk::VectorImage< short, 3 > DwiImageType; - try{ - itk::ImageRegionIterator< DwiImageType > it1(dwi1, dwi1->GetLargestPossibleRegion()); - itk::ImageRegionIterator< DwiImageType > it2(dwi2, dwi2->GetLargestPossibleRegion()); - int count = 0; - while(!it1.IsAtEnd()) - { - if (it1.Get()!=it2.Get()) - { - if (count<10) - { - MITK_INFO << "**************************************"; - MITK_INFO << "Test value: " << it1.GetIndex() << ":" << it1.Get(); - MITK_INFO << "Ref. value: " << it2.GetIndex() << ":" << it2.Get(); - } - out = false; - count++; - } - ++it1; - ++it2; - } - if (count>=10) - MITK_INFO << "Skipping errors."; - MITK_INFO << "Errors detected: " << count; - } - catch(...) - { - return false; - } - return out; + FiberfoxParameters parameters; + parameters.LoadParameters(GetTestDataFilePath("DiffusionImaging/Fiberfox/params/param2.ffp")); + m_Parameters.push_back(parameters); + m_RefImages.push_back(dynamic_cast(mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/Fiberfox/params/param2.dwi"))[0].GetPointer())); } - void StartSimulation(FiberfoxParameters parameters, mitk::Image::Pointer refImage, std::string out) { - itk::TractsToDWIImageFilter< short >::Pointer tractsToDwiFilter = itk::TractsToDWIImageFilter< short >::New(); - tractsToDwiFilter->SetUseConstantRandSeed(true); - tractsToDwiFilter->SetParameters(parameters); - tractsToDwiFilter->SetFiberBundle(m_FiberBundle); - tractsToDwiFilter->Update(); - - mitk::Image::Pointer testImage = mitk::GrabItkImageMemory( tractsToDwiFilter->GetOutput() ); - testImage->GetPropertyList()->ReplaceProperty( mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str(), mitk::GradientDirectionsProperty::New( parameters.m_SignalGen.GetGradientDirections() ) ); - testImage->SetProperty( mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str(), mitk::FloatProperty::New( parameters.m_SignalGen.m_Bvalue ) ); - - mitk::DiffusionPropertyHelper propertyHelper( testImage ); - propertyHelper.InitializeImage(); - - if (refImage.IsNotNull()) - { - if( static_cast( refImage->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() )->GetGradientDirectionsContainer().IsNotNull() ) - { - ItkDwiType::Pointer itkTestImagePointer = ItkDwiType::New(); - mitk::CastToItkImage(testImage, itkTestImagePointer); - ItkDwiType::Pointer itkRefImagePointer = ItkDwiType::New(); - mitk::CastToItkImage(refImage, itkRefImagePointer); - - bool cond = CompareDwi(itkTestImagePointer, itkRefImagePointer); - if (!cond) - { - MITK_INFO << "Saving test image to " << mitk::IOUtil::GetTempPath(); - mitk::IOUtil::Save(testImage, mitk::IOUtil::GetTempPath()+out); - } - CPPUNIT_ASSERT_MESSAGE("Simulated images should be equal", cond); - } - } + FiberfoxParameters parameters; + parameters.LoadParameters(GetTestDataFilePath("DiffusionImaging/Fiberfox/params/param3.ffp")); + m_Parameters.push_back(parameters); + m_RefImages.push_back(dynamic_cast(mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/Fiberfox/params/param3.dwi"))[0].GetPointer())); } - void Test0() { - StartSimulation(m_Parameters.at(0), m_RefImages.at(0), "param1.dwi"); + FiberfoxParameters parameters; + parameters.LoadParameters(GetTestDataFilePath("DiffusionImaging/Fiberfox/params/param4.ffp")); + m_Parameters.push_back(parameters); + m_RefImages.push_back(dynamic_cast(mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/Fiberfox/params/param4.dwi"))[0].GetPointer())); } - void Test1() { - StartSimulation(m_Parameters.at(1), m_RefImages.at(1), "param2.dwi"); + FiberfoxParameters parameters; + parameters.LoadParameters(GetTestDataFilePath("DiffusionImaging/Fiberfox/params/param5.ffp")); + m_Parameters.push_back(parameters); + m_RefImages.push_back(dynamic_cast(mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/Fiberfox/params/param5.dwi"))[0].GetPointer())); } - void Test2() { - StartSimulation(m_Parameters.at(2), m_RefImages.at(2), "param3.dwi"); + FiberfoxParameters parameters; + parameters.LoadParameters(GetTestDataFilePath("DiffusionImaging/Fiberfox/params/param6.ffp")); + m_Parameters.push_back(parameters); + m_RefImages.push_back(dynamic_cast(mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/Fiberfox/params/param6.dwi"))[0].GetPointer())); } - void Test3() { - StartSimulation(m_Parameters.at(3), m_RefImages.at(3), "param4.dwi"); + FiberfoxParameters parameters; + parameters.LoadParameters(GetTestDataFilePath("DiffusionImaging/Fiberfox/params/param7.ffp")); + m_Parameters.push_back(parameters); + m_RefImages.push_back(dynamic_cast(mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/Fiberfox/params/param7.dwi"))[0].GetPointer())); } - void Test4() { - StartSimulation(m_Parameters.at(4), m_RefImages.at(4), "param5.dwi"); + FiberfoxParameters parameters; + parameters.LoadParameters(GetTestDataFilePath("DiffusionImaging/Fiberfox/params/param8.ffp")); + m_Parameters.push_back(parameters); + m_RefImages.push_back(dynamic_cast(mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/Fiberfox/params/param8.dwi"))[0].GetPointer())); } - void Test5() { - StartSimulation(m_Parameters.at(5), m_RefImages.at(5), "param6.dwi"); + FiberfoxParameters parameters; + parameters.LoadParameters(GetTestDataFilePath("DiffusionImaging/Fiberfox/params/param9.ffp")); + m_Parameters.push_back(parameters); + m_RefImages.push_back(dynamic_cast(mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/Fiberfox/params/param9.dwi"))[0].GetPointer())); } + } + + void tearDown() override + { + + } + + bool CompareDwi(itk::VectorImage< short, 3 >* dwi1, itk::VectorImage< short, 3 >* dwi2) + { + bool out = true; + typedef itk::VectorImage< short, 3 > DwiImageType; + try{ + itk::ImageRegionIterator< DwiImageType > it1(dwi1, dwi1->GetLargestPossibleRegion()); + itk::ImageRegionIterator< DwiImageType > it2(dwi2, dwi2->GetLargestPossibleRegion()); + int count = 0; + while(!it1.IsAtEnd()) + { + for (unsigned int i=0; iGetVectorLength(); ++i) + { + short d = abs(it1.Get()[i]-it2.Get()[i]); - void Test6() - { - StartSimulation(m_Parameters.at(6), m_RefImages.at(6), "param7.dwi"); + if (d>0) + { + if (count<10) + { + MITK_INFO << "**************************************"; + MITK_INFO << "Test value: " << it1.GetIndex() << ":" << it1.Get()[i]; + MITK_INFO << "Ref. value: " << it2.GetIndex() << ":" << it2.Get()[i]; + } + out = false; + count++; + } + } + ++it1; + ++it2; + } + if (count>=10) + MITK_INFO << "Skipping errors."; + MITK_INFO << "Errors detected: " << count; } - - void Test7() + catch(...) { - StartSimulation(m_Parameters.at(7), m_RefImages.at(7), "param8.dwi"); + return false; } + return out; + } + + void StartSimulation(FiberfoxParameters parameters, mitk::Image::Pointer refImage, std::string out) + { + itk::TractsToDWIImageFilter< short >::Pointer tractsToDwiFilter = itk::TractsToDWIImageFilter< short >::New(); + tractsToDwiFilter->SetUseConstantRandSeed(true); + tractsToDwiFilter->SetParameters(parameters); + tractsToDwiFilter->SetFiberBundle(m_FiberBundle); + tractsToDwiFilter->Update(); + + mitk::Image::Pointer testImage = mitk::GrabItkImageMemory( tractsToDwiFilter->GetOutput() ); + testImage->GetPropertyList()->ReplaceProperty( mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str(), mitk::GradientDirectionsProperty::New( parameters.m_SignalGen.GetGradientDirections() ) ); + testImage->SetProperty( mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str(), mitk::FloatProperty::New( parameters.m_SignalGen.m_Bvalue ) ); - void Test8() + mitk::DiffusionPropertyHelper propertyHelper( testImage ); + propertyHelper.InitializeImage(); + + if (refImage.IsNotNull()) { - StartSimulation(m_Parameters.at(8), m_RefImages.at(8), "param9.dwi"); + if( static_cast( refImage->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() )->GetGradientDirectionsContainer().IsNotNull() ) + { + ItkDwiType::Pointer itkTestImagePointer = ItkDwiType::New(); + mitk::CastToItkImage(testImage, itkTestImagePointer); + ItkDwiType::Pointer itkRefImagePointer = ItkDwiType::New(); + mitk::CastToItkImage(refImage, itkRefImagePointer); + + bool cond = CompareDwi(itkTestImagePointer, itkRefImagePointer); + if (!cond) + { + MITK_INFO << "Saving test image to " << mitk::IOUtil::GetTempPath(); + mitk::IOUtil::Save(testImage, mitk::IOUtil::GetTempPath()+out); + } + CPPUNIT_ASSERT_MESSAGE("Simulated images should be equal", cond); + } } + } + + void Test0() + { + StartSimulation(m_Parameters.at(0), m_RefImages.at(0), "param1.dwi"); + } + + void Test1() + { + StartSimulation(m_Parameters.at(1), m_RefImages.at(1), "param2.dwi"); + } + + void Test2() + { + StartSimulation(m_Parameters.at(2), m_RefImages.at(2), "param3.dwi"); + } + + void Test3() + { + StartSimulation(m_Parameters.at(3), m_RefImages.at(3), "param4.dwi"); + } + + void Test4() + { + StartSimulation(m_Parameters.at(4), m_RefImages.at(4), "param5.dwi"); + } + + void Test5() + { + StartSimulation(m_Parameters.at(5), m_RefImages.at(5), "param6.dwi"); + } + + void Test6() + { + StartSimulation(m_Parameters.at(6), m_RefImages.at(6), "param7.dwi"); + } + + void Test7() + { + StartSimulation(m_Parameters.at(7), m_RefImages.at(7), "param8.dwi"); + } + + void Test8() + { + StartSimulation(m_Parameters.at(8), m_RefImages.at(8), "param9.dwi"); + } }; MITK_TEST_SUITE_REGISTRATION(mitkFiberfoxSignalGeneration)