diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorReconstructionWithEigenvalueCorrectionFilter.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorReconstructionWithEigenvalueCorrectionFilter.h
index b9259e967c..fcaf9ddbee 100644
--- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorReconstructionWithEigenvalueCorrectionFilter.h
+++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorReconstructionWithEigenvalueCorrectionFilter.h
@@ -1,209 +1,209 @@
 /*===================================================================
 
 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 _itk_TensorReconstructionWithEigenvalueCorrectionFilter_h_
 #define _itk_TensorReconstructionWithEigenvalueCorrectionFilter_h_
 
 #include "itkImageToImageFilter.h"
 #include <itkDiffusionTensor3D.h>
 #include <itkVectorImage.h>
 #include <vnl/algo/vnl_symmetric_eigensystem.h>
 
 #include <math.h>
 #include <mitkDiffusionImage.h>
 
 
 
 
 namespace itk
 {
 
   template <class TDiffusionPixelType, class TTensorPixelType>
   class TensorReconstructionWithEigenvalueCorrectionFilter
     : public ImageToImageFilter< itk::Image< TDiffusionPixelType, 3 >, itk::Image<itk::DiffusionTensor3D<TTensorPixelType>,3> >
   {
 
   public:
 
 
     typedef itk::VectorImage<short, 3>  ImageType;
 
 
 
     typedef TensorReconstructionWithEigenvalueCorrectionFilter                  Self;
     typedef SmartPointer<Self>                          Pointer;
     typedef SmartPointer<const Self>                    ConstPointer;
 
 
 
     typedef ImageToImageFilter< Image< TDiffusionPixelType, 3>,
       Image< DiffusionTensor3D< TTensorPixelType >, 3 > >
       Superclass;
 
     /** Method for creation through the object factory. */
     itkNewMacro(Self)
 
     /** Runtime information support. */
     itkTypeMacro(TensorReconstructionWithEigenvalueCorrectionFilter, ImageToImageFilter)
 
     typedef TDiffusionPixelType                       ReferencePixelType;
     typedef TDiffusionPixelType                       GradientPixelType;
     typedef DiffusionTensor3D< TTensorPixelType >     TensorPixelType;
 
 
 
     typedef Image< TensorPixelType, 3 >              TensorImageType;
     typedef TensorImageType                          OutputImageType;
     typedef typename Superclass::OutputImageRegionType
       OutputImageRegionType;
 
 
     /** An alternative typedef defining one (of the many) gradient images.
     * It will be assumed that the vectorImage has the same dimension as the
     * Reference image and a vector length parameter of \c n (number of
     * gradient directions) */
     typedef VectorImage< GradientPixelType, 3 >      GradientImagesType;
     typedef typename GradientImagesType::PixelType         GradientVectorType;
 
 
 
     /** Holds each magnetic field gradient used to acquire one DWImage */
     typedef vnl_vector_fixed< double, 3 >            GradientDirectionType;
 
     /** Container to hold gradient directions of the 'n' DW measurements */
     typedef VectorContainer< unsigned int,
       GradientDirectionType >                  GradientDirectionContainerType;
 
 
 
     /** Another set method to add a gradient directions and its corresponding
     * image. The image here is a VectorImage. The user is expected to pass the
     * gradient directions in a container. The ith element of the container
     * corresponds to the gradient direction of the ith component image the
     * VectorImage.  For the baseline image, a vector of all zeros
     * should be set. */
     void SetGradientImage( GradientDirectionContainerType *,
       const GradientImagesType *image);
 
 
 
     itkSetMacro( BValue, TTensorPixelType)
     itkSetMacro( B0Threshold, float)
 
     itkGetMacro(PseudoInverse, vnl_matrix<double>)
     itkGetMacro(H, vnl_matrix<double>)
     itkGetMacro(BVec, vnl_vector<double>)
     itkGetMacro(B0Mask, vnl_vector<short>)
 
 
     ImageType::Pointer GetCorrectedDiffusionVolumes()
     {
       return m_CorrectedDiffusionVolumes;
     }
 
     itk::Image<short, 3>::Pointer GetMask()
     {
       return m_MaskImage;
     }
 
 
   protected:
 
     TensorReconstructionWithEigenvalueCorrectionFilter();
     ~TensorReconstructionWithEigenvalueCorrectionFilter() {}
 
 
     void GenerateData();
 
 
     typedef enum
     {
       GradientIsInASingleImage = 1,
       GradientIsInManyImages,
       Else
     } GradientImageTypeEnumeration;
 
 
 
   private:
 
     double CheckNeighbours(int x, int y, int z,int f, itk::Size<3> size, itk::Image<short, 3> ::Pointer mask,itk::VectorImage<short, 3>::Pointer corrected_diffusion_temp);
 
     void CalculateAttenuation(vnl_vector<double> org_data, vnl_vector<double> &atten,int nof,int numberb0);
 
     //void CalculateTensor(vnl_vector<double> atten, vnl_vector<double> &tensor,int nof,int numberb0);
 
