Commit d8ef5aef authored by Andi Braimllari's avatar Andi Braimllari Committed by Tobias Lasser
Browse files

add an operator for the shearlet transform

parent 79b8f158
Pipeline #819190 passed with stages
in 12 minutes and 28 seconds
......@@ -202,7 +202,7 @@ with section("markup"):
with section("lint"):
# a list of lint codes to disable
disabled_codes = ["C0103", "C0111", "C0304", "C0306", "E1120", "R0912", "C0113", "C0301", "C0303", "C0305", "C0307", "W0105"]
disabled_codes = ["C0103", "C0111", "C0304", "C0306", "E1120", "R0912", "R0915", "C0113", "C0301", "C0303", "C0305", "C0307", "W0105"]
# regular expression pattern describing valid function names
function_pattern = '[0-9a-z_]+'
......
Shearlets
---------
Shearlets are a multiscale framework which provide an optimally sparse approximations of multivariate data governed by
anisotropic features [1]. They are constructed by applying three operations, translation, parabolic scaling, and
shearing, to mother functions. Here, we offer band-limited cone-adapted [2] shearlets, meaning that we have compact
support in the Fourier domain. For further reading regarding the exact functions used, please refer to the FFST
paper [3].
### Decompose and reconstruct an image using shearlets
We will now briefly explain the `ShearletTransform` class and go through an example.
We start by generating a modified 2D Shepp-Logan phantom and create a `ShearletTransform` object. Afterwards, we apply
the direct shearlet transform to the image, which generates a `DataContainer<real_t>` of shape (256, 256, 61). The last
dimension corresponds to the oversampling factor of the operation. Here, we generate 61 layers. Finally, we apply the
inverse shearlet transform to retrieve the reconstruction of the phantom.
```c++
// generate 2D phantom
IndexVector_t size(2);
size << 256, 256;
auto phantom = PhantomGenerator<real_t>::createModifiedSheppLogan(size);
ShearletTransform<real_t, real_t> shearletTransform(size);
Logger::get("Info")->info("Applying shearlet transform");
DataContainer<real_t> shearletCoefficients = shearletTransform.apply(phantom);
Logger::get("Info")->info("Applying inverse shearlet transform");
DataContainer<real_t> reconstruction = shearletTransform.applyAdjoint(shearletCoefficients);
// write the reconstruction out
EDF::write(reconstruction, "2dreconstruction_shearlet.edf");
```
Note that `apply` works only on two-dimensional objects, e.g. on a grayscale image.
By default, the `apply` method will return real values, and the imaginary parts cut out. If we want to the complex
numbers as-generated, then we can simply specify `ShearletTransform<std::complex<real_t>, real_t>`, and go from there.
As a subclass of the `LinearOperator` class, the `ShearletTransform` does not eagerly perform computations. Therefore,
the spectra are stored right after they've been generated on the first run of the `apply`/`applyAdjoint` calls. This
avoids redundant computation.
Given that the spectra are only related to the shape of the image and number of scales, one can reuse such an object
depending on the context.
After applying the direct and inverse shearlet transform, the image reconstruction is generated. From here we can do a
side-by-side comparison of the original phantom and its reconstruction,
![Modified Shepp-Logan phantom](./images/shearlets_2dphantom.png)
![Modified Shepp-Logan phantom shearlet reconstruction](./images/shearlets_2dreconstruction.png)
### References
[1] Gitta Kutyniok and Demetrio Labate, *Introduction to Shearlets*, https://www.math.uh.edu/~dlabate/SHBookIntro.pdf.
[2] Kanghui Guo, Gitta Kutyniok, and Demetrio Labate, *Sparse Multidimensional Representations using Anisotropic
Dilation and Shear Operators*, ISBN 0-0-9728482-x-x, https://www.math.uh.edu/~dlabate/Athens.pdf.
[3] Sören Häuser, Gabriele Steidl, *Fast Finite Shearlet Transform: a tutorial*, https://arxiv.org/pdf/1202.1773.pdf.
......@@ -28,3 +28,8 @@ Dictionary
===================
.. doxygenclass:: elsa::Dictionary
ShearletTransform
=================
.. doxygenclass:: elsa::ShearletTransform
......@@ -6,6 +6,7 @@ set(MODULE_HEADERS
Utilities/Badge.hpp
Utilities/DataContainerFormatter.hpp
Utilities/FormatConfig.h
Utilities/Math.hpp
Utilities/Statistics.hpp
Utilities/TypeCasts.hpp
Descriptors/DataDescriptor.h
......
#pragma once
#include "elsaDefines.h"
namespace elsa
{
/// proposed in Y. Meyer, Oscillating Patterns in Image Processing and Nonlinear Evolution
/// Equations. AMS, 2001
template <typename data_t>
data_t meyerFunction(data_t x)
{
if (x < 0) {
return 0;
} else if (0 <= x && x <= 1) {
return 35 * std::pow(x, 4) - 84 * std::pow(x, 5) + 70 * std::pow(x, 6)
- 20 * std::pow(x, 7);
} else {
return 1;
}
}
namespace shearlet
{
/// defined in Sören Häuser and Gabriele Steidl, Fast Finite Shearlet Transform: a
/// tutorial, 2014
template <typename data_t>
data_t b(data_t w)
{
if (1 <= std::abs(w) && std::abs(w) <= 2) {
return std::sin(pi<data_t> / 2.0 * meyerFunction(std::abs(w) - 1));
} else if (2 < std::abs(w) && std::abs(w) <= 4) {
return std::cos(pi<data_t> / 2.0 * meyerFunction(1.0 / 2 * std::abs(w) - 1));
} else {
return 0;
}
}
/// defined in Sören Häuser and Gabriele Steidl, Fast Finite Shearlet Transform: a
/// tutorial, 2014
template <typename data_t>
data_t phi(data_t w)
{
if (std::abs(w) <= 1.0 / 2) {
return 1;
} else if (1.0 / 2 < std::abs(w) && std::abs(w) < 1) {
return std::cos(pi<data_t> / 2.0 * meyerFunction(2 * std::abs(w) - 1));
} else {
return 0;
}
}
/// defined in Sören Häuser and Gabriele Steidl, Fast Finite Shearlet Transform: a
/// tutorial, 2014
template <typename data_t>
data_t phiHat(data_t w, data_t h)
{
if (std::abs(h) <= std::abs(w)) {
return phi(w);
} else {
return phi(h);
}
}
/// defined in Sören Häuser and Gabriele Steidl, Fast Finite Shearlet Transform: a
/// tutorial, 2014
template <typename data_t>
data_t psiHat1(data_t w)
{
return std::sqrt(std::pow(b(2 * w), 2) + std::pow(b(w), 2));
}
/// defined in Sören Häuser and Gabriele Steidl, Fast Finite Shearlet Transform: a
/// tutorial, 2014
template <typename data_t>
data_t psiHat2(data_t w)
{
if (w <= 0) {
return std::sqrt(meyerFunction(1 + w));
} else {
return std::sqrt(meyerFunction(1 - w));
}
}
/// defined in Sören Häuser and Gabriele Steidl, Fast Finite Shearlet Transform: a
/// tutorial, 2014
template <typename data_t>
data_t psiHat(data_t w, data_t h)
{
if (w == 0) {
return 0;
} else {
return psiHat1(w) * psiHat2(h / w);
}
}
} // namespace shearlet
} // namespace elsa
......@@ -47,6 +47,7 @@
#include "Scaling.h"
#include "FiniteDifferences.h"
#include "BlockLinearOperator.h"
#include "ShearletTransform.h"
// Proximity Operator headers
#include "ProximityOperator.h"
......
......@@ -5,7 +5,8 @@ set(MODULE_HEADERS
FiniteDifferences.h
FourierTransform.h
BlockLinearOperator.h
Dictionary.h)
Dictionary.h
ShearletTransform.h)
# list all the code files of the module
set(MODULE_SOURCES
......@@ -14,7 +15,8 @@ set(MODULE_SOURCES
FiniteDifferences.cpp
FourierTransform.cpp
BlockLinearOperator.cpp
Dictionary.cpp)
Dictionary.cpp
ShearletTransform.cpp)
list(APPEND MODULE_PUBLIC_DEPS elsa_core elsa_logging)
list(APPEND MODULE_PRIVATE_DEPS)
......
#include "ShearletTransform.h"
#include "FourierTransform.h"
#include "VolumeDescriptor.h"
#include "Timer.h"
#include "Math.hpp"
namespace elsa
{
template <typename ret_t, typename data_t>
ShearletTransform<ret_t, data_t>::ShearletTransform(IndexVector_t spatialDimensions)
: ShearletTransform(spatialDimensions[0], spatialDimensions[1])
{
if (spatialDimensions.size() != 2) {
throw LogicError("ShearletTransform: Only 2D shape supported");
}
}
template <typename ret_t, typename data_t>
ShearletTransform<ret_t, data_t>::ShearletTransform(index_t width, index_t height)
: ShearletTransform(width, height, calculateNumOfScales(width, height))
{
}
template <typename ret_t, typename data_t>
ShearletTransform<ret_t, data_t>::ShearletTransform(index_t width, index_t height,
index_t numOfScales)
: ShearletTransform(width, height, numOfScales, std::nullopt)
{
}
template <typename ret_t, typename data_t>
ShearletTransform<ret_t, data_t>::ShearletTransform(
index_t width, index_t height, index_t numOfScales,
std::optional<DataContainer<data_t>> spectra)
: LinearOperator<ret_t>(
VolumeDescriptor{{width, height}},
VolumeDescriptor{{width, height, calculateNumOfLayers(numOfScales)}}),
_spectra{spectra},
_width{width},
_height{height},
_numOfScales{numOfScales},
_numOfLayers{calculateNumOfLayers(numOfScales)}
{
if (width < 0 || height < 0) {
throw LogicError("ShearletTransform: negative width/height were provided");
}
if (numOfScales < 0) {
throw LogicError("ShearletTransform: negative number of scales was provided");
}
}
// TODO implement sumByAxis in DataContainer and remove me
template <typename ret_t, typename data_t>
DataContainer<elsa::complex<data_t>> ShearletTransform<ret_t, data_t>::sumByLastAxis(
DataContainer<elsa::complex<data_t>> dc) const
{
auto coeffsPerDim = dc.getDataDescriptor().getNumberOfCoefficientsPerDimension();
index_t width = coeffsPerDim[0];
index_t height = coeffsPerDim[1];
index_t layers = coeffsPerDim[2];
DataContainer<elsa::complex<data_t>> summedDC(VolumeDescriptor{{width, height}});
for (index_t j = 0; j < width; j++) {
for (index_t k = 0; k < height; k++) {
elsa::complex<data_t> currValue = 0;
for (index_t i = 0; i < layers; i++) {
currValue += dc(j, k, i);
}
summedDC(j, k) = currValue;
}
}
return summedDC;
}
template <typename ret_t, typename data_t>
void ShearletTransform<ret_t, data_t>::applyImpl(const DataContainer<ret_t>& x,
DataContainer<ret_t>& Ax) const
{
Timer timeguard("ShearletTransform", "apply");
if (_width != this->getDomainDescriptor().getNumberOfCoefficientsPerDimension()[0]
|| _height != this->getDomainDescriptor().getNumberOfCoefficientsPerDimension()[1]) {
throw InvalidArgumentError("ShearletTransform: Width and height of the input do not "
"match to that of this shearlet system");
}
Logger::get("ShearletTransform")
->info("Running the shearlet transform on a 2D signal of shape ({}, {}), on {} "
"scales with an oversampling factor of {} and {} spectra",
_width, _height, _numOfScales, _numOfLayers,
isSpectraComputed() ? "precomputed" : "non-precomputed");
if (!isSpectraComputed()) {
computeSpectra();
}
FourierTransform<elsa::complex<data_t>> fourierTransform(x.getDataDescriptor());
DataContainer<elsa::complex<data_t>> fftImg = fourierTransform.apply(x.asComplex());
for (index_t i = 0; i < getNumOfLayers(); i++) {
DataContainer<elsa::complex<data_t>> temp =
getSpectra().slice(i).viewAs(x.getDataDescriptor()).asComplex() * fftImg;
if constexpr (isComplex<ret_t>) {
Ax.slice(i) = fourierTransform.applyAdjoint(temp);
} else {
Ax.slice(i) = real(fourierTransform.applyAdjoint(temp));
}
}
}
template <typename ret_t, typename data_t>
void ShearletTransform<ret_t, data_t>::applyAdjointImpl(const DataContainer<ret_t>& y,
DataContainer<ret_t>& Aty) const
{
Timer timeguard("ShearletTransform", "applyAdjoint");
if (_width != this->getDomainDescriptor().getNumberOfCoefficientsPerDimension()[0]
|| _height != this->getDomainDescriptor().getNumberOfCoefficientsPerDimension()[1]) {
throw InvalidArgumentError("ShearletTransform: Width and height of the input do not "
"match to that of this shearlet system");
}
Logger::get("ShearletTransform")
->info("Running the inverse shearlet transform on a 2D signal of shape ({}, {}), "
"on {} "
"scales with an oversampling factor of {} and {} spectra",
_width, _height, _numOfScales, _numOfLayers,
isSpectraComputed() ? "precomputed" : "non-precomputed");
if (!isSpectraComputed()) {
computeSpectra();
}
FourierTransform<elsa::complex<data_t>> fourierTransform(Aty.getDataDescriptor());
DataContainer<elsa::complex<data_t>> intermRes(y.getDataDescriptor());
for (index_t i = 0; i < getNumOfLayers(); i++) {
DataContainer<elsa::complex<data_t>> temp =
fourierTransform.apply(y.slice(i).viewAs(Aty.getDataDescriptor()).asComplex())
* getSpectra().slice(i).viewAs(Aty.getDataDescriptor()).asComplex();
intermRes.slice(i) = fourierTransform.applyAdjoint(temp);
}
if constexpr (isComplex<ret_t>) {
Aty = sumByLastAxis(intermRes);
} else {
Aty = real(sumByLastAxis(intermRes));
}
}
template <typename ret_t, typename data_t>
void ShearletTransform<ret_t, data_t>::computeSpectra() const
{
if (isSpectraComputed()) {
Logger::get("ShearletTransform")->warn("Spectra have already been computed!");
}
_spectra = DataContainer<data_t>(VolumeDescriptor{{_width, _height, _numOfLayers}});
index_t i = 0;
_computeSpectraAtLowFreq(i);
for (index_t j = 0; j < _numOfScales; j++) {
auto twoPowJ = static_cast<index_t>(std::pow(2, j));
_computeSpectraAtSeamLines(i, j, -twoPowJ);
for (auto k = -twoPowJ + 1; k < twoPowJ; k++) {
_computeSpectraAtConicRegions(i, j, k);
}
_computeSpectraAtSeamLines(i, j, twoPowJ);
}
assert(i == _numOfLayers
&& "ShearletTransform: The layers of the spectra were indexed wrong");
}
template <typename ret_t, typename data_t>
void ShearletTransform<ret_t, data_t>::_computeSpectraAtLowFreq(index_t& i) const
{
DataContainer<data_t> sectionZero(VolumeDescriptor{{_width, _height}});
sectionZero = 0;
auto negativeHalfWidth = static_cast<index_t>(-std::floor(_width / 2.0));
auto halfWidth = static_cast<index_t>(std::ceil(_width / 2.0));
auto negativeHalfHeight = static_cast<index_t>(-std::floor(_height / 2.0));
auto halfHeight = static_cast<index_t>(std::ceil(_height / 2.0));
// TODO attempt to refactor the negative indexing
for (auto w = negativeHalfWidth; w < halfWidth; w++) {
for (auto h = negativeHalfHeight; h < halfHeight; h++) {
sectionZero(w < 0 ? w + _width : w, h < 0 ? h + _height : h) =
shearlet::phiHat<data_t>(static_cast<data_t>(w), static_cast<data_t>(h));
}
}
_spectra.value().slice(i++) = sectionZero;
}
template <typename ret_t, typename data_t>
void ShearletTransform<ret_t, data_t>::_computeSpectraAtConicRegions(index_t& i, index_t j,
index_t k) const
{
DataContainer<data_t> sectionh(VolumeDescriptor{{_width, _height}});
sectionh = 0;
DataContainer<data_t> sectionv(VolumeDescriptor{{_width, _height}});
sectionv = 0;
auto negativeHalfWidth = static_cast<index_t>(-std::floor(_width / 2.0));
auto halfWidth = static_cast<index_t>(std::ceil(_width / 2.0));
auto negativeHalfHeight = static_cast<index_t>(-std::floor(_height / 2.0));
auto halfHeight = static_cast<index_t>(std::ceil(_height / 2.0));
// TODO attempt to refactor the negative indexing
for (auto w = negativeHalfWidth; w < halfWidth; w++) {
for (auto h = negativeHalfHeight; h < halfHeight; h++) {
if (std::abs(h) <= std::abs(w)) {
sectionh(w < 0 ? w + _width : w, h < 0 ? h + _height : h) =
shearlet::psiHat<data_t>(std::pow(4, -j) * w,
std::pow(4, -j) * k * w + std::pow(2, -j) * h);
} else {
sectionv(w < 0 ? w + _width : w, h < 0 ? h + _height : h) =
shearlet::psiHat<data_t>(std::pow(4, -j) * h,
std::pow(4, -j) * k * h + std::pow(2, -j) * w);
}
}
}
_spectra.value().slice(i++) = sectionh;
_spectra.value().slice(i++) = sectionv;
}
template <typename ret_t, typename data_t>
void ShearletTransform<ret_t, data_t>::_computeSpectraAtSeamLines(index_t& i, index_t j,
index_t k) const
{
DataContainer<data_t> sectionhxv(VolumeDescriptor{{_width, _height}});
sectionhxv = 0;
auto negativeHalfWidth = static_cast<index_t>(-std::floor(_width / 2.0));
auto halfWidth = static_cast<index_t>(std::ceil(_width / 2.0));
auto negativeHalfHeight = static_cast<index_t>(-std::floor(_height / 2.0));
auto halfHeight = static_cast<index_t>(std::ceil(_height / 2.0));
// TODO attempt to refactor the negative indexing
for (auto w = negativeHalfWidth; w < halfWidth; w++) {
for (auto h = negativeHalfHeight; h < halfHeight; h++) {
if (std::abs(h) <= std::abs(w)) {
sectionhxv(w < 0 ? w + _width : w, h < 0 ? h + _height : h) =
shearlet::psiHat<data_t>(std::pow(4, -j) * w,
std::pow(4, -j) * k * w + std::pow(2, -j) * h);
} else {
sectionhxv(w < 0 ? w + _width : w, h < 0 ? h + _height : h) =
shearlet::psiHat<data_t>(std::pow(4, -j) * h,
std::pow(4, -j) * k * h + std::pow(2, -j) * w);
}
}
}
_spectra.value().slice(i++) = sectionhxv;
}
template <typename ret_t, typename data_t>
auto ShearletTransform<ret_t, data_t>::getSpectra() const -> DataContainer<data_t>
{
if (!_spectra.has_value()) {
throw LogicError(std::string("ShearletTransform: the spectra is not yet computed"));
}
return _spectra.value();
}
template <typename ret_t, typename data_t>
bool ShearletTransform<ret_t, data_t>::isSpectraComputed() const
{
return _spectra.has_value();
}
template <typename ret_t, typename data_t>
index_t ShearletTransform<ret_t, data_t>::calculateNumOfScales(index_t width, index_t height)
{
return static_cast<index_t>(std::log2(std::max(width, height)) / 2.0);
}
template <typename ret_t, typename data_t>
index_t ShearletTransform<ret_t, data_t>::calculateNumOfLayers(index_t width, index_t height)
{
return static_cast<index_t>(std::pow(2, (calculateNumOfScales(width, height) + 2)) - 3);
}
template <typename ret_t, typename data_t>
index_t ShearletTransform<ret_t, data_t>::calculateNumOfLayers(index_t numOfScales)
{
return static_cast<index_t>(std::pow(2, numOfScales + 2) - 3);
}
template <typename ret_t, typename data_t>
auto ShearletTransform<ret_t, data_t>::getWidth() const -> index_t
{
return _width;
}
template <typename ret_t, typename data_t>
auto ShearletTransform<ret_t, data_t>::getHeight() const -> index_t
{
return _height;
}
template <typename ret_t, typename data_t>
auto ShearletTransform<ret_t, data_t>::getNumOfLayers() const -> index_t
{
return _numOfLayers;
}
template <typename ret_t, typename data_t>
ShearletTransform<ret_t, data_t>* ShearletTransform<ret_t, data_t>::cloneImpl() const
{
return new ShearletTransform<ret_t, data_t>(_width, _height, _numOfScales, _spectra);
}
template <typename ret_t, typename data_t>
bool ShearletTransform<ret_t, data_t>::isEqual(const LinearOperator<ret_t>& other) const
{
if (!LinearOperator<ret_t>::isEqual(other))
return false;
auto otherST = downcast_safe<ShearletTransform<ret_t, data_t>>(&other);
if (!otherST)
return false;
if (_width != otherST->_width)
return false;
if (_height != otherST->_height)
return false;
if (_numOfScales != otherST->_numOfScales)
return false;
return true;
}
// ------------------------------------------
// explicit template instantiation
template class ShearletTransform<float, float>;
template class ShearletTransform<elsa::complex<float>, float>;
template class ShearletTransform<double, double>;
template class ShearletTransform<elsa::complex<double>, double>;
} // namespace elsa
#pragma once
#include "LinearOperator.h"
namespace elsa
{
/**
* @brief Class representing a (regular) Cone-Adapted Discrete Shearlet Transform
*
* @author Andi Braimllari - initial code
*
* @tparam data_t data type for the domain and range of the problem, defaulting to real_t
*
* ShearletTransform represents a band-limited (compact support in Fourier domain)
* representation system. It oversamples a 2D signal of (W, H) to (W, H, L). Most of the
* computation is taken for the spectra, which is stored after the first run. It only
* handles signals with one channel, e.g. grayscale images. Increasing the number of scales
* will increase precision.