From f9e8952216d066d50eea1f2b21f25619375bb5d5 Mon Sep 17 00:00:00 2001 From: Samuel Farrens Date: Mon, 14 Sep 2020 10:48:42 +0200 Subject: [PATCH] BigMac support for AppleClang (#33) * set macOS FFTW flags * fixed gitignore bug * added NFFT fix * updated fftw and nfft builds * updated travis * added Dockerfile and updated readme * added build for mr2d1d_stat * cleaned up CMakeLists * changed lib order * updated .gitignore * added nfft fix for linux * set macOS FFTW flags * added build for mr2d1d_stat * added new libraries * updated fftw build instructions * resolved issue with FFTW headers * fixed library def --- .gitignore | 37 +- .travis.yml | 75 +-- CMakeLists.txt | 332 +++++------- Dockerfile | 22 + README.md | 106 ++-- modules/BuildFFTW.cmake | 63 +++ modules/BuildFFTW3.cmake | 29 -- modules/BuildNFFT.cmake | 51 +- modules/fftw.cmd | 15 - modules/functions.cmake | 43 ++ modules/nfft.cmd | 14 - modules/nfft_fftw.cmd | 15 - src/libsparse1d/SB_Filter1D.h | 104 ++-- src/libsparse2d/MW_Deconv.cc | 687 +++++++++++++++++++++++++ src/libsparse2d/MW_Deconv.h | 122 +++++ src/libsparse2d/MW_Filter.cc | 518 +++++++++++++++++++ src/libsparse2d/MW_Filter.h | 153 ++++++ src/libtools/CMem.cc | 231 +++++++++ src/libtools/CMem.h | 87 ++++ src/libtools/IM_Block2D.cc | 187 +++++++ src/libtools/IM_Block2D.h | 69 +++ src/libtools/IM_DCT.cc | 592 +++++++++++++++++++++ src/libtools/IM_DCT.h | 76 +++ src/{ => msvst}/msvst_2d1d.cc | 0 src/{ => msvst}/msvst_iwt2d.cc | 0 src/{ => msvst}/msvst_iwt2d_coupled.cc | 0 src/{ => msvst}/msvst_uwt2d.cc | 0 src/{ => unused}/cb_filter.cc | 0 src/{ => unused}/cb_mca.cc | 0 src/{ => unused}/mr1d_deconv.cc | 0 30 files changed, 3210 insertions(+), 418 deletions(-) create mode 100644 Dockerfile create mode 100644 modules/BuildFFTW.cmake delete mode 100644 modules/BuildFFTW3.cmake delete mode 100755 modules/fftw.cmd create mode 100644 modules/functions.cmake delete mode 100755 modules/nfft.cmd delete mode 100755 modules/nfft_fftw.cmd create mode 100644 src/libsparse2d/MW_Deconv.cc create mode 100644 src/libsparse2d/MW_Deconv.h create mode 100644 src/libsparse2d/MW_Filter.cc create mode 100644 src/libsparse2d/MW_Filter.h create mode 100644 src/libtools/CMem.cc create mode 100644 src/libtools/CMem.h create mode 100644 src/libtools/IM_Block2D.cc create mode 100644 src/libtools/IM_Block2D.h create mode 100644 src/libtools/IM_DCT.cc create mode 100644 src/libtools/IM_DCT.h rename src/{ => msvst}/msvst_2d1d.cc (100%) rename src/{ => msvst}/msvst_iwt2d.cc (100%) rename src/{ => msvst}/msvst_iwt2d_coupled.cc (100%) rename src/{ => msvst}/msvst_uwt2d.cc (100%) rename src/{ => unused}/cb_filter.cc (100%) rename src/{ => unused}/cb_mca.cc (100%) rename src/{ => unused}/mr1d_deconv.cc (100%) diff --git a/.gitignore b/.gitignore index 4ab7804..dce5ed3 100644 --- a/.gitignore +++ b/.gitignore @@ -1,2 +1,35 @@ -# Ignore build files -*build* +# Prerequisites +*.d + +# Compiled Object files +*.slo +*.lo +*.o +*.obj + +# Precompiled Headers +*.gch +*.pch + +# Compiled Dynamic libraries +*.so +*.dylib +*.dll + +# Fortran module files +*.mod +*.smod + +# Compiled Static libraries +*.lai +*.la +*.a +*.lib + +# Executables +*.exe +*.out +*.app + +# Ignore build directory +build diff --git a/.travis.yml b/.travis.yml index d0c238f..111bd19 100644 --- a/.travis.yml +++ b/.travis.yml @@ -1,55 +1,62 @@ # System Set-Up language: cpp -# install dependencies -addons: - apt: - packages: - - cmake - - libcfitsio3-dev - homebrew: - packages: &macos_packages - - cmake - - pkgconfig - - cfitsio - update: true +branches: + only: + - master -matrix: +jobs: include: - os: linux + dist: focal + compiler: gcc + name: "Ubuntu Focal" + env: + - CC=gcc-9 + - CXX=g++-9 + addons: + apt: + sources: + - sourceline: 'ppa:ubuntu-toolchain-r/test' + packages: + - gcc-9 + - g++-9 + - cmake + - libcfitsio-dev - os: osx - name: "macOS with GCC" - osx_image: xcode11.3 + osx_image: xcode12 + language: shell + name: "macOS 10.15 GCC" env: - - MATRIX_EVAL="CC=gcc-9 && CXX=g++-9" + - CC=gcc-9 + - CXX=g++-9 addons: homebrew: packages: - - gcc - - *macos_packages + - cmake + - cfitsio + - gcc@9 + update: true - os: osx - name: "macOS with AppleClang" - osx_image: xcode11.3 + osx_image: xcode12 + language: shell + name: "macOS 10.15 AppleClang" addons: homebrew: + taps: sfarrens/sf packages: - - libomp - - *macos_packages + - cmake + - cfitsio + - bigmac + update: true before_install: - - eval "${MATRIX_EVAL}" - -# set up installation -install: - mkdir build - cd build - - cmake .. -DRUN_TESTS=OFF -# build package and install -script: - - make install - # - make test +install: + - cmake .. -DBUILD_MSVST=ON + - make -# notification settings -notification: - email: false +script: + - make test diff --git a/CMakeLists.txt b/CMakeLists.txt index c67df8a..cbfc1d7 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -3,212 +3,140 @@ ################# cmake_minimum_required(VERSION 3.0.0) +project(sparse2d) + +# Include modules in CMake module path +list(APPEND CMAKE_MODULE_PATH "${CMAKE_SOURCE_DIR}/modules/") + +# Load CMake tools include(ExternalProject) include(FindPkgConfig) -list(APPEND CMAKE_MODULE_PATH "${CMAKE_SOURCE_DIR}/modules/") +# Load custom CMake functions +include(functions) + +# Set build type +if(CMAKE_BUILD_TYPE STREQUAL "") + set(CMAKE_BUILD_TYPE "Release") +endif() +message(STATUS "CMake Build Type: ${CMAKE_BUILD_TYPE}") + +# ----------------- # +# Find Dependencies # +# ----------------- # + +# Locate OpenMP +if("${CMAKE_CXX_COMPILER_ID}" STREQUAL "AppleClang") + find_package(BigMac 0.0.1 REQUIRED) +else() + find_package(OpenMP REQUIRED) +endif() + +# Locate CFITSIO +find_cfitsio() + +# Set build location for external modules +set(MODULE_BUILD_DIR "${CMAKE_CURRENT_BINARY_DIR}/module_build/") + +# Set compilation flags +set(CMAKE_CXX_FLAGS + "${CMAKE_CXX_FLAGS} -DNO_DISP_IO -fPIC -Wno-write-strings\ + ${OpenMP_CXX_FLAGS} ${BigMac_NOWARN}" +) +set(CMAKE_CXX_FLAGS_RELEASE + "${CMAKE_CXX_FLAGS_RELEASE} -g0 -fomit-frame-pointer" +) +set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -O0") + +# ---- # +# FFTW # +# -----# + +# Optional use of FFTW +option(USE_FFTW "Use FFTW for Fourier transformation" OFF) + +# Optional build of FFTW +option(BUILD_FFTW "Build FFTW libraries" OFF) +include(BuildFFTW) + +# ---- # +# NFFT # +# -----# + +# Optional build of NFFT +option(BUILD_NFFT "Download and build NFFT" OFF) +include(BuildNFFT) + +# ------------ # +# Sparse2D STD # +# -------------# + +# Build standard Sparse2D libs +set(sparse2d_libs mga2d sparse3d sparse2d sparse1d tools) +foreach(library ${sparse2d_libs}) + build_lib(${library}) +endforeach() + +# Find all Sparse2D targets +find_targets(sparse2d_targets src cc) + +# Build binaries +foreach(program ${sparse2d_targets}) + build_bin(${program} "${sparse2d_libs}" src cc) +endforeach(program) + +# Install libraries +install(TARGETS ${sparse2d_libs} DESTINATION lib EXPORT sparse2d_libs) +install(EXPORT sparse2d_libs DESTINATION lib/cmake) + +# Install binaraies +install(TARGETS ${sparse2d_targets} DESTINATION bin) + +# ----- # +# MSVST # +# ------# + +option(BUILD_MSVST "Build MSVST library" OFF) +message(STATUS "MSVST Build: ${BUILD_MSVST}") + +if(BUILD_MSVST) + + # Build MSVST lib + build_lib(msvst) + set(msvst_libs msvst tools) + + # Find all MSVST targets + find_targets(msvst_targets src/msvst cc) + + # Build binaries + foreach(program ${msvst_targets}) + build_bin(${program} "${msvst_libs}" src/msvst cc) + endforeach(program) -project(sparse2d) + # Install binaraies + install(TARGETS ${msvst_targets} DESTINATION bin) + +endif() + +# ---------- # +# UNIT TESTS # +# -----------# + +option(RUN_TESTS "Build and run unit tests" ON) +message(STATUS "Run Tests: ${RUN_TESTS}") + +if(RUN_TESTS) + + # Add unit tests + enable_testing() + + # Find all unit test targets + find_targets(unit_tests tests cpp) - # Set build type - if(CMAKE_BUILD_TYPE STREQUAL "") - set(CMAKE_BUILD_TYPE "Release") - endif() - message(STATUS "CMake Build Type: ${CMAKE_BUILD_TYPE}") - - # Set the default installation path - if(CMAKE_INSTALL_PREFIX_INITIALIZED_TO_DEFAULT) - set(CMAKE_INSTALL_PREFIX "${CMAKE_CURRENT_BINARY_DIR}/sparse2d" - CACHE PATH "Setting new default path" FORCE) - endif() - message(STATUS "Installing to: ${CMAKE_INSTALL_PREFIX}") - - if("${CMAKE_CXX_COMPILER_ID}" STREQUAL "AppleClang") - set(OpenMP_INCLUDE_PATH "/usr/local/include") - set(OpenMP_LIB_PATH "/usr/local/lib") - set(OpenMP_CXX_FLAGS "-Xpreprocessor -fopenmp -lomp") - set(OpenMP_CXX_LIB_NAMES "omp") - set(OpenMP_CXX_LIBRARIES "${OpenMP_LIB_PATH}/libomp.dylib") - include_directories(${OpenMP_INCLUDE_PATH}) - link_directories(${OpenMP_LIB_PATH}) - else() - find_package(OpenMP REQUIRED) - endif() - - # Locate CFITSIO using pkg-config or use command line arguments to configure - # CFITSIO - if(CFITSIO_LIBRARIES STREQUAL "" OR NOT DEFINED CFITSIO_LIBRARIES OR - CFITSIO_LIBRARY_DIRS STREQUAL "" OR NOT DEFINED CFITSIO_LIBRARY_DIRS OR - CFITSIO_INCLUDE_DIRS STREQUAL "" OR NOT DEFINED CFITSIO_INCLUDE_DIRS) - pkg_check_modules(CFITSIO REQUIRED cfitsio) - else() - message(STATUS "Use manually configured cfitsio") - message(STATUS " includes: ${CFITSIO_INCLUDE_DIRS}") - message(STATUS " libs: ${CFITSIO_LIBRARY_DIRS}") - message(STATUS " flags: ${CFITSIO_LIBRARIES}") - endif() - include_directories(${CFITSIO_INCLUDE_DIRS}) - link_directories(${CFITSIO_LIBRARY_DIRS}) - - # Build location for external modules - set(MODULE_BUILD_DIR "${CMAKE_CURRENT_BINARY_DIR}/module_build/") - - # Optional build of FFTW3 - option(USE_FFTW "Use FFTW for Fourier transformation" OFF) - include(BuildFFTW3) - - # Optional build of NFFT - option(COMPILE_NFFT "Get and compile NFFT" OFF) - include(BuildNFFT) - - # Set compilation flags Mac OSX or any other system - if("${CMAKE_CXX_COMPILER_ID}" STREQUAL "AppleClang") - - set(CMAKE_CXX_FLAGS_RELEASE "-DNO_DISP_IO -g0 -O2 -fPIC \ --fomit-frame-pointer ${OpenMP_CXX_FLAGS} -Wno-write-strings -DNDEBUG \ -${FFTW_FLAGS}") - set(CMAKE_CXX_FLAGS_DEBUG "-DNO_DISP_IO -g -fPIC ${OpenMP_CXX_FLAGS} \ --Wno-write-strings ${FFTW_CXX_FLAGS}") - - else() - - set(CMAKE_CXX_FLAGS_RELEASE "-DNO_DISP_IO -g0 -O2 -fPIC \ --fomit-frame-pointer ${OpenMP_CXX_FLAGS} -Wno-write-strings -DNDEBUG \ -${FFTW_FLAGS}" - ) - set(CMAKE_CXX_FLAGS_DEBUG "-DNO_DISP_IO -g -fPIC ${OpenMP_CXX_FLAGS} \ --Wno-write-strings ${FFTW_CXX_FLAGS}" - ) - - message(STATUS "Using CXX Flags for ${CMAKE_SYSTEM}") - - endif() - - # Build tools library - FILE(GLOB src_lib1 "${PROJECT_SOURCE_DIR}/src/libtools/*.cc") - include_directories("${PROJECT_SOURCE_DIR}/src/libtools") - add_library(tools STATIC ${src_lib1}) - target_link_libraries(tools ${CFITSIO_LIBRARIES} ${FFTW_LD_FLAGS}) - if(USE_FFTW) - add_dependencies(tools fftw3) - endif(USE_FFTW) - - # Build sparse1d library - FILE(GLOB src_lib2 "${PROJECT_SOURCE_DIR}/src/libsparse1d/*.cc") - include_directories("${PROJECT_SOURCE_DIR}/src/libsparse1d") - add_library(sparse1d STATIC ${src_lib2}) - target_link_libraries(sparse1d ${CFITSIO_LIBRARIES} ${FFTW_LD_FLAGS}) - if(USE_FFTW) - add_dependencies(sparse1d fftw3) - endif(USE_FFTW) - - # Build sparse2d library - FILE(GLOB src_lib2 "${PROJECT_SOURCE_DIR}/src/libsparse2d/*.cc") - include_directories("${PROJECT_SOURCE_DIR}/src/libsparse2d") - add_library(sparse2d STATIC ${src_lib2}) - target_link_libraries(sparse2d ${CFITSIO_LIBRARIES} ${FFTW_LD_FLAGS}) - if(USE_FFTW) - add_dependencies(sparse2d fftw3) - endif(USE_FFTW) - - option(SPARSE3D "Build Sparse3D library" ON) - if(SPARSE3D) - # Build sparse3d library - FILE(GLOB src_lib2 "${PROJECT_SOURCE_DIR}/src/libsparse3d/*.cc") - include_directories("${PROJECT_SOURCE_DIR}/src/libsparse3d") - add_library(sparse3d STATIC ${src_lib2}) - target_link_libraries(sparse3d ${CFITSIO_LIBRARIES} ${FFTW_LD_FLAGS}) - if(USE_FFTW) - add_dependencies(sparse3d fftw3) - endif(USE_FFTW) - endif(SPARSE3D) - message(STATUS "Sparse3D Build: ${SPARSE3D}") - - # Build mga2d library - FILE(GLOB src_mgalib2 "${PROJECT_SOURCE_DIR}/src/libmga2d/*.cc") - include_directories("${PROJECT_SOURCE_DIR}/src/libmga2d") - add_library(mga2d STATIC ${src_mgalib2}) - target_link_libraries(mga2d ${CFITSIO_LIBRARIES} ${FFTW_LD_FLAGS}) - if(USE_FFTW) - add_dependencies(mga2d fftw3) - endif(USE_FFTW) - - option(BUILD_MSVST "Build MSVST library" OFF) - if(BUILD_MSVST) - FILE(GLOB msvst_src_lib "${PROJECT_SOURCE_DIR}/src/libmsvst/*.cc") - include_directories("${PROJECT_SOURCE_DIR}/src/libmsvst") - add_library(msvst STATIC ${msvst_src_lib} ) - endif(BUILD_MSVST) - message(STATUS "MSVST Build: ${BUILD_MSVST}") - - # Compile and link executables - set(BINMR2D mr_transform mr_recons mr_filter mr_deconv mr_info cur_contrast - cur_deconv cur_filter cur_stat cur_trans im1d_convert im1d_info im1d_stf - im1d_tfreq mr1d_detect mr1d_trans im1d_der_loggamma im1d_period im1d_tend - mr1d_filter im_isospec) - if(SPARSE3D) - set(BINMR2D ${BINMR2D} mr3d_trans mr3d_filter mr3d_recons mr3d_stat - mr2d1d_trans) - endif(SPARSE3D) - - foreach(program ${BINMR2D}) - add_executable(${program} ${PROJECT_SOURCE_DIR}/src/${program}.cc) - target_link_libraries(${program} mga2d sparse2d sparse1d tools) - if(SPARSE3D) - target_link_libraries(${program} mga2d sparse3d sparse2d sparse1d tools) - endif(SPARSE3D) + # Build binaries + foreach(program ${unit_tests}) + build_bin(${program} "${sparse2d_libs}" tests cpp) + add_test(${program} ${program}) endforeach(program) - if(BUILD_MSVST) - set(MSVST msvst_uwt2d msvst_iwt2d msvst_2d1d msvst_iwt2d_coupled) - foreach(program ${MSVST}) - add_executable(${program} ${PROJECT_SOURCE_DIR}/src/${program}.cc) - target_link_libraries(${program} msvst tools) - endforeach(program) - FILE(GLOB inc_lib "${PROJECT_SOURCE_DIR}/src/libmsvst/*.h") - INSTALL(FILES ${inc_lib} DESTINATION include/sparse2d) - INSTALL(TARGETS ${MSVST} DESTINATION bin) - endif(BUILD_MSVST) - - # Install headers - FILE(GLOB inc_lib "${PROJECT_SOURCE_DIR}/src/libsparse1d/*.h") - INSTALL(FILES ${inc_lib} DESTINATION include/sparse2d) - FILE(GLOB inc_lib "${PROJECT_SOURCE_DIR}/src/libsparse2d/*.h") - INSTALL(FILES ${inc_lib} DESTINATION include/sparse2d) - FILE(GLOB inc_lib "${PROJECT_SOURCE_DIR}/src/libmga2d/*.h") - INSTALL(FILES ${inc_lib} DESTINATION include/sparse2d) - FILE(GLOB inc_lib "${PROJECT_SOURCE_DIR}/src/libtools/*.h") - INSTALL(FILES ${inc_lib} DESTINATION include/sparse2d) - if(SPARSE3D) - FILE(GLOB inc_lib "${PROJECT_SOURCE_DIR}/src/libsparse3d/*.h") - INSTALL(FILES ${inc_lib} DESTINATION include/sparse2d) - endif(SPARSE3D) - - # Install library - INSTALL(TARGETS sparse1d sparse2d mga2d tools DESTINATION lib) - if(SPARSE3D) - INSTALL(TARGETS sparse3d DESTINATION lib) - endif(SPARSE3D) - - # install executables - INSTALL(TARGETS ${BINMR2D} DESTINATION bin) - - # ---------- # - # UNIT TESTS # - # -----------# - - option(RUN_TESTS "Build and run unit tests" ON) - if(RUN_TESTS) - # Add unit tests - enable_testing() - set(UNIT_TESTS mr_transform_test) - - foreach(program ${UNIT_TESTS}) - add_executable(${program} ${PROJECT_SOURCE_DIR}/tests/${program}.cpp) - target_link_libraries(${program} mga2d sparse2d sparse1d tools) - if(SPARSE3D) - target_link_libraries(${program} mga2d sparse3d sparse2d sparse1d tools) - endif(SPARSE3D) - add_test(${program} ${program}) - endforeach(program) - endif(RUN_TESTS) +endif() diff --git a/Dockerfile b/Dockerfile new file mode 100644 index 0000000..c0c02f2 --- /dev/null +++ b/Dockerfile @@ -0,0 +1,22 @@ +FROM ubuntu + +LABEL Description="Sparse2D" + +ARG DEBIAN_FRONTEND=noninteractive +ARG CC=gcc-9 +ARG CXX=g++-9 + +RUN apt-get update && \ + apt-get install -y autoconf automake libtool pkg-config && \ + apt-get install -y gcc-9 g++-9 && \ + apt-get install -y git cmake libcfitsio-dev && \ + apt-get clean + +RUN cd home && \ + git clone https://github.com/CosmoStat/Sparse2D && \ + cd Sparse2D && \ + mkdir build && \ + cd build && \ + cmake .. -DBUILD-MSVST=ON -DUSE_FFTW=ON -DBUILD-NFFT=ON && \ + make && \ + make install diff --git a/README.md b/README.md index d585971..01c392b 100644 --- a/README.md +++ b/README.md @@ -1,85 +1,119 @@ # Sparse2D [![Build Status](https://travis-ci.org/CosmoStat/Sparse2D.svg?branch=master)](https://travis-ci.org/CosmoStat/Sparse2D) +[![cpp](https://img.shields.io/badge/language-C%2B%2B-red)](https://isocpp.org/std/the-standard) Sparse2D provides an array of sparsity based tools and a convenient C++ library for performing various wavelet tranforms. -This package is part of the iSAP suite, available at the [CosmoStat webpage](http://www.cosmostat.org/software/isap) +This package is part of the iSAP suite, available on the [CosmoStat website](http://www.cosmostat.org/software/isap). + +Python bindings to Sparse2D are provided in [PySAP](https://github.com/CEA-COSMIC/pysap). PySAP handles Sparse2D installation internally. ## Installation instructions ### Prerequisites -The following softwares are required: - - - C++ compiler (gcc strongly recommended) - - Recent version of CFITSIO ( >V3.31) needs to be first installed. - - CMake (http://www.cmake.org) and pkg-config + - C/C++ compiler + - [CMake](http://www.cmake.org) + - [CFITSIO](https://heasarc.gsfc.nasa.gov/fitsio/) ( >V3.31) + - [pkg-config](https://www.freedesktop.org/wiki/Software/pkg-config/) + - [BigMac](https://github.com/sfarrens/bigmac) (For macOS AppleClang compiler) -Please use a package management tool to properly install cfistio and the other dependencies. On linux you can use apt (Ubuntu). -On Mac OSX we recommend homebrew, the command "brew install cfitsio" will install the cfitsio package. +Please use a package management tool to properly install cfistio and the other dependencies on linux (*e.g.* `apt` on Ubuntu or `brew` on macOS). -### Building from Source +### Docker -To configure the sparse2d toolbox for your machine type: +If you have [Docker](https://www.docker.com/) installed, you can pull the latest build of the Sparse2D image from [Docker Hub](TBD). ```bash - $ mkdir build - $ cd build - $ cmake .. +$ docker pull cosmostat/sparse2d ``` -Once CMake is done, to create the toolbox and compile all the C++ files, type: +No further installation is required. + +To run this image on data in your current working directory, simply run: ```bash - $ make - $ make install +$ docker run -v ${PWD}:/workdir -it cosmostat/sparse2d /bin/bash -c "cd workdir && " ``` -The library, headers and executables will by default be installed in a directory -called `sparse2d` inside the current build directory. The default install directory -can be changed by running: +where `` is one of the Sparse2D binaries. The reference to `${PWD}` can be replaced by the path to any directory on your system and options can be passed to `` inside the double quotes. + +### Homebrew + +Sparse2D can be built on macOS using [Homebrew](https://brew.sh/). ```bash - $ cmake .. -DCMAKE_INSTALL_PREFIX:PATH= +$ brew tap cosmostat/science +$ brew install sparse2d +``` + +The Homebrew formula handles all of the required dependencies. + +### Building from Source + +#### Basic Sparse2D + +Clone the repository: + ``` +$ git clone https://github.com/CosmoStat/Sparse2D.git +$ cd Sparse2D +``` + +Create a build directory for Sparse2D: -You may want to add the bin directory to your path for easy access to the sparse2d executables. +```bash +$ mkdir build +$ cd build +``` -Sparse2D can also be built using [FFTW3](http://www.fftw.org/) and/or -[nFFT](https://www-user.tu-chemnitz.de/~potts/nfft/) with the following options: +Build Sparse2D: ```bash - $ cmake .. -DUSE_FFTW=ON -DCOMPILE_NFFT=ON +$ cmake .. +$ make +$ make install ``` -Finally, if you wish to build using a compiler other than the default on your -system (*e.g.* gcc on macOS) you can do so as follows: +#### MSVST + +Additional MSVST binaries can be build with the following option: ```bash - $ CC=gcc CXX=g++ cmake .. +$ cmake .. -DBUILD_MSVST=ON ``` -### macOS +#### FFTW -Sparse2D can be built on macOS using [Homebrew](https://brew.sh/). Before doing so, you should make sure you have a C compiler that supports OpenMP (note that the native clang does not) and cmake. These requirements can also be installed with Homebrew as follows +Sparse2D can also be built using pre-installed [FFTW](http://www.fftw.org/) libraries with the following option: ```bash - $ brew install gcc cmake +$ cmake .. -DUSE_FFTW=ON ``` -Once these packages have been installed, Sparse2D can be installed with +Alternatively, Sparse2D can build FFTW libraries from source as follows: ```bash - $ brew tap cosmostat/science - $ brew install sparse2d +$ cmake .. -DUSE_FFTW=ON -DBUILD_FFTW=ON ``` -By default this package builds with Homebrew gcc, however you can also build with an alternative compiler as follows +#### nFFT + +Sparse2D can also build [nFFT](https://github.com/NFFT/nfft) libraries with the following options: + +```bash +$ cmake .. -DBUILD_FFTW=ON -DBUILD_NFFT=ON +``` + +#### Non-default Compiler + +Finally, if you wish to build using a compiler other than the default on your +system (*e.g.* `gcc` on macOS) you can do so as follows: ```bash - $ brew install sparse2d --cc=COMPILER +$ CC=gcc CXX=g++ cmake .. ``` ## Usage -The two main executables of the package are **mr_transform** and **mr_filter**, see the instructions in the README/ folder for usage examples +The two main executables of the package are `mr_transform` and `mr_filter`, see the instructions in the [README](./README) folder for usage examples diff --git a/modules/BuildFFTW.cmake b/modules/BuildFFTW.cmake new file mode 100644 index 0000000..f8ccd1f --- /dev/null +++ b/modules/BuildFFTW.cmake @@ -0,0 +1,63 @@ +#=============# +# Build FFTW # +#=============# + +# Set FFTW Version +set(FFTWVersion 3.3.8) +set(FFTWHash 8aac833c943d8e90d51b697b27d4384d) + +# Find FFTW +if(USE_FFTW AND NOT BUILD_FFTW) + find_package(FFTW3 REQUIRED) + message(STATUS "Checking for module 'fftw3'") + message(STATUS " Found fftw3, version ${FFTW3_VERSION}") +elseif(USE_FFTW AND BUILD_FFTW) + set(FFTW3_INCLUDE_DIRS "${MODULE_BUILD_DIR}include") + set(FFTW3_LIBRARY_DIRS "${MODULE_BUILD_DIR}lib") +endif() + +# Set FFTW flags +if(USE_FFTW) + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DUSE_FFTW") + set(FFTW_LD_FLAGS "-lfftw3f_omp -lfftw3_omp -lfftw3f -lfftw3 -lm") +endif() + +# Build FFTW +if(BUILD_FFTW) + + # Set C/C++ compiler and flags + set(FFTW_COMPILE + CC=${CMAKE_C_COMPILER} + CXX=${CMAKE_CXX_COMPILER} + ${BigMac_FFTW}) + + # Set FFTW configuration flags + set(FFTW_CONFIG_FLAGS + --prefix=${MODULE_BUILD_DIR} + --enable-maintainer-mode + --enable-shared + --enable-threads + --enable-sse2 + --enable-openmp) + + # Download and build FFTW + ExternalProject_Add(fftw + URL http://www.fftw.org/fftw-${FFTWVersion}.tar.gz + URL_HASH MD5=${FFTWHash} + SOURCE_DIR ${MODULE_BUILD_DIR}fftw + BUILD_IN_SOURCE 1 + CONFIGURE_COMMAND ./configure ${FFTW_COMPILE} ${FFTW_CONFIG_FLAGS} + COMMAND make && make install + COMMAND ./configure ${FFTW_COMPILE} ${FFTW_CONFIG_FLAGS} --enable-single + COMMAND make && make install + COMMAND make clean + ) + +endif() + +# Source FFTW headers and libraries +include_directories(${FFTW3_INCLUDE_DIRS}) +link_directories(${FFTW3_LIBRARY_DIRS}) + +message(STATUS "Use FFTW: ${USE_FFTW}") +message(STATUS "FFTW Build: ${BUILD_FFTW}") diff --git a/modules/BuildFFTW3.cmake b/modules/BuildFFTW3.cmake deleted file mode 100644 index 5e9bbcf..0000000 --- a/modules/BuildFFTW3.cmake +++ /dev/null @@ -1,29 +0,0 @@ -#=============# -# Build FFTW3 # -#=============# - -set(FFTWVersion 3.3.8) - -if(USE_FFTW) - - ExternalProject_Add(fftw3 - URL http://www.fftw.org/fftw-${FFTWVersion}.tar.gz - CONFIGURE_COMMAND "" - BUILD_COMMAND "${CMAKE_SOURCE_DIR}/modules/fftw.cmd" - SOURCE_DIR "${MODULE_BUILD_DIR}/fftw" - INSTALL_COMMAND "" - BUILD_IN_SOURCE 1 - ) - - set(FFTW_CXX_FLAGS "-DUSE_FFTW") - set(FFTW_LD_FLAGS "-L ${MODULE_BUILD_DIR}/fftw/lib -lfftw3f_omp -lfftw3_omp \ --lfftw3f -lfftw3 -lm") - -else(USE_FFTW) - - set(FFTW_CXX_FLAGS "") - set(FFTW_LD_FLAGS "") - -endif(USE_FFTW) - -message(STATUS "FFTW3 Build: ${USE_FFTW}") diff --git a/modules/BuildNFFT.cmake b/modules/BuildNFFT.cmake index 4279f1a..1aed6d9 100644 --- a/modules/BuildNFFT.cmake +++ b/modules/BuildNFFT.cmake @@ -1,26 +1,37 @@ -#==============# -# Compile nFFT # -#==============# +#============# +# Build nFFT # +#============# -set(NFFTVersion 3.4.0) +# Set NFFT Version +set(NFFTVersion 3.5.1) +set(NFFTHash 8d53164d7cd85ad77e1bd03e36c4a99ef73c77f640e527db816cdc3fcb43d6aa) -if(COMPILE_NFFT) +if(BUILD_NFFT) - if(USE_FFTW) - set(BUILD_CMD nfft_fftw.cmd) - else(USE_FFTW) - set(BUILD_CMD nfft.cmd) - endif(USE_FFTW) + # Set C/C++ compiler and flags + set(NFFT_COMPILE + CC=${CMAKE_C_COMPILER} + CXX=${CMAKE_CXX_COMPILER} + ${BigMac_NFFT}) - ExternalProject_Add(nfft - URL https://github.com/NFFT/nfft/archive/${NFFTVersion}.tar.gz - CONFIGURE_COMMAND "" - BUILD_COMMAND "${CMAKE_SOURCE_DIR}/modules/${BUILD_CMD}" - SOURCE_DIR "${MODULE_BUILD_DIR}/nfft" - INSTALL_COMMAND "" - BUILD_IN_SOURCE 1 - ) + # Set NFFT configuration flags + set(NFFT_CONFIG_FLAGS + --prefix=${MODULE_BUILD_DIR} + --with-fftw3=${MODULE_BUILD_DIR} + --enable-openmp) -endif(COMPILE_NFFT) + # Download and build NFFT + ExternalProject_Add(nfft + URL https://github.com/NFFT/nfft/archive/${NFFTVersion}.tar.gz + URL_HASH SHA256=${NFFTHash} + SOURCE_DIR ${MODULE_BUILD_DIR}nfft + BUILD_IN_SOURCE 1 + CONFIGURE_COMMAND ./bootstrap.sh && + ./configure ${NFFT_COMPILE} ${NFFT_CONFIG_FLAGS} + BUILD_COMMAND export LD_LIBRARY_PATH=${MODULE_BUILD_DIR}lib && make -j8 + INSTALL_COMMAND make install + ) -message(STATUS "nFFT Compilation: ${COMPILE_NFFT}") +endif(BUILD_NFFT) + +message(STATUS "nFFT Build: ${BUILD_NFFT}") diff --git a/modules/fftw.cmd b/modules/fftw.cmd deleted file mode 100755 index aa632d3..0000000 --- a/modules/fftw.cmd +++ /dev/null @@ -1,15 +0,0 @@ -#! /bin/bash -dest_dir=`pwd` -touch ChangeLog -rm -rf autom4te.cache -autoreconf --verbose --install --symlink --force -autoreconf --verbose --install --symlink --force -autoreconf --verbose --install --symlink --force -rm -f config.cache -./configure --prefix=$dest_dir --enable-maintainer-mode --enable-shared --enable-sse2 --enable-openmp -make -make install -./configure --prefix=$dest_dir --enable-maintainer-mode --enable-shared --enable-single --enable-sse2 --enable-openmp -make -make install -make distclean diff --git a/modules/functions.cmake b/modules/functions.cmake new file mode 100644 index 0000000..51cac09 --- /dev/null +++ b/modules/functions.cmake @@ -0,0 +1,43 @@ +# Locate CFITSIO using pkg-config or use command line arguments to configure +# CFITSIO +function(find_cfitsio) + if(CFITSIO_LIBRARIES STREQUAL "" OR NOT DEFINED CFITSIO_LIBRARIES OR + CFITSIO_LIBRARY_DIRS STREQUAL "" OR NOT DEFINED CFITSIO_LIBRARY_DIRS OR + CFITSIO_INCLUDE_DIRS STREQUAL "" OR NOT DEFINED CFITSIO_INCLUDE_DIRS) + pkg_check_modules(CFITSIO REQUIRED cfitsio) + else() + message(STATUS "Use manually configured cfitsio") + message(STATUS " includes: ${CFITSIO_INCLUDE_DIRS}") + message(STATUS " libs: ${CFITSIO_LIBRARY_DIRS}") + message(STATUS " flags: ${CFITSIO_LIBRARIES}") + endif() + include_directories(${CFITSIO_INCLUDE_DIRS}) + link_directories(${CFITSIO_LIBRARY_DIRS}) +endfunction() + +# Extract target names from source files +function(find_targets targets target_path ext) + file(GLOB src_targets "${PROJECT_SOURCE_DIR}/${target_path}/*.${ext}") + list(TRANSFORM src_targets REPLACE "${PROJECT_SOURCE_DIR}/${target_path}/" "") + list(TRANSFORM src_targets REPLACE ".${ext}" "") + set(${targets} "${src_targets}" PARENT_SCOPE) +endfunction() + +# Build library +function(build_lib library) + file(GLOB src_${library} "${PROJECT_SOURCE_DIR}/src/lib${library}/*.cc") + file(GLOB inc_${library} "${PROJECT_SOURCE_DIR}/src/lib${library}/*.h") + include_directories("${PROJECT_SOURCE_DIR}/src/lib${library}") + add_library(${library} STATIC ${src_${library}}) + target_link_libraries(${library} ${CFITSIO_LIBRARIES} ${FFTW_LD_FLAGS}) + if(BUILD_FFTW) + add_dependencies(${library} fftw) + endif() + INSTALL(FILES ${inc_${library}} DESTINATION include/sparse2d) +endfunction() + +# Build binary +function(build_bin program libs target_path ext) + add_executable(${program} "${PROJECT_SOURCE_DIR}/${target_path}/${program}.${ext}") + target_link_libraries(${program} ${CFITSIO_LIBRARIES} ${FFTW_LD_FLAGS} ${libs}) +endfunction() diff --git a/modules/nfft.cmd b/modules/nfft.cmd deleted file mode 100755 index f5197eb..0000000 --- a/modules/nfft.cmd +++ /dev/null @@ -1,14 +0,0 @@ -#! /bin/bash -dest_dir=`pwd` -alias libtoolize=$(type -p glibtoolize libtoolize | head -1) -touch ChangeLog -rm -rf autom4te.cache -libtoolize -autoreconf --verbose --install --force -autoreconf --verbose --install --force -autoreconf --verbose --install --force -rm -f config.cache -./configure --prefix=$dest_dir --enable-openmp -make -j8 -make install -make distclean diff --git a/modules/nfft_fftw.cmd b/modules/nfft_fftw.cmd deleted file mode 100755 index 5a15ec5..0000000 --- a/modules/nfft_fftw.cmd +++ /dev/null @@ -1,15 +0,0 @@ -#! /bin/bash -dest_dir=`pwd` -alias libtoolize=$(type -p glibtoolize libtoolize | head -1) -touch ChangeLog -rm -rf autom4te.cache -libtoolize -autoreconf --verbose --install --force -autoreconf --verbose --install --force -autoreconf --verbose --install --force -rm -f config.cache -./configure --prefix=$dest_dir --with-fftw3=$dest_dir --enable-openmp -make -j8 -make install -make distclean -rm include/config.h.in~ diff --git a/src/libsparse1d/SB_Filter1D.h b/src/libsparse1d/SB_Filter1D.h index 10a5c4e..ab52029 100644 --- a/src/libsparse1d/SB_Filter1D.h +++ b/src/libsparse1d/SB_Filter1D.h @@ -8,22 +8,22 @@ ** ** Author: Jean-Luc Starck ** -** Date: 5/03/2000 -** +** Date: 5/03/2000 +** ** File: SB_Filter1D.h ** ** Modification history : ** ******************************************************************************* ** -** DESCRIPTION -** ----------- -** -** PARAMETRES -** ---------- -** -** RESULTS -** ------- +** DESCRIPTION +** ----------- +** +** PARAMETRES +** ---------- +** +** RESULTS +** ------- ** ** ******************************************************************************/ @@ -55,7 +55,7 @@ class SubBand1D { Bool SubSample_H_Even; Bool SubSample_G_Odd; - int DistPix; // distance between two adjacent pixels + int DistPix; // distance between two adjacent pixels // normally defaulted to 1. // only used with the decimated transform SubBand1D (){setBorder(I_MIRROR);DistPix=1; @@ -68,7 +68,7 @@ class SubBand1D { case I_CONT: test_index_function = test_index_cont; break; case I_MIRROR: test_index_function = test_index_mirror; break; case I_PERIOD: test_index_function = test_index_period; break; - case I_ZERO: + case I_ZERO: default: test_index_function = test_index_zero; break; @@ -77,27 +77,27 @@ class SubBand1D { inline int test_index(int i, int N) { return (*test_index_function)( i, N ); - } - + } + virtual void transform (int , float *, float *, float *){}; // (int N, float *High, float *Low, float *Det) // sub-band decomposiiton with decimation // N = number of pixels in the input signal High // High is decomposed in two sub signals Low and Det // Low = low resolution signal: its size is (N+1)/2 - // Det = Detail information: its size is N /2 - + // Det = Detail information: its size is N /2 + virtual void recons (int , float *, float *, float *){}; // (int N, float *Low, float *Det, float *High) // sub-band reconstruction with decimation // inverse transform: the output High (size N) is reconstructed // from Low (size (N+1)/2) and Det (size N/2) - - + + virtual void noise_transform (int , float *, float *, float *){}; - // (int N, float *High, float *Low, float *Det) - - + // (int N, float *High, float *Low, float *Det) + + virtual void transform (int , float *, float *, float *, int){}; // (int N, float *High, float *Low, float *Det, int Step) // sub-band decomposiiton without decimation @@ -105,19 +105,21 @@ class SubBand1D { // High is decomposed in two sub signals Low and Det // Low = low resolution signal: its size is N // Det = Detail information: its size is N - + virtual void recons (int , float *, float *, float *, int){}; // (int N, float *Low, float *Det, float *High, int Step) // sub-band reconstruction without decimation // inverse transform: the output High (size N) is reconstructed // from Low (size N) and Det (size N) - + virtual ~SubBand1D(){}; + type_border Border; // Border management type + protected: TestIndexFunction test_index_function; - type_border Border; // Border management type + }; /*********************************************************************** @@ -127,12 +129,12 @@ class SubBand1D { class SubBandFilter: public SubBand1D { - int NormCoef; // Norm = 2 for L1 normalization + int NormCoef; // Norm = 2 for L1 normalization // Norm = 1 for L2 normalization type_sb_filter TypeFilter; // type of bands float *H0,*G0,*H1,*G1; // filter banks int Size_H0, Size_H1, Size_G0, Size_G1; // filter sizes - int Start_H0, Start_H1, Start_G0, Start_G1; // First position of the + int Start_H0, Start_H1, Start_G0, Start_G1; // First position of the // filter sb_type_norm TypeNorm; void reset_param(); @@ -153,7 +155,7 @@ class SubBandFilter: public SubBand1D // Decimated SubBand decomposition void convol_h0 (int N, float *Input, float *Output); // convolved a signal by H0 filter - void convol_g0 (int N, float *Input, float *Output); + void convol_g0 (int N, float *Input, float *Output); // convolved a signal by G0 filter void noise_convol_h0 (int N, float *Input, float *Output); // convolved a signal by H0^2 filter @@ -170,12 +172,12 @@ class SubBandFilter: public SubBand1D // SignalOut = output smooth signal of size (N+1)/2 // DetailOut = output wavelet coefficient of size N/2 // the DistPix field is used in order to know the distance to use - // between two consecutive pixels. i.e if a given scale + // between two consecutive pixels. i.e if a given scale // has been undecimated, the pixel to consider must separated by a // a distance of 2. void recons (int N, float *SignalIn, float *DetailIn, float *SignalOut); - // reconstruct a signal for the low resolution part and the + // reconstruct a signal for the low resolution part and the // associated wavelet coefficient // N = ouput signal size // SignalIn = input signal of size(N+1)/2 @@ -188,7 +190,7 @@ class SubBandFilter: public SubBand1D // DetailOut = output signal coefficient of size N/2 = SignalIn convolved with G0^2 // Undecimated SubBand decomposition - void convol_filter(int N, float *Input, float *Output, + void convol_filter(int N, float *Input, float *Output, float *F, int SizeFilter, int Start_Filter, int Step); // convolution routine used for the transformation // Input and Output signals have the same size N @@ -197,8 +199,8 @@ class SubBandFilter: public SubBand1D // Start_Filter = low range index of the filter // Step = distance between two cofficients ("trous" size). - void rec_convol_filter(int N, float *Input, float *Output, - float *F, int SizeFilter, int Start_Filter, int Step); + void rec_convol_filter(int N, float *Input, float *Output, + float *F, int SizeFilter, int Start_Filter, int Step); // convolution routine used for the reconstruction // Input and Output signals have the same size N // F = pointer to filter value @@ -214,26 +216,26 @@ class SubBandFilter: public SubBand1D // Step = distance between two cofficients ("trous" size). void recons (int N, float *SignalIn, float *DetailIn, float *SignalOut, int Step); - // reconstruct a signal for the low undecimated resolution part and the + // reconstruct a signal for the low undecimated resolution part and the // associated unbdecimated wavelet coefficient // N = input-ouput signal size // SignalIn = input signal of size N // DetailIn = input wavelet coefficient N // SignalOut = outpout signal of size N - + void transform (fltarray& Sig_in, fltarray& Sig_out, fltarray* Smooth=NULL); //decimated subband transform // Sig_in = input signal of size N // Sig_out= output detail signal of size N // Smooth= output smooth signal of size N - + void recons (fltarray& Sig_in, fltarray& Sig_out, fltarray& Smooth); - // reconstruct a signal for the low decimated resolution part and the + // reconstruct a signal for the low decimated resolution part and the // associated decimated wavelet coefficient // Sig_in = input signal of size N // Sig_out = outpout signal of size N // Smooth = input smooth wavelet coefficient N - + ~SubBandFilter(); }; @@ -246,7 +248,7 @@ class G_Minmax: public SubBand1D { public: void transform(int N, float *High, float *Low, float *Det); void recons(int N, float *Low, float *Det, float *High); - + void transform(int N, float *High, float *Low, float *Det, int Step); void recons(int N, float *Low, float *Det, float *High, int Step); G_Minmax(){Border = I_CONT;} @@ -277,7 +279,7 @@ class Lifting: public SubBand1D { void recons(int N, float *Low, float *Det, float *High); void transform(int N, float *High, float *Low, float *Det, int Step); void recons(int N, float *Low, float *Det, float *High, int Step); - + type_lift TypeTrans; Lifting(){TypeTrans=DEF_LIFT;} Lifting(type_lift TL) {TypeTrans=TL;} @@ -289,7 +291,7 @@ class Lifting: public SubBand1D { ************************************************************************/ // U_B3SPLINE: H = B3spline, G=Id-H, Ht=H, Gt=Id+H -// U_B3SPLINE_2: H = B3spline, G=Id-H, Gt=G, HT =2*Id-H +// U_B3SPLINE_2: H = B3spline, G=Id-H, Gt=G, HT =2*Id-H #define NBR_UNDEC_FILTER 4 enum type_undec_filter {U_B3SPLINE,U_B3SPLINE_2,U_B2SPLINE,U_HAAR_B3S_POS,U_HAAR_B3S}; @@ -332,7 +334,7 @@ class UndecSubBandFilter: public SubBand1D // inverse transform: the output High (size N) is reconstructed // from Low (size N) and Det (size N) - ~UndecSubBandFilter() { if (FilterH != NULL) delete [] FilterH; FilterH=NULL; + ~UndecSubBandFilter() { if (FilterH != NULL) delete [] FilterH; FilterH=NULL; if (FilterG != NULL) delete [] FilterG; FilterG=NULL; if (FilterTG != NULL) delete [] FilterTG; FilterTG=NULL; if ((TypeUnderFilter == U_HAAR_B3S) && (FilterTH != NULL)) delete [] FilterTH; @@ -362,7 +364,7 @@ class HALF_1D_WT { int NbrUndecimatedScale; public: HALF_1D_WT() {Ptr_SB1D = NULL; NbrUndecimatedScale=0;} - HALF_1D_WT(SubBand1D &SB1D, int NbrUndec) + HALF_1D_WT(SubBand1D &SB1D, int NbrUndec) {Ptr_SB1D = &SB1D;NbrUndecimatedScale=NbrUndec;} void alloc(SubBand1D &SB1D, int NbrUndec) {Ptr_SB1D = &SB1D;NbrUndecimatedScale=NbrUndec;} @@ -379,12 +381,12 @@ class PAVE_1D_WT { SubBand1D *Ptr_SB1D; public: PAVE_1D_WT(SubBand1D &SB1D) {Ptr_SB1D = &SB1D;}; - float getAbsMaxTransf (fltarray & WT_Trans, float SigmaNoise, Bool OnlyPositivDetect=False, Bool Verbose=False); + float getAbsMaxTransf (fltarray & WT_Trans, float SigmaNoise, Bool OnlyPositivDetect=False, Bool Verbose=False); void transform (fltarray &SignalIn, fltarray &Transf_out, int Nbr_Plan); - void transform (fltarray &SignalIn, fltarray &Transf_out, + void transform (fltarray &SignalIn, fltarray &Transf_out, fltarray &Smooth, int step); void recons (fltarray &Transf_in, fltarray &Sig_Out, int Nbr_Plan); - void recons (fltarray &Transf_in, fltarray &Sig_Out, + void recons (fltarray &Transf_in, fltarray &Sig_Out, fltarray &Smooth, int step); ~PAVE_1D_WT(){Ptr_SB1D = NULL;} }; @@ -397,17 +399,17 @@ class PAVE_1D_WT { class HALF_DECIMATED_1D_WT { private: SubBandFilter *Ptr_SB1D; - + void set_tabdec(int NumUndec, Bool * & TabDec, int Nbr_Plan); void transform (fltarray& Signal, fltarray* TabTrans, Bool *TabDec); float update (float CoefSol, float Threshold, float SoftLevel); bool is_on_support (float CoefSol, float Threshold, float SoftLevel); void recons (fltarray* TabTrans, fltarray& Signal, Bool *TabDec); - + public: int NbScale; int NbrUndecimatedScale; - type_border Bord; + type_border Bord; float SigmaNoise; bool OnlyPositivDetect; bool SuppressIsolatedPixel; @@ -415,13 +417,13 @@ class HALF_DECIMATED_1D_WT { bool Verbose; bool Write; bool UseNormL1; - float NSigmaSoft; - + float NSigmaSoft; + HALF_DECIMATED_1D_WT (SubBandFilter &SB1D) {Ptr_SB1D = &SB1D;} int alloc (fltarray* & TabTrans, int Nx, int Nbr_Plan, Bool *TabDec); int alloc (fltarray* & TabTrans, int Nx, int Nbr_Plan, int NbrUndecimatedScale); void free (fltarray* TabTrans, int Nbr_Plan); - + void transform (fltarray& Signal, fltarray* TabTrans); void KillScaleNotUsed (fltarray* WT_Trans, int FirstDetectScale); void KillLastScale (fltarray* WT_Trans); diff --git a/src/libsparse2d/MW_Deconv.cc b/src/libsparse2d/MW_Deconv.cc new file mode 100644 index 0000000..62247eb --- /dev/null +++ b/src/libsparse2d/MW_Deconv.cc @@ -0,0 +1,687 @@ +/****************************************************************************** +** Copyright (C) 1998 by CEA +******************************************************************************* +** +** UNIT +** +** Version: 1.0 +** +** Author: Jean-Luc Starck +** +** Date: 23/12/98 +** +** File: MW_Deconv.cc +** +******************************************************************************* +** +** DESCRIPTION Image Deconvolution using the multiscale entropy +** ----------- +** +** +******************************************************************************/ + +#include "IM_Obj.h" +#include "IM_Math.h" +// #include "IM_IO.h" +#include "MR_Obj.h" +#include "MR_NoiseModel.h" +#include "CErf.h" +#include "CMem.h" +#include "IM_Deconv.h" +#include "MR_Deconv.h" +#include "MW_Deconv.h" +#include "MW_Filter.h" + + +/*****************************************************************************/ + +void MWDeconv::mw_grad(Ifloat & Gradient) +// Apply the multiscale entropy to the multiresolution data set MR_Obj +{ + float Alpha=0.,Alphab; + float Sigma; + int Nbr_band = MR_Obj.nbr_band(); + int b,i,j; + + // apply the multiscale entropy on each band + MR_Obj.transform(Obj); + + for (b=0; b < Nbr_band-1; b++) + { + // for all wavelet coefficients + for (i=0; i < MR_Obj.size_band_nl(b); i++) + for (j=0; j < MR_Obj.size_band_nc(b); j++) + { + float Val = MR_Obj(b,i,j); + Sigma = ModelData->sigma(b,i,j); + Alphab = RegulParam*TabAlpha(b); + //Modify Lambda using the SNR of the data + switch(TypeRegul) + { + case DEC_NO_REGUL: + Alpha = 0.; + break; + case DEC_REGUL_NO_PROTECT: + Alpha = Alphab; + break; + case DEC_REGUL_PROBA: + Alpha = Alphab*MR_Data(b,i,j); + break; + case DEC_REGUL_SUPPORT: + //if ((*ModelData)(b,i,j) == True) Alpha = 0.; + //else + Alpha = Alphab* MR_RegulObj(b,i,j); + break; + case DEC_REGUL_PROBA_SUPPORT: + //if ((*ModelData)(b,i,j) == True) Alpha = 0.; + //else Alpha = Alphab*MR_Data(b,i,j); + Alpha = Alphab * MR_Data(b,i,j) * MR_RegulObj(b,i,j); + break; + } + // correct the wavelet coefficients + if (DecMethod == DEC_MR_MEM_NOISE) + { + int Sign = (Val >= 0) ? 1: -1; + MR_Obj(b,i,j) = MemWObj.grad_hn(ABS(Val)/Sigma)*Sigma*Sign*Alpha; + } + else MR_Obj(b,i,j) = Val*Alpha; + } + + } + MR_Obj.band(Nbr_band-1).init(); + MR_Obj.recons(Buff); + Gradient -= Buff; + + +} + +/*********************************************************************/ + +float MWDeconv::mw_find_optim(Ifloat &Gradient) +// MR_Obj contains the wavelet transform of the current solution +{ + int b,i,j; + int Nbr_band = MR_Obj.nbr_band(); + float hcroise,hgrad, coef_num,coef_den, Wg, ValRet; + + MultiResol MR_Grad(Resi.nl(), Resi.nc(), MR_Obj.nbr_scale(), + MR_Obj.Type_Transform, "MR_Grad"); + MR_Grad.transform(Gradient); + MR_Obj.transform(Obj); + psf_convol (Gradient, Psf_cf, Buff); + coef_den = energy(Buff); + Buff *= Resi; + coef_num = flux(Buff); + + hcroise=0; + hgrad=0; + for (b = 0; b < Nbr_band-1; b++) + for (i = 0; i < MR_Grad.size_band_nl(b); i++) + for (j = 0; j < MR_Grad.size_band_nc(b); j++) + { + Wg = (TypeRegul == DEC_REGUL_PROBA) ? MR_Data(b,i,j): 1.; + Wg *= TabAlpha(b)*RegulParam; + hcroise += MR_Obj(b,i,j) * MR_Grad(b,i,j) * Wg; + hgrad += MR_Grad(b,i,j) * MR_Grad(b,i,j)* Wg; + } + coef_num -= hcroise; + coef_den += hgrad; + cout << "Coef num = " << coef_num << " Den = " << coef_den; + cout << " param = " << coef_num / coef_den << endl; + ValRet = coef_num / coef_den; + return ValRet; +} + + +/***************************************************************************/ + +void MWDeconv::compute_mem_resi() +{ + int b,i,j; + + psf_convol (Obj, Psf_cf, Resi); + for (i=0; i< Nl*Nc; i++) Resi(i) = Imag(i) - Resi(i); + SigmaResi = sigma(Resi); + + switch (TypeWeight) + { + case DEC_WEIGHT_PROBA: + MR_Obj.transform(Resi); + for (b = 0; b < NbrBand-1; b++) + for (i = 0; i < MR_Obj.size_band_nl(b); i++) + for (j = 0; j < MR_Obj.size_band_nc(b); j++) + MR_Obj(b,i,j) *= (1. - MR_Data(b,i,j)); + MR_Obj.recons(Resi); + break; + case DEC_WEIGHT_SUPPORT: + MR_Obj.transform(Resi); + ModelData->threshold(MR_Obj); + MR_Obj.recons(Resi); + break; + case DEC_NO_WEIGHT: + break; + default: + cerr << "Error: unknown data weighting method ... " << endl; + exit(-1); + break; + } +} + +/****************************************************************************/ + +float MWDeconv::sigma_obj(float Sigma) +// The routine estimate the regularization paramter per band +// 1) Create a random Gaussian noise: N +// 2) Convolve it by the PSF: I = P*N +// 3) Take its wavelet transform: W = WT(I) +// 4) Calculate the standard deviation in each band j: TabSigma[j] +// 5) MaxSigma = max(TabSigma) +// 6) The regularisation parameter at each band j is given by +// TabAlpha[j] = Sigma_j / TabSigma[i] +// 7) Renormalisation: TabAlpha[j] = TabAlpha[j] / MaxSigma +{ + int b; + Ifloat Dat(Nl, Nc, "noise"); + TabSigmaObj.alloc(NbrBand); + im_noise_gaussian (Dat, 1.); + psf_convol_conj (Dat, Psf_cf); + MR_Obj.transform(Dat); + + for (b = 0; b < NbrBand-1; b++) + TabSigmaObj (b) = Sigma * MR_Obj.band_norm(b) / sigma( MR_Obj.band(b)); + + float Max = TabSigmaObj.max(); + for (b = 0; b < NbrBand-1; b++) TabAlpha(b) = TabSigmaObj(b) / Max; + + if (Verbose == True) + { + for (b = 0; b < NbrBand-1; b++) + { + cout << "Band " << b << " alpha = " << TabAlpha(b) << endl; + // cout << " Sigma obj = " << TabSigmaObj(b) << " P2 = " + // << TabSigmaObj(b)*TabSigmaObj(b)<< endl; + } + } + return (Sigma / sigma(Dat)); +} + +/****************************************************************************/ + +void MWDeconv::get_mr_support_first_guess(type_deconv Deconv, int Niter) +{ + // save real parameters + type_deconv RealDecMethod = DecMethod; + int RealMaxIter = MaxIter; + float RealIterCvg = IterCvg; + Bool RealVerb = Verbose; + Bool RealOptim = OptimParam; + // int b,i,j; + //fltarray TabNsig[NbrBand]; + + // Change some paramters + IterCvg = 1.; + DecMethod = Deconv; + MaxIter = Niter; + Verbose = False; + OptimParam = False; + //float NSig = 5.; + //for (b=0; b < NbrBand; b++) TabNsig[b] = (ModelData->NSigma)[b]; + //for (b=0; b < NbrBand; b++) (ModelData->NSigma)[b] = NSig; + // run the deconvolution + im_iter_deconv(); + + // restore real parameters + DecMethod = RealDecMethod; + MaxIter = RealMaxIter; + IterCvg = RealIterCvg; + Verbose = RealVerb; + OptimParam = RealOptim; + //for (b=0; b < NbrBand; b++) (ModelData->NSigma)[b] = TabNsig[b]; + + if (Verbose == True) + { + INFO_X(Obj, "MR Van-Citter Sol: "); + } +} + +/****************************************************************************/ + +void MWDeconv::find_support_obj() +{ + int b,j,i; + + float Sig = sigma(Obj); + + MR_Obj.transform(Obj); + // io_write_ima_float("sol1.fits", Obj); + for (b = 0; b < NbrBand-1; b++) + { + int Scale = MR_Obj.band_to_scale(b); + int Band = (int) pow ((double) 2., (double)(Scale+1)); + + Buff.resize(MR_Obj.size_band_nl(b), MR_Obj.size_band_nc(b)); + float Level = NSigmaObj*Sig*MR_Obj.band_norm(b); + for (i = 0; i < MR_Obj.size_band_nl(b); i++) + for (j = 0; j < MR_Obj.size_band_nc(b); j++) + { + if (ABS(MR_Obj(b,i,j)) > Level) Buff(i,j) = 0.; + else Buff(i,j) = 1; + } + smooth_bspline(Buff, MR_RegulObj.band(b), I_ZERO, Scale); + for (i = Band; i < MR_Obj.size_band_nl(b)-Band; i++) + for (j = Band; j < MR_Obj.size_band_nc(b)-Band; j++) + if (ABS(MR_Obj(b,i,j)) > Level) MR_RegulObj(b,i,j) = 0.; + } + Buff.resize(Nl,Nc); + // MR_RegulObj.write("sup.mr"); +} + +/****************************************************************************/ + +void MWDeconv::mw_vaguelet(Bool UseSupport, Ifloat *FirstGuess, Ifloat *ICF) +{ + int b,i,j; + CMemWave MemWObj; + Bool RecEachIter = False; + Bool DataSNR = False; + + WaveFilterResi = True; + NbrBand = ModelData->nbr_band(); + Nscale = ModelData->nbr_scale(); + Transform= ModelData->type_trans(); + TabAlpha.alloc(NbrBand); + init_deconv(FirstGuess, ICF); + + Ifloat NoiseSimu(Nl,Nc,"simu"); + va_inverse_data(NoiseSimu); + // io_write_ima_float("xx_noise.fits", NoiseSimu); + + Resi = Obj; +// if (UseSupport == True) +// { +// MR_Data.transform(Obj); +// if (MR_Data.Set_Transform != TRANSF_MALLAT) +// { +// type_filter Filter = (MR_Data.Set_Transform == TRANSF_MALLAT) ? +// FILTER_THRESHOLD: FILTER_ITER_THRESHOLD; +// MRFiltering CFilter(*ModelData, Filter); +// NoiseSimu=Obj; +// CFilter.Sup_Set = False; +// CFilter.Verbose = Verbose; +// CFilter.KillLastScale = KillLastScale; +// CFilter.mr_support_filter(NoiseSimu, MR_Data, Resi); +// } +// else +// { +// ModelData->threshold(MR_Data); +// MR_Data.recons(Resi); +// } +// Obj -= Resi; +// } + + // io_write_ima_float("xx_inv.fits", Obj); + for (b = 0; b < NbrBand-1; b++) + TabAlpha(b) = (RegulParam <= 0) ? 1.: RegulParam; + + // MR_Data.transform(Obj); + // apply the filtering + Bool DataEdge = False; + MultiResol *MR_Model = NULL; + float FilCvgParam= DEF_MEM_FILT_CONGER; + int FilNbrIter = DEF_MEM_FILT_MAX_ITER; + + MRNoiseModel MD(NOISE_CORREL, Nl, Nc, Nscale, Transform); + MD.RmsMap = NoiseSimu; + for (b=0; b < NbrBand; b++) MD.NSigma[b]= ModelData->NSigma[b]; + for (b=0; b < NbrBand; b++) MD.TabEps[b]= ModelData->TabEps[b]; + MD.OnlyPositivDetect = ModelData->OnlyPositivDetect; + MD.SupIsol = ModelData->SupIsol; + + MD.model(Obj, MR_Data); + // MD.threshold(MR_Data); + mw_fm3(MR_Data, MD, TabAlpha, MemWObj, DataSNR, + DataEdge, MR_Model, FilCvgParam, FilNbrIter, Verbose, + RecEachIter, PositivConstraint, UseSupport); + // MR_Data.write("xx.mr"); + int Lb = MR_Data.nbr_band()-1; + if (KillLastScale == True) MR_Data.band(Lb).init(); + + MR_Data.recons(Obj); + // Resi -= Obj; + // io_write_ima_float("xx_resi.fits", Resi); + // if (UseSupport == True) Obj += Resi; + + if (PositivConstraint == True) threshold(Obj); + if (KillLastScale == True) + { + MultiResol MR_Resi(Nl, Nc, MD.nbr_scale(), MD.type_trans(), "MR_Data"); + for (int Iter=0; Iter < 10; Iter++) + { + MR_Resi.transform(Obj); + for (b = 0; b < MR_Data.nbr_band()-1; b++) + for (i = 0; i < MR_Data.size_band_nl(b); i++) + for (j = 0; j < MR_Data.size_band_nc(b); j++) + MR_Resi(b,i,j) = MR_Data(b,i,j) - MR_Resi(b,i,j); + MR_Resi.band(Lb).init(); + MR_Resi.recons(Resi); + Obj += Resi; + if (PositivConstraint == True) threshold(Obj); + } + } + + compute_resi(); + if (GaussConv == True) convol_gauss(Obj, Fwhm); + else if (UseICF == True) psf_convol(Obj, Ima_ICF, Obj); +} + +/****************************************************************************/ + +void MWDeconv::mw_deconv(Ifloat *FirstGuess, Ifloat *ICF) +{ + int b,i,Iter=0; + Bool Stop = False; + float OldFonc=1e+15, Fonc, Delta; + float Cvgparam = IterCvg; +// if ((OptimParam == True) && (DecMethod == DEC_MR_MEM_NOISE)) +// { +// cerr << "Optimization is not implemented when using noise entropy ... " << endl; +// exit(-1); +// } + + + // initialization + // ---------------- + WaveFilterResi = True; + NbrBand = ModelData->nbr_band(); + Nscale = ModelData->nbr_scale(); + Transform= ModelData->type_trans(); + TabAlpha.alloc(NbrBand); + // Border = I_CONT; + init_deconv(FirstGuess,ICF); + NbIter_SupportCalc = DEF_NB_ITER_SUP_CAL; + if ((RegulParam <= 0.) || (TypeRegul == DEC_NO_REGUL)) + { + RegulParam = 0.; + TypeRegul = DEC_NO_REGUL; + } + + + int Np = Nl*Nc; + Ifloat Gradient(Nl, Nc, "Info"); + Buff.alloc(Nl, Nc, "Buffer"); + MR_Obj.alloc(Nl, Nc, Nscale, Transform, "MR_Obj"); + MR_Obj.Border = Border; + if ((TypeRegul == DEC_REGUL_PROBA_SUPPORT) || + (TypeRegul == DEC_REGUL_SUPPORT)) + { + MR_RegulObj.alloc(Nl, Nc, Nscale, Transform, "MR_Obj"); + MR_RegulObj.Border = Border; + } + Noise_Ima = ModelData->SigmaNoise; + + // Noise_Obj = sigma_obj(Noise_Ima); + for (b = 0; b < NbrBand-1; b++) TabAlpha(b) = 1.; + float FluxImag=0.; + + if (NormFlux == True) FluxImag = flux(Imag); + + // Find the first guess using the regularized Van Cittert Method + if (GetMR_FirstGuess == True) get_mr_support_first_guess(DEC_MR_GRADIENT,50); + + if ((TypeWeight == DEC_WEIGHT_PROBA) || (TypeRegul == DEC_REGUL_PROBA) + || (TypeRegul == DEC_REGUL_PROBA_SUPPORT)) + { + ModelData->prob_noise(MR_Data); + } + + + if (Verbose == True) + { + cout << "Image size = " << Nl << " " << Nc << endl; + cout << "Noise_Ima = " << Noise_Ima << endl; + // cout << "Sigma obj = " << Noise_Obj << endl; + } + + // Residual estimation + compute_mem_resi(); + + // if (Verbose == True) cout << "start iterate" << endl; + do + { + // Calculate the multiplicate term in Gradient's iteration + // io_write_ima_float("xx_resi.fits", Resi); + psf_convol_conj (Resi, Psf_cf, Gradient); + + // calculates the gradient + if (TypeRegul != DEC_NO_REGUL) mw_grad(Gradient); + + // Functional gradient + //io_write_ima_float("xx_grad.fits", Gradient); + //if ((Verbose == True) && (Iter == MaxIter-1)) + //{ + // i = 670; + // cout << "P(10,10) = " << Gradient(i) << endl; + //} + //io_write_ima_float("xx_fgrad.fits", Gradient); + + // Supress the last scale in the Gradient + if (KillLastScale == True) + { + int Nbr_band = MR_Obj.nbr_band(); + MR_Obj.transform(Gradient); + MR_Obj.band(Nbr_band-1).init(); + MR_Obj.recons(Gradient); + } + + // convergence parameter + if (OptimParam == True) + { + // Cvgparam = mw_find_optim(Gradient); + // PB: mw_find_optim does not work correctly + // it replaces by a X2 optimization + + // Xi2 optimization + psf_convol(Gradient, Psf_cf, Buff); + Cvgparam = flux(Buff*Resi)/energy(Buff); + if (Cvgparam > 10) Cvgparam = 10.; + } + + // calculate the next object estimation with positivity constraint + for (i = 0; i < Np; i++) Obj(i) += Cvgparam * Gradient(i); + + if (PositivConstraint == True) threshold(Obj); + if (NormFlux == True) + { + float FluxObj = flux(Obj); + FluxObj = FluxImag / FluxObj; + for (i= 0; i < Np; i++) Obj(i) *= FluxObj; + } + + // New Residual estimation + compute_mem_resi(); + + Iter ++; + Fonc = SigmaResi / Noise_Ima; + Delta = (Iter > 1) ? OldFonc - Fonc : Fonc; + if (Iter > MaxIter) Stop = True; + else if ((ABS(Delta) < EpsCvg) && (EpsCvg > 0)) Stop = True; + OldFonc = Fonc; + + if ((Verbose == True) || (Stop == True)) + { + if (Iter % 2 == 0) + { +// cout.width(8); // see p343 c++ stroustrup +// cout.fill(' '); +// cout.setf(ios::right,ios::adjustfield); +// cout.setf(ios::scientific,ios::floatfield); +// cout.precision(4); + if (OptimParam == True) + printf("%d: Sigma(grad) = %f, Delta = %f, Sigma(resi) = %f, Cvgparam = %f \n", + Iter, Fonc, Delta, SigmaResi, Cvgparam); + else printf("%d: Sigma(grad) = %f, Delta = %f, Sigma(resi) = %f \n", + Iter, Fonc, Delta, SigmaResi); + } + } + + // Get new multiresolution support + if ((Iter % NbIter_SupportCalc == 0) && + ((TypeRegul == DEC_REGUL_SUPPORT) || + (TypeRegul == DEC_REGUL_PROBA_SUPPORT))) + { + find_support_obj (); + } + + } while (Stop == False); + + psf_convol (Obj, Psf_cf, Resi); + for (i=0; i< Nl*Nc; i++) Resi(i) = Imag(i) - Resi(i); + + if (GaussConv == True) convol_gauss(Obj, Fwhm); + else if (UseICF == True) psf_convol(Obj, Ima_ICF, Obj); + + if (Verbose == True) + { + INFO_X(Obj, "Solution: "); + INFO_X(Resi, "Resi: "); + } +} + +/*****************************************************************************/ + +inline float laplacian_val(Ifloat &Obj, int i, int j) +{ + type_border type=I_CONT; + float Val = - Obj(i,j) + 0.25*( Obj(i+1,j,type) + + Obj(i-1,j,type) + Obj(i,j+1,type) + Obj(i,j-1,type)); + return Val; +} + +/*****************************************************************************/ + +// static void laplacien_regul(Ifloat &Obj, Ifloat &Grad, float Regul) +// { +// int Nl = Obj.nl(); +// int Nc = Obj.nc(); +// int i,j; +// Ifloat Aux(Nl,Nc,"Aux"); +// +// for (i = 0; i < Nl; i ++) +// for (j = 0; j < Nc; j ++) Aux(i,j) = laplacian_val(Obj,i,j); +// +// for (i = 0; i < Nl; i ++) +// for (j = 0; j < Nc; j ++) +// { +// Grad(i,j) -= Regul*laplacian_val(Aux,i,j); +// } +// } + +/****************************************************************************/ +/* +// TRY to build two solutions: one free of noise and the second +// containing the part contminated by the noise +void MWDeconv::mw_deconv(Ifloat *FirstGuess, Ifloat *ICF) +{ + Ifloat Gradient(Nl, Nc, "Info"); + float Func, OldFunc=1e10, Sigma, Delta; + int i,Iter=0; + Bool Stop = False; + float Cvgparam = IterCvg; + + // initialization + // ---------------- + WaveFilterResi = True; + NbrBand = ModelData->nbr_band(); + Nscale = ModelData->nbr_scale(); + Transform= ModelData->type_trans(); + TabAlpha.alloc(NbrBand); + // Border = I_CONT; + init_deconv(FirstGuess,ICF); + int Np = Nl*Nc; + NbIter_SupportCalc = DEF_NB_ITER_SUP_CAL; + if ((RegulParam <= 0.) || (TypeRegul == DEC_NO_REGUL)) + { + RegulParam = 0.; + TypeRegul = DEC_NO_REGUL; + } + + Buff.alloc(Nl, Nc, "Buffer"); + Noise_Ima = ModelData->SigmaNoise; + float Sigma2N = Noise_Ima*Noise_Ima*Np; + if (NormFlux == True) FluxImag = flux(Imag); + else FluxImag=0.; + + // Residual estimation + compute_resi(); + // INFO(Resi, "RESI"); + // start iterations + // cout << "Start iterate " << endl; + // INFO(Obj, "obj"); + ObjN.alloc(Nl,Nc,"objn"); + ObjS.alloc(Nl,Nc,"objn"); + cout << "Cvgparam = " << Cvgparam<< endl; + do + { + // compute the gradient + Gradient = Resi; + MR_Data.transform (Resi); + ModelData->threshold(MR_Data); + MR_Data.recons (Resi); + Gradient -= Resi; + //INFO(Resi, "Resi sig"); + //INFO(Gradient, "Resi non sig"); + psf_convol_conj (Resi, Psf_cf); + //INFO(Resi, "Resi Conv"); + for (i = 0; i < Np; i++) ObjS(i) += Cvgparam *Resi(i); + //INFO(ObjS, "ObjS"); + + psf_convol_conj (Gradient, Psf_cf, Resi); + if (RegulParam > 0) + { + laplacien_regul(ObjN, Resi, RegulParam); + for (i = 0; i < Np; i++) ObjN(i) += Cvgparam*Resi(i); + // INFO(ObjN, "ObjN"); + for (i = 0; i < Np; i++) Obj(i) = ObjS(i) + ObjN(i); + } + else Obj = ObjS; + + if (PositivConstraint == True) threshold(Obj); + if (NormFlux == True) + { + float FluxObj = flux(Obj); + FluxObj = FluxImag / FluxObj; + for (i= 0; i < Np; i++) Obj(i) *= FluxObj; + } + + // next residual + compute_resi(); + Sigma = sigma(Resi) / Sigma2N; + Func = Sigma / Noise_Ima; + Delta = (Iter > 1) ? OldFunc - Func : Func; + if (Iter > MaxIter) Stop = True; + else if ((ABS(Delta) < EpsCvg) && (EpsCvg > 0)) Stop = True; + OldFunc = Func; + + if ((Verbose == True) && ((Iter >= 0) && (Iter % 1 == 0))) + { + Func = fonctional() / (float) Np; + printf("%d: Sigma = %f, Func = %f\n", Iter, Sigma, Func); + } + + Iter ++; + if (Iter > MaxIter) Stop = True; + else if (EpsCvg > FLOAT_EPSILON) + { + if (OldFunc < Sigma) Stop = True; + if ((ABS(OldFunc - Sigma)) < EpsCvg) Stop = True; + } + OldFunc = Sigma; + } while (Stop == False); + + if (Verbose == True) + { + INFO(Obj, "Solution: "); + INFO(Resi, "Resi: "); + } +} +*/ diff --git a/src/libsparse2d/MW_Deconv.h b/src/libsparse2d/MW_Deconv.h new file mode 100644 index 0000000..5a68e51 --- /dev/null +++ b/src/libsparse2d/MW_Deconv.h @@ -0,0 +1,122 @@ +/****************************************************************************** +** Copyright (C) 1998 by CEA +******************************************************************************* +** +** UNIT +** +** Version: 1.0 +** +** Author: Jean-Luc Starck +** +** Date: 23/12/98 +** +** File: MW_Deconv.h +** +******************************************************************************* +** +** DESCRIPTION +** ----------- +** +** PARAMETRES +** ---------- +** +** +** RESULTS +** ------- +** +** +******************************************************************************/ + + +#ifndef __MWDECONV__ +#define __MWDECONV__ + +const int DEF_DEC_MAX_ITER = 500; +const int DEF_NB_ITER_SUP_CAL = 5; + +// different type of data weighting +const int DEC_NBR_WEIGHT = 3; +enum type_mem_weight {DEC_NO_WEIGHT, DEC_WEIGHT_PROBA, DEC_WEIGHT_SUPPORT}; +const type_mem_weight DEF_DEC_WEIGHT = DEC_WEIGHT_SUPPORT; + +// different type of regularization +const int DEC_NBR_REGUL = 5; +enum type_mem_regul {DEC_NO_REGUL, DEC_REGUL_NO_PROTECT, DEC_REGUL_PROBA, + DEC_REGUL_SUPPORT,DEC_REGUL_PROBA_SUPPORT}; + +const type_mem_regul DEF_DEC_REGUL = DEC_REGUL_SUPPORT; + +const type_deconv DEF_DEC_ENTROP = DEC_MR_MEM; +const type_deconv DEF_DEC_ENTROP_GAUSS = DEC_MR_MEM_NOISE; + +class MWDeconv: public MRDeconv +{ + float sigma_obj(float SigmaNoise); + // create a noisy map of standard deviation SigmaNoise + // convolve it by the PSF + // return the standard deviation of convolved map + + void init_mem_param() + { + RegulParam = 0; + MaxIter = DEF_DEC_MAX_ITER; + DecMethod = DEF_DEC_ENTROP; + Nscale = 0; + NbrBand = 0; + Transform = DEFAULT_TRANSFORM; + TypeWeight = DEC_NO_WEIGHT; + TypeRegul = DEC_REGUL_SUPPORT; + WaveFilterResi=True; + GetMR_FirstGuess=True; + } + + // Apply first another regularized iterative method + // in order to improve the convergence. + void get_mr_support_first_guess(type_deconv Deconv=DEC_MR_CITTERT, int Niter=10); + + CMemWave MemWObj; // Gradient class + void compute_mem_resi(); // compute the residual + void mw_grad(Ifloat & Gradient); // compute the gradient + float mw_find_optim(Ifloat &Gradient); // optimize the convergence paramter + void find_support_obj(); // find the multiresolution support + // of the object + + MultiResol MR_Obj; // Multiresolution buffer used in some calculation + MultiResol MR_RegulObj; // Multiresolution buffer used in some calculation + int NbIter_SupportCalc; // Number of iterations between two reestimations + // of the multiresolution support (def is 5). + Ifloat Buff; // Buffer image used during some calcualtion + Ifloat ObjN; + Ifloat ObjS; + float Noise_Obj; // Noise in the object space + fltarray TabSigmaObj; // object noise estimation per scale + fltarray TabAlpha; // Regularization parameter per scale + + int Nscale; + int NbrBand; + type_transform Transform; // type of transform + float SigmaResi ; // residual standard deviation + + public: + type_mem_regul TypeRegul; // type of regularization weighting + type_mem_weight TypeWeight; // type of weighting applied to the data + Bool GetMR_FirstGuess; // regularized Van Cittert method is + // first apply to improve the convergence + float NSigmaObj; // Nsigma for object multiresolution + // support detection + + MWDeconv(): MRDeconv() { init_mem_param();} + + void mw_deconv(Ifloat *FirstGuess, Ifloat *ICF); + // deconvolution by the multiscale entropy method + + void mw_vaguelet(Bool DataSNR, Ifloat *FirstGuess, Ifloat *ICF); + + ~MWDeconv(){}; + }; + +#endif + + + + diff --git a/src/libsparse2d/MW_Filter.cc b/src/libsparse2d/MW_Filter.cc new file mode 100644 index 0000000..471d74a --- /dev/null +++ b/src/libsparse2d/MW_Filter.cc @@ -0,0 +1,518 @@ +/****************************************************************************** +** Copyright (C) 1998 by CEA +******************************************************************************* +** +** UNIT +** +** Version: 1.0 +** +** Author: Jean-Luc Starck +** +** Date: 98/01/12 +** +** File: MW_Filter.cc +** +******************************************************************************* +** +** DESCRIPTION Image filtering using the multiscale entropy +** ----------- +** +** +******************************************************************************/ + +#include "IM_Obj.h" +// #include "IM_IO.h" +#include "MR_Obj.h" +#include "MR_NoiseModel.h" +#include "CErf.h" +#include "CMem.h" +#include "MW_Filter.h" + + +/****************************************************************************/ + +void mw_filter(MultiResol &MR_Data, MRNoiseModel & NoiseModel, + int TypeFilter, fltarray &TabAlpha, Bool DataSNR, + Bool DataModel, MultiResol *MR_Model, + float ConvgParam, int MaxIter, Bool PositivConst, Bool Verbose, Bool UseSupport) +/* Apply the multiscale entropy to the multiresolution data set MR_Data + MR_Data = in-out: multiresolution transform of an image + NoiseModel = in: noise modeling of the image in the wavelet space + TypeFilter = in: type of filtering + = DEF_MEM_ALPHA_CST ==> alpha is fixed + = DEF_MEM_ALPHA_OPT ==> alpha is globally optimized + = DEF_MEM_ALPHA_OPT_SCALE ==> alpha is optimized at each scale + + TabAlpha = in-out: TabAlpha(b) gives the alpha regularization parameter for + the band b + DataModel = in: true if we have a model + MR_Model = in: pointer to a multiresolution model or NULL (if none) + DataSNR = in: true if alpha parameter are modified using the SNR + of the data + UseSupport= in: true if alpha parameter are modified using the multiresolution + support + Warning: works well, when each wavelet can be modelized by a Gaussian +*/ +{ + CMemWave MemWObj; + + switch (TypeFilter) + { + case DEF_MEM_ALPHA_CST: + mw_fm1(MR_Data, NoiseModel,TabAlpha, MemWObj, DataModel, MR_Model, DataSNR, UseSupport); + break; + case DEF_MEM_ALPHA_OPT: + mw_fm2( MR_Data, NoiseModel, TabAlpha, MemWObj, DataSNR, + DataModel, MR_Model, ConvgParam, MaxIter, PositivConst, Verbose, UseSupport); + break; + case DEF_MEM_ALPHA_OPT_SCALE: + mw_fm3( MR_Data, NoiseModel, TabAlpha, MemWObj, DataSNR, + DataModel, MR_Model, ConvgParam, MaxIter, Verbose, + True, PositivConst, UseSupport); + break; + default: + cerr << "Error: unknow type of filtering ... " << endl; + exit(-1); + break; + } +} + +/****************************************************************************/ + +void mw_fm1(MultiResol & MR_Data, MRNoiseModel & NoiseModel, + fltarray &TabAlpha, CMemWave & MemWObj, + Bool UseModel, MultiResol *MR_Model, + Bool UseDataSNR, Bool UseSupport) +/* Apply the multiscale entropy to the multiresolution data set MR_Data + MR_Data = in-out: multiresolution transform of an image + NoiseModel = in: noise modeling of the image in the wavelet space + TabAlpha = in: TabAlpha(b) gives the alpha regularization parameter for + the band b + MemWObj = in: class for entropy filtering + UseModel = in: true if we have a model + MR_Model = in: pointer to a multiresolution model or NULL (if none) + UseDataSNR = in: true if alpha parameter are modified using the SNR + of the data + UseSupport= in: true if alpha parameter are modified using the multiresolution + support +*/ + +{ + float AlphaP,Val; + int b,i,j; + float Model=0.; + + // apply the multiscale entropy on each band + for (b = 0; b < MR_Data.nbr_band()-1; b++) + { + int s = MR_Data.band_to_scale(b); + + // for all wavelet coefficients + for (i = 0; i < MR_Data.size_band_nl(b); i++) + for (j = 0; j < MR_Data.size_band_nc(b); j++) + { + if (s < NoiseModel.FirstDectectScale) MR_Data(b, i,j) = 0; + else + { + float Sigma = NoiseModel.sigma(b,i,j); + float RegulParam = TabAlpha(b); + Val = MR_Data(b,i,j); + + // use the input model for the wavelet coefficients + if ((UseModel != True) || (MR_Model == NULL)) Model = 0.; + else Model = (*MR_Model)(b,i,j); + + // Modify Alpha using the SNR of the data + if ((UseDataSNR == True) + && ((NoiseModel.OnlyPositivDetect == False) || (Val >= 0))) + { + // SNR = MAX [ ABS(data) / (N*Sigma), 1 ] + AlphaP = ABS(Val / (Sigma* NoiseModel.NSigma[s])); + if (AlphaP > 1) AlphaP = 1; + + // if AlphaP = 0, J = h_n ==> Solution = model + // when RegulParam < 0, the filtering routine set to + // the model the solution + if (AlphaP < FLOAT_EPSILON) RegulParam = -1; + else RegulParam *= (1-AlphaP)/AlphaP; + } + else if ((UseSupport == True) + && (NoiseModel(b,i,j) == True)) RegulParam = 0; + + // correct the wavelet coefficient + if (RegulParam > 0) MR_Data(b, i,j) = MemWObj.filter(Val, RegulParam, Sigma, Model); + else MR_Data(b, i,j) = Val; + +// if ((b==0) && (i==59) && (j==47)) +// { +// if (NoiseModel(b,i,j) == True) cout << "Support coef ... " << endl; +// if (UseSupport == True) cout << "Use Support" << endl; +// cout << "Data = " << MR_Data(b, i,j) << " Corr = " << Val << endl; +// } + } + } + } +} + +/***************************************/ + + +void mw_fm2(MultiResol &MR_Data, MRNoiseModel & NoiseModel, + fltarray &TabAlpha, CMemWave & MemWObj, Bool DataSNR, Bool DataModel, MultiResol *MR_Model, + float ConvgParam, int MaxIter, Bool PositivConst, Bool Verbose, Bool UseSupport) +/* Apply the multiscale entropy to the multiresolution data set MR_Data + The Parameter Alpha is estimated iteratively (by dichotomy) in order to have + a residu compatible with the noise. + + MR_Data = in-out: multiresolution transform of an image + NoiseModel = in: noise modeling of the image in the wavelet space + TabAlpha = in: TabAlpha(b) gives the alpha regularization parameter for + the band b + MemWObj = in: class for entropy filtering + UseModel = in: true if we have a model + MR_Model = in: pointer to a multiresolution model or NULL (if none) + UseDataSNR = in: true if alpha parameter are modified using the SNR + of the data + ConvgParam = in: Convergence parameter + MaxIter = in: maximum number of iterations + PositivConst = in: apply the positivy constraint + UseSupport= in: true if alpha parameter are modified using the multiresolution + support +*/ +{ + int Nbr_Band = NoiseModel.nbr_band()-1; + int Nbr_Plan = NoiseModel.nbr_scale(); + int Nl = MR_Data.size_ima_nl(); + int Nc = MR_Data.size_ima_nc(); + type_transform Transform = NoiseModel.type_trans(); + MultiResol MR_Sol(Nl,Nc,Nbr_Plan,Transform,"MR_Sol"); + Ifloat Data(Nl,Nc, "data"); + Ifloat Ima(Nl,Nc, "Ima"); + double RegulMin, RegulMax, SigmaNoise, Delta; + int b,Iter=0; + float RegulParam; + float AlphaUser = TabAlpha(0); + + // we need the raw data in order to be able to calculate the residual + // image: resi = data - solution + MR_Data.recons(Data); + + // initialization + MR_Sol.band(Nbr_Band) = MR_Data.band(Nbr_Band); + RegulMin=0.0; // minimum alpha value + RegulMax=DEF_MAX_MEM_ALPHA; // maximum alpha value + + cout << "Sigma Noise= " << NoiseModel.SigmaNoise << endl; + do + { + Iter ++; + // copy the raw wavelet coef. into the solution array + for (b = 0; b < Nbr_Band; b++) MR_Sol.band(b) = MR_Data.band(b); + + // RegulParam = alpha parameter new calculation + RegulParam = (RegulMin + RegulMax)/2; + + // put it in TabAlpha + for (b = 0; b < Nbr_Band; b++) TabAlpha(b) = RegulParam; + + // correct the wavelet coefficients + mw_fm1(MR_Sol, NoiseModel,TabAlpha, MemWObj, DataModel, MR_Model, False, False); + + // reconstruct the filtered image + MR_Sol.recons(Ima); + + // Apply positivuty constraint + if (PositivConst == True) threshold(Ima); + + // SigmaNoise = sigma(residual) + SigmaNoise = sigma(Data-Ima); + + // modify RegulMin,RegulMax by dichotomy + if (SigmaNoise > NoiseModel.SigmaNoise) RegulMax=RegulParam; + else RegulMin = RegulParam; + + // distance between sigma(residual) and sigma(noise) + Delta = ABS(SigmaNoise - NoiseModel.SigmaNoise); + + if (Verbose == True) + cout << Iter << ": Sigma Resi = " << SigmaNoise << " Regul. Param = " << RegulParam << endl; + + // convergence test + } while( (Delta>ConvgParam) && (Iter < MaxIter)); + + // the calculate alpha is + RegulParam = (RegulMin + RegulMax)/2; + for (b = 0; b < Nbr_Band; b++) TabAlpha(b) = RegulParam*AlphaUser; + + if (Verbose == True) + cout << "Optimal Regul. Param = " << RegulParam << endl; + + // final correction of the wavelet coefficients + mw_fm1(MR_Data, NoiseModel,TabAlpha, MemWObj, DataModel, MR_Model,DataSNR,UseSupport); +} + +/***************************************/ + +void mw_fm3(MultiResol &MR_Data, MRNoiseModel & NoiseModel, + fltarray &TabAlpha, CMemWave & MemWObj, Bool DataSNR, Bool DataModel, MultiResol *MR_Model, + float ConvgParam, int MaxIter, Bool Verbose, + Bool RecEachIter, Bool PosConstraint, Bool UseSupport) +/* Apply the multiscale entropy to the multiresolution data set MR_Data + The Parameter Alpha is estimated iteratively at scale separately + (by dichotomy) in order to have a residu compatible with the noise + at all the scales + + MR_Data = in-out: multiresolution transform of an image + NoiseModel = in: noise modeling of the image in the wavelet space + TabAlpha = in: TabAlpha(b) gives the alpha regularization parameter for + the band b + MemWObj = in: class for entropy filtering + UseModel = in: true if we have a model + MR_Model = in: pointer to a multiresolution model or NULL (if none) + UseDataSNR = in: true if alpha parameter are modified using the SNR + of the data + ConvgParam = in: Convergence parameter + MaxIter = in: maximum number of iterations + PositivConst = in: apply the positivy constraint + UseSupport= in: true if alpha parameter are modified using the multiresolution + support +*/ +{ + int Nbr_Band = NoiseModel.nbr_band()-1; + int Nbr_Plan = NoiseModel.nbr_scale(); + int Nl = MR_Data.size_ima_nl(); + int Nc = MR_Data.size_ima_nc(); + type_transform Transform = NoiseModel.type_trans(); + MultiResol MR_Sol(Nl,Nc,Nbr_Plan,Transform,"MR_Sol"); + Ifloat ImaSol(Nl,Nc, "Ima"); + double SigmaNoise; + int b,i,j, Iter=0; + fltarray RegulMin( Nbr_Band ); + fltarray RegulMax( Nbr_Band ); + fltarray TabDelta( Nbr_Band ); + float AlphaUser = TabAlpha(0); + + // initialization + MR_Sol.band(Nbr_Band) = MR_Data.band(Nbr_Band); + + // initialization + for (b = 0; b < Nbr_Band; b++) + { + RegulMin(b) = 0; // minimum alpha value + RegulMax(b) = DEF_MAX_MEM_ALPHA; // maximum alpha value + if (MR_Data.band_to_scale(b) == 0) TabAlpha(b)= 5; + else TabAlpha(b)= 1.; + TabDelta(b) = 0.; + } + if (Verbose == True) + { + if (RecEachIter == True) cout << "Reconstruction at each iteration ... " << endl; + if (UseSupport == True) cout << "Use the multiresolution support ... " << endl; + if (DataSNR == True) cout << "Protect high coeff from the regularization ... " << endl; + cout << " AlphaUser = " << AlphaUser << endl; + } + + do + { + Iter ++; + if (Verbose == True) cout << "Iter " << Iter << endl; + + // copy the raw wavelet coef. into the solution array + for (b = 0; b <= Nbr_Band; b++) MR_Sol.band(b) = MR_Data.band(b); + + // RegulParam = alpha parameter new calculation + for (b = 0; b < Nbr_Band; b++) + TabAlpha(b) = (RegulMin(b) + RegulMax(b))/2; + + // correct the wavelet coefficients + mw_fm1(MR_Sol, NoiseModel,TabAlpha, MemWObj, DataModel, MR_Model, DataSNR, UseSupport); + + // transform the result, apply the positivity contrain + // and reconstruct + if (RecEachIter == True) + { + MR_Sol.recons(ImaSol); + INFO_X(ImaSol, "RES"); + if (PosConstraint == True) threshold(ImaSol); + MR_Sol.transform(ImaSol); + // MR_Sol.band(Nbr_Band) = MR_Data.band(Nbr_Band); + } + + // new estimation of the Alpha paramters + for (b = 0; b < Nbr_Band; b++) + { + int Nlb = MR_Data.size_band_nl(b); + int Ncb = MR_Data.size_band_nc(b); + SigmaNoise = 0.; + long Nc = 0; + for (i = 0; i < Nlb; i++) + for (j = 0; j < Ncb; j++) + { + if ((UseSupport == False ) || + ((UseSupport == True) && (NoiseModel(b,i,j) == False))) + { + double Resi = MR_Data(b,i,j) - MR_Sol(b,i,j); + double SigmaCoef = NoiseModel.sigma(b,i,j); + double ExpectVariance = SigmaCoef*SigmaCoef; + if (ExpectVariance > FLOAT_EPSILON) + SigmaNoise += (Resi*Resi) / ExpectVariance; + else SigmaNoise += 1 + Resi*Resi; + Nc++; + } + } + if (Nc > 0) SigmaNoise = sqrt(SigmaNoise/(double) (Nc)); + else SigmaNoise = 1.; + if (SigmaNoise >= 1) RegulMax (b) = TabAlpha(b); + else RegulMin(b) = TabAlpha(b); + + TabDelta(b) = ABS(SigmaNoise - 1.); + if (Verbose == True) + cout << " band " << b+1 << " Normalized sigma " << SigmaNoise << " Alpha = " << TabAlpha(b) << " Delta = " << TabDelta(b) << endl; + } + // convergence test + } while( (TabDelta.max() > ConvgParam) && (Iter < MaxIter)); + + // the calculate alpha is + for (b = 0; b < Nbr_Band; b++) + TabAlpha(b) = (RegulMin(b) + RegulMax(b))/2*AlphaUser; + + if (Verbose == True) + for (b = 0; b < Nbr_Band; b++) + cout << "band " << b+1 << " Optimal Regul. Param = " << TabAlpha(b) << endl; + + // final correction of the wavelet coefficients + mw_fm1(MR_Data, NoiseModel, TabAlpha, MemWObj, DataModel, MR_Model,DataSNR,UseSupport); + MR_Data.recons(ImaSol); +} + +/***************************************/ + +// same as fm3 but alpha parameter is not estimated by dichotomy, but +// fix step gradient +// void mw_fm3_bis(MultiResol &MR_Data, MRNoiseModel & NoiseModel, +// fltarray &TabAlpha, CMemWave &MemWObj, Bool DataSNR, +// Bool DataModel, MultiResol *MR_Model, +// float ConvgParam, int MaxIter) +// /* Apply the multiscale entropy to the multiresolution data set MR_Data. +// The alpha parameter is different at each scale, and all alpha are +// estimated iteratively in order to have a residu compatible with the noise +// at all the scale. +// +// MR_Data = in-out: multiresolution transform of an image +// NoiseModel = in: noise modeling of the image in the wavelet space +// TabAlpha = in: TabAlpha(b) gives the alpha regularization parameter for +// the band b +// MemWObj = in: class for entropy filtering +// UseModel = in: true if we have a model +// MR_Model = in: pointer to a multiresolution model or NULL (if none) +// UseDataSNR = in: true if alpha parameter are modified using the SNR +// of the data +// ConvgParam = in: Convergence parameter +// MaxIter = in: maximum number of iterations +// +// minimization of alpha by fixed step gradient +// ccl: dichotomy seem better adapted. +// */ +// { +// int Nbr_Band = MR_Data.nbr_band()-1; +// int Nl = MR_Data.size_ima_nl(); +// int Nc = MR_Data.size_ima_nc(); +// type_transform Transform = NoiseModel.type_trans(); +// int Nbr_Plan = MR_Data.nbr_scale(); +// MultiResol MR_Sol(Nl,Nc,Nbr_Plan,Transform,"MR_Sol"); +// MultiResol MR_Resi(Nl,Nc,Nbr_Plan,Transform,"MR_Resi"); +// +// Ifloat Data(Nl,Nc, "data"); +// Ifloat Ima(Nl,Nc, "data"); +// double SigmaNoise; +// int b,i,j,Iter=0; +// fltarray TabDelta( Nbr_Band ); +// fltarray TabStep( Nbr_Band ); +// +// // we need the raw data in order to be able to calculate the residual +// // image: resi = data - solution +// MR_Data.recons(Data); +// +// // initialization +// for (b = 0; b < Nbr_Band; b++) +// { +// if (MR_Data.band_to_scale(b) == 0) TabAlpha(b)= 5; +// else TabAlpha(b)= 1.; +// TabDelta(b) = 0.; +// } +// +// +// +// do +// { +// Iter ++; +// cout << "Iter " << Iter << endl; +// +// // copy the raw wavelet coef. into the solution array +// for (b = 0; b <= Nbr_Band; b++) MR_Sol.band(b) = MR_Data.band(b); +// +// // correct the wavelet coefficients +// mw_fm1(MR_Sol, NoiseModel,TabAlpha, MemWObj, DataModel, MR_Model,DataSNR); +// +// // calculate the residual image +// //MR_Sol.recons(Ima); +// //Ima = Data - Ima; +// +// // wavelet transform of the residual +// //MR_Resi.transform(Ima); +// +// // new estimation of the Alpha paramters +// for (b = 0; b < Nbr_Band; b++) +// { +// int Nlb = MR_Data.size_band_nl(b); +// int Ncb = MR_Data.size_band_nc(b); +// SigmaNoise = 0.; +// +// // standard deviation of the residual in the band b +// // and Alpha paramter new estimation +// float SumNum = 0.; +// // float SumDen = 0.; +// for (i = 0; i < Nlb; i++) +// for (j = 0; j < Ncb; j++) +// { +// float dhs,dhn,Coef = MR_Sol(b,i,j); +// float Resi= MR_Data(b,i,j) - Coef; +// //float Resi= MR_Resi(b,i,j); +// float SigmaCoef = NoiseModel.sigma(b,i,j); +// float ExpectVariance = SigmaCoef*SigmaCoef; +// SigmaNoise += (Resi*Resi) / ExpectVariance; +// +// dhn = MemWObj.grad_hn(Coef); +// dhs = MemWObj.grad_hs(Resi); +// SumNum += dhn/ExpectVariance * Resi; +// //SumNum += dhn/ExpectVariance*(-Resi-dhs); +// //SumDen += dhn*dhn/ExpectVariance; +// } +// SigmaNoise = sqrt(SigmaNoise/(Nlb*Ncb)); +// if (SigmaNoise >= 1) TabAlpha(b) -= 0.2*SumNum/(Nlb*Ncb); +// else TabAlpha(b) += 0.2*SumNum/(Nlb*Ncb); +// if ( TabAlpha(b) <= 0.) TabAlpha(b) = 0.; +// //if (SumDen > FLOAT_EPSILON) SumNum /= SumDen; +// //else TabAlpha(b) = 1.; +// //if (SumNum <= 0.) TabAlpha(b) /= 2.; +// //else TabAlpha(b) = SumNum; +// +// // distance between sigma(residual) and sigma(noise) +// TabDelta(b) = ABS(SigmaNoise - 1.); +// +// cout << " band " << b+1 << " Normalized sigma " << SigmaNoise << " Alpha = " << TabAlpha(b) << " Delta = " << TabDelta(b) << endl; +// } +// +// // convergence test +// } while( (TabDelta.max() > ConvgParam) && (Iter < MaxIter)); +// +// +// +// for (b = 0; b < Nbr_Band; b++) +// cout << "band " << b+1 << " Optimal Regul. Param = " << TabAlpha(b) << endl; +// +// // final correction of the wavelet coefficients +// mw_fm1(MR_Data, NoiseModel,TabAlpha, MemWObj, DataModel, MR_Model,DataSNR); +// } + +/****************************************************************************/ diff --git a/src/libsparse2d/MW_Filter.h b/src/libsparse2d/MW_Filter.h new file mode 100644 index 0000000..e36c123 --- /dev/null +++ b/src/libsparse2d/MW_Filter.h @@ -0,0 +1,153 @@ +/****************************************************************************** +** Copyright (C) 1995 by CEA +******************************************************************************* +** +** UNIT +** +** Version: 3.1 +** +** Author: Jean-Luc Starck +** +** Date: 96/05/02 +** +** File: MR_Filter.h +** +******************************************************************************* +** +** DESCRIPTION +** ----------- +** +** PARAMETRES +** ---------- +** +** +** RESULTS +** ------- +** +** +******************************************************************************/ + + +#ifndef __MWFILTER__ +#define __MWFILTER__ + +const int DEF_MEM_ALPHA_CST = 1; +const int DEF_MEM_ALPHA_OPT = 2; +const int DEF_MEM_ALPHA_OPT_SCALE = 3; +const float DEF_MAX_MEM_ALPHA = 100.0; + +const int DEF_MEM_FILT_MAX_ITER = 20; +const float DEF_MEM_FILT_CONGER = 0.01; + +/****************************************************************************/ + +void mw_filter(MultiResol &MR_Data, MRNoiseModel & NoiseModel, + int TypeFilter, fltarray &TabAlpha, Bool DataSNR=False, + Bool DataModel=False, MultiResol *MR_Model=NULL, + float ConvgParam = DEF_MEM_FILT_CONGER, int MaxIter=DEF_MEM_FILT_MAX_ITER, + Bool PositivConst=False, Bool Verbose=False, Bool UseSupport=False); +/* Apply the multiscale entropy to the multiresolution data set MR_Data + MR_Data = in-out: multiresolution transform of an image + NoiseModel = in: noise modeling of the image in the wavelet space + TypeFilter = in: type of filtering + = DEF_MEM_ALPHA_CST ==> alpha is fixed + = DEF_MEM_ALPHA_OPT ==> alpha is globally optimized + = DEF_MEM_ALPHA_OPT_SCALE ==> alpha is optimized at each scale + + TabAlpha = in-out: TabAlpha(b) gives the alpha regularization parameter for + the band b + DataModel = in: true if we have a model + MR_Model = in: pointer to a multiresolution model or NULL (if none) + DataSNR = in: true if alpha parameter are modified using the SNR + of the data + PositivConst = in: apply the positivy constraint + UseSupport= in: true if alpha parameter are modified using the multiresolution + support + Warning: works well, when each wavelet can be modelized by a Gaussian +*/ + +/****************************************************************************/ + +void mw_fm1(MultiResol & MR_Data, MRNoiseModel & NoiseModel, + fltarray &TabAlpha, CMemWave & MemWObj, + Bool UseModel, MultiResol *MR_Model, Bool UseDataSNR, Bool UseSupport=False); +/* Apply the multiscale entropy to the multiresolution data set MR_Data + MR_Data = in-out: multiresolution transform of an image + NoiseModel = in: noise modeling of the image in the wavelet space + TabAlpha = in: TabAlpha(b) gives the alpha regularization parameter for + the band b + MemWObj = in: class for entropy filtering + UseModel = in: true if we have a model + MR_Model = in: pointer to a multiresolution model or NULL (if none) + UseDataSNR = in: true if alpha parameter are modified using the SNR + of the data + UseSupport= in: true if alpha parameter are modified using the multiresolution + support + +*/ + + +/****************************************************************************/ + +void mw_fm2(MultiResol &MR_Data, MRNoiseModel & NoiseModel, + fltarray &TabAlpha, CMemWave & MemWObj, Bool DataSNR, Bool DataModel, MultiResol *MR_Model, + float ConvgParam, int MaxIter, Bool PositivConst=False, Bool Verbose=False, Bool UseSupport=False); +/* Apply the multiscale entropy to the multiresolution data set MR_Data + The Parameter Alpha is estimated iteratively (by dichotomy) in order to have + a residu compatible with the noise. + + MR_Data = in-out: multiresolution transform of an image + NoiseModel = in: noise modeling of the image in the wavelet space + TabAlpha = in: TabAlpha(b) gives the alpha regularization parameter for + the band b + MemWObj = in: class for entropy filtering + UseModel = in: true if we have a model + MR_Model = in: pointer to a multiresolution model or NULL (if none) + UseDataSNR = in: true if alpha parameter are modified using the SNR + of the data + ConvgParam = in: Convergence parameter + MaxIter = in: maximum number of iterations + UseSupport= in: true if alpha parameter are modified using the multiresolution + support +*/ + +/****************************************************************************/ + + +void mw_fm3(MultiResol &MR_Data, MRNoiseModel & NoiseModel, + fltarray &TabAlpha, CMemWave & MemWObj, Bool DataSNR, Bool DataModel, MultiResol *MR_Model, + float ConvgParam, int MaxIter, Bool Verbose=False, + Bool RecEachIter=False, Bool PosConstraint=False, Bool UseSupport=False); +/* Apply the multiscale entropy to the multiresolution data set MR_Data + The Parameter Alpha is estimated iteratively at scale separately + (by dichotomy) in order to have a residu compatible with the noise + at all the scales + + MR_Data = in-out: multiresolution transform of an image + NoiseModel = in: noise modeling of the image in the wavelet space + TabAlpha = in: TabAlpha(b) gives the alpha regularization parameter for + the band b + MemWObj = in: class for entropy filtering + UseModel = in: true if we have a model + MR_Model = in: pointer to a multiresolution model or NULL (if none) + UseDataSNR = in: true if alpha parameter are modified using the SNR + of the data + ConvgParam = in: Convergence parameter + MaxIter = in: maximum number of iterations + RecEachIter = in: if true, a reconstruction, followed by a transform + is performed at each iteration. This allows to resolve + the problem of non-orthogonality of the coefficient. + Set this option for orthogonal transform has no sense. + PosConstraint: in: if set and if RecEachIter is set, the positivity contraint + is applied. + UseSupport= in: true if alpha parameter are modified using the multiresolution + support +*/ + +/****************************************************************************/ + +#endif + + + + diff --git a/src/libtools/CMem.cc b/src/libtools/CMem.cc new file mode 100644 index 0000000..7f2285f --- /dev/null +++ b/src/libtools/CMem.cc @@ -0,0 +1,231 @@ +/****************************************************************************** +** Copyright (C) 1998 by CEA +******************************************************************************* +** +** UNIT +** +** Version: 1.0 +** +** Author: Jean-Luc Starck +** +** Date: 98/01/12 +** +** File: CMem.cc +** +******************************************************************************* +** +** DESCRIPTION class MEM for signal restoration +** ----------- +** +******************************************************************************/ + +#include "GlobalInc.h" +#include "CErf.h" +#include "CMem.h" + +/****************************************************************************/ + +// calculate the derivative of hn in case of Gaussian noise +// hn = x/Sigma^2*erfc(x/(sqrt(2)sigma)) + +// sqrt(2/PI)/sigma[1-exp(-x^2/(2sigma^2))] +double CMemWave::grad_hn_sig1(double Val, CERF & ErTab) +{ + int Sign = (Val >= 0) ? 1: -1; + double Coef = Val*Sign; + Coef = Coef*ErTab.get_erfc(Coef/C2) + C1 * (1-exp(-(Coef*Coef)/2.)); + return Coef*Sign; +} + +/****************************************************************************/ + +// calculate the derivative of hs in case of Gaussian noise +// hs = -x/Sigma^2*erf(x/(sqrt(2)sigma)) + +// sqrt(2/PI)/sigma[1-exp(-x^2/(2sigma^2))] +double CMemWave::grad_hs_sig1(double Val, CERF & ErTab) +{ + int Sign = (Val >= 0) ? 1: -1; + double Coef = Val*Sign; + Coef = -Coef*ErTab.get_erf(Coef/C2) + C1 * (1-exp(-(Coef*Coef)/2.)); + return Coef*Sign; +} + +/***************************************/ + + // initialization +CMemWave::CMemWave() +{ + double Val=0.; + int i=0; + CERF ErTab; + + MaxIter= DEF_MEM_MAX_ITER; + ConvgParam= DEF_MEM_CVG; + C1 = sqrt(2./PI); + C2 = sqrt(2.); + Step = DEF_MEM_STEP; + MaxVal = DEF_MEM_MAX; + + Np = (int) (MaxVal / Step + 1.5); + TabHsGauss = new double [Np]; + TabHnGauss = new double [Np]; + + while (Val < MaxVal) + { + TabHsGauss[i] = grad_hs_sig1(Val, ErTab); + TabHnGauss[i++] = grad_hn_sig1(Val, ErTab); + Val += Step; + } +} + +/***************************************/ + +float CMemWave::grad_hs(float Val) +{ + int Ind; + double ValRet = 0.; + + if (Val >= MaxVal) ValRet = TabHsGauss[Np-1] - Val; + else if (Val > 0) + { + Ind = (int) (Val / Step + 0.5); + ValRet = TabHsGauss[Ind]; + } + return ValRet; +} + +/***************************************/ + +float CMemWave::grad_hn(float Val) +{ + int Ind; + double ValRet = 0.; + + if (Val >= MaxVal) ValRet = TabHnGauss[Np-1]; + else if (Val > 0) + { + Ind = (int) (Val / Step + 0.5); + ValRet = TabHnGauss[Ind]; + } + return ValRet; +} + +/***************************************/ + +// calculate by dichotomy method the solution wich minimize: +// J = hs(y-x) + Alpha hn(x) +// It is obtained when: grad_hs(y-x) = Alpha grad_hn(x) +float CMemWave::filter (float CoefDat, float Alpha, float Sig, float Model) +{ + float ValRet=0; + float Sigma=Sig; + if (Sigma < FLOAT_EPSILON) Sigma = FLOAT_EPSILON ; + + if (ABS(Alpha) < FLOAT_EPSILON) ValRet = CoefDat; + else if (Alpha < 0) ValRet = Model; + else + { + double ValDat = ABS(CoefDat-Model) / Sigma; + double CoefMax = ValDat+3; + double Coef,Diff,CoefMin = 0; + int Iter=0; + do + { + Iter++; + Coef = (CoefMin+CoefMax)/2.; + Diff = grad_hs(ValDat-Coef) + Alpha * grad_hn(Coef); + if (Diff > 0.) CoefMax = Coef; + else CoefMin = Coef; + } while (( CoefMax-CoefMin > ConvgParam) && (Iter < MaxIter)); + Coef = (CoefDat >= Model) ? Coef : -Coef; + ValRet = Coef*Sigma+Model; + } + return ValRet; + } + +/***************************************/ + + +// calculate by dichotomy method the solution wich minimize: +// J = hs(y-x) + Alpha hn(x) +// It is obtained when: grad_hs(y-x) = Alpha grad_hn(x) +float CMemWave::filter (float CoefDat, float Alpha, float SigS, float SigN, float Model) +{ + float ValRet=0; + float Sigma=SigS; + float SigmaN=SigN; + + if (Sigma < FLOAT_EPSILON) Sigma = FLOAT_EPSILON ; + if (SigmaN < FLOAT_EPSILON) SigmaN = FLOAT_EPSILON ; + + if (ABS(Alpha) < FLOAT_EPSILON) ValRet = CoefDat; + else if (Alpha < 0) ValRet = Model; + else + { + double ValDat = ABS(CoefDat-Model); + double CoefMax = ValDat+3*MAX(Sigma,SigmaN); + double Coef,Diff,CoefMin = 0; + int Iter=0; + do + { + Iter++; + Coef = (CoefMin+CoefMax)/2.; + Diff = Sigma*grad_hs( (ValDat-Coef)/Sigma) + Alpha * SigmaN * grad_hn(Coef/SigmaN); + if (Diff > 0.) CoefMax = Coef; + else CoefMin = Coef; + } while (( (CoefMax-CoefMin)/Sigma > ConvgParam) && (Iter < MaxIter)); + Coef = (CoefDat >= Model) ? Coef : -Coef; + ValRet = Coef+Model; + } + return ValRet; + } + /***************************************/ + +void CMemWave::test_tab(fltarray &TabTest, float NSig, Bool UseDataSNR) +{ + int k,U=1000; + fltarray Tab(U,8); + float AlphaP, Rg=1., RegulParam; + float TabReg[7]; + TabReg[0] = 0.; + TabReg[1] = 0.1; + TabReg[2] = 0.5; + TabReg[3] = 1.; + TabReg[4] = 2.; + TabReg[5] = 5.; + TabReg[6] = 10.; + + for (int i=0; i< U; i++) + { + Tab(i,0) = 0.01*i; + if (UseDataSNR == True) + { + AlphaP = ABS( Tab(i,0) / NSig); + if (AlphaP > 1) AlphaP = 1; + if (AlphaP < FLOAT_EPSILON) Rg = -1; + else Rg = (1-AlphaP)/AlphaP; + } + for (k=1; k < 7; k++) + { + RegulParam = TabReg[k]*Rg; + if (ABS(RegulParam) < FLOAT_EPSILON) Tab(i,k)= Tab(i,0); + else if (RegulParam < 0) Tab(i,k) = 0; + else Tab(i, k) = filter( Tab(i,0), RegulParam); + } + + // cout << "Val = " << i*0.5 << " Dicho = " << filter(0.01*i, RegulParam) << endl; + } + TabTest = Tab; + } + +/***************************************/ + +CMemWave::~CMemWave() +{ + if (TabHnGauss != NULL) delete [] TabHnGauss; + if ( TabHsGauss != NULL) delete [] TabHsGauss; + Np = 0; +} + +/***************************************/ + + diff --git a/src/libtools/CMem.h b/src/libtools/CMem.h new file mode 100644 index 0000000..49c8b13 --- /dev/null +++ b/src/libtools/CMem.h @@ -0,0 +1,87 @@ +/****************************************************************************** +** Copyright (C) 1998 by CEA +******************************************************************************* +** +** UNIT +** +** Version: 1.0 +** +** Author: Jean-Luc Starck +** +** Date: 98/01/12 +** +** File: Cerf.h +** +******************************************************************************* +** +** DESCRIPTION: class MEM for signal restoration +** ----------- +** +******************************************************************************/ + +#ifndef _CMEM_H_ +#define _CMEM_H_ + +const float DEF_MEM_STEP = 0.01; +const float DEF_MEM_MAX = 5.; +const int DEF_MEM_MAX_ITER = 100; +const float DEF_MEM_CVG = 0.001; + +class CMemWave { + double C1; // internal constant = sqrt(2/PI) + double C2; // internal constant = sqrt(2) + + + double *TabHsGauss; // table for the derivative of Hs + double *TabHnGauss; // table for the derivative of Hn + int Np; // size of the tables + + float Step; // Step used for the probability integration + float MaxVal; // Maximum value pre-calculated + + double grad_hn_sig1(double Val, CERF & ErTab); + // calculate the derivative of hn in case of Gaussian noise (sigma=1) + // hn = x/Sigma^2*erfc(x/(sqrt(2)sigma)) + + // sqrt(2/PI)/sigma[1-exp(-x^2/(2sigma^2))] + + double grad_hs_sig1(double Val, CERF & ErTab); + // calculate the derivative of hs in case of Gaussian noise (sigma=1) + // hs = -x/Sigma^2*erf(x/(sqrt(2)sigma)) + + // sqrt(2/PI)/sigma[1-exp(-x^2/(2sigma^2))] + + public: + int MaxIter; // maximum number iteration for the dichotomy + float ConvgParam; // convergence parameter for the dichotomy + + CMemWave(); // initialization: set the table + + float grad_hs(float Val); // calculate the derivative of hs at Val + + float grad_hn(float Val); // calculate the derivative of hn at Val + + float filter (float CoefDat, float Alpha, float Sigma=1., float Model=0.); + // calculate by the dichotomy method the solution wich minimize: + // J = hs(y-x) + Alpha hn(x) + // It is obtained when: grad_hs(y-x) = Alpha grad_hn(x) + + float filter (float CoefDat, float Alpha, float SigmaS, float SigmaN, float Model); + // calculate by the dichotomy method the solution wich minimize: + // J = hs(y-x) + Alpha hn(x) + // It is obtained when: grad_hs(y-x) = Alpha grad_hn(x) + + void test_tab(fltarray &TabTest, float Nsig=3., Bool UseDataSNR=False); + // set a table TabTest(8,1000) + // TabTest(*,i) represents the filtering of a signal S, with + // S(k) = 0.01 * k : S is in [0,10[, with a step = 0.01 + // i = 0 ==> alpha = 0 + // i = 1 ==> alpha = 0.1 + // i = 2 ==> alpha = 0.5 + // i = 3 ==> alpha = 1 + // i = 4 ==> alpha = 2 + // i = 5 ==> alpha = 5 + // i = 6 ==> alpha = 10.; + + ~CMemWave(); +}; + +#endif diff --git a/src/libtools/IM_Block2D.cc b/src/libtools/IM_Block2D.cc new file mode 100644 index 0000000..511a8a0 --- /dev/null +++ b/src/libtools/IM_Block2D.cc @@ -0,0 +1,187 @@ +/****************************************************************************** +** Copyright (C) 2003 by CEA +******************************************************************************* +** +** UNIT +** +** Version: 1.0 +** +** Author: Jean-Luc Starck +** +** Date: 28/03/2003 +** +** File: IM_Block2D.cc +** +** Modification history: +** +******************************************************************************* +** +** DESCRIPTION 2D block manadgement +** ----------- +** +******************************************************************************/ + +#include "IM_Block2D.h" + +/***********************************************************************/ + +void Block2D::alloc(int Nlc, int Ncc, int ParamBlockSize, float Overlapping) +{ + Overlap = Overlapping; + if(Overlap>0) BlockOverlap=True; + alloc(Nlc, Ncc, ParamBlockSize); +} + +/***********************************************************************/ + +void Block2D::alloc(int Nlc, int Ncc, int ParamBlockSize) +{ + int MaxXY; + Nl=Nlc; + Nc=Ncc; + BlockSize = ParamBlockSize; + + + if ((BlockOverlap == True) && (Overlap==0)) Overlap=0.5; + MaxXY = MAX(Nl,Nc); + + if ((BlockSize==0) || (BlockSize==MaxXY)) + { + Nlb = Ncb = 1; + BlockSize = MaxXY; + BlockOverlap=False; + Overlap = 0; + OverlapSize = 0; + } + else + { + BlockImaSize = ceil(BlockSize*(1-Overlap)); + OverlapSize = BlockSize-BlockImaSize; + + // cout << "BlockSize = " << BlockSize << endl; + + // cout << "BlockImaSize = " << BlockImaSize << endl; + int ix=Nl-((Nl-OverlapSize)/BlockImaSize)*BlockImaSize; + Nlb = ((Nl-OverlapSize)/BlockImaSize) + (ix!=OverlapSize); + int iy=Nc-((Nc-OverlapSize)/BlockImaSize)*BlockImaSize; + Ncb = ((Nc-OverlapSize)/BlockImaSize) + (iy!=OverlapSize); + + // Nlb = (Nl % BlockImaSize == 0) ? Nl / BlockImaSize: Nl / BlockImaSize + 1; + // Ncb = (Nc % BlockImaSize == 0) ? Nc / BlockImaSize: Nc / BlockImaSize + 1; + } + + NbrBlock = Nlb*Ncb; + // if (BlockOverlap == False) BlockImaSize = 0; + if (Verbose == True) + { + printf( "\nima size = (%2d,%2d), BlockSize = %2d, BlockImaSize = %2d, OverlapSize = %2d\n", Nl,Nc,BlockSize, BlockImaSize, OverlapSize); + printf( "Nbr Blocks = (%2d,%2d)\n", Nlb,Ncb); + printf( "NbrBlocks*BlockImaSize + OverlapSize = (%2d,%2d)\n", Nlb*BlockImaSize+OverlapSize, Ncb*BlockImaSize+OverlapSize); + } +} + +/***********************************************************************/ + +double Block2D::get_weight(int PosPix, int CurrentBlock, int MaxBlockNumber) +{ + double ValReturn; + double ValL=1.,ValR=1.; + + if((PosPix < OverlapSize) && (CurrentBlock > 0)) + ValL = (double) pow(sin((double) (PosPix+1) + / (double) OverlapSize * PI/2.), 2.); + + if((PosPix > BlockImaSize-1) && (CurrentBlock < MaxBlockNumber-1)) + ValR = (double) pow(cos((double) (PosPix-(BlockImaSize-1)) + / (double) OverlapSize * PI/2.), 2.); + + ValReturn = ValL*ValR; + + return ValReturn; +} + +/*************************************************************************/ + + +void Block2D::get_block_ima(int Bi, int Bj, Ifloat &Ima, Ifloat &ImaBlock, Bool Weight) +// Extract the block (Bi,Bj) from Ima and put it in ImaBlock +{ + int k,l; + int Depi = BlockImaSize*Bi; + int Depj = BlockImaSize*Bj; + + if ((WeightFirst == False) || (Weight == False)) + { + // cout << "get_block_ima NO Weight " << endl; + for (k = 0; k < BlockSize; k++) + for (l = 0; l < BlockSize; l++) + { + int Indk = test_index_mirror(Depi+k,Nl); + int Indl = test_index_mirror(Depj+l,Nc); + ImaBlock(k,l) = Ima(Indk,Indl); + } + } + else + { + // cout << "get_block_ima Weight " << endl; + for (k = 0; k < BlockSize; k++) + for (l = 0; l < BlockSize; l++) + { + int Indk = test_index_mirror(Depi+k,Nl); + int Indl = test_index_mirror(Depj+l,Nc); + ImaBlock(k,l) = Ima(Indk,Indl) * get_weight(k, Bi, Nlb) * get_weight(l, Bj, Ncb); + } + } + // cout << "BLOCK: " << Bi <<","<= 0) && (Depi+k < Nl) + && (Depj+l >= 0) && (Depj+l < Nc)) + Ima(Depi+k,Depj+l) = ImaBlock(k,l); + // cout << "BLOCK: " << Bi <<","<= 0) && (Depi+k < Nl) + && (Depj+l >= 0) && (Depj+l < Nc)) + Ima(Depi+k,Depj+l) += ImaBlock(k,l)* + get_weight(k, Bi, Nlb) * get_weight(l, Bj, Ncb); + } + else + { + // cout << "add_block_ima NO Weight " << endl; + for (k = 0; k < BlockSize; k++) + for (l = 0; l < BlockSize; l++) + if ((Depi+k >= 0) && (Depi+k < Nl) + && (Depj+l >= 0) && (Depj+l < Nc)) + Ima(Depi+k,Depj+l) += ImaBlock(k,l); + } +} + +/****************************************************************************/ diff --git a/src/libtools/IM_Block2D.h b/src/libtools/IM_Block2D.h new file mode 100644 index 0000000..fe7655d --- /dev/null +++ b/src/libtools/IM_Block2D.h @@ -0,0 +1,69 @@ +/****************************************************************************** +** Copyright (C) 2003 by CEA +******************************************************************************* +** +** UNIT +** +** Version: 1.0 +** +** Author: Jean-Luc Starck +** +** Date: 28/03/2003 +** +** File: IM_Block2D.h +** +** Modification history : +** +******************************************************************************* +** +** DESCRIPTION : CLASS Definition for 2D block manadgement +** ----------- +******************************************************************************/ + +#ifndef _BLOCK2D_H_ +#define _BLOCK2D_H_ + +#include "IM_Obj.h" + +/***********************************************************************/ + +class Block2D { + private: + int Nc, Nl; + int BlockImaSize; // Ima block size without taking account + // the overlapping aera + int BlockSize; // Ima block size taking account + // the overlapping aera + int Ncb,Nlb; // Number of blocks in the x,y directions + int NbrBlock; + void reset_param() {Nc=Nl=0;BlockImaSize=BlockSize=NbrBlock; + Ncb=Nlb=0;BlockOverlap=False; + Verbose=False;WeightFirst=False; Overlap=0.5;OverlapSize=0;} + + public: + Bool Verbose; + Bool BlockOverlap; + Bool WeightFirst; + float Overlap; + int OverlapSize; + inline int nl() { return Nl;} + inline int nc() { return Nc;} + inline int nbr_block_nl() { return Nlb;} + inline int nbr_block_nc() { return Ncb;} + inline int nbr_block() { return NbrBlock;} + inline int block_size() { return BlockSize;} + inline int l_nl() {return BlockSize*Nlb;} + inline int l_nc() {return BlockSize*Ncb;} + Block2D() {reset_param();} + void alloc(int Nlc, int Ncc, int ParamBlockSize); + void alloc(int Nlc, int Ncc, int ParamBlockSize, float Overlaping); + void get_block_ima(int Bi, int Bj, Ifloat &Ima, Ifloat &ImaBlock, Bool Weight=False); + void put_block_ima(int Bi, int Bj, Ifloat &Ima, Ifloat &ImaBlock); + double get_weight(int PosPix, int CurrentBlock, int MaxBlockNumber); + void add_block_ima(int Bi, int Bj, Ifloat &Ima, Ifloat &ImaBlock, Bool Weight=False); + ~Block2D() {reset_param();} +}; + +/***********************************************************************/ + +#endif diff --git a/src/libtools/IM_DCT.cc b/src/libtools/IM_DCT.cc new file mode 100644 index 0000000..38f77ee --- /dev/null +++ b/src/libtools/IM_DCT.cc @@ -0,0 +1,592 @@ +/****************************************************************************** +** Copyright (C) 2000 by CEA +******************************************************************************* +** +** UNIT +** +** Version: 1.0 +** +** Author: Jean-Luc Starck +** +** Date: 10/02/00 +** +** File: IM_DCT.cc +** +** +** +******************************************************************************/ + + +#include"IM_DCT.h" + +extern void ddct2d (int, int, int, double **, double *, int *, double *); +extern void ddst2d (int, int, int, double **, double *, int *, double *); +#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr + +void four1(float data[], unsigned long nn, int isign) +{ + unsigned long n,mmax,m,j,istep,i; + double wtemp,wr,wpr,wpi,wi,theta; + float tempr,tempi; + + n=nn << 1; + j=1; + for (i=1;i i) { + SWAP(data[j],data[i]); + SWAP(data[j+1],data[i+1]); + } + m=n >> 1; + while (m >= 2 && j > m) { + j -= m; + m >>= 1; + } + j += m; + } + mmax=2; + while (n > mmax) { + istep=mmax << 1; + theta=isign*(6.28318530717959/mmax); + wtemp=sin(0.5*theta); + wpr = -2.0*wtemp*wtemp; + wpi=sin(theta); + wr=1.0; + wi=0.0; + for (m=1;m>1); + if (isign == 1) { + c2 = -0.5; + four1(data,n>>1,1); + } else { + c2=0.5; + theta = -theta; + } + wtemp=sin(0.5*theta); + wpr = -2.0*wtemp*wtemp; + wpi=sin(theta); + wr=1.0+wpr; + wi=wpi; + np3=n+3; + for (i=2;i<=(n>>2);i++) { + i4=1+(i3=np3-(i2=1+(i1=i+i-1))); + h1r=c1*(data[i1]+data[i3]); + h1i=c1*(data[i2]-data[i4]); + h2r = -c2*(data[i2]+data[i4]); + h2i=c2*(data[i1]-data[i3]); + data[i1]=h1r+wr*h2r-wi*h2i; + data[i2]=h1i+wr*h2i+wi*h2r; + data[i3]=h1r-wr*h2r+wi*h2i; + data[i4] = -h1i+wr*h2i+wi*h2r; + wr=(wtemp=wr)*wpr-wi*wpi+wr; + wi=wi*wpr+wtemp*wpi+wi; + } + if (isign == 1) { + data[1] = (h1r=data[1])+data[2]; + data[2] = h1r-data[2]; + } else { + data[1]=c1*((h1r=data[1])+data[2]); + data[2]=c1*(h1r-data[2]); + four1(data,n>>1,-1); + } +} + +void cosft2(float y[], int n, int isign) +{ + void realft(float data[], unsigned long n, int isign); + int i; + float sum,sum1,y1,y2,ytemp; + double theta,wi=0.0,wi1,wpi,wpr,wr=1.0,wr1,wtemp; + + theta=0.5*PI/n; + wr1=cos(theta); + wi1=sin(theta); + wpr = -2.0*wi1*wi1; + wpi=sin(2.0*theta); + if (isign == 1) { + for (i=1;i<=n/2;i++) { + y1=0.5*(y[i]+y[n-i+1]); + y2=wi1*(y[i]-y[n-i+1]); + y[i]=y1+y2; + y[n-i+1]=y1-y2; + wr1=(wtemp=wr1)*wpr-wi1*wpi+wr1; + wi1=wi1*wpr+wtemp*wpi+wi1; + } + realft(y,n,1); + for (i=3;i<=n;i+=2) { + wr=(wtemp=wr)*wpr-wi*wpi+wr; + wi=wi*wpr+wtemp*wpi+wi; + y1=y[i]*wr-y[i+1]*wi; + y2=y[i+1]*wr+y[i]*wi; + y[i]=y1; + y[i+1]=y2; + } + sum=0.5*y[2]; + for (i=n;i>=2;i-=2) { + sum1=sum; + sum += y[i]; + y[i]=sum1; + } + } else if (isign == -1) { + ytemp=y[n]; + for (i=n;i>=4;i-=2) y[i]=y[i-2]-y[i]; + y[2]=2.0*ytemp; + for (i=3;i<=n;i+=2) { + wr=(wtemp=wr)*wpr-wi*wpi+wr; + wi=wi*wpr+wtemp*wpi+wi; + y1=y[i]*wr+y[i+1]*wi; + y2=y[i+1]*wr-y[i]*wi; + y[i]=y1; + y[i+1]=y2; + } + realft(y,n,-1); + for (i=1;i<=n/2;i++) { + y1=y[i]+y[n-i+1]; + y2=(0.5/wi1)*(y[i]-y[n-i+1]); + y[i]=0.5*(y1+y2); + y[n-i+1]=0.5*(y1-y2); + wr1=(wtemp=wr1)*wpr-wi1*wpi+wr1; + wi1=wi1*wpr+wtemp*wpi+wi1; + } + } +} + +/******************************************************************************/ + +void im_dct(Ifloat &Ima, Ifloat &Trans, Bool Reverse) +{ + int Nl = Ima.nl(); + int Nc = Ima.nc(); + int i,j; + float Norm = 2./Nc; + int isign = (Reverse == False) ? 1: -1; + float *Ptr = Trans.buffer() - 1; + fltarray Tab; + Trans = Ima; + + if (Reverse == False) + { + for (i = 0; i < Nl; i++) + { + cosft2(Ptr, Nc, isign); + Ptr += Nc; + } + Tab.alloc(Nl); + for (j=0; j < Nc; j++) + { + for (i = 0; i < Nl; i++) Tab(i) = Trans(i,j); + Ptr = Tab.buffer() - 1; + cosft2(Ptr, Nl, isign); + for (i = 0; i < Nl; i++) Trans(i,j) = Tab(i); + } + } + else + { + Tab.alloc(Nl); + for (j=0; j < Nc; j++) + { + for (i = 0; i < Nl; i++) Tab(i) = Ima(i,j); + Ptr = Tab.buffer() - 1; + cosft2(Ptr, Nl, isign); + for (i = 0; i < Nl; i++) Trans(i,j) = Tab(i); + } + Ptr = Trans.buffer() - 1; + for (i = 0; i < Nl; i++) + { + cosft2(Ptr, Nc, isign); + Ptr += Nc; + } + } + for (i = 0; i < Nl; i++) + for (j=0; j < Nc; j++) Trans(i,j) *= Norm; +} + +/******************************************************************************/ + +// +// void im_dctold(Ifloat &Ima, Ifloat &Trans, Bool Reverse) +// { +// int Nl = Ima.nl(); +// int Nc = Ima.nc(); +// int Dir = (Reverse == False) ? 1: -1; +// double **a, *t, *w, *dd; +// int i,j; +// // double Norm = 4. / (double)(Nl*Nc); +// double Norm = 1.; +// +// a = new double * [Nl]; +// dd = new double [Nl*Nc]; +// t = new double [2 * Nl]; +// int n = MAX(Nl, Nc/2); +// int *ip = new int [2 + (int) sqrt(n + 0.5)]; +// n = MAX(Nl * 5 / 4, Nc * 5 / 4) + Nc / 4; +// w = new double [n]; +// ip[0] = 0; +// +// for (i=0; i < Nl; i++) a[i] = &dd[i*Nc]; +// for (i=0; i < Nl; i++) +// for (j=0; j < Nc; j++) a[i][j] = (double) Ima(i,j); +// +// ddct2d (Nl, Nc, Dir, a, t, ip, w); +// +// if (Reverse == True) +// { +// for (i=0; i < Nl; i++) +// for (j=0; j < Nc; j++) Trans(i,j) = (float) (a[i][j]*Norm); +// } +// else +// for (i=0; i < Nl; i++) +// for (j=0; j < Nc; j++) Trans(i,j) = (float) a[i][j]; +// +// +// delete [] dd; +// delete [] a; +// delete [] w; +// delete [] t; +// delete [] ip; +// } + +/******************************************************************************/ + +void im_dst(Ifloat &Ima, Ifloat &Trans, Bool Reverse) +{ + int Nl = Ima.nl(); + int Nc = Ima.nc(); + int Dir = (Reverse == False) ? 1: -1; + double **a, *t, *w, *dd; + int i,j; + double Norm = 4. / (double)(Nl*Nc); + + a = new double * [Nl]; + dd = new double [Nl*Nc]; + t = new double [2 * Nl]; + int n = MAX(Nl, Nc/2); + int *ip = new int [2 + (int) sqrt(n + 0.5)]; + n = MAX(Nl * 5 / 4, Nc * 5 / 4) + Nc / 4; + w = new double [n]; + ip[0] = 0; + + for (i=0; i < Nl; i++) a[i] = &dd[i*Nc]; + for (i=0; i < Nl; i++) + for (j=0; j < Nc; j++) a[i][j] = (double) Ima(i,j); + ddst2d (Nl, Nc, Dir, a, t, ip, w); + + if (Reverse == True) + { + for (i=0; i < Nl; i++) + for (j=0; j < Nc; j++) Trans(i,j) = (float) (a[i][j]*Norm); + } + else + for (i=0; i < Nl; i++) + for (j=0; j < Nc; j++) Trans(i,j) = (float) a[i][j]; + + + delete [] dd; + delete [] a; + delete [] w; + delete [] t; + delete [] ip; +} + +/******************************************************************************/ + +void im_blockdct(Ifloat &Ima, Ifloat &Trans, int BlockSize, Bool Overlap, Bool WeightFirst) +{ + int Nl = Ima.nl(); + int Nc = Ima.nc(); + int i,j,Nlt, Nct, BS; + Block2D B2D,B2DTrans; + Ifloat ImaBlock, Block, TransBlock; + // B2D.Verbose = True; + B2D.BlockOverlap = Overlap; + B2D.WeightFirst = WeightFirst; + // B2DTrans.Verbose = True; + + B2D.alloc(Nl, Nc, BlockSize); + BS = B2D.block_size(); + Nlt = B2D.nbr_block_nl() * BS; + Nct = B2D.nbr_block_nc() * BS; + + if ((Trans.nl() != Nlt) || (Trans.nc() != Nct)) Trans.resize(Nlt,Nct); + ImaBlock.alloc(BS,BS,"ImaBlock"); + TransBlock.alloc(BS,BS,"ImaTransBlock"); + B2DTrans.alloc(Nlt, Nct, BlockSize); + + for (i = 0; i < B2D.nbr_block_nl(); i++) + for (j = 0; j < B2D.nbr_block_nc(); j++) + { + if (WeightFirst == False) B2D.get_block_ima(i, j, Ima, ImaBlock); + else B2D.get_block_ima(i, j, Ima, ImaBlock); + im_dct(ImaBlock,TransBlock); + B2DTrans.put_block_ima(i, j, Trans, TransBlock); + } +} + +/******************************************************************************/ + +void im_blockdct(Ifloat &Ima, Ifloat &Trans, Block2D & B2D, Block2D & B2DTrans) +{ + int i,j,Nlt, Nct, BS; + Ifloat ImaBlock, Block, TransBlock; + BS = B2D.block_size(); + Nlt = B2D.l_nl(); + Nct = B2D.l_nc(); + + if ((Trans.nl() != Nlt) || (Trans.nc() != Nct)) Trans.resize(Nlt,Nct); + ImaBlock.alloc(BS,BS,"ImaBlock"); + TransBlock.alloc(BS,BS,"ImaTransBlock"); + for (i = 0; i < B2D.nbr_block_nl(); i++) + for (j = 0; j < B2D.nbr_block_nc(); j++) + { + B2D.get_block_ima(i, j, Ima, ImaBlock); + im_dct(ImaBlock,TransBlock); + B2DTrans.put_block_ima(i, j, Trans, TransBlock); + } +} + +/******************************************************************************/ + +void im_invblockdct(Ifloat &Trans, Ifloat &Ima, int BlockSize, Bool Overlap, Bool WeightFirst) +{ + int Nl = Trans.nl(); + int Nc = Trans.nc(); + if (Overlap == True) + { + Nl /=2; + Nc /= 2; + } + im_invblockdct(Trans,Ima, BlockSize, Nl, Nc, Overlap, WeightFirst); +} + +/******************************************************************************/ + +void im_invblockdct(Ifloat &Trans, Ifloat &Ima, Block2D & B2D, Block2D & B2DTrans) +{ + int i,j, Nlt, Nct, BS; + Ifloat ImaBlock, Block, TransBlock; + int Nl = B2D.nl(); + int Nc = B2D.nc(); + BS = B2D.block_size(); + Nlt = B2D.nbr_block_nl() * BS; + Nct = B2D.nbr_block_nc() * BS; + if ((Trans.nl() != Nlt) || (Trans.nc() != Nct)) + { + cout << "Error in im_invblockdct: Nlt = " << Nlt << " Nct = " << Nct << endl; + cout << " Expected values are Nlt = " << Trans.nl() << " Nct = " << Trans.nc() << endl; + exit(-1); + } + if ((Ima.nl() != Nl) || (Ima.nc() != Nc)) Ima.resize(Nl,Nc); + ImaBlock.alloc(BS,BS,"ImaBlock"); + TransBlock.alloc(BS,BS,"ImaTransBlock"); + Ima.init(); + for (i = 0; i < B2D.nbr_block_nl(); i++) + for (j = 0; j < B2D.nbr_block_nc(); j++) + { + B2DTrans.get_block_ima(i, j, Trans, TransBlock); + im_dct(TransBlock,ImaBlock,True); + B2D.add_block_ima(i, j, Ima, ImaBlock); + } +} + +/******************************************************************************/ + +void im_invblockdct(Ifloat &Trans, Ifloat &Ima, int BlockSize, int Nl, int Nc, Bool Overlap, Bool WeightFirst) +{ + int i,j, Nlt, Nct, BS; + Block2D B2D,B2DTrans; + Ifloat ImaBlock, Block, TransBlock; + + // B2D.Verbose = True; + B2D.BlockOverlap = Overlap; + B2D.WeightFirst = WeightFirst; + B2D.alloc(Nl, Nc, BlockSize); + BS = B2D.block_size(); + Nlt = B2D.nbr_block_nl() * BS; + Nct = B2D.nbr_block_nc() * BS; + if ((Trans.nl() != Nlt) || (Trans.nc() != Nct)) + { + cout << "Error in im_invblockdct: Nlt = " << Nlt << " Nct = " << Nct << endl; + cout << " Expected values are Nlt = " << Trans.nl() << " Nct = " << Trans.nc() << endl; + exit(-1); + } + if ((Ima.nl() != Nl) || (Ima.nc() != Nc)) Ima.resize(Nl,Nc); + ImaBlock.alloc(BS,BS,"ImaBlock"); + TransBlock.alloc(BS,BS,"ImaTransBlock"); + B2DTrans.alloc(Nlt, Nct, BlockSize); + + Ima.init(); + for (i = 0; i < B2D.nbr_block_nl(); i++) + for (j = 0; j < B2D.nbr_block_nc(); j++) + { + B2DTrans.get_block_ima(i, j, Trans, TransBlock); + im_dct(TransBlock,ImaBlock,True); + B2D.add_block_ima(i, j, Ima, ImaBlock); + } +} + +/******************************************************************************/ + +void LOCAL_DCT2D::reset() +{ + Nl = Nc = Nlt = Nlt = 0; + BlockSize= 0; + Overlap= False; +} + +/******************************************************************************/ + +void LOCAL_DCT2D::alloc_from_trans(int Nli, int Nci, int BS, Bool Overlapping, Bool WeightF) +{ + int Nl1 = Nli; + int Nc1 = Nci; + if (Overlapping == True) + { + Nl1 /= 2; + Nc1 /= 2; + } + alloc(Nl1, Nc1, BS, Overlapping, WeightF); +} + +/******************************************************************************/ + +void LOCAL_DCT2D::alloc(int Nli, int Nci, int BS, Bool Overlapping, Bool WeightF) +{ + Nl= Nli; + Nc = Nci; + BlockSize = BS; + Overlap = Overlapping; + + if ((Nl == 0) || (Nc == 0)) + { + cout << "Error: LOCAL_DCT2D class is allocated for an image of zero dimension ... " << endl; + cout << " Image size is: Nl = " << Nl << " Nc = " << Nc << endl; + exit(-1); + } + + + if ((BlockSize > MIN(Nl,Nc)) || (BlockSize <= 0)) BlockSize = MIN(Nl,Nc); + if (is_power_of_2(BlockSize) == False) + { + cout << "Error: block size must be power of two ... " << endl; + cout << " BlockSize = " << BlockSize << endl; + exit(-1); + } + B2D.BlockOverlap = Overlapping; + B2D.WeightFirst = WeightF; + B2D.alloc(Nl, Nc, BlockSize); + BlockSize = B2D.block_size(); + Nlt = B2D.nbr_block_nl() * BlockSize; + Nct = B2D.nbr_block_nc() * BlockSize; + B2DTrans.alloc(Nlt, Nct, BlockSize); + DCTIma.alloc(Nlt, Nct, "DCTIma"); +} + +/******************************************************************************/ + +void LOCAL_DCT2D::transform(Ifloat &Ima) +{ + if ((Nl == 0) || (Nc == 0)) + { + cout << "Error: LOCAL_DCT2D class is not allocated ... " << endl; + exit(-1); + } + + if ((Ima.nl() != Nl) || (Ima.nc() != Nc)) + { + cout << "Error: Bad image size in LOCAL_DCT2D::transform ... " << endl; + cout << " Image size is: Nl = " << Ima.nl() << " Nc = " << Ima.nc() << endl; + cout << " Expected size is: Nl = " << Nl << " Nc = " << Nc << endl; + exit(-1); + } + im_blockdct(Ima, DCTIma, B2D, B2DTrans); +} + +/******************************************************************************/ + +void LOCAL_DCT2D::recons(Ifloat &Ima) +{ + if ((Nl == 0) || (Nc == 0)) + { + cout << "Error: LOCAL_DCT2D class is not allocated ... " << endl; + exit(-1); + } + + if ((Ima.nl() != Nl) || (Ima.nc() != Nc)) Ima.resize(Nl,Nc); + im_invblockdct(DCTIma, Ima, B2D, B2DTrans); +} + +/******************************************************************************/ + +void LOCAL_DCT2D::set_dct(Ifloat &DCTData) +{ + if ((DCTData.nl() != Nlt) || (DCTData.nc() != Nct)) + { + cout << "Error: Bad image size in LOCAL_DCT2D::set_dct ... " << endl; + cout << " Image size is: Nl = " << DCTData.nl() << " Nc = " << DCTData.nc() << endl; + cout << " Expected size is: Nl = " << Nlt << " Nc = " << Nct << endl; + exit(-1); + } + DCTIma = DCTData; +} + +/******************************************************************************/ + +float LOCAL_DCT2D::max_without_zerofreq() +{ + int bi,bj; + int BS = B2DTrans.block_size(); + Ifloat BlockCosIma(BS, BS, "blockima"); + float Max=0; + float MaxBlock; + for (bi = 0; bi < B2DTrans.nbr_block_nl(); bi++) + for (bj = 0; bj < B2DTrans.nbr_block_nc(); bj++) + { + B2DTrans.get_block_ima(bi,bj, DCTIma, BlockCosIma); + BlockCosIma(0,0) = 0.; + MaxBlock = BlockCosIma.max(); + if (Max < MaxBlock) Max = MaxBlock; + } + return Max; +} + +/******************************************************************************/ + +float LOCAL_DCT2D::maxfabs_without_zerofreq() +{ + int bi,bj; + int BS = B2DTrans.block_size(); + Ifloat BlockCosIma(BS, BS, "blockima"); + float Max=0; + float MaxBlock; + for (bi = 0; bi < B2DTrans.nbr_block_nl(); bi++) + for (bj = 0; bj < B2DTrans.nbr_block_nc(); bj++) + { + B2DTrans.get_block_ima(bi,bj, DCTIma, BlockCosIma); + BlockCosIma(0,0) = 0.; + MaxBlock = BlockCosIma.maxfabs(); + if (ABS(Max) < ABS(MaxBlock)) Max = MaxBlock; + } + return Max; +} diff --git a/src/libtools/IM_DCT.h b/src/libtools/IM_DCT.h new file mode 100644 index 0000000..8651afb --- /dev/null +++ b/src/libtools/IM_DCT.h @@ -0,0 +1,76 @@ +/****************************************************************************** +** Copyright (C) 2000 by CEA +******************************************************************************* +** +** UNIT +** +** Version: 1.0 +** +** Author: Jean-Luc Starck +** +** Date: 10/02/00 +** +** File: IM_DCT.h +** +** Modification history : +** +******************************************************************************/ + +#ifndef _DCT_H_ +#define _DCT_H_ + +#include"IM_Obj.h" +#include"IM_Block2D.h" + +void ddct16x16s(int isgn, double **a); +void ddct8x8s(int isgn, double **a); + +void ddst2d(int n1, int n2, int isgn, double **a, double *t,int *ip, double *w); +// ddst2d: Discrete Sine Transform +void ddct2d(int n1, int n2, int isgn, double **a, double *t,int *ip, double *w); +// ddct2d: Discrete Cosine Transform +void rdft2d(int n1, int n2, int isgn, double **a, double *t,int *ip, double *w); +// rdft2d: Real Discrete Fourier Transform +void cdft2d(int n1, int n2, int isgn, double **a, double *t,int *ip, double *w); +// cdft2d: Complex Discrete Fourier Transform + + +void im_dct(Ifloat &Ima, Ifloat &Trans, Bool Reverse=False); +void im_dst(Ifloat &Ima, Ifloat &Trans, Bool Reverse=False); +void im_blockdct(Ifloat &Ima, Ifloat &Trans, int BlockSize, Bool Overlap=False, Bool WeightF=False); +void im_blockdct(Ifloat &Ima, Ifloat &Trans, Block2D & B2D, Block2D & B2DTrans); +void im_invblockdct(Ifloat &Trans, Ifloat &Ima, int BlockSize, int Nl, int Nc, Bool Overlap=False, Bool WeightF=False); +void im_invblockdct(Ifloat &Trans, Ifloat &Ima, int BlockSize, Bool Overlap=False, Bool WeightF=False); +void im_invblockdct(Ifloat &Trans, Ifloat &Ima, Block2D & B2D, Block2D & B2DTrans); + +class LOCAL_DCT2D { + Bool Overlap; + int BlockSize; + int Nl, Nc; + int Nlt, Nct; + void reset(); + public: + Block2D B2D, B2DTrans; + Ifloat DCTIma; + LOCAL_DCT2D() {reset();} + ~LOCAL_DCT2D() {}; + void alloc(int Nl, int Nc, int Block_Size, Bool Overlapping=False, Bool WeightF=False); + void alloc_from_trans(int Nli, int Nci, int BS, Bool Overlapping=False, Bool WeightF=False); + int nl_ima() {return Nl;} + int nc_ima() {return Nc;} + int nl() {return Nlt;} + int nc() {return Nct;} + void set_dct(Ifloat &DCTData); + inline float & operator() (int i, int j) {return DCTIma(i,j);} + Ifloat & trans_ima() { return DCTIma;} + float * buffer() const { return DCTIma.buffer();} + int block_size() {return BlockSize;} + void transform (Ifloat &Ima); + void recons (Ifloat &Ima); + float max() { return DCTIma.max();} + float max_without_zerofreq(); + float maxfabs_without_zerofreq(); + inline float norm() { return 1.;} +}; + +#endif diff --git a/src/msvst_2d1d.cc b/src/msvst/msvst_2d1d.cc similarity index 100% rename from src/msvst_2d1d.cc rename to src/msvst/msvst_2d1d.cc diff --git a/src/msvst_iwt2d.cc b/src/msvst/msvst_iwt2d.cc similarity index 100% rename from src/msvst_iwt2d.cc rename to src/msvst/msvst_iwt2d.cc diff --git a/src/msvst_iwt2d_coupled.cc b/src/msvst/msvst_iwt2d_coupled.cc similarity index 100% rename from src/msvst_iwt2d_coupled.cc rename to src/msvst/msvst_iwt2d_coupled.cc diff --git a/src/msvst_uwt2d.cc b/src/msvst/msvst_uwt2d.cc similarity index 100% rename from src/msvst_uwt2d.cc rename to src/msvst/msvst_uwt2d.cc diff --git a/src/cb_filter.cc b/src/unused/cb_filter.cc similarity index 100% rename from src/cb_filter.cc rename to src/unused/cb_filter.cc diff --git a/src/cb_mca.cc b/src/unused/cb_mca.cc similarity index 100% rename from src/cb_mca.cc rename to src/unused/cb_mca.cc diff --git a/src/mr1d_deconv.cc b/src/unused/mr1d_deconv.cc similarity index 100% rename from src/mr1d_deconv.cc rename to src/unused/mr1d_deconv.cc