stable version - todo: check if given img has only 1 dimension (gray scale)

segment-graph.h

@@ -48,7 +48,18 @@ bool operator<(const edge &a, const edge &b) {
 universe *segment_graph(int num_vertices, int num_edges, edge *edges, 
 			float c) { 
   // sort edges by weight
-  std::sort(edges, edges + num_edges);
+  // note: if weights occure more then once, this might lead to non-deterministic (non-reproducable) behaviour, since
+  // "Equivalent elements are not guaranteed to keep their original relative order"
+//   std::sort(edges, edges + num_edges);
+  
+/*   for (int i = 0; i < 5; i++) {
+      int a = edges[i].a;
+      int b = edges[i].b;
+      mexPrintf("edges before %d: %d <-> %d with weight %f\n",i, a, b, edges[i].w);      
+    } */  
+  
+  // adaptation: use stable_sort instead, which will keep the relative position of equal elements :)
+  std::stable_sort(edges, edges + num_edges);
 
   // make a disjoint-set forest
   universe *u = new universe(num_vertices);
@@ -57,6 +68,20 @@ universe *segment_graph(int num_vertices, int num_edges, edge *edges,
   float *threshold = new float[num_vertices];
   for (int i = 0; i < num_vertices; i++)
     threshold[i] = THRESHOLD(1,c);
+  
+    //that's the same as with find, and also non-reproducable
+//    for (int i = 0; i < 5; i++) {
+//       int a = edges[i].a;
+//       int b = edges[i].b;
+//       mexPrintf("edges after %d: %d <-> %d with weight %f\n",i, a, b, edges[i].w);    
+//     } 
+  
+//      for (int i = 0; i < 5; i++) {
+//       int a = u->find(edges[i].a);
+//       int b = u->find(edges[i].b);
+//       mexPrintf("components %d: %d <-> %d\n",i, a, b);      
+//     } 
+  // this is already non-reproducable
 
   // for each edge, in non-decreasing weight order...
   for (int i = 0; i < num_edges; i++) {
@@ -74,6 +99,7 @@ universe *segment_graph(int num_vertices, int num_edges, edge *edges,
       }
     }
   }
+    
 
   // free up
   delete threshold;

segment-image-labelOutput.h

@@ -2,12 +2,16 @@
 #define SEGMENT_IMAGE_LABELOUTPUT
 
 #include <cstdlib>
+#include <map>
+#include <limits>
 #include "./image.h"
 #include "./misc.h"
 #include "./filter.h"
 #include "./segment-graph.h"
 #include "./segment-image.h"
 
+#include <mex.h>
+
 /*
  * Segment an image
  *
@@ -49,6 +53,10 @@ image<int> *segment_image_labelOutput(image<rgb> *im, float sigma, float c, int
   for (int y = 0; y < height; y++) {
     for (int x = 0; x < width; x++) {
       if (x < width-1) {
+// 	if ( num == 0)
+// 	{
+// 	  mexPrintf("r: %f , b: %f g: %f\n", imRef(smooth_r, x, y), imRef(smooth_g, x, y), imRef(smooth_b, x, y));
+// 	}
 	edges[num].a = y * width + x;
 	edges[num].b = y * width + (x+1);
 	edges[num].w = diff(smooth_r, smooth_g, smooth_b, x, y, x+1, y);
@@ -80,11 +88,27 @@ image<int> *segment_image_labelOutput(image<rgb> *im, float sigma, float c, int
   delete smooth_r;
   delete smooth_g;
   delete smooth_b;
+  
+//      for (int i = 0; i < 5; i++) {
+//       int a = edges[i].a;
+//       int b = edges[i].b;
+//       mexPrintf("components %d: %d <-> %d\n",i, a, b);      
+//     }   
+    // so far, everything is always reproducable
 
   // segment
   universe *u = segment_graph(width*height, num, edges, c);
   
-  // post process small components
+//   // post process small components
+//   mexPrintf("components before PP: %d\n",num);  
+  
+
+//     for (int i = 0; i < 10; i++) {
+//       int a = u->find(edges[i].a);
+//       int b = u->find(edges[i].b);
+//       mexPrintf("components %d: %d <-> %d\n",i, a, b);      
+//     }
+  
   for (int i = 0; i < num; i++) {
     int a = u->find(edges[i].a);
     int b = u->find(edges[i].b);
@@ -92,14 +116,39 @@ image<int> *segment_image_labelOutput(image<rgb> *im, float sigma, float c, int
       u->join(a, b);
   }
   delete [] edges;
+  
+  //NOTE since u->num_sets() is an int, we are not able to handle more than 
+  // max_int segments ( and therefore also not imgs with more pixels)
   *num_ccs = u->num_sets();
+  
+//   mexPrintf("components after PP: %d\n",*num_ccs); 
 
   image<int> *output = new image<int>(width, height);
+ 
+  //how many different regions do we finally have, and which of them are the corresponding indices?
+  std::map<int,int> regionLabels;
+  int idx ( 0 );
   
-  for (int y = 0; y < height; y++) {
-    for (int x = 0; x < width; x++) {
+  for (int y = 0; y < height; y++)
+  {
+    for (int x = 0; x < width; x++)
+    {
+      int comp = u->find(y * width + x);
+      if ( regionLabels.find( comp ) == regionLabels.end() )
+      {
+	regionLabels.insert( std::pair<int,int>(comp,idx) );
+	idx++;
+      }
+    }
+  }  
+  
+  // insert region indices into img, map indices from current ones to consecutive ones
+  for (int y = 0; y < height; y++)
+  {
+    for (int x = 0; x < width; x++)
+    {
       int comp = u->find(y * width + x);
-      imRef(output, x, y) = comp;
+      imRef(output, x, y) = regionLabels.find(comp)->second;
     }
   }  
  

segmentFelzenszwalb.cpp

@@ -10,10 +10,18 @@
 #include "./segment-image.h"
 #include "./segment-image-labelOutput.h"
 
+#ifndef MIN
+#define MIN(x, y) (((x) < (y)) ? (x) : (y))	/* min and max value
+macros */
+#endif
+#ifndef MAX
+#define MAX(x, y) (((x) > (y)) ? (x) : (y))
+#endif
+
 void mexFunction(int nOutput, mxArray *pOutput[], /* Output variables */
   int nInput, const mxArray *pInput[]) /* Input variables */
 { 
-    
+   
    /* -----------------------------------------------------------------
    *                         Check the arguments
    * -------------------------------------------------------------- */
@@ -172,13 +180,17 @@ void mexFunction(int nOutput, mxArray *pOutput[], /* Output variables */
    image<int> * imgResult  = NULL;
    
    if ( computeColorOutput )
+   {
      imgResultRGB = segment_image(imgInput, sigma, k, minSize, &num_ccs); 
+   }
    else
+   {
     imgResult = segment_image_labelOutput(imgInput, sigma, k, minSize, &num_ccs); 
+   }
        
    if ( verbose )
      mexPrintf("segmentation done\n");   
-  
+   
    
    /* -----------------------------------------------------------------
    *                      Conversion to Matlab structures
@@ -186,11 +198,19 @@ void mexFunction(int nOutput, mxArray *pOutput[], /* Output variables */
    
    if ( nOutput == 0 )
    {
-     if ( computeColorOutput )
-       savePPM( imgResultRGB, output );
-    else
-      save_image( imgResult, output );
-     mexPrintf("save results\n"); 
+     //was a filename given to store the image in?
+     if (nInput >=  6) 
+     {
+      if ( computeColorOutput )
+	savePPM( imgResultRGB, output );
+      else
+	save_image( imgResult, output );
+      mexPrintf("save results\n"); 
+     }
+     else
+     {
+       mexPrintf("ATTENTION -- Neither output variables nor filename to store the result to given!");
+     }
    }
    else
    {
@@ -205,9 +225,13 @@ void mexFunction(int nOutput, mxArray *pOutput[], /* Output variables */
 	int height (imgResultRGB->height() );       
 	/* keep in mind that matlab stores images  height × width × color, whereas C does it the other way round*/
 	int dims[] = {height, width,3};
-	pOutput[0] = mxCreateNumericArray (3, dims, mxUINT8_CLASS, mxREAL);
-	unsigned char *out1; /* pointer to output 1 */
-	out1 = (unsigned char *)mxGetPr( pOutput[0] ); /* pointer to output 1 */
+	pOutput[0] = mxCreateNumericArray (3, dims, mxINT32_CLASS, mxREAL);
+	//unsigned char *out1; /* pointer to output 1 */
+	//out1 = (unsigned char *)mxGetPr( pOutput[0] ); /* pointer to output 1 */
+	//unsigned char *out1; /* pointer to output 1 */
+	//out1 = (unsigned char *)mxGetPr( pOutput[0] ); /* pointer to output 1 */
+        int *out1; /* pointer to output 1 */
+        out1 = (int *)mxGetData( pOutput[0] ); /* pointer to output 1 *	
 
 	      
 	/* start with RED channel*/     
@@ -217,7 +241,7 @@ void mexFunction(int nOutput, mxArray *pOutput[], /* Output variables */
 	  uint rowOffset ( x*height );
 	  for ( uint y = 0; y < height; y++)
 	  {
-	      out1[rowOffset + y ] = (double) (imRef(imgResultRGB, x, y)).r;
+	      out1[rowOffset + y ] = (int) (imRef(imgResultRGB, x, y)).r;
 	  }
 	}
 	
@@ -228,7 +252,7 @@ void mexFunction(int nOutput, mxArray *pOutput[], /* Output variables */
 	  uint rowOffset ( x*height );
 	    for ( uint y = 0; y < height; y++)
 	    {
-	      out1[channelOffsetG + rowOffset + y ] = (double) (imRef(imgResultRGB, x, y)).g;
+	      out1[channelOffsetG + rowOffset + y ] = (int) (imRef(imgResultRGB, x, y)).g;
 	    }
 	}
 	
@@ -239,7 +263,7 @@ void mexFunction(int nOutput, mxArray *pOutput[], /* Output variables */
 	  uint rowOffset ( x*height );
 	    for ( uint y = 0; y < height; y++)
 	    {
-	      out1[channelOffsetB + rowOffset + y ] = (double) (imRef(imgResultRGB, x, y)).b;
+	      out1[channelOffsetB + rowOffset + y ] = (int) (imRef(imgResultRGB, x, y)).b;
 	    }
 	}
      }
@@ -248,29 +272,68 @@ void mexFunction(int nOutput, mxArray *pOutput[], /* Output variables */
 	int width ( imgResult->width() );
 	int height (imgResult->height() ); 
         int dims[] = {height, width};
-       	pOutput[0] = mxCreateNumericArray (2, dims, mxUINT8_CLASS, mxREAL);
-	unsigned char *out1; /* pointer to output 1 */
-	out1 = (unsigned char *)mxGetPr( pOutput[0] ); /* pointer to output 1 */
 	
-       	/* keep in mind that matlab stores images col by width, whereas C does it the other way round*/
-	for ( uint x = 0; x < width; x++)
+	if (num_ccs < 255)
 	{
-	  uint rowOffset ( x*height );
-	  for ( uint y = 0; y < height; y++)
+	  pOutput[0] = mxCreateNumericArray (2, dims, mxUINT8_CLASS, mxREAL);
+	  unsigned char *out1; /* pointer to output 1 */
+	  out1 = (unsigned char *)mxGetPr( pOutput[0] ); /* pointer to output 1 */
+	  
+	  /* keep in mind that matlab stores images col by width, whereas C does it the other way round*/
+	  for ( uint x = 0; x < width; x++)
+	  {
+	    uint rowOffset ( x*height );
+	    for ( uint y = 0; y < height; y++)
+	    {
+		out1[rowOffset + y ] = (unsigned char ) (imRef(imgResult, x, y));
+	    }
+	  }	  
+	}
+	else
+	{
+	  pOutput[0] = mxCreateNumericArray (2, dims, mxINT32_CLASS, mxREAL);
+	  int *out1; /* pointer to output 1 */
+	  out1 = (int *)mxGetData( pOutput[0] ); /* pointer to output 1 *
+	  
+	  /* keep in mind that matlab stores images col by width, whereas C does it the other way round*/
+	  for ( uint x = 0; x < width; x++)
 	  {
-	      out1[rowOffset + y ] = (uint) (imRef(imgResult, x, y));
+	    uint rowOffset ( x*height );
+	    for ( uint y = 0; y < height; y++)
+	    {
+		out1[rowOffset + y ] = (int) (imRef(imgResult, x, y));
+	    }
 	  }
 	}
-     }
+    }
+     
      
      /* return number of segments*/
      if ( nOutput >= 2 )
      {
        int dims[] = {1};
-       pOutput[1] = mxCreateNumericArray(1,dims,mxINT8_CLASS,mxREAL);
-       unsigned char *out2; /* pointer to output 2 */
-       out2 = (unsigned char *)mxGetPr( pOutput[1] ); /* pointer to output 2 */      
-       out2[0] = num_ccs;
+       if (num_ccs < 255)
+       {
+	pOutput[1] = mxCreateNumericArray(1,dims,mxUINT8_CLASS,mxREAL);
+	unsigned char *out2; /* pointer to output 2 */
+	out2 = (unsigned char *)mxGetPr( pOutput[1] ); /* pointer to output 2 */
+	
+	out2[0] = num_ccs;
+       }
+       else
+       {
+	pOutput[1] = mxCreateNumericArray(1,dims,mxINT32_CLASS,mxREAL);
+	//unsigned char *out2; /* pointer to output 2 */
+	//out2 = (unsigned char *)mxGetPr( pOutput[1] ); /* pointer to output 2 */
+	int *out2; /* pointer to output 2 */
+	out2 = (int *) mxGetData( pOutput[1] ); /* pointer to output 2 */
+	
+	out2[0] = num_ccs;
+
+       }
+       
+       if ( verbose ) 
+	mexPrintf("number of components: %i", num_ccs);
      }
 
     /* done */