June 2004 - Olena-patches

oln 10.191: slow Gaussian test.
by Niels van Vliet 13 Jun '04

13 Jun '04

Index: olena/ChangeLog from Niels Van Vliet <niels(a)lrde.epita.fr> * olena/tests/convol/tests/slow_gaussian: New file. Index: olena/tests/convol/tests/slow_gaussian --- olena/tests/convol/tests/slow_gaussian Sun, 13 Jun 2004 13:20:26 +0200 van-vl_n () +++ olena/tests/convol/tests/slow_gaussian Wed, 09 Jun 2004 20:04:10 +0200 van-vl_n (oln/r/11_slow_gauss 644) @@ -0,0 +1,81 @@ +// -*- c++ -*- +// Copyright (C) 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, 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. + +#include "data.hh" +#include "check.hh" + +#include <oln/basics2d.hh> +#include <oln/convol/slow_gaussian.hh> +#include <ntg/int.hh> +#include <oln/utils/md5.hh> +using namespace oln; +using namespace ntg; + +int main() +{ + // Image1d + { + image1d<int_u8> l(4); + l(0) = 1; + l(1) = 6; + l(2) = 64; + l(3) = 0; + + l = convol::slow::gaussian(l, 0.5, 3); + if (l(0) != 2 || l(1) != 12 || l(2) != 51 || l(3) != 7) + { + std::cout << "FAIL" << std::endl; + return true; + } + } + + // Image2d + { + oln::utils::key::value_type data_key[] = + {0x95, 0xe3, 0x63, 0xe0, 0xef, 0x21, 0xcd, 0x34, 0xa2, + 0x15, 0xa3, 0x4e, 0x82, 0x37, 0x40, 0x70}; + oln::utils::key key(data_key); + + image2d<int_u8> lena(rdata("lena.pgm")); + + float_d sigma = 2.5; + if (oln::utils::md5(convol::slow::gaussian(lena, sigma, 3)) != key) + { + std::cout << "FAIL" << std::endl; + return true; + } + } + + // Image3d. Result not tested. + { + image3d<int_u8> i(10, 10, 10); + convol::slow::gaussian(i, 0.2, 3); + } + + std::cout << "OK" << std::endl; +}

1 0

oln 10.190: Behavior Const and Slow Gaussian.
by Niels van Vliet 09 Jun '04

09 Jun '04

