Compare revisions

Jakob Maljaars · Jakob Maljaars · df8a2004 · df8a2004 · df8a2004 · df8a2004
--- a/.clang-format
+++ b/.clang-format
+---
+Language:        Cpp
+# BasedOnStyle:  Google
+AccessModifierOffset: -1
+AlignAfterOpenBracket: Align
+AlignConsecutiveMacros: false
+AlignConsecutiveAssignments: false
+AlignConsecutiveDeclarations: false
+AlignEscapedNewlines: Left
+AlignOperands:   true
+AlignTrailingComments: true
+AllowAllArgumentsOnNextLine: true
+AllowAllConstructorInitializersOnNextLine: true
+AllowAllParametersOfDeclarationOnNextLine: true
+AllowShortBlocksOnASingleLine: false
+AllowShortCaseLabelsOnASingleLine: false
+AllowShortFunctionsOnASingleLine: All
+AllowShortLambdasOnASingleLine: All
+AllowShortIfStatementsOnASingleLine: WithoutElse
+AllowShortLoopsOnASingleLine: true
+AlwaysBreakAfterDefinitionReturnType: None
+AlwaysBreakAfterReturnType: None
+AlwaysBreakBeforeMultilineStrings: true
+AlwaysBreakTemplateDeclarations: Yes
+BinPackArguments: true
+BinPackParameters: true
+BraceWrapping:
+  AfterCaseLabel:  false
+  AfterClass:      false
+  AfterControlStatement: false
+  AfterEnum:       false
+  AfterFunction:   false
+  AfterNamespace:  false
+  AfterObjCDeclaration: false
+  AfterStruct:     false
+  AfterUnion:      false
+  AfterExternBlock: false
+  BeforeCatch:     false
+  BeforeElse:      false
+  IndentBraces:    false
+  SplitEmptyFunction: true
+  SplitEmptyRecord: true
+  SplitEmptyNamespace: true
+BreakBeforeBinaryOperators: None
+BreakBeforeBraces: Attach
+BreakBeforeInheritanceComma: false
+BreakInheritanceList: BeforeColon
+BreakBeforeTernaryOperators: true
+BreakConstructorInitializersBeforeComma: false
+BreakConstructorInitializers: BeforeColon
+BreakAfterJavaFieldAnnotations: false
+BreakStringLiterals: true
+ColumnLimit:     80
+CommentPragmas:  '^ IWYU pragma:'
+CompactNamespaces: false
+ConstructorInitializerAllOnOneLineOrOnePerLine: true
+ConstructorInitializerIndentWidth: 4
+ContinuationIndentWidth: 4
+Cpp11BracedListStyle: true
+DerivePointerAlignment: false
+DisableFormat:   false
+ExperimentalAutoDetectBinPacking: false
+FixNamespaceComments: true
+ForEachMacros:
+  - foreach
+  - Q_FOREACH
+  - BOOST_FOREACH
+IncludeBlocks:   Regroup
+IncludeCategories:
+  - Regex:           '^<ext/.*\.h>'
+    Priority:        2
+  - Regex:           '^<.*\.h>'
+    Priority:        1
+  - Regex:           '^<.*'
+    Priority:        2
+  - Regex:           '.*'
+    Priority:        3
+IncludeIsMainRegex: '([-_](test|unittest))?$'
+IndentCaseLabels: true
+IndentPPDirectives: None
+IndentWidth:     2
+IndentWrappedFunctionNames: false
+JavaScriptQuotes: Leave
+JavaScriptWrapImports: true
+KeepEmptyLinesAtTheStartOfBlocks: false
+MacroBlockBegin: ''
+MacroBlockEnd:   ''
+MaxEmptyLinesToKeep: 1
+NamespaceIndentation: None
+ObjCBinPackProtocolList: Never
+ObjCBlockIndentWidth: 2
+ObjCSpaceAfterProperty: false
+ObjCSpaceBeforeProtocolList: true
+PenaltyBreakAssignment: 2
+PenaltyBreakBeforeFirstCallParameter: 1
+PenaltyBreakComment: 300
+PenaltyBreakFirstLessLess: 120
+PenaltyBreakString: 1000
+PenaltyBreakTemplateDeclaration: 10
+PenaltyExcessCharacter: 1000000
+PenaltyReturnTypeOnItsOwnLine: 200
+PointerAlignment: Left
+RawStringFormats:
+  - Language:        Cpp
+    Delimiters:
+      - cc
+      - CC
+      - cpp
+      - Cpp
+      - CPP
+      - 'c++'
+      - 'C++'
+    CanonicalDelimiter: ''
+    BasedOnStyle:    google
+  - Language:        TextProto
+    Delimiters:
+      - pb
+      - PB
+      - proto
+      - PROTO
+    EnclosingFunctions:
+      - EqualsProto
+      - EquivToProto
+      - PARSE_PARTIAL_TEXT_PROTO
+      - PARSE_TEST_PROTO
+      - PARSE_TEXT_PROTO
+      - ParseTextOrDie
+      - ParseTextProtoOrDie
+    CanonicalDelimiter: ''
+    BasedOnStyle:    google
+ReflowComments:  true
+SortIncludes:    false
+SortUsingDeclarations: true
+SpaceAfterCStyleCast: false
+SpaceAfterLogicalNot: false
+SpaceAfterTemplateKeyword: true
+SpaceBeforeAssignmentOperators: true
+SpaceBeforeCpp11BracedList: false
+SpaceBeforeCtorInitializerColon: true
+SpaceBeforeInheritanceColon: true
+SpaceBeforeParens: ControlStatements
+SpaceBeforeRangeBasedForLoopColon: true
+SpaceInEmptyParentheses: false
+SpacesBeforeTrailingComments: 2
+SpacesInAngles:  false
+SpacesInContainerLiterals: true
+SpacesInCStyleCastParentheses: false
+SpacesInParentheses: false
+SpacesInSquareBrackets: false
+Standard:        c++17
+StatementMacros:
+  - Q_UNUSED
+  - QT_REQUIRE_VERSION
+TabWidth:        8
+UseTab:          Never
+...
--- a/.gitlab-ci.astron.yml
+++ b/.gitlab-ci.astron.yml
+# Copyright (C) 2022 ASTRON (Netherlands Institute for Radio Astronomy)
+# SPDX-License-Identifier: GPL-3.0-or-later
+
+# This file contains the pipelines that run on the Astron repository of Radler,
+# which is at https://git.astron.nl/RD/radler
+
+include: .gitlab-ci.common.yml
\ No newline at end of file
--- a/.gitlab-ci.common.yml
+++ b/.gitlab-ci.common.yml
+# Copyright (C) 2022 ASTRON (Netherlands Institute for Radio Astronomy)
+# SPDX-License-Identifier: GPL-3.0-or-later
+
+workflow:
+  rules:
+    # don't create a pipeline if its a commit pipeline, on a branch and that branch has open merge requests (we will get a MR build instead)
+    - if: $CI_PIPELINE_SOURCE == "push" && $CI_COMMIT_BRANCH && $CI_OPEN_MERGE_REQUESTS
+      when: never
+    - when: always
+
+stages:
+  - versioning
+  - prepare
+  - linting
+  - build
+  - test
+
+# The 'IMAGE' variables allow reusing docker images between different pipelines.
+# See https://confluence.skatelescope.org/display/SE/Caching+Docker+images+using+GitLab+CI+registry
+versioning:
+  stage: versioning
+  image: bitnami/git
+  script:
+    # Unshallowing ensures that 'git log' works
+    # - git fetch --unshallow
+    - git fetch --depth=1000
+    - echo BASE_IMAGE_2004=${CI_REGISTRY_IMAGE}/base_2004:$(git log -n 1 --pretty=format:%H -- docker/ubuntu_20_04_base) > versions.env
+    - echo BASE_IMAGE_2204=${CI_REGISTRY_IMAGE}/base_2204:$(git log -n 1 --pretty=format:%H -- docker/ubuntu_22_04_base) >> versions.env
+    - cat versions.env
+  artifacts:
+    reports:
+      dotenv: versions.env
+
+.prepare:
+  stage: prepare
+  needs: ["versioning"]
+  image: docker:20.10
+  services:
+    - docker:20.10-dind
+  before_script:
+    - echo $CI_REGISTRY_PASSWORD | docker login -u $CI_REGISTRY_USER --password-stdin $CI_REGISTRY
+  script:
+    - |
+      if ! docker manifest inspect $DOCKER_IMAGE > /dev/null; then
+        docker build $DOCKER_BUILD_ARG --tag $DOCKER_IMAGE -f $DOCKER_FILE .
+        docker push $DOCKER_IMAGE
+      fi
+  # Skip the job if there are no changes to the Docker file. This shortcut only
+  # works for push and merge request jobs.
+  # A manual pipeline run will thus create missing docker images.
+  rules:
+    - changes:
+      - $DOCKER_FILE
+
+prepare-base-2004:
+  extends: .prepare
+  variables:
+    DOCKER_IMAGE: $BASE_IMAGE_2004
+    DOCKER_FILE: docker/ubuntu_20_04_base
+
+prepare-base-2204:
+  extends: .prepare
+  variables:
+    DOCKER_IMAGE: $BASE_IMAGE_2204
+    DOCKER_FILE: docker/ubuntu_22_04_base
+
+.needs-2004:
+  needs:
+    - job: versioning
+    - job: prepare-base-2004
+      optional: true
+  image: $BASE_IMAGE_2004
+
+.needs-2204:
+  needs:
+    - job: versioning
+    - job: prepare-base-2204
+      optional: true
+  image: $BASE_IMAGE_2204
+
+clang-format:
+  extends: .needs-2004
+  stage: linting
+  script:
+    # Explicitly checking out submodules might become redundant
+    # once git fetch --unshallow works.
+    - git submodule update --init --recursive --checkout --depth 1
+    - ./scripts/run-format.sh
+
+.build:
+  stage: build
+  script:
+    - cmake --version
+    - mkdir build && cd build
+    - cmake -DBUILD_TESTING=ON -DCMAKE_INSTALL_PREFIX=.. -DCMAKE_CXX_FLAGS="-coverage" -DCMAKE_EXE_LINKER_FLAGS="-coverage" ..
+    - make -j`nproc`
+    - make install
+  artifacts:
+    paths:
+      - build
+
+build-2004:
+  extends: [".needs-2004",".build"]
+
+build-2204:
+  extends: [".needs-2204",".build"]
+
+.test:
+  stage: test
+  script:
+     - cd build/
+     - ctest -j`nproc` --output-on-failure -T test
+
+test-2004:
+  extends: .test
+  needs: ["versioning","build-2004"]
+  image: $BASE_IMAGE_2004
+  after_script:
+    - gcovr -j`nproc` -r .. -e '.*/external/.*' -e '.*/CompilerIdCXX/.*' -e '.*/test/.*'
+
+test-2204:
+  extends: .test
+  needs: ["versioning","build-2204"]
+  image: $BASE_IMAGE_2204
+  after_script:
+    - gcovr -j`nproc` -r .. -e '.*/external/.*' -e '.*/CompilerIdCXX/.*' -e '.*/test/.*'
--- a/.gitmodules
+++ b/.gitmodules
+[submodule "external/aocommon"]
+	path = external/aocommon
+	url = https://gitlab.com/aroffringa/aocommon.git
+[submodule "external/schaapcommon"]
+	path = external/schaapcommon
+	url = https://git.astron.nl/RD/schaapcommon.git
+[submodule "external/pybind11"]
+	path = external/pybind11
+	url = https://github.com/pybind/pybind11.git
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
+# Copyright (C) 2022 ASTRON (Netherlands Institute for Radio Astronomy)
+# SPDX-License-Identifier: GPL-3.0-or-later
+
+cmake_minimum_required(VERSION 3.8)
+
+# When Radler is compiled as an ExternalProject inside another project, set this
+# option to On. See, e.g., the wsclean CMake file for an example.
+option(COMPILE_AS_EXTERNAL_PROJECT OFF)
+
+set(RADLER_VERSION 0.0.0)
+project(radler VERSION ${RADLER_VERSION})
+if(RADLER_VERSION MATCHES "^([0-9]+)\\.([0-9]+)\\.([0-9]+)")
+  set(RADLER_VERSION_MAJOR "${CMAKE_MATCH_1}")
+  set(RADLER_VERSION_MINOR "${CMAKE_MATCH_2}")
+  set(RADLER_VERSION_PATCH "${CMAKE_MATCH_3}")
+else()
+  message(FATAL_ERROR "Failed to parse RADLER_VERSION='${RADLER_VERSION}'")
+endif()
+
+if(POLICY CMP0076)
+  cmake_policy(SET CMP0076 NEW)
+endif()
+
+set(CMAKE_CXX_STANDARD 17)
+set(CMAKE_CXX_STANDARD_REQUIRED YES)
+set(CMAKE_CXX_EXTENSIONS OFF)
+
+include(CheckCXXCompilerFlag)
+if(PORTABLE)
+  if(DEFINED TARGET_CPU)
+    message(WARNING "You have selected to build PORTABLE binaries. "
+                    "TARGET_CPU settings will be ignored.")
+    unset(TARGET_CPU CACHE)
+  endif()
+endif()
+
+if(NOT COMPILE_AS_EXTERNAL_PROJECT)
+  # Include submodules
+  set(ExternalSubmoduleDirectories aocommon pybind11 schaapcommon)
+  foreach(ExternalSubmodule ${ExternalSubmoduleDirectories})
+    if(NOT EXISTS ${CMAKE_SOURCE_DIR}/external/${ExternalSubmodule})
+      message(
+        FATAL_ERROR
+          "The external submodule '${ExternalSubmodule}' is missing in the external/ subdirectory. "
+          "This is likely the result of downloading a git tarball without submodules. "
+          "This is not supported: git tarballs do not provide the required versioning "
+          "information for the submodules. Please perform a git clone of this repository."
+      )
+    endif()
+  endforeach()
+
+  # Find and include git submodules
+  find_package(Git QUIET)
+  if(GIT_FOUND AND EXISTS "${PROJECT_SOURCE_DIR}/.git")
+    # Update submodules as needed
+    option(GIT_SUBMODULE "Check submodules during build" ON)
+    if(GIT_SUBMODULE)
+      message(STATUS "Submodule update")
+      execute_process(
+        COMMAND ${GIT_EXECUTABLE} submodule update --init --recursive --checkout
+                --depth 1
+        WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}
+        RESULT_VARIABLE GIT_SUBMOD_RESULT)
+      if(NOT GIT_SUBMOD_RESULT EQUAL "0")
+        message(
+          FATAL_ERROR
+            "git submodule update --init failed with ${GIT_SUBMOD_RESULT}, please checkout submodules"
+        )
+      endif()
+    endif()
+  endif()
+endif()
+
+if(COMPILE_AS_EXTERNAL_PROJECT)
+  message(
+    STATUS "Radler is compiled as an external project within another project.")
+  if(NOT DEFINED AOCOMMON_INCLUDE_DIR)
+    message(
+      FATAL_ERROR
+        "AOCOMMON_INCLUDE_DIR not specified. Please add -DAOCOMMON_INCLUDE_DIR to the CMAKE_ARGS."
+    )
+  endif()
+  if(NOT DEFINED SCHAAPCOMMON_SOURCE_DIR)
+    message(
+      FATAL_ERROR
+        "SCHAAPCOMMON_SOURCE_DIR not specified. Please add -DSCHAAPCOMMON_SOURCE_DIR to the CMAKE_ARGS."
+    )
+  endif()
+  if(NOT DEFINED PYBIND11_SOURCE_DIR)
+    message(
+      FATAL_ERROR
+        "PYBIND11_SOURCE_DIR not specified. Please add -DPYBIND11_SOURCE_DIR to the CMAKE_ARGS."
+    )
+  endif()
+else()
+  set(AOCOMMON_INCLUDE_DIR ${CMAKE_SOURCE_DIR}/external/aocommon/include)
+  set(SCHAAPCOMMON_SOURCE_DIR ${CMAKE_SOURCE_DIR}/external/schaapcommon)
+  set(PYBIND11_SOURCE_DIR ${CMAKE_SOURCE_DIR}/external/pybind11)
+endif()
+
+set(CMAKE_MODULE_PATH ${CMAKE_SOURCE_DIR}/cmake)
+
+# Boost date_time is needed in aocommon
+find_package(
+  Boost
+  COMPONENTS date_time
+  REQUIRED)
+
+# Threads
+find_package(Threads REQUIRED)
+
+# Find and include HDF5
+find_package(
+  HDF5
+  COMPONENTS C CXX
+  REQUIRED)
+add_definitions(${HDF5_DEFINITIONS} -DH5_USE_110_API)
+
+# Casacore
+set(CASACORE_MAKE_REQUIRED_EXTERNALS_OPTIONAL TRUE)
+find_package(Casacore REQUIRED COMPONENTS fits casa ms tables measures)
+
+# CFitsio
+find_package(CFITSIO REQUIRED)
+
+# Python3
+find_package(PythonLibs 3 REQUIRED)
+find_package(PythonInterp REQUIRED)
+message(STATUS "Using python version ${PYTHON_VERSION_STRING}")
+
+if(COMPILE_AS_EXTERNAL_PROJECT)
+  add_subdirectory(${SCHAAPCOMMON_SOURCE_DIR}
+                   ${PROJECT_BINARY_DIR}/schaapcommon)
+  add_subdirectory(${PYBIND11_SOURCE_DIR} ${PROJECT_BINARY_DIR}/pybind11)
+else()
+  add_subdirectory(${SCHAAPCOMMON_SOURCE_DIR})
+  add_subdirectory(${PYBIND11_SOURCE_DIR})
+endif()
+target_include_directories(schaapcommon SYSTEM PUBLIC ${AOCOMMON_INCLUDE_DIR})
+
+set(RADLER_TARGET_INCLUDE_DIRS
+    ${AOCOMMON_INCLUDE_DIR}
+    ${Boost_INCLUDE_DIRS}
+    ${CASACORE_INCLUDE_DIRS}
+    ${CFITSIO_INCLUDE_DIR}
+    ${HDF5_INCLUDE_DIRS}
+    ${pybind11_INCLUDE_DIRS}
+    ${SCHAAPCOMMON_SOURCE_DIR}/include)
+
+set(RADLER_TARGET_LIBS
+    ${Boost_DATE_TIME_LIBRARY} ${CMAKE_THREAD_LIBS_INIT} ${CASACORE_LIBRARIES}
+    ${CFITSIO_LIBRARY} pybind11::embed schaapcommon)
+
+# Source directories
+add_subdirectory(cpp)
+
+# Compile tests
+if(NOT ${COMPILE_AS_EXTERNAL_PROJECT} AND ${BUILD_TESTING})
+  include(CTest)
+  find_package(
+    Boost
+    COMPONENTS unit_test_framework
+    REQUIRED)
+  add_subdirectory(cpp/test)
+  add_subdirectory(cpp/algorithms/test)
+  add_subdirectory(cpp/math/test)
+
+  # demo targets are excluded from a default build
+  add_subdirectory(cpp/demo)
+endif()
--- a/LICENSE
+++ b/LICENSE
--- a/README.md
+++ b/README.md
-# Radio Astronomical DeconvolvER
+# Radio Astronomical Deconvolution Library

