diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkProjectionFilter.txx b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkProjectionFilter.txx index deb59b26b4..39228e4f76 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkProjectionFilter.txx +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkProjectionFilter.txx @@ -1,202 +1,202 @@ /*=================================================================== 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 _itkProjectionFilter_txx #define _itkProjectionFilter_txx #include "itkProjectionFilter.h" #include "mitkProgressBar.h" -#include +//#include #define SEARCHSIGMA 10 /* length in linear voxel dimensions */ #define MAXSEARCHLENGTH (3*SEARCHSIGMA) namespace itk { ProjectionFilter::ProjectionFilter() { } ProjectionFilter::~ProjectionFilter() { } void ProjectionFilter::Project() { // Contains only code for the projection of FA data. The original FSL code contains some extra lines // For projection of other measurements than FA mitk::ProgressBar::GetInstance()->AddStepsToDo( 3 ); Float4DImageType::Pointer data_4d_projected = Float4DImageType::New(); data_4d_projected->SetRegions(m_AllFA->GetRequestedRegion()); data_4d_projected->SetDirection(m_AllFA->GetDirection()); data_4d_projected->SetSpacing(m_AllFA->GetSpacing()); data_4d_projected->SetOrigin(m_AllFA->GetOrigin()); data_4d_projected->Allocate(); data_4d_projected->FillBuffer(0.0); Float4DImageType::SizeType size = m_AllFA->GetRequestedRegion().GetSize(); for(int t=0; tGetPixel(ix) != 0) { VectorImageType::PixelType dir = m_Directions->GetPixel(ix); short maxvalX=0, maxvalY=0, maxvalZ=0; Float4DImageType::IndexType ix4d; ix4d[0]=x; ix4d[1]=y; ix4d[2]=z; ix4d[3]=t; float maxval = m_AllFA->GetPixel(ix4d); float maxval_weighted = maxval; float exponentfactor = -0.5 * (dir[0]*dir[0]+dir[1]*dir[1]+dir[2]*dir[2]) / (float)(SEARCHSIGMA*SEARCHSIGMA); // No tubular structure here if(m_Tube->GetPixel(ix) == 0) { for(int iters=0;iters<2;iters++) { float distance=0; for(int d=1;d=size[0] && dy<=size[1] && dz<=size[2]) { d=MAXSEARCHLENGTH; } else if(m_DistanceMap->GetPixel(ix3d)>=distance) { float distanceweight = exp(d * d * exponentfactor); distance = m_DistanceMap->GetPixel(ix3d); ix4d[0]=x+dir[0]*D; ix4d[1]=y+dir[1]*D; ix4d[2]=z+dir[2]*D; ix4d[3]=t; if(distanceweight*m_AllFA->GetPixel(ix4d)>maxval_weighted) { maxval = m_AllFA->GetPixel(ix4d); maxval_weighted = maxval*distanceweight; maxvalX = dir[0]*D; maxvalY = dir[1]*D; maxvalZ = dir[2]*D; } } else{ d=MAXSEARCHLENGTH; } } } } // Tubular structure else { for(int dy=-MAXSEARCHLENGTH; dy<=MAXSEARCHLENGTH;dy++) { for(int dx=-MAXSEARCHLENGTH; dx<=MAXSEARCHLENGTH; dx++) { float distanceweight = exp(-0.5 * (dx*dx+dy*dy) / (float)(SEARCHSIGMA*SEARCHSIGMA) ); float r=sqrt((float)(dx*dx+dy*dy)); if (r>0) { int allok=1; for(float rr=1; rr<=r+0.1; rr++) /* search outwards from centre to current voxel - test that distancemap always increasing */ { int dx1=round(rr*dx/r); int dy1=round(rr*dy/r); int dx2=round((rr+1)*dx/r); int dy2=round((rr+1)*dy/r); RealImageType::IndexType ix1, ix2; ix1[0]=x+dx1; ix1[1]=y+dy1; ix1[2]=z; ix2[0]=x+dx2; ix2[1]=y+dy2; ix2[2]=z; if(m_DistanceMap->GetPixel(ix1) > m_DistanceMap->GetPixel(ix2) ) { allok=0; } } ix4d[0]=x+dx; ix4d[1]=y+dy, ix4d[2]=z; ix4d[3]=t; if( allok && (distanceweight * m_AllFA->GetPixel(ix4d) > maxval_weighted) ) { maxval = m_AllFA->GetPixel(ix4d); maxval_weighted = maxval * distanceweight; maxvalX=dx; maxvalY=dy; maxvalZ=0; } } } } } ix4d[0]=x; ix4d[1]=y; ix4d[2]=z; ix4d[3]=t; data_4d_projected->SetPixel(ix4d, maxval); } } } } } m_Projections = data_4d_projected; } } #endif // _itkProjectionFilter_txx