diff --git a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkPreprocessingView.cpp b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkPreprocessingView.cpp
index 38a10a8072..c17dfc7dbc 100644
--- a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkPreprocessingView.cpp
+++ b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkPreprocessingView.cpp
@@ -1,467 +1,466 @@
 /*=========================================================================
 
 Program:   Medical Imaging & Interaction Toolkit
 Module:    $RCSfile$
 Language:  C++
 Date:      $Date: 2009-05-28 17:19:30 +0200 (Do, 28 Mai 2009) $
 Version:   $Revision: 17495 $
 
 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.
 
 =========================================================================*/
 
 //#define MBILOG_ENABLE_DEBUG
 
 #include "QmitkPreprocessingView.h"
 #include "mitkDiffusionImagingConfigure.h"
 
 // qt includes
 #include <QMessageBox>
 
 // itk includes
 #include "itkTimeProbe.h"
 
 // mitk includes
 #include "mitkProgressBar.h"
 #include "mitkStatusBar.h"
 
 #include "mitkNodePredicateDataType.h"
 #include "QmitkDataStorageComboBox.h"
 #include "QmitkStdMultiWidget.h"
 
 #include "itkVectorContainer.h"
 
 #include "mitkProperties.h"
 #include "mitkVtkResliceInterpolationProperty.h"
 #include "mitkLookupTable.h"
 #include "mitkLookupTableProperty.h"
 #include "mitkTransferFunction.h"
 #include "mitkTransferFunctionProperty.h"
 #include "mitkDataNodeObject.h"
 #include "mitkOdfNormalizationMethodProperty.h"
 #include "mitkOdfScaleByProperty.h"
 
 #include "itkB0ImageExtractionImageFilter.h"
 #include "itkBrainMaskExtractionImageFilter.h"
 #include "itkCastImageFilter.h"
 
 #include "berryIStructuredSelection.h"
 #include "berryIWorkbenchWindow.h"
 #include "berryISelectionService.h"
 
 #include <QTableWidgetItem>
 #include <QTableWidget>
 
 const std::string QmitkPreprocessingView::VIEW_ID =
 "org.mitk.views.preprocessing";
 
 #define DI_INFO MITK_INFO("DiffusionImaging")
 
 
 typedef float TTensorPixelType;
 
 using namespace berry;
 
 struct PrpSelListener : ISelectionListener
 {
 
   berryObjectMacro(PrpSelListener);
 
   PrpSelListener(QmitkPreprocessingView* view)
   {
     m_View = view;
   }
 
   void DoSelectionChanged(ISelection::ConstPointer selection)
   {
     // save current selection in member variable
     m_View->m_CurrentSelection = selection.Cast<const IStructuredSelection>();
 
     // do something with the selected items
     if(m_View->m_CurrentSelection)
     {
       bool foundDwiVolume = false;
       m_View->m_DiffusionImage = NULL;
 
       // iterate selection
       for (IStructuredSelection::iterator i = m_View->m_CurrentSelection->Begin();
         i != m_View->m_CurrentSelection->End(); ++i)
       {
 
         // extract datatree node
         if (mitk::DataNodeObject::Pointer nodeObj = i->Cast<mitk::DataNodeObject>())
         {
           mitk::DataNode::Pointer node = nodeObj->GetDataNode();
 
           // only look at interesting types
           if(QString("DiffusionImage").compare(node->GetData()->GetNameOfClass())==0)
           {
             foundDwiVolume = true;
             m_View->m_DiffusionImage = dynamic_cast<mitk::DiffusionImage<DiffusionPixelType>*>(node->GetData());
           }
         }
       }
 
       m_View->m_Controls->m_ButtonBrainMask->setEnabled(foundDwiVolume);
       m_View->m_Controls->m_ButtonAverageGradients->setEnabled(foundDwiVolume);
       m_View->m_Controls->m_ButtonExtractB0->setEnabled(foundDwiVolume);
       m_View->m_Controls->m_ModifyMeasurementFrame->setEnabled(foundDwiVolume);
       m_View->m_Controls->m_MeasurementFrameTable->setEnabled(foundDwiVolume);
 
       if (foundDwiVolume)
       {
         vnl_matrix_fixed< double, 3, 3 > mf = m_View->m_DiffusionImage->GetMeasurementFrame();
         for (int r=0; r<3; r++)
           for (int c=0; c<3; c++)
           {
             QTableWidgetItem* item = m_View->m_Controls->m_MeasurementFrameTable->item(r,c);
             delete item;
             item = new QTableWidgetItem();
             item->setTextAlignment(Qt::AlignCenter | Qt::AlignVCenter);
             item->setText(QString::number(mf.get(r,c)));
             m_View->m_Controls->m_MeasurementFrameTable->setItem(r,c,item);
           }
         m_View->m_Controls->m_GradientsLabel->setText(QString::number(m_View->m_DiffusionImage->GetNumDirections()));
 
         if (m_View->m_DiffusionImage->IsMultiBval())
           m_View->m_Controls->m_BvalLabel->setText("Acquisition with multiple b-values!");
         else
           m_View->m_Controls->m_BvalLabel->setText(QString::number(m_View->m_DiffusionImage->GetB_Value()));
       }
       else
       {
         for (int r=0; r<3; r++)
           for (int c=0; c<3; c++)
           {
             QTableWidgetItem* item = m_View->m_Controls->m_MeasurementFrameTable->item(r,c);
             delete item;
             item = new QTableWidgetItem();
             m_View->m_Controls->m_MeasurementFrameTable->setItem(r,c,item);
           }
         m_View->m_Controls->m_GradientsLabel->setText("-");
         m_View->m_Controls->m_BvalLabel->setText("-");
       }
     }
   }
 
   void SelectionChanged(IWorkbenchPart::Pointer part, ISelection::ConstPointer selection)
   {
     // check, if selection comes from datamanager
     if (part)
     {
       QString partname(part->GetPartName().c_str());
       if(partname.compare("Datamanager")==0)
       {
 
         // apply selection
         DoSelectionChanged(selection);
 
       }
     }
   }
 
   QmitkPreprocessingView* m_View;
 };
 
 QmitkPreprocessingView::QmitkPreprocessingView()
 : QmitkFunctionality(),
 m_Controls(NULL),
 m_MultiWidget(NULL),
 m_DiffusionImage(NULL)
 {
 }
 
 QmitkPreprocessingView::QmitkPreprocessingView(const QmitkPreprocessingView& other)
 {
   Q_UNUSED(other)
   throw std::runtime_error("Copy constructor not implemented");
 }
 
 QmitkPreprocessingView::~QmitkPreprocessingView()
 {
   this->GetSite()->GetWorkbenchWindow()->GetSelectionService()->RemovePostSelectionListener(/*"org.mitk.views.datamanager",*/ m_SelListener);
 }
 
 void QmitkPreprocessingView::CreateQtPartControl(QWidget *parent)
 {
   if (!m_Controls)
   {
     // create GUI widgets
     m_Controls = new Ui::QmitkPreprocessingViewControls;
     m_Controls->setupUi(parent);
     this->CreateConnections();
 
     m_Controls->m_MeasurementFrameTable->horizontalHeader()->setResizeMode(QHeaderView::Stretch);
     m_Controls->m_MeasurementFrameTable->verticalHeader()->setResizeMode(QHeaderView::Stretch);
   }
 
   m_SelListener = berry::ISelectionListener::Pointer(new PrpSelListener(this));
   this->GetSite()->GetWorkbenchWindow()->GetSelectionService()->AddPostSelectionListener(/*"org.mitk.views.datamanager",*/ m_SelListener);
   berry::ISelection::ConstPointer sel(
     this->GetSite()->GetWorkbenchWindow()->GetSelectionService()->GetSelection("org.mitk.views.datamanager"));
   m_CurrentSelection = sel.Cast<const IStructuredSelection>();
   m_SelListener.Cast<PrpSelListener>()->DoSelectionChanged(sel);
 }
 
 void QmitkPreprocessingView::StdMultiWidgetAvailable (QmitkStdMultiWidget &stdMultiWidget)
 {
   m_MultiWidget = &stdMultiWidget;
 }
 
 void QmitkPreprocessingView::StdMultiWidgetNotAvailable()
 {
   m_MultiWidget = NULL;
 }
 
 void QmitkPreprocessingView::CreateConnections()
 {
   if ( m_Controls )
   {
     connect( (QObject*)(m_Controls->m_ButtonAverageGradients), SIGNAL(clicked()), this, SLOT(AverageGradients()) );
     connect( (QObject*)(m_Controls->m_ButtonExtractB0), SIGNAL(clicked()), this, SLOT(ExtractB0()) );
     connect( (QObject*)(m_Controls->m_ButtonBrainMask), SIGNAL(clicked()), this, SLOT(BrainMask()) );
     connect( (QObject*)(m_Controls->m_ModifyMeasurementFrame), SIGNAL(clicked()), this, SLOT(ApplyMesurementFrame()) );
   }
 }
 
 void QmitkPreprocessingView::Activated()
 {
   QmitkFunctionality::Activated();
 
   berry::ISelection::ConstPointer sel(
     this->GetSite()->GetWorkbenchWindow()->GetSelectionService()->GetSelection("org.mitk.views.datamanager"));
   m_CurrentSelection = sel.Cast<const IStructuredSelection>();
   m_SelListener.Cast<PrpSelListener>()->DoSelectionChanged(sel);
 }
 
 void QmitkPreprocessingView::Deactivated()
 {
   QmitkFunctionality::Deactivated();
 }
 
 void QmitkPreprocessingView::ApplyMesurementFrame()
 {
   if (m_DiffusionImage.IsNull())
     return;
   vnl_matrix_fixed< double, 3, 3 > mf;
   for (int r=0; r<3; r++)
     for (int c=0; c<3; c++)
     {
       QTableWidgetItem* item = m_Controls->m_MeasurementFrameTable->item(r,c);
       if (!item)
         return;
       mf[r][c] = item->text().toDouble();
     }
   m_DiffusionImage->SetMeasurementFrame(mf);
 }
 
 void QmitkPreprocessingView::ExtractB0()
 {
 
   if (m_CurrentSelection)
   {
     mitk::DataStorage::SetOfObjects::Pointer set =
       mitk::DataStorage::SetOfObjects::New();
 
     int at = 0;
     for (IStructuredSelection::iterator i = m_CurrentSelection->Begin();
       i != m_CurrentSelection->End();
       ++i)
     {
 
       if (mitk::DataNodeObject::Pointer nodeObj = i->Cast<mitk::DataNodeObject>())
       {
         mitk::DataNode::Pointer node = nodeObj->GetDataNode();
         if(QString("DiffusionImage").compare(node->GetData()->GetNameOfClass())==0)
         {
           set->InsertElement(at++, node);
         }
       }
     }
 
     DoExtractB0(set);
 
   }
 }
 
 void QmitkPreprocessingView::DoExtractB0
 (mitk::DataStorage::SetOfObjects::Pointer inImages)
 {
   typedef mitk::DiffusionImage<DiffusionPixelType>              DiffusionImageType;
   typedef DiffusionImageType::GradientDirectionContainerType    GradientContainerType;
 
   int nrFiles = inImages->size();
   if (!nrFiles) return;
 
   mitk::DataStorage::SetOfObjects::const_iterator itemiter( inImages->begin() );
   mitk::DataStorage::SetOfObjects::const_iterator itemiterend( inImages->end() );
 
   std::vector<mitk::DataNode::Pointer> nodes;
   while ( itemiter != itemiterend ) // for all items
   {
 
     DiffusionImageType* vols =
       static_cast<DiffusionImageType*>(
       (*itemiter)->GetData());
 
     std::string nodename;
     (*itemiter)->GetStringProperty("name", nodename);
 
     // Extract image using found index
     typedef itk::B0ImageExtractionImageFilter<short, short> FilterType;
     FilterType::Pointer filter = FilterType::New();
     filter->SetInput(vols->GetVectorImage());
     filter->SetDirections(vols->GetDirections());
     filter->Update();
 
     mitk::Image::Pointer mitkImage = mitk::Image::New();
     mitkImage->InitializeByItk( filter->GetOutput() );
     mitkImage->SetVolume( filter->GetOutput()->GetBufferPointer() );
     mitk::DataNode::Pointer node=mitk::DataNode::New();
     node->SetData( mitkImage );
     node->SetProperty( "name", mitk::StringProperty::New(nodename + "_B0"));
 
     GetDefaultDataStorage()->Add(node);
 
     ++itemiter;
   }
 
 
 }
 
 void QmitkPreprocessingView::AverageGradients()
 {
 
   if (m_CurrentSelection)
   {
     mitk::DataStorage::SetOfObjects::Pointer set =
       mitk::DataStorage::SetOfObjects::New();
 
     int at = 0;
     for (IStructuredSelection::iterator i = m_CurrentSelection->Begin();
       i != m_CurrentSelection->End();
       ++i)
     {
 
       if (mitk::DataNodeObject::Pointer nodeObj = i->Cast<mitk::DataNodeObject>())
       {
         mitk::DataNode::Pointer node = nodeObj->GetDataNode();
         if(QString("DiffusionImage").compare(node->GetData()->GetNameOfClass())==0)
         {
           set->InsertElement(at++, node);
         }
       }
     }
 
     DoAverageGradients(set);
 
   }
 
 }
 
 void QmitkPreprocessingView::DoAverageGradients
 (mitk::DataStorage::SetOfObjects::Pointer inImages)
 {
   int nrFiles = inImages->size();
   if (!nrFiles) return;
 
   mitk::DataStorage::SetOfObjects::const_iterator itemiter( inImages->begin() );
   mitk::DataStorage::SetOfObjects::const_iterator itemiterend( inImages->end() );
 
   std::vector<mitk::DataNode::Pointer> nodes;
   while ( itemiter != itemiterend ) // for all items
   {
 
     mitk::DiffusionImage<DiffusionPixelType>* vols =
       static_cast<mitk::DiffusionImage<DiffusionPixelType>*>(
       (*itemiter)->GetData());
 
     vols->AverageRedundantGradients(m_Controls->m_Blur->value());
 
     ++itemiter;
   }
 }
 
 void QmitkPreprocessingView::BrainMask()
 {
 
   if (m_CurrentSelection)
   {
     mitk::DataStorage::SetOfObjects::Pointer set =
       mitk::DataStorage::SetOfObjects::New();
 
     int at = 0;
     for (IStructuredSelection::iterator i = m_CurrentSelection->Begin();
       i != m_CurrentSelection->End();
       ++i)
     {
 
       if (mitk::DataNodeObject::Pointer nodeObj = i->Cast<mitk::DataNodeObject>())
       {
         mitk::DataNode::Pointer node = nodeObj->GetDataNode();
         if(QString("DiffusionImage").compare(node->GetData()->GetNameOfClass())==0)
         {
           set->InsertElement(at++, node);
         }
       }
     }
 
     DoBrainMask(set);
 
   }
 
 }
 
 void QmitkPreprocessingView::DoBrainMask
 (mitk::DataStorage::SetOfObjects::Pointer inImages)
 {
   int nrFiles = inImages->size();
   if (!nrFiles) return;
 
   mitk::DataStorage::SetOfObjects::const_iterator itemiter( inImages->begin() );
   mitk::DataStorage::SetOfObjects::const_iterator itemiterend( inImages->end() );
 
-  std::vector<mitk::DataNode::Pointer> nodes;
   while ( itemiter != itemiterend ) // for all items
   {
 
     mitk::DiffusionImage<DiffusionPixelType>* vols =
       static_cast<mitk::DiffusionImage<DiffusionPixelType>*>(
       (*itemiter)->GetData());
 
     std::string nodename;
     (*itemiter)->GetStringProperty("name", nodename);
 
     // Extract image using found index
     typedef itk::B0ImageExtractionImageFilter<short, short> FilterType;
     FilterType::Pointer filter = FilterType::New();
     filter->SetInput(vols->GetVectorImage());
     filter->SetDirections(vols->GetDirections());
 
     typedef itk::CastImageFilter<itk::Image<short,3>, itk::Image<unsigned short,3> > CastFilterType;
     CastFilterType::Pointer castfilter = CastFilterType::New();
     castfilter->SetInput(filter->GetOutput());
 
     typedef itk::BrainMaskExtractionImageFilter<unsigned char> MaskFilterType;
     MaskFilterType::Pointer maskfilter = MaskFilterType::New();
     maskfilter->SetInput(castfilter->GetOutput());
     maskfilter->Update();
 
     mitk::Image::Pointer mitkImage = mitk::Image::New();
     mitkImage->InitializeByItk( maskfilter->GetOutput() );
     mitkImage->SetVolume( maskfilter->GetOutput()->GetBufferPointer() );
     mitk::DataNode::Pointer node=mitk::DataNode::New();
     node->SetData( mitkImage );
     node->SetProperty( "name", mitk::StringProperty::New(nodename + "_Mask"));
 
     GetDefaultDataStorage()->Add(node);
 
     ++itemiter;
   }
 }
