Tentative de débuggage du calcul des descripteurs de fourrier

tests/examples/test-descripteurs.cpp

@@ -0,0 +1,38 @@
+#include <math.hpp>
+
+int main(int argc, char** argv) {
+  std::string imagename = "";
+  int seuil = 25;
+  int cmax = 10;
+  if (argc > 2) {
+    imagename = argv[1];
+    seuil = atoi(argv[2]);
+  }
+
+  cv::namedWindow("Image", CV_WINDOW_AUTOSIZE);
+  cv::namedWindow("Binaire", CV_WINDOW_AUTOSIZE);
+  cv::namedWindow("Contour", CV_WINDOW_AUTOSIZE);
+
+  cv::Mat image = cv::imread(imagename, CV_LOAD_IMAGE_COLOR);
+  cv::Mat binaire(image.rows, image.cols, CV_8UC1);
+  cv::Mat contour_image(image.rows, image.cols, CV_8UC1);
+
+  std::vector<std::vector<cv::Point>> contours;
+  std::vector<std::vector<cv::Point>> contrs;
+  std::vector<cv::Vec4i> hierarchy;
+
+  math::filter(image, binaire, seuil);
+  cv::findContours(binaire, contours, hierarchy, CV_RETR_EXTERNAL, CV_CHAIN_APPROX_NONE);
+
+  int index = math::max_cont(contours);
+  contrs.push_back(contours[index]);
+  contrs.push_back(math::simplify_contour(contrs[0], cmax));
+  cv::drawContours(contour_image, contrs, -1, 255);
+
+  imshow("Image", image);
+  imshow("Binaire", binaire);
+  imshow("Contour", contour_image);
+
+  cv::waitKey(0);
+  return 0;
+}

tests/examples/test-fft.cpp

@@ -27,6 +27,8 @@ int main(int argc, char** argv) {
     s.push_back(math::complex(std::sin(2*math::pi()*f0*float(i)/fe), 0));
   }
 
+  std::cout << math::mean(s) << std::endl;
+
   math::csignal tfd;
   clock_t begin = std::clock();
 

tests/src/math.hpp

