https://svn.lrde.epita.fr/svn/oln/branches/cleanup-2008/milena
Index: ChangeLog
from Nicolas Ballas <ballas(a)lrde.epita.fr>
Add directional filters for the gaussian and its derivatives.
* tests/linear/gaussian.cc: Update tests.
* mln/linear/gaussian.hh: Add directional filters.
mln/linear/gaussian.hh | 293 ++++++++++++++++++++++++++++++++++++++++++++---
tests/linear/gaussian.cc | 3
2 files changed, 281 insertions(+), 15 deletions(-)
Index: tests/linear/gaussian.cc
--- tests/linear/gaussian.cc (revision 2826)
+++ tests/linear/gaussian.cc (working copy)
@@ -44,6 +44,8 @@
#include "tests/data.hh"
+#include <mln/trace/all.hh>
+
int main()
@@ -54,6 +56,5 @@
io::pgm::load(lena, MLN_IMG_DIR "/lena.pgm");
image2d<value::int_u8> out = linear::gaussian(lena, 5.1f);
-
io::pgm::save(out, "out.pgm");
}
Index: mln/linear/gaussian.hh
--- mln/linear/gaussian.hh (revision 2826)
+++ mln/linear/gaussian.hh (working copy)
@@ -43,6 +43,8 @@
# include <mln/geom/ncols.hh>
# include <mln/geom/nrows.hh>
+# include <mln/trait/image/props.hh>
+
# include <mln/extension/adjust_fill.hh>
# include <mln/level/stretch.hh>
@@ -51,6 +53,8 @@
# include <vector>
# include <cmath>
+# include <mln/debug/println.hh>
+
namespace mln
{
@@ -247,7 +251,6 @@
}
// Combine results from causal and non-causal parts.
-
current = start;
for (int i = 0; i < len; ++i)
{
@@ -331,8 +334,10 @@
float scale0 = 1 + exp2b0 - 2 * cos0 * expb0;
scale0 *= scale0 * scale0;
- float scaleA = a1 * sin0 * expb0 * (1 + expb0 * (2 * cos0 * (1 + exp2b0) + exp2b0 - 6))
+
- a0 * expb0 * (2 * expb0 * (2 - cos0 * cos0) * (1 - exp2b0) - cos0 * (1 - exp2b0
* exp2b0));
+ float scaleA = a1 * sin0 * expb0 *
+ (1 + expb0 * (2 * cos0 * (1 + exp2b0) + exp2b0 - 6)) +
+ a0 * expb0 * (2 * expb0 * (2 - cos0 * cos0) *
+ (1 - exp2b0) - cos0 * (1 - exp2b0 * exp2b0));
float expb1 = exp(b1);
float exp2b1 = exp(2.0 * b1);
@@ -340,8 +345,10 @@
float scale1 = 1 + exp2b1 - 2 * cos1 * expb1;
scale1 *= scale1 * scale1;
- float scaleC = c1 * sin1 * expb1 * (1 + expb1 * (2 * cos1 * (1 + exp2b1) + exp2b1 - 6))
+
- c0 * expb1 * (2 * expb1 * (2 - cos1 * cos1) * (1 - exp2b1) - cos1 * (1 - exp2b1
* exp2b1));
+ float scaleC = c1 * sin1 * expb1 *
+ (1 + expb1 * (2 * cos1 * (1 + exp2b1) + exp2b1 - 6)) +
+ c0 * expb1 * (2 * expb1 * (2 - cos1 * cos1) *
+ (1 - exp2b1) - cos1 * (1 - exp2b1 * exp2b1));
float sumA = scaleA / scale0;
float sumC = scaleC / scale1;
@@ -350,29 +357,114 @@
return (sumA + sumC);
}
+
+ template <class I, class F>
+ inline
+ void
+ generic_filter_(trait::image::dimension::one_d,
+ Image<I>& img_, const F& coef, int dir)
+ {
+ I& img = exact(img_);
+ mln_precondition(dir < I::site::dim);
+
+
+ recursivefilter_<mln_value(I)>(img, coef,
+ point1d(-img.border()),
+ point1d(geom::ninds(img) - 1 +
+ img.border()),
+ geom::ninds(img) + 2 * img.border(),
+ dpoint1d(1));
+ }
+
template <class I, class F>
inline
void
- generic_filter_(Image<I>& img_, const F& coef)
+ generic_filter_(trait::image::dimension::two_d,
+ Image<I>& img_, const F& coef, int dir)
{
I& img = exact(img_);
+ mln_precondition(dir < I::site::dim);
+
+
+ if (dir == 0)
+ {
// Apply on rows.
for (unsigned j = 0; j < geom::ncols(img); ++j)
recursivefilter_< mln_value(I) >(img, coef,
point2d(-img.border(), j),
- point2d(geom::nrows(img) - 1 + img.border(), j),
+ point2d(geom::nrows(img) - 1 +
+ img.border(), j),
geom::nrows(img) + 2 * img.border(),
dpoint2d(1, 0));
+ }
+ if (dir == 1)
+ {
// Apply on columns.
for (unsigned i = 0; i < geom::nrows(img); ++i)
recursivefilter_< mln_value(I) >(img, coef,
point2d(i, -img.border()),
- point2d(i, geom::ncols(img) - 1 + img.border()),
+ point2d(i, geom::ncols(img) - 1 +
+ img.border()),
geom::ncols(img) + 2 * img.border(),
dpoint2d(0, 1));
}
+ }
+
+ template <class I, class F>
+ inline
+ void
+ generic_filter_(trait::image::dimension::three_d,
+ Image<I>& img_, const F& coef, int dir)
+ {
+ I& img = exact(img_);
+ mln_precondition(dir < I::site::dim);
+
+ if (dir == 0)
+ {
+ // Apply on slices.
+ for (unsigned j = 0; j < geom::nrows(img); ++j)
+ for (unsigned k = 0; k < geom::ncols(img); ++k)
+ recursivefilter_<mln_value(I)>(img, coef,
+ point3d(-img.border(), j , k),
+ point3d(geom::nslis(img) - 1 +
+ img.border(), j, k),
+ geom::nslis(img) + 2 *
+ img.border(),
+ dpoint3d(1, 0, 0));
+ }
+
+
+ if (dir == 1)
+ {
+ // Apply on rows.
+ for (unsigned i = 0; i < geom::nslis(img); ++i)
+ for (unsigned k = 0; k < geom::ncols(img); ++k)
+ recursivefilter_<mln_value(I)>(img, coef,
+ point3d(i, -img.border(), k),
+ point3d(i, geom::nrows(img) - 1 +
+ img.border(), k),
+ geom::nrows(img) + 2 *
+ img.border(),
+ dpoint3d(0, 1, 0));
+ }
+
+ if (dir == 2)
+ {
+ // Apply on columns.
+ for (unsigned i = 0; i < geom::nslis(img); ++i)
+ for (unsigned j = 0; j < geom::nrows(img); ++i)
+ recursivefilter_<mln_value(I)>(img, coef,
+ point3d(i, j, -img.border()),
+ point3d(i, j, geom::ncols(img) -
+ 1 + img.border()),
+ geom::ncols(img) + 2 *
+ img.border(),
+ dpoint3d(0, 0, 1));
+ }
+ }
+
template <class I, class F, class O>
@@ -391,15 +483,42 @@
// On tiny sigma, Derich algorithm doesn't work.
// It is the same thing that to convolve with a Dirac.
if (sigma > 0.006)
- generic_filter_(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::stretch(work_img, out);
+ for (int i = 0; i < I::site::dim; ++i)
+ generic_filter_(mln_trait_image_dimension(I)(),
+ work_img, coef, i);
+
+ // FIXME deal with overflow problem
+ // for instance, when we paste a float image into a int_u8 images.
+ level::paste(work_img, out);
+ }
+
+ template <class I, class F, class O>
+ inline
+ void
+ generic_filter_common_(trait::value::nature::scalar,
+ const Image<I>& in,
+ const F& coef,
+ float sigma,
+ Image<O>& out,
+ int dir)
+ {
+ mln_ch_value(O, float) work_img(exact(in).domain());
+ level::paste(in, work_img);
+ extension::adjust_fill(work_img, 4, 0);
+
+ // On tiny sigma, Derich algorithm doesn't work.
+ // It is the same thing that to convolve with a Dirac.
+ if (sigma > 0.006)
+ generic_filter_(mln_trait_image_dimension(I)(),
+ work_img, coef, dir);
+
+ // FIXME deal with overflow problem
+ // for instance, when we paste a float image into a int_u8 images.
level::paste(work_img, out);
}
+
template <class I, class F, class O>
inline
void
@@ -417,13 +536,145 @@
// On tiny sigma, Derich algorithm doesn't work.
// It is the same thing that to convolve with a Dirac.
if (sigma > 0.006)
- generic_filter_(out, coef);
+ for (int i = 0; i < I::site::dim; ++i)
+ generic_filter_(mln_trait_image_dimension(I)(),
+ out, coef, i);
+ }
+
+ template <class I, class F, class O>
+ inline
+ void
+ generic_filter_common_(trait::value::nature::vectorial,
+ const Image<I>& in,
+ const F& coef,
+ float sigma,
+ Image<O>& out,
+ int dir)
+ {
+ // typedef algebra::vec<3, float> vec3f;
+ // mln_ch_value(O, vec3f) work_img(exact(in).domain());
+ // FIXME : paste does not work (rgb8 -> vec3f).
+ level::paste(in, out);
+
+ // On tiny sigma, Derich algorithm doesn't work.
+ // It is the same thing that to convolve with a Dirac.
+ if (sigma > 0.006)
+ generic_filter_(mln_trait_image_dimension(I)(),
+ out, coef, dir);
}
} // end of namespace mln::linear::impl
+
// Facade.
+
+ /*! Apply an approximated gaussian filter of \p sigma on \p input.
+ * on a specific direction \p dir
+ * if \p dir = 0, the filter is applied on the first image dimension.
+ * if \p dir = 1, the filter is applied on the second image dimension.
+ * And so on...
+ *
+ * \pre input.has_data
+ * \pre dir < dimension(input)
+ */
+ template <class I>
+ inline
+ mln_concrete(I)
+ gaussian(const Image<I>& input, float sigma, int dir)
+ {
+ mln_precondition(exact(input).has_data());
+ mln_precondition(dir < I::site::dim);
+
+ mln_concrete(I) output;
+ initialize(output, input);
+ impl::recursivefilter_coef_ coef(1.68f, 3.735f,
+ 1.783f, 1.723f,
+ -0.6803f, -0.2598f,
+ 0.6318f, 1.997f,
+ sigma, impl::gaussian_norm_coef_);
+
+ impl::generic_filter_common_(mln_trait_value_nature(mln_value(I))(),
+ input, coef, sigma, output, dir);
+ return output;
+ }
+
+
+ /*! Apply an approximated first derivative gaussian filter of \p sigma on
+ * \p input.
+ * on a specific direction \p dir
+ * if \p dir = 0, the filter is applied on the first image dimension.
+ * if \p dir = 1, the filter is applied on the second image dimension.
+ * And so on...
+ *
+ * \pre input.has_data
+ * \pre dir < dimension(input)
+ */
+ template <class I>
+ inline
+ mln_concrete(I)
+ gaussian_1st_derivative(const Image<I>& input, float sigma, int dir)
+ {
+ mln_precondition(exact(input).has_data());
+ mln_precondition(dir < I::site::dim);
+
+ mln_concrete(I) output;
+ initialize(output, input);
+
+ impl::recursivefilter_coef_
+ coef(-0.6472f, -4.531f,
+ 1.527f, 1.516f,
+ 0.6494f, 0.9557f,
+ 0.6719f, 2.072f,
+ sigma, impl::gaussian_1st_deriv_coef_norm_);
+
+ impl::generic_filter_common_(mln_trait_value_nature(mln_value(I))(),
+ input, coef, sigma, output, dir);
+ return output;
+ }
+
+ /*! Apply an approximated second derivative gaussian filter of \p sigma on
+ * \p input.
+ * on a specific direction \p dir
+ * if \p dir = 0, the filter is applied on the first image dimension.
+ * if \p dir = 1, the filter is applied on the second image dimension.
+ * And so on...
+ *
+ * \pre input.has_data
+ * \pre dir < dimension(input)
+ */
+ template <class I>
+ inline
+ mln_concrete(I)
+ gaussian_2nd_derivative(const Image<I>& input, float sigma, int dir)
+ {
+ mln_precondition(exact(input).has_data());
+ mln_precondition(dir < I::site::dim);
+
+ mln_concrete(I) output;
+ initialize(output, input);
+
+ impl::recursivefilter_coef_
+ coef(-1.331f, 3.661f,
+ 1.24f, 1.314f,
+ 0.3225f, -1.738f,
+ 0.748f, 2.166f,
+ sigma, impl::gaussian_2nd_deriv_coef_norm_);
+
+ impl::generic_filter_common_(mln_trait_value_nature(mln_value(I))(),
+ input, coef, sigma, output, dir);
+ return output;
+ }
+
+
+
+
+
+ /*! Apply an approximated gaussian filter of \p sigma on \p input.
+ * This filter is applied in all the input image direction.
+ *
+ * \pre input.has_data
+ */
template <class I>
inline
mln_concrete(I)
@@ -446,6 +697,13 @@
return output;
}
+
+ /*! Apply an approximated first derivative gaussian filter of \p sigma on
+ * \p input
+ * This filter is applied in all the input image direction.
+ *
+ * \pre input.has_data
+ */
template <class I>
inline
mln_concrete(I)
@@ -467,6 +725,13 @@
return output;
}
+
+ /*! Apply an approximated second derivative gaussian filter of \p sigma on
+ * \p input
+ * This filter is applied in all the input image direction.
+ *
+ * \pre input.has_data
+ */
template <class I>
inline
mln_concrete(I)