Index: olena/ChangeLog from Niels Van Vliet <niels(a)lrde.epita.fr> * olena/oln/core/abstract/behavior.hh (adapt_border): Change constness of the image modified to const. * olena/oln/core/behavior.hh: Likewise. * olena/oln/convol/slow_gaussian.hh: Add behavior. Remove parameter mu. Add doc and example. Add radius_factor. * olena/oln/convol/slow_gaussian.hxx: Clean code. Remove parameter mu. Add radius_factor. * olena/oln/morpho/gradient.inc: Fix doc. Index: olena/oln/morpho/gradient.inc --- olena/oln/morpho/gradient.inc Sun, 14 Mar 2004 19:03:34 +0100 van-vl_n (oln/43_gradient.i 1.16 600) +++ olena/oln/morpho/gradient.inc Wed, 09 Jun 2004 10:48:12 +0200 van-vl_n (oln/43_gradient.i 1.16 600) @@ -71,7 +71,7 @@ ** IMG_OUT "oln_morpho_beucher_gradient.pbm"); ** return 0; ** } -** \encode +** \endcode ** ** \image html lena128_pgm.png ** \image latex lena128_pgm.png Index: olena/oln/core/abstract/behavior.hh --- olena/oln/core/abstract/behavior.hh Thu, 11 Mar 2004 17:12:19 +0100 thivol_d (oln/j/46_behavior.h 1.3 600) +++ olena/oln/core/abstract/behavior.hh Wed, 09 Jun 2004 19:01:52 +0200 van-vl_n (oln/j/46_behavior.h 1.3 600) @@ -58,7 +58,7 @@ ** Adapt the border of an image regarding the kind of behavior wanted. */ template <class I> - void adapt_border(oln::abstract::image &im, coord border_size) const + void adapt_border(const oln::abstract::image &im, coord border_size) const { mlc_dispatch(adapt_border)(im, border_size); }; Index: olena/oln/core/behavior.hh --- olena/oln/core/behavior.hh Fri, 12 Mar 2004 20:17:58 +0100 thivol_d (oln/j/47_behavior.h 1.3 600) +++ olena/oln/core/behavior.hh Wed, 09 Jun 2004 19:03:10 +0200 van-vl_n (oln/j/47_behavior.h 1.3 600) @@ -47,7 +47,8 @@ typedef mlc_exact_vt_type(self_type, Exact) exact_type; template <class I> - void adapt_border_impl(oln::abstract::image &im, coord border_size) const + void adapt_border_impl(const oln::abstract::image &im, + coord border_size) const { im.border_adapt_mirror(border_size); }; @@ -75,9 +76,11 @@ }; template <class I> - void adapt_border_impl(abstract::image &im, coord border_size) const + void adapt_border_impl(const abstract::image &im, + coord border_size) const { - im.border_adapt_assign(border_size, ntg::cast::force<oln_value_type(I)>(value_)); + im.border_adapt_assign(border_size, + ntg::cast::force<oln_value_type(I)>(value_)); }; protected: @@ -101,7 +104,8 @@ typedef mlc_exact_vt_type(self_type, Exact) exact_type; template <class I> - void adapt_border_impl(abstract::image &im, coord border_size) const + void adapt_border_impl(const abstract::image &im, + coord border_size) const { im.border_adapt_copy(border_size); }; Index: olena/oln/convol/slow_gaussian.hh --- olena/oln/convol/slow_gaussian.hh Wed, 12 May 2004 23:10:18 +0200 van-vl_n (oln/r/7_slow_gauss 1.1 600) +++ olena/oln/convol/slow_gaussian.hh Wed, 09 Jun 2004 19:14:55 +0200 van-vl_n (oln/r/7_slow_gauss 1.1 600) @@ -31,6 +31,7 @@ # include <oln/core/image.hh> # include <ntg/float.hh> # include "convolution.hh" +# include <oln/core/behavior.hh> namespace oln { namespace convol { @@ -45,33 +46,67 @@ { }; - /*! Gaussian filter. + /*! Slow Gaussian filter. ** - ** \todo FIXME: this class has been done on hurry. ** \arg in input image. ** \arg sigma sigma. - ** \arg mu point corresponding to the center. + ** \arg radius_factor the kernel of the convolution is done with a + ** kernel null at a distance of sigma * radius_factor. + ** \arg behavior Object to know how to work on borders. By default + ** mirror is used. + ** \code + ** #include <oln/basics2d.hh> + ** #include <oln/convol/slow_gaussian.hh> + ** #include <ntg/int.hh> + ** + ** using namespace oln; + ** using namespace ntg; + ** + ** int main() + ** { + ** image2d<int_u8> lena(IMG_IN "lena.pgm"); + ** + ** float_d sigma = 2.5; + ** lena = convol::slow::gaussian(lena, sigma, 3); + ** save(lena, IMG_OUT "oln_convol_slow_gaussian.pgm"); + ** } + ** \endcode + ** \image html lena_pgm.png + ** \image latex lena_pgm.png + ** => + ** \image html oln_convol_slow_gaussian.pgm + ** \image latex oln_convol_slow_gaussian.pgm */ - template <class I> + + template <class I, class BE> oln_concrete_type(I) gaussian(const oln::abstract::non_vectorial_image & in, const ntg::float_d sigma, - const oln_point_type(I) & mu = oln_point_type(I)()) + const ntg::float_d radius_factor, + const oln::abstract::behavior<BE> &behavior) { - static const typename gaussian_kernel<I::dim>::ret k - = gaussian_kernel<I::dim>::kernel_norm(sigma, mu); - in.border_adapt_mirror(k.delta()); - //FIXME: we should use fast convol. - //FIXME: This is hugely. - typename mute<I, double>::ret im = oln::convol::slow::convolve<double> - (in, - gaussian_kernel<I::dim>::kernel(sigma, mu)); + const typename gaussian_kernel<I::dim>::ret k + = gaussian_kernel<I::dim>::kernel(sigma, radius_factor); + behavior.adapt_border(in, coord(k.delta())); + + typename mute<I, ntg::float_d>::ret im = + oln::convol::slow::convolve<ntg::float_d> (in, k); oln_concrete_type(I) out(in.size()); oln_iter_type(I) it(out); for_all(it) out[it] = oln_value_type(I)(im[it]); return out; } + + //! Gaussian filter, with default borders. + template <class I> + oln_concrete_type(I) + gaussian(const oln::abstract::non_vectorial_image & in, + const ntg::float_d sigma, + const ntg::float_d radius_factor) + { + return gaussian(in, sigma, radius_factor, mirror_behavior<>()); + } }//End namespace slow }// End namespace convol }// End namespace oln Index: olena/oln/convol/slow_gaussian.hxx --- olena/oln/convol/slow_gaussian.hxx Wed, 12 May 2004 23:10:18 +0200 van-vl_n (oln/r/8_slow_gauss 1.1 600) +++ olena/oln/convol/slow_gaussian.hxx Wed, 09 Jun 2004 18:31:54 +0200 van-vl_n (oln/r/8_slow_gauss 1.1 600) @@ -32,9 +32,6 @@ # include <oln/core/w_window2d.hh> # include <oln/core/w_window3d.hh> -//! SIGMA_COEF times sigma is the size of the kernel. -# define SIGMA_COEF 3 - namespace oln { namespace convol { namespace slow { @@ -44,14 +41,15 @@ and w_windownd does not supprot concrete_type, etc */ template<class T> - T + inline T normalise(const T &in) { T w = in; - double sum = 0; + ntg::float_d sum = 0.; for (unsigned i = 0; i < in.card(); ++i) sum += in.w(i); sum = 1. / sum; + assert(finite(sum)); for (unsigned i = 0; i < in.card(); ++i) w.set(w.dp(i), in.w(i) * sum); return w; @@ -63,26 +61,38 @@ template<> struct gaussian_kernel<1> { - typedef oln::w_window1d<double> ret; + typedef oln::w_window1d<ntg::float_d> ret; enum {dim = 1}; - static w_window1d<double> - kernel_norm(double sigma, const point1d & mu = point1d(0)) - { - return internal::normalise(kernel(sigma, mu)); - } - - w_window1d<double> - static kernel(double sigma, const point1d & mu = point1d(0)) - { - assert(sigma >= 0); - w_window1d<double> w; - - const int size = int(sigma * SIGMA_COEF); - const double inv_sigma_sqrt_2pi = 1. / (sigma * sqrt(M_PI * 2.)); + //! Return a discrete Gaussian kernel, that is equal to zero + // at a distance of radius_factor * sigma. + static w_window1d<ntg::float_d> + kernel(ntg::float_d sigma, + ntg::float_d radius_factor) + { + precondition(sigma > 0.); + precondition(radius_factor >= 0.); + return internal::normalise(kernel_values(sigma, radius_factor)); + } + + //! Return values of the Gaussian kernel at discrete points. + // + // Note: the integral is not equal to 1 (discrete grid and 0 outside + // radius_factor * sigma). You should use kernel(sigma, radius_factor). + static w_window1d<ntg::float_d> + kernel_values(ntg::float_d sigma, + ntg::float_d radius_factor) + { + precondition(sigma > 0.); + precondition(radius_factor >= 0.); + + w_window1d<ntg::float_d> w; + + const int size = int(sigma * radius_factor); + const ntg::float_d inv_sigma_sqrt_2pi + = 1. / (sigma * sqrt(M_PI * 2.)); for (int x = -size; x <= +size; ++x) - w.add(x, inv_sigma_sqrt_2pi * exp( -(x - mu.nth(0)) * (x - mu.nth(0)) / (2 * sigma * sigma)) - ); + w.add(x, inv_sigma_sqrt_2pi * exp(-x * x /(2 * sigma * sigma))); return w; } }; @@ -91,27 +101,42 @@ template<> struct gaussian_kernel<2> { - typedef w_window2d<double> ret; + typedef w_window2d<ntg::float_d> ret; enum {dim = 2}; - static w_window2d<double> - kernel_norm(double sigma, const point2d & mu = point2d(0, 0)) - { - return internal::normalise(kernel(sigma, mu)); - } - w_window2d<double> - static kernel(double sigma, const point2d & mu = point2d(0, 0)) - { - assert(sigma >= 0); - w_window2d<double> w; - const int size = int(sigma * SIGMA_COEF); - const double inv_sigma_sigma_pi_2 = 1. / (sigma * sigma * M_PI * 2.); + //! Return a discrete Gaussian kernel, that is equal to zero at a + // distance of radius_factor * sigma. + static w_window2d<ntg::float_d> + kernel(ntg::float_d sigma, + ntg::float_d radius_factor) + { + precondition(sigma > 0.); + precondition(radius_factor >= 0.); + + return internal::normalise(kernel_values(sigma, radius_factor)); + } + + //! Return values of the Gaussian kernel at discrete points. + // + // Note: the integral is not equal to 1 (discrete grid and 0 outside + // radius_factor * sigma). You should use kernel(sigma, radius_factor). + static w_window2d<ntg::float_d> + kernel_values(ntg::float_d sigma, + ntg::float_d radius_factor) + { + precondition(sigma > 0.); + precondition(radius_factor >= 0.); + w_window2d<ntg::float_d> w; + + const int size = int(sigma * radius_factor); + const ntg::float_d inv_sigma_sigma_pi_2 = 1. / (sigma * sigma * + M_PI * 2.); for (int x = -size; x <= +size; ++x) for (int y = -size; y <= +size; ++y) + if (x * x + y * y <= size * size) w.add(x, y, inv_sigma_sigma_pi_2 - * exp(- ((x - mu.nth(0)) * (x - mu.nth(0)) + (y - mu.nth(1)) * (y - mu.nth(1)) - ) + * exp(- (x * x + y * y) / ( 2. * sigma * sigma) ) ); @@ -123,31 +148,47 @@ template<> struct gaussian_kernel<3> { - typedef w_window3d<double> ret; + typedef w_window3d<ntg::float_d> ret; enum {dim = 3}; - static w_window3d<double> - kernel_norm(double sigma, const point3d & mu = point3d(0, 0, 0)) + //! Return a discrete Gaussian kernel, that is equal to zero at a + // distance of radius_factor * sigma. + static w_window3d<ntg::float_d> + kernel(ntg::float_d sigma, + ntg::float_d radius_factor) { - return internal::normalise(kernel(sigma, mu)); + precondition(sigma > 0.); + precondition(radius_factor >= 0.); + + return internal::normalise(kernel_values(sigma, radius_factor)); } - static w_window3d<double> - kernel(double sigma, const point3d & mu = point3d(0, 0, 0)) + + + //! Return values of the Gaussian kernel at discrete points. + // + // Note: the integral is not equal to 1 (discrete grid and 0 outside + // radius_factor * sigma). You should use kernel(sigma, radius_factor). + static w_window3d<ntg::float_d> + kernel_values(ntg::float_d sigma, + ntg::float_d radius_factor) { - assert(sigma >= 0); - w_window3d<double> w; + precondition(sigma > 0.); + precondition(radius_factor >= 0.); + + w_window3d<ntg::float_d> w; - const int size = int(sigma * SIGMA_COEF); - const double k = 1. / (sigma * sigma * sigma * sqrt((M_PI * 2.) * (M_PI * 2.) * (M_PI * 2.))); + const int size = int(sigma * radius_factor); + const ntg::float_d k = 1. / (sigma * sigma * sigma * + sqrt((M_PI * 2.) * (M_PI * 2.) * + (M_PI * 2.))); for (int x = -size; x <= +size; ++x) for (int y = -size; y <= +size; ++y) for (int z = -size; z <= +size; ++z) + if (x * x + y * y + z * z <= size) w.add(x, y, z, k * - exp(-((x - mu.nth(0)) * (x - mu.nth(0)) + - (y - mu.nth(1)) * (y - mu.nth(1)) + - (z - mu.nth(2)) * (z - mu.nth(2)) - ) / (2. * sigma * sigma))); + exp(-(x * x + y * y + z * z) + ) / (2. * sigma * sigma)); return w; } };

