vtkImagePCAFilter.h

Go to the documentation of this file.
00001 /*=auto=========================================================================
00002 
00003 (c) Copyright 2001 Massachusetts Institute of Technology
00004 
00005 Permission is hereby granted, without payment, to copy, modify, display
00006 and distribute this software and its documentation, if any, for any purpose,
00007 provided that the above copyright notice and the following three paragraphs
00008 appear on all copies of this software.  Use of this software constitutes
00009 acceptance of these terms and conditions.
00010 
00011 IN NO EVENT SHALL MIT BE LIABLE TO ANY PARTY FOR DIRECT, INDIRECT, SPECIAL,
00012 INCIDENTAL, OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE OF THIS SOFTWARE
00013 AND ITS DOCUMENTATION, EVEN IF MIT HAS BEEN ADVISED OF THE POSSIBILITY OF
00014 SUCH DAMAGE.
00015 
00016 MIT SPECIFICALLY DISCLAIMS ANY EXPRESS OR IMPLIED WARRANTIES INCLUDING,
00017 BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY, FITNESS FOR
00018 A PARTICULAR PURPOSE, AND NON-INFRINGEMENT.
00019 
00020 THE SOFTWARE IS PROVIDED "AS IS."  MIT HAS NO OBLIGATION TO PROVIDE
00021 MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.
00022 
00023 =========================================================================auto=*/
00024 // .NAME vtkImagePCAFilter
00025 // This is a generic class for vtkImagePCATraining and vtkImagePCAAnalysis
00026 #ifndef __vtkImagePCAFilter_h
00027 #define __vtkImagePCAFilter_h
00028 
00029 #include "vtkAtlasCreatorCxxModuleWin32Header.h"
00030 
00031 // #include "vtkSlicer.h"
00032 #include "vtkImageMultipleInputFilter.h"
00033 #include "vtkFloatArray.h"
00034 
00035 class VTK_ATLASCREATORCXXMODULE_EXPORT vtkImagePCAFilter : public vtkImageMultipleInputFilter
00036 {
00037   public:
00038   // -----------------------------------------------------
00039   // Genral Functions for the filter
00040   // -----------------------------------------------------
00041   static vtkImagePCAFilter *New();
00042   vtkTypeMacro(vtkImagePCAFilter,vtkObject);
00043   void PrintSelf(ostream& os, vtkIndent indent);
00044 
00045   void SetInputIndex(int index, vtkImageData *image) {this->SetInput(index,image);}
00046 
00047   vtkGetMacro(NumberOfEigenModes, int);
00048 
00049   vtkGetMacro(XDim, int);
00050   vtkGetMacro(YDim, int);
00051   vtkGetMacro(ZDim, int);
00052 
00053   float GetEigenValue(int index) {
00054     if (index < this->NumberOfEigenModes) return this->EigenValues[index];
00055     else {
00056       vtkWarningMacro( << "GetEigenValue: Index exceeds dimension of EigenValue" );
00057       return 0;
00058     }
00059   }
00060   void SetEigenValue(int index, float value) {
00061      if (index < this->NumberOfEigenModes) this->EigenValues[index] = value;
00062      else vtkWarningMacro( << "SetEigenValue: Index (" << index << ") exceeds Number of EigenModes ("<< this->NumberOfEigenModes << ") !" );
00063   }
00064 
00065   void TransfereLabelmapIntoPCADistanceMap(vtkImageData *inData, int Label, float MaxDistanceValue, float BoundaryValue, vtkImageData *outData);
00066 
00067   // Checks the Input if everything is set correctly
00068   int CheckInput(vtkImageData **InData);
00069   int CheckInput(vtkImageData *InData);
00070 
00071 protected:
00072   vtkImagePCAFilter();
00073   vtkImagePCAFilter(const vtkImagePCAFilter&) {};
00074   ~vtkImagePCAFilter();
00075 
00076   void operator=(const vtkImagePCAFilter&) {};
00077   void ThreadedExecute(vtkImageData **inData, vtkImageData *outData,int outExt[6], int id){};
00078 
00079   void DeleteEigenValues();
00080 
00081   void SetNumberOfEigenModes(int init) {
00082     this->DeleteEigenValues();
00083     this->NumberOfEigenModes =  init;
00084     this->EigenValues = new float[init];
00085     memset(this->EigenValues,0,sizeof(float)*init);
00086   }
00087 
00088   int CheckEigenValues() {
00089     int i = 0;
00090     while ((i < this->NumberOfEigenModes) && (this->EigenValues[i] > 0.0)) i++;
00091     return (i < this->NumberOfEigenModes ? 1 : 0);
00092   }
00093 
00094   int NumberOfEigenModes;
00095   float *EigenValues;
00096 
00097   // Save size so that everytime it is the same
00098   int XDim;
00099   int YDim;
00100   int ZDim;
00101   int InputScalarType;
00102 
00103 };
00104 
00105 //BTX
00106 //------------------------------------------------------------------------
00107 // Matrix ops. Some taken from vtkThinPlateSplineTransform.cxx
00108 static inline double** NewMatrix(int rows, int cols)
00109 {
00110   double *matrix = new double[rows*cols];
00111   double **m = new double *[rows];
00112   for(int i = 0; i < rows; i++) {
00113     m[i] = &matrix[i*cols];
00114   }
00115   return m;
00116 }
00117 
00118 //------------------------------------------------------------------------
00119 static inline void DeleteMatrix(double **m)
00120 {
00121   delete [] *m;
00122   delete [] m;
00123 }
00124 
00125 //------------------------------------------------------------------------
00126 template <class TMatrix1>
00127 static inline void MatrixMultiply(TMatrix1 **a, double **b, double **c,
00128                                   int arows, int acols,
00129                                   int brows, int bcols)
00130 {
00131   if(acols != brows)  {
00132     return; // acols must equal br otherwise we can't proceed
00133   }
00134 
00135   // c must have size arows*bcols (we assume this)
00136 
00137   for(int i = 0; i < arows; i++) {
00138     for(int j = 0; j < bcols; j++) {
00139       c[i][j] = 0.0;
00140       for(int k = 0; k < acols; k++) {
00141         c[i][j] += double(a[i][k])*b[k][j];
00142       }
00143     }
00144   }
00145 }
00146 
00147 
00148 //------------------------------------------------------------------------
00149 // Subtracting the mean column from the observation matrix is equal
00150 // to subtracting the mean shape from all shapes.
00151 // The mean column is equal to the Procrustes mean (it is also returned)
00152 template <class TMatrix1, class TMatrix2>
00153 static inline void SubtractMeanColumn(TMatrix1 **m, TMatrix2 *mean, int rows, int cols)
00154 {
00155   int r,c;
00156   double csum;
00157   for (r = 0; r < rows; r++) {
00158     // What does it mean also done at Simon
00159     csum = 0.0F;
00160     for (c = 0; c < cols; c++) {
00161       csum += double(m[r][c]);
00162     }
00163     // calculate average value of row
00164     csum /= cols;
00165 
00166     // Mean shape vector is updated
00167     mean[r] =  TMatrix2(csum);
00168 
00169     // average value is subtracted from all elements in the row
00170     for (c = 0; c < cols; c++) {
00171       m[r][c] -=  TMatrix1(csum);
00172     }
00173   }
00174 }
00175 
00176 //------------------------------------------------------------------------
00177 // Normalise all columns to have length 1
00178 // meaning that all eigenvectors are normalised
00179 static inline void NormaliseColumns(double **m, int rows, int cols)
00180 {
00181   for (int c = 0; c < cols; c++) {
00182     double cl = 0;
00183     for (int r = 0; r < rows; r++) {
00184       cl += m[r][c] * m[r][c];
00185     }
00186     cl = sqrt(cl);
00187 
00188     // If cl == 0 something is rotten, dont do anything now
00189     if (cl != 0) {
00190       for (int r = 0; r < rows; r++) {
00191         m[r][c] /= cl;
00192       }
00193     }
00194   }
00195 }
00196 
00197 //------------------------------------------------------------------------
00198 // Here it is assumed that a rows >> a cols
00199 // Output matrix is [a cols X a cols]
00200 static inline void SmallCovarianceMatrix(double **a, double **c,
00201                                          int arows, int acols)
00202 {
00203   const int s = acols;
00204 
00205   // c must have size acols*acols (we assume this)
00206   for(int i = 0; i < acols; i++) {
00207     for(int j = 0; j < acols; j++) {
00208       // Use symmetry
00209       if (i <= j) {
00210         c[i][j] = 0.0;
00211         for(int k = 0; k < arows; k++) {
00212           c[i][j] += a[k][i]*a[k][j];
00213         }
00214         c[i][j] /= (s-1);
00215         c[j][i] = c[i][j];
00216       }
00217     }
00218   }
00219 }
00220 
00221 //------------------------------------------------------------------------
00222 static inline double* NewVector(int length)
00223 {
00224   double *vec = new double[length];
00225   return vec;
00226 }
00227 
00228 //------------------------------------------------------------------------
00229 static inline void DeleteVector(double *v)
00230 {
00231   delete [] v;
00232 }
00233 //ETX
00234 #endif
00235 
00236 
00237 
00238 
00239 
00240 
00241 
00242 
00243 
00244 
00245 

Generated on 6 Apr 2011 for Slicer3 by  doxygen 1.6.1