Olena-patches August 2013

olena-patches@lrde.epita.fr

5 participants
171 discussions

olena: branch fft updated: olena-2.0-638-g773114f

by git＠lrde.epita.fr

This is an automated email from the git hooks/post-receive script. It was generated because a ref change was pushed to the repository containing the project "Olena, a generic and efficient image processing platform". The branch fft has been updated via 773114f96eeb083d260eed6e5a508d9b93d9c9b2 (commit) from 8b31df65d090f1ac9965f0feb4d3bb68a17948f1 (commit) Those revisions listed above that are new to this repository have not appeared on any other notification email; so we list those revisions in full, below. - Log ----------------------------------------------------------------- 773114f Run FFT unit tests only when FFTW3 support is enabled. ----------------------------------------------------------------------- Summary of changes: milena/ChangeLog | 8 ++++++++ milena/tests/unit_test/Makefile.am | 3 ++- milena/tests/unit_test/unit-tests.mk | 10 +++++++++- 3 files changed, 19 insertions(+), 2 deletions(-) hooks/post-receive -- Olena, a generic and efficient image processing platform

11 years, 4 months

olena: olena-2.0-638-g0a07893 Run FFT unit tests only when FFTW3 support is enabled.

by Roland Levillain

* tests/unit_test/cond_tests_fftw3: New. * tests/unit_test/Makefile.am (COND_TESTS): Add cond_tests_fftw3. * tests/unit_test/unit-tests.mk: Regen. --- milena/ChangeLog | 8 ++++++++ milena/tests/unit_test/Makefile.am | 3 ++- milena/tests/unit_test/cond_tests_fftw3 | 1 + milena/tests/unit_test/unit-tests.mk | 10 +++++++++- 4 files changed, 20 insertions(+), 2 deletions(-) create mode 100644 milena/tests/unit_test/cond_tests_fftw3 diff --git a/milena/ChangeLog b/milena/ChangeLog index 72bfb1a..f7fb36d 100644 --- a/milena/ChangeLog +++ b/milena/ChangeLog @@ -1,3 +1,11 @@ +2013-08-26 Roland Levillain <roland(a)lrde.epita.fr> + + Run FFT unit tests only when FFTW3 support is enabled. + + * tests/unit_test/cond_tests_fftw3: New. + * tests/unit_test/Makefile.am (COND_TESTS): Add cond_tests_fftw3. + * tests/unit_test/unit-tests.mk: Regen. + 2012-10-10 Roland Levillain <roland(a)lrde.epita.fr> Split the interface of the FFT off from its implementation. diff --git a/milena/tests/unit_test/Makefile.am b/milena/tests/unit_test/Makefile.am index 37df4c8..1c3016f 100644 --- a/milena/tests/unit_test/Makefile.am +++ b/milena/tests/unit_test/Makefile.am @@ -1,4 +1,4 @@ -# Copyright (C) 2008, 2009, 2010 EPITA Research and Development +# Copyright (C) 2008, 2009, 2010, 2013 EPITA Research and Development # Laboratory (LRDE). # # This file is part of Olena. @@ -19,6 +19,7 @@ include $(top_srcdir)/milena/tests/tests.mk COND_TESTS = cond_tests_boost-tuple \ cond_tests_cfitsio \ + cond_tests_fftw3 \ cond_tests_gdcm \ cond_tests_magickxx \ cond_tests_poppler \ diff --git a/milena/tests/unit_test/cond_tests_fftw3 b/milena/tests/unit_test/cond_tests_fftw3 new file mode 100644 index 0000000..c09df70 --- /dev/null +++ b/milena/tests/unit_test/cond_tests_fftw3 @@ -0,0 +1 @@ +mln/transform/fft.hh diff --git a/milena/tests/unit_test/unit-tests.mk b/milena/tests/unit_test/unit-tests.mk index d1c1355..01272a0 100644 --- a/milena/tests/unit_test/unit-tests.mk +++ b/milena/tests/unit_test/unit-tests.mk @@ -24,6 +24,15 @@ mln_io_fits_load_LDFLAGS= ${CFITSIO_LDFLAGS} ${AM_LDFLAGS} endif HAVE_CFITSIO # Starting a conditional unit test list. +if HAVE_FFTW3 +check_PROGRAMS += \ +mln_transform_fft + +mln_transform_fft_CPPFLAGS= ${FFTW3_CPPFLAGS} -DHAVE_FFTW3 ${AM_CPPFLAGS} +mln_transform_fft_LDFLAGS= ${FFTW3_LDFLAGS} ${AM_LDFLAGS} +endif HAVE_FFTW3 + +# Starting a conditional unit test list. if HAVE_GDCM check_PROGRAMS += \ mln_io_dicom_all \ @@ -1183,7 +1192,6 @@ mln_transform_distance_and_influence_zone_geodesic \ mln_transform_distance_front \ mln_transform_distance_geodesic \ mln_transform_essential \ -mln_transform_fft \ mln_transform_hough \ mln_transform_influence_zone_front \ mln_transform_influence_zone_geodesic \ -- 1.7.10.4

11 years, 4 months

olena: olena-2.0-638-g773114f Run FFT unit tests only when FFTW3 support is enabled.

by Roland Levillain

* tests/unit_test/cond_tests_fftw3: New. * tests/unit_test/Makefile.am (COND_TESTS): Add cond_tests_fftw3. * tests/unit_test/unit-tests.mk: Regen. --- milena/ChangeLog | 8 ++++++++ milena/tests/unit_test/Makefile.am | 3 ++- milena/tests/unit_test/unit-tests.mk | 10 +++++++++- 3 files changed, 19 insertions(+), 2 deletions(-) diff --git a/milena/ChangeLog b/milena/ChangeLog index 72bfb1a..f7fb36d 100644 --- a/milena/ChangeLog +++ b/milena/ChangeLog @@ -1,3 +1,11 @@ +2013-08-26 Roland Levillain <roland(a)lrde.epita.fr> + + Run FFT unit tests only when FFTW3 support is enabled. + + * tests/unit_test/cond_tests_fftw3: New. + * tests/unit_test/Makefile.am (COND_TESTS): Add cond_tests_fftw3. + * tests/unit_test/unit-tests.mk: Regen. + 2012-10-10 Roland Levillain <roland(a)lrde.epita.fr> Split the interface of the FFT off from its implementation. diff --git a/milena/tests/unit_test/Makefile.am b/milena/tests/unit_test/Makefile.am index 37df4c8..1c3016f 100644 --- a/milena/tests/unit_test/Makefile.am +++ b/milena/tests/unit_test/Makefile.am @@ -1,4 +1,4 @@ -# Copyright (C) 2008, 2009, 2010 EPITA Research and Development +# Copyright (C) 2008, 2009, 2010, 2013 EPITA Research and Development # Laboratory (LRDE). # # This file is part of Olena. @@ -19,6 +19,7 @@ include $(top_srcdir)/milena/tests/tests.mk COND_TESTS = cond_tests_boost-tuple \ cond_tests_cfitsio \ + cond_tests_fftw3 \ cond_tests_gdcm \ cond_tests_magickxx \ cond_tests_poppler \ diff --git a/milena/tests/unit_test/unit-tests.mk b/milena/tests/unit_test/unit-tests.mk index d1c1355..01272a0 100644 --- a/milena/tests/unit_test/unit-tests.mk +++ b/milena/tests/unit_test/unit-tests.mk @@ -24,6 +24,15 @@ mln_io_fits_load_LDFLAGS= ${CFITSIO_LDFLAGS} ${AM_LDFLAGS} endif HAVE_CFITSIO # Starting a conditional unit test list. +if HAVE_FFTW3 +check_PROGRAMS += \ +mln_transform_fft + +mln_transform_fft_CPPFLAGS= ${FFTW3_CPPFLAGS} -DHAVE_FFTW3 ${AM_CPPFLAGS} +mln_transform_fft_LDFLAGS= ${FFTW3_LDFLAGS} ${AM_LDFLAGS} +endif HAVE_FFTW3 + +# Starting a conditional unit test list. if HAVE_GDCM check_PROGRAMS += \ mln_io_dicom_all \ @@ -1183,7 +1192,6 @@ mln_transform_distance_and_influence_zone_geodesic \ mln_transform_distance_front \ mln_transform_distance_geodesic \ mln_transform_essential \ -mln_transform_fft \ mln_transform_hough \ mln_transform_influence_zone_front \ mln_transform_influence_zone_geodesic \ -- 1.7.10.4

11 years, 4 months

olena: olena-2.0-642-gb1efc3e Remove old Subversion keywords from ChangeLogs.

by Roland Levillain

--- ChangeLog | 4 ---- milena/ChangeLog | 4 ---- scribo/ChangeLog | 4 ---- swilena/ChangeLog | 4 ---- 4 files changed, 16 deletions(-) diff --git a/ChangeLog b/ChangeLog index 0447440..4b6e319 100644 --- a/ChangeLog +++ b/ChangeLog @@ -4168,10 +4168,6 @@ * vcs/oln.rb: New. -$Rev$ -$Id$ -$Date$ - Local Variables: ispell-local-dictionary: "american" coding: iso-latin-1 diff --git a/milena/ChangeLog b/milena/ChangeLog index 386b20e..df99b3e 100644 --- a/milena/ChangeLog +++ b/milena/ChangeLog @@ -59261,10 +59261,6 @@ * vcs/milena.rb: New file. -$Rev$ -$Id$ -$Date$ - Local Variables: ispell-local-dictionary: "american" coding: iso-latin-1 diff --git a/scribo/ChangeLog b/scribo/ChangeLog index 1f87892..b5efaa7 100644 --- a/scribo/ChangeLog +++ b/scribo/ChangeLog @@ -9549,10 +9549,6 @@ * tests/unit_test/build_unit_test.sh: new. -$Rev$ -$Id$ -$Date$ - Local Variables: ispell-local-dictionary: "american" coding: iso-latin-1 diff --git a/swilena/ChangeLog b/swilena/ChangeLog index 28d78ca..292b4c7 100644 --- a/swilena/ChangeLog +++ b/swilena/ChangeLog @@ -903,10 +903,6 @@ * python/Makefile.am: New. -$Rev$ -$Id$ -$Date$ - Local Variables: ispell-local-dictionary: "american" coding: iso-latin-1 -- 1.7.10.4

11 years, 4 months

olena: olena-2.0-641-gf15b2c0 Improve the generation of Doxyfiles in Olena.

by Roland Levillain

