added glcm computation

image/MultiChannelImage3DT.h

@@ -107,6 +107,14 @@ public:
   /** calc integral image */
   void calcIntegral( uint channel = 0 );
   
+  /** 
+   * @brief calculate the gray level co-occurence matrix of a channel
+   * @param max Co-occurence Matrix
+   * @param dis displacement vector
+   * @param channel channel
+   */
+  void calcGLCM( const std::vector<int> dis, uint channel = 0 );
+  
   /**
    * @brief calculate the integral value in the volume given by upper left front corner and lower right back corner, including out of boundary check
    * @warning make sure that the given channel is an integral 3d image

image/MultiChannelImage3DT.tcc

@@ -622,6 +622,73 @@ void MultiChannelImage3DT<P>::calcIntegral( uint channel )
   }
 }
 
+template<class P>
+void MultiChannelImage3DT<P>::calcGLCM( const std::vector<int> dis, uint channel )
+{
+  assert( channel < numChannels );
+  assert( data[channel] != NULL );
+  assert( dis.size() > 2 );
+  
+  P min, max;
+  statistics( min, max, channel );
+  
+  long k = 0;
+  int wsize = 2;    // window size = 2*wsize + 1
+  
+  int tmp[xsize][ysize][zsize];
+  std::vector<std::vector<double> > mat;
+  
+  for (int z = 0; z < zsize; z++ )
+  {
+    for (int y = 0; y < ysize; y++ )
+    {
+      for (int x = 0; x < xsize; x++)
+      {
+        mat = std::vector<std::vector<double> > ( max+1, std::vector<double> ( max+1, 0.0 ) );
+        for (int u=-wsize; u<=wsize; u++)
+        {
+          for (int v=-wsize; v<=wsize; v++)
+          {
+            for (int w=-wsize; w<=wsize; w++)
+            {
+              if ( (u+x+dis[0]>=0) && (u+x+dis[0]<xsize)
+                && (v+y+dis[1]>=0) && (v+y+dis[1]<ysize)
+                && (w+z+dis[2]>=0) && (w+z+dis[2]<zsize) )
+              {
+                P val = get( x+u, y+v, z+w, channel );
+                P dval = get( x+u+dis[0], y+v+dis[1], z+w+dis[2], channel );
+                mat[val][dval]++;
+                mat[dval][val]++;   // symmetry
+                k += 2;
+              }
+            }
+          }
+        }
+        double ASM = 0.0;
+        for (int i = 0; i <= max; i++)
+        {
+          for (int j = 0; j <= max; j++)
+          {
+            mat[i][j] /= k;
+            ASM += mat[i][j]*mat[i][j];
+          }
+        }
+        tmp[x][y][z] = (int) (ASM * get(x,y,z,channel) );
+      }
+    }
+  }
+  for (int z = 0; z < zsize; z++ )
+  {
+    for (int y = 0; y < ysize; y++ )
+    {
+      for (int x = 0; x < xsize; x++)
+      {
+        set(x,y,z, tmp[x][y][z], channel);
+      }
+    }
+  }
+}
+
 template<class P>
 P MultiChannelImage3DT<P>::getIntegralValue(int ulfx, int ulfy, int ulfz, int lrbx, int lrby, int lrbz, int channel)
 {

progs/testSemanticSegmentation.cpp

@@ -153,6 +153,7 @@ int main ( int argc, char **argv )
     NICE::MultiChannelImage3DT<double> probabilities;
     NICE::MultiChannelImage3DT<double> imgData;
     semseg->make3DImage ( filelist, imgData );
+    cout << segresult.width() << " " << segresult.height() << " " << segresult.channels() << endl;
     semseg->semanticseg ( imgData, segresult, probabilities, filelist );
 
     fprintf ( stderr, "testSemanticSegmentation: Segmentation finished !\n" );