diff --git a/Modules/DiffusionImaging/Algorithms/itkBrainMaskExtractionImageFilter.txx b/Modules/DiffusionImaging/Algorithms/itkBrainMaskExtractionImageFilter.txx
index 3d1f1da9fc..80a009fe07 100644
--- a/Modules/DiffusionImaging/Algorithms/itkBrainMaskExtractionImageFilter.txx
+++ b/Modules/DiffusionImaging/Algorithms/itkBrainMaskExtractionImageFilter.txx
@@ -1,436 +1,437 @@
 /*=========================================================================
 
 Program:   Insight Segmentation & Registration Toolkit
 Language:  C++
 Date:      $Date: 2006-07-19 15:11:41 $
 Version:   $Revision: 1.11 $
 
 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 
+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 __itkBrainMaskExtractionImageFilter_txx
 #define __itkBrainMaskExtractionImageFilter_txx
 
 #include "itkBrainMaskExtractionImageFilter.h"
 
 #include <itkBinaryThresholdImageFilter.h>
 #include <itkBinaryErodeImageFilter.h>
 #include <itkBinaryDilateImageFilter.h>
 #include <itkVotingBinaryHoleFillingImageFilter.h>
 #include <itkVotingBinaryIterativeHoleFillingImageFilter.h>
 #include <itkBinaryBallStructuringElement.h>
 #include <itkBinaryCrossStructuringElement.h>
 #include <itkAndImageFilter.h>
-#include <itkRecursiveGaussianImageFilter.h> 
+#include <itkRecursiveGaussianImageFilter.h>
 #include <itkMaskImageFilter.h>
 
 namespace itk {
 
   template< class TOutputImagePixelType >
     BrainMaskExtractionImageFilter< TOutputImagePixelType >
     ::BrainMaskExtractionImageFilter()
   {
     // At least 1 inputs is necessary for a vector image.
     // For images added one at a time we need at least six
-    this->SetNumberOfRequiredInputs( 1 ); 
+    this->SetNumberOfRequiredInputs( 1 );
   }
 
   template< class TOutputImagePixelType >
     void BrainMaskExtractionImageFilter< TOutputImagePixelType >
     ::GenerateData()
   {
 
     // gaussian filter
     typename itk::RecursiveGaussianImageFilter<InputImageType,InputImageType>::Pointer gaussian
       = itk::RecursiveGaussianImageFilter<InputImageType,InputImageType>::New();
     gaussian->SetInput( this->GetInput(0) );
     gaussian->SetSigma( 1.0 );
 
     try
     {
       gaussian->Update();
     }
     catch( itk::ExceptionObject &e)
     {
       std::cerr << e;
       return;
     }
 
     // threshold the image
-    typename itk::BinaryThresholdImageFilter<InputImageType,OutputImageType>::Pointer threshold = 
+    typename itk::BinaryThresholdImageFilter<InputImageType,OutputImageType>::Pointer threshold =
       itk::BinaryThresholdImageFilter<InputImageType,OutputImageType>::New();
 
     threshold->SetInput( gaussian->GetOutput() );
 
     int seuil = static_cast<int>( ComputeHistogram( gaussian->GetOutput() ) );
     threshold->SetLowerThreshold( seuil );
 
     std::cout << "Thresholding..." << std::flush;
     try
     {
       threshold->Update();
     }
     catch( itk::ExceptionObject &e)
     {
       std::cerr << e;
       return;
     }
     std::cout << "Done." << std::endl;
 
 #ifdef DEBUG_ME
     {
       WriterType::Pointer writer = WriterType::New();
       writer->SetInput( threshold->GetOutput() );
       writer->SetFileName( "AfterThreshold.hdr" );
       writer->Update();
     }
 #endif
 
     // erode to remove background noise
     typedef itk::BinaryBallStructuringElement<int, 3> StructuralElementType;
     StructuralElementType ball;
 
-    typename itk::BinaryErodeImageFilter<OutputImageType,OutputImageType,StructuralElementType>::Pointer erode = 
+    typename itk::BinaryErodeImageFilter<OutputImageType,OutputImageType,StructuralElementType>::Pointer erode =
       itk::BinaryErodeImageFilter<OutputImageType,OutputImageType,StructuralElementType>::New();
 
     ball.SetRadius( 3 );
 
     erode->SetInput( threshold->GetOutput() );
     erode->SetKernel( ball );
 
     std::cout << "Eroding..." << std::flush;
     try
     {
       erode->Update();
     }
     catch( itk::ExceptionObject &e)
     {
       std::cerr << e;
       return;
     }
     std::cout << "Done." << std::endl;
 
 #ifdef DEBUG_ME
     {
       WriterType::Pointer writer = WriterType::New();
       writer->SetInput( erode->GetOutput() );
       writer->SetFileName( "AfterErode.hdr" );
       writer->Update();
     }
 #endif
 
 
     typedef BinaryCrossStructuringElement<int, 3> CrossType;
 
     typedef BinaryDilateImageFilter<OutputImageType,OutputImageType,CrossType> DilateFilterType;
     typedef AndImageFilter<OutputImageType,OutputImageType,OutputImageType> AndFilterType;
 
     typename OutputImageType::Pointer M0 = threshold->GetOutput();
     typename OutputImageType::Pointer Mn = OutputImageType::New();
     Mn->SetRegions( M0->GetLargestPossibleRegion() );
     Mn->SetSpacing( M0->GetSpacing() );
     Mn->SetOrigin( M0->GetOrigin() );
+    Mn->SetDirection( M0->GetDirection() );
     Mn->Allocate();
 
     typename OutputImageType::Pointer Mnplus1 = erode->GetOutput();
 
 
     CrossType cross;
     cross.SetRadius( 1 );
     //unsigned long rad[3]={3,3,3};
 
     //ball2.SetRadius( rad );
 
     std::cout << "Conditional reconstruction..." << std::flush;
     int iter = 0;
     do
     {
       std::cout << "Iteration: " << iter++ << std::endl;
       CopyImage( Mn, Mnplus1);
 
       typename DilateFilterType::Pointer dilater = DilateFilterType::New();
       dilater->SetInput( Mn );
       dilater->SetKernel( cross );
 
       try
       {
         dilater->Update();
       }
       catch( itk::ExceptionObject &e)
       {
         std::cerr << e;
         return;
       }
 
       typename AndFilterType::Pointer andfilter = AndFilterType::New();
       andfilter->SetInput(0, M0);
       andfilter->SetInput(1, dilater->GetOutput() );
 
       try
       {
         andfilter->Update();
       }
       catch( itk::ExceptionObject &e)
       {
         std::cerr << e;
         return;
       }
 
       Mnplus1 = andfilter->GetOutput();
 
       /*
       #ifdef DEBUG_ME
       {
       WriterType::Pointer writer = WriterType::New();
       writer->SetInput( andfilter->GetOutput() );
       char filename[512];
       sprintf( filename, "CondReconstruction_iter_%d.hdr", iter);
       writer->SetFileName( filename );
       writer->Update();
       }
       #endif*/
 
     } while( !CompareImages( Mn, Mnplus1) );
 
     std::cout << "Done." << std::endl;
 
 #ifdef DEBUG_ME
     {
       WriterType::Pointer writer = WriterType::New();
       writer->SetInput( Mn );
       writer->SetFileName( "AfterCondReconstruction.hdr" );
       writer->Update();
     }
 #endif
     // now fill the holes
 
     typename itk::VotingBinaryIterativeHoleFillingImageFilter< OutputImageType >::Pointer filler =
       itk::VotingBinaryIterativeHoleFillingImageFilter< OutputImageType >::New();
     filler->SetInput( Mn );
     filler->SetMaximumNumberOfIterations (1000);
 
     std::cout << "Filling the holes..." << std::flush;
     try
     {
       filler->Update();
     }
     catch( itk::ExceptionObject &e)
     {
       std::cerr << e;
       return;
     }
     std::cout << "Done." << std::endl;
 
