oln 10.56: Gaussian correction.

Index: integre/ChangeLog from Giovanni Palma <giovanni@lrde.epita.fr> * ntg/utils/cast.hh: Make force use a single cast instead of delegate its work to an unsafe type. Index: olena/ChangeLog from Giovanni Palma <giovanni@lrde.epita.fr> * oln/convol/fast_gaussian.hxx: Make the gaussian take care of the borders. Index: olena/oln/convol/fast_gaussian.hxx --- olena/oln/convol/fast_gaussian.hxx Thu, 07 Aug 2003 02:08:21 +0200 david (oln/25_fast_gauss 1.7.1.8.1.4 640) +++ olena/oln/convol/fast_gaussian.hxx Wed, 28 Jan 2004 15:34:38 +0100 palma_g (oln/25_fast_gauss 1.7.1.8.1.4 640) @@ -159,9 +159,9 @@ // Apply on columns. recursivefilter_<float>(img, coef, - oln_point_type(I)(0), - oln_point_type(I)(img.ncols() - 1), - img.ncols(), + oln_point_type(I)(-img.border()), + oln_point_type(I)(img.ncols() - 1 + img.border()), + img.ncols() + 2 * img.border(), oln_dpoint_type(I)(1)); } }; @@ -175,21 +175,20 @@ void doit(abstract::image_with_dim<2, I>& img, const F& coef) { - // Apply on rows. for (coord j = 0; j < img.ncols(); ++j) recursivefilter_<float>(img, coef, - oln_point_type(I)(0, j), - oln_point_type(I)(img.nrows() - 1, j), - img.nrows(), + oln_point_type(I)(-img.border(), j), + oln_point_type(I)(img.nrows() - 1 + img.border(), j), + img.nrows() + 2 * img.border(), oln_dpoint_type(I)(1, 0)); // Apply on columns. for (coord i = 0; i < img.nrows(); ++i) recursivefilter_<float>(img, coef, - oln_point_type(I)(i, 0), - oln_point_type(I)(i, img.ncols() - 1), - img.ncols(), + oln_point_type(I)(i, -img.border()), + oln_point_type(I)(i, img.ncols() - 1 + img.border()), + img.ncols() + 2 * img.border(), oln_dpoint_type(I)(0, 1)); } }; @@ -206,27 +205,27 @@ for (coord j = 0; j < img.nrows(); ++j) for (coord k = 0; k < img.ncols(); ++k) recursivefilter_<float>(img, coef, - oln_point_type(I)(0, j, k), - oln_point_type(I)(img.nslices() - 1, j, k), - img.ncols(), + oln_point_type(I)(-img.border(), j, k), + oln_point_type(I)(img.nslices() - 1 + img.border(), j, k), + img.ncols() + 2 * img.border(), oln_dpoint_type(I)(1, 0, 0)); // Apply on rows. for (coord i = 0; i < img.nslices(); ++i) for (coord k = 0; k < img.ncols(); ++k) recursivefilter_<float>(img, coef, - oln_point_type(I)(i, 0, k), - oln_point_type(I)(i, img.nrows() - 1, k), - img.nrows(), + oln_point_type(I)(i, -img.border(), k), + oln_point_type(I)(i, img.nrows() - 1 + img.border(), k), + img.nrows() + 2 * img.border(), oln_dpoint_type(I)(0, 1, 0)); // Apply on columns. for (coord i = 0; i < img.nslices(); ++i) for (coord j = 0; j < img.nrows(); ++j) recursivefilter_<float>(img, coef, - oln_point_type(I)(i, j, 0), - oln_point_type(I)(i, j, img.ncols() - 1), - img.ncols(), + oln_point_type(I)(i, j, -img.border()), + oln_point_type(I)(i, j, img.ncols() - 1 + img.border()), + img.ncols() + 2 * img.border(), oln_dpoint_type(I)(0, 0, 1)); } }; @@ -235,7 +234,8 @@ typename mute<I, typename convoutput<C,B,oln_value_type(I)>::ret>::ret gaussian_common_(const convert::abstract::conversion<C,B>& c, const abstract::image<I>& in, - const F& coef) + const F& coef, + ntg::float_s sigma) { typename mute<I, ntg::float_s>::ret work_img(in.size()); @@ -243,6 +243,12 @@ for_all(it) work_img[it] = ntg::cast::force<ntg::float_s>(in[it]); + /* FIXME: relation between sigma and the border shouldn't + be linear, so when sigma is big enougth, the signal may + be parasitized by the non signal values. + */ + work_img.border_adapt_mirror(ntg::cast::round<coord>(5 * sigma)); + gaussian_<I::dim>::doit(work_img, coef); /* Convert the result image to the user-requested datatype. @@ -250,6 +256,7 @@ user expects a ntg::float_s image. */ typename mute<I, typename convoutput<C,B,oln_value_type(I)>::ret>::ret out_img(in.size()); + for_all(it) out_img[it] = c(work_img[it]); @@ -271,7 +278,7 @@ sigma, internal::recursivefilter_coef_<float>::DericheGaussian); - return internal::gaussian_common_(c, in, coef); + return internal::gaussian_common_(c, in, coef, sigma); } template <class C, class B, class I> @@ -288,7 +295,7 @@ internal::recursivefilter_coef_<float> ::DericheGaussianFirstDerivative); - return internal::gaussian_common_(c, in, coef); + return internal::gaussian_common_(c, in, coef, sigma); } template <class C, class B, class I> @@ -305,7 +312,7 @@ internal::recursivefilter_coef_<float> ::DericheGaussianSecondDerivative); - return internal::gaussian_common_(c, in, coef); + return internal::gaussian_common_(c, in, coef, sigma); } } // fast Index: integre/ntg/utils/cast.hh --- integre/ntg/utils/cast.hh Thu, 27 Nov 2003 11:26:27 +0100 burrus_n (oln/i/26_cast.hh 1.3.1.11 640) +++ integre/ntg/utils/cast.hh Wed, 28 Jan 2004 16:22:43 +0100 palma_g (oln/i/26_cast.hh 1.3.1.11 640) @@ -70,12 +70,14 @@ | force | `------*/ + // a cast is performed instead of a constructor (unsafe one) call + // because this last one may not be available with the good + // signature. template<class Tdest, class Tsrc> inline const Tdest force(const Tsrc& val) { - ntg_unsafe_type(Tdest) tmp (val); - return tmp; + return (Tdest)val; } /*------. @@ -115,13 +117,13 @@ { // FIXME: this code seems out of date. -#if 0 +#if 1 // KLUDGE: Cast the rounded value to Tdest::value_t before - // returning it as Tdest. Otherwise g++-3.0 complains there + // returning it as Tdest. Otherwise g++-3.* complains there // is no Tdest constructor taking a float argument. return (ntg_storage_type(Tdest)) round(val.exact()); #endif - return round(val.exact()); + // return round(val.exact()); } }; @@ -132,13 +134,13 @@ doit(const float_s& val) { // FIXME: this code seems out of date. -#if 0 +#if 1 // KLUDGE: Cast the rounded value to Tdest::value_t before - // returning it as Tdest. Otherwise g++-3.0 complains there + // returning it as Tdest. Otherwise g++-3.* complains there // is no Tdest constructor taking a float argument. - // return (ntg_storage_type(Tdest)) roundf(val); + return (ntg_storage_type(Tdest)) roundf(val); #endif - return roundf(val); + // return roundf(val); } }; -- Giovanni Palma EPITA - promo 2005 - membre d'EpX - LRDE Mob. : +33 (0)6 60 97 31 74