-    void CorrectDiffusionImage(int nof,int numberb0,itk::Size<3> size,itk::VectorImage<short, 3>::Pointer corrected_diffusion_temp,itk::VectorImage<short, 3>::Pointer corrected_diffusion,itk::Image<short, 3>::Pointer mask,vnl_vector< double> pixel_max,vnl_vector< double> pixel_min);
+    void CorrectDiffusionImage(int nof,int numberb0,itk::Size<3> size,itk::VectorImage<short, 3>::Pointer corrected_diffusion,itk::Image<short, 3>::Pointer mask,vnl_vector< double> pixel_max,vnl_vector< double> pixel_min);
 
     void GenerateTensorImage(int nof,int numberb0,itk::Size<3> size,itk::VectorImage<short, 3>::Pointer corrected_diffusion,itk::Image<short, 3>::Pointer mask,double what_mask,itk::Image< itk::DiffusionTensor3D<float>, 3 >::Pointer tensorImg );
 
     void DeepCopyTensorImage(itk::Image< itk::DiffusionTensor3D<float>, 3 >::Pointer tensorImg, itk::Image< itk::DiffusionTensor3D<float>, 3 >::Pointer temp_tensorImg);
 
     void DeepCopyDiffusionImage(itk::VectorImage<short, 3>::Pointer corrected_diffusion, itk::VectorImage<short, 3>::Pointer corrected_diffusion_temp,int nof);
 
 
     void TurnMask( itk::Size<3> size, itk::Image<short, 3>::Pointer mask, double previous_mask, double set_mask);
 
     double CheckNegatives ( itk::Size<3> size, itk::Image<short, 3>::Pointer mask, itk::Image< itk::DiffusionTensor3D<float>, 3 >::Pointer tensorImg );
 
 
     /** Gradient image was specified in a single image or in multiple images */
     GradientImageTypeEnumeration                      m_GradientImageTypeEnumeration;
 
     /** Number of gradient measurements */
     unsigned int                                      m_NumberOfGradientDirections;
 
     /** container to hold gradient directions */
     GradientDirectionContainerType::Pointer           m_GradientDirectionContainer;
 
 
     /** b-value */
     TTensorPixelType                                  m_BValue;
 
     /** Number of baseline images */
     unsigned int                                      m_NumberOfBaselineImages;
 
     ImageType::Pointer                                m_CorrectedDiffusionVolumes;
 
     float                                             m_B0Threshold;
 
 
 
     itk::Image<short, 3>::Pointer m_MaskImage;
     vnl_matrix<double> m_PseudoInverse;
     vnl_matrix<double> m_H;
     vnl_vector<double> m_BVec;
 
     /** m_B0Mask indicates whether a volume identified by an index is B0-weighted or not */
     vnl_vector<short> m_B0Mask;
 
     typename GradientImagesType::Pointer m_GradientImagePointer;
 
   };
 
 
 
 
 
 
 } // end of namespace
 
 
 #ifndef ITK_MANUAL_INSTANTIATION
 #include "itkTensorReconstructionWithEigenvalueCorrectionFilter.txx"
 #endif
 
 
 #endif
diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorReconstructionWithEigenvalueCorrectionFilter.txx b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorReconstructionWithEigenvalueCorrectionFilter.txx
index 45f7c987fe..e4fea36144 100644
--- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorReconstructionWithEigenvalueCorrectionFilter.txx
+++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorReconstructionWithEigenvalueCorrectionFilter.txx
@@ -1,965 +1,990 @@
 /*===================================================================
 
 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 _itk_TensorReconstructionWithEigenvalueCorrectionFilter_txx_
 #define _itk_TensorReconstructionWithEigenvalueCorrectioFiltern_txx_
 #endif
 #include "itkImageRegionConstIterator.h"
 #include <math.h>
 #include <mitkImageWriter.h>
 #include "itkImageFileWriter.h"
 #include "itkImage.h"
 #include "itkImageRegionIterator.h"
 
 
 
 namespace itk
 {
 
 
   template <class TDiffusionPixelType, class TTensorPixelType>
   TensorReconstructionWithEigenvalueCorrectionFilter<TDiffusionPixelType, TTensorPixelType>
     ::TensorReconstructionWithEigenvalueCorrectionFilter()
   {
     m_B0Threshold = 50.0;
   }
 
   template <class TDiffusionPixelType, class TTensorPixelType>
   void
   TensorReconstructionWithEigenvalueCorrectionFilter<TDiffusionPixelType, TTensorPixelType>
   ::GenerateData ()
   {
 
     // GradientImagesType input
     m_GradientImagePointer = static_cast< GradientImagesType * >(
       this->ProcessObject::GetInput(0) );
 
     typename GradientImagesType::SizeType size = m_GradientImagePointer->GetLargestPossibleRegion().GetSize();
 
     // number of volumes
     int nof = m_GradientDirectionContainer->Size();
 
     // determine the number of b-zero values
     int numberb0=0;
     for(int i=0; i<nof; i++)
     {
       vnl_vector_fixed <double, 3 > vec = m_GradientDirectionContainer->ElementAt(i);
 
       // due to roundings, the values are not always exactly zero
       if(vec[0]<0.0001 && vec[1]<0.0001 && vec[2]<0.0001 && vec[0]>-0.0001&& vec[1]>-0.0001 && vec[2]>-0.0001)
       {
         numberb0++;
       }
     }
 
 
-    // Matrix to store al diffusion encoding gradients
+    // Matrix to store all diffusion encoding gradients
     vnl_matrix<double> directions(nof-numberb0,3);
 
     m_B0Mask.set_size(nof);
 
 
     int cnt=0;
     for(int i=0; i<nof; i++)
     {
       vnl_vector_fixed <double, 3 > vec = m_GradientDirectionContainer->ElementAt(i);
 
       if(vec[0]<0.0001 && vec[1]<0.0001 && vec[2]<0.0001 && vec[0]>-0.001&& vec[1]>-0.001 && vec[2]>-0.001)
       {
         // the diffusion encoding gradient is approximately zero, wo we are dealing with a non-diffusion weighted volume
         m_B0Mask[i]=1;
       }
       else
       {
         // dealing with a diffusion weighted volume
         m_B0Mask[i]=0;
 
         // set the diffusion encoding gradient to the directions matrix
         directions[cnt][0] = vec[0];
         directions[cnt][1] = vec[1];
         directions[cnt][2] = vec[2];
 
         cnt++;
       }
     }
 
+    // looking for maximal norm among gradients.
+    // The norm is calculated with use of spectral radius theorem- based on determination of eigenvalue.
+
     vnl_matrix<double> dirsTimesDirsTrans = directions*directions.transpose();
 
     vnl_vector< double> diagonal(nof-numberb0);
     vnl_vector< double> b_vec(nof-numberb0);
     vnl_vector< double> temporary(3);
 
+
     for (int i=0;i<nof-numberb0;i++)
     {
       diagonal[i]=dirsTimesDirsTrans[i][i];
     }
 
     double max_diagonal = diagonal.max_value();
 
 
+    // normalization: free-water method assumes that directions matrix norm 1 is equal to b=1000
+
     directions=directions*sqrt(m_BValue/1000.0)*(1.0/sqrt(max_diagonal));
 
+    //calculation of norms and stoaring them in vector
 
     dirsTimesDirsTrans = directions*directions.transpose();
 
     for (int i=0;i<nof-numberb0;i++)
     {
       b_vec[i]= dirsTimesDirsTrans[i][i];
     }
 
+    // Calculation of so called design matrix that is used to obtain expected signal.
+
     vnl_matrix<double> H(nof-numberb0, 6);
     vnl_matrix<double> H_org(nof-numberb0, 6);
     vnl_vector<double> pre_tensor(9);
 
-    int etbt[6] = { 0, 4, 8, 1, 5, 2 };
+    //H is matrix that containes covariances for directions. It is stored twice because its original value is needed later
+    // while H is changed
+
+    int etbt[6] = { 0, 4, 8, 1, 5, 2 };// tensor order
 
     for (int i = 0; i < nof-numberb0; i++)
     {
       for (int j = 0; j < 3; j++)
       {
         temporary[j] = -directions[i][j];
       }
       for (int j = 0; j < 3; j++)
       {
        for (int k = 0; k < 3; k++)
        {
          pre_tensor[k + 3 * j] = temporary[k] * directions[i][j];
        }
       }
       for (int j = 0; j < 6; j++)
       {
         H[i][j] = pre_tensor[etbt[j]];
       }
       for (int j = 0; j < 3; j++)
       {
         H[i][3 + j] *= 2.0;
       }
     }
 
     H_org=H;
 
+    // calculation of inverse matrix by means of pseudoinverse
+
     vnl_matrix<double> inputtopseudoinverse=H.transpose()*H;
     vnl_symmetric_eigensystem<double> eig( inputtopseudoinverse);
     vnl_matrix<double> pseudoInverse = eig.pinverse()*H.transpose();
 
 
 
-
-
-
-
     ImageType::Pointer corrected_diffusion_temp = ImageType::New();
 
     typedef itk::VariableLengthVector<short> VariableVectorType;
     VariableVectorType variableLengthVector;
     variableLengthVector.SetSize(nof);
 
 
     typedef itk::VariableLengthVector<short> VariableVectorType;
     VariableVectorType corrected_single;
     corrected_single.SetSize(nof-1);
 
 
     typedef itk::Image<short, 3> MaskImageType;
     MaskImageType::Pointer mask = MaskImageType::New();
     mask->SetRegions(m_GradientImagePointer->GetLargestPossibleRegion().GetSize());
     mask->SetSpacing(m_GradientImagePointer->GetSpacing());
     mask->SetOrigin(m_GradientImagePointer->GetOrigin());
     mask->Allocate();
 
+    // Image thresholding: For every voxel mean B0 image is calculated and then voxels of mean B0 less than assumed
+    //treshold are excluded from the dataset with use of defined mask image. 1 in mask voxel has mean B0 > assumed treshold.
 
     int mask_cnt=0;
     for(int x=0;x<size[0];x++)
     {
       for(int y=0;y<size[1];y++)
       {
         for(int z=0;z<size[2];z++)
         {
           double mean_b=0.0;
 
           itk::Index<3> ix = {{x,y,z}};
 
 
           GradientVectorType pixel = m_GradientImagePointer->GetPixel(ix);
           for (int i=0;i<nof;i++)
           {
             if(m_B0Mask[i]==1)
             {
               mean_b = mean_b + pixel[i];
             }
           }
           mean_b=mean_b/numberb0;
           if(mean_b > m_B0Threshold)
           {
             mask->SetPixel(ix, 1);
             mask_cnt++;
           }
           else
           {
             mask->SetPixel(ix, 0);
           }
 
         }
       }
     }
 
 
-    //create a copy of the original image
+    //create a copy of the original image- it is then used in pre-processing methods
+
     m_CorrectedDiffusionVolumes = ImageType::New();
     m_CorrectedDiffusionVolumes->SetRegions(size);
     m_CorrectedDiffusionVolumes->SetSpacing(m_GradientImagePointer->GetSpacing());
     m_CorrectedDiffusionVolumes->SetOrigin(m_GradientImagePointer->GetOrigin());
     m_CorrectedDiffusionVolumes->SetVectorLength(nof);
     m_CorrectedDiffusionVolumes->Allocate();
 
 
     for ( int x=0;x<size[0];x++)
     {
       for ( int y=0;y<size[1];y++)
       {
         for ( int z=0;z<size[2];z++)
         {
           itk::Index<3> ix = {{x,y,z}};
 
           GradientVectorType p = m_GradientImagePointer->GetPixel(ix);
 
           m_CorrectedDiffusionVolumes->SetPixel(ix,p);
 
         }
       }
     }
 
 
 
-   //Sometimes the gradient voxels may contain negative values ( even if B0 voxel is > = 50 ). This must be corrected
+   //Sometimes the gradient voxels may contain negative values ( even if B0 voxel is > = 50 ). This must be corrected by smoothing DWI
+   //Smoothing is done by aproximation of negative voxel value by its correct ( positive) 27-th neighborhood.
 
 
     double mask_val=0.0;
     vnl_vector<double> org_vec(nof-numberb0);
 
     int counter_corrected =0;
 
     for ( int x=0;x<size[0];x++)
     {
       for ( int y=0;y<size[1];y++)
       {
         for ( int z=0;z<size[2];z++)
         {
           itk::Index<3> ix = {{x,y,z}};
 
           mask_val = mask->GetPixel(ix);
           GradientVectorType pixel2 = m_CorrectedDiffusionVolumes->GetPixel(ix);
 
           for (int i=0;i<nof;i++)
           {
             org_vec[i]=pixel2[i];
           }
 
           if(mask_val >0) // we are dooing this only if the voxels are in the mask
           {
 
           for( int f=0;f<nof;f++)
           {
             if(org_vec[f] <= 0) // if in one of the gradients it is 0 or les than 0
             {
               org_vec[f] = CheckNeighbours(x,y,z,f,size,mask,m_CorrectedDiffusionVolumes);
               counter_corrected++;
 
               // If method returned 0 it means that there is no valid ( correct) neighborhood. It is done outside the function because
               //later the same method is called in different concept when this is not important
             }
           }
 
           for (int i=0;i<nof;i++)
           {
             pixel2[i]=org_vec[i];
           }
 
           m_CorrectedDiffusionVolumes->SetPixel(ix, pixel2);
 
 
           }
         }
       }
     }
 
-
+    //Declaration of tensor image that is used for estimation of free water map
 
     typedef itk::Image< itk::DiffusionTensor3D<float>, 3 > TensorImageType;
     TensorImageType::Pointer tensorImg = TensorImageType::New();
     tensorImg->SetRegions(m_GradientImagePointer->GetLargestPossibleRegion().GetSize());
     tensorImg->SetSpacing(m_GradientImagePointer->GetSpacing());
     tensorImg->SetOrigin(m_GradientImagePointer->GetOrigin());
     tensorImg->Allocate();
 
 
     TensorImageType::Pointer temp_tensorImg = TensorImageType::New();
 
+    //Deep copy into memmory tensor image. It is done because temporary tensor image is needed for the pre-processing methods.
 
     DeepCopyTensorImage(tensorImg,temp_tensorImg);
 
-    vnl_vector< double> pixel_max(nof-numberb0);
-    vnl_vector< double> pixel_min(nof-numberb0);
+    //Declaration of vectors that contains values of toohigh and toolow atenuation for each gradient. Attenuation is only calculated for
+    //non B0 images so nof-numberb0.
+
+    vnl_vector< double> pixel_max(nof-numberb0);// to high attenuation
+    vnl_vector< double> pixel_min(nof-numberb0);//to low attenuation
+
+    // to high and to low attenuation is calculated with use of highest allowed =5 and lowest allowed =0.01 diffusion coefficient
 
     for (int i=0;i<nof-numberb0;i++)
     {
-      pixel_max[i]=exp(-b_vec[i]*0.01);// here some values for low and high diffusivity were pre define as 0.01 for low and 5 for high
+      pixel_max[i]=exp(-b_vec[i]*0.01);
       pixel_min[i]= exp(-b_vec[i]*5);
     }
 
     m_PseudoInverse = pseudoInverse;
     m_H = H_org;
     m_BVec=b_vec;
 
-
+    // in pre-processing we are dealing with 3 types of voxels: voxels excluded = 0 in mask, voxels correct =1 and voxels under correction
+    //= 2. During pre processing most of voxels should be switched from 2 to 1.
+    //To smoothly edal with the mask the Turn mask method is defined that takes previous value of mask and new value of mask.
 
     double what_mask=0;
+
+    // what amsk is a variable declared to simplify tensor calculaton. Tensors are obtained only for voxels that have a mask value bigger
+    //than what mask. Sometimes it is 1 sometimes 2.
+
+    // initialization of required variables
     double set_mask=2.0;
     double previous_mask=0;
 
     double old_number_negative_eigs=0;
     double new_number_negative_eigs=0;
 
-
     bool  stil_correcting = true;
 
-    TurnMask(size,mask,previous_mask,set_mask);// simply defining all possible tensors as with negative eigenvalues
+    TurnMask(size,mask,previous_mask,set_mask);
+    // simply defining all possible tensors as with negative eigenvalues
 
     what_mask=1.0;// we want to calculate all the tensors
 
 
+    //Pre-processing is done in iterations as long as next iteration does not introduce higher number of bad voxels.
+    //In such a case smoothed DWI should be last one that has smaller or equal number of bad voxels in comparison to the
+    //previous iteration. To obtain this temporary DWI image must be stored in memory.
+
     DeepCopyDiffusionImage(m_CorrectedDiffusionVolumes,corrected_diffusion_temp,nof);
 
-    GenerateTensorImage(nof,numberb0,size,corrected_diffusion_temp,mask,what_mask,tensorImg);// generating of the tensors
+    // generating of initial tensor image
+    GenerateTensorImage(nof,numberb0,size,corrected_diffusion_temp,mask,what_mask,tensorImg);
+
 
+    // checking how many tensors has problems, this is working only for mask =2
+    old_number_negative_eigs = CheckNegatives (size,mask,tensorImg);
 
-    old_number_negative_eigs = CheckNegatives (size,mask,tensorImg);// checking how many tensors has problems, this is working only for mask =2
 
-    std::cout << "Number of negative eigenvalues: " << old_number_negative_eigs << std::endl;// info for Thomas: Debug stuff - to be removed
+    //info for the user printed in the consol-debug information to be removed in the future
+    std::cout << "Number of negative eigenvalues: " << old_number_negative_eigs << std::endl;
 
-    CorrectDiffusionImage(nof,numberb0,size,corrected_diffusion_temp,corrected_diffusion_temp,mask,pixel_max,pixel_min);
 
+    //Smoothing DWI - method is described when it is defined
+    CorrectDiffusionImage(nof,numberb0,size,corrected_diffusion_temp,mask,pixel_max,pixel_min);
 
 
+    //pre-processing loop-most important part of the algorithm. The iterative  is repeated as long as correction still helds. The process
+    //might be described in 3 steps: correct diffusion, generate tensor, check how may negative eigenvalues are there.
 
     while (stil_correcting == true)
     {
-
-      std::cout << "Number of negative eigenvalues: " << old_number_negative_eigs << std::endl;// info for Thomas: Debug stuff - to be removed
+      //info for the user printed in the consol-debug information to be removed in the future
+      std::cout << "Number of negative eigenvalues: " << old_number_negative_eigs << std::endl;
 
       GenerateTensorImage(nof,numberb0,size,corrected_diffusion_temp,mask,what_mask,tensorImg);
 
       new_number_negative_eigs = CheckNegatives (size,mask,tensorImg);
 
       if(new_number_negative_eigs<old_number_negative_eigs)
       {
+        // if the correction does not introduce new negative eigenvalues and corrects some of negative detected n previous iteration
+        // smoothed DWI is used to substitute DWI from previous iteration
         stil_correcting=true;
         old_number_negative_eigs=new_number_negative_eigs;
         DeepCopyDiffusionImage(corrected_diffusion_temp,m_CorrectedDiffusionVolumes,nof);
       }
 
       else
       {
         stil_correcting=false;
       }
 
 
-      CorrectDiffusionImage(nof,numberb0,size,corrected_diffusion_temp,corrected_diffusion_temp,mask,pixel_max,pixel_min);
+      CorrectDiffusionImage(nof,numberb0,size,corrected_diffusion_temp,mask,pixel_max,pixel_min);
 
 
     }
 
+    //Changing the mask value for voxels that are still not corrcted. Original idea is to take them into analysis even though
+    // eigenvalues are still negative
     TurnMask(size, mask,1,1);
 
+    //Generation of final pre-processed tensor image that might be used as an input for FWE method
     GenerateTensorImage(nof,numberb0,size,m_CorrectedDiffusionVolumes,mask,what_mask,tensorImg);
 
     m_MaskImage = mask;
 
     this->SetNthOutput(0, tensorImg);
 
   }
 
 
 
   template <class TDiffusionPixelType, class TTensorPixelType>
   void
     TensorReconstructionWithEigenvalueCorrectionFilter<TDiffusionPixelType, TTensorPixelType>
     ::SetGradientImage( GradientDirectionContainerType *gradientDirection,
     const GradientImagesType *gradientImage )
   {
-    // Make sure crazy users did not call both AddGradientImage and
-    // SetGradientImage
+
     if( m_GradientImageTypeEnumeration == GradientIsInManyImages )
     {
       itkExceptionMacro( << "Cannot call both methods:"
         << "AddGradientImage and SetGradientImage. Please call only one of them.");
     }
 
     this->m_GradientDirectionContainer = gradientDirection;
 
 
     unsigned int numImages = gradientDirection->Size();
     this->m_NumberOfBaselineImages = 0;
 
     this->m_NumberOfGradientDirections = numImages - this->m_NumberOfBaselineImages;
 
     // ensure that the gradient image we received has as many components as
     // the number of gradient directions
     if( gradientImage->GetVectorLength() != this->m_NumberOfBaselineImages + this->m_NumberOfGradientDirections )
     {
       itkExceptionMacro( << this->m_NumberOfGradientDirections << " gradients + " << this->m_NumberOfBaselineImages
         << "baselines = " << this->m_NumberOfGradientDirections + this->m_NumberOfBaselineImages
         << " directions specified but image has " << gradientImage->GetVectorLength()
         << " components.");
     }
 
     this->ProcessObject::SetNthInput( 0,
       const_cast< GradientImagesType* >(gradientImage) );
     m_GradientImageTypeEnumeration = GradientIsInASingleImage;
   }
 
   template <class TDiffusionPixelType, class TTensorPixelType>
   double
   TensorReconstructionWithEigenvalueCorrectionFilter<TDiffusionPixelType, TTensorPixelType>
   ::CheckNeighbours(int x, int y, int z,int f, itk::Size<3> size, itk::Image<short, 3>::Pointer mask, itk::VectorImage<short, 3>::Pointer corrected_diffusion_temp)
   {
+    // method is used for finding a new value for the voxel with use of its 27 neighborhood. To perform such a smoothing correct voxels are
+    // counted an arithmetical mean is calculated and stored as a new value for the voxel. If there is no proper neigborhood voxel is turned
+    // to the value of 0.
 
     // Definition of neighbourhood avoiding crossing the image boundaries
     int x_max=size[0];
     int y_max=size[1];
     int z_max=size[2];
 
     double back_x=std::max(0,x-1);
     double back_y=std::max(0,y-1);
     double back_z=std::max(0,z-1);
 
     double forth_x=std::min((x+1),x_max);
     double forth_y=std::min((y+1),y_max);
     double forth_z=std::min((z+1),z_max);
 
 
     double tempsum=0;
     itk::Index<3> ix;
     double temp_number=0;
     double temp_mask=0;
 
     for(int i=back_x; i<=forth_x; i++)
     {
       for (int j=back_y; j<=forth_y; j++)
       {
         for (int k=back_z; k<=forth_z; k++)
         {
 
           ix[0] = i;
           ix[1] = j;
           ix[2] = k;
 
           temp_mask=mask->GetPixel(ix);
 
-          //if(temp_mask > 0 && temp_mask < 2 )
-          //{
-              //GradientVectorType p = m_GradientImagePointer->GetPixel(ix);
+
           GradientVectorType p = corrected_diffusion_temp->GetPixel(ix);
 
           double test= p[f];
 
-          if (test > 0.0 )// hmm this must be here becaus the method is used in multiple ocasions. Sometiems we may deal with negative values
+          if (test > 0.0 )// taking only positive values and counting them
           {
             if(!(i==x && j==y && k== z))
             {
-                tempsum=tempsum+p[f];// sum for calculation of mean
-                temp_number++;// number for calculation of mean
+                tempsum=tempsum+p[f];
+                temp_number++;
             }
 
           }
 
 
         }
       }
-    }// end of size
-
+    }// end of jumping through 27th neighbours
 
+    //getting back to the original position of the voxel
 
     ix[0] = x;ix[1] = y;ix[2] = z;
 
     if (temp_number <= 0.0)
     {
-      //ix[0] = x;ix[1] = y;ix[2] = z;
-      //GradientVectorType p = m_GradientImagePointer->GetPixel(ix);
-      //tempsum=p[f];
       tempsum=0;
       mask->SetPixel(ix,0);
     }
     else
     {
       tempsum=tempsum/temp_number;
 
     }
 
 
-    return tempsum;
+    return tempsum;// smoothed value of voxel
   }
 
 
   template <class TDiffusionPixelType, class TTensorPixelType>
   void
   TensorReconstructionWithEigenvalueCorrectionFilter<TDiffusionPixelType, TTensorPixelType>
   ::CalculateAttenuation(vnl_vector<double> org_data,vnl_vector<double> &atten,int nof, int numberb0)
   {
     double mean_b=0.0;
 
     for (int i=0;i<nof;i++)
     {
       if(m_B0Mask[i]>0)
       {
         double o_d=org_data[i];
         mean_b=mean_b+org_data[i];
       }
 
     }
     mean_b=mean_b/numberb0;
     int cnt=0;
     for (int i=0;i<nof;i++)
     {
       if(m_B0Mask[i]==0)
       {
         //if(org_data[i]<0.001){org_data[i]=0.01;}
         atten[cnt]=org_data[i]/mean_b;
         cnt++;
       }
     }
   }
 
-
-
-  /*
-  template <class TDiffusionPixelType, class TTensorPixelType>
-  void
-  TensorReconstructionWithEigenvalueCorrectionFilter<TDiffusionPixelType, TTensorPixelType>
-  ::CalculateTensor(vnl_vector<double> atten,vnl_vector<double> &tensor, int nof,int numberb0)
-  {
-    for (int i=0;i<nof-numberb0;i++)
-    {
-
-      if(atten[i]<=0.0)
-      {
-          atten[i]==0.1;
-      }
-     std::cout << "problem";
-      atten[i]=log((double)atten[i]);
-    }
-    tensor = m_PseudoInverse*atten;
-
-
-  }
-
-*/
-
-
   template <class TDiffusionPixelType, class TTensorPixelType>
   double
   TensorReconstructionWithEigenvalueCorrectionFilter<TDiffusionPixelType, TTensorPixelType>
   ::CheckNegatives ( itk::Size<3> size, itk::Image<short, 3>::Pointer mask, itk::Image< itk::DiffusionTensor3D<float>, 3 >::Pointer tensorImg )
   {
 
       // The method was created to simplif the flow of negative eigenvalue correction process. The method itself just return the number
       // of voxels (tensors) with negative eigenvalues. Then if the voxel was previously bad ( mask=2 ) but it is not bad anymore mask is
       //changed to 1.
 
+      // declaration of important structures and variables
       double badvoxels=0;
       double pixel=0;
       itk::Index<3> ix;
       itk::DiffusionTensor3D<float> ten;
       vnl_matrix<double> temp_tensor(3,3);
-      vnl_vector<double> eigen_vals(3);// declaration of important structures and variables
+      vnl_vector<double> eigen_vals(3);
       vnl_vector<double> tensor (6);
 
-
+      // for every pixel from the image
       for (int x=0;x<size[0];x++)
       {
 
         for (int y=0;y<size[1];y++)
         {
 
           for (int z=0;z<size[2];z++)
           {
-              // for every pixel
+
 
               ix[0] = x; ix[1] = y; ix[2] = z;
 
               pixel = mask->GetPixel(ix);
 
+              // but only if previously marked as bad one-negative eigen value
+
               if(pixel > 1)
-              {// but only if previously marked as bad or potentially bad one
+              {
 
               ten = tensorImg->GetPixel(ix);
 
+              // changing order from tensor structure in MITK into vnl structure 3x3 symmetric tensor matrix
+
               tensor[0] = ten(0,0);
               tensor[3] = ten(0,1);
               tensor[5] = ten(0,2);
               tensor[1] = ten(1,1);
               tensor[4] = ten(1,2);
               tensor[2] = ten(2,2);
 
               temp_tensor[0][0]= tensor[0]; temp_tensor[1][0]= tensor[3]; temp_tensor[2][0]= tensor[5];
               temp_tensor[0][1]= tensor[3]; temp_tensor[1][1]= tensor[1]; temp_tensor[2][1]= tensor[4];
               temp_tensor[0][2]= tensor[5]; temp_tensor[1][2]= tensor[4]; temp_tensor[2][2]= tensor[2];
 
+              //checking negativity of tensor eigenvalues
 
-              vnl_symmetric_eigensystem<double> eigen_tensor(temp_tensor);//creating an eigensystem out of 3x3 symmetric tensor matrix
+              vnl_symmetric_eigensystem<double> eigen_tensor(temp_tensor);
 
 
               eigen_vals[0]=eigen_tensor.get_eigenvalue(0);
               eigen_vals[1]=eigen_tensor.get_eigenvalue(1);
-              eigen_vals[2]=eigen_tensor.get_eigenvalue(2);// calculating eigenvalues for current tensor
+              eigen_vals[2]=eigen_tensor.get_eigenvalue(2);
+
+              //comparison to 0.01 instead of 0 was proposed by O.Pasternak
 
               if( eigen_vals[0]>0.01 && eigen_vals[1]>0.01 && eigen_vals[2]>0.01)
               {
                 mask->SetPixel(ix,1);
 
               }
               else
               {
-                badvoxels++;// increasing a number of detected bad
+                badvoxels++;
               }
 
               }
 
           }
          }
        }
 
-      double ret = badvoxels;//returning just the number of bad voxels
+      double ret = badvoxels;
 
       return ret;
 
-  }// end of void check negativity
+  }
 
 
   template <class TDiffusionPixelType, class TTensorPixelType>
   void
   TensorReconstructionWithEigenvalueCorrectionFilter<TDiffusionPixelType, TTensorPixelType>
