Line data Source code
1 : #include "ShearletTransform.h"
2 : #include "FourierTransform.h"
3 : #include "VolumeDescriptor.h"
4 : #include "Timer.h"
5 : #include "Math.hpp"
6 :
7 : namespace elsa
8 : {
9 : template <typename ret_t, typename data_t>
10 : ShearletTransform<ret_t, data_t>::ShearletTransform(IndexVector_t spatialDimensions)
11 : : ShearletTransform(spatialDimensions[0], spatialDimensions[1])
12 2 : {
13 2 : if (spatialDimensions.size() != 2) {
14 0 : throw LogicError("ShearletTransform: Only 2D shape supported");
15 0 : }
16 2 : }
17 :
18 : template <typename ret_t, typename data_t>
19 : ShearletTransform<ret_t, data_t>::ShearletTransform(index_t width, index_t height)
20 : : ShearletTransform(width, height, calculateNumOfScales(width, height))
21 6 : {
22 6 : }
23 :
24 : template <typename ret_t, typename data_t>
25 : ShearletTransform<ret_t, data_t>::ShearletTransform(index_t width, index_t height,
26 : index_t numOfScales)
27 : : ShearletTransform(width, height, numOfScales, std::nullopt)
28 12 : {
29 12 : }
30 :
31 : template <typename ret_t, typename data_t>
32 : ShearletTransform<ret_t, data_t>::ShearletTransform(
33 : index_t width, index_t height, index_t numOfScales,
34 : std::optional<DataContainer<data_t>> spectra)
35 : : LinearOperator<ret_t>(
36 : VolumeDescriptor{{width, height}},
37 : VolumeDescriptor{{width, height, calculateNumOfLayers(numOfScales)}}),
38 : _spectra{spectra},
39 : _width{width},
40 : _height{height},
41 : _numOfScales{numOfScales},
42 : _numOfLayers{calculateNumOfLayers(numOfScales)}
43 14 : {
44 14 : if (width < 0 || height < 0) {
45 0 : throw LogicError("ShearletTransform: negative width/height were provided");
46 0 : }
47 14 : if (numOfScales < 0) {
48 0 : throw LogicError("ShearletTransform: negative number of scales was provided");
49 0 : }
50 14 : }
51 :
52 : // TODO implement sumByAxis in DataContainer and remove me
53 : template <typename ret_t, typename data_t>
54 : DataContainer<elsa::complex<data_t>> ShearletTransform<ret_t, data_t>::sumByLastAxis(
55 : DataContainer<elsa::complex<data_t>> dc) const
56 2 : {
57 2 : auto coeffsPerDim = dc.getDataDescriptor().getNumberOfCoefficientsPerDimension();
58 2 : index_t width = coeffsPerDim[0];
59 2 : index_t height = coeffsPerDim[1];
60 2 : index_t layers = coeffsPerDim[2];
61 2 : DataContainer<elsa::complex<data_t>> summedDC(VolumeDescriptor{{width, height}});
62 :
63 34 : for (index_t j = 0; j < width; j++) {
64 544 : for (index_t k = 0; k < height; k++) {
65 512 : elsa::complex<data_t> currValue = 0;
66 31744 : for (index_t i = 0; i < layers; i++) {
67 31232 : currValue += dc(j, k, i);
68 31232 : }
69 512 : summedDC(j, k) = currValue;
70 512 : }
71 32 : }
72 :
73 2 : return summedDC;
74 2 : }
75 :
76 : template <typename ret_t, typename data_t>
77 : void ShearletTransform<ret_t, data_t>::applyImpl(const DataContainer<ret_t>& x,
78 : DataContainer<ret_t>& Ax) const
79 2 : {
80 2 : Timer timeguard("ShearletTransform", "apply");
81 :
82 2 : if (_width != this->getDomainDescriptor().getNumberOfCoefficientsPerDimension()[0]
83 2 : || _height != this->getDomainDescriptor().getNumberOfCoefficientsPerDimension()[1]) {
84 0 : throw InvalidArgumentError("ShearletTransform: Width and height of the input do not "
85 0 : "match to that of this shearlet system");
86 0 : }
87 :
88 2 : Logger::get("ShearletTransform")
89 2 : ->info("Running the shearlet transform on a 2D signal of shape ({}, {}), on {} "
90 2 : "scales with an oversampling factor of {} and {} spectra",
91 2 : _width, _height, _numOfScales, _numOfLayers,
92 2 : isSpectraComputed() ? "precomputed" : "non-precomputed");
93 :
94 2 : if (!isSpectraComputed()) {
95 2 : computeSpectra();
96 2 : }
97 :
98 2 : FourierTransform<elsa::complex<data_t>> fourierTransform(x.getDataDescriptor());
99 :
100 2 : DataContainer<elsa::complex<data_t>> fftImg = fourierTransform.apply(x.asComplex());
101 :
102 124 : for (index_t i = 0; i < getNumOfLayers(); i++) {
103 122 : DataContainer<elsa::complex<data_t>> temp =
104 122 : getSpectra().slice(i).viewAs(x.getDataDescriptor()).asComplex() * fftImg;
105 122 : if constexpr (isComplex<ret_t>) {
106 0 : Ax.slice(i) = fourierTransform.applyAdjoint(temp);
107 0 : } else {
108 0 : Ax.slice(i) = real(fourierTransform.applyAdjoint(temp));
109 0 : }
110 122 : }
111 2 : }
112 :
113 : template <typename ret_t, typename data_t>
114 : void ShearletTransform<ret_t, data_t>::applyAdjointImpl(const DataContainer<ret_t>& y,
115 : DataContainer<ret_t>& Aty) const
116 2 : {
117 2 : Timer timeguard("ShearletTransform", "applyAdjoint");
118 :
119 2 : if (_width != this->getDomainDescriptor().getNumberOfCoefficientsPerDimension()[0]
120 2 : || _height != this->getDomainDescriptor().getNumberOfCoefficientsPerDimension()[1]) {
121 0 : throw InvalidArgumentError("ShearletTransform: Width and height of the input do not "
122 0 : "match to that of this shearlet system");
123 0 : }
124 :
125 2 : Logger::get("ShearletTransform")
126 2 : ->info("Running the inverse shearlet transform on a 2D signal of shape ({}, {}), "
127 2 : "on {} "
128 2 : "scales with an oversampling factor of {} and {} spectra",
129 2 : _width, _height, _numOfScales, _numOfLayers,
130 2 : isSpectraComputed() ? "precomputed" : "non-precomputed");
131 :
132 2 : if (!isSpectraComputed()) {
133 0 : computeSpectra();
134 0 : }
135 :
136 2 : FourierTransform<elsa::complex<data_t>> fourierTransform(Aty.getDataDescriptor());
137 :
138 2 : DataContainer<elsa::complex<data_t>> intermRes(y.getDataDescriptor());
139 :
140 124 : for (index_t i = 0; i < getNumOfLayers(); i++) {
141 122 : DataContainer<elsa::complex<data_t>> temp =
142 122 : fourierTransform.apply(y.slice(i).viewAs(Aty.getDataDescriptor()).asComplex())
143 122 : * getSpectra().slice(i).viewAs(Aty.getDataDescriptor()).asComplex();
144 122 : intermRes.slice(i) = fourierTransform.applyAdjoint(temp);
145 122 : }
146 :
147 2 : if constexpr (isComplex<ret_t>) {
148 0 : Aty = sumByLastAxis(intermRes);
149 0 : } else {
150 0 : Aty = real(sumByLastAxis(intermRes));
151 0 : }
152 2 : }
153 :
154 : template <typename ret_t, typename data_t>
155 : void ShearletTransform<ret_t, data_t>::computeSpectra() const
156 6 : {
157 6 : if (isSpectraComputed()) {
158 0 : Logger::get("ShearletTransform")->warn("Spectra have already been computed!");
159 0 : }
160 :
161 6 : _spectra = DataContainer<data_t>(VolumeDescriptor{{_width, _height, _numOfLayers}});
162 :
163 6 : _computeSpectraAtLowFreq();
164 :
165 30 : for (index_t j = 0; j < _numOfScales; j++) {
166 24 : auto twoPowJ = static_cast<index_t>(std::pow(2, j));
167 24 : auto shearletsAtJ = static_cast<index_t>(std::pow(2, j + 2));
168 24 : index_t shearletsUpUntilJ = shearletsAtJ - 3;
169 24 : index_t index = 1;
170 :
171 24 : _computeSpectraAtSeamLines(j, -twoPowJ, shearletsUpUntilJ + twoPowJ);
172 180 : for (auto k = -twoPowJ + 1; k < twoPowJ; k++) {
173 : // modulo instead of remainder for negative numbers is needed here, therefore doing
174 : // "((a % b) + b) % b" instead of "a % b"
175 156 : index_t modIndex =
176 156 : (((twoPowJ - index + 1) % shearletsAtJ) + shearletsAtJ) % shearletsAtJ;
177 156 : if (modIndex == 0) {
178 18 : modIndex = shearletsAtJ - 1;
179 138 : } else {
180 138 : --modIndex;
181 138 : }
182 :
183 156 : _computeSpectraAtConicRegions(j, k, shearletsUpUntilJ + modIndex,
184 156 : shearletsUpUntilJ + twoPowJ + index);
185 156 : ++index;
186 156 : }
187 24 : _computeSpectraAtSeamLines(j, twoPowJ, shearletsUpUntilJ + twoPowJ + index);
188 24 : }
189 6 : }
190 :
191 : template <typename ret_t, typename data_t>
192 : void ShearletTransform<ret_t, data_t>::_computeSpectraAtLowFreq() const
193 6 : {
194 6 : DataContainer<data_t> sectionZero(VolumeDescriptor{{_width, _height}});
195 6 : sectionZero = 0;
196 :
197 6 : auto shape = getShapeFractions();
198 :
199 : // TODO attempt to refactor the negative indexing
200 166 : for (auto w = shape.negativeHalfWidth; w < shape.halfWidth; w++) {
201 4768 : for (auto h = shape.negativeHalfHeight; h < shape.halfHeight; h++) {
202 4608 : sectionZero(w < 0 ? w + _width : w, h < 0 ? h + _height : h) =
203 4608 : shearlet::phiHat<data_t>(static_cast<data_t>(w), static_cast<data_t>(h));
204 4608 : }
205 160 : }
206 :
207 6 : _spectra.value().slice(0) = sectionZero;
208 6 : }
209 :
210 : template <typename ret_t, typename data_t>
211 : void ShearletTransform<ret_t, data_t>::_computeSpectraAtConicRegions(index_t j, index_t k,
212 : index_t hSliceIndex,
213 : index_t vSliceIndex) const
214 156 : {
215 156 : DataContainer<data_t> sectionh(VolumeDescriptor{{_width, _height}});
216 156 : sectionh = 0;
217 156 : DataContainer<data_t> sectionv(VolumeDescriptor{{_width, _height}});
218 156 : sectionv = 0;
219 :
220 156 : auto shape = getShapeFractions();
221 156 : auto jr = static_cast<data_t>(j);
222 156 : auto kr = static_cast<data_t>(k);
223 :
224 : // TODO attempt to refactor the negative indexing
225 4316 : for (auto w = shape.negativeHalfWidth; w < shape.halfWidth; w++) {
226 4160 : auto wr = static_cast<data_t>(w);
227 123968 : for (auto h = shape.negativeHalfHeight; h < shape.halfHeight; h++) {
228 119808 : auto hr = static_cast<data_t>(h);
229 119808 : if (std::abs(h) <= std::abs(w)) {
230 63908 : sectionh(w < 0 ? w + _width : w, h < 0 ? h + _height : h) =
231 63908 : shearlet::psiHat<data_t>(std::pow(4.f, -jr) * wr,
232 63908 : std::pow(4.f, -jr) * kr * wr
233 63908 : + std::pow(2.f, -jr) * hr);
234 63908 : } else {
235 55900 : sectionv(w < 0 ? w + _width : w, h < 0 ? h + _height : h) =
236 55900 : shearlet::psiHat<data_t>(std::pow(4.f, -jr) * hr,
237 55900 : std::pow(4.f, -jr) * kr * hr
238 55900 : + std::pow(2.f, -jr) * wr);
239 55900 : }
240 119808 : }
241 4160 : }
242 :
243 156 : _spectra.value().slice(hSliceIndex) = sectionh;
244 156 : _spectra.value().slice(vSliceIndex) = sectionv;
245 156 : }
246 :
247 : template <typename ret_t, typename data_t>
248 : void ShearletTransform<ret_t, data_t>::_computeSpectraAtSeamLines(index_t j, index_t k,
249 : index_t hxvSliceIndex) const
250 48 : {
251 48 : DataContainer<data_t> sectionhxv(VolumeDescriptor{{_width, _height}});
252 48 : sectionhxv = 0;
253 :
254 48 : auto shape = getShapeFractions();
255 48 : auto jr = static_cast<data_t>(j);
256 48 : auto kr = static_cast<data_t>(k);
257 :
258 : // TODO attempt to refactor the negative indexing
259 1328 : for (auto w = shape.negativeHalfWidth; w < shape.halfWidth; w++) {
260 1280 : auto wr = static_cast<data_t>(w);
261 38144 : for (auto h = shape.negativeHalfHeight; h < shape.halfHeight; h++) {
262 36864 : auto hr = static_cast<data_t>(h);
263 36864 : if (std::abs(h) <= std::abs(w)) {
264 19664 : sectionhxv(w < 0 ? w + _width : w, h < 0 ? h + _height : h) =
265 19664 : shearlet::psiHat<data_t>(std::pow(4.f, -jr) * wr,
266 19664 : std::pow(4.f, -jr) * kr * wr
267 19664 : + std::pow(2.f, -jr) * hr);
268 19664 : } else {
269 17200 : sectionhxv(w < 0 ? w + _width : w, h < 0 ? h + _height : h) =
270 17200 : shearlet::psiHat<data_t>(std::pow(4.f, -jr) * hr,
271 17200 : std::pow(4.f, -jr) * kr * hr
272 17200 : + std::pow(2.f, -jr) * wr);
273 17200 : }
274 36864 : }
275 1280 : }
276 :
277 48 : _spectra.value().slice(hxvSliceIndex) = sectionhxv;
278 48 : }
279 :
280 : /**
281 : * helper function to calculate input data fractions.
282 : */
283 : template <typename ret_t, typename data_t>
284 : auto ShearletTransform<ret_t, data_t>::getShapeFractions() const -> shape_fractions
285 210 : {
286 210 : shape_fractions ret;
287 210 : auto width = static_cast<real_t>(_width);
288 210 : auto height = static_cast<real_t>(_height);
289 :
290 210 : ret.negativeHalfWidth = static_cast<index_t>(-std::floor(width / 2.0));
291 210 : ret.halfWidth = static_cast<index_t>(std::ceil(width / 2.0));
292 210 : ret.negativeHalfHeight = static_cast<index_t>(-std::floor(height / 2.0));
293 210 : ret.halfHeight = static_cast<index_t>(std::ceil(height / 2.0));
294 :
295 210 : return ret;
296 210 : }
297 :
298 : template <typename ret_t, typename data_t>
299 : auto ShearletTransform<ret_t, data_t>::getSpectra() const -> DataContainer<data_t>
300 248 : {
301 248 : if (!_spectra.has_value()) {
302 2 : throw LogicError(std::string("ShearletTransform: the spectra is not yet computed"));
303 2 : }
304 246 : return _spectra.value();
305 246 : }
306 :
307 : template <typename ret_t, typename data_t>
308 : bool ShearletTransform<ret_t, data_t>::isSpectraComputed() const
309 16 : {
310 16 : return _spectra.has_value();
311 16 : }
312 :
313 : template <typename ret_t, typename data_t>
314 : index_t ShearletTransform<ret_t, data_t>::calculateNumOfScales(index_t width, index_t height)
315 6 : {
316 6 : return static_cast<index_t>(std::log2(std::max(width, height)) / 2.0);
317 6 : }
318 :
319 : template <typename ret_t, typename data_t>
320 : index_t ShearletTransform<ret_t, data_t>::calculateNumOfLayers(index_t width, index_t height)
321 0 : {
322 0 : return static_cast<index_t>(std::pow(2, (calculateNumOfScales(width, height) + 2)) - 3);
323 0 : }
324 :
325 : template <typename ret_t, typename data_t>
326 : index_t ShearletTransform<ret_t, data_t>::calculateNumOfLayers(index_t numOfScales)
327 28 : {
328 28 : return static_cast<index_t>(std::pow(2, numOfScales + 2) - 3);
329 28 : }
330 :
331 : template <typename ret_t, typename data_t>
332 : auto ShearletTransform<ret_t, data_t>::getWidth() const -> index_t
333 2 : {
334 2 : return _width;
335 2 : }
336 :
337 : template <typename ret_t, typename data_t>
338 : auto ShearletTransform<ret_t, data_t>::getHeight() const -> index_t
339 2 : {
340 2 : return _height;
341 2 : }
342 :
343 : template <typename ret_t, typename data_t>
344 : auto ShearletTransform<ret_t, data_t>::getNumOfLayers() const -> index_t
345 250 : {
346 250 : return _numOfLayers;
347 250 : }
348 :
349 : template <typename ret_t, typename data_t>
350 : ShearletTransform<ret_t, data_t>* ShearletTransform<ret_t, data_t>::cloneImpl() const
351 2 : {
352 2 : return new ShearletTransform<ret_t, data_t>(_width, _height, _numOfScales, _spectra);
353 2 : }
354 :
355 : template <typename ret_t, typename data_t>
356 : bool ShearletTransform<ret_t, data_t>::isEqual(const LinearOperator<ret_t>& other) const
357 2 : {
358 2 : if (!LinearOperator<ret_t>::isEqual(other))
359 0 : return false;
360 :
361 2 : auto otherST = downcast_safe<ShearletTransform<ret_t, data_t>>(&other);
362 :
363 2 : if (!otherST)
364 0 : return false;
365 :
366 2 : if (_width != otherST->_width)
367 0 : return false;
368 :
369 2 : if (_height != otherST->_height)
370 0 : return false;
371 :
372 2 : if (_numOfScales != otherST->_numOfScales)
373 0 : return false;
374 :
375 2 : return true;
376 2 : }
377 :
378 : // ------------------------------------------
379 : // explicit template instantiation
380 : template class ShearletTransform<float, float>;
381 : template class ShearletTransform<elsa::complex<float>, float>;
382 : template class ShearletTransform<double, double>;
383 : template class ShearletTransform<elsa::complex<double>, double>;
384 : } // namespace elsa
|