Giovanni Palma <giovanni@lrde.epita.fr> writes:
Index: integre/ChangeLog from Giovanni Palma <giovanni@lrde.epita.fr>
* ntg/utils/cast.hh: Make force use a single cast instead of delegate its work to an unsafe type.
[...]
Index: integre/ntg/utils/cast.hh --- integre/ntg/utils/cast.hh Thu, 27 Nov 2003 11:26:27 +0100 burrus_n (oln/i/26_cast.hh 1.3.1.11 640) +++ integre/ntg/utils/cast.hh Wed, 28 Jan 2004 16:22:43 +0100 palma_g (oln/i/26_cast.hh 1.3.1.11 640) @@ -70,12 +70,14 @@ | force | `------*/
+ // a cast is performed instead of a constructor (unsafe one) call + // because this last one may not be available with the good + // signature. template<class Tdest, class Tsrc> inline const Tdest force(const Tsrc& val) { - ntg_unsafe_type(Tdest) tmp (val); - return tmp; + return (Tdest)val; }
Ce code pose pb dans quels cas ? Le pb en enlevant le passage par un type unsafe, c'est qu il peut y avoir un check dynamique. L interet de cast::force est justement d assurer qu aucune verification ne sera faite. Si ce n'est pas le cas, autant utiliser directement le cast habituel sans passer par cast::force.

Nicolas Burrus <burrus_n@epita.fr> writes:
Giovanni Palma <giovanni@lrde.epita.fr> writes:
Index: integre/ChangeLog from Giovanni Palma <giovanni@lrde.epita.fr>
* ntg/utils/cast.hh: Make force use a single cast instead of delegate its work to an unsafe type.
[...]
Index: integre/ntg/utils/cast.hh --- integre/ntg/utils/cast.hh Thu, 27 Nov 2003 11:26:27 +0100 burrus_n (oln/i/26_cast.hh 1.3.1.11 640) +++ integre/ntg/utils/cast.hh Wed, 28 Jan 2004 16:22:43 +0100 palma_g (oln/i/26_cast.hh 1.3.1.11 640) @@ -70,12 +70,14 @@ | force | `------*/
+ // a cast is performed instead of a constructor (unsafe one) call + // because this last one may not be available with the good + // signature. template<class Tdest, class Tsrc> inline const Tdest force(const Tsrc& val) { - ntg_unsafe_type(Tdest) tmp (val); - return tmp; + return (Tdest)val; }
Ce code pose pb dans quels cas ? Le pb en enlevant le passage par un type unsafe, c'est qu il peut y avoir un check dynamique. L interet de cast::force est justement d assurer qu aucune verification ne sera faite. Si ce n'est pas le cas, autant utiliser directement le cast habituel sans passer par cast::force.
Le probleme rencontre concerne tout ce qui est dans /doc/demo pour la gaussienne. Concretement, lors du passage d'une image de float a une image de u_int, avec un behavior force, le passage par un type unsafe tente d'appeler unsafe::check qui genere des warning lors de la compil. Ces derniers viennent du fait que g++-3.* ne veut pas generer un int a partir d'un float sans faire de cast explicite. Ce qui m'a amener a faire la modif en ce point, c'est que quand on est dans le behavior force, le code d'origine nous ramenait un un comportement unsafe (par l'appel de cast::force), or il est specifie dans les commentaires, que ces derniers ne doivent pas se comporter de la meme maniere. Je suis ouvert a toute remarque car mon interpretation de la chose n'est peut etre pas la bonne. -- Giovanni Palma EPITA - promo 2005 - membre d'EpX - LRDE Mob. : +33 (0)6 60 97 31 74

Index: integre/ChangeLog from Giovanni Palma <giovanni@lrde.epita.fr>
* ntg/utils/cast.hh: Make force use a single cast instead of delegate its work to an unsafe type.
Index: olena/ChangeLog from Giovanni Palma <giovanni@lrde.epita.fr>
* oln/convol/fast_gaussian.hxx: Make the gaussian take care of the borders.
Existe-t-il un test qui permette de contrôler ça ?

Akim Demaille <akim@epita.fr> writes:
Index: integre/ChangeLog from Giovanni Palma <giovanni@lrde.epita.fr>
* ntg/utils/cast.hh: Make force use a single cast instead of delegate its work to an unsafe type.
Index: olena/ChangeLog from Giovanni Palma <giovanni@lrde.epita.fr>
* oln/convol/fast_gaussian.hxx: Make the gaussian take care of the borders.
Existe-t-il un test qui permette de contrôler ça ?
Pas vraiment, dans la mesure où il est difficile de comparer avec un autre algo (olena n'en comporte qu'un seul). -- Giovanni Palma EPITA - promo 2005 - membre d'EpX - LRDE Mob. : +33 (0)6 60 97 31 74

* oln/convol/fast_gaussian.hxx: Make the gaussian take care of the borders.
Existe-t-il un test qui permette de contrôler ça ?
Pas vraiment, dans la mesure où il est difficile de comparer avec un autre algo (olena n'en comporte qu'un seul).
N'y a-t-il des "invariants" qu'on puisse trouver et vérifier ?

Akim Demaille <akim@epita.fr> writes:
* oln/convol/fast_gaussian.hxx: Make the gaussian take care of the borders.
Existe-t-il un test qui permette de contrôler ça ?
Pas vraiment, dans la mesure où il est difficile de comparer avec un autre algo (olena n'en comporte qu'un seul).
N'y a-t-il des "invariants" qu'on puisse trouver et vérifier ?
Bah c'est un algo d'approximation de la convolution gaussienne par des filtres récursifs, on ne peut donc pas le comparer à l'algo exact de manière naïve. Par contre, si je me souviens bien, dans l'article, il y a une majoration de l'erreur : on peut donc faire la différence des deux images et vérifier. L'article, c'est : Recursively Implementing the Gaussian and its Derivatives Mise en Oeuvre Récursive de la Gaussienne et de ses Dérivées Rachid Deriche Inria http://citeseer.nj.nec.com/deriche93recursively.html A l'époque, on m'avait dit de tester sur un dirac que la masse était conservée mais l'approximation ne le garantie plus sur ce genre de cas pathologiques. Voila si ca peut aider.

Akim Demaille wrote:
Index: integre/ChangeLog from Giovanni Palma <giovanni@lrde.epita.fr>
* ntg/utils/cast.hh: Make force use a single cast instead of delegate its work to an unsafe type.
Index: olena/ChangeLog from Giovanni Palma <giovanni@lrde.epita.fr>
* oln/convol/fast_gaussian.hxx: Make the gaussian take care of the borders.
Existe-t-il un test qui permette de contrôler ça ?
il n'y a rien à "contrôler" la modification de giov rend le calcul "plus juste" qu'il ne l'était, c'est tout
participants (5)
-
Akim Demaille
-
Giovanni Palma
-
Nicolas Burrus
-
Thierry Géraud
-
yann@lrde.epita.fr