-    typename OutputImageType::Pointer outputImage = 
+    typename OutputImageType::Pointer outputImage =
       static_cast< OutputImageType * >(this->ProcessObject::GetOutput(0));
 
     outputImage->SetSpacing( filler->GetOutput()->GetSpacing() );   // Set the image spacing
     outputImage->SetOrigin( filler->GetOutput()->GetOrigin() );     // Set the image origin
-    outputImage->SetDirection( filler->GetOutput()->GetDirection() );  // Set the image direction
     outputImage->SetLargestPossibleRegion( filler->GetOutput()->GetLargestPossibleRegion());
     outputImage->SetBufferedRegion( filler->GetOutput()->GetLargestPossibleRegion() );
+    outputImage->SetDirection( filler->GetOutput()->GetDirection() );
     outputImage->Allocate();
 
     itk::ImageRegionIterator<OutputImageType> itIn( filler->GetOutput(), filler->GetOutput()->GetLargestPossibleRegion() );
     itk::ImageRegionIterator<OutputImageType> itOut( outputImage, outputImage->GetLargestPossibleRegion() );
 
     while( !itIn.IsAtEnd() )
     {
         itOut.Set(itIn.Get());
         ++itIn;
         ++itOut;
     }
 
   }
 
 
   template< class TOutputImagePixelType >
     void BrainMaskExtractionImageFilter< TOutputImagePixelType >
     ::CopyImage( typename OutputImageType::Pointer target, typename OutputImageType::Pointer source)
   {
     itk::ImageRegionConstIterator<OutputImageType> itIn( source, source->GetLargestPossibleRegion() );
     itk::ImageRegionIterator<OutputImageType> itOut( target, target->GetLargestPossibleRegion() );
 
     while( !itOut.IsAtEnd() )
     {
       itOut.Set( itIn.Get() );
       ++itIn;
       ++itOut;
     }
   }
 
 
   template< class TOutputImagePixelType >
     bool BrainMaskExtractionImageFilter< TOutputImagePixelType >
     ::CompareImages( typename OutputImageType::Pointer im1, typename OutputImageType::Pointer im2)
   {
     itk::ImageRegionConstIterator<OutputImageType> itIn( im1, im1->GetLargestPossibleRegion() );
     itk::ImageRegionConstIterator<OutputImageType> itOut( im2, im2->GetLargestPossibleRegion() );
 
     while( !itOut.IsAtEnd() )
     {
       if( itOut.Value() != itIn.Value() )
       {
         return false;
       }
       ++itOut;
       ++itIn;
     }
 
     return true;
   }
 
 
   template< class TOutputImagePixelType >
     int BrainMaskExtractionImageFilter< TOutputImagePixelType >
     ::ComputeHistogram( typename InputImageType::Pointer image)
   {
     // IMPORTANT: IMAGE MUST BE UNSIGNED SHORT
     int N=65535;
     int* histogram = new int[N];
     for( int i=0; i<N; i++)
     {
       histogram[i] = 0;
     }
 
     itk::ImageRegionConstIterator<InputImageType> itIn( image, image->GetLargestPossibleRegion() );
 
     long totVoxels = 0;
 
     int max = -1;
     int min = 9999999;
     while( !itIn.IsAtEnd() )
     {
       histogram[ (int)(itIn.Value()) ]++;
 
       if( itIn.Value()>max )
       {
         max = itIn.Value();
       }
 
       if( itIn.Value()<min )
       {
         min = itIn.Value();
       }
       ++itIn;
       ++totVoxels;
     }
 
 
     //int EPS = 1;
     int seuil = 0;
     //int newseuil = (max + min)/2;
 
     N = max;
 
     double V = 0.0;
     int S = 0;
 
     double mean = 0.0;
     for( int i=min; i<=max; i++)
     {
       mean += (double)(i) * (double)histogram[i];
     }
     mean /= (double)totVoxels;
 
     //std::cout << "Min/Max/Mean: " << min << "/" << max << "/" << mean << std::endl;
 
     //while( abs(newseuil - seuil)>EPS )
     for( seuil = min; seuil<=max; seuil++)
     {
       //seuil = newseuil;
 
       // compute the classes:
       double mean0 = 0.0;
       double mean1 = 0.0;
       //double std0 = 0.0;
       //double std1 = 0.0;
 
       double num0 = 0.0;
       double num1 = 0.0;
       for( int i=min; i<seuil; i++)
       {
         //mean0 += (double)(histogram[i])/(double)(totVoxels) * (double)(i);
         mean0 += (double)histogram[i] * (double)i;
         num0  += (double)histogram[i];
       }
 
       for( int i=seuil; i<max; i++)
       {
         //mean1 += (double)(histogram[i])/(double)(totVoxels) * (double)(i);
         mean1 += (double)histogram[i] * (double)i;
         num1  += (double)histogram[i];
       }
 
       if( num0 )
       {
         mean0/=num0;
       }
       if( num1 )
       {
         mean1/=num1;
       }
 
       /*
       for( int i=0; i<seuil; i++)
       {
       std0 += (double)histogram[i]/(double)totVoxels* ((double)i - mean0)*((double)i - mean0);
       //std0 += (double)histogram[i]/num0* ((double)i - mean0)*((double)i - mean0);
       }
 
       for( int i=seuil; i<65536; i++)
       {
       std1 += (double)histogram[i]/(double)totVoxels* ((double)i - mean1)* ((double)i - mean1);
       //std1 += (double)histogram[i]/num1* ((double)i - mean1)* ((double)i - mean1);
       }
       */
 
 
       //std0 = sqrt(std0);
       //std1 = sqrt(std1);
 
       //std::cout << "Classe 0: " << mean0 << " " << std0 << std::endl;
       //std::cout << "Classe 1: " << mean1 << " " << std1 << std::endl;
 
       //std::cout << "Classe 0: " << mean0 << std::endl;
       //std::cout << "Classe 1: " << mean1 << std::endl;
 
       //double alpha = std0*std1/(std0+std1);
 
       //std::cout << "alpha: " << alpha << " " << std0 << " " << std1 << std::endl;
-      //newseuil = (int) (alpha*(mean0/std0 + mean1/std1));		
+      //newseuil = (int) (alpha*(mean0/std0 + mean1/std1));
 
       //std::cout << "New threshold: " << newseuil << std::endl;
 
       /*
       if( newseuil > 65535 || newseuil<0)
       {
       std::cerr << "Error: threshold is too low or high, exiting" << std::endl;
       return -1;
       }*/
 
       double Vm = num0 * (mean0 - mean)*(mean0 - mean) + num1*(mean1 - mean)*(mean1 - mean);
       if( Vm > V )
       {
         V = Vm;
         S = seuil;
         //std::cout << "New seuil: " << S << std::endl;
         //getchar();
       }
 
     }
 
     delete [] histogram;
 
     std::cout << "Seuil: " << S << std::endl;
 
     return S;
   }
 }
 
 #endif // __itkBrainMaskExtractionImageFilter_txx