URL: https://svn.lrde.epita.fr/svn/oln/trunk/milena ChangeLog: 2007-09-05 Simon Nivault <simon.nivault@lrde.epita.fr> Add fast-gaussian * sandbox/nivault/fast_gaussian.hh: New. * sandbox/nivault/tests/test.cc: New. * sandbox/nivault/tests/test: New. * sandbox/nivault/tests: New. fast_gaussian.hh | 288 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ tests/test | 4 tests/test.cc | 42 ++++++++ 3 files changed, 334 insertions(+) Index: trunk/milena/sandbox/nivault/tests/test =================================================================== --- trunk/milena/sandbox/nivault/tests/test (revision 0) +++ trunk/milena/sandbox/nivault/tests/test (revision 1069) @@ -0,0 +1,4 @@ +#!/bin/bash + +g++ -I../../.. -I../.. -I.. test.cc -O1 +./a.out Property changes on: trunk/milena/sandbox/nivault/tests/test ___________________________________________________________________ Name: svn:executable + * Index: trunk/milena/sandbox/nivault/tests/test.cc =================================================================== --- trunk/milena/sandbox/nivault/tests/test.cc (revision 0) +++ trunk/milena/sandbox/nivault/tests/test.cc (revision 1069) @@ -0,0 +1,42 @@ +// Copyright (C) 2007 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, 51 Franklin Street, Fifth Floor, +// 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. + +#include <mln/core/image2d_b.hh> +#include <mln/io/load_pgm.hh> +#include <mln/io/save_pgm.hh> + +#include <fast_gaussian.hh> + +int main() +{ + mln::image2d_b<mln::value::int_u<8u> > im1 = mln::io::load_pgm("../../../img/lena.pgm"); + mln::image2d_b<float> im2; + + mln::linear::gaussian(im1, 0.2, im2); + + // mln::io::save_pgm(im2, "gausslena.pgm"); + } Index: trunk/milena/sandbox/nivault/fast_gaussian.hh =================================================================== --- trunk/milena/sandbox/nivault/fast_gaussian.hh (revision 0) +++ trunk/milena/sandbox/nivault/fast_gaussian.hh (revision 1069) @@ -0,0 +1,288 @@ +// Copyright (C) 2001, 2002, 2003, 2004 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, Inc., 51 Franklin Street, Fifth Floor, +// Boston, MA 02110-1301, 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_CONVOL_FAST_GAUSSIAN_FILTER_HH__ +# define OLENA_CONVOL_FAST_GAUSSIAN_FILTER_HH__ + +#include <mln/core/concept/image.hh> +#include <mln/level/paste.hh> +#include <mln/geom/ncols.hh> + + +namespace mln +{ + + namespace linear + { + + + struct recursivefilter_coef_ + { + + /*! + ** \brief Constructor. + */ + recursivefilter_coef_(float a0, float a1, + float b0, float b1, + float c0, float c1, + float w0, float w1, + float s); + std::vector<float> n, d, nm, dm; + float sumA, sumC; + }; + + recursivefilter_coef_::recursivefilter_coef_(float a0, float a1, + float b0, float b1, + float c0, float c1, + float w0, float w1, + float s) + { + n.reserve(5); + d.reserve(5); + nm.reserve(5); + dm.reserve(5); + + b0 /= s; + b1 /= s; + w0 /= s; + w1 /= s; + + float sin0 = sin(w0); + float sin1 = sin(w1); + float cos0 = cos(w0); + float cos1 = cos(w1); + + sumA = + (2.0 * a1 * exp( b0 ) * cos0 * cos0 - a0 * sin0 * exp( 2.0 * b0 ) + + a0 * sin0 - 2.0 * a1 * exp( b0 )) / + (( 2.0 * cos0 * exp( b0 ) - exp( 2.0 * b0 ) - 1 ) * sin0); + + sumC = + (2.0 * c1 * exp( b1 ) * cos1 * cos1 - c0 * sin1 * exp( 2.0 * b1 ) + + c0 * sin1 - 2.0 * c1 * exp( b1 )) + / (( 2.0 * cos1 * exp( b1 ) - exp( 2.0 * b1 ) - 1 ) * sin1); + + a0 /= (sumA + sumC); + a1 /= (sumA + sumC); + c0 /= (sumA + sumC); + c1 /= (sumA + sumC); + + n[3] = + exp( -b1 - 2*b0 ) * (c1 * sin1 - cos1 * c0) + + exp( -b0 - 2*b1 ) * (a1 * sin0 - cos0 * a0); + n[2] = + 2 * exp(-b0 - b1) * ((a0 + c0) * cos1 * cos0 - + cos1 * a1 * sin0 - + cos0 * c1 * sin1) + + c0 * exp(-2 * b0) + a0 * exp(-2 * b1); + n[1] = + exp(-b1) * (c1 * sin1 - (c0 + 2*a0) * cos1) + + exp(-b0) * (a1 * sin0 - (2*c0 + a0) * cos0); + n[0] = + a0 + c0; + + d[4] = exp(-2 * b0 - 2 * b1); + d[3] = + -2 * cos0 * exp(-b0 - 2*b1) - + 2 * cos1 * exp(-b1 - 2*b0); + d[2] = + 4 * cos1 * cos0 * exp(-b0 - b1) + + exp(-2*b1) + exp(-2*b0); + d[1] = + -2*exp(-b1) * cos1 - 2 * exp(-b0) * cos0; + + for (unsigned i = 1; i <= 3; ++i) + { + dm[i] = d[i]; + nm[i] = n[i] - d[i] * n[0]; + } + dm[4] = d[4]; + nm[4] = -d[4] * n[0]; + } + + + + template <class WorkType, class I> + void + recursivefilter_(I& ima, + const recursivefilter_coef_& c, + const mln_point(I)& start, + const mln_point(I)& finish, + int len, + const mln_dpoint(I)& d) + { + std::vector<WorkType> tmp1(len); + std::vector<WorkType> tmp2(len); + + // The fourth degree approximation implies to have a special + // look on the four first points we consider that there is + // no signal before 0 (to be discussed) + + // -- + // Causal part + + tmp1[0] = + c.n[0]*ima(start); + + tmp1[1] = + c.n[0]*ima(start + d) + + c.n[1]*ima(start) + - c.d[1]*tmp1[0]; + + tmp1[2] = + c.n[0]*ima(start + d + d) + + c.n[1]*ima(start + d) + + c.n[2]*ima(start) + - c.d[1]*tmp1[1] + - c.d[2]*tmp1[0]; + + tmp1[3] = + c.n[0]*ima(start + d + d + d) + + c.n[1]*ima(start + d + d) + + c.n[2]*ima(start + d) + + c.n[3]*ima(start) + - c.d[1]*tmp1[2] - c.d[2]*tmp1[1] + - c.d[3]*tmp1[0]; + + mln_point(I) current(start + d + d + d + d); + for (mln_coord(I) i = 4; i < len; ++i) + { + tmp1[i] = + c.n[0]*ima(current) + + c.n[1]*ima(current - d) + + c.n[2]*ima(current - d - d) + + c.n[3]*ima(current - d - d - d) + - c.d[1]*tmp1[i - 1] - c.d[2]*tmp1[i - 2] + - c.d[3]*tmp1[i - 3] - c.d[4]*tmp1[i - 4]; + current += d; + } + + // Non causal part + + tmp2[len - 1] = 0; + + tmp2[len - 2] = + c.nm[1]*ima(finish); + + tmp2[len - 3] = + c.nm[1]*ima(finish - d) + + c.nm[2]*ima(finish) + - c.dm[1]*tmp2[len-2]; + + tmp2[len - 4] = + c.nm[1]*ima(finish - d - d) + + c.nm[2]*ima(finish - d) + + c.nm[3]*ima(finish) + - c.dm[1]*tmp2[len-3] + - c.dm[2]*tmp2[len-2]; + + current = finish - d - d - d ; + + for (int i = len - 5; i >= 0; --i) + { + tmp2[i] = + c.nm[1]*ima(current) + + c.nm[2]*ima(current + d) + + c.nm[3]*ima(current + d + d) + + c.nm[4]*ima(current + d + d + d) + - c.dm[1]*tmp2[i+1] - c.dm[2]*tmp2[i+2] + - c.dm[3]*tmp2[i+3] - c.dm[4]*tmp2[i+4]; + current -= d; + } + + // Combine results from causal and non-causal parts. + + current = start; + for (int i = 0; i < len; ++i) + { + ima(current) = tmp1[i] + tmp2[i]; + current += d; + } + } + + + template <class I, class F> + void + gaussian_(Image<I>& img_, const F& coef) + { + I& img = exact(img_); + typedef mln_point(I) P; + // Apply on columns. + + recursivefilter_<float>(img, coef, + make::point2d(0, - img.border()), // FIXME + make::point2d(0, geom::ncols(img) - 1 + img.border()), + geom::ncols(img) + 2 * img.border(), + make::dpoint2d(0, 1)); + } + + + template <class I, class F, class O> + void + gaussian_common_(const Image<I>& in, + const F& coef, + float sigma, + Image<O>& out) + { + mln_ch_value(O, float) work_img(exact(in).domain()); + level::paste(in, work_img); + + // On tiny sigma, Derich algorithm doesn't work. + // It is the same thing that to convolve with a Dirac. + if (sigma > 0.006) + gaussian_(work_img, coef); + /* Convert the result image to the user-requested datatype. + FIXME: We are making an unnecessary copy in case the + user expects a ntg::float_s image. */ + level::paste(work_img, out); + } + + + + + template <class I, class O> + void + gaussian(const Image<I>& in, float sigma, + Image<O>& out) + { + recursivefilter_coef_ + coef(1.68f, 3.735f, + 1.783f, 1.723f, + -0.6803f, -0.2598f, + 0.6318f, 1.997f, + sigma); + + gaussian_common_(in, coef, sigma, out); + } + + } + +} + + + +#endif // OLENA_CONVOL_FAST_GAUSSIAN_FILTER_HH__