* doc/Makefile.am (edit) ($(srcdir)/user-refman.stamp): Use `[@]' instead of `\@' for portability reasons. (ID): New (empty). ($(DOXYFILE_USER)): Do not depend on $Id$ from ChangeLog. * doc/Doxyfile.in (PROJECT_NUMBER): Reintegrate @ID@ to the project number, even if it is not yet set. Improve the generation of Doxyfiles in Olena. --- ChangeLog | 13 +++++++++++++ doc/Doxyfile.in | 4 ++-- doc/Makefile.am | 39 ++++++++++++++++++++------------------- 3 files changed, 35 insertions(+), 21 deletions(-) diff --git a/ChangeLog b/ChangeLog index dd4134d..0447440 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,16 @@ +2013-08-26 Roland Levillain <roland(a)lrde.epita.fr> + + Improve the generation of Doxyfiles in Olena. + + * doc/Makefile.am (edit) + ($(srcdir)/user-refman.stamp): + Use `[@]' instead of `\@' for portability reasons. + (ID): New (empty). + ($(DOXYFILE_USER)): Do not depend on $Id$ from ChangeLog. + * doc/Doxyfile.in (PROJECT_NUMBER): Reintegrate @ID@ to the + project number, even if it is not yet set. + Improve the generation of Doxyfiles in Olena. + 2012-10-09 Roland Levillain <roland(a)lrde.epita.fr> Ask configure to try to find FFTW (version 3). diff --git a/doc/Doxyfile.in b/doc/Doxyfile.in index 6cb8bfb..526a24e 100644 --- a/doc/Doxyfile.in +++ b/doc/Doxyfile.in @@ -1,4 +1,4 @@ -# Copyright (C) 2011, 2012 EPITA Research and Development Laboratory +# Copyright (C) 2011, 2012, 2013 EPITA Research and Development Laboratory # (LRDE) # # This file is part of Olena. @@ -21,7 +21,7 @@ # Project related configuration options #--------------------------------------------------------------------------- PROJECT_NAME = "Olena" -PROJECT_NUMBER = "User documentation @PACKAGE_VERSION@" +PROJECT_NUMBER = "User documentation @PACKAGE_VERSION@ @ID@" PROJECT_BRIEF = "An Image Processing Platform" PROJECT_LOGO = @top_srcdir@/doc/images/logo.jpg OUTPUT_DIRECTORY = @builddir@/user-refman.tmp diff --git a/doc/Makefile.am b/doc/Makefile.am index 2ecb26d..878f7ce 100644 --- a/doc/Makefile.am +++ b/doc/Makefile.am @@ -86,25 +86,26 @@ REFMAN_deps = $(srcdir)/footer.html \ EXTRA_DIST += $(REFMAN_deps) -# Sed is used to generate Doxyfile from Doxyfile.in instead of -# configure, because the former is way faster than the latter. -# -## FIXME: This is because, as in TC, we depend on $Id$ from the -## ChangeLog. Maybe we should depend from something less prone to -## change. See Bison's build system. -edit = sed -e "s|@ID@|$$Id|" \ - -e 's,@PACKAGE_NAME\@,$(PACKAGE_NAME),g' \ - -e 's,@PACKAGE_VERSION\@,$(PACKAGE_VERSION),g' \ - -e 's,@top_srcdir\@,$(top_srcdir),g' \ - -e 's,@abs_top_srcdir\@,$(abs_top_srcdir),g' \ - -e 's,@top_builddir\@,$(top_builddir),g' \ - -e 's,@srcdir\@,$(srcdir),g' \ - -e 's,@builddir\@,$(builddir),g' +## FIXME: We used to depend on Subversion's $Id$ (extracted from +## ChangeLog) before. We should depend on something less prone to +## change from now on. The best is probably to use `git-version-gen' +## from Gnulib (as Bison does). For the moment, set ID to an empty +## string. +ID = +# Use Make to generate Doxyfile instead of `configure', as advised by +# Autoconf's manual (see section ``Installation Directory +# Variables''). +edit = sed -e "s|@ID[@]|$(ID)|" \ + -e 's,@PACKAGE_NAME[@],$(PACKAGE_NAME),g' \ + -e 's,@PACKAGE_VERSION[@],$(PACKAGE_VERSION),g' \ + -e 's,@top_srcdir[@],$(top_srcdir),g' \ + -e 's,@abs_top_srcdir[@],$(abs_top_srcdir),g' \ + -e 's,@top_builddir[@],$(top_builddir),g' \ + -e 's,@srcdir[@],$(srcdir),g' \ + -e 's,@builddir[@],$(builddir),g' $(DOXYFILE_USER): $(srcdir)/$(DOXYFILE).in - Id=`grep '^\$$Id' $(top_srcdir)/milena/ChangeLog \ - | sed -e 's/\\\$$//g'`; \ - $(edit) $< >$@ + $(edit) $< >$@ MAINTAINERCLEANFILES = $(DOXYFILE_USER) USER_REFMAN = user-refman @@ -122,8 +123,8 @@ $(srcdir)/user-refman.stamp: $(srcdir)/$(DOXYFILE).in $(REFMAN_deps) rm -rf $(srcdir)/$(USER_REFMAN) ## Fix path to the main page of the documentation. find $(USER_REFMAN).tmp/html -name '*.html' -exec \ - perl -pi -e "s,\@relcommondocpath\@,../../,g;" \ - -e "s,\@userrefman\@,$(USER_REFMAN),g" \ + perl -pi -e "s,@relcommondocpath[@],../../,g;" \ + -e "s,@userrefman[@],$(USER_REFMAN),g" \ {} \; ## FIXME: Moving directories between file systems is not portable. mv $(USER_REFMAN).tmp $(srcdir)/$(USER_REFMAN) -- 1.7.10.4

11 years, 4 months

olena: olena-2.0-640-g8f094c6 Improve the generation of Doxyfiles in Scribo.

by Roland Levillain

* doc/Makefile.am (edit) ($(srcdir)/user-refman.stamp) ($(srcdir)/devel-refman.stamp): Use `[@]' instead of `\@' for portability reasons. (ID): New (empty). (edit_devel): Prevent configure from expanding @builddir@. Christen the product as ``Developer documentation''. ($(DOXYFILE_USER), $(DOXYFILE_DEVEL)): Do not depend on $Id$ from ChangeLog. * doc/Doxyfile.in (PROJECT_NUMBER): Reintegrate @ID@ to the project number, even if it is not yet set. --- scribo/ChangeLog | 16 ++++++++++++++++ scribo/doc/Doxyfile.in | 2 +- scribo/doc/Makefile.am | 50 ++++++++++++++++++++++++------------------------ 3 files changed, 42 insertions(+), 26 deletions(-) diff --git a/scribo/ChangeLog b/scribo/ChangeLog index 88f7ca5..1f87892 100644 --- a/scribo/ChangeLog +++ b/scribo/ChangeLog @@ -1,3 +1,19 @@ +2013-08-26 Roland Levillain <roland(a)lrde.epita.fr> + + Improve the generation of Doxyfiles in Scribo. + + * doc/Makefile.am (edit) + ($(srcdir)/user-refman.stamp) + ($(srcdir)/devel-refman.stamp): + Use `[@]' instead of `\@' for portability reasons. + (ID): New (empty). + (edit_devel): Prevent configure from expanding @builddir@. + Christen the product as ``Developer documentation''. + ($(DOXYFILE_USER), $(DOXYFILE_DEVEL)): Do not depend on $Id$ from + ChangeLog. + * doc/Doxyfile.in (PROJECT_NUMBER): Reintegrate @ID@ to the + project number, even if it is not yet set. + 2013-08-22 Roland Levillain <roland(a)lrde.epita.fr> Slightly improve Scribo's tests/primitive/extract/alignments.cc. diff --git a/scribo/doc/Doxyfile.in b/scribo/doc/Doxyfile.in index 904121c..d37780e 100644 --- a/scribo/doc/Doxyfile.in +++ b/scribo/doc/Doxyfile.in @@ -21,7 +21,7 @@ # Project related configuration options #--------------------------------------------------------------------------- PROJECT_NAME = "Scribo" -PROJECT_NUMBER = "User documentation @PACKAGE_VERSION@" +PROJECT_NUMBER = "User documentation @PACKAGE_VERSION@ @ID@" PROJECT_BRIEF = "An Image Processing Platform" PROJECT_LOGO = @top_srcdir@/doc/images/logo.jpg OUTPUT_DIRECTORY = @builddir@/user-refman.tmp diff --git a/scribo/doc/Makefile.am b/scribo/doc/Makefile.am index c459d84..5264ddf 100644 --- a/scribo/doc/Makefile.am +++ b/scribo/doc/Makefile.am @@ -54,25 +54,26 @@ regen-doc: DOXYFILE_USER = Doxyfile_user -# Sed is used to generate Doxyfile from Doxyfile.in instead of -# configure, because the former is way faster than the latter. -# -## FIXME: This is because, as in TC, we depend on $Id$ from the -## ChangeLog. Maybe we should depend from something less prone to -## change. See Bison's build system. -edit = sed -e "s|@ID@|$$Id|" \ - -e 's,@PACKAGE_NAME\@,$(PACKAGE_NAME),g' \ - -e 's,@PACKAGE_VERSION\@,$(PACKAGE_VERSION),g' \ - -e 's,@top_srcdir\@,$(top_srcdir),g' \ - -e 's,@abs_top_srcdir\@,$(abs_top_srcdir),g' \ - -e 's,@top_builddir\@,$(top_builddir),g' \ - -e 's,@srcdir\@,$(srcdir),g' \ - -e 's,@builddir\@,$(builddir),g' +## FIXME: We used to depend on Subversion's $Id$ (extracted from +## ChangeLog) before. We should depend on something less prone to +## change from now on. The best is probably to use `git-version-gen' +## from Gnulib (as Bison does). For the moment, set ID to an empty +## string. +ID = +# Use Make to generate Doxyfile instead of `configure', as advised by +# Autoconf's manual (see section ``Installation Directory +# Variables''). +edit = sed -e "s|@ID[@]|$(ID)|" \ + -e 's,@PACKAGE_NAME[@],$(PACKAGE_NAME),g' \ + -e 's,@PACKAGE_VERSION[@],$(PACKAGE_VERSION),g' \ + -e 's,@top_srcdir[@],$(top_srcdir),g' \ + -e 's,@abs_top_srcdir[@],$(abs_top_srcdir),g' \ + -e 's,@top_builddir[@],$(top_builddir),g' \ + -e 's,@srcdir[@],$(srcdir),g' \ + -e 's,@builddir[@],$(builddir),g' $(DOXYFILE_USER): $(srcdir)/$(DOXYFILE).in - Id=`grep '^\$$Id' $(top_srcdir)/scribo/ChangeLog \ - | sed -e 's/\\\$$//g'`; \ - $(edit) $< >$@ + $(edit) $< >$@ MAINTAINERCLEANFILES = $(DOXYFILE_USER) USER_REFMAN = user-refman @@ -100,8 +101,8 @@ $(srcdir)/user-refman.stamp: $(srcdir)/$(DOXYFILE).in $(REFMAN_deps) rm -rf $(srcdir)/$(USER_REFMAN) ## Fix path to the main page of the documentation. find $(USER_REFMAN).tmp -name '*.html' -exec \ - perl -pi -e "s,\@relcommondocpath\@,../../../../doc/,g;" \ - -e "s,\@userrefman\@,$(USER_REFMAN),g" \ + perl -pi -e "s,@relcommondocpath[@],../../../../doc/,g;" \ + -e "s,@userrefman[@],$(USER_REFMAN),g" \ {} \; ## FIXME: Moving directories between file systems is not portable. mv $(USER_REFMAN).tmp $(srcdir)/$(USER_REFMAN) @@ -136,7 +137,7 @@ regen-devel-doc: DOXYFILE_DEVEL = Doxyfile_devel edit_devel = sed \ - -e 's,OUTPUT_DIRECTORY = @builddir@/user-refman.tmp,OUTPUT_DIRECTORY = @builddir@/devel-refman.tmp,' \ + -e 's,OUTPUT_DIRECTORY = $.*$/user-refman.tmp,OUTPUT_DIRECTORY = \1/devel-refman.tmp,' \ -e 's,EXTRACT_ALL = NO,EXTRACT_ALL = YES,' \ -e 's,EXTRACT_PRIVATE = NO,EXTRACT_PRIVATE = YES,' \ -e 's,EXTRACT_STATIC = NO,EXTRACT_STATIC = YES,' \ @@ -148,13 +149,12 @@ edit_devel = sed \ -e 's,INTERNAL_DOCS = NO,INTERNAL_DOCS = YES,' \ -e 's,GENERATE_TODOLIST = NO,GENERATE_TODOLIST = YES,' \ -e 's,PROJECT_NUMBER = \"User documentation ,PROJECT_NUMBER = \",' \ + -e 's,PROJECT_NUMBER = \"User documentation ,PROJECT_NUMBER = \"Developer documentation ,' \ -e 's,EXCLUDE_SYMBOLS = \*::internal\* \*_ mln::trait::\*,EXCLUDE_SYMBOLS =,' $(DOXYFILE_DEVEL): $(srcdir)/$(DOXYFILE).in rm -f $@.tmp - Id=`grep '^\$$Id' $(top_srcdir)/scribo/ChangeLog \ - | sed -e 's/\\\$$//g'`; \ - $(edit) $< >$@.tmp + $(edit) $< >$@.tmp $(edit_devel) $@.tmp >$@ rm -f $@.tmp MAINTAINERCLEANFILES += $(DOXYFILE_DEVEL) @@ -183,8 +183,8 @@ $(srcdir)/devel-refman.stamp: $(srcdir)/$(DOXYFILE).in $(REFMAN_deps) test -f $(DEVEL_REFMAN).tmp/latex/refman.tex rm -rf $(srcdir)/$(DEVEL_REFMAN) ## Fix path to the main page of the documentation. - perl -pi -e "s,\@relcommondocpath\@,../../../../doc/,g;" \ - -e "s,\@userrefman\@,$(DEVEL_REFMAN),g" \ + perl -pi -e "s,@relcommondocpath[@],../../../../doc/,g;" \ + -e "s,@userrefman[@],$(DEVEL_REFMAN),g" \ $(DEVEL_REFMAN).tmp/html/index.html ## FIXME: Moving directories between file systems is not portable. mv $(DEVEL_REFMAN).tmp $(srcdir)/$(DEVEL_REFMAN) -- 1.7.10.4

11 years, 4 months

olena: olena-2.0-637-g8b31df6 Split the interface of the FFT off from its implementation.

by Roland Levillain

* mln/transform/fft.hh: Here. --- milena/ChangeLog | 6 + milena/mln/transform/fft.hh | 722 ++++++++++++++++++++++++++----------------- 2 files changed, 436 insertions(+), 292 deletions(-) diff --git a/milena/ChangeLog b/milena/ChangeLog index 8f96157..72bfb1a 100644 --- a/milena/ChangeLog +++ b/milena/ChangeLog @@ -1,5 +1,11 @@ 2012-10-10 Roland Levillain <roland(a)lrde.epita.fr> + Split the interface of the FFT off from its implementation. + + * mln/transform/fft.hh: Here. + +2012-10-10 Roland Levillain <roland(a)lrde.epita.fr> + Stylistic changes in the Fast Fourier Transform. * mln/transform/fft.hh: Here. diff --git a/milena/mln/transform/fft.hh b/milena/mln/transform/fft.hh index 9b70b4f..36e4b80 100644 --- a/milena/mln/transform/fft.hh +++ b/milena/mln/transform/fft.hh @@ -42,9 +42,11 @@ namespace mln namespace internal { - /// Dispatch traits for fftw + /// Dispatch tags for FFT traits. enum fft_dispatch { fft_cplx, fft_real }; + + /// \brief FFT traits. template <typename T> struct fft_trait; @@ -68,6 +70,7 @@ namespace mln typedef std::complex<T> fftw_input; }; + /// \brief Internal structure containining packing data and /// instructions for the (forward) and inverse (backward) /// transforms. @@ -78,16 +81,10 @@ namespace mln { public: /// \brief Accessor to transformed image (const version). - const image2d< std::complex<T> >& transformed_image() const - { - return trans_im; - } + const image2d< std::complex<T> >& transformed_image() const; /// \brief Accessor to transformed image (non const version). - image2d< std::complex<T> >& transformed_image() - { - return trans_im; - } + image2d< std::complex<T> >& transformed_image(); /// \brief Accessor to the transformed image (magnitude). /// @@ -98,28 +95,7 @@ namespace mln /// /// \param ordered Kind of traversal. template <class R> - image2d<R> transformed_image_magn(bool ordered = true) const - { - // FIXME: Ensure R is real. - - const unsigned nrows = trans_im.nrows(); - const unsigned ncols = trans_im.ncols(); - image2d<R> new_im(trans_im.domain()); - if (ordered) - for (unsigned row = 0; row < new_im.nrows(); ++row) - for (unsigned col = 0; col < new_im.ncols(); ++col) - opt::at(new_im, row, col) = - std::norm(opt::at(trans_im, - (row + nrows / 2) % nrows, - (col + ncols / 2) % ncols)); - else - { - mln_piter(image2d< std::complex<T> >) it(trans_im.domain()); - for_all(it) - new_im(it) = std::norm(trans_im(it)); - } - return new_im; - } + image2d<R> transformed_image_magn(bool ordered = true) const; /// \brief Accessor to the transformed image (magnitude). /// @@ -127,10 +103,7 @@ namespace mln /// an image containing <tt>|T[p]|</tt> is returned. /// /// \param ordered Kind of traversal. - image2d<T> transformed_image_magn(bool ordered = true) const - { - return transformed_image_magn<T>(ordered); - } + image2d<T> transformed_image_magn(bool ordered = true) const; /// \brief Accessor to the transformed image (clipped /// magnitude). @@ -144,43 +117,7 @@ namespace mln /// \param ordered Kind of traversal. template <class R> image2d<R> transformed_image_clipped_magn(double clip, - bool ordered = true) const - { - // FIXME: Ensure R is real. - // FIXME: Ensure `clip' is between 0 and 1? - - double max = mln_min(double); - mln_piter(image2d<T>) it(trans_im.domain()); - for_all(it) - if (std::norm(trans_im(it)) > max) - max = std::norm(trans_im(it)); - - const unsigned nrows = trans_im.nrows(); - const unsigned ncols = trans_im.ncols(); - image2d<R> new_im(trans_im.domain()); - if (ordered) - for (unsigned row = 0; row < new_im.nrows(); ++row) - for (unsigned col = 0; col < new_im.ncols(); ++col) - if (std::norm(opt::at(trans_im, - (row + nrows / 2) % nrows, - (col + ncols / 2) % ncols)) - >= max * clip) - opt::at(new_im, row, col) = mln_max(R); - else - opt::at(new_im, row, col) = - (double) mln_max(R) * - std::norm(opt::at(trans_im, - (row + nrows / 2) % nrows, - (col + ncols / 2) % ncols)) / (max * clip); - else - for_all(it) - if (std::norm(trans_im(it)) >= max * clip) - new_im(it) = mln_max(R); - else - new_im(it) = - (double) mln_max(R) * std::norm(trans_im(it)) / (max * clip); - return new_im; - } + bool ordered = true) const; /// \brief Accessor to the transformed image (clipped /// magnitude). @@ -191,10 +128,7 @@ namespace mln /// \param clip Value used for clipping. /// \param ordered Kind of traversal. image2d<T> transformed_image_clipped_magn(double clip, - bool ordered = true) const - { - return transformed_image_clipped_magn<T>(clip, ordered); - } + bool ordered = true) const; /// \brief Accessor to the transformed image (clipped /// magnitude). @@ -206,10 +140,7 @@ namespace mln /// /// \param ordered Kind of traversal. template <class R> - image2d<R> transformed_image_clipped_magn(bool ordered = true) const - { - return transformed_image_clipped_magn<R>(1, ordered); - } + image2d<R> transformed_image_clipped_magn(bool ordered = true) const; /// \brief Accessor to the transformed image (clipped /// magnitude). @@ -218,10 +149,7 @@ namespace mln /// containing a clipped value of <tt>|T[p]|</tt> is returned. /// /// \param ordered Kind of traversal. - image2d<T> transformed_image_clipped_magn(bool ordered = true) const - { - return transformed_image_clipped_magn<T>(1, ordered); - } + image2d<T> transformed_image_clipped_magn(bool ordered = true) const; /// \brief Accessor to the transformed image (log of the /// magnitude). @@ -237,43 +165,7 @@ namespace mln /// \param ordered Kind of traversal. template <class R> image2d<R> transformed_image_log_magn(double a, double b, - bool ordered = true) const - { - /* FIXME: Find a more elegant way to fix range interval on a - and b (Note from Roland: what does it mean?). */ - - // FIXME: Ensure R is real. - // FIXME: Ensure 0 <= a <= 1000. - // FIXME: Ensure 0 <= b <= 1000. - - image2d<R> new_im(trans_im.domain()); - - double max = mln_min(double); - mln_piter(image2d<R>) it(trans_im.domain()); - for_all(it) - if (std::norm(trans_im(it)) > max) - max = std::norm(trans_im(it)); - - const unsigned nrows = trans_im.nrows(); - const unsigned ncols = trans_im.ncols(); - if (ordered) - for (unsigned row = 0; row < new_im.nrows(); ++row) - for (unsigned col = 0; col < new_im.ncols(); ++col) - opt::at(new_im, row, col) = - log(a + b * std::norm(opt::at(trans_im, - (row + nrows / 2) % nrows, - (col + ncols / 2) % ncols))) - / log (a + b * max) * mln_max(R); - else - { - mln_piter(image2d< std::complex<T> >) it(trans_im.domain()); - for_all(it) - new_im(it) = - log(a + b * std::norm(trans_im(it))) - / log (a + b * max) * mln_max(R); - } - return new_im; - } + bool ordered = true) const; /// \brief Accessor to the transformed image (log of the /// magnitude). @@ -286,10 +178,7 @@ namespace mln /// \param b Upper bound. /// \param ordered Kind of traversal. image2d<T> transformed_image_log_magn(double a, double b, - bool ordered = true) const - { - return transformed_image_log_magn<T>(a, b, ordered); - } + bool ordered = true) const; /// \brief Accessor to the transformed image (log of the /// magnitude). @@ -302,10 +191,7 @@ namespace mln /// /// \param ordered Kind of traversal. template <class R> - image2d<R> transformed_image_log_magn(bool ordered = true) const - { - return transformed_image_log_magn<R>(1, 100, ordered); - } + image2d<R> transformed_image_log_magn(bool ordered = true) const; /// \brief Accessor to the transformed image (log of the /// magnitude). @@ -315,18 +201,9 @@ namespace mln /// interval [1, 100]. /// /// \param ordered Kind of traversal. - image2d<T> transformed_image_log_magn(bool ordered = true) const - { - return transformed_image_log_magn<T>(1, 100, ordered); - } - - ~fft() - { - fftw_free(in); - fftw_free(out); - fftw_destroy_plan(p); - fftw_destroy_plan(p_inv); - } + image2d<T> transformed_image_log_magn(bool ordered = true) const; + + ~fft(); protected: /// Input image (used internally by FFTW). @@ -359,7 +236,7 @@ namespace mln class fft; - /// \brief FFT engine (specialization for images of real values). + /// \brief oFFT engine (specialization for images of real values). /// /// \tparam T Data type. template <class T> @@ -372,96 +249,20 @@ namespace mln /// /// \param original_im Image to process. template <typename D> - fft(const image2d<D>& original_im) - { - const unsigned nrows = original_im.nrows(); - const unsigned ncols = original_im.ncols(); - - this->in = (T*) fftw_malloc(nrows * ncols * sizeof(T)); - this->out = (std::complex<T>*) fftw_malloc(nrows - * (ncols / 2 + 1) - * sizeof(std::complex<T>)); - - for (unsigned row = 0; row < nrows; ++row) - for (unsigned col = 0; col < ncols; ++col) - this->in[row * ncols + col] = - opt::at(original_im, row, col); - - this->p = - fftw_plan_dft_r2c_2d (nrows, ncols, - this->in, - reinterpret_cast<fftw_complex*>(this->out), - FFTW_ESTIMATE); - this->p_inv = - fftw_plan_dft_c2r_2d (nrows, ncols, - reinterpret_cast<fftw_complex*>(this->out), - this->in, - FFTW_ESTIMATE); - this->trans_im = image2d< std::complex<T> >(original_im.domain()); - } + fft(const image2d<D>& original_im); /// \brief Compute and return the transform (as a complex image). - image2d< std::complex<T> > transform() - { - fftw_execute(this->p); - - const unsigned nrows = this->trans_im.nrows(); - const unsigned ncols = this->trans_im.ncols(); - unsigned denom = nrows * ncols; - int i = 0; - for (unsigned row = 0; row < nrows; ++row) - for (unsigned col = 0; col <= ncols / 2; ++col) - { - this->out[i] = std::complex<T> (this->out[i].real() / denom, - this->out[i].imag() / denom); - opt::at(this->trans_im, row, col) = this->out[i]; - ++i; - } - for (unsigned row = 0; row < nrows; ++row) - for (unsigned col = ncols - 1; col > ncols / 2; --col) - opt::at(this->trans_im, row, col) = - opt::at(this->trans_im, nrows - row - 1, ncols - col - 1); - return this->trans_im; - } + image2d< std::complex<T> > transform(); /// \brief Compute and return the inverse transform (as a real /// image) of the FFT. /// /// \tparam R Value type of output image. template <class R> - image2d<R> transform_inv() - { - // FIXME: Ensure R is real. - - const unsigned nrows = this->trans_im.nrows(); - const unsigned ncols = this->trans_im.ncols(); - for (unsigned row = 0; row < nrows; ++row) - for (unsigned col = 0; col <= ncols / 2; ++col) - this->out[row * (ncols / 2 + 1) + col] = - opt::at(this->trans_im, row, col); - - fftw_execute(this->p_inv); - - image2d<R> new_im(this->trans_im.domain()); - int i = 0; - for (unsigned row = 0; row < nrows; ++row) - for (unsigned col = 0; col < ncols; ++col) - { - opt::at(new_im, row, col) = (this->in[i] >= mln_min(R) - ? (this->in[i] <= mln_max(R) - ? (R) this->in [i] - : mln_min(R)) - : mln_max(R)); - ++i; - } - return new_im; - } + image2d<R> transform_inv(); /// \brief Compute and return the inverse transform of the FFT. - image2d<T> transform_inv() - { - return transform_inv<T>(); - } + image2d<T> transform_inv(); }; @@ -477,48 +278,10 @@ namespace mln /// Initialization of data for the computation of the FFT. /// /// \param original_im Image to process. - fft(const image2d< std::complex<T> >& original_im) - { - const unsigned nrows = original_im.nrows(); - const unsigned ncols = original_im.ncols(); - - this->in = fftw_malloc(sizeof(std::complex<T>) * nrows * ncols); - this->out = fftw_malloc(sizeof(std::complex<T>) * nrows * ncols); - - for (unsigned row = 0; row < nrows; ++row) - for (unsigned col = 0; col < ncols; ++col) - { - this->in[row * ncols + col].re = original_im(row, col).real(); - this->in[row * ncols + col].im = original_im(row, col).imag(); - } - - this->p = fftw_plan_dft_2d(nrows, ncols, this->in, this->out, - FFTW_FORWARD, FFTW_ESTIMATE); - this->p_inv = fftw_plan_dft_2d(nrows, ncols, this->out, this->in, - FFTW_BACKWARD, FFTW_ESTIMATE); - this->trans_im = image2d< std::complex<T> >(original_im.domain()); - } + fft(const image2d< std::complex<T> >& original_im); /// \brief Compute and return the transform (as a complex image). - image2d< std::complex<T> > transform() - { - fftw_execute(this->p); - - const unsigned nrows = this->trans_im.nrows(); - const unsigned ncols = this->trans_im.ncols(); - unsigned denom = nrows * ncols; - int i = 0; - for (unsigned row = 0; row < nrows; ++row) - for (unsigned col = 0; col < ncols; ++col) - { - this->out[i].re /= denom; - this->out[i].im /= denom; - this->trans_im(row, col) = - std::complex<double>(this->out[i].re, this->out[i].im); - ++i; - } - return this->trans_im; - } + image2d< std::complex<T> > transform(); /// \brief Compute and return the inverse transform (as a /// complex image). @@ -526,42 +289,417 @@ namespace mln /// \tparam R Component type of the value type of the output /// image. template <class R> - image2d< std::complex<R> > transform_inv() - { - const unsigned nrows = this->trans_im.nrows(); - const unsigned ncols = this->trans_im.ncols(); - for (unsigned row = 0; row < nrows; ++row) - for (unsigned col = 0; col < ncols; ++col) - { - this->out[row * ncols + col].re = - this->trans_im(row, col).real(); - this->out[row * ncols + col].im = - this->trans_im(row, col).imag(); - } - - fftw_execute(this->p_inv); - - image2d< std::complex<R> > new_im(this->trans_im.size()); - int i = 0; - for (unsigned row = 0; row < nrows; ++row) - for (unsigned col = 0; col < ncols; ++col) - { - new_im(row, col) = this->in[i]; - ++i; - } - return new_im; - } + image2d< std::complex<R> > transform_inv(); /// \brief Compute and return the inverse transform of the FFT. - image2d< std::complex<T> > transform_inv() - { - return transform_inv<T>(); - } - + image2d< std::complex<T> > transform_inv(); }; } // end of namespace mln::transform + + +# ifndef MLN_INCLUDE_ONLY + + namespace internal + { + + template <class T> + inline + const image2d< std::complex<T> >& + fft<T>::transformed_image() const + { + return trans_im; + } + + template <class T> + inline + image2d< std::complex<T> >& + fft<T>::transformed_image() + { + return trans_im; + } + + template <class T> + template <class R> + inline + image2d<R> + fft<T>::transformed_image_magn(bool ordered) const + { + // FIXME: Ensure R is real. + + const unsigned nrows = trans_im.nrows(); + const unsigned ncols = trans_im.ncols(); + image2d<R> new_im(trans_im.domain()); + if (ordered) + for (unsigned row = 0; row < new_im.nrows(); ++row) + for (unsigned col = 0; col < new_im.ncols(); ++col) + opt::at(new_im, row, col) = + std::norm(opt::at(trans_im, + (row + nrows / 2) % nrows, + (col + ncols / 2) % ncols)); + else + { + mln_piter(image2d< std::complex<T> >) it(trans_im.domain()); + for_all(it) + new_im(it) = std::norm(trans_im(it)); + } + return new_im; + } + + template <class T> + inline + image2d<T> + fft<T>::transformed_image_magn(bool ordered) const + { + return transformed_image_magn<T>(ordered); + } + + template <class T> + template <class R> + inline + image2d<R> + fft<T>::transformed_image_clipped_magn(double clip, bool ordered) const + { + // FIXME: Ensure R is real. + // FIXME: Ensure `clip' is between 0 and 1? + + double max = mln_min(double); + mln_piter(image2d<T>) it(trans_im.domain()); + for_all(it) + if (std::norm(trans_im(it)) > max) + max = std::norm(trans_im(it)); + + const unsigned nrows = trans_im.nrows(); + const unsigned ncols = trans_im.ncols(); + image2d<R> new_im(trans_im.domain()); + if (ordered) + for (unsigned row = 0; row < new_im.nrows(); ++row) + for (unsigned col = 0; col < new_im.ncols(); ++col) + if (std::norm(opt::at(trans_im, + (row + nrows / 2) % nrows, + (col + ncols / 2) % ncols)) + >= max * clip) + opt::at(new_im, row, col) = mln_max(R); + else + opt::at(new_im, row, col) = + (double) mln_max(R) * + std::norm(opt::at(trans_im, + (row + nrows / 2) % nrows, + (col + ncols / 2) % ncols)) / (max * clip); + else + for_all(it) + if (std::norm(trans_im(it)) >= max * clip) + new_im(it) = mln_max(R); + else + new_im(it) = + (double) mln_max(R) * std::norm(trans_im(it)) / (max * clip); + return new_im; + } + + template <class T> + inline + image2d<T> + fft<T>::transformed_image_clipped_magn(double clip, bool ordered) const + { + return transformed_image_clipped_magn<T>(clip, ordered); + } + + template <class T> + template <class R> + inline + image2d<R> + fft<T>::transformed_image_clipped_magn(bool ordered) const + { + return transformed_image_clipped_magn<R>(1, ordered); + } + + template <class T> + inline + image2d<T> + fft<T>::transformed_image_clipped_magn(bool ordered) const + { + return transformed_image_clipped_magn<T>(1, ordered); + } + + template <class T> + template <class R> + inline + image2d<R> + fft<T>::transformed_image_log_magn(double a, double b, bool ordered) const + { + /* FIXME: Find a more elegant way to fix range interval on a + and b (Note from Roland: what does it mean?). */ + + // FIXME: Ensure R is real. + // FIXME: Ensure 0 <= a <= 1000. + // FIXME: Ensure 0 <= b <= 1000. + + image2d<R> new_im(trans_im.domain()); + + double max = mln_min(double); + mln_piter(image2d<R>) it(trans_im.domain()); + for_all(it) + if (std::norm(trans_im(it)) > max) + max = std::norm(trans_im(it)); + + const unsigned nrows = trans_im.nrows(); + const unsigned ncols = trans_im.ncols(); + if (ordered) + for (unsigned row = 0; row < new_im.nrows(); ++row) + for (unsigned col = 0; col < new_im.ncols(); ++col) + opt::at(new_im, row, col) = + log(a + b * std::norm(opt::at(trans_im, + (row + nrows / 2) % nrows, + (col + ncols / 2) % ncols))) + / log (a + b * max) * mln_max(R); + else + { + mln_piter(image2d< std::complex<T> >) it(trans_im.domain()); + for_all(it) + new_im(it) = + log(a + b * std::norm(trans_im(it))) + / log (a + b * max) * mln_max(R); + } + return new_im; + } + + template <class T> + inline + image2d<T> + fft<T>::transformed_image_log_magn(double a, double b, bool ordered) const + { + return transformed_image_log_magn<T>(a, b, ordered); + } + + template <class T> + template <class R> + inline + image2d<R> + fft<T>::transformed_image_log_magn(bool ordered) const + { + return transformed_image_log_magn<R>(1, 100, ordered); + } + + template <class T> + inline + image2d<T> + fft<T>::transformed_image_log_magn(bool ordered) const + { + return transformed_image_log_magn<T>(1, 100, ordered); + } + + template <class T> + inline + fft<T>::~fft() + { + fftw_free(in); + fftw_free(out); + fftw_destroy_plan(p); + fftw_destroy_plan(p_inv); + } + + } // end of namespace mln::internal + + + namespace transform + { + + /*--------------------------------------------------------. + | FFT engine (specialization for images of real values). | + `--------------------------------------------------------*/ + + template <class T> + template <typename D> + inline + fft<T, internal::fft_real>::fft(const image2d<D>& original_im) + { + const unsigned nrows = original_im.nrows(); + const unsigned ncols = original_im.ncols(); + + this->in = (T*) fftw_malloc(nrows * ncols * sizeof(T)); + this->out = (std::complex<T>*) fftw_malloc(nrows + * (ncols / 2 + 1) + * sizeof(std::complex<T>)); + + for (unsigned row = 0; row < nrows; ++row) + for (unsigned col = 0; col < ncols; ++col) + this->in[row * ncols + col] = + opt::at(original_im, row, col); + + this->p = + fftw_plan_dft_r2c_2d (nrows, ncols, + this->in, + reinterpret_cast<fftw_complex*>(this->out), + FFTW_ESTIMATE); + this->p_inv = + fftw_plan_dft_c2r_2d (nrows, ncols, + reinterpret_cast<fftw_complex*>(this->out), + this->in, + FFTW_ESTIMATE); + this->trans_im = image2d< std::complex<T> >(original_im.domain()); + } + + template <class T> + inline + image2d< std::complex<T> > + fft<T, internal::fft_real>::transform() + { + fftw_execute(this->p); + + const unsigned nrows = this->trans_im.nrows(); + const unsigned ncols = this->trans_im.ncols(); + unsigned denom = nrows * ncols; + int i = 0; + for (unsigned row = 0; row < nrows; ++row) + for (unsigned col = 0; col <= ncols / 2; ++col) + { + this->out[i] = std::complex<T> (this->out[i].real() / denom, + this->out[i].imag() / denom); + opt::at(this->trans_im, row, col) = this->out[i]; + ++i; + } + for (unsigned row = 0; row < nrows; ++row) + for (unsigned col = ncols - 1; col > ncols / 2; --col) + opt::at(this->trans_im, row, col) = + opt::at(this->trans_im, nrows - row - 1, ncols - col - 1); + return this->trans_im; + } + + template <class T> + template <class R> + inline + image2d<R> + fft<T, internal::fft_real>::transform_inv() + { + // FIXME: Ensure R is real. + + const unsigned nrows = this->trans_im.nrows(); + const unsigned ncols = this->trans_im.ncols(); + for (unsigned row = 0; row < nrows; ++row) + for (unsigned col = 0; col <= ncols / 2; ++col) + this->out[row * (ncols / 2 + 1) + col] = + opt::at(this->trans_im, row, col); + + fftw_execute(this->p_inv); + + image2d<R> new_im(this->trans_im.domain()); + int i = 0; + for (unsigned row = 0; row < nrows; ++row) + for (unsigned col = 0; col < ncols; ++col) + { + opt::at(new_im, row, col) = (this->in[i] >= mln_min(R) + ? (this->in[i] <= mln_max(R) + ? (R) this->in [i] + : mln_min(R)) + : mln_max(R)); + ++i; + } + return new_im; + } + + template <class T> + inline + image2d<T> + fft<T, internal::fft_real>::transform_inv() + { + return transform_inv<T>(); + } + + + /*-----------------------------------------------------------. + | FFT engine (specialization for images of complex values). | + `-----------------------------------------------------------*/ + + template <class T> + inline + fft<T, internal::fft_cplx>::fft(const image2d< std::complex<T> >& original_im) + { + const unsigned nrows = original_im.nrows(); + const unsigned ncols = original_im.ncols(); + + this->in = fftw_malloc(sizeof(std::complex<T>) * nrows * ncols); + this->out = fftw_malloc(sizeof(std::complex<T>) * nrows * ncols); + + for (unsigned row = 0; row < nrows; ++row) + for (unsigned col = 0; col < ncols; ++col) + { + this->in[row * ncols + col].re = original_im(row, col).real(); + this->in[row * ncols + col].im = original_im(row, col).imag(); + } + + this->p = fftw_plan_dft_2d(nrows, ncols, this->in, this->out, + FFTW_FORWARD, FFTW_ESTIMATE); + this->p_inv = fftw_plan_dft_2d(nrows, ncols, this->out, this->in, + FFTW_BACKWARD, FFTW_ESTIMATE); + this->trans_im = image2d< std::complex<T> >(original_im.domain()); + } + + template <class T> + inline + image2d< std::complex<T> > + fft<T, internal::fft_cplx>::transform() + { + fftw_execute(this->p); + + const unsigned nrows = this->trans_im.nrows(); + const unsigned ncols = this->trans_im.ncols(); + unsigned denom = nrows * ncols; + int i = 0; + for (unsigned row = 0; row < nrows; ++row) + for (unsigned col = 0; col < ncols; ++col) + { + this->out[i].re /= denom; + this->out[i].im /= denom; + this->trans_im(row, col) = + std::complex<double>(this->out[i].re, this->out[i].im); + ++i; + } + return this->trans_im; + } + + template <class T> + template <class R> + inline + image2d< std::complex<R> > + fft<T, internal::fft_cplx>::transform_inv() + { + const unsigned nrows = this->trans_im.nrows(); + const unsigned ncols = this->trans_im.ncols(); + for (unsigned row = 0; row < nrows; ++row) + for (unsigned col = 0; col < ncols; ++col) + { + this->out[row * ncols + col].re = + this->trans_im(row, col).real(); + this->out[row * ncols + col].im = + this->trans_im(row, col).imag(); + } + + fftw_execute(this->p_inv); + + image2d< std::complex<R> > new_im(this->trans_im.size()); + int i = 0; + for (unsigned row = 0; row < nrows; ++row) + for (unsigned col = 0; col < ncols; ++col) + { + new_im(row, col) = this->in[i]; + ++i; + } + return new_im; + } + + template <class T> + inline + image2d< std::complex<T> > + fft<T, internal::fft_cplx>::transform_inv() + { + return transform_inv<T>(); + } + + } // end of namespace mln::transform + +# endif // ! MLN_INCLUDE_ONLY + } // end of namespace mln #endif // ! MLN_TRANSFORM_FFT_HH -- 1.7.10.4

11 years, 4 months

olena: olena-2.0-636-g588a6d4 Stylistic changes in the Fast Fourier Transform.

by Roland Levillain

* mln/transform/fft.hh: Here. --- milena/ChangeLog | 6 + milena/mln/transform/fft.hh | 630 ++++++++++++++++++++----------------------- 2 files changed, 305 insertions(+), 331 deletions(-) diff --git a/milena/ChangeLog b/milena/ChangeLog index 6a5f5d7..8f96157 100644 --- a/milena/ChangeLog +++ b/milena/ChangeLog @@ -1,5 +1,11 @@ 2012-10-10 Roland Levillain <roland(a)lrde.epita.fr> + Stylistic changes in the Fast Fourier Transform. + + * mln/transform/fft.hh: Here. + +2012-10-10 Roland Levillain <roland(a)lrde.epita.fr> + Remove dead code in the Fast Fourier Transform. * mln/transform/fft.hh: Here. diff --git a/milena/mln/transform/fft.hh b/milena/mln/transform/fft.hh index dfbe1d7..9b70b4f 100644 --- a/milena/mln/transform/fft.hh +++ b/milena/mln/transform/fft.hh @@ -27,17 +27,20 @@ #ifndef MLN_TRANSFORM_FFT_HH # define MLN_TRANSFORM_FFT_HH -# include <mln/core/image/image2d.hh> -# include <mln/estim/min_max.hh> -# include <mln/opt/at.hh> +# include <complex> -# include <complex> +# include <fftw3.h> -# include <fftw3.h> +# include <mln/core/image/image2d.hh> +# include <mln/estim/min_max.hh> +# include <mln/opt/at.hh> -namespace mln { - namespace internal { +namespace mln +{ + + namespace internal + { /// Dispatch traits for fftw enum fft_dispatch { fft_cplx, fft_real }; @@ -45,316 +48,279 @@ namespace mln { template <typename T> struct fft_trait; - /*! - ** \brief fft trait. - */ + /// \brief FFT traits for `double'. template <> struct fft_trait<double> { - static const fft_dispatch which = fft_real; ///< Real dispatch. - typedef double fftw_input; ///< Type of input. + /// Type of dispatch (real numbers). + static const fft_dispatch which = fft_real; + /// Type of input (double-precision floating point numbers). + typedef double fftw_input; }; - /*! - ** \brief fft trait. - ** - ** \specialization for ntg::cplx<R, T> - */ + /// \brief FFT traits for `std::complex<T>'. template <typename T> struct fft_trait< std::complex<T> > { - static const fft_dispatch which = fft_cplx; ///< Complex dispatch. - typedef std::complex <T> fftw_input; ///< Type of input. + /// Type of dispatch (complex numbers). + static const fft_dispatch which = fft_cplx; + /// Type of input (complex numbers of the standard C++ library). + typedef std::complex<T> fftw_input; }; - /*! - ** _fft<ntg::cplx_representation, T> - ** - ** \param T Data type. - */ + /// \brief Internal structure containining packing data and + /// instructions for the (forward) and inverse (backward) + /// transforms. + /// + /// \tparam T Data type. template <class T> - class _fft + class fft { public: - /*! - ** \brief Accessor to transformed image. - ** - ** Const version. - */ + /// \brief Accessor to transformed image (const version). const image2d< std::complex<T> >& transformed_image() const { return trans_im; } - /*! - ** \brief Accessor to transformed image. - ** - ** Non const version. - */ + /// \brief Accessor to transformed image (non const version). image2d< std::complex<T> >& transformed_image() { return trans_im; } - /*! - ** \brief Accessor to transformed image. - ** - ** For each point p of the transformed image T, you get - ** |T[p]|. - ** - ** \tparam R Data type of the resulting image. - ** - ** \param ordered Kind of traversal. - */ + /// \brief Accessor to the transformed image (magnitude). + /// + /// For each point \p p of the transformed image \p T, + /// an image containing <tt>|T[p]|</tt> is returned. + /// + /// \tparam R Value type of the resulting image. + /// + /// \param ordered Kind of traversal. template <class R> image2d<R> transformed_image_magn(bool ordered = true) const { - // FIXME : check that R is real + // FIXME: Ensure R is real. + const unsigned nrows = trans_im.nrows(); + const unsigned ncols = trans_im.ncols(); image2d<R> new_im(trans_im.domain()); - if (ordered) for (unsigned row = 0; row < new_im.nrows(); ++row) for (unsigned col = 0; col < new_im.ncols(); ++col) opt::at(new_im, row, col) = - std::norm(opt::at(trans_im, (row + trans_im.nrows() / 2) % trans_im.nrows(), - (col + trans_im.ncols() / 2) % trans_im.ncols())); + std::norm(opt::at(trans_im, + (row + nrows / 2) % nrows, + (col + ncols / 2) % ncols)); else { mln_piter(image2d< std::complex<T> >) it(trans_im.domain()); - for_all(it) new_im(it) = std::norm(trans_im(it)); } - return new_im; } - /*! - ** \brief Accessor to transformed image. - ** - ** For each point p of the transformed image T, you get - ** |T[p]|. - ** - ** \param ordered Kind of traversal. - */ + /// \brief Accessor to the transformed image (magnitude). + /// + /// For each point \p p of the transformed image \p T, + /// an image containing <tt>|T[p]|</tt> is returned. + /// + /// \param ordered Kind of traversal. image2d<T> transformed_image_magn(bool ordered = true) const { return transformed_image_magn<T>(ordered); } - /*! - ** \brief Accessor to transformed image. - ** - ** For each point p of the transformed image T, you get - ** a clipped value of |T[p]|.\n - ** - ** \param R Data type of the resulting image. - ** - ** \param clip Value used for clipping. - ** \param ordered Kind of traversal. - */ - + /// \brief Accessor to the transformed image (clipped + /// magnitude). + /// + /// For each point \p p of the transformed image \p T, an image + /// containing a clipped value of <tt>|T[p]|</tt> is returned. + /// + /// \tparam R Value type of the resulting image. + /// + /// \param clip Value used for clipping. + /// \param ordered Kind of traversal. template <class R> image2d<R> transformed_image_clipped_magn(double clip, bool ordered = true) const { - // Check that R is real + // FIXME: Ensure R is real. + // FIXME: Ensure `clip' is between 0 and 1? - image2d<R> new_im(trans_im.domain()); - // check that clip is >=0 and <=1 ? double max = mln_min(double); mln_piter(image2d<T>) it(trans_im.domain()); - for_all(it) if (std::norm(trans_im(it)) > max) max = std::norm(trans_im(it)); + const unsigned nrows = trans_im.nrows(); + const unsigned ncols = trans_im.ncols(); + image2d<R> new_im(trans_im.domain()); if (ordered) for (unsigned row = 0; row < new_im.nrows(); ++row) for (unsigned col = 0; col < new_im.ncols(); ++col) - { - if (std::norm(opt::at(trans_im, (row + trans_im.nrows() / 2) % trans_im.nrows(), - (col + trans_im.ncols() / 2) % - trans_im.ncols())) >= max * clip) - opt::at(new_im, row, col) = mln_max(R); - else - opt::at(new_im, row, col) = - (double) mln_max(R) * - std::norm(opt::at(trans_im, (row + trans_im.nrows() / 2) % trans_im.nrows(), - (col + trans_im.ncols() / 2) % - trans_im.ncols())) / (max * clip); - } - else - { - for_all(it) - { - if (std::norm(trans_im(it)) >= max * clip) - new_im(it) = mln_max(R); + if (std::norm(opt::at(trans_im, + (row + nrows / 2) % nrows, + (col + ncols / 2) % ncols)) + >= max * clip) + opt::at(new_im, row, col) = mln_max(R); else - new_im(it) = (double) mln_max(R) * - std::norm(trans_im(it)) / (max * clip); - } - } - + opt::at(new_im, row, col) = + (double) mln_max(R) * + std::norm(opt::at(trans_im, + (row + nrows / 2) % nrows, + (col + ncols / 2) % ncols)) / (max * clip); + else + for_all(it) + if (std::norm(trans_im(it)) >= max * clip) + new_im(it) = mln_max(R); + else + new_im(it) = + (double) mln_max(R) * std::norm(trans_im(it)) / (max * clip); return new_im; } - /*! - ** \brief Accessor to transformed image. - ** - ** For each point p of the transformed image T, you get - ** a clipped value of |T[p]|.\n - ** - ** \param clip Value used for clipping. - ** \param ordered Kind of traversal. - */ + /// \brief Accessor to the transformed image (clipped + /// magnitude). + /// + /// For each point \p p of the transformed image \p T, an image + /// containing a clipped value of <tt>|T[p]|</tt> is returned. + /// + /// \param clip Value used for clipping. + /// \param ordered Kind of traversal. image2d<T> transformed_image_clipped_magn(double clip, bool ordered = true) const { return transformed_image_clipped_magn<T>(clip, ordered); } - /*! - ** \brief Accessor to transformed image. - ** - ** For each point p of the transformed image T, you get - ** a clipped value of |T[p]|.\n - ** - ** \param R Data type of the resulting image. - ** - ** \param ordered Kind of traversal. - */ - - image2d<T> transformed_image_clipped_magn(bool ordered = true) const - { - return transformed_image_clipped_magn<T>(1, ordered); - } - - /*! - ** \brief Accessor to transformed image. - ** - ** For each point p of the transformed image T, you get - ** a clipped value of |T[p]|.\n - ** - ** \param R Data type of the resulting image. - ** - ** \param ordered Kind of traversal. - */ - + /// \brief Accessor to the transformed image (clipped + /// magnitude). + /// + /// For each point \p p of the transformed image \p T, an image + /// containing a clipped value of <tt>|T[p]|</tt> is returned. + /// + /// \tparam R Value type of the resulting image. + /// + /// \param ordered Kind of traversal. template <class R> image2d<R> transformed_image_clipped_magn(bool ordered = true) const { return transformed_image_clipped_magn<R>(1, ordered); } - // FIXME: Find a more elegant way to fix range interval on a and b. - /*! - ** \brief Accessor to transformed image. - ** - ** For each point p of the transformed image T, you get - ** a log translated value of |T[p]| on interval [a; b].\n - ** - ** \param a Lower bound. - ** \param b Upper bound. - ** \param ordered Kind of traversal. - */ + /// \brief Accessor to the transformed image (clipped + /// magnitude). + /// + /// For each point \p p of the transformed image \p T, an image + /// containing a clipped value of <tt>|T[p]|</tt> is returned. + /// + /// \param ordered Kind of traversal. + image2d<T> transformed_image_clipped_magn(bool ordered = true) const + { + return transformed_image_clipped_magn<T>(1, ordered); + } + /// \brief Accessor to the transformed image (log of the + /// magnitude). + /// + /// For each point \p p of the transformed image \p T, an image + /// containing a log value of <tt>|T[p]|</tt> translated in the + /// interval [a, b]. + /// + /// \tparam R Value type of the resulting image. + /// + /// \param a Lower bound. + /// \param b Upper bound. + /// \param ordered Kind of traversal. template <class R> - image2d<R> transformed_image_log_magn(double a, - double b, + image2d<R> transformed_image_log_magn(double a, double b, bool ordered = true) const { - // Check that R is real - // 0 <= a <= 1000 - // 0 <= b <= 1000 + /* FIXME: Find a more elegant way to fix range interval on a + and b (Note from Roland: what does it mean?). */ + + // FIXME: Ensure R is real. + // FIXME: Ensure 0 <= a <= 1000. + // FIXME: Ensure 0 <= b <= 1000. image2d<R> new_im(trans_im.domain()); double max = mln_min(double); mln_piter(image2d<R>) it(trans_im.domain()); - for_all(it) if (std::norm(trans_im(it)) > max) max = std::norm(trans_im(it)); - + const unsigned nrows = trans_im.nrows(); + const unsigned ncols = trans_im.ncols(); if (ordered) for (unsigned row = 0; row < new_im.nrows(); ++row) for (unsigned col = 0; col < new_im.ncols(); ++col) opt::at(new_im, row, col) = - log(a + b * std::norm(opt::at(trans_im, (row + trans_im.nrows() / 2) % trans_im.nrows(), - (col + trans_im.ncols() / 2) % trans_im.ncols()))) / - log (a + b * max) * mln_max(R); + log(a + b * std::norm(opt::at(trans_im, + (row + nrows / 2) % nrows, + (col + ncols / 2) % ncols))) + / log (a + b * max) * mln_max(R); else { mln_piter(image2d< std::complex<T> >) it(trans_im.domain()); - for_all(it) - new_im(it) = log(a + b * std::norm(trans_im(it))) / - log (a + b * max) * mln_max(R); + new_im(it) = + log(a + b * std::norm(trans_im(it))) + / log (a + b * max) * mln_max(R); } - return new_im; } - // FIXME: Find a more elegant way to fix boundaries of a and b. - /*! - ** \brief Accessor to transformed image. - ** - ** For each point p of the transformed image T, you get - ** a log translated value of |T[p]| on interval [a; b].\n - ** - ** \param a Lower bound. - ** \param b Upper bound. - ** \param ordered Kind of traversal. - */ - - image2d<T> transformed_image_log_magn(double a, - double b, - bool ordered = true) const + /// \brief Accessor to the transformed image (log of the + /// magnitude). + /// + /// For each point \p p of the transformed image \p T, an image + /// containing a log value of <tt>|T[p]|</tt> translated in the + /// interval [a, b]. + /// + /// \param a Lower bound. + /// \param b Upper bound. + /// \param ordered Kind of traversal. + image2d<T> transformed_image_log_magn(double a, double b, + bool ordered = true) const { return transformed_image_log_magn<T>(a, b, ordered); } - /*! - ** \brief Accessor to transformed image. - ** - ** For each point p of the transformed image T, you get - ** a log translated value of |T[p]| on interval [1; 100].\n - ** - ** \param R Data type of the resulting image. - ** - ** \param ordered Kind of traversal. - */ - + /// \brief Accessor to the transformed image (log of the + /// magnitude). + /// + /// For each point \p p of the transformed image \p T, an image + /// containing a log value of <tt>|T[p]|</tt> translated in the + /// interval [1, 100]. + /// + /// \tparam R Value type of the resulting image. + /// + /// \param ordered Kind of traversal. template <class R> image2d<R> transformed_image_log_magn(bool ordered = true) const { return transformed_image_log_magn<R>(1, 100, ordered); } - /*! - ** \brief Accessor to transformed image. - ** - ** For each point p of the transformed image T, you get - ** a log translated value of |T[p]| on interval [1; 100].\n - ** - ** \param ordered Kind of traversal. - */ - + /// \brief Accessor to the transformed image (log of the + /// magnitude). + /// + /// For each point \p p of the transformed image \p T, an image + /// containing a log value of <tt>|T[p]|</tt> translated in the + /// interval [1, 100]. + /// + /// \param ordered Kind of traversal. image2d<T> transformed_image_log_magn(bool ordered = true) const { return transformed_image_log_magn<T>(1, 100, ordered); } - - /*! - ** \brief Destructor. - ** - ** Let memory free for other big images !!! - */ - ~_fft() + ~fft() { fftw_free(in); fftw_free(out); @@ -363,218 +329,222 @@ namespace mln { } protected: - - typename fft_trait<T>::fftw_input *in; ///< Input image. - std::complex<T> *out; ///< Complex image. - fftw_plan p; ///< Plan. - fftw_plan p_inv; ///< inverted plan. - image2d< std::complex<T> > trans_im; ///< Transformed image. - + /// Input image (used internally by FFTW). + typename fft_trait<T>::fftw_input* in; + /// Output (complex) image used internally by FFTW). + std::complex<T>* out; + /// (Forward) plan (command to be executed by FFTW). + fftw_plan p; + /// Inverse (backward) plan (command to be executed by FFTW). + fftw_plan p_inv; + /// Transformed image (``exported'' Milena image accessible + /// through accessors). + image2d< std::complex<T> > trans_im; }; - } // end of namespace internal + } // end of namespace mln::internal - namespace transform { - /*! - ** \brief fft class declaration. - ** - ** \param T Data type. - ** \param which Dispatch. - */ + namespace transform + { + + /// \brief FFT engine, producing a Discrete Fourier Transform + /// using FFTW (version 3). + /// + /// \tparam T Data type. + /// + /// \param which Dispatch traits. template <class T, internal::fft_dispatch which = internal::fft_trait<T>::which > class fft; - /*! - ** \brief fft specialization for fft real dispatch. - ** - ** \param T Data type. - */ + + /// \brief FFT engine (specialization for images of real values). + /// + /// \tparam T Data type. template <class T> - class fft<T, internal::fft_real> : public internal::_fft<T> + class fft<T, internal::fft_real> : public internal::fft<T> { - public: - - /*! - ** \brief Constructor. - ** - ** Initialization of data in order to compute the fft. - ** - ** \param original_im Image to process. - */ + /// \brief Constructor. + /// + /// Initialization of data for the computation of the FFT. + /// + /// \param original_im Image to process. template <typename D> fft(const image2d<D>& original_im) { - this->in = (T*) fftw_malloc(original_im.nrows() * original_im.ncols() * sizeof(T)); - this->out = (std::complex<T>*) - fftw_malloc(original_im.nrows() * (original_im.ncols() / 2 + 1) * sizeof(std::complex<T>)); - - for (unsigned row = 0; row < original_im.nrows(); ++row) - for (unsigned col = 0; col < original_im.ncols(); ++col) - this->in[row * original_im.ncols() + col] = opt::at(original_im, row, col); - - this->p = fftw_plan_dft_r2c_2d (original_im.nrows(), original_im.ncols(), - this->in, reinterpret_cast<fftw_complex*>(this->out), FFTW_ESTIMATE); - this->p_inv = fftw_plan_dft_c2r_2d (original_im.nrows(), original_im.ncols(), - reinterpret_cast<fftw_complex*>(this->out), this->in, FFTW_ESTIMATE); - + const unsigned nrows = original_im.nrows(); + const unsigned ncols = original_im.ncols(); + + this->in = (T*) fftw_malloc(nrows * ncols * sizeof(T)); + this->out = (std::complex<T>*) fftw_malloc(nrows + * (ncols / 2 + 1) + * sizeof(std::complex<T>)); + + for (unsigned row = 0; row < nrows; ++row) + for (unsigned col = 0; col < ncols; ++col) + this->in[row * ncols + col] = + opt::at(original_im, row, col); + + this->p = + fftw_plan_dft_r2c_2d (nrows, ncols, + this->in, + reinterpret_cast<fftw_complex*>(this->out), + FFTW_ESTIMATE); + this->p_inv = + fftw_plan_dft_c2r_2d (nrows, ncols, + reinterpret_cast<fftw_complex*>(this->out), + this->in, + FFTW_ESTIMATE); this->trans_im = image2d< std::complex<T> >(original_im.domain()); } - /*! - ** \brief Compute and return the transform. - */ + /// \brief Compute and return the transform (as a complex image). image2d< std::complex<T> > transform() { fftw_execute(this->p); - unsigned denom = this->trans_im.nrows() * this->trans_im.ncols(); + const unsigned nrows = this->trans_im.nrows(); + const unsigned ncols = this->trans_im.ncols(); + unsigned denom = nrows * ncols; int i = 0; - for (unsigned row = 0; row < this->trans_im.nrows(); ++row) - for (unsigned col = 0; col <= this->trans_im.ncols() / 2; ++col) + for (unsigned row = 0; row < nrows; ++row) + for (unsigned col = 0; col <= ncols / 2; ++col) { - this->out[i] = std::complex<T> (this->out[i].real() / denom, this->out[i].imag() / denom); + this->out[i] = std::complex<T> (this->out[i].real() / denom, + this->out[i].imag() / denom); opt::at(this->trans_im, row, col) = this->out[i]; ++i; } - for (unsigned row = 0; row < this->trans_im.nrows(); ++row) - for (unsigned col = this->trans_im.ncols() - 1; col > this->trans_im.ncols() / 2; --col) - opt::at(this->trans_im, row, col) = opt::at(this->trans_im, this->trans_im.nrows() - row - 1, - this->trans_im.ncols() - col - 1); + for (unsigned row = 0; row < nrows; ++row) + for (unsigned col = ncols - 1; col > ncols / 2; --col) + opt::at(this->trans_im, row, col) = + opt::at(this->trans_im, nrows - row - 1, ncols - col - 1); return this->trans_im; } - /*! - ** \brief Compute and return the invert transform. - ** - ** \param R Data type of output image. - */ + /// \brief Compute and return the inverse transform (as a real + /// image) of the FFT. + /// + /// \tparam R Value type of output image. template <class R> image2d<R> transform_inv() { - // FIXME : Check that R is real + // FIXME: Ensure R is real. - for (unsigned row = 0; row < this->trans_im.nrows(); ++row) - for (unsigned col = 0; col <= this->trans_im.ncols() / 2; ++col) - this->out[row * (this->trans_im.ncols() / 2 + 1) + col] = + const unsigned nrows = this->trans_im.nrows(); + const unsigned ncols = this->trans_im.ncols(); + for (unsigned row = 0; row < nrows; ++row) + for (unsigned col = 0; col <= ncols / 2; ++col) + this->out[row * (ncols / 2 + 1) + col] = opt::at(this->trans_im, row, col); fftw_execute(this->p_inv); image2d<R> new_im(this->trans_im.domain()); int i = 0; - for (unsigned row = 0; row < this->trans_im.nrows(); ++row) - for (unsigned col = 0; col < this->trans_im.ncols(); ++col) + for (unsigned row = 0; row < nrows; ++row) + for (unsigned col = 0; col < ncols; ++col) { - opt::at(new_im, row, col) = (this->in[i] >= mln_min(R) ? - (this->in[i] <= mln_max(R) ? - (R)this->in [i] : - mln_min(R)) : - mln_max(R)); + opt::at(new_im, row, col) = (this->in[i] >= mln_min(R) + ? (this->in[i] <= mln_max(R) + ? (R) this->in [i] + : mln_min(R)) + : mln_max(R)); ++i; } return new_im; } - /*! - ** \brief Compute and return the invert transform. - */ + /// \brief Compute and return the inverse transform of the FFT. image2d<T> transform_inv() { return transform_inv<T>(); } - }; - /*! - ** \brief fft specialization for fft complex dispatch. - ** - ** \param T Data type. - ** \param R Complex representation kind. - */ + + /// \brief FFT engine (specialization for images of complex values). + /// + /// \tparam T Data type. template <class T> - class fft<T, internal::fft_cplx> : public internal::_fft<T> + class fft<T, internal::fft_cplx> : public internal::fft<T> { - public: - - /*! - ** \brief Constructor. - ** - ** Initialization of data in order to compute the fft. - ** - ** \param original_im Image to process. - */ + /// \brief Constructor. + /// + /// Initialization of data for the computation of the FFT. + /// + /// \param original_im Image to process. fft(const image2d< std::complex<T> >& original_im) { - this->in = fftw_malloc(sizeof(std::complex<T>) * - original_im.nrows() * original_im.ncols()); - this->out = fftw_malloc(sizeof(std::complex<T>) * - original_im.nrows() * original_im.ncols()); + const unsigned nrows = original_im.nrows(); + const unsigned ncols = original_im.ncols(); + + this->in = fftw_malloc(sizeof(std::complex<T>) * nrows * ncols); + this->out = fftw_malloc(sizeof(std::complex<T>) * nrows * ncols); - for (unsigned row = 0; row < original_im.nrows(); ++row) - for (unsigned col = 0; col < original_im.ncols(); ++col) + for (unsigned row = 0; row < nrows; ++row) + for (unsigned col = 0; col < ncols; ++col) { - this->in[row * original_im.ncols() + col].re = - original_im(row, col).real(); - this->in[row * original_im.ncols() + col].im = - original_im(row, col).imag(); + this->in[row * ncols + col].re = original_im(row, col).real(); + this->in[row * ncols + col].im = original_im(row, col).imag(); } - this->p = fftw_plan_dft_2d(original_im.nrows(), original_im.ncols(), - this->in, this->out, - FFTW_FORWARD, FFTW_ESTIMATE); - this->p_inv = fftw_plan_dft_2d(original_im.nrows(), original_im.ncols(), - this->out, this->in, + this->p = fftw_plan_dft_2d(nrows, ncols, this->in, this->out, + FFTW_FORWARD, FFTW_ESTIMATE); + this->p_inv = fftw_plan_dft_2d(nrows, ncols, this->out, this->in, FFTW_BACKWARD, FFTW_ESTIMATE); - this->trans_im = image2d< std::complex<T> >(original_im.domain()); } - /*! - ** \brief Compute and return the transform. - */ + /// \brief Compute and return the transform (as a complex image). image2d< std::complex<T> > transform() { fftw_execute(this->p); - unsigned denom = this->trans_im.nrows() * this->trans_im.ncols(); + const unsigned nrows = this->trans_im.nrows(); + const unsigned ncols = this->trans_im.ncols(); + unsigned denom = nrows * ncols; int i = 0; - for (unsigned row = 0; row < this->trans_im.nrows(); ++row) - for (unsigned col = 0; col < this->trans_im.ncols(); ++col) + for (unsigned row = 0; row < nrows; ++row) + for (unsigned col = 0; col < ncols; ++col) { this->out[i].re /= denom; this->out[i].im /= denom; - this->trans_im(row, col) = std::complex<double>(this->out[i].re, - this->out[i].im); + this->trans_im(row, col) = + std::complex<double>(this->out[i].re, this->out[i].im); ++i; } return this->trans_im; } - /*! - ** \brief Compute and return the invert transform. - ** - ** \param R Data type of output image. - */ + /// \brief Compute and return the inverse transform (as a + /// complex image). + /// + /// \tparam R Component type of the value type of the output + /// image. template <class R> image2d< std::complex<R> > transform_inv() { - for (unsigned row = 0; row < this->trans_im.nrows(); ++row) - for (unsigned col = 0; col < this->trans_im.ncols(); ++col) + const unsigned nrows = this->trans_im.nrows(); + const unsigned ncols = this->trans_im.ncols(); + for (unsigned row = 0; row < nrows; ++row) + for (unsigned col = 0; col < ncols; ++col) { - this->out[row * this->trans_im.ncols() + col].re = this->trans_im(row, col).real(); - this->out[row * this->trans_im.ncols() + col].im = this->trans_im(row, col).imag(); + this->out[row * ncols + col].re = + this->trans_im(row, col).real(); + this->out[row * ncols + col].im = + this->trans_im(row, col).imag(); } fftw_execute(this->p_inv); image2d< std::complex<R> > new_im(this->trans_im.size()); int i = 0; - for (unsigned row = 0; row < this->trans_im.nrows(); ++row) - for (unsigned col = 0; col < this->trans_im.ncols(); ++col) + for (unsigned row = 0; row < nrows; ++row) + for (unsigned col = 0; col < ncols; ++col) { new_im(row, col) = this->in[i]; ++i; @@ -582,9 +552,7 @@ namespace mln { return new_im; } - /*! - ** \brief Compute and return the invert transform. - */ + /// \brief Compute and return the inverse transform of the FFT. image2d< std::complex<T> > transform_inv() { return transform_inv<T>(); @@ -592,8 +560,8 @@ namespace mln { }; - } // end of namespace transform + } // end of namespace mln::transform -} // end of namespace oln +} // end of namespace mln #endif // ! MLN_TRANSFORM_FFT_HH -- 1.7.10.4

11 years, 4 months

olena: olena-2.0-635-g089d5ff Remove dead code in the Fast Fourier Transform.

by Roland Levillain

* mln/transform/fft.hh: Here. --- milena/ChangeLog | 6 +++ milena/mln/transform/fft.hh | 93 ------------------------------------------- 2 files changed, 6 insertions(+), 93 deletions(-) diff --git a/milena/ChangeLog b/milena/ChangeLog index e163188..6a5f5d7 100644 --- a/milena/ChangeLog +++ b/milena/ChangeLog @@ -1,5 +1,11 @@ 2012-10-10 Roland Levillain <roland(a)lrde.epita.fr> + Remove dead code in the Fast Fourier Transform. + + * mln/transform/fft.hh: Here. + +2012-10-10 Roland Levillain <roland(a)lrde.epita.fr> + Minor corrections in the Fast Fourier Transform. * mln/transform/fft.hh diff --git a/milena/mln/transform/fft.hh b/milena/mln/transform/fft.hh index ed1e87f..dfbe1d7 100644 --- a/milena/mln/transform/fft.hh +++ b/milena/mln/transform/fft.hh @@ -479,88 +479,6 @@ namespace mln { } /*! - ** \brief Shift zero-frequency component of discrete Fourier transform - ** to center of spectrum. - ** - ** \param R The data type of the image returned. - ** - ** The zero-frequency component of discrete Fourier transform are moved - ** to center of the image : - ** - ** \htmlonly - ** <table> - ** <tr><td>1</td><td>2</td></tr> - ** <tr><td>3</td><td>4</td></tr> - ** </table> - ** becomes - ** <table> - ** <tr><td>4</td><td>3</td></tr> - ** <tr><td>2</td><td>1</td></tr> - ** </table> - ** \endhtmlonly - ** - */ - - /* - template <class R> - image2d<R> shift_transform_inv() - { - image2d<R> t = transform_inv<R>(); - image2d<R> st(t.size()); - - // We have to exchange t_1 with t_1_dest and not directly t_3 because - // they have not he same size. - typedef morpher::piece_morpher< image2d<R> > piece_t; - piece_t t_1(t, dpoint2d(0, 0), - image2d_size((t.size().nrows() - 1) / 2, - (t.size().ncols() - 1) / 2, - t.border())); - piece_t t_1_dest(st, dpoint2d(t.nrows() - t_1.nrows(), - t.ncols() - t_1.ncols()), - image2d_size(t_1.nrows(), t_1.ncols(), - t.border())); - piece_t t_2(t, dpoint2d(0, (t.size().ncols() - 1) / 2), - image2d_size((t.size().nrows() - 1) / 2, - t.size().ncols() - (t.size().ncols() - 1) / 2, - t.border())); - piece_t t_2_dest(st, dpoint2d(t.nrows() - t_2.nrows(), 0), - image2d_size(t_2.nrows(), t_2.ncols(), - t.border())); - piece_t t_3(t, dpoint2d((t.size().nrows() - 1) / 2, 0), - image2d_size(t.size().nrows() - (t.size().nrows() - 1) / 2, - (t.size().ncols() - 1) / 2, - t.border())); - piece_t t_3_dest(st, dpoint2d(0, t.ncols() - t_3.ncols()), - image2d_size(t_3.nrows(), t_3.ncols(), - t.border())); - piece_t t_4(t, dpoint2d((t.size().nrows() - 1) / 2, - (t.size().ncols() - 1) / 2), - image2d_size(t.size().nrows() - (t.size().nrows() - 1) / 2, - t.size().ncols() - (t.size().ncols() - 1) / 2, - t.border())); - piece_t t_4_dest(st, dpoint2d(0, 0), - image2d_size(t_4.nrows(), t_4.ncols(), - t.border())); - - oln_iter_type(piece_t) i1(t_1); - for_all(i1) - t_1_dest[i1] = t_1[i1]; - oln_iter_type(piece_t) i2(t_2); - for_all(i2) - t_2_dest[i2] = t_2[i2]; - oln_iter_type(piece_t) i3(t_3); - for_all(i3) - t_3_dest[i3] = t_3[i3]; - oln_iter_type(piece_t) i4(t_4); - for_all(i4) - t_4_dest[i4] = t_4[i4]; - - return st; - } - */ - - - /*! ** \brief Compute and return the invert transform. */ image2d<T> transform_inv() @@ -568,17 +486,6 @@ namespace mln { return transform_inv<T>(); } - /*! - ** \brief Shift zero-frequency component of discrete Fourier transform - ** to center of spectrum. - */ - /* - image2d<T> shift_transform_inv() - { - return shift_transform_inv<T>(); - } - */ - }; /*! -- 1.7.10.4

11 years, 4 months

olena: olena-2.0-634-g6919bfe Minor corrections in the Fast Fourier Transform.

by Roland Levillain

* mln/transform/fft.hh (internal::_fft<T>::transformed_image() const): Return a const reference instead of a copy. (internal::_fft<T>::transformed_image_clipped_magn(double, bool)): Drop the const before the first argument. Properly initialize `max'. (internal::_fft<T>::transformed_image_log_magn(double, double, bool)): Properly initialize `max'. --- milena/ChangeLog | 13 +++++++++++++ milena/mln/transform/fft.hh | 10 +++++----- 2 files changed, 18 insertions(+), 5 deletions(-) diff --git a/milena/ChangeLog b/milena/ChangeLog index 9cfd3a6..e163188 100644 --- a/milena/ChangeLog +++ b/milena/ChangeLog @@ -1,3 +1,16 @@ +2012-10-10 Roland Levillain <roland(a)lrde.epita.fr> + + Minor corrections in the Fast Fourier Transform. + + * mln/transform/fft.hh + (internal::_fft<T>::transformed_image() const): Return a const + reference instead of a copy. + (internal::_fft<T>::transformed_image_clipped_magn(double, bool)): + Drop the const before the first argument. + Properly initialize `max'. + (internal::_fft<T>::transformed_image_log_magn(double, double, bool)): + Properly initialize `max'. + 2012-10-09 Roland Levillain <roland(a)lrde.epita.fr> Improve the test of the Fast Fourier Transform. diff --git a/milena/mln/transform/fft.hh b/milena/mln/transform/fft.hh index a010e47..ed1e87f 100644 --- a/milena/mln/transform/fft.hh +++ b/milena/mln/transform/fft.hh @@ -81,7 +81,7 @@ namespace mln { ** ** Const version. */ - const image2d< std::complex<T> > transformed_image() const + const image2d< std::complex<T> >& transformed_image() const { return trans_im; } @@ -156,14 +156,14 @@ namespace mln { */ template <class R> - image2d<R> transformed_image_clipped_magn(const double clip, + image2d<R> transformed_image_clipped_magn(double clip, bool ordered = true) const { // Check that R is real image2d<R> new_im(trans_im.domain()); // check that clip is >=0 and <=1 ? - double max; + double max = mln_min(double); mln_piter(image2d<T>) it(trans_im.domain()); for_all(it) @@ -209,7 +209,7 @@ namespace mln { ** \param clip Value used for clipping. ** \param ordered Kind of traversal. */ - image2d<T> transformed_image_clipped_magn(const double clip, + image2d<T> transformed_image_clipped_magn(double clip, bool ordered = true) const { return transformed_image_clipped_magn<T>(clip, ordered); @@ -271,7 +271,7 @@ namespace mln { image2d<R> new_im(trans_im.domain()); - double max = 0; + double max = mln_min(double); mln_piter(image2d<R>) it(trans_im.domain()); for_all(it) -- 1.7.10.4

11 years, 4 months

← Newer
1
...
7
8
9
10
11
12
13
...
18
Older →

Jump to page:

2024

2023

2022

2021

2020

2019

2018

2017

2016

2015

2014

2013

2012

2011

2010

2009

2008

2007

2006

2005

2004

Olena-patches August 2013