[Olena-patches] proto-1.0 451: Add a closed snake algorithm.

5 Apr 2006

https://svn/svn/oln/prototypes/proto-1.0/olena

Index: ChangeLog
from  Nicolas Widynski  <nicolas.widynski@lrde.epita.fr>

	Add a closed snake algorithm.

	* oln/appli/snakes: New. Directory.
	* oln/appli/snakes/snakes.hh: New. Snake algorithm.
				      2 forces : gradient and snake perimeter.

 snakes.hh |  413 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 413 insertions(+)

Index: oln/appli/snakes/snakes.hh

--- oln/appli/snakes/snakes.hh	(revision 0)
+++ oln/appli/snakes/snakes.hh	(revision 0)
@@ -0,0 +1,413 @@
+// Copyright (C) 2001-2006  EPITA Research and Development Laboratory
+//
+// This file is part of the Olena Library.  This library is free
+// software; you can redistribute it and/or modify it under the terms
+// of the GNU General Public License version 2 as published by the
+// Free Software Foundation.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+// General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this library; see the file COPYING.  If not, write to
+// the Free Software Foundation, 59 Temple Place - Suite 330, Boston,
+// MA 02111-1307, USA.
+//
+// As a special exception, you may use this file as part of a free
+// software library without restriction.  Specifically, if other files
+// instantiate templates or use macros or inline functions from this
+// file, or you compile this file and link it with other files to
+// produce an executable, this file does not by itself cause the
+// resulting executable to be covered by the GNU General Public
+// License.  This exception does not however invalidate any other
+// reasons why the executable file might be covered by the GNU General
+// Public License.
+
+#ifndef OLENA_SNAKES_HH
+# define OLENA_SNAKES_HH
+
+# include <mlc/is_a.hh>
+# include <oln/basics2d.hh>
+# include <ntg/all.hh>
+# include <oln/morpho/gradient_morpho.hh>
+# include <ntg/int.hh>
+# include <vector>
+
+#define sqr(u) ((u)*(u))
+const float PI = 3.14159265;
+
+
+namespace oln {
+
+  namespace snakes {
+
+    namespace impl {
+
+      inline
+      float euclidian_distance(int x1, int x2, int y1, int y2) // sqr
+      {
+	return sqr(x1 - x2) + sqr(y1 - y2);
+      }
+
+      template <typename T>
+      void verify_integrity_(const image2d<T>& grad,
+			     int &x1,
+			     int &x2,
+			     int &y1,
+			     int &y2)
+      {
+	if (x1 < 0)
+	  x1 = 0;
+	if (x2 < 0)
+	  x2 = 0;
+	if (x1 > grad.nrows())
+	  x1 = grad.nrows() - 1;
+	if (x2 > grad.nrows())
+	  x2 = grad.nrows() - 1;
+
+	if (y1 < 0)
+	  y1 = 0;
+	if (y2 < 0)
+	  y2 = 0;
+	if (y1 > grad.ncols())
+	  y1 = grad.ncols() - 1;
+	if (y2 > grad.ncols())
+	  y2 = grad.ncols() - 1;
+      }
+
+
+      std::vector<point2d> draw_line_(const point2d& p,
+				      const point2d& q)
+      {
+	std::vector<point2d> res;
+
+	point2d p_cur = p;
+	if (q.col() == p.col()) // cas vertical
+	  {
+	    int sens = q.row() > p.row() ? 1 : -1;
+	    while (p_cur != q)
+	      {
+		double pc_p = euclidian_distance(p_cur.row(), p.row(), p_cur.col(), p.col());
+		double pc_q = euclidian_distance(p_cur.row(), q.row(), p_cur.col(), q.col());
+		double p_q = euclidian_distance(p.row(), q.row(), p.col(), q.col());
+
+		res.push_back(p_cur);
+		p_cur.row() = p_cur.row() + sens;
+	      }
+	  }
+	else
+	  {
+	    if (q.row() == p.row()) // cas horizontal
+	      {
+		int sens = q.col() > p.col() ? 1 : -1;
+		while (p_cur != q)
+		  {
+		    double pc_p = euclidian_distance(p_cur.row(), p.row(), p_cur.col(), p.col());
+		    double pc_q = euclidian_distance(p_cur.row(), q.row(), p_cur.col(), q.col());
+		    double p_q = euclidian_distance(p.row(), q.row(), p.col(), q.col());
+
+		    res.push_back(p_cur);
+		    p_cur.col() = p_cur.col() + sens;
+		  }
+	      }
+	    else // cas diagonal
+	      {
+		long decalage = 0;
+
+		long nom = q.row() - p.row();
+		long denom = q.col() - p.col();
+		int signe = (nom * denom) > 0 ? 1 : -1;
+		nom = nom * signe;
+		denom = denom * signe;
+
+		int vsens = q.row() > p.row() ? 1 : -1;
+		int hsens = q.col() > p.col() ? 1 : -1;
+
+		int deboule = 0;
+
+		while (p_cur != q)
+		  {
+		    double pc_p = euclidian_distance(p_cur.row(), p.row(), p_cur.col(), p.col());
+		    double pc_q = euclidian_distance(p_cur.row(), q.row(), p_cur.col(), q.col());
+		    double p_q = euclidian_distance(p.row(), q.row(), p.col(), q.col());
+
+		    res.push_back(p_cur);
+
+		    if (abs(decalage) >= abs(denom))
+		      {
+			p_cur.row() = p_cur.row() + vsens;
+			decalage -= vsens * denom;
+		      }
+		    else
+		      {
+			p_cur.col() = p_cur.col() + hsens;
+			decalage += nom * hsens;
+		      }
+		  }
+	      }
+	  }
+	return res;
+      }
+
+
+
+
+      template <typename T>
+      std::vector<point2d> compute_normal_points_(const image2d<T>& grad,
+						  const point2d& p1,
+						  const point2d& p2,
+						  const point2d& p3,
+						  unsigned int dmax)
+      {
+	float x_c = p2.row();
+	float y_c = p2.col();
+	float a_p1_p3 = (p1.col() - p3.col()) != 0 ? (float)(((float)(p1.row() - p3.row())) / (p1.col() - p3.col())) : 0;
+	float a = -a_p1_p3;
+	float b = p2.col() - a * p2.row(); // use of - slope of p1 et p3
+
+	int x1 = 0;
+	int y1 = 0;
+	int x2 = 0;
+	int y2 = 0;
+
+	if (a != 0)
+	  {
+	    float A = 1 + sqr(a);
+	    float B = -2 * (x_c + a * y_c - a * b);
+	    float C = sqr(x_c) + sqr((b - y_c)) - sqr(dmax);
+
+	    float delta = sqr(B) - (4 * A * C);
+
+	    if (delta > 0)
+	      {
+		x1 = (int)(roundf((-B + sqrt(delta)) / (2 * A)));
+		x2 = (int)(roundf((-B - sqrt(delta)) / (2 * A)));
+
+		if (x1 == x2)
+		  {
+		    float tmp1 = roundf((-B + sqrt(delta)) / (2 * A));
+		    float tmp2 = roundf((-B - sqrt(delta)) / (2 * A));
+		    if (tmp1 > tmp2)
+		      x1++;
+		    else
+		      x2++;
+		  }
+
+		y1 = (int)(a * x1 + b);
+		y2 = (int)(a * x2 + b);
+	      }
+	  }
+	else
+	  {
+	    x1 = p2.row() - dmax;
+	    y1 = p2.col();
+	    x2 = p2.row() + dmax;
+	    y2 = p2.col();
+	  }
+	verify_integrity_(grad, x1, x2, y1, y2);
+
+	return draw_line_(point2d(x1, y1), point2d(x2, y2));
+      }
+
+
+      float angle_force(const point2d& p1, const point2d& p2, const point2d& p3)
+      { // al kashi
+	float b = euclidian_distance(p1.row(), p2.row(), p1.col(), p2.col());
+	float c = euclidian_distance(p3.row(), p2.row(), p3.col(), p2.col());
+	float a = euclidian_distance(p1.row(), p3.row(), p1.col(), p3.col());
+	float angle =  acos((b + c - a) / (2 * sqrt(b) * sqrt(c))); // rad
+	angle = (angle * 180) / PI; //deg
+	return angle;
+      }
+
+
+
+      template <typename T>
+      std::vector<point2d> gen_snakes_(const image2d<T>& grad,
+				       const std::vector<point2d>& v,
+				       unsigned int fen,
+				       float lambda,
+				       unsigned int dmax)
+      {
+	std::vector<point2d>::const_iterator it = v.begin();
+	std::vector<point2d> res;
+
+	for (; (it + 2) != v.end(); it++)
+	  {
+	    point2d p1 = *it;
+	    point2d p2 = *(it + 1);
+	    point2d p3 = *(it + 2);
+	    std::vector<point2d> v_normal = compute_normal_points_(grad, p1, p2, p3, dmax);
+	    std::vector<point2d>::const_iterator it2 = v_normal.begin();
+
+	    float max = -1;
+	    point2d n_pt(0,0);
+
+	    for (; it2 != v_normal.end(); it2++)
+	      {
+		float force = grad[*it2];
+		if (max < force)
+		  {
+		    max = force;
+		    n_pt = *it2;
+		  }
+	      }
+	    res.push_back(n_pt);
+	  }
+
+	res.push_back(res[0]);
+	res.push_back(res[1]);
+
+
+	if (lambda != 0)
+	  {
+	    bool conv = false;
+
+	    while (conv == false)
+	      {
+		conv = true;
+		float max = 0;
+		std::vector<point2d>::iterator it3 = res.begin();
+		std::vector<point2d>::iterator p_s;
+		point2d n_pt(0,0);
+		for (; (it3 + 2) != res.end(); it3++)
+		  {
+		    point2d p1 = *it3;
+		    point2d p2 = *(it3 + 1);
+		    point2d p3 = *(it3 + 2);
+
+		    std::vector<point2d> v_normal = compute_normal_points_(grad, p1, p2, p3, dmax);
+		    std::vector<point2d>::const_iterator it2 = v_normal.begin();
+
+		    float angle_ref = angle_force(p1, p2, p3);
+
+		    float force_ref = grad[p2] - (180 - angle_ref) * lambda;
+
+		    for (; it2 != v_normal.end(); it2++)
+		      {
+			float angle = angle_force(p1, *it2, p3);
+			float force = grad[*it2] - (180 - angle) * lambda;
+
+			if (p2.col() != (*it2).col() and p2.row() != (*it2).row() and force - force_ref > max)
+			  {
+			    max = force - force_ref;
+			    p_s = it3 + 1;
+			    n_pt = *it2;
+			  }
+		      }
+		  }
+
+		if (max > 0)
+		  {
+		    conv = false;
+		    if (p_s == (res.begin() + 1)) // 2nd
+		      {
+			(*(res.end() - 1)).col() = n_pt.col();
+			(*(res.end() - 1)).row() = n_pt.row();
+		      }
+		    else
+		      if (p_s == (res.end() - 2)) // last
+			{
+			  (*(res.begin())).col() = n_pt.col();
+			  (*(res.begin())).row() = n_pt.row();
+			}
+		    (*p_s).col() = n_pt.col();
+		    (*p_s).row() = n_pt.row();
+		  }
+	      }
+	  }
+
+	return res;
+      }
+
+
+      std::vector<point2d> regen_snaxels(std::vector<point2d>& v)
+      {
+	std::vector<point2d> res;
+	std::vector<point2d> tmp;
+	std::vector<point2d>::iterator it = v.begin();
+	std::vector<point2d>::iterator it2;
+
+	for (it = v.begin(); (it + 1) != v.end(); it++)
+	  {
+	    std::vector<point2d> tmp2 = draw_line_(*it, *(it+1));
+	    for (it2 = tmp2.begin(); it2 != tmp2.end(); it2++)
+	      tmp.push_back(*it2);
+	  }
+
+	int nb_pts = v.size() * 2;
+	int cpt = 1;
+
+	for (it2 = tmp.begin(); it2 != tmp.end(); it2++, cpt++)
+	  if ((cpt % (tmp.size() / nb_pts)) == 0)
+	    res.push_back(*it2);
+
+	return res;
+      }
+
+      template <typename T>
+      void clean_ima(image2d<T>& ima)
+      {
+	for (int j = 0; j < ima.size().ncols(); j++)
+	  {
+	    ima(-1, j) = ima(0, j);
+	    ima(ima.size().nrows(), j) = ima(ima.size().nrows() - 1, j);
+	  }
+
+	for (int j = 0; j < ima.size().nrows(); j++)
+	  {
+	    ima(j, -1) = ima(j, 0);
+	    ima(j, ima.size().ncols()) = ima(j, ima.size().ncols() - 1);
+	  }
+      }
+
+
+      template <typename T>
+      std::vector<point2d>
+      snakes_(const image2d<T>& input,
+	      unsigned int fen,
+	      unsigned int nb_gen,
+	      float lambda,
+	      unsigned int dmax)
+      {
+	image2d<T> grad = morpho::gradient_morpho(input, win_c4p());
+	clean_ima(grad);
+
+	std::vector<point2d> res;
+
+	// build init points + 2 (close snake)
+
+	int i = 0;
+	while (i < nb_gen)
+	  {
+	    res = gen_snakes_(grad, res, fen, lambda, dmax);
+	    i++;
+	    if (i < nb_gen)
+	      res = regen_snaxels(res);
+	  }
+
+	return res;
+      }
+
+
+    } // end of namespace oln::impl::snakes
+
+    template <typename T>
+    std::vector<point2d>
+    snakes(const image2d<T>& input,
+	   unsigned int fen,
+	   unsigned int nb_gen,
+	   float lambda,
+	   unsigned int dmax)
+    {
+      return impl::snakes_(input, fen, nb_gen, lambda, dmax);
+    }
+
+  } // end of namespace oln::appli
+
+} // end of namespace oln
+
+
+#endif // ! OLENA_SNAKES_HH

    