1 0

oln 10.189: Change functions to functors
by Damien Thivolle 08 Jun '04

08 Jun '04

Index: olena/ChangeLog from Damien Thivolle <damien(a)lrde.epita.fr> * oln/transforms/shapelets.hh: Change functions to functors, bug fixes. Index: olena/oln/transforms/shapelets.hh --- olena/oln/transforms/shapelets.hh Sun, 06 Jun 2004 23:17:15 +0200 thivol_d (oln/r/10_shapelets. 1.1 600) +++ olena/oln/transforms/shapelets.hh Mon, 07 Jun 2004 19:39:27 +0200 thivol_d (oln/r/10_shapelets. 1.2 600) @@ -7,71 +7,89 @@ #include <oln/core/point2d.hh> #include <oln/level/fill.hh> #include <vector> +#include <bits/stl_function.h> namespace oln { - inline double - hermite(int n, double x) + namespace internal + { + template <class T> + struct hermite : public std::binary_function<T, int, T> + { + T operator()(unsigned n, T x) { - assert(n >= 0); if (n == 0) return 1; if (n == 1) return 2. * x; - return 2. * (x * hermite(n - 1, x) - (n - 1) * hermite(n - 2, x)); + return 2. * x * (*this)(n - 1, x) - 2. * double(n - 1) * (*this)(n - 2, x); } + }; - inline double - fact(double n) + template <class T> + struct fact : public std::unary_function<T, T> + { + T operator() (T n) { precondition(n >= 0); + T res = 1; - double res = n; - - if (n <= 1) - return 1; - while (n > 1) + while (n > 0) { res *= n; n--; } return res; } + }; - inline double - beta(int n, double x) + template <class T> + struct shapelets_beta : public std::binary_function<T, int, T> + { + T operator()(int n, T x) { assert(n >= 0); - const double c = sqrt(pow(2., n) * sqrt(M_PI) * fact(n)); + const double c = sqrt(pow(2.0, n) * sqrt(M_PI) * T(fact<int>()(n))); - assert(finite(hermite(n, x) * exp(x * x / -2.0) / c)); - return hermite(n, x) * exp(x * x / -2.0) / c; + return hermite<T>()(n, x) * exp(x * x / -2.0) / c; } + }; + template <unsigned Dim, class T> + struct shapelets_basis + { + }; - - inline double - shapelet_basis_1d(int n, point1d p, double b) + template <class T> + struct shapelets_basis<1, T> + { + T operator() (int n, point1d p, T b) { assert(n >= 0 && b >= 0); - return beta(n, p.col() / b) / sqrt(b); + return shapelets_beta<T>(n, p.col() / b) / sqrt(b); } + }; - double - shapelet_basis_2d(int m, int n, point2d p, double b) + template <class T> + struct shapelets_basis<2, T> + { + T operator() (int m, int n, T x, T y, T b) { assert(m >= 0 && n >= 0 && b >= 0); - return beta(m, p.col() / b) * beta(n, p.row() / b) / b; + return shapelets_beta<T>()(m, x / b) * shapelets_beta<T>()(n, y / b) / b; } + }; + //FIXME: ima should be non_vectorial template <class I> std::vector<double> - to_shapelet_basis_2d(const oln::abstract::image_with_dim<2, I>& ima, + shapelets(const oln::abstract::image_with_dim<2, I>& ima, int m, int n, double b) { assert(m >= 0 && n >= 0 && b >= 0); + shapelets_basis<2, double> func; // Center of the image const int col = ima.ncols() / 2; @@ -90,7 +108,7 @@ for_all(it) { s += (ima[it] - ntg_max_val(oln_value_type(I)) / 2) * - shapelet_basis_2d(k, l, point2d(it.row() - row, it.col() - col), b); + func(k, l, double(it.row() - row), double(it.col() - col), b); } res[k * n + l] = s; } @@ -120,18 +138,20 @@ // Add the value at the point for_all(it) resf[it] += vec[k * n + l] * - shapelet_basis_2d(k, l, point2d(it.row() - nrows / 2, it.col() - ncols / 2), b); + shapelets_basis<2, double>()(k, l, double(it.row() - nrows / 2), double(it.col() - ncols / 2), b); image2d<D> res(oln::image2d_size(nrows, ncols, 0)); for_all(it) { - if ( !(resf[it] >= 0 && resf[it] < 256)) + if ( !(resf[it] >= 0 && resf[it] <= 255)) { - std::cout << "err:" << resf[it] << std::endl; + std::cout << "- "; + std::cout << "err:" << resf[it]; // min = (resf[it] < min)? resf[it]:min; // max = (resf[it] >max)? resf[it]:max; res[it] = resf[it] < 0 ? 0 : ntg_max_val(D); + std::cout << " -" << std::endl; } else res[it] = resf[it]; @@ -142,3 +162,4 @@ } } +} -- Damien Thivolle damien.thivolle(a)lrde.epita.fr

