EMLocalInterface.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 EMLocalInterface
00025 
00026 #ifndef _EMLOCALINTERFACE_H_INCLUDED
00027 #define _EMLOCALINTERFACE_H_INCLUDED 1
00028 
00029 #include "vtkEMSegment.h"
00030 #include "vtkMultiThreader.h"
00031 
00032 // Defines the maximum number of threads used throughout this program.
00033 // This is very helpful when validating results from machines with different cpu number
00034 // EMLOCALSEGMENTER_MAX_MULTI_THREAD = -1 => no restriction
00035  
00036 #define EMLOCALSEGMENTER_MAX_MULTI_THREAD -1
00037 
00038 //--------------------------------------------------------------------
00039 // Registration Parameters 
00040 //--------------------------------------------------------------------
00041 // Parameters needed to define cost function 
00042 #define EMSEGMENT_REGISTRATION_DISABLED     0
00043 #define EMSEGMENT_REGISTRATION_APPLY        1 
00044 #define EMSEGMENT_REGISTRATION_GLOBAL_ONLY  2
00045 #define EMSEGMENT_REGISTRATION_CLASS_ONLY   3
00046 #define EMSEGMENT_REGISTRATION_SIMULTANEOUS 4
00047 #define EMSEGMENT_REGISTRATION_SEQUENTIAL   5
00048  
00049 #define EMSEGMENT_REGISTRATION_INTERPOLATION_LINEAR 1
00050 #define EMSEGMENT_REGISTRATION_INTERPOLATION_NEIGHBOUR 2 
00051 
00052 #define EMSEGMENT_REGISTRATION_SIMPLEX 1 
00053 #define EMSEGMENT_REGISTRATION_POWELL 2 
00054 
00055 #define EMSEGMENT_STOP_FIXED    0 
00056 #define EMSEGMENT_STOP_LABELMAP 1 
00057 #define EMSEGMENT_STOP_WEIGHTS  2 
00058 
00059 #define EMSEGMENT_PCASHAPE_DEPENDENT 0
00060 #define EMSEGMENT_PCASHAPE_INDEPENDENT 1
00061 #define EMSEGMENT_PCASHAPE_APPLY 2
00062 
00063 
00064 //--------------------------------------------------------------------
00065 // Hierachy / SuperClass specific parameters 
00066 //--------------------------------------------------------------------
00067 //BTX  
00068 class VTK_EMSEGMENT_EXPORT EMLocal_Hierarchical_Class_Parameters {
00069  public: 
00070   int           NumClasses;
00071   int           NumTotalTypeCLASS;
00072   int*          NumChildClasses;
00073   // Do not put it here even though it makees sense bc in MF changes for each job
00074   // void          ** ProbDataPtr;
00075   int           *ProbDataIncY; 
00076   int           *ProbDataIncZ; 
00077   float         *ProbDataWeight;
00078   float         *ProbDataMinusWeight;
00079   int           ProbDataType;
00080   double        **LogMu;
00081   double        ***InvLogCov;
00082   double        *InvSqrtDetLogCov;
00083   double        *TissueProbability;
00084   int           *VirtualNumInputImages;
00085   double        ***MrfParams;
00086   EMLocal_Hierarchical_Class_Parameters(); 
00087   ~EMLocal_Hierarchical_Class_Parameters() {} 
00088   void Copy(EMLocal_Hierarchical_Class_Parameters init) ; 
00089 }; 
00090 
00091 //  Utility function used by shape.cxx and EMPrivatRegistration - should be put in .h file
00092 inline int EMLocalInterface_InterpolationNearestNeighbourVoxelIndex(float col, float row, float slice, int dataIncY,int dataIncZ, int* Image_Length) {
00093   /* AVOID floor() overhead */
00094   /* Is floor call necessary ?
00095      i = (int)floor(row + .5);
00096      j = (int)floor(col + .5);
00097      k = (int)floor(slice + .5);
00098   */
00099   int i = (int)((row < 0) ? (row - 0.5) : (row + 0.5));
00100   int j = (int)((col < 0) ? (col - 0.5) : (col + 0.5));
00101   int k = (int)((slice < 0) ? (slice - 0.5) : (slice + 0.5));
00102 
00103   int MaxIndexCol   = Image_Length[0] - 1;
00104   int MaxIndexRow   = Image_Length[1] - 1;
00105   int MaxIndexSlice = Image_Length[2] - 1;
00106 
00107   if (i < 0) {i = 0;}
00108   else if (i > MaxIndexRow) i = MaxIndexRow;  
00109   if (j < 0) j = 0;
00110   else if (j > MaxIndexCol) j =  MaxIndexCol;
00111   if (k < 0) k = 0;
00112   else if (k > MaxIndexSlice) k = MaxIndexSlice;
00113 
00114   int colsize   = Image_Length[0] +dataIncY; 
00115   int slicesize = Image_Length[1]*colsize + dataIncZ;
00116 
00117   return i*colsize + j + k*slicesize;
00118 }    
00119 
00120 template <class T>
00121 inline double  EMLocalInterface_Interpolation(
00122     float col, float row,   float slice,
00123     int ncol,  int nrow,  int nslice,
00124     T* data, int dataIncY, int dataIncZ, int InterpolationType, int* Image_Length ) {
00125 
00126 
00127   int j = (int)floor(col);
00128   int i = (int)floor(row);
00129   int k = (int)floor(slice);
00130 
00131   /* Decide if we can do trilinear interpolation */
00132   // ------------------------------------------------
00133   // from NearestNeighbourInterpolation 
00134   // ------------------------------------------------
00135 
00136   if ((InterpolationType == EMSEGMENT_REGISTRATION_INTERPOLATION_NEIGHBOUR) || (i < 0) || (j < 0) || (k < 0) || (i >= (nrow - 1)) || (j >= (ncol - 1)) || ((k >= (nslice - 1) ) && (nslice != 1))) {
00137     int VoxelJump = EMLocalInterface_InterpolationNearestNeighbourVoxelIndex(col, row, slice, dataIncY,dataIncZ, Image_Length);
00138     return double(data[VoxelJump]);
00139   }
00140 
00141   // ------------------------------------------------
00142   // From LinearInterpolation
00143   // ------------------------------------------------
00144   int index;
00145 
00146   int colsize   = ncol+dataIncY; 
00147   int slicesize = nrow*colsize + dataIncZ;
00148 
00149   double t;
00150   double tt;
00151   double u;
00152   double uu;
00153   double v;
00154   double vv;
00155 
00156   t = double(row - i);
00157   tt = double(1.0 - t);
00158   u = double(col - j);
00159   uu = double(1.0 - u);
00160   v = double(slice - k);
00161   vv = double(1.0 - v);
00162 
00163   if (k >= (nslice - 1) && (nslice == 1)) {
00164     /* we have 2D data - so do a restricted interpolation */
00165 
00166     /* So off slice will make no difference */
00167     v = 0.0;
00168     vv = 1.0; 
00169     /* So we don't dereference an illegal location */
00170     slicesize = 0;
00171   }
00172 
00173 /*
00174 fprintf(stdout,"row %g col %g slice %g ", row, col, slice);
00175 fprintf(stdout,"(i,j,k) (%d, %d, %d)\n", i, j, k);
00176 fprintf(stdout, "t %g tt %g u %g uu %g v %g vv %g\n", t, tt, u, uu, v, vv);
00177 */
00178 
00179   index = i*colsize + j + k*slicesize;
00180 /*
00181 fprintf(stdout, "centre index: %d lastindex: %d\n", index, lastindex);
00182 fprintf(stdout, "data indexed: %d %d %d %d %d %d %d\n", index + ncol,
00183 index + 1, index + 1 + ncol, index + slicesize, index + ncol + slicesize, 
00184 index + 1 + slicesize, index + 1 + ncol + slicesize);
00185 */
00186 
00187 
00188  return double(tt*uu*vv*double(data[index])     + t*uu*vv*double(data[index + colsize]) + 
00189      tt*u*vv*double(data[index + 1])            + t*u*vv*double(data[index + 1 + colsize]) + 
00190      tt*uu*v*double(data[index + slicesize])    + t*uu*v*double(data[index + colsize + slicesize]) +
00191      tt*u*v*double(data[index + 1 + slicesize]) + t*u*v*double(data[index + 1 + colsize + slicesize]));
00192 
00193 }
00194 
00195 
00196 template <class T>
00197 inline void  EMLocalInterface_Interpolation(float col, float row, float slice, int ncol,  int nrow, int nslice, T* data, int dataIncY, int dataIncZ, 
00198                                               int InterpolationType, int* Image_Length, double &result) {
00199      result = EMLocalInterface_Interpolation(col, row, slice, ncol, nrow, nslice, data, dataIncY, dataIncZ, InterpolationType, Image_Length);
00200 }
00201 
00202 // This interpolation descirbes the following:
00203 // ccol, crow, cslice are in image space 
00204 // The function calculates the coordinates in target space (in our case the atlas space)
00205 inline void EMLocalInterface_findCoordInTargetOfMatchingSourceCentreTarget(
00206     float *invRot, float *invTran,
00207     int ccol, int crow,   int cslice,
00208     float &ocol,  float &orow,  float &oslice,
00209     float Image_MidX, float Image_MidY, float Image_MidZ) {
00210 
00211   float mmcoords[3];
00212 
00213   // Kilian: Scale has to be applied this way so that the center stays the same
00214   mmcoords[0] = float(ccol)   - Image_MidX;
00215   mmcoords[1] = float(crow)   - Image_MidY;
00216   mmcoords[2] = float(cslice) - Image_MidZ;
00217   ocol    = invRot[0]*mmcoords[0] +
00218                      invRot[1]*mmcoords[1]+                     
00219                      invRot[2]*mmcoords[2] + invTran[0] + Image_MidX;
00220   orow    = invRot[3]*mmcoords[0] +
00221                      invRot[4]*mmcoords[1] +
00222                      invRot[5]*mmcoords[2] + invTran[1] + Image_MidY;
00223   oslice  = invRot[6]*mmcoords[0] +
00224                      invRot[7]*mmcoords[1] +
00225                      invRot[8]*mmcoords[2] + invTran[2] + Image_MidZ;
00226 
00227 }
00228 
00229 
00230 // ==============================================
00231 // Shape Based Functions 
00232 // This is currnelty without scaling of the PCA parameters 
00233 // Voxel jump is needed so we can do registration easily
00234 
00235 inline float  EMLocalInterface_CalcDistanceMap(const double *ShapeParameter, float **PCAEigenVectorsPtr, float *MeanShape, int DimEigenVector, int VoxelJump) {
00236   float term = *(MeanShape + VoxelJump); 
00237   for (int l = 0 ; l < DimEigenVector; l++) term += (*(PCAEigenVectorsPtr[l] + VoxelJump)) * ShapeParameter[l];
00238   return  term;
00239 }
00240 
00241 inline float  EMLocalInterface_CalcDistanceMap(const float *ShapeParameter, float **PCAEigenVectorsPtr, float *MeanShape, int DimEigenVector, int VoxelJump) {
00242   float term = *(MeanShape + VoxelJump); 
00243   for (int l = 0 ; l < DimEigenVector; l++) term += (*(PCAEigenVectorsPtr[l] + VoxelJump)) * ShapeParameter[l];
00244   return  term;
00245 }
00246 
00247 inline int EMLocalInterface_DefineMultiThreadJump(int *Start, int BoundaryMaxX, int BoundaryMaxY, int DataIncY,  int DataIncZ) {
00248   int LengthOfXDim = BoundaryMaxX + DataIncY;
00249   int LengthOfYDim = LengthOfXDim*(BoundaryMaxY) + DataIncZ;  
00250   return Start[0] + Start[1] * LengthOfXDim + LengthOfYDim *Start[2];
00251 }
00252 
00253 inline int EMLocalInterface_GetDefaultNumberOfThreads(int DisableFlag) {
00254   if (DisableFlag) return 1;
00255   int result = vtkMultiThreader::GetGlobalDefaultNumberOfThreads();
00256   if ((EMLOCALSEGMENTER_MAX_MULTI_THREAD > 0 ) && (result >  EMLOCALSEGMENTER_MAX_MULTI_THREAD)) return EMLOCALSEGMENTER_MAX_MULTI_THREAD;
00257   return result;
00258 }
00259 #endif

Generated on 6 Apr 2011 for Slicer3 by  doxygen 1.6.1