diff --git a/Modules/DiffusionImaging/Algorithms/itkDistanceMapFilter.h b/Modules/DiffusionImaging/Algorithms/itkDistanceMapFilter.h
new file mode 100644
index 0000000000..40713ee910
--- /dev/null
+++ b/Modules/DiffusionImaging/Algorithms/itkDistanceMapFilter.h
@@ -0,0 +1,103 @@
+/*=========================================================================
+
+Program:   Medical Imaging & Interaction Toolkit
+Module:    $RCSfile$
+Language:  C++
+Date:      $Date: 2007-03-19 18:26:32 +0100 (Mo, 19 Mär 2007) $
+Version:   $Revision: 9819 $
+
+Copyright (c) German Cancer Research Center, Division of Medical and
+Biological Informatics. All rights reserved.
+See MITKCopyright.txt or http://www.mitk.org/copyright.html for details.
+
+This software is distributed WITHOUT ANY WARRANTY; without even
+the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+PURPOSE.  See the above copyright notices for more information.
+
+=========================================================================*/
+
+#ifndef ITKDISTANCEMAPFILTER_H_
+#define ITKDISTANCEMAPFILTER_H_
+
+#include "itkImageToImageFilter.h"
+#include "itkImage.h"
+#include "mitkImage.h"
+
+
+namespace itk
+{
+template < class TInputImage, class TOutputImage >
+class DistanceMapFilter : public ImageToImageFilter<TInputImage, TOutputImage>
+{
+	/*!
+	\brief itkDistanceMapFilter 
+	
+	\brief Creates a distance map from a FA skeleton image.
+	
+	\sa itkImageToImageFilter
+	
+	
+	\verbatim
+	Last contributor: $Author: vanbrugg $
+	\endverbatim
+	*/
+
+public:
+
+	/** Typedef for input ImageType. */
+  typedef TInputImage  InputImageType;
+
+  /** Typedef for input imageType Pointer. */
+  typedef typename InputImageType::Pointer InputImagePointer;
+
+  /** Typedef for output ImageType. */
+  typedef TOutputImage OutputImageType;
+
+  /** Typedef for input imageType Pointer. */
+  typedef typename OutputImageType::Pointer OutputImagePointer;
+
+public:
+
+	/** */
+  typedef DistanceMapFilter Self;
+  
+  /** Superclass */
+  typedef ImageToImageFilter<InputImageType, OutputImageType> Superclass;
+  
+  /** Smart Pointer */
+  typedef SmartPointer<Self> Pointer;
+  
+  /** Smart Pointer */
+  typedef SmartPointer<const Self> ConstPointer;
+
+	/** */
+  itkNewMacro( Self);
+
+	/** Generate Data. The image will be divided into a number of pieces, a number of threads 
+  will be spawned and Threaded GenerateData() will be called in each thread. */
+  virtual void GenerateData();
+
+ 
+protected:
+
+	/** Constructor */
+  DistanceMapFilter();
+  
+  /** Destructor */
+  virtual ~DistanceMapFilter();
+
+  
+  
+protected:
+ 	
+  	
+  
+};
+
+}
+
+#ifndef ITK_MANUAL_INSTANTIATION
+#include "itkDistanceMapFilter.txx"
+#endif
+
+#endif 
diff --git a/Modules/DiffusionImaging/Algorithms/itkDistanceMapFilter.txx b/Modules/DiffusionImaging/Algorithms/itkDistanceMapFilter.txx
new file mode 100644
index 0000000000..3d29cf2814
--- /dev/null
+++ b/Modules/DiffusionImaging/Algorithms/itkDistanceMapFilter.txx
@@ -0,0 +1,87 @@
+/*=========================================================================
+
+Program:   Medical Imaging & Interaction Toolkit
+Module:    $RCSfile$
+Language:  C++
+Date:      $ $
+Version:   $ $
+
+Copyright (c) German Cancer Research Center, Division of Medical and
+Biological Informatics. All rights reserved.
+See MITKCopyright.txt or http://www.mitk.org/copyright.html for details.
+
+This software is distributed WITHOUT ANY WARRANTY; without even
+the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+PURPOSE.  See the above copyright notices for more information.
+
+=========================================================================*/
+#ifndef _itkDistanceMapFilter_txx
+#define _itkDistanceMapFilter_txx
+
+#include "itkDistanceMapFilter.h"
+
+#include "mitkProgressBar.h"
+#include <itkSignedMaurerDistanceMapImageFilter.h>
+
+
+namespace itk
+{
+
+  template< class TInputImage, class TOutputImage >
+  DistanceMapFilter<TInputImage, TOutputImage>::DistanceMapFilter()
+  {
+   
+  }
+
+  template< class TInputImage, class TOutputImage >
+  DistanceMapFilter<TInputImage, TOutputImage>::~DistanceMapFilter()
+  {
+    
+
+  }
+
+ 
+
+ 
+  template< class TInputImage, class TOutputImage >
+  void DistanceMapFilter<TInputImage, TOutputImage>::GenerateData()
+  {
+    //----------------------------------------------------------------------//
+    //  Progress bar                                                        //
+    //----------------------------------------------------------------------//
+    //mitk::ProgressBar::GetInstance()->AddStepsToDo( 3 );
+
+   
+
+    typedef itk::SignedMaurerDistanceMapImageFilter<InputImageType, OutputImageType> DistanceFilterType;
+    typename DistanceFilterType::Pointer dFilter = DistanceFilterType::New();
+    dFilter->SetInput(this->GetInput());
+    dFilter->SetUseImageSpacing(true);
+    dFilter->SetSquaredDistance(false);
+    dFilter->SetInsideIsPositive(true);
+    dFilter->Update();
+
+    typename OutputImageType::Pointer outputImg = dFilter->GetOutput();
+
+    typedef itk::ImageRegionIterator<OutputImageType> ImageIteratorType;
+
+    ImageIteratorType outIt(outputImg, outputImg->GetRequestedRegion());
+    outIt.GoToBegin();
+    while(!outIt.IsAtEnd())
+    {
+      typename OutputImageType::PixelType p = outIt.Get();
+      p *= -1;
+      outIt.Set(p);
+      ++outIt;
+    }
+
+    Superclass::SetNthOutput( 0, outputImg );
+
+  }
+
+
+
+  
+  
+}
+#endif // _itkDistanceMapFilter_txx
\ No newline at end of file
diff --git a/Modules/DiffusionImaging/Algorithms/itkProjectionFilter.h b/Modules/DiffusionImaging/Algorithms/itkProjectionFilter.h
new file mode 100644
index 0000000000..f924d200fa
--- /dev/null
+++ b/Modules/DiffusionImaging/Algorithms/itkProjectionFilter.h
@@ -0,0 +1,131 @@
+/*=========================================================================
+
+Program:   Medical Imaging & Interaction Toolkit
+Module:    $RCSfile$
+Language:  C++
+Date:      $Date: 2007-03-19 18:26:32 +0100 (Mo, 19 Mär 2007) $
+Version:   $Revision: 9819 $
+
+Copyright (c) German Cancer Research Center, Division of Medical and
+Biological Informatics. All rights reserved.
+See MITKCopyright.txt or http://www.mitk.org/copyright.html for details.
+
+This software is distributed WITHOUT ANY WARRANTY; without even
+the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+PURPOSE.  See the above copyright notices for more information.
+
+=========================================================================*/
+
+#ifndef ITKPROJECTIONFILTER_H_
+#define ITKPROJECTIONFILTER_H_
+
+#include "itkObject.h"
+#include "itkImage.h"
+#include "mitkImage.h"
+#include "mitkTbssImage.h"
+
+namespace itk
+{
+class ProjectionFilter : public Object
+{
+	/*!
+	\brief itkProjectionFilter 
+	
+	\brief Performs the Projection Step of the TBSS algorithm from Smith et al. 2006.
+	
+	
+	*/
+
+public:	
+
+  typedef itk::Image<float, 3> RealImageType;
+
+  typedef itk::CovariantVector<int,3> VectorType;
+
+  typedef itk::Image<VectorType, 3> VectorImageType;
+
+  typedef itk::Image<char,3> CharImageType;
+
+  typedef itk::Image<float, 4> Float4DImageType;
+
+ // typedef itk::VectorImageType<float, 3> SubjectVolumesType;
+
+
+public:
+
+	/** */
+  typedef ProjectionFilter Self;
+  
+  /** Superclass */
+  typedef Object Superclass;
+  
+  /** Smart Pointer */
+  typedef SmartPointer<Self> Pointer;
+  
+  /** Smart Pointer */
+  typedef SmartPointer<const Self> ConstPointer;
+
+	/** */
+  itkNewMacro( Self);
+
+	/** Generate Data. The image will be divided into a number of pieces, a number of threads 
+  will be spawned and Threaded GenerateData() will be called in each thread. */
+  void Project();
+
+  itkSetMacro(DistanceMap, RealImageType::Pointer);
+
+  itkSetMacro(Directions, VectorImageType::Pointer);
+
+  itkSetMacro(Skeleton, CharImageType::Pointer);
+
+  itkSetMacro(Tube, CharImageType::Pointer);
+
+  itkSetMacro(AllFA, Float4DImageType::Pointer);
+
+  itkGetMacro(Projections, Float4DImageType::Pointer);
+
+  itkGetMacro(Projected, RealImageType::Pointer);
+
+
+ 
+protected:
+
+	/** Constructor */
+  ProjectionFilter();
+  
+  /** Destructor */
+  virtual ~ProjectionFilter();
+
+  RealImageType::Pointer m_DistanceMap;
+
+  VectorImageType::Pointer m_Directions;
+
+  CharImageType::Pointer m_Skeleton;
+
+  CharImageType::Pointer m_Tube;  
+
+  Float4DImageType::Pointer m_Projections;
+
+  Float4DImageType::Pointer m_AllFA;
+  
+  RealImageType::Pointer m_Projected;
+
+  int round(float x) 
+  { 
+    if (x>0.0) return ((int) (x+0.5));
+    else       return ((int) (x-0.5));
+  }
+  
+protected:
+ 	
+  	
+  
+};
+
+}
+
+#ifndef ITK_MANUAL_INSTANTIATION
+#include "itkProjectionFilter.txx"
+#endif
+
+#endif 
diff --git a/Modules/DiffusionImaging/Algorithms/itkProjectionFilter.txx b/Modules/DiffusionImaging/Algorithms/itkProjectionFilter.txx
new file mode 100644
index 0000000000..05e487dd3f
--- /dev/null
+++ b/Modules/DiffusionImaging/Algorithms/itkProjectionFilter.txx
@@ -0,0 +1,204 @@
+/*=========================================================================
+
+Program:   Medical Imaging & Interaction Toolkit
+Module:    $RCSfile$
+Language:  C++
+Date:      $ $
+Version:   $ $
+
+Copyright (c) German Cancer Research Center, Division of Medical and
+Biological Informatics. All rights reserved.
+See MITKCopyright.txt or http://www.mitk.org/copyright.html for details.
+
+This software is distributed WITHOUT ANY WARRANTY; without even
+the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+PURPOSE.  See the above copyright notices for more information.
+
+=========================================================================*/
+#ifndef _itkProjectionFilter_txx
+#define _itkProjectionFilter_txx
+
+#include "itkProjectionFilter.h"
+
+#include "mitkProgressBar.h"
+#include <itkSignedMaurerDistanceMapImageFilter.h>
+
+#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; t<size[3]; t++)
+    {
+      for(int z=1; z<size[2]-1; z++)
+      {
+        for(int y=1; y<size[1]-1; y++)
+        {
+
+          std::cout << "x";
+          for(int x=1; x<size[0]-1; x++)
+          {           
+
+            VectorImageType::IndexType ix; 
+            ix[0] = x; ix[1] = y; ix[2] = z;
+
+            if(m_Skeleton->GetPixel(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), maxval_weighted = maxval,
+                exponentfactor = -0.5 * (dir[0]*dir[0]+dir[1]*dir[1]+dir[2]*dir[2]) / (float)(SEARCHSIGMA*SEARCHSIGMA);
+
+
+              // No cingulum here
+              if(m_Tube->GetPixel(ix) == 0)
+              {
+                for(int iters=0;iters<2;iters++)
+                {
+                  float distance=0;
+
+                  for(int d=1;d<MAXSEARCHLENGTH;d++)
+                  {
+                    int D=d;
+                    if (iters==1) D=-d;
+                  
+                    FloatImageType::IndexType ix3d;
+                    int dx = x+dir[0]*D, dy = y+dir[1]*D, dz = z+dir[2]*D;
+                    ix3d[0] = dx; ix3d[1] = dy; ix3d[2] = dz;
+                    
+                    if(dx<0 || dy<0 || dz<0
+                      || dx>=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; 
+                    }       
+
+                    
+                  
+                  } // endfor(int d=1;d<MAXSEARCHLENGTH;d++)
+                
+                } // endfor(int iters=0;iters<2;iters++) 
+              } // endif (m_Tube->GetPixel(ix) == 0)
+
+
+              // Cingulum here
+              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;
+                      }
+
+
+                    } //endif(r>0)
+
+                  } //endfor(int xxx=-MAXSEARCHLENGTH; xxx<=MAXSEARCHLENGTH; xxx++)
+                } // endfor(int dy=-MAXSEARCHLENGTH; dy<=MAXSEARCHLENGTH;y++)
+              } // endelse 
+
+              ix4d[0]=x; ix4d[1]=y; ix4d[2]=z; ix4d[3]=t; 
+              data_4d_projected->SetPixel(ix4d, maxval);
+
+            } // endif(m_Skeleton->GetPixel(x) != 0) 
+          } // endfor(int x=1; x<size[0]-1; x++)
+        } // endfor(int y=1; y<size[1]-1; y++)
+      } // endfor(int z=1; z<size[2]-1; z++)
+    } // endfor(int t=0; t<size[3]; t++)
+  
+    m_Projections = data_4d_projected;
+
+   
+
+  } // GenerateData()
+
+
+
+  
+  
+}
+#endif // _itkProjectionFilter_txx
diff --git a/Modules/DiffusionImaging/Algorithms/itkShortestPathCostFunctionTbss.txx b/Modules/DiffusionImaging/Algorithms/itkShortestPathCostFunctionTbss.txx
new file mode 100644
index 0000000000..32896d160d
--- /dev/null
+++ b/Modules/DiffusionImaging/Algorithms/itkShortestPathCostFunctionTbss.txx
@@ -0,0 +1,74 @@
+#ifndef __itkShortestPathCostFunctionTbss_txx
+#define __itkShortestPathCostFunctionTbss_txx
+
+#include "itkShortestPathCostFunctionTbss.h"
+#include <limits>
+
+
+
+namespace itk
+{ 
+
+  // Constructor
+  template<class TInputImageType>
+  ShortestPathCostFunctionTbss<TInputImageType>
+    ::ShortestPathCostFunctionTbss()
+  {
+  }
+
+
+ template<class TInputImageType>
+  double ShortestPathCostFunctionTbss<TInputImageType>
+    ::GetCost(IndexType p1 ,IndexType  p2)
+  {
+    // Very simple Metric:
+    // The squared difference of both values is defined as cost
+
+    double a,b,c;    
+    ConstIteratorType it( this->m_Image, this->m_Image->GetRequestedRegion());
+    
+    it.SetIndex(p1);
+    a = it.Get();
+    it.SetIndex(p2);
+    b = it.Get();
+
+
+    if(a==b)
+    {     
+      double dxSqt = (p1[0]-p2[0]) * (p1[0]-p2[0]);// * 1000;
+      double dySqt = (p1[1]-p2[1]) * (p1[1]-p2[1]);
+      double dzSqt = (p1[2]-p2[2]) * (p1[2]-p2[2]);
+      c = sqrt(dxSqt + dySqt + dzSqt);
+    }
+    else if(p1[0] != p2[0])
+    {
+      c=std::numeric_limits<double>::max( );
+    }
+    else
+    {
+      c=std::numeric_limits<double>::max( );
+    }
+
+    return c; 
+  }
+
+
+  template<class TInputImageType>
+  void ShortestPathCostFunctionTbss<TInputImageType>
+    ::Initialize()
+  {
+  }
+
+  template<class TInputImageType>
+  double ShortestPathCostFunctionMitral<TInputImageType>
+    ::GetMinCost()
+  {
+    return 1.0;
+  }
+
+
+
+
+} // end namespace itk
+
+#endif // __itkShortestPathCostFunctionSimple_txx
diff --git a/Modules/DiffusionImaging/Algorithms/itkSkeletonizationFilter.h b/Modules/DiffusionImaging/Algorithms/itkSkeletonizationFilter.h
new file mode 100644
index 0000000000..6323bb0f7b
--- /dev/null
+++ b/Modules/DiffusionImaging/Algorithms/itkSkeletonizationFilter.h
@@ -0,0 +1,118 @@
+/*=========================================================================
+
+Program:   Medical Imaging & Interaction Toolkit
+Module:    $RCSfile$
+Language:  C++
+Date:      $Date: 2007-03-19 18:26:32 +0100 (Mo, 19 Mär 2007) $
+Version:   $Revision: 9819 $
+
+Copyright (c) German Cancer Research Center, Division of Medical and
+Biological Informatics. All rights reserved.
+See MITKCopyright.txt or http://www.mitk.org/copyright.html for details.
+
+This software is distributed WITHOUT ANY WARRANTY; without even
+the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+PURPOSE.  See the above copyright notices for more information.
+
+=========================================================================*/
+
+#ifndef ITKSKELETONIZATIONFILTER_H_
+#define ITKSKELETONIZATIONFILTER_H_
+
+#include "itkImageToImageFilter.h"
+#include "itkImage.h"
+#include "mitkImage.h"
+#include <itkVectorImage.h>
+
+namespace itk
+{
+template < class TInputImage, class TOutputImage >
+class SkeletonizationFilter : public ImageToImageFilter<TInputImage, TOutputImage>
+{
+	/*!
+	\brief itkSkeletonizationFilter
+	
+	\brief Skeletonization algorithm from TBSS (Smith et al. 2006)
+	
+	\sa itkImageToImageFilter
+	
+	
+	*/
+
+public:
+
+	/** Typedef for input ImageType. */
+  typedef TInputImage  InputImageType;
+
+  typedef itk::CovariantVector<int,3> VectorType;
+
+  typedef itk::Image<VectorType, 3> VectorImageType;
+
+
+  /** Typedef for output ImageType. */
+  typedef TOutputImage OutputImageType;  
+
+  typedef itk::VectorImage<int, 3> GradientImageType;
+
+
+	/** */
+  typedef SkeletonizationFilter Self;
+  
+  /** Superclass */
+  typedef ImageToImageFilter<InputImageType, OutputImageType> Superclass;
+  
+  /** Smart Pointer */
+  typedef SmartPointer<Self> Pointer;
+  
+  /** Smart Pointer */
+  typedef SmartPointer<const Self> ConstPointer;
+
+	/** */
+  itkNewMacro( Self);
+
+	/** Generate Data. The image will be divided into a number of pieces, a number of threads 
+  will be spawned and Threaded GenerateData() will be called in each thread. */
+  virtual void GenerateData();
+
+  
+  GradientImageType::Pointer GetDirectionImage();
+
+
+
+
+  
+
+protected:
+
+	/** Constructor */
+  SkeletonizationFilter();
+  
+  /** Destructor */
+  virtual ~SkeletonizationFilter();
+
+  void CalculatePerpendicularDirections();
+
+  VectorImageType::Pointer m_DirectionImage;
+
+  //FloatVectorImageType::Pointer m_FixedDirImage;
+
+  int round(float x) 
+  { 
+    if (x>0.0) return ((int) (x+0.5));
+    else       return ((int) (x-0.5));
+  }
+  
+  
+protected:
+ 	
+  	
+  
+};
+
+}
+
+#ifndef ITK_MANUAL_INSTANTIATION
+#include "itkSkeletonizationFilter.txx"
+#endif
+
+#endif
diff --git a/Modules/DiffusionImaging/Algorithms/itkSkeletonizationFilter.txx b/Modules/DiffusionImaging/Algorithms/itkSkeletonizationFilter.txx
new file mode 100644
index 0000000000..01da799365
--- /dev/null
+++ b/Modules/DiffusionImaging/Algorithms/itkSkeletonizationFilter.txx
@@ -0,0 +1,314 @@
+/*=========================================================================
+
+Program:   Medical Imaging & Interaction Toolkit
+Module:    $RCSfile$
+Language:  C++
+Date:      $ $
+Version:   $ $
+
+Copyright (c) German Cancer Research Center, Division of Medical and
+Biological Informatics. All rights reserved.
+See MITKCopyright.txt or http://www.mitk.org/copyright.html for details.
+
+This software is distributed WITHOUT ANY WARRANTY; without even
+the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+PURPOSE.  See the above copyright notices for more information.
+
+=========================================================================*/
+#ifndef _itkSkeletonizationFilter_txx
+#define _itkSkeletonizationFilter_txx
+
+#include "itkSkeletonizationFilter.h"
+#include "mitkProgressBar.h"
+#include <limits>
+#include <time.h> 
+
+namespace itk
+{
+
+  template< class TInputImage, class TOutputImage >
+  SkeletonizationFilter<TInputImage, TOutputImage>::SkeletonizationFilter()
+  {
+    m_DirectionImage = VectorImageType::New();
+  }
+
+  template< class TInputImage, class TOutputImage >
+  SkeletonizationFilter<TInputImage, TOutputImage>::~SkeletonizationFilter()
+  {    
+
+  }
+
+ 
+
+ 
+  template< class TInputImage, class TOutputImage >
+  void SkeletonizationFilter<TInputImage, TOutputImage>::GenerateData()
+  {
+    //----------------------------------------------------------------------//
+    //  Progress bar                                                        //
+    //----------------------------------------------------------------------//
+    mitk::ProgressBar::GetInstance()->AddStepsToDo( 3 );    
+
+    std::cout << "Skeletonize" << std::endl;
+    CalculatePerpendicularDirections();
+  }
+
+
+
+  template< class TInputImage, class TOutputImage >
+  void SkeletonizationFilter<TInputImage, TOutputImage>::CalculatePerpendicularDirections()
+  {      
+
+    const InputImageType* faImage = this->GetInput();
+    typename InputImageType::SizeType size = faImage->GetRequestedRegion().GetSize();
+
+    //typename RealImageType::SizeType size = m_FaImage->GetRequestedRegion().GetSize();
+
+    m_DirectionImage->SetRegions(faImage->GetRequestedRegion());
+    m_DirectionImage->SetDirection(faImage->GetDirection());
+    m_DirectionImage->SetSpacing(faImage->GetSpacing());    
+    m_DirectionImage->SetOrigin(faImage->GetOrigin());
+    m_DirectionImage->Allocate();
+    m_DirectionImage->FillBuffer(0.0);
+    
+
+    for(int z=1; z<size[2]-1; z++) for(int y=1; y<size[1]-1; y++) for(int x=1; x<size[0]-1; x++)
+    { 
+      typename InputImageType::IndexType ix;
+      ix[0]=x; ix[1]=y; ix[2]=z;
+      float theval = faImage->GetPixel(ix);
+
+      if(theval != 0)
+      {      
+
+        /* Calculate point of gravity. We will consider each 3x3x3 neighbourhood as a unit cube. The center
+         * point of each voxel will be a multiplicative of 1/6. The center of the unit cube is 3/6 = 1/2/
+         */
+        
+        float cogX = 0.0; float cogY = 0.0; float cogZ = 0.0; float sum = 0.0; float l;
+        int vecX = 0; int vecY = 0; int vecZ = 0;
+        
+        for(int dz=-1; dz<=1; dz++) for(int dy=-1; dy<=1; dy++) for(int dx=-1; dx<=1;dx++)
+        {
+          typename InputImageType::IndexType p;
+          p[0] = x+dx; p[1] = y+dy; p[2] = z+dz;
+          float mass = faImage->GetPixel(p);
+
+          sum += mass;
+          cogX += (float)dx*mass; cogY += (float)dy*mass; cogZ += (float)dz*mass;           
+        }
+
+        cogX /= sum; cogY /= sum; cogZ /= sum;
+        l = sqrt(cogX*cogX + cogY*cogY + cogZ*cogZ);
+
+        if (l > 0.1) /* is CofG far enough away from centre voxel? */
+        {
+          vecX = std::max(std::min(round(cogX/l),1),-1);
+          vecY = std::max(std::min(round(cogY/l),1),-1);
+          vecZ = std::max(std::min(round(cogZ/l),1),-1);          
+        }
+        else
+          // Find direction of max curvature
+        { 
+          
+          float maxcost=0, centreval=2*theval;
+          for(int zz=0; zz<=1; zz++) // note - starts at zero as we're only searching half the voxels 
+          {
+            for(int yy=-1; yy<=1; yy++)
+            {
+              for(int xx=-1; xx<=1; xx++)
+              {
+                if ( (zz==1) || (yy==1) || ((yy==0)&&(xx==1)) )
+                {
+                  float weighting = pow( (float)(xx*xx+yy*yy+zz*zz) , (float)-0.7 ); // power is arbitrary: maybe test other functions here 
+                  
+
+                  typename InputImageType::IndexType i,j;
+                  i[0] = x+xx; i[1] = y+yy; i[2] = z+zz;
+                  j[0] = x-xx; j[1] = y-yy; j[2] = z-zz;
+                  float cost = weighting * ( centreval 
+                    - (float)faImage->GetPixel(i)
+                    - (float)faImage->GetPixel(j));
+
+                  if (cost>maxcost)
+                  {
+                    maxcost=cost;
+                    vecX=xx;
+                    vecY=yy;
+                    vecZ=zz;
+                  }
+                }
+              }
+            }
+          }         
+        }
+
+        VectorType vec;
+        vec[0] = vecX; vec[1] = vecY; vec[2]=vecZ;
+        m_DirectionImage->SetPixel(ix, vec);        
+      
+      }     
+    }
+
+    mitk::ProgressBar::GetInstance()->Progress();
+
+
+    // Smooth m_DirectionImage and store in directionSmoothed by finding the
+    // mode in a 3*3 neighbourhoud
+    VectorImageType::Pointer directionSmoothed = VectorImageType::New();
+    directionSmoothed->SetRegions(faImage->GetRequestedRegion());
+    directionSmoothed->SetDirection(faImage->GetDirection());
+    directionSmoothed->SetSpacing(faImage->GetSpacing());    
+    directionSmoothed->SetOrigin(faImage->GetOrigin());
+    directionSmoothed->Allocate();
+
+    VectorImageType::PixelType p;
+    p[0]=0; p[1]=0; p[2]=0; 
+    directionSmoothed->FillBuffer(p);
+
+  
+
+    for(int z=1; z<size[2]-1; z++) for(int y=1; y<size[1]-1; y++) for(int x=1; x<size[0]-1; x++)
+    {   
+      VectorImageType::IndexType ix;
+      ix[0]=x; ix[1]=y; ix[2]=z;  
+
+      // Find the vector that occured most
+      int* localsum = new int[27];
+      int localmax=0, xxx, yyy, zzz;
+
+      for(int zz=0; zz<27; zz++) localsum[zz]=0;
+
+      for(int zz=-1; zz<=1; zz++) for(int yy=-1; yy<=1; yy++) for(int xx=-1; xx<=1; xx++)
+      {
+        VectorImageType::IndexType i;
+        i[0] = x+xx; i[1] = y+yy; i[2] = z+zz;
+        VectorType v = m_DirectionImage->GetPixel(i);
+        xxx = v[0];
+        yyy = v[1];
+        zzz = v[2];
+       
+        localsum[(1+zzz)*9+(1+yyy)*3+1+xxx]++;
+        localsum[(1-zzz)*9+(1-yyy)*3+1-xxx]++;
+      }
+
+      for(int zz=-1; zz<=1; zz++) for(int yy=-1; yy<=1; yy++) for(int xx=-1; xx<=1; xx++)
+      {
+        if (localsum[(1+zz)*9+(1+yy)*3+1+xx]>localmax)
+        {
+          localmax=localsum[(1+zz)*9+(1+yy)*3+1+xx];
+          VectorType v;
+          v[0] = xx; v[1] = yy; v[2] = zz;
+          directionSmoothed->SetPixel(ix, v);
+        }
+      }
+
+      delete localsum;
+      
+    }
+
+    m_DirectionImage = directionSmoothed;
+
+    mitk::ProgressBar::GetInstance()->Progress();
+
+    // Do non-max-suppression in the direction of perp and set as output of the filter
+    typename OutputImageType::Pointer outputImg = OutputImageType::New();
+    outputImg->SetRegions(faImage->GetRequestedRegion());
+    outputImg->SetDirection(faImage->GetDirection());
+    outputImg->SetSpacing(faImage->GetSpacing());    
+    outputImg->SetOrigin(faImage->GetOrigin());
+    outputImg->Allocate();
+    outputImg->FillBuffer(0.0);  
+
+    for(int z=1; z<size[2]-1; z++) for(int y=1; y<size[1]-1; y++) for(int x=1; x<size[0]-1; x++)
+    {
+
+      typename InputImageType::IndexType ix;
+      ix[0]=x; ix[1]=y; ix[2]=z;
+      
+      float theval = faImage->GetPixel(ix);
+      VectorType v = directionSmoothed->GetPixel(ix);
+      
+      typename VectorImageType::IndexType i;
+
+      i[0] = x-v[0]; i[1] = y-v[1]; i[2] = z-v[2];
+      float min = faImage->GetPixel(i);
+
+      i[0] = x+v[0]; i[1] = y+v[1]; i[2] = z+v[2];
+      float plus = faImage->GetPixel(i);
+
+      i[0] = x-2*v[0]; i[1] = y-2*v[1]; i[2] = z-2*v[2];
+      float minmin = faImage->GetPixel(i);
+
+      i[0] = x+2*v[0]; i[1] = y+2*v[1]; i[2] = z+2*v[2];
+      float plusplus = faImage->GetPixel(i);
+
+      if( ((v[0]!=0) || (v[1]!=0) || (v[2]!=0)) &&
+           theval >= plus && theval >  min && theval >= plusplus && theval >  minmin  )
+      {
+        outputImg->SetPixel(ix, theval);
+      }
+     
+    }
+   
+    
+
+    Superclass::SetNthOutput( 0, outputImg );
+    mitk::ProgressBar::GetInstance()->Progress();
+   
+
+  }
+
+  
+  template< class TInputImage, class TOutputImage >
+  itk::VectorImage<int, 3>::Pointer SkeletonizationFilter<TInputImage, TOutputImage>::GetDirectionImage()
+  {
+    GradientImageType::Pointer gradImg = GradientImageType::New();
+
+    if(m_DirectionImage.IsNotNull())
+    {
+      gradImg->SetSpacing(m_DirectionImage->GetSpacing());
+      gradImg->SetOrigin(m_DirectionImage->GetOrigin());
+      gradImg->SetDirection(m_DirectionImage->GetDirection());
+      gradImg->SetRegions(m_DirectionImage->GetLargestPossibleRegion().GetSize());
+      gradImg->SetVectorLength(3);
+      gradImg->Allocate();
+
+
+      VectorImageType::SizeType size = m_DirectionImage->GetLargestPossibleRegion().GetSize();
+
+      for(int i=0; i<size[0]; i++)
+      {
+        for(int j=0; j<size[1]; j++)
+        {
+          for(int k=0; k<size[2]; k++)
+          {
+            itk::Index<3> ix;
+            ix[0] = i;
+            ix[1] = j;
+            ix[2] = k;
+
+            VectorType vec = m_DirectionImage->GetPixel(ix);
+
+            itk::VariableLengthVector<int> pixel;
+            pixel.SetSize(3);
+            pixel.SetElement(0, vec.GetElement(0));
+            pixel.SetElement(1, vec.GetElement(1));
+            pixel.SetElement(2, vec.GetElement(2));
+
+            gradImg->SetPixel(ix, pixel);
+
+          }
+        }
+      }
+
+    }
+
+
+
+
+    return gradImg;
+  }
+  
+}
+#endif // _itkSkeletonizationFilter_txx
diff --git a/Modules/DiffusionImaging/Algorithms/mitkTractAnalyzer.cpp b/Modules/DiffusionImaging/Algorithms/mitkTractAnalyzer.cpp
new file mode 100644
index 0000000000..d5e9be0f60
--- /dev/null
+++ b/Modules/DiffusionImaging/Algorithms/mitkTractAnalyzer.cpp
@@ -0,0 +1,162 @@
+
+#ifndef __mitkTractAnalyzer_cpp
+#define __mitkTractAnalyzer_cpp
+
+
+
+#include <mitkLogMacros.h>
+#include <mitkTractAnalyzer.h>
+
+#include <itkImageRegionIteratorWithIndex.h>
+
+
+#include <itkShortestPathCostFunctionTbss.h>
+#include <itkShortestPathImageFilter.h>
+
+#include <iostream>
+#include <fstream>
+
+
+using namespace std;
+
+namespace mitk {
+
+  TractAnalyzer::TractAnalyzer()
+  {
+
+  }
+
+
+  void TractAnalyzer::BuildGraph(itk::Index<3> startPoint, itk::Index<3> endPoint)
+  {
+
+    typedef itk::ShortestPathImageFilter<FloatImageType, CharImageType> ShortestPathFilterType;
+    typedef itk::ShortestPathCostFunctionTbss<FloatImageType> CostFunctionType;
+
+
+    FloatImageType::Pointer meanSkeleton;
+
+    mitk::CastToItkImage(m_InputImage, meanSkeleton);
+
+  // Only use the mitk image
+
+
+
+    if(meanSkeleton)
+    {
+      CostFunctionType::Pointer costFunction = CostFunctionType::New();
+      costFunction->SetImage(meanSkeleton);
+      costFunction->SetStartIndex(startPoint);
+      costFunction->SetEndIndex(endPoint);
+      costFunction->SetThreshold(m_Threshold);
+
+      ShortestPathFilterType::Pointer pathFinder = ShortestPathFilterType::New();
+      pathFinder->SetCostFunction(costFunction);
+      pathFinder->SetFullNeighborsMode(true);
+      //pathFinder->SetCalcMode(ShortestPathFilterType::A_STAR);
+      pathFinder->SetInput(meanSkeleton);
+      pathFinder->SetStartIndex(startPoint);
+      pathFinder->SetEndIndex(endPoint);
+      pathFinder->Update();
+
+      m_Path = pathFinder->GetVectorPath();
+
+
+      m_RoiImg = pathFinder->GetOutput();
+    }
+
+
+
+  }
+
+  void TractAnalyzer::MeasureRoi()
+  {
+
+    // Output two types
+    ProjectionsImageType::SizeType size = m_Projections->GetLargestPossibleRegion().GetSize();
+
+    std::ofstream file(m_FileName.c_str());
+
+    std::vector<std::string> individuals;
+    for(int i=0; i<m_Groups.size(); i++)
+    {
+      std::pair<std::string, int> group = m_Groups[i];
+      for(int j=0; j<group.second; j++)
+      {
+        individuals.push_back(group.first);
+      }
+    }
+
+
+    // Position ID
+
+    file << "position ";
+
+    for(int i=0; i<m_Roi.size(); i++)
+    {
+      file << "pos" << i << " ";
+    }
+
+    file << "\n";
+
+    // For every subect
+    for(int j=0; j<size[3]; j++)
+    {
+
+      file << individuals[j] << " ";
+
+      // For every point on the current ROI
+      for(int k=0; k<m_Roi.size(); k++)
+      {
+        itk::Index<3> ix = m_Roi[k];
+        itk::Index<4> ix4;
+        ix4[0] = ix[0];
+        ix4[1] = ix[1];
+        ix4[2] = ix[2];
+        ix4[3] = j;
+
+        float f = m_Projections->GetPixel(ix4);
+
+        file << f << " ";
+
+      }
+
+      file << "\n";
+
+    }
+
+    file.close();
+
+
+
+    // Write the long format output
+    std::ofstream fileLong(m_FileNameLong.c_str());
+
+    fileLong << "ID " << "group " << "position " << "value\n";
+
+    for(int i=0; i<m_Roi.size(); i++)
+    {
+      for (int j=0; j<individuals.size(); j++)
+      {
+        itk::Index<3> ix = m_Roi[i];
+        itk::Index<4> ix4;
+        ix4[0] = ix[0];
+        ix4[1] = ix[1];
+        ix4[2] = ix[2];
+        ix4[3] = j;
+
+        float f = m_Projections->GetPixel(ix4);
+
+        fileLong << "ID" << j << " " << individuals[j] << " " << "pos" << i << " "<< f << "\n";
+      }
+    }
+
+    fileLong.close();
+
+
+  }
+
+
+
+}
+#endif
diff --git a/Modules/DiffusionImaging/Algorithms/mitkTractAnalyzer.h b/Modules/DiffusionImaging/Algorithms/mitkTractAnalyzer.h
new file mode 100644
index 0000000000..6fbc00d524
--- /dev/null
+++ b/Modules/DiffusionImaging/Algorithms/mitkTractAnalyzer.h
@@ -0,0 +1,144 @@
+/*=========================================================================
+
+  Program:   Insight Segmentation & Registration Toolkit
+  Module:    $RCSfile: itkDiffusionTensor3DReconstructionImageFilter.h,v $
+  Language:  C++
+  Date:      $Date: 2006-03-27 17:01:06 $
+  Version:   $Revision: 1.12 $
+
+  Copyright (c) Insight Software Consortium. All rights reserved.
+  See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
+
+     This software is distributed WITHOUT ANY WARRANTY; without even 
+     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
+     PURPOSE.  See the above copyright notices for more information.
+
+=========================================================================*/
+#ifndef __mitkTractAnalyzer_h_
+#define __mitkTractAnalyzer_h_
+
+#include "MitkDiffusionImagingMBIExports.h"
+#include <itkImage.h>
+#include <mitkTbssImage.h>
+
+
+namespace mitk{
+/** \class TractAnalyzer
+ */
+
+
+class MitkDiffusionImagingMBI_EXPORT TractAnalyzer
+{
+
+public:
+
+
+  TractAnalyzer();
+  ~TractAnalyzer() {};
+  
+  typedef itk::Image<char,3> CharImageType;
+  typedef itk::Image<float,3> FloatImageType;
+  typedef itk::Image<float,4> ProjectionsImageType;
+  typedef itk::VectorImage<float, 3> VectorImageType;
+  
+  
+
+  /*
+  void SetSkeleton(CharImageType::Pointer skeleton)
+  {
+    m_Skeleton = skeleton;
+  }
+
+  void SetMeanSkeleton(FloatImageType::Pointer i)
+  {
+    m_MeanSkeleton = i;
+  }*/
+
+
+  void SetTbssImage(mitk::TbssImage::Pointer tbssImg)
+  {
+    m_TbssImage = tbssImg;
+  }
+
+  void SetProjections(ProjectionsImageType::Pointer projections)
+  {
+    m_Projections = projections;
+  }
+
+  void BuildGraph(itk::Index<3> startPoint, itk::Index<3> endPoint);
+
+  std::vector< itk::Index<3> > GetPath()
+  {
+    return m_Path;
+  }
+
+  void SetFileName(std::string fname)
+  {
+    m_FileName = fname;
+  }
+
+  void SetFileNameLong(std::string fname)
+  {
+    m_FileNameLong = fname;
+  }
+
+  void SetRoi(std::vector< itk::Index<3> > roi)
+  {
+    m_Roi = roi;
+  }
+
+  CharImageType::Pointer GetRoiImage()
+  {
+    return m_RoiImg;
+  }
+
+  void SetGroups(std::vector< std::pair<std::string, int> > groups)
+  {
+    m_Groups = groups;
+  }
+
+  void MeasureRoi();
+
+  void SetInputImage(mitk::Image::Pointer inputImage)
+  {
+    m_InputImage = inputImage;
+  }
+
+
+  void SetThreshold(double threshold)
+  {
+    m_Threshold = threshold;
+  }
+
+
+protected:
+  
+
+  //CharImageType::Pointer m_Skeleton;
+  CharImageType::Pointer m_RoiImg;
+  ProjectionsImageType::Pointer m_Projections;
+  //FloatImageType::Pointer m_MeanSkeleton;
+  mitk::TbssImage::Pointer m_TbssImage;
+
+  mitk::Image::Pointer m_InputImage;
+
+  double m_Threshold;
+
+  std::vector< itk::Index<3> > m_Path;
+
+  std::string m_FileName;
+
+  std::string m_FileNameLong; // For the regression analysis 'friendly' file
+
+  std::vector< std::pair<std::string, int> > m_Groups;
+
+  std::vector< itk::Index<3> > m_Roi;
+
+private:
+
+};
+
+}
+
+#endif //__itkTractAnalyzer_h_
+