4 3

oln 10.188: Add shapelets tranforms
by Damien Thivolle 08 Jun '04

08 Jun '04

Index: ChangeLog from Damien Thivolle <damien(a)lrde.epita.fr> * config/oln.m4: Add detection of Doxygen 1.3.7 Index: olena/ChangeLog from Damien Thivolle <damien(a)lrde.epita.fr> * oln/transforms/shapelets.hh: Add file, it contains the shapelets transforms implementation for olena. * oln/makefile.src: Add oln/transforms/shapelets.hh to the distribution's files. Index: config/oln.m4 --- config/oln.m4 Wed, 14 Apr 2004 20:42:25 +0200 thivol_d (oln/j/15_oln.m4 1.38.1.5 640) +++ config/oln.m4 Sun, 06 Jun 2004 23:02:13 +0200 thivol_d (oln/j/15_oln.m4 1.38.1.5 640) @@ -348,7 +348,7 @@ [dnl AC_CACHE_CHECK([for doxygen >= 1.3.4], [ac_cv_doxygen_version], - [if $DOXYGEN --version 2>&1 | grep ['1.3.[456]'] > /dev/null 2>&1; then + [if $DOXYGEN --version 2>&1 | grep ['1.3.[4567]'] > /dev/null 2>&1; then ac_cv_doxygen_version=recent else ac_cv_doxygen_version=old Index: olena/oln/makefile.src --- olena/oln/makefile.src Tue, 13 Apr 2004 13:42:44 +0200 palma_g (oln/r/4_makefile.s 1.1 644) +++ olena/oln/makefile.src Sun, 06 Jun 2004 23:11:13 +0200 thivol_d (oln/r/4_makefile.s 1.1 644) @@ -212,6 +212,7 @@ transforms/fft.hh \ transforms/dwt.hh \ transforms/wavelet_coeffs.hh \ + transforms/shapelets.hh \ utils/copy.hh \ utils/key.hh \ utils/key.hxx \ Index: olena/oln/transforms/shapelets.hh --- olena/oln/transforms/shapelets.hh Sun, 06 Jun 2004 23:15:43 +0200 thivol_d () +++ olena/oln/transforms/shapelets.hh Sun, 06 Jun 2004 23:07:28 +0200 thivol_d (oln/r/10_shapelets. 644) @@ -0,0 +1,144 @@ +#include <ntg/int.hh> +#include <oln/basics2d.hh> + + +#include <oln/core/image.hh> +#include <oln/core/point1d.hh> +#include <oln/core/point2d.hh> +#include <oln/level/fill.hh> +#include <vector> + +namespace oln +{ + inline double + hermite(int n, double x) + { + assert(n >= 0); + if (n == 0) + return 1; + if (n == 1) + return 2. * x; + return 2. * (x * hermite(n - 1, x) - (n - 1) * hermite(n - 2, x)); + } + + inline double + fact(double n) + { + precondition(n >= 0); + + double res = n; + + if (n <= 1) + return 1; + while (n > 1) + { + res *= n; + n--; + } + return res; + } + + inline double + beta(int n, double x) + { + assert(n >= 0); + const double c = sqrt(pow(2., n) * sqrt(M_PI) * fact(n)); + + assert(finite(hermite(n, x) * exp(x * x / -2.0) / c)); + return hermite(n, x) * exp(x * x / -2.0) / c; + } + + + + inline double + shapelet_basis_1d(int n, point1d p, double b) + { + assert(n >= 0 && b >= 0); + return beta(n, p.col() / b) / sqrt(b); + } + + double + shapelet_basis_2d(int m, int n, point2d p, double b) + { + assert(m >= 0 && n >= 0 && b >= 0); + return beta(m, p.col() / b) * beta(n, p.row() / b) / b; + } + + //FIXME: ima should be non_vectorial + template <class I> + std::vector<double> + to_shapelet_basis_2d(const oln::abstract::image_with_dim<2, I>& ima, + int m, int n, + double b) + { + assert(m >= 0 && n >= 0 && b >= 0); + + // Center of the image + const int col = ima.ncols() / 2; + const int row = ima.nrows() / 2; + + std::vector<double> res(m * n); + + // For all vector in the base of shapelets + for (int k = 0; k < m; k++) + for (int l = 0; l < n; l++) + // res[k * n + l] = scalar product of the image and the function + { + double s = 0.0; + + oln_iter_type(I) it(ima); + for_all(it) + { + s += (ima[it] - ntg_max_val(oln_value_type(I)) / 2) * + shapelet_basis_2d(k, l, point2d(it.row() - row, it.col() - col), b); + } + res[k * n + l] = s; + } + return res; + } + + + template <class D> + image2d<D> + unshapelets(const std::vector<double> &vec, + int m, int n, double b, + coord nrows, coord ncols) + { + precondition(nrows >= 1 && ncols >= 1 && + m >= 0 && n >= 0 && b >= 0); + + + typedef image2d<double> tmp_type; + tmp_type resf(oln::image2d_size(nrows, ncols, 0)); + + oln_iter_type(tmp_type) it(resf); + + oln::level::fill(resf, ntg_max_val(D) / 2 ); + // For each vector + for (int k = 0; k < m; k++) + for (int l = 0; l < n; l++) + // Add the value at the point + for_all(it) + resf[it] += vec[k * n + l] * + shapelet_basis_2d(k, l, point2d(it.row() - nrows / 2, it.col() - ncols / 2), b); + + image2d<D> res(oln::image2d_size(nrows, ncols, 0)); + + for_all(it) + { + if ( !(resf[it] >= 0 && resf[it] < 256)) + { + std::cout << "err:" << resf[it] << std::endl; +// min = (resf[it] < min)? resf[it]:min; +// max = (resf[it] >max)? resf[it]:max; + res[it] = resf[it] < 0 ? 0 : ntg_max_val(D); + } + else + res[it] = resf[it]; + } + // std::cout << "min:" << min + // << " max:" << max << std::endl; + return res; + } + +} -- Damien Thivolle damien.thivolle(a)lrde.epita.fr

3 4