Figs bugs, add window option to emvisi2_video

- fix a memory issue (allocation too small)
- fix some numerical instability

Author: jpilet

emvisi2.cpp

@@ -105,7 +105,9 @@ void EMVisi2::run(int nbIter, float smooth_amount, float smooth_threshold)
 		reset_gaussians();
 	}
 	for (int i=0;i<nbIter;i++) {
-		iterate();
+		if (!iterate()) {
+                  break;
+                }
 	}
 	if (smooth_amount>0)
 		smooth(smooth_amount, smooth_threshold);
@@ -127,7 +129,7 @@ void EMVisi2::reset_gaussians()
 		weights[i] = 1.0f/(NB_GAUSSIANS+1);
 }
 
-void EMVisi2::iterate()
+bool EMVisi2::iterate()
 {
 	char str[256];
 	uniform_resp=0;
@@ -171,12 +173,20 @@ void EMVisi2::iterate()
 	for (int i=0; i<NB_VISI_GAUSSIANS; i++) {
 		visi_g[i].compute_sigma();
 		totN += visi_g[i].n;
+                if (visi_g[i].n < 1e-2f) {
+                  recycle = false;
+                  return false;
+                }
 	}
 
 	for (int i=0; i<NB_OCCL_GAUSSIANS; i++) {
 		occl_g[i].compute_sigma();
 		totN += occl_g[i].n;
-	}
+                if (occl_g[i].n < 1e-2f) {
+                  recycle = false;
+                  return false;
+                }
+        }
 
 	totN += uniform_resp;
 	weights[NB_GAUSSIANS] = uniform_resp/totN;
@@ -199,6 +209,7 @@ void EMVisi2::iterate()
 		sprintf(str, "proba%02d.png", iteration );
 		scale_save(str, proba, -1, -1);
 	}
+        return true;
 }
 
 float EMVisi2::process_pixel(const float *rgb, const float *frgb, const float dl, const float nccv, const float ncch, float *proba, float *visi_proba)
@@ -511,7 +522,7 @@ void EMVisi2::smooth(float amount, float threshold) {
 	}
 
 	// allocate the graph. Note: worst case memory scenario.
-	int n_nodes= gc_mask.cols * gc_mask.rows;
+	int n_nodes= (gc_mask.cols + 2) * (gc_mask.rows + 2);
 	int n_edges = 2*((wa.cols)*(wa.rows-1) + (wa.cols-1)*wa.rows);
 	FGraph g(n_nodes, n_edges);
 	int *nodesid = new int[n_nodes];
@@ -546,6 +557,8 @@ void EMVisi2::smooth(float amount, float threshold) {
 				if (m[x] == rval) {
 					// add a new node
                                         int next_node = g.addVtx();
+                                        assert(row_id >= nodesid);
+                                        assert(row_id < nodesid + n_nodes);
 					*row_id = next_node;
 
 					// terminal weights
@@ -601,7 +614,7 @@ void EMVisi2::smooth(float amount, float threshold) {
 
 			for (int x=rect.x; x<rect.x+rect.width; x++) {
 				if (*row_id >= 0) {
-					p[x] = !g.inSourceSegment(*row_id);
+					p[x] = (g.inSourceSegment(*row_id) ? 0 : 1);
 				}
 				row_id++;
 			}

emvisi2.h

@@ -90,7 +90,7 @@ class EMVisi2 {
 	bool init();
         int setModel(const cv::Mat im1, const cv::Mat mask = cv::Mat());
 	int setTarget(cv::Mat target);
-	void iterate();
+	bool iterate();
 	void smooth(float amount, float threshold);
 	void run(int nbIter=3, float smooth_amount=2.5, float smooth_threshold=.0001);
 	float process_pixel(const float *rgb, const float *frgb, const float dl, const float nccv, const float ncch, float *proba, float *visi_proba);

emvisi2_video.cpp

@@ -61,6 +61,7 @@ int main(int argc, char *argv[])
   cv::Mat background;
   list<string> images;
   path destinationFolder("out");
+  bool window = false;
 
   // parse command line
   for (int narg=1; narg<argc; ++narg) {
@@ -85,6 +86,8 @@ int main(int argc, char *argv[])
 
     if (strcmp(argv[narg],"-v")==0) {
       emv.save_images=true;
+    } else if (strcmp(argv[narg],"-w")==0) {
+      window = true;
     } else if (strcmp(argv[narg],"-s")==0) {
       smooth = 0.001f;
     } else if (argv[narg][0] == '-') {
@@ -111,6 +114,8 @@ int main(int argc, char *argv[])
 
   cout << "done in " << timer.duration() << " ms.\n";
 
+  int playing = 0;
+
   for (list<string>::const_iterator it = images.begin();
        it != images.end(); ++it) {
     Mat frame = loadImage(*it, -1);
@@ -128,6 +133,21 @@ int main(int argc, char *argv[])
     cout << "computed in " << timer.duration() << " ms.\n";
 
     save_proba(destinationFilename(*it, destinationFolder).c_str(),emv.proba);
+
+    if (window) {
+      std::vector<cv::Mat> channels;
+      cv::split(frame, channels);
+      emv.proba.convertTo(channels[1], CV_8UC1, 255, 0);
+      cv::Mat result;
+      cv::merge(channels, result);
+
+      imshow("Segmentation", result);
+      int k = cv::waitKey(playing);
+      switch (k) {
+        case ' ': playing = (playing + 1) & 1; break;
+        case 'p': --it; --it; break;
+      }
+    }
   }
   return 0;
 }

imstat.h

@@ -232,14 +232,19 @@ class MultiGaussian3 : public Accumulator<T> {
 			+t14*sigma[2][2]
 			+t15+t15;
 
+		float epsilon = 1e-30;
+
+                if (isnan(vsv)) {
+                  return epsilon;
+                }
+
 		assert(!isnan(vsv));
 		assert(!isnan(one_over_sq_det));
 		if (vsv<0) vsv=0;
 
 		// this constant is 1/(2Pi)^(3/2)
 		float p = .0634936359342409f*one_over_sq_det*exp_table(-.5f*vsv);
 
-		float epsilon = 1e-30;
 		if (p<epsilon) p=epsilon;
 
 		assert(!isnan(p));