-RADLER is a library providing functionality for deconvolving astronomical images. RADLER evolved as a stand-alone library from the w-stacking clean (WSClean) imager https://gitlab.com/aroffringa/wsclean.
+Radler is a library providing functionality for deconvolving astronomical images. Radler evolved as a stand-alone library from the w-stacking clean (WSClean) imager https://gitlab.com/aroffringa/wsclean.
--- a/cmake/FindCFITSIO.cmake
+++ b/cmake/FindCFITSIO.cmake
+# * Try to find CFITSIO. Variables used by this module: CFITSIO_ROOT_DIR     -
+#   CFITSIO root directory Variables defined by this module: CFITSIO_FOUND -
+#   system has CFITSIO CFITSIO_INCLUDE_DIR  - the CFITSIO include directory
+#   (cached) CFITSIO_INCLUDE_DIRS - the CFITSIO include directories (identical
+#   to CFITSIO_INCLUDE_DIR) CFITSIO_LIBRARY      - the CFITSIO library (cached)
+#   CFITSIO_LIBRARIES    - the CFITSIO libraries (identical to CFITSIO_LIBRARY)
+#   CFITSIO_VERSION_STRING the found version of CFITSIO, padded to 3 digits
+
+# Copyright (C) 2022 ASTRON (Netherlands Institute for Radio Astronomy)
+# SPDX-License-Identifier: GPL-3.0-or-later
+
+if(NOT CFITSIO_FOUND)
+
+  find_path(
+    CFITSIO_INCLUDE_DIR fitsio.h
+    HINTS ${CFITSIO_ROOT_DIR}
+    PATH_SUFFIXES include include/cfitsio include/libcfitsio0)
+
+  if(CFITSIO_INCLUDE_DIR)
+    file(READ "${CFITSIO_INCLUDE_DIR}/fitsio.h" CFITSIO_H)
+    set(CFITSIO_VERSION_REGEX
+        ".*#define CFITSIO_VERSION[^0-9]*([0-9]+)\\.([0-9]+).*")
+    if("${CFITSIO_H}" MATCHES ${CFITSIO_VERSION_REGEX})
+      # Pad CFITSIO minor version to three digit because 3.181 is older than
+      # 3.35
+      string(REGEX REPLACE ${CFITSIO_VERSION_REGEX} "\\1.\\200"
+                           CFITSIO_VERSION_STRING "${CFITSIO_H}")
+      string(SUBSTRING ${CFITSIO_VERSION_STRING} 0 5 CFITSIO_VERSION_STRING)
+      string(REGEX REPLACE "^([0-9]+)[.]([0-9]+)" "\\1" CFITSIO_VERSION_MAJOR
+                           ${CFITSIO_VERSION_STRING})
+      # CFITSIO_VERSION_MINOR will contain 80 for 3.08, 181 for 3.181 and 200
+      # for 3.2
+      string(REGEX REPLACE "^([0-9]+)[.]0*([0-9]+)" "\\2" CFITSIO_VERSION_MINOR
+                           ${CFITSIO_VERSION_STRING})
+    else()
+      set(CFITSIO_VERSION_STRING "Unknown")
+    endif()
+  endif(CFITSIO_INCLUDE_DIR)
+
+  find_library(
+    CFITSIO_LIBRARY cfitsio
+    HINTS ${CFITSIO_ROOT_DIR}
+    PATH_SUFFIXES lib)
+  find_library(M_LIBRARY m)
+  mark_as_advanced(CFITSIO_INCLUDE_DIR CFITSIO_LIBRARY M_LIBRARY)
+
+  include(FindPackageHandleStandardArgs)
+  find_package_handle_standard_args(
+    CFITSIO
+    REQUIRED_VARS CFITSIO_LIBRARY M_LIBRARY CFITSIO_INCLUDE_DIR
+    VERSION_VAR CFITSIO_VERSION_STRING)
+
+  set(CFITSIO_INCLUDE_DIRS ${CFITSIO_INCLUDE_DIR})
+  set(CFITSIO_LIBRARIES ${CFITSIO_LIBRARY} ${M_LIBRARY})
+
+endif(NOT CFITSIO_FOUND)
--- a/cmake/FindCasacore.cmake
+++ b/cmake/FindCasacore.cmake
+# * Try to find Casacore include dirs and libraries Usage: find_package(Casacore
+#   [REQUIRED] [COMPONENTS components...]) Valid components are: casa,
+#   coordinates, derivedmscal, fits, images, lattices, meas, measures, mirlib,
+#   ms, msfits, python, scimath, scimath_f, tables
+#
+# Note that most components are dependent on other (more basic) components. In
+# that case, it suffices to specify the "top-level" components; dependent
+# components will be searched for automatically.
+#
+# The dependency tree can be generated using the script get_casacore_deps.sh.
+# For this, you need to have a complete casacore installation, built with shared
+# libraries, at your disposal.
+#
+# The dependencies in this macro were generated against casacore release 1.7.0.
+#
+# Variables used by this module: CASACORE_ROOT_DIR         - Casacore root
+# directory.
+#
+# Variables defined by this module: CASACORE_FOUND            - System has
+# Casacore, which means that the include dir was found, as well as all libraries
+# specified (not cached) CASACORE_INCLUDE_DIR      - Casacore include directory
+# (cached) CASACORE_INCLUDE_DIRS     - Casacore include directories (not cached)
+# identical to CASACORE_INCLUDE_DIR CASACORE_LIBRARIES        - The Casacore
+# libraries (not cached) CASA_${COMPONENT}_LIBRARY - The absolute path of
+# Casacore library "component" (cached) HAVE_AIPSPP               - True if
+# system has Casacore (cached) for backward compatibility with AIPS++
+# HAVE_CASACORE             - True if system has Casacore (cached) identical to
+# CASACORE_FOUND TAQL_EXECUTABLE           - The absolute path of the TaQL
+# executable (cached)
+#
+# ATTENTION: The component names need to be in lower case, just as the casacore
+# library names. However, the CMake variables use all upper case.
+
+# Copyright (C) 2009 ASTRON (Netherlands Institute for Radio Astronomy) P.O.Box
+# 2, 7990 AA Dwingeloo, The Netherlands
+#
+# This file is part of the LOFAR software suite. The LOFAR software suite is
+# free software: you can redistribute it and/or modify it under the terms of the
+# GNU General Public License as published by the Free Software Foundation,
+# either version 3 of the License, or (at your option) any later version.
+#
+# The LOFAR software suite is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+# or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for
+# more details.
+#
+# You should have received a copy of the GNU General Public License along with
+# the LOFAR software suite. If not, see <http://www.gnu.org/licenses/>.
+#
+# $Id: FindCasacore.cmake 31487 2015-04-16 11:28:17Z dijkema $
+
+# * casacore_resolve_dependencies(_result)
+#
+# Resolve the Casacore library dependencies for the given components. The list
+# of dependent libraries will be returned in the variable result. It is sorted
+# from least dependent to most dependent library, so it can be directly fed to
+# the linker.
+#
+# Usage: casacore_resolve_dependencies(result components...)
+#
+macro(casacore_resolve_dependencies _result)
+  set(${_result} ${ARGN})
+  set(_index 0)
+  # Do a breadth-first search through the dependency graph; append to the result
+  # list the dependent components for each item in that list. Duplicates will be
+  # removed later.
+  while(1)
+    list(LENGTH ${_result} _length)
+    if(NOT _index LESS _length)
+      break()
+    endif(NOT _index LESS _length)
+    list(GET ${_result} ${_index} item)
+    list(APPEND ${_result} ${Casacore_${item}_DEPENDENCIES})
+    math(EXPR _index "${_index}+1")
+  endwhile(1)
+  # Remove all duplicates in the current result list, while retaining only the
+  # last of each duplicate.
+  list(REVERSE ${_result})
+  list(REMOVE_DUPLICATES ${_result})
+  list(REVERSE ${_result})
+endmacro(casacore_resolve_dependencies _result)
+
+# * casacore_find_library(_name)
+#
+# Search for the library ${_name}. If library is found, add it to
+# CASACORE_LIBRARIES; if not, add ${_name} to CASACORE_MISSING_COMPONENTS and
+# set CASACORE_FOUND to false.
+#
+# Usage: casacore_find_library(name)
+#
+macro(casacore_find_library _name)
+  string(TOUPPER ${_name} _NAME)
+  find_library(
+    ${_NAME}_LIBRARY ${_name}
+    HINTS ${CASACORE_ROOT_DIR}
+    PATH_SUFFIXES lib)
+  mark_as_advanced(${_NAME}_LIBRARY)
+  if(${_NAME}_LIBRARY)
+    list(APPEND CASACORE_LIBRARIES ${${_NAME}_LIBRARY})
+  else(${_NAME}_LIBRARY)
+    set(CASACORE_FOUND FALSE)
+    list(APPEND CASACORE_MISSING_COMPONENTS ${_name})
+  endif(${_NAME}_LIBRARY)
+endmacro(casacore_find_library _name)
+
+# * casacore_find_package(_name)
+#
+# Search for the package ${_name}. If the package is found, add the contents of
+# ${_name}_INCLUDE_DIRS to CASACORE_INCLUDE_DIRS and ${_name}_LIBRARIES to
+# CASACORE_LIBRARIES.
+#
+# If Casacore itself is required, then, strictly speaking, the packages it
+# requires must be present. However, when linking against static libraries they
+# may not be needed. One can override the REQUIRED setting by switching
+# CASACORE_MAKE_REQUIRED_EXTERNALS_OPTIONAL to ON. Beware that this might cause
+# compile and/or link errors.
+#
+# Usage: casacore_find_package(name [REQUIRED])
+#
+macro(casacore_find_package _name)
+  if("${ARGN}" MATCHES "^REQUIRED$"
+     AND Casacore_FIND_REQUIRED
+     AND NOT CASACORE_MAKE_REQUIRED_EXTERNALS_OPTIONAL)
+    find_package(${_name} REQUIRED)
+  else()
+    find_package(${_name})
+  endif()
+  if(${_name}_FOUND)
+    list(APPEND CASACORE_INCLUDE_DIRS ${${_name}_INCLUDE_DIRS})
+    list(APPEND CASACORE_LIBRARIES ${${_name}_LIBRARIES})
+  endif(${_name}_FOUND)
+endmacro(casacore_find_package _name)
+
+# Define the Casacore components.
+set(Casacore_components
+    casa
+    coordinates
+    derivedmscal
+    fits
+    images
+    lattices
+    meas
+    measures
+    mirlib
+    ms
+    msfits
+    python
+    scimath
+    scimath_f
+    tables)
+
+# Define the Casacore components' inter-dependencies.
+set(Casacore_casa_DEPENDENCIES)
+set(Casacore_coordinates_DEPENDENCIES fits measures casa)
+set(Casacore_derivedmscal_DEPENDENCIES ms measures tables casa)
+set(Casacore_fits_DEPENDENCIES measures tables casa)
+set(Casacore_images_DEPENDENCIES
+    mirlib
+    lattices
+    coordinates
+    fits
+    measures
+    scimath
+    tables
+    casa)
+set(Casacore_lattices_DEPENDENCIES tables scimath casa)
+set(Casacore_meas_DEPENDENCIES measures tables casa)
+set(Casacore_measures_DEPENDENCIES tables casa)
+set(Casacore_mirlib_DEPENDENCIES)
+set(Casacore_ms_DEPENDENCIES measures scimath tables casa)
+set(Casacore_msfits_DEPENDENCIES ms fits measures tables casa)
+set(Casacore_python_DEPENDENCIES casa)
+set(Casacore_scimath_DEPENDENCIES scimath_f casa)
+set(Casacore_scimath_f_DEPENDENCIES)
+set(Casacore_tables_DEPENDENCIES casa)
+
+# Initialize variables.
+set(CASACORE_FOUND FALSE)
+set(CASACORE_DEFINITIONS)
+set(CASACORE_LIBRARIES)
+set(CASACORE_MISSING_COMPONENTS)
+
+# Search for the header file first.
+if(NOT CASACORE_INCLUDE_DIR)
+  find_path(
+    CASACORE_INCLUDE_DIR casacore/casa/aips.h
+    HINTS ${CASACORE_ROOT_DIR}
+    PATH_SUFFIXES include)
+  mark_as_advanced(CASACORE_INCLUDE_DIR)
+endif(NOT CASACORE_INCLUDE_DIR)
+
+# Fallback for systems that have old casacore installed in directory not called
+# 'casacore' This fallback can be removed once we move to casacore 2.0 which
+# always puts headers in 'casacore'
+if(NOT CASACORE_INCLUDE_DIR)
+  find_path(
+    CASACORE_INCLUDE_DIR casa/aips.h
+    HINTS ${CASACORE_ROOT_DIR}
+    PATH_SUFFIXES include)
+  mark_as_advanced(CASACORE_INCLUDE_DIR)
+endif(NOT CASACORE_INCLUDE_DIR)
+
+if(NOT CASACORE_INCLUDE_DIR)
+  set(CASACORE_ERROR_MESSAGE
+      "Casacore: unable to find the header file casa/aips.h.\nPlease set CASACORE_ROOT_DIR to the root directory containing Casacore."
+  )
+else(NOT CASACORE_INCLUDE_DIR)
+  # We've found the header file; let's continue.
+  set(CASACORE_FOUND TRUE)
+  # Note that new Casacore uses #include<casacore/casa/...>, while LOFAR still
+  # uses #include<casa/...>. Hence use both in -I path.
+  set(CASACORE_INCLUDE_DIRS ${CASACORE_INCLUDE_DIR}
+                            ${CASACORE_INCLUDE_DIR}/casacore)
+
+  # Search for some often used binaries.
+  find_program(TAQL_EXECUTABLE taql HINTS ${CASACORE_ROOT_DIR}/bin)
+  mark_as_advanced(TAQL_EXECUTABLE)
+
+  # If the user specified components explicity, use that list; otherwise we'll
+  # assume that the user wants to use all components.
+  if(NOT Casacore_FIND_COMPONENTS)
+    set(Casacore_FIND_COMPONENTS ${Casacore_components})
+  endif(NOT Casacore_FIND_COMPONENTS)
+
+  # Get a list of all dependent Casacore libraries that need to be found.
+  casacore_resolve_dependencies(_find_components ${Casacore_FIND_COMPONENTS})
+
+  # Find the library for each component, and handle external dependencies
+  foreach(_comp ${_find_components})
+    casacore_find_library(casa_${_comp})
+    if(${_comp} STREQUAL casa)
+      casacore_find_package(HDF5)
+      casacore_find_library(m)
+      list(APPEND CASACORE_LIBRARIES ${CMAKE_DL_LIBS})
+    elseif(${_comp} STREQUAL coordinates)
+      casacore_find_package(WCSLIB REQUIRED)
+    elseif(${_comp} STREQUAL fits)
+      casacore_find_package(CFITSIO REQUIRED)
+    elseif(${_comp} STREQUAL scimath_f)
+      casacore_find_package(LAPACK REQUIRED)
+    endif(${_comp} STREQUAL casa)
+  endforeach(_comp ${_find_components})
+endif(NOT CASACORE_INCLUDE_DIR)
+
+# Set HAVE_CASACORE; and HAVE_AIPSPP (for backward compatibility with AIPS++).
+if(CASACORE_FOUND)
+  set(HAVE_CASACORE
+      TRUE
+      CACHE INTERNAL "Define if Casacore is installed")
+  set(HAVE_AIPSPP
+      TRUE
+      CACHE INTERNAL "Define if AIPS++/Casacore is installed")
+endif(CASACORE_FOUND)
+
+# Compose diagnostic message if not all necessary components were found.
+if(CASACORE_MISSING_COMPONENTS)
+  set(CASACORE_ERROR_MESSAGE
+      "Casacore: the following components could not be found:\n     ${CASACORE_MISSING_COMPONENTS}"
+  )
+endif(CASACORE_MISSING_COMPONENTS)
+
+# Print diagnostics.
+if(CASACORE_FOUND)
+  if(NOT Casacore_FIND_QUIETLY)
+    message(STATUS "Found the following Casacore components: ")
+    foreach(_comp ${_find_components})
+      string(TOUPPER casa_${_comp} _COMP)
+      message(STATUS "  ${_comp}: ${${_COMP}_LIBRARY}")
+    endforeach(_comp ${_find_components})
+  endif(NOT Casacore_FIND_QUIETLY)
+else(CASACORE_FOUND)
+  if(Casacore_FIND_REQUIRED)
+    message(FATAL_ERROR "${CASACORE_ERROR_MESSAGE}")
+  else(Casacore_FIND_REQUIRED)
+    message(STATUS "${CASACORE_ERROR_MESSAGE}")
+  endif(Casacore_FIND_REQUIRED)
+endif(CASACORE_FOUND)
--- a/cmake/config/radler-config-version.cmake.in
+++ b/cmake/config/radler-config-version.cmake.in
+# Copyright (C) 2022 ASTRON (Netherlands Institute for Radio Astronomy)
+# SPDX-License-Identifier: GPL-3.0-or-later
+
+set(PACKAGE_VERSION "@RADLER_VERSION@")
+
+# Check whether the requested PACKAGE_FIND_VERSION is compatible
+if("${PACKAGE_VERSION}" VERSION_LESS "${PACKAGE_FIND_VERSION}")
+  set(PACKAGE_VERSION_COMPATIBLE FALSE)
+else()
+  set(PACKAGE_VERSION_COMPATIBLE TRUE)
+  if ("${PACKAGE_VERSION}" VERSION_EQUAL "${PACKAGE_FIND_VERSION}")
+    set(PACKAGE_VERSION_EXACT TRUE)
+  endif()
+endif()
\ No newline at end of file
--- a/cmake/config/radler-config.cmake.in
+++ b/cmake/config/radler-config.cmake.in
+# Copyright (C) 2022 ASTRON (Netherlands Institute for Radio Astronomy)
+# SPDX-License-Identifier: GPL-3.0-or-later
+#
+#  Config file for the radler library, it sets the following variables
+#
+# - RADLER_FOUND
+# - RADLER_ROOT_DIR
+# - RADLER_INCLUDE_DIR
+# - RADLER_INCLUDE_DIRS (equals RADLER_INCLUDE_DIR)
+# - RADLER_LIB_PATH
+# - RADLER_LIB
+# - RADLER_VERSION[_MAJOR/_MINOR/_PATCH]
+
+# Compute paths
+get_filename_component(_RADLER_CMAKE_DIR "${CMAKE_CURRENT_LIST_FILE}" PATH)
+get_filename_component(_RADLER_CMAKE_DIR_ABS "${_RADLER_CMAKE_DIR}" ABSOLUTE)
+get_filename_component(_RADLER_ROOT_DIR "${_RADLER_CMAKE_DIR_ABS}/../.." ABSOLUTE)
+
+# Use FORCE to override cached variables
+set(RADLER_ROOT_DIR "${_RADLER_ROOT_DIR}"
+    CACHE PATH "Radler root (prefix) directory" FORCE)
+
+set(RADLER_INCLUDE_DIR "${RADLER_ROOT_DIR}/include"
+    CACHE PATH "Radler include directory" FORCE)
+
+set(RADLER_INCLUDE_DIRS ${RADLER_INCLUDE_DIR})
+
+set(RADLER_LIB_PATH "${RADLER_ROOT_DIR}/lib"
+    CACHE PATH "Radler library directory" FORCE)
+
+find_library(RADLER_LIB radler PATH ${RADLER_LIB_PATH} NO_DEFAULT_PATH
+             DOC "Radler library directory")
+message(STATUS "Found Radler @RADLER_VERSION@.")
+message(STATUS "  Radler include dir: ${RADLER_INCLUDE_DIR}")
+message(STATUS "  Radler lib: ${RADLER_LIB}")
+
+set(RADLER_VERSION "@RADLER_VERSION@")
+set(RADLER_VERSION_MAJOR @RADLER_VERSION_MAJOR@)
+set(RADLER_VERSION_MINOR @RADLER_VERSION_MINOR@)
+set(RADLER_VERSION_PATCH @RADLER_VERSION_PATCH@)
+set(RADLER_FOUND 1)
+
+unset(_RADLER_ROOT_DIR)
+unset(_RADLER_CMAKE_DIR)
+unset(_RADLER_CMAKE_DIR_ABS)
--- a/cmake/unittest.cmake
+++ b/cmake/unittest.cmake
+# Copyright (C) 2022 ASTRON (Netherlands Institute for Radio Astronomy)
+# SPDX-License-Identifier: GPL-3.0-or-later
+
+# Generic function for adding a radler unittest. It creates a 'test' that builds
+# the unittest and a test that runs it. Arguments:
+#
+# * First argument: Module name, e.g., 'math'.
+# * Next arguments: Source file names.
+#
+# Return value:
+#
+# * Sets TEST_NAME to the unit test name in the parent scope.
+
+function(add_unittest MODULE_NAME)
+  set(TEST_NAME "unittests_${MODULE_NAME}")
+  set(TEST_NAME
+      ${TEST_NAME}
+      PARENT_SCOPE)
+  set(FILENAMES ${ARGN})
+
+  # Add boost dynamic link flag for all test files.
+  # https://www.boost.org/doc/libs/1_66_0/libs/test/doc/html/boost_test/usage_variants.html
+  # Without this flag, linking is incorrect and boost performs duplicate
+  # delete() calls after running all tests, in the cleanup phase.
+  set_source_files_properties(${FILENAMES} PROPERTIES COMPILE_DEFINITIONS
+                                                      "BOOST_TEST_DYN_LINK")
+
+  add_executable(${TEST_NAME} ${FILENAMES})
+  target_link_libraries(${TEST_NAME} radler
+                        ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY})
+  target_include_directories(
+    ${TEST_NAME} PRIVATE $<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/cpp>)
+
+  # Add test for automatically (re)building the test if needed.
+  add_test(build_${TEST_NAME} ${CMAKE_COMMAND} --build ${CMAKE_BINARY_DIR}
+           --target ${TEST_NAME})
+  set_tests_properties(build_${TEST_NAME} PROPERTIES FIXTURES_SETUP
+                                                     ${TEST_NAME})
+
+  add_test(NAME ${TEST_NAME} COMMAND ${TEST_NAME} -f JUNIT -k ${TEST_NAME}.xml
+                                     --catch_system_error=yes)
+  set_tests_properties(${TEST_NAME} PROPERTIES LABELS unit FIXTURES_REQUIRED
+                                               ${TEST_NAME})
+
+endfunction()
--- a/cpp/CMakeLists.txt
+++ b/cpp/CMakeLists.txt
+# Copyright (C) 2022 ASTRON (Netherlands Institute for Radio Astronomy)
+# SPDX-License-Identifier: GPL-3.0-or-later
+
+# CMake function to keep directory structure when installing headers.
+function(install_headers_with_directory HEADER_LIST)
+  foreach(HEADER ${HEADER_LIST})
+    string(REGEX MATCH ".*\/" DIR ${HEADER})
+    install(FILES ${HEADER}
+            DESTINATION ${CMAKE_INSTALL_PREFIX}/include/radler/${DIR})
+  endforeach(HEADER)
+endfunction(install_headers_with_directory)
+
+add_library(${PROJECT_NAME} SHARED "")
+
+set(RADLER_FILES
+    radler.cc
+    deconvolution_table.cc
+    component_list.cc
+    image_set.cc
+    algorithms/deconvolution_algorithm.cc
+    algorithms/generic_clean.cc
+    algorithms/iuwt_deconvolution_algorithm.cc
+    # algorithms/ls_deconvolution.cc // TODO: Complete or remove this class.
+    algorithms/more_sane.cc
+    algorithms/multiscale_algorithm.cc
+    algorithms/parallel_deconvolution.cc
+    algorithms/python_deconvolution.cc
+    algorithms/simple_clean.cc
+    algorithms/subminor_loop.cc
+    algorithms/threaded_deconvolution_tools.cc
+    algorithms/iuwt/image_analysis.cc
+    algorithms/iuwt/iuwt_decomposition.cc
+    algorithms/iuwt/iuwt_mask.cc
+    algorithms/multiscale/multiscale_transforms.cc
+    math/peak_finder.cc
+    math/rms_image.cc
+    utils/casa_mask_reader.cc)
+
+# A number of files perform the 'core' high-performance floating point
+# operations. In these files, NaNs are avoided and thus -ffast-math is allowed.
+# Note that visibilities can be NaN hence this can not be turned on for all
+# files.
+set_source_files_properties(
+  image_set.cpp
+  algorithms/generic_clean.cpp
+  algorithms/multiscale_algorithm.cpp
+  algorithms/threaded_deconvolution_tools.cpp
+  algorithms/simple_clean.cpp
+  algorithms/subminor_loop.cpp
+  algorithms/multiscale/multiscale_transforms.cpp
+  PROPERTIES COMPILE_FLAGS -ffast-math)
+
+# Using pybind11 requires using -fvisibility=hidden. See
+# https://pybind11.readthedocs.io/en/stable/faq.html
+set_source_files_properties(algorithms/python_deconvolution.cpp
+                            PROPERTIES COMPILE_FLAGS -fvisibility=hidden)
+
+target_sources(${PROJECT_NAME} PRIVATE ${RADLER_FILES})
+target_link_libraries(${PROJECT_NAME} ${RADLER_TARGET_LIBS})
+target_include_directories(${PROJECT_NAME} SYSTEM
+                           PUBLIC ${RADLER_TARGET_INCLUDE_DIRS})
+
+# Allows including the paths relative to base algorithm, in line with Google
+# style
+# https://google.github.io/styleguide/cppguide.html#Names_and_Order_of_Includes
+target_include_directories(
+  ${PROJECT_NAME} PRIVATE $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}>)
+
+target_compile_options(radler PRIVATE -O3 -Wall -Wzero-as-null-pointer-constant)
+
+if(NOT PORTABLE)
+  if(DEFINED TARGET_CPU)
+    target_compile_options(radler BEFORE PRIVATE -march=${TARGET_CPU})
+  else()
+    check_cxx_compiler_flag("-march=native" COMPILER_HAS_MARCH_NATIVE)
+    if(COMPILER_HAS_MARCH_NATIVE)
+      target_compile_options(radler BEFORE PRIVATE -march=native)
+    else()
+      message(
+        WARNING "The compiler doesn't support -march=native for your CPU.")
+    endif()
+  endif()
+endif()
+
+if(NOT COMPILE_AS_EXTERNAL_PROJECT)
+  include(GNUInstallDirs)
+  install(TARGETS ${PROJECT_NAME}
+          LIBRARY DESTINATION ${CMAKE_INSTALL_FULL_LIBDIR})
+else()
+  install(TARGETS ${PROJECT_NAME}
+          LIBRARY DESTINATION ${CMAKE_INSTALL_PREFIX}/lib)
+endif()
+
+set(PUBLIC_HEADERS
+    component_list.h
+    radler.h
+    deconvolution_settings.h
+    deconvolution_table.h
+    deconvolution_table_entry.h
+    image_set.h
+    algorithms/multiscale/multiscale_transforms.h)
+
+install_headers_with_directory("${PUBLIC_HEADERS}")
+
+if(NOT COMPILE_AS_EXTERNAL_PROJECT)
+  configure_file(${PROJECT_SOURCE_DIR}/cmake/config/radler-config.cmake.in
+                 ${PROJECT_BINARY_DIR}/CMakeFiles/radler-config.cmake @ONLY)
+  configure_file(
+    ${PROJECT_SOURCE_DIR}/cmake/config/radler-config-version.cmake.in
+    ${PROJECT_BINARY_DIR}/CMakeFiles/radler-config-version.cmake @ONLY)
+
+  # Install configuration files
+  install(FILES ${PROJECT_BINARY_DIR}/CMakeFiles/radler-config.cmake
+                ${PROJECT_BINARY_DIR}/CMakeFiles/radler-config-version.cmake
+          DESTINATION ${CMAKE_INSTALL_PREFIX}/share/radler)
+endif()
--- a/cpp/algorithms/deconvolution_algorithm.cc
+++ b/cpp/algorithms/deconvolution_algorithm.cc
+// Copyright (C) 2022 ASTRON (Netherlands Institute for Radio Astronomy)
+// SPDX-License-Identifier: GPL-3.0-or-later
+
+#include "algorithms/deconvolution_algorithm.h"
+
+#include <aocommon/system.h>
+
+namespace radler::algorithms {
+
+DeconvolutionAlgorithm::DeconvolutionAlgorithm()
+    : _threshold(0.0),
+      _majorIterThreshold(0.0),
+      _gain(0.1),
+      _mGain(1.0),
+      _cleanBorderRatio(0.05),
+      _maxIter(500),
+      _iterationNumber(0),
+      _threadCount(aocommon::system::ProcessorCount()),
+      _allowNegativeComponents(true),
+      _stopOnNegativeComponent(false),
+      _cleanMask(nullptr),
+      _logReceiver(nullptr),
+      _spectralFitter(schaapcommon::fitters::SpectralFittingMode::NoFitting,
+                      0) {}
+
+void DeconvolutionAlgorithm::ResizeImage(float* dest, size_t newWidth,
+                                         size_t newHeight, const float* source,
+                                         size_t width, size_t height) {
+  size_t srcStartX = (width - newWidth) / 2,
+         srcStartY = (height - newHeight) / 2;
+  for (size_t y = 0; y != newHeight; ++y) {
+    float* destPtr = dest + y * newWidth;
+    const float* srcPtr = source + (y + srcStartY) * width + srcStartX;
+    std::copy_n(srcPtr, newWidth, destPtr);
+  }
+}
+
+void DeconvolutionAlgorithm::RemoveNaNsInPSF(float* psf, size_t width,
+                                             size_t height) {
+  float* endPtr = psf + width * height;
+  while (psf != endPtr) {
+    if (!std::isfinite(*psf)) *psf = 0.0;
+    ++psf;
+  }
+}
+
+void DeconvolutionAlgorithm::PerformSpectralFit(float* values, size_t x,
+                                                size_t y) const {
+  _spectralFitter.FitAndEvaluate(values, x, y, _fittingScratch);
+}
+}  // namespace radler::algorithms
\ No newline at end of file
--- a/cpp/algorithms/deconvolution_algorithm.h
+++ b/cpp/algorithms/deconvolution_algorithm.h
+// Copyright (C) 2022 ASTRON (Netherlands Institute for Radio Astronomy)
+// SPDX-License-Identifier: GPL-3.0-or-later
+
+#ifndef RADLER_ALGORITHMS_DECONVOLUTION_ALGORITHM_H_
+#define RADLER_ALGORITHMS_DECONVOLUTION_ALGORITHM_H_
+
+#include <string>
+#include <cmath>
+
+#include <aocommon/image.h>
+#include <aocommon/logger.h>
+#include <aocommon/polarization.h>
+#include <aocommon/uvector.h>
+
+#include <schaapcommon/fitters/spectralfitter.h>
+
+#include "image_set.h"
+
+namespace radler::algorithms {
+
+class DeconvolutionAlgorithm {
+ public:
+  virtual ~DeconvolutionAlgorithm() {}
+
+  virtual float ExecuteMajorIteration(
+      ImageSet& dataImage, ImageSet& modelImage,
+      const std::vector<aocommon::Image>& psfImages,
+      bool& reachedMajorThreshold) = 0;
+
+  virtual std::unique_ptr<DeconvolutionAlgorithm> Clone() const = 0;
+
+  void SetMaxNIter(size_t nIter) { _maxIter = nIter; }
+
+  void SetThreshold(float threshold) { _threshold = threshold; }
+
+  void SetMajorIterThreshold(float mThreshold) {
+    _majorIterThreshold = mThreshold;
+  }
+
+  void SetGain(float gain) { _gain = gain; }
+
+  void SetMGain(float mGain) { _mGain = mGain; }
+
+  void SetAllowNegativeComponents(bool allowNegativeComponents) {
+    _allowNegativeComponents = allowNegativeComponents;
+  }
+
+  void SetStopOnNegativeComponents(bool stopOnNegative) {
+    _stopOnNegativeComponent = stopOnNegative;
+  }
+
+  void SetCleanBorderRatio(float borderRatio) {
+    _cleanBorderRatio = borderRatio;
+  }
+
+  void SetThreadCount(size_t threadCount) { _threadCount = threadCount; }
+
+  void SetLogReceiver(aocommon::LogReceiver& receiver) {
+    _logReceiver = &receiver;
+  }
+
+  size_t MaxNIter() const { return _maxIter; }
+  float Threshold() const { return _threshold; }
+  float MajorIterThreshold() const { return _majorIterThreshold; }
+  float Gain() const { return _gain; }
+  float MGain() const { return _mGain; }
+  float CleanBorderRatio() const { return _cleanBorderRatio; }
+  bool AllowNegativeComponents() const { return _allowNegativeComponents; }
+  bool StopOnNegativeComponents() const { return _stopOnNegativeComponent; }
+
+  void SetCleanMask(const bool* cleanMask) { _cleanMask = cleanMask; }
+
+  size_t IterationNumber() const { return _iterationNumber; }
+
+  void SetIterationNumber(size_t iterationNumber) {
+    _iterationNumber = iterationNumber;
+  }
+
+  static void ResizeImage(float* dest, size_t newWidth, size_t newHeight,
+                          const float* source, size_t width, size_t height);
+
+  static void RemoveNaNsInPSF(float* psf, size_t width, size_t height);
+
+  void CopyConfigFrom(const DeconvolutionAlgorithm& source) {
+    _threshold = source._threshold;
+    _gain = source._gain;
+    _mGain = source._mGain;
+    _cleanBorderRatio = source._cleanBorderRatio;
+    _maxIter = source._maxIter;
+    // skip _iterationNumber
+    _allowNegativeComponents = source._allowNegativeComponents;
+    _stopOnNegativeComponent = source._stopOnNegativeComponent;
+    _cleanMask = source._cleanMask;
+    _spectralFitter = source._spectralFitter;
+  }
+
+  void SetSpectralFittingMode(schaapcommon::fitters::SpectralFittingMode mode,
+                              size_t nTerms) {
+    _spectralFitter.SetMode(mode, nTerms);
+  }
+
+  void SetSpectrallyForcedImages(std::vector<aocommon::Image>&& images) {
+    _spectralFitter.SetForcedImages(std::move(images));
+  }
+
+  void InitializeFrequencies(const aocommon::UVector<double>& frequencies,
+                             const aocommon::UVector<float>& weights) {
+    _spectralFitter.SetFrequencies(frequencies.data(), weights.data(),
+                                   frequencies.size());
+  }
+
+  const schaapcommon::fitters::SpectralFitter& Fitter() const {
+    return _spectralFitter;
+  }
+
+  void SetRMSFactorImage(aocommon::Image&& image) {
+    _rmsFactorImage = std::move(image);
+  }
+  const aocommon::Image& RMSFactorImage() const { return _rmsFactorImage; }
+
+ protected:
+  DeconvolutionAlgorithm();
+
+  DeconvolutionAlgorithm(const DeconvolutionAlgorithm&) = default;
+  DeconvolutionAlgorithm& operator=(const DeconvolutionAlgorithm&) = default;
+
+  void PerformSpectralFit(float* values, size_t x, size_t y) const;
+
+  float _threshold, _majorIterThreshold, _gain, _mGain, _cleanBorderRatio;
+  size_t _maxIter, _iterationNumber, _threadCount;
+  bool _allowNegativeComponents, _stopOnNegativeComponent;
+  const bool* _cleanMask;
+  aocommon::Image _rmsFactorImage;
+  mutable std::vector<float> _fittingScratch;
+
+  aocommon::LogReceiver* _logReceiver;
+
+  schaapcommon::fitters::SpectralFitter _spectralFitter;
+};
+}  // namespace radler::algorithms
+#endif  // RADLER_ALGORITHMS_DECONVOLUTION_ALGORITHM_H_
--- a/cpp/algorithms/generic_clean.cc
+++ b/cpp/algorithms/generic_clean.cc
+// Copyright (C) 2022 ASTRON (Netherlands Institute for Radio Astronomy)
+// SPDX-License-Identifier: GPL-3.0-or-later
+
+#include "algorithms/generic_clean.h"
+
+#include <aocommon/image.h>
+#include <aocommon/lane.h>
+#include <aocommon/units/fluxdensity.h>
+
+#include "algorithms/subminor_loop.h"
+#include "algorithms/threaded_deconvolution_tools.h"
+#include "math/peak_finder.h"
+
+using aocommon::units::FluxDensity;
+
+namespace radler::algorithms {
+namespace {
+std::string peakDescription(const aocommon::Image& image, size_t x, size_t y) {
+  std::ostringstream str;
+  const size_t index = x + y * image.Width();
+  const float peak = image[index];
+  str << FluxDensity::ToNiceString(peak) << " at " << x << "," << y;
+  return str.str();
+}
+}  // namespace
+
+GenericClean::GenericClean(bool useSubMinorOptimization)
+    : _convolutionPadding(1.1),
+      _useSubMinorOptimization(useSubMinorOptimization) {}
+
+float GenericClean::ExecuteMajorIteration(
+    ImageSet& dirtySet, ImageSet& modelSet,
+    const std::vector<aocommon::Image>& psfs, bool& reachedMajorThreshold) {
+  const size_t width = dirtySet.Width();
+  const size_t height = dirtySet.Height();
+  const size_t iterationCounterAtStart = _iterationNumber;
+  if (_stopOnNegativeComponent) _allowNegativeComponents = true;
+  _convolutionWidth = ceil(_convolutionPadding * width);
+  _convolutionHeight = ceil(_convolutionPadding * height);
+  if (_convolutionWidth % 2 != 0) ++_convolutionWidth;
+  if (_convolutionHeight % 2 != 0) ++_convolutionHeight;
+
+  aocommon::Image integrated(width, height);
+  aocommon::Image scratchA(_convolutionWidth, _convolutionHeight);
+  aocommon::Image scratchB(_convolutionWidth, _convolutionHeight);
+  dirtySet.GetLinearIntegrated(integrated);
+  size_t componentX = 0;
+  size_t componentY = 0;
+  std::optional<float> maxValue =
+      findPeak(integrated, scratchA.Data(), componentX, componentY);
+  if (!maxValue) {
+    _logReceiver->Info << "No peak found.\n";
+    reachedMajorThreshold = false;
+    return 0.0;
+  }
+  _logReceiver->Info << "Initial peak: "
+                     << peakDescription(integrated, componentX, componentY)
+                     << '\n';
+  float firstThreshold = this->_threshold;
+  float majorIterThreshold = std::max<float>(
+      MajorIterThreshold(), std::fabs(*maxValue) * (1.0 - this->_mGain));
+  if (majorIterThreshold > firstThreshold) {
+    firstThreshold = majorIterThreshold;
+    _logReceiver->Info << "Next major iteration at: "
+                       << FluxDensity::ToNiceString(majorIterThreshold) << '\n';
+  } else if (this->_mGain != 1.0) {
+    _logReceiver->Info
+        << "Major iteration threshold reached global threshold of "
+        << FluxDensity::ToNiceString(this->_threshold)
+        << ": final major iteration.\n";
+  }
+
+  if (_useSubMinorOptimization) {
+    size_t startIteration = _iterationNumber;
+    SubMinorLoop subMinorLoop(width, height, _convolutionWidth,
+                              _convolutionHeight, *_logReceiver);
+    subMinorLoop.SetIterationInfo(_iterationNumber, MaxNIter());
+    subMinorLoop.SetThreshold(firstThreshold, firstThreshold * 0.99);
+    subMinorLoop.SetGain(Gain());
+    subMinorLoop.SetAllowNegativeComponents(AllowNegativeComponents());
+    subMinorLoop.SetStopOnNegativeComponent(StopOnNegativeComponents());
+    subMinorLoop.SetSpectralFitter(&Fitter());
+    if (!_rmsFactorImage.Empty())
+      subMinorLoop.SetRMSFactorImage(_rmsFactorImage);
+    if (_cleanMask) subMinorLoop.SetMask(_cleanMask);
+    const size_t horBorderSize = std::round(width * CleanBorderRatio());
+    const size_t vertBorderSize = std::round(height * CleanBorderRatio());
+    subMinorLoop.SetCleanBorders(horBorderSize, vertBorderSize);
+    subMinorLoop.SetThreadCount(_threadCount);
+
+    maxValue = subMinorLoop.Run(dirtySet, psfs);
+
+    _iterationNumber = subMinorLoop.CurrentIteration();
+
+    _logReceiver->Info
+        << "Performed " << _iterationNumber << " iterations in total, "
+        << (_iterationNumber - startIteration)
+        << " in this major iteration with sub-minor optimization.\n";
+
+    for (size_t imageIndex = 0; imageIndex != dirtySet.size(); ++imageIndex) {
+      // TODO this can be multi-threaded if each thread has its own temporaries
+      const aocommon::Image& psf = psfs[dirtySet.PSFIndex(imageIndex)];
+      subMinorLoop.CorrectResidualDirty(scratchA.Data(), scratchB.Data(),
+                                        integrated.Data(), imageIndex,
+                                        dirtySet.Data(imageIndex), psf.Data());
+
+      subMinorLoop.GetFullIndividualModel(imageIndex, scratchA.Data());
+      float* model = modelSet.Data(imageIndex);
+      for (size_t i = 0; i != width * height; ++i)
+        model[i] += scratchA.Data()[i];
+    }
+  } else {
+    ThreadedDeconvolutionTools tools(_threadCount);
+    size_t peakIndex = componentX + componentY * width;
+
+    aocommon::UVector<float> peakValues(dirtySet.size());
+
+    while (maxValue && fabs(*maxValue) > firstThreshold &&
+           this->_iterationNumber < this->_maxIter &&
+           !(maxValue < 0.0f && this->_stopOnNegativeComponent)) {
+      if (this->_iterationNumber <= 10 ||
+          (this->_iterationNumber <= 100 && this->_iterationNumber % 10 == 0) ||
+          (this->_iterationNumber <= 1000 &&
+           this->_iterationNumber % 100 == 0) ||
+          this->_iterationNumber % 1000 == 0)
+        _logReceiver->Info << "Iteration " << this->_iterationNumber << ": "
+                           << peakDescription(integrated, componentX,
+                                              componentY)
+                           << '\n';
+
+      for (size_t i = 0; i != dirtySet.size(); ++i)
+        peakValues[i] = dirtySet[i][peakIndex];
+
+      PerformSpectralFit(peakValues.data(), componentX, componentY);
+
+      for (size_t i = 0; i != dirtySet.size(); ++i) {
+        peakValues[i] *= this->_gain;
+        modelSet.Data(i)[peakIndex] += peakValues[i];
+
+        size_t psfIndex = dirtySet.PSFIndex(i);
+
+        tools.SubtractImage(dirtySet.Data(i), psfs[psfIndex], componentX,
+                            componentY, peakValues[i]);
+      }
+
+      dirtySet.GetSquareIntegrated(integrated, scratchA);
+      maxValue = findPeak(integrated, scratchA.Data(), componentX, componentY);
+
+      peakIndex = componentX + componentY * width;
+
+      ++this->_iterationNumber;
+    }
+  }
+  if (maxValue) {
+    _logReceiver->Info << "Stopped on peak "
+                       << FluxDensity::ToNiceString(*maxValue) << ", because ";
+    bool maxIterReached = _iterationNumber >= MaxNIter(),
+         finalThresholdReached =
+             std::fabs(*maxValue) <= _threshold || maxValue == 0.0f,
+         negativeReached = maxValue < 0.0f && this->_stopOnNegativeComponent,
+         mgainReached = std::fabs(*maxValue) <= majorIterThreshold,
+         didWork = (_iterationNumber - iterationCounterAtStart) != 0;
+
+    if (maxIterReached)
+      _logReceiver->Info << "maximum number of iterations was reached.\n";
+    else if (finalThresholdReached)
+      _logReceiver->Info << "the threshold was reached.\n";
+    else if (negativeReached)
+      _logReceiver->Info << "a negative component was found.\n";
+    else if (!didWork)
+      _logReceiver->Info << "no iterations could be performed.\n";
+    else
+      _logReceiver->Info << "the minor-loop threshold was reached. Continuing "
+                            "cleaning after inversion/prediction round.\n";
+    reachedMajorThreshold =
+        mgainReached && didWork && !negativeReached && !finalThresholdReached;
+    return *maxValue;
+  } else {
+    _logReceiver->Info << "Deconvolution aborted.\n";
+    reachedMajorThreshold = false;
+    return 0.0;
+  }
+}
+
+std::optional<float> GenericClean::findPeak(const aocommon::Image& image,
+                                            float* scratch_buffer, size_t& x,
+                                            size_t& y) {
+  const float* actual_image = image.Data();
+  if (!_rmsFactorImage.Empty()) {
+    std::copy_n(image.Data(), image.Size(), scratch_buffer);
+    for (size_t i = 0; i != image.Size(); ++i) {
+      scratch_buffer[i] *= _rmsFactorImage[i];
+    }
+    actual_image = scratch_buffer;
+  }
+
+  if (_cleanMask == nullptr) {
+    return math::PeakFinder::Find(actual_image, image.Width(), image.Height(),
+                                  x, y, _allowNegativeComponents, 0,
+                                  image.Height(), _cleanBorderRatio);
+  } else {
+    return math::PeakFinder::FindWithMask(
+        actual_image, image.Width(), image.Height(), x, y,
+        _allowNegativeComponents, 0, image.Height(), _cleanMask,
+        _cleanBorderRatio);
+  }
+}
+}  // namespace radler::algorithms
\ No newline at end of file
--- a/cpp/algorithms/generic_clean.h
+++ b/cpp/algorithms/generic_clean.h
+// Copyright (C) 2022 ASTRON (Netherlands Institute for Radio Astronomy)
+// SPDX-License-Identifier: GPL-3.0-or-later
+
+#ifndef RADLER_GENERIC_CLEAN_H_
+#define RADLER_GENERIC_CLEAN_H_
+
+#include <optional>
+
+#include <aocommon/uvector.h>
+
+#include "image_set.h"
+#include "algorithms/deconvolution_algorithm.h"
+#include "algorithms/simple_clean.h"
+
+namespace radler::algorithms {
+/**
+ * This class implements a generalized version of Högbom clean. It performs a
+ * single-channel or joined cleaning, depending on the number of images
+ * provided. It can use a Clark-like optimization to speed up the cleaning. When
+ * multiple frequencies are provided, it can perform spectral fitting.
+ */
+class GenericClean : public DeconvolutionAlgorithm {
+ public:
+  explicit GenericClean(bool useSubMinorOptimization);
+
+  float ExecuteMajorIteration(ImageSet& dirtySet, ImageSet& modelSet,
+                              const std::vector<aocommon::Image>& psfs,
+                              bool& reachedMajorThreshold) final override;
+
+  virtual std::unique_ptr<DeconvolutionAlgorithm> Clone() const final override {
+    return std::unique_ptr<DeconvolutionAlgorithm>(new GenericClean(*this));
+  }
+
+ private:
+  size_t _convolutionWidth;
+  size_t _convolutionHeight;
+  const float _convolutionPadding;
+  bool _useSubMinorOptimization;
+
+  // Scratch buffer should at least accomodate space for image.Size() floats
+  // and is only used to avoid unnecessary memory allocations.
+  std::optional<float> findPeak(const aocommon::Image& image,
+                                float* scratch_buffer, size_t& x, size_t& y);
+};
+}  // namespace radler::algorithms
+#endif
--- a/cpp/algorithms/iuwt/image_analysis.cc
+++ b/cpp/algorithms/iuwt/image_analysis.cc
+// Copyright (C) 2022 ASTRON (Netherlands Institute for Radio Astronomy)
+// SPDX-License-Identifier: GPL-3.0-or-later
+
+#include "algorithms/iuwt/image_analysis.h"
+
+#include <stack>
+
+namespace radler::algorithms::iuwt {
+
+bool ImageAnalysis::IsHighestOnScale0(const IUWTDecomposition& iuwt,
+                                      IUWTMask& markedMask, size_t& x,
+                                      size_t& y, size_t endScale,
+                                      float& highestScale0) {
+  const size_t width = iuwt.Width(), height = iuwt.Height();
+  Component component(x, y, 0);
+  std::stack<Component> todo;
+  todo.push(component);
+  markedMask[0][x + y * width] = false;
+  highestScale0 = iuwt[0][x + y * width];
+  float highestOtherScales = 0.0;
+  while (!todo.empty()) {
+    Component c = todo.top();
+    todo.pop();
+    size_t index = c.x + c.y * width;
+    if (c.scale == 0) {
+      if (exceedsThreshold(iuwt[0][index], highestScale0)) {
+        highestScale0 = iuwt[0][index];
+        x = c.x;
+        y = c.y;
+      }
+    } else {
+      if (exceedsThreshold(iuwt[c.scale][index], highestOtherScales))
+        highestOtherScales = iuwt[c.scale][index];
+    }
+    if (c.x > 0) {
+      if (markedMask[c.scale][index - 1]) {
+        markedMask[c.scale][index - 1] = false;
+        todo.push(Component(c.x - 1, c.y, c.scale));
+      }
+    }
+    if (c.x < width - 1) {
+      if (markedMask[c.scale][index + 1]) {
+        markedMask[c.scale][index + 1] = false;
+        todo.push(Component(c.x + 1, c.y, c.scale));
+      }
+    }
+    if (c.y > 0) {
+      if (markedMask[c.scale][index - width]) {
+        markedMask[c.scale][index - width] = false;
+        todo.push(Component(c.x, c.y - 1, c.scale));
+      }
+    }
+    if (c.y < height - 1) {
+      if (markedMask[c.scale][index + width]) {
+        markedMask[c.scale][index + width] = false;
+        todo.push(Component(c.x, c.y + 1, c.scale));
+      }
+    }
+    if (c.scale > int(0)) {
+      if (markedMask[c.scale - 1][index]) {
+        markedMask[c.scale - 1][index] = false;
+        todo.push(Component(c.x, c.y, c.scale - 1));
+      }
+    }
+    if (c.scale < int(endScale) - 1) {
+      if (markedMask[c.scale + 1][index]) {
+        markedMask[c.scale + 1][index] = false;
+        todo.push(Component(c.x, c.y, c.scale + 1));
+      }
+    }
+  }
+  return std::fabs(highestScale0) > std::fabs(highestOtherScales);
+}
+
+void ImageAnalysis::Floodfill(const IUWTDecomposition& iuwt, IUWTMask& mask,
+                              const aocommon::UVector<float>& thresholds,
+                              size_t minScale, size_t endScale,
+                              const Component& component, float cleanBorder,
+                              size_t& areaSize) {
+  const size_t width = iuwt.Width(), height = iuwt.Height();
+  size_t xBorder = cleanBorder * width;
+  size_t yBorder = cleanBorder * height;
+  size_t minX = xBorder, maxX = width - xBorder;
+  size_t minY = yBorder, maxY = height - yBorder;
+
+  areaSize = 0;
+  endScale = std::min<size_t>(endScale, iuwt.NScales());
+  std::stack<Component> todo;
+  todo.push(component);
+  mask[component.scale][component.x + component.y * width] = true;
+  while (!todo.empty()) {
+    Component c = todo.top();
+    ++areaSize;
+    todo.pop();
+    size_t index = c.x + c.y * width;
+    if (c.x > minX) {
+      if (exceedsThreshold(iuwt[c.scale][index - 1], thresholds[c.scale]) &&
+          !mask[c.scale][index - 1]) {
+        mask[c.scale][index - 1] = true;
+        todo.push(Component(c.x - 1, c.y, c.scale));
+      }
+    }
+    if (c.x < maxX - 1) {
+      if (exceedsThreshold(iuwt[c.scale][index + 1], thresholds[c.scale]) &&
+          !mask[c.scale][index + 1]) {
+        mask[c.scale][index + 1] = true;
+        todo.push(Component(c.x + 1, c.y, c.scale));
+      }
+    }
+    if (c.y > minY) {
+      if (exceedsThreshold(iuwt[c.scale][index - width], thresholds[c.scale]) &&
+          !mask[c.scale][index - width]) {
+        mask[c.scale][index - width] = true;
+        todo.push(Component(c.x, c.y - 1, c.scale));
+      }
+    }
+    if (c.y < maxY - 1) {
+      if (exceedsThreshold(iuwt[c.scale][index + width], thresholds[c.scale]) &&
+          !mask[c.scale][index + width]) {
+        mask[c.scale][index + width] = true;
+        todo.push(Component(c.x, c.y + 1, c.scale));
+      }
+    }
+    if (c.scale > int(minScale)) {
+      if (exceedsThreshold(iuwt[c.scale - 1][index], thresholds[c.scale - 1]) &&
+          !mask[c.scale - 1][index]) {
+        mask[c.scale - 1][index] = true;
+        todo.push(Component(c.x, c.y, c.scale - 1));
+      }
+    }
+    if (c.scale < int(endScale) - 1) {
+      if (exceedsThreshold(iuwt[c.scale + 1][index], thresholds[c.scale + 1]) &&
+          !mask[c.scale + 1][index]) {
+        mask[c.scale + 1][index] = true;
+        todo.push(Component(c.x, c.y, c.scale + 1));
+      }
+    }
+  }
+}
+
+void ImageAnalysis::MaskedFloodfill(const IUWTDecomposition& iuwt,
+                                    IUWTMask& mask,
+                                    const aocommon::UVector<float>& thresholds,
+                                    size_t minScale, size_t endScale,
+                                    const Component& component,
+                                    float cleanBorder, const bool* priorMask,
+                                    size_t& areaSize) {
+  const size_t width = iuwt.Width(), height = iuwt.Height();
+  size_t xBorder = cleanBorder * width;
+  size_t yBorder = cleanBorder * height;
+  size_t minX = xBorder, maxX = width - xBorder;
+  size_t minY = yBorder, maxY = height - yBorder;
+
+  areaSize = 0;
+  endScale = std::min<size_t>(endScale, iuwt.NScales());
+  std::stack<Component> todo;
+  todo.push(component);
+  mask[component.scale][component.x + component.y * width] = true;
+  while (!todo.empty()) {
+    Component c = todo.top();
+    ++areaSize;
+    todo.pop();
+    size_t index = c.x + c.y * width;
+    if (c.x > minX) {
+      if (exceedsThreshold(iuwt[c.scale][index - 1], thresholds[c.scale]) &&
+          !mask[c.scale][index - 1] && priorMask[index - 1]) {
+        mask[c.scale][index - 1] = true;
+        todo.push(Component(c.x - 1, c.y, c.scale));
+      }
+    }
+    if (c.x < maxX - 1) {
+      if (exceedsThreshold(iuwt[c.scale][index + 1], thresholds[c.scale]) &&
+          !mask[c.scale][index + 1] && priorMask[index + 1]) {
+        mask[c.scale][index + 1] = true;
+        todo.push(Component(c.x + 1, c.y, c.scale));
+      }
+    }
+    if (c.y > minY) {
+      if (exceedsThreshold(iuwt[c.scale][index - width], thresholds[c.scale]) &&
+          !mask[c.scale][index - width] && priorMask[index - width]) {
+        mask[c.scale][index - width] = true;
+        todo.push(Component(c.x, c.y - 1, c.scale));
+      }
+    }
+    if (c.y < maxY - 1) {
+      if (exceedsThreshold(iuwt[c.scale][index + width], thresholds[c.scale]) &&
+          !mask[c.scale][index + width] && priorMask[index + width]) {
+        mask[c.scale][index + width] = true;
+        todo.push(Component(c.x, c.y + 1, c.scale));
+      }
+    }
+    if (c.scale > int(minScale)) {
+      if (exceedsThreshold(iuwt[c.scale - 1][index], thresholds[c.scale - 1]) &&
+          !mask[c.scale - 1][index] && priorMask[index]) {
+        mask[c.scale - 1][index] = true;
+        todo.push(Component(c.x, c.y, c.scale - 1));
+      }
+    }
+    if (c.scale < int(endScale) - 1) {
+      if (exceedsThreshold(iuwt[c.scale + 1][index], thresholds[c.scale + 1]) &&
+          !mask[c.scale + 1][index] && priorMask[index]) {
+        mask[c.scale + 1][index] = true;
+        todo.push(Component(c.x, c.y, c.scale + 1));
+      }
+    }
+  }
+}
+
+void ImageAnalysis::SelectStructures(const IUWTDecomposition& iuwt,
+                                     IUWTMask& mask,
+                                     const aocommon::UVector<float>& thresholds,
+                                     size_t minScale, size_t endScale,
+                                     float cleanBorder, const bool* priorMask,
+                                     size_t& areaSize) {
+  const size_t width = iuwt.Width(), height = iuwt.Height();
+  const size_t xBorder = cleanBorder * width, yBorder = cleanBorder * height,
+               minX = xBorder, maxX = width - xBorder, minY = yBorder,
+               maxY = height - yBorder;
+
+  areaSize = 0;
+
+  for (size_t scale = minScale; scale != endScale; ++scale) {
+    for (size_t y = minY; y != maxY; ++y) {
+      for (size_t x = minX; x != maxX; ++x) {
+        size_t index = x + y * width;
+        bool isInPriorMask = (priorMask == nullptr) || priorMask[index];
+        if (exceedsThreshold(iuwt[scale][index], thresholds[scale]) &&
+            !mask[scale][index] && isInPriorMask) {
+          Component component(x, y, scale);
+          size_t subAreaSize = 0;
+          if (priorMask == nullptr)
+            Floodfill(iuwt, mask, thresholds, minScale, endScale, component,
+                      cleanBorder, subAreaSize);
+          else
+            MaskedFloodfill(iuwt, mask, thresholds, minScale, endScale,
+                            component, cleanBorder, priorMask, subAreaSize);
+          areaSize += subAreaSize;
+        }
+      }
+    }
+  }
+}
+
+void ImageAnalysis::FloodFill2D(const float* image, bool* mask, float threshold,
+                                const ImageAnalysis::Component2D& component,
+                                size_t width, size_t height, size_t& areaSize) {
+  areaSize = 0;
+  std::stack<Component2D> todo;
+  todo.push(component);
+  mask[component.x + component.y * width] = true;
+  while (!todo.empty()) {
+    Component2D c = todo.top();
+    ++areaSize;
+    todo.pop();
+    size_t index = c.x + c.y * width;
+    if (c.x > 0) {
+      if (exceedsThreshold(image[index - 1], threshold) && !mask[index - 1]) {
+        mask[index - 1] = true;
+        todo.push(Component2D(c.x - 1, c.y));
+      }
+    }
+    if (c.x < width - 1) {
+      if (exceedsThreshold(image[index + 1], threshold) && !mask[index + 1]) {
+        mask[index + 1] = true;
+        todo.push(Component2D(c.x + 1, c.y));
+      }
+    }
+    if (c.y > 0) {
+      if (exceedsThreshold(image[index - width], threshold) &&
+          !mask[index - width]) {
+        mask[index - width] = true;
+        todo.push(Component2D(c.x, c.y - 1));
+      }
+    }
+    if (c.y < height - 1) {
+      if (exceedsThreshold(image[index + width], threshold) &&
+          !mask[index + width]) {
+        mask[index + width] = true;
+        todo.push(Component2D(c.x, c.y + 1));
+      }
+    }
+  }
+}
+
+void ImageAnalysis::FloodFill2D(const float* image, bool* mask, float threshold,
+                                const ImageAnalysis::Component2D& component,
+                                size_t width, size_t height,
+                                std::vector<Component2D>& area) {
+  area.clear();
+  std::stack<Component2D> todo;
+  todo.push(component);
+  mask[component.x + component.y * width] = true;
+  while (!todo.empty()) {
+    Component2D c = todo.top();
+    area.push_back(c);
+    todo.pop();
+    size_t index = c.x + c.y * width;
+    if (c.x > 0) {
+      if (exceedsThresholdAbs(image[index - 1], threshold) &&
+          !mask[index - 1]) {
+        mask[index - 1] = true;
+        todo.push(Component2D(c.x - 1, c.y));
+      }
+    }
+    if (c.x < width - 1) {
+      if (exceedsThresholdAbs(image[index + 1], threshold) &&
+          !mask[index + 1]) {
+        mask[index + 1] = true;
+        todo.push(Component2D(c.x + 1, c.y));
+      }
+    }
+    if (c.y > 0) {
+      if (exceedsThresholdAbs(image[index - width], threshold) &&
+          !mask[index - width]) {
+        mask[index - width] = true;
+        todo.push(Component2D(c.x, c.y - 1));
+      }
+    }
+    if (c.y < height - 1) {
+      if (exceedsThresholdAbs(image[index + width], threshold) &&
+          !mask[index + width]) {
+        mask[index + width] = true;
+        todo.push(Component2D(c.x, c.y + 1));
+      }
+    }
+  }
+}
+}  // namespace radler::algorithms::iuwt
\ No newline at end of file
--- a/cpp/algorithms/iuwt/image_analysis.h
+++ b/cpp/algorithms/iuwt/image_analysis.h
+// Copyright (C) 2022 ASTRON (Netherlands Institute for Radio Astronomy)
+// SPDX-License-Identifier: GPL-3.0-or-later
+
+#ifndef RADLER_ALGORITHMS_IUWT_IMAGE_ANALYSIS_H_
+#define RADLER_ALGORITHMS_IUWT_IMAGE_ANALYSIS_H_
+
+#include "algorithms/iuwt/iuwt_decomposition.h"
+
+namespace radler::algorithms::iuwt {
+
+class ImageAnalysis {
+ public:
+  struct Component {
+    Component(size_t _x, size_t _y, int _scale) : x(_x), y(_y), scale(_scale) {}
+
+    std::string ToString() const {
+      std::ostringstream str;
+      str << x << ',' << y << ", scale " << scale;
+      return str.str();
+    }
+
+    size_t x, y;
+    int scale;
+  };
+
+  struct Component2D {
+    Component2D(size_t _x, size_t _y) : x(_x), y(_y) {}
+
+    std::string ToString() const {
+      std::ostringstream str;
+      str << x << ',' << y;
+      return str.str();
+    }
+
+    size_t x, y;
+  };
+
+  static void SelectStructures(const IUWTDecomposition& iuwt, IUWTMask& mask,
+                               const aocommon::UVector<float>& thresholds,
+                               size_t minScale, size_t endScale,
+                               float cleanBorder, const bool* priorMask,
+                               size_t& areaSize);
+
+  static bool IsHighestOnScale0(const IUWTDecomposition& iuwt,
+                                IUWTMask& markedMask, size_t& x, size_t& y,
+                                size_t endScale, float& highestScale0);
+
+  static void Floodfill(const IUWTDecomposition& iuwt, IUWTMask& mask,
+                        const aocommon::UVector<float>& thresholds,
+                        size_t minScale, size_t endScale,
+                        const Component& component, float cleanBorder,
+                        size_t& areaSize);
+
+  static void MaskedFloodfill(const IUWTDecomposition& iuwt, IUWTMask& mask,
+                              const aocommon::UVector<float>& thresholds,
+                              size_t minScale, size_t endScale,
+                              const Component& component, float cleanBorder,
+                              const bool* priorMask, size_t& areaSize);
+
+  static void FloodFill2D(const float* image, bool* mask, float threshold,
+                          const Component2D& component, size_t width,
+                          size_t height, size_t& areaSize);
+
+  /**
+   * Exactly like above, but now collecting the components in the
+   * area vector, instead of returning just the area size.
+   */
+  static void FloodFill2D(const float* image, bool* mask, float threshold,
+                          const Component2D& component, size_t width,
+                          size_t height, std::vector<Component2D>& area);
+
+ private:
+  static bool exceedsThreshold(float val, float threshold) {
+    if (threshold >= 0.0)
+      return val > threshold;
+    else
+      return val < threshold || val > -threshold;
+  }
+  static bool exceedsThresholdAbs(float val, float threshold) {
+    return std::fabs(val) > threshold;
+  }
+};
+}  // namespace radler::algorithms::iuwt
+#endif  // RADLER_ALGORITHMS_IUWT_IMAGE_ANALYSIS_H_
--- a/cpp/algorithms/iuwt/iuwt_decomposition.cc
+++ b/cpp/algorithms/iuwt/iuwt_decomposition.cc
+// Copyright (C) 2022 ASTRON (Netherlands Institute for Radio Astronomy)
+// SPDX-License-Identifier: GPL-3.0-or-later
+
+#include "algorithms/iuwt/iuwt_decomposition.h"
+
+using aocommon::Image;
+
+namespace radler::algorithms::iuwt {
+
+void IUWTDecomposition::DecomposeMT(aocommon::StaticFor<size_t>& loop,
+                                    const float* input, float* scratch,
+                                    bool includeLargest) {
+  Image& i1(_scales.back().Coefficients());
+  i1 = Image(_width, _height);
+
+  // The first iteration of the loop, unrolled, so that we don't have to
+  // copy the input into i0.
+  Image& coefficients0 = _scales[0].Coefficients();
+  coefficients0 = Image(_width, _height);
+  convolveMT(loop, i1.Data(), input, scratch, _width, _height, 1);
+  convolveMT(loop, coefficients0.Data(), i1.Data(), scratch, _width, _height,
+             1);
+
+  // coefficients = i0 - i2
+  differenceMT(loop, coefficients0.Data(), input, coefficients0.Data(), _width,
+               _height);
+
+  // i0 = i1;
+  Image i0(i1);
+
+  for (int scale = 1; scale != int(_scaleCount); ++scale) {
+    Image& coefficients = _scales[scale].Coefficients();
+    coefficients = Image(_width, _height);
+    convolveMT(loop, i1.Data(), i0.Data(), scratch, _width, _height, scale + 1);
+    convolveMT(loop, coefficients.Data(), i1.Data(), scratch, _width, _height,
+               scale + 1);
+
+    // coefficients = i0 - i2
+    differenceMT(loop, coefficients.Data(), i0.Data(), coefficients.Data(),
+                 _width, _height);
+
+    // i0 = i1;
+    if (scale + 1 != int(_scaleCount)) {
+      std::copy_n(i1.Data(), _width * _height, i0.Data());
+    }
+  }
+
+  // The largest (residual) scales are in i1, but since the
+  // largest scale is aliased to i1, it's already stored there.
+  // Hence we can skip this:
+  // if(includeLargest)
+  //	_scales.back().Coefficients() = i1;
+
+  // Do free the memory of the largest scale if it is not necessary:
+  if (!includeLargest) _scales.back().Coefficients().Reset();
+}
+
+void IUWTDecomposition::convolveMT(aocommon::StaticFor<size_t>& loop,
+                                   float* output, const float* image,
+                                   float* scratch, size_t width, size_t height,
+                                   int scale) {
+  loop.Run(0, height, [&](size_t y_start, size_t y_end) {
+    convolveHorizontalPartial(scratch, image, width, y_start, y_end, scale);
+  });
+
+  loop.Run(0, width, [&](size_t x_start, size_t x_end) {
+    convolveVerticalPartialFast(output, scratch, width, height, x_start, x_end,
+                                scale);
+  });
+}
+
+void IUWTDecomposition::differenceMT(aocommon::StaticFor<size_t>& loop,
+                                     float* dest, const float* lhs,
+                                     const float* rhs, size_t width,
+                                     size_t height) {
+  loop.Run(0, height, [&](size_t y_start, size_t y_end) {
+    differencePartial(dest, lhs, rhs, width, y_start, y_end);
+  });
+}
+
+void IUWTDecomposition::convolveHorizontalFast(float* output,
+                                               const float* image, size_t width,
+                                               size_t height, int scale) {
+  const size_t H_SIZE = 5;
+  const float h[H_SIZE] = {1.0 / 16.0, 4.0 / 16.0, 6.0 / 16.0, 4.0 / 16.0,
+                           1.0 / 16.0};
+  int scaleDist = (1 << scale);
+  int dist[H_SIZE];
+  size_t minX[H_SIZE], maxX[H_SIZE];
+  for (int hIndex = 0; hIndex != H_SIZE; ++hIndex) {
+    int hShift = hIndex - H_SIZE / 2;
+    dist[hIndex] = (scaleDist - 1) * hShift;
+    minX[hIndex] = std::max<int>(0, -dist[hIndex]);
+    maxX[hIndex] = std::min<int>(width, width - dist[hIndex]);
+  }
+  for (size_t y = 0; y != height; ++y) {
+    float* outputPtr = &output[y * width];
+    const float* inputPtr = &image[y * width];
+
+    for (size_t x = 0; x != minX[1]; ++x) {
+      outputPtr[x] = inputPtr[x + dist[2]] * h[2] +
+                     inputPtr[x + dist[3]] * h[3] +
+                     inputPtr[x + dist[4]] * h[4];
+    }
+
+    for (size_t x = minX[1]; x != minX[0]; ++x) {
+      outputPtr[x] =
+          inputPtr[x + dist[2]] * h[2] + inputPtr[x + dist[1]] * h[1] +
+          inputPtr[x + dist[3]] * h[3] + inputPtr[x + dist[4]] * h[4];
+    }
+
+    for (size_t x = minX[0]; x != maxX[4]; ++x) {
+      outputPtr[x] =
+          inputPtr[x + dist[2]] * h[2] + inputPtr[x + dist[1]] * h[1] +
+          inputPtr[x + dist[0]] * h[0] + inputPtr[x + dist[3]] * h[3] +
+          inputPtr[x + dist[4]] * h[4];
+    }
+
+    for (size_t x = maxX[4]; x != maxX[3]; ++x) {
+      outputPtr[x] =
+          inputPtr[x + dist[2]] * h[2] + inputPtr[x + dist[1]] * h[1] +
+          inputPtr[x + dist[0]] * h[0] + inputPtr[x + dist[3]] * h[3];
+    }
+
+    for (size_t x = maxX[3]; x != width; ++x) {
+      outputPtr[x] = inputPtr[x + dist[2]] * h[2] +
+                     inputPtr[x + dist[1]] * h[1] +
+                     inputPtr[x + dist[0]] * h[0];
+    }
+  }
+}
+
+// This version is not as fast as the one below.
+void IUWTDecomposition::convolveVerticalPartialFastFailed(
+    float* output, const float* image, size_t width, size_t height,
+    size_t startX, size_t endX, int scale) {
+  const size_t H_SIZE = 5;
+  const float h[H_SIZE] = {1.0 / 16.0, 4.0 / 16.0, 6.0 / 16.0, 4.0 / 16.0,
+                           1.0 / 16.0};
+  int scaleDist = (1 << scale);
+  for (size_t y = 0; y < height; ++y) {
+    float* outputPtr = &output[y * width];
+    const float* inputPtr = &image[y * width];
+    for (size_t x = startX; x != endX; ++x) outputPtr[x] = inputPtr[x] * h[2];
+  }
+  for (int hIndex = 0; hIndex != H_SIZE; ++hIndex) {
+    if (hIndex != 2) {
+      int hShift = hIndex - H_SIZE / 2;
+      int dist = (scaleDist - 1) * hShift;
+      size_t minY = std::max<int>(0, -dist),
+             maxY = std::min<int>(height, height - dist);
+      for (size_t y = minY; y < maxY; ++y) {
+        float* outputPtr = &output[y * width];
+        const float* inputPtr = &image[(y + dist) * width];
+        for (size_t x = startX; x != endX; ++x)
+          outputPtr[x] += inputPtr[x] * h[hIndex];
+      }
+    }
+  }
+}
+
+void IUWTDecomposition::convolveVerticalPartialFast(float* output,
+                                                    const float* image,
+                                                    size_t width, size_t height,
+                                                    size_t startX, size_t endX,
+                                                    int scale) {
+  const size_t H_SIZE = 5;
+  const float h[H_SIZE] = {1.0 / 16.0, 4.0 / 16.0, 6.0 / 16.0, 4.0 / 16.0,
+                           1.0 / 16.0};
+  int scaleDist = (1 << scale);
+  int dist[H_SIZE];
+  size_t minY[H_SIZE], maxY[H_SIZE];
+  for (int hIndex = 0; hIndex != H_SIZE; ++hIndex) {
+    int hShift = hIndex - H_SIZE / 2;
+    dist[hIndex] = (scaleDist - 1) * hShift;
+    minY[hIndex] = std::max<int>(0, -dist[hIndex]);
+    maxY[hIndex] = std::min<int>(height, height - dist[hIndex]);
+  }
+
+  for (size_t y = 0; y != minY[1]; ++y) {
+    float* outputPtr = &output[y * width];
+    const float* inputPtr2 = &image[(y + dist[2]) * width];
+    const float* inputPtr3 = &image[(y + dist[3]) * width];
+    const float* inputPtr4 = &image[(y + dist[4]) * width];
+    for (size_t x = startX; x != endX; ++x) {
+      outputPtr[x] =
+          inputPtr2[x] * h[2] + inputPtr3[x] * h[3] + inputPtr4[x] * h[4];
+    }
+  }
+
+  for (size_t y = minY[1]; y != minY[0]; ++y) {
+    float* outputPtr = &output[y * width];
+    const float* inputPtr1 = &image[(y + dist[1]) * width];
+    const float* inputPtr2 = &image[(y + dist[2]) * width];
+    const float* inputPtr3 = &image[(y + dist[3]) * width];
+    const float* inputPtr4 = &image[(y + dist[4]) * width];
+    for (size_t x = startX; x != endX; ++x) {
+      outputPtr[x] = inputPtr1[x] * h[1] + inputPtr2[x] * h[2] +
+                     inputPtr3[x] * h[3] + inputPtr4[x] * h[4];
+    }
+  }
+
+  for (size_t y = minY[0]; y != maxY[4]; ++y) {
+    float* outputPtr = &output[y * width];
+    const float* inputPtr0 = &image[(y + dist[0]) * width];
+    const float* inputPtr1 = &image[(y + dist[1]) * width];
+    const float* inputPtr2 = &image[(y + dist[2]) * width];
+    const float* inputPtr3 = &image[(y + dist[3]) * width];
+    const float* inputPtr4 = &image[(y + dist[4]) * width];
+    for (size_t x = startX; x != endX; ++x) {
+      outputPtr[x] = inputPtr0[x] * h[0] + inputPtr1[x] * h[1] +
+                     inputPtr2[x] * h[2] + inputPtr3[x] * h[3] +
+                     inputPtr4[x] * h[4];
+    }
+  }
+
+  for (size_t y = maxY[4]; y != maxY[3]; ++y) {
+    float* outputPtr = &output[y * width];
+    const float* inputPtr0 = &image[(y + dist[0]) * width];
+    const float* inputPtr1 = &image[(y + dist[1]) * width];
+    const float* inputPtr2 = &image[(y + dist[2]) * width];
+    const float* inputPtr3 = &image[(y + dist[3]) * width];
+    for (size_t x = startX; x != endX; ++x) {
+      outputPtr[x] = inputPtr0[x] * h[0] + inputPtr1[x] * h[1] +
+                     inputPtr2[x] * h[2] + inputPtr3[x] * h[3];
+    }
+  }
+
+  for (size_t y = maxY[3]; y != height; ++y) {
+    float* outputPtr = &output[y * width];
+    const float* inputPtr0 = &image[(y + dist[0]) * width];
+    const float* inputPtr1 = &image[(y + dist[1]) * width];
+    const float* inputPtr2 = &image[(y + dist[2]) * width];
+    for (size_t x = startX; x != endX; ++x) {
+      outputPtr[x] =
+          inputPtr0[x] * h[0] + inputPtr1[x] * h[1] + inputPtr2[x] * h[2];
+    }
+  }
+}
+}  // namespace radler::algorithms::iuwt
\ No newline at end of file
No results found