-  ::CorrectDiffusionImage(int nof,int numberb0,itk::Size<3> size,itk::VectorImage<short, 3>::Pointer corrected_diffusion_temp,itk::VectorImage<short, 3>::Pointer corrected_diffusion,itk::Image<short, 3>::Pointer mask,vnl_vector< double> pixel_max,vnl_vector< double> pixel_min)
+  ::CorrectDiffusionImage(int nof,int numberb0,itk::Size<3> size,itk::VectorImage<short, 3>::Pointer corrected_diffusion,itk::Image<short, 3>::Pointer mask,vnl_vector< double> pixel_max,vnl_vector< double> pixel_min)
   {
-      // in this method the diffusion image temporray is corrected with use of information from updated mask. The diffusion image is tempora
-      //-ry while we can't use corrected voxel to correct other voxel in the same iteration.
+      // in this method the voxels that has tensor negative eigenvalues are smoothed. Smoothing is done on DWI image.For the voxel
+      //detected as bad one, B0 image is smoothed obligatory. All other gradient images are smoothed only when value of attenuation
+      //is out of declared bounds for too high or too low attenuation.
 
+     // declaration of important variables
 
       itk::Index<3> ix;
       vnl_vector<double> org_data(nof-numberb0);
       vnl_vector<double> atten(nof-numberb0);
-      double cnt_atten=0;// declaration of important variables
+      double cnt_atten=0;
 
       for (int z=0;z<size[2];z++)
       {
 
         for (int x=0;x<size[0];x++)
         {
 
           for (int y=0;y<size[1];y++)
           {
-              // for every pixel
+
 
               ix[0] = x; ix[1] = y; ix[2] = z;
 
 
               if(mask->GetPixel(ix) > 1.0)
               {
                   GradientVectorType  pt = corrected_diffusion->GetPixel(ix);
 
                   for (int i=0;i<nof;i++)
                   {
                     org_data[i]=pt[i];
                   }
 
-                  //CalculateAttenuation(org_data,atten,nof,numberb0);
 
                   double mean_b=0.0;
 
                   for (int i=0;i<nof;i++)
                   {
                     if(m_B0Mask[i]>0)
                     {
                       mean_b=mean_b+org_data[i];
                     }
 
                    }
                   mean_b=mean_b/numberb0;
                   int cnt=0;
                   for (int i=0;i<nof;i++)
                   {
                     if(m_B0Mask[i]==0)
                     {
-                      //if(org_data[i]<0.001){org_data[i]=0.01;}
+
                       atten[cnt]=org_data[i]/mean_b;
                       cnt++;
                     }
                   }
 
                   cnt_atten=0;
 
+                  //smoothing certain gradient images taht are out of declared constraints
+
                   for (int f=0;f<nof;f++)
                   {
                     if(m_B0Mask[f]==0)
                     {
 
                         if(atten[cnt_atten]<pixel_min[cnt_atten] || atten[cnt_atten]> pixel_max[cnt_atten])
                     {
-                        org_data[f] = CheckNeighbours(x,y,z,f,size,mask,corrected_diffusion_temp);
+                        org_data[f] = CheckNeighbours(x,y,z,f,size,mask,corrected_diffusion);
 
                     }
 
 
                     cnt_atten++;
 
                    }
+                    //smoothing B0
                     if(m_B0Mask[f]==1)
                     {
 
-
-                        org_data[f] = CheckNeighbours(x,y,z,f,size,mask,corrected_diffusion_temp);
-
+                        org_data[f] = CheckNeighbours(x,y,z,f,size,mask,corrected_diffusion);
 
                     }
 
 
-                  }//end for
-
-
+                  }
 
                   for (int i=0;i<nof;i++)
                   {
                     pt[i]=org_data[i];
                   }
 
                   corrected_diffusion_temp->SetPixel(ix, pt);
 
               }
               else
               {
                   GradientVectorType  pt = corrected_diffusion->GetPixel(ix);
-                  corrected_diffusion_temp->SetPixel(ix, pt);
+                  corrected_diffusion->SetPixel(ix, pt);
 
               }
            }
         }
        }
 
-      std::cout << "break point";
+
 
   }
 
   template <class TDiffusionPixelType, class TTensorPixelType>
   void
   TensorReconstructionWithEigenvalueCorrectionFilter<TDiffusionPixelType, TTensorPixelType>
   ::GenerateTensorImage(int nof,int numberb0,itk::Size<3> size,itk::VectorImage<short, 3>::Pointer corrected_diffusion,itk::Image<short, 3>::Pointer mask,double what_mask,itk::Image< itk::DiffusionTensor3D<float>, 3 >::Pointer tensorImg)
   {
       // in this method the whole tensor image is updated with a tensors for defined voxels ( defined by a value of mask);
 
 
       itk::Index<3> ix;
       vnl_vector<double> org_data(nof-numberb0);
       vnl_vector<double> atten(nof-numberb0);
       vnl_vector<double> tensor(6);
       itk::DiffusionTensor3D<float> ten;
       double mask_val=0;
 
 
       for (int x=0;x<size[0];x++)
       {
 
         for (int y=0;y<size[1];y++)
         {
 
          for (int z=0;z<size[2];z++)
           {
-              // for every pixel
+
 
               ix[0] = x; ix[1] = y; ix[2] = z;
 
               mask_val= mask->GetPixel(ix);
 
+              //Tensors are calculated only for voxels above theshold for B0 image.
+
               if( mask_val > 0.0 )
               {
+
+                 // calculation of attenuation with use of gradient image and  and mean B0 image
                  GradientVectorType pt = corrected_diffusion->GetPixel(ix);
 
                   for (int i=0;i<nof;i++)
                   {
                     org_data[i]=pt[i];
                   }
 
                   double mean_b=0.0;
 
                   for (int i=0;i<nof;i++)
                   {
                     if(m_B0Mask[i]>0)
                     {
                       double o_d=org_data[i];
                       mean_b=mean_b+org_data[i];
                     }
 
                   }
                   mean_b=mean_b/numberb0;
                   int cnt=0;
                   for (int i=0;i<nof;i++)
                   {
                     if (org_data[i]<= 0)
 
                     {
                         org_data[i]=0.1;
                     }
                     if(m_B0Mask[i]==0)
                     {
                       atten[cnt]=org_data[i]/mean_b;
                       cnt++;
                     }
                   }
 
 
 
 
                   for (int i=0;i<nof-numberb0;i++)
                   {
 
-                    if(atten[i]==0.0)
-                   std::cout << "problem";
                     atten[i]=log((double)atten[i]);
                   }
+
+                  // Calculation of tensor with use of previously calculated inverse of design matrix and attenuation
                   tensor = m_PseudoInverse*atten;
 
                   ten(0,0) = tensor[0];
                   ten(0,1) = tensor[3];
                   ten(0,2) = tensor[5];
                   ten(1,1) = tensor[1];
                   ten(1,2) = tensor[4];
                   ten(2,2) = tensor[2];
 
                   tensorImg->SetPixel(ix, ten);
 
 
 
               }
+              // for voxels with mask value 0 - tensor is simply 0 ( outside brain value)
               else if (mask_val < 1.0)
               {
                   ten(0,0) = 0;
                   ten(0,1) = 0;
                   ten(0,2) = 0;
                   ten(1,1) = 0;
                   ten(1,2) = 0;
                   ten(2,2) = 0;
                   tensorImg->SetPixel(ix, ten);
               }
 
-
-
-
            }
         }
        }
 
 
 
   }// end of Generate Tensor
 
 
   template <class TDiffusionPixelType, class TTensorPixelType>
   void
   TensorReconstructionWithEigenvalueCorrectionFilter<TDiffusionPixelType, TTensorPixelType>
   ::TurnMask( itk::Size<3> size, itk::Image<short, 3>::Pointer mask, double previous_mask, double set_mask)
   {
-    // in this method voxels in the mask that poses a certain value are substituded by other value.// in this method voxels in the mask that poses a certain value are substituded by other value.
+    // The method changes voxels in the mask that poses a certain value with other value.
+
     itk::Index<3> ix;
     double temp_mask_value=0;
 
     for(int x=0;x<size[0];x++)
     {
       for(int y=0;y<size[1];y++)
       {
         for(int z=0;z<size[2];z++)
         {
           ix[0] = x; ix[1] = y; ix[2] = z;
           temp_mask_value=mask->GetPixel(ix);
 
           if(temp_mask_value>previous_mask)
           {
             mask->SetPixel(ix,set_mask);
           }
         }
       }
     }
   }
 
 
 
   template <class TDiffusionPixelType, class TTensorPixelType>
   void
   TensorReconstructionWithEigenvalueCorrectionFilter<TDiffusionPixelType, TTensorPixelType>
   ::DeepCopyDiffusionImage(itk::VectorImage<short, 3>::Pointer corrected_diffusion, itk::VectorImage<short, 3>::Pointer corrected_diffusion_temp,int nof)
   {
     corrected_diffusion_temp->SetSpacing(corrected_diffusion->GetSpacing());
     corrected_diffusion_temp->SetOrigin(corrected_diffusion->GetOrigin());
     corrected_diffusion_temp->SetVectorLength(nof);
     corrected_diffusion_temp->SetRegions(corrected_diffusion->GetLargestPossibleRegion());
     corrected_diffusion_temp->Allocate();
 
     itk::ImageRegionConstIterator<ImageType> inputIterator(corrected_diffusion, corrected_diffusion->GetLargestPossibleRegion());
     itk::ImageRegionIterator<ImageType> outputIterator(corrected_diffusion_temp, corrected_diffusion_temp->GetLargestPossibleRegion());
 
     inputIterator.GoToBegin();
     outputIterator.GoToBegin();
 
     while(!inputIterator.IsAtEnd())
       {
       outputIterator.Set(inputIterator.Get());
       ++inputIterator;
       ++outputIterator;
       }
   }
 
   template <class TDiffusionPixelType, class TTensorPixelType>
   void
   TensorReconstructionWithEigenvalueCorrectionFilter<TDiffusionPixelType, TTensorPixelType>
   ::DeepCopyTensorImage(itk::Image< itk::DiffusionTensor3D<float>, 3 >::Pointer tensorImg, itk::Image< itk::DiffusionTensor3D<float>, 3 >::Pointer temp_tensorImg)
   {
 
       temp_tensorImg->SetSpacing(tensorImg->GetSpacing());
       temp_tensorImg->SetOrigin(tensorImg->GetOrigin());
       temp_tensorImg->SetRegions(tensorImg->GetLargestPossibleRegion());
       temp_tensorImg->Allocate();
 
     itk::ImageRegionConstIterator<TensorImageType> inputIterator(tensorImg, tensorImg->GetLargestPossibleRegion());
     itk::ImageRegionIterator<TensorImageType> outputIterator(temp_tensorImg, temp_tensorImg->GetLargestPossibleRegion());
 
     inputIterator.GoToBegin();
     outputIterator.GoToBegin();
 
     while(!inputIterator.IsAtEnd())
       {
       outputIterator.Set(inputIterator.Get());
       ++inputIterator;
       ++outputIterator;
       }
   }
 
 
 
 } // end of namespace