@@ -13,13 +13,34 @@ namespace math {
   using contour = std::vector<cv::Point>;
   constexpr double pi() {return std::atan(1)*4;}
 
+  int filter(const cv::Mat& img, cv::Mat output, int seuil) {
+  	bool detect = false;
+  	uchar R, G, B;
+  	int rows = img.rows;
+  	int cols = img.cols;
+  	int dim = img.channels();
+  	int indexNB;
+  
+  	for (int index=0,indexNB=0;index<dim*rows*cols;index+=dim,indexNB++) {
+  		detect=0;
+  		B = img.data[index    ];
+  		G = img.data[index + 1];
+  		R = img.data[index + 2];
+  
+  		if ((R>G) && (R>B))
+  			if (((R-B)>=seuil) || ((R-G)>=seuil))
+  				detect=1;
+  		if (detect==1)
+  			output.data[indexNB]=255;
+  		else
+  			output.data[indexNB]=0;
+  	}
+  }
+
   csignal cont2sig(const contour& cont) {
     csignal sig;
-    auto sig_it = sig.begin();
+    for (auto p: cont) {
-    auto cont_it = cont.begin();
+      sig.push_back(complex(p.x, p.y));
-
-    for (auto cont_it = cont.begin(); cont_it != cont.end(); ++cont_it) {
-      *(sig_it++) = complex((*cont_it).x, (*cont_it).y);
     }
     return sig;
   };
@@ -35,7 +56,7 @@ namespace math {
   csignal diff(const csignal& input, complex mean) {
     csignal res;
     for (auto x: input) {
-      res.push_back(x - mean);
+      res.push_back(x-mean);
     }
     return res;
   }
@@ -79,6 +100,8 @@ namespace math {
     int opt_size;
     if (N < input.size()) {
       opt_size = 1 << (int)std::ceil(std::log(input.size())/std::log(2));
+    } else if (N==0){
+      opt_size = input.size();
     } else {
       opt_size = 1 << (int)std::ceil(std::log(N)/std::log(2));
     }
@@ -89,30 +112,85 @@ namespace math {
     return fft_rec(sig);
   };
 
-  contour coef2cont(const csignal& tfd, complex mean, int size, int cmax) {
+  void operator*=(csignal& sig, complex& m) {
-    contour cont;
+    for(auto x: sig) {
-    auto tf_it = tfd.begin();
+      x *= m;
-    auto cont_it = cont.begin();
+    }
-    int kmin = tfd.size()/2 - cmax;
+  }
-    int kmax = tfd.size()/2 + cmax;
 
-    for (int m=0; m<tfd.size(); ++m) {
+  void operator*=(csignal& sig, complex&& m) {
+    for(auto x: sig) {
+      x *= m;
+    }
+  }
+
+  void operator/=(csignal& sig, complex& m) {
+    for(auto x: sig) {
+      x /= m;
+    }
+  }
+
+  void operator/=(csignal& sig, complex&& m) {
+    for(auto x: sig) {
+      x /= m;
+    }
+  }
+
+  csignal extract(const csignal& tfd, int cmax) {
+    csignal res;
+    for (int k=0; k<cmax; ++k) {
+      res.push_back(tfd[tfd.size() - cmax + k]);
+    }
+    for (int k=cmax; k<2*cmax; ++k) {
+      res.push_back(tfd[k]);
+    }
+    return res;
+  }
+
+  contour sig2cont(const csignal& sig) {
+    contour res;
+    for (auto x: sig) {
+      res.push_back(cv::Point(x.real(), x.imag()));
+    }
+    return res;
+  }
+
+  csignal desc2sig(const csignal& desc, complex mean, int N, int kmin) {
+    csignal cont;
+    auto desc_it = desc.begin();
+
+    for (int m=0; m<N; ++m) {
       complex sum;
-      for (int k=kmin; k<kmax; ++k) {
+      for (int k=0; k<desc.size(); ++k) {
-        sum += tfd[k]*std::exp(complex(0, 2*pi()*k*m/tfd.size()));
+        sum += desc[k]*std::exp(complex(0, 2*pi()*(k+kmin)*m/N));
       }
-      complex zm = mean + sum;
+      cont.push_back(mean + sum);
-      *(cont_it++) = cv::Point(zm.real(), zm.imag());
     }
     return cont;
   };
 
   contour simplify_contour(const contour& cont, int cmax) {
-    contour res;
     csignal z = cont2sig(cont);
     complex zm = mean(z);
     csignal tfd = fft(diff(z, zm));
-    return coef2cont(tfd, zm, 0, cmax);
+    tfd /= tfd.size();
+    csignal desc = extract(tfd, cmax);
+
+    if (std::abs(desc[desc.size()-1]) > std::abs(desc[0])) {
+      std::reverse(desc.begin(), desc.end());
+    }
+
+    double phy = std::arg(desc[desc.size()-1]*desc[0])/2;
+    desc *= std::exp(complex(0, -phy));
+    double theta = std::arg(desc[0]);
+
+    for (int k=0; k<desc.size(); ++k) {
+      desc[k] *= std::exp(complex(0, -theta*k));
+    }
+    desc /= desc[0];
+
+    csignal sig = desc2sig(desc, zm, z.size(), cmax);
+    return sig2cont(sig);
   };
 
   int max_cont(const std::vector<contour>& contours) {

tests/src/traitement.cpp

@@ -28,15 +28,14 @@ int filter(const cv::Mat& img, cv::Mat output, int seuil) {
 }
 
 int main(int argc, char** argv) {
-	int seuil=0;
+	int seuil=80;
+	int cmax = 10;
+	int N = 5000;
 
 	if (argc > 1) {
 		seuil = atol(argv[1]);
 	}
 
-	int cmax = 10;
-	int N = 5000;
-
 	char detect;
 	cv::VideoCapture cap(0); 
 	if(!cap.isOpened())  
@@ -44,6 +43,7 @@ int main(int argc, char** argv) {
 	cv::namedWindow("Image",1);
 	cv::namedWindow("Detection",1);
 	cv::namedWindow("Contours",1);
+	cv::namedWindow("New Contours",1);
 	while(true) {
 		int X, Y, DIM, index;
 		unsigned int numc;
@@ -66,15 +66,27 @@ int main(int argc, char** argv) {
 		cv::imshow("Detection", binaire);
 		cv::findContours(binaire, contours, hierarchy, CV_RETR_EXTERNAL, CV_CHAIN_APPROX_NONE);
 		cv::Mat Dessin = cv::Mat::zeros(X,Y, CV_8UC1);
+		cv::Mat new_contour_image = cv::Mat::zeros(X,Y, CV_8UC1);
 
-    int id = math::max_cont(contours);
+    if (contours.size() > 0) {
-		std::cout << contours[id].size() << std::endl;
+		  std::vector<std::vector<cv::Point>> contrs;
 
-		std::vector<cv::Point> new_cont = math::simplify_contour(contours[id], cmax);
+      int id = math::max_cont(contours);
+      contrs.push_back(contours[id]);
 
-		cv::drawContours(Dessin, contours, id, 255);
+      std::cout << "Number of countours: "
-		cv::drawContours(Dessin, new_cont, id, 255);
+                << contours.size()
+                << "; Index of biggest contour: "
+                << id
+                << std::endl;
+
+      contrs.push_back(math::simplify_contour(contrs[0], cmax));
+
+      cv::drawContours(Dessin, contrs, 0, 255);
+      cv::drawContours(new_contour_image, contrs, 1, 255);
+    }
 		cv::imshow("Contours", Dessin);
+		cv::imshow("New Contours", new_contour_image);
 		if(cv::waitKey(30) == 27) {
 			break;
 		}