1 /*
2  * Copyright 2021 The Android Open Source Project
3  *
4  * Licensed under the Apache License, Version 2.0 (the "License");
5  * you may not use this file except in compliance with the License.
6  * You may obtain a copy of the License at
7  *
8  *      http://www.apache.org/licenses/LICENSE-2.0
9  *
10  * Unless required by applicable law or agreed to in writing, software
11  * distributed under the License is distributed on an "AS IS" BASIS,
12  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13  * See the License for the specific language governing permissions and
14  * limitations under the License.
15  */
16 
17 #pragma once
18 
19 #include <cmath> // for sqrt() and ceil()
20 #include <memory>
21 #include <vector>
22 
23 #include "array2d.h"
24 #include "fast_dct.h"
25 
26 #include "fft.h"   // external FFT routines
27 #include "fft2d.h" // external 2D FFT routines
28 
29 namespace android {
30 
31 template <typename Torig, typename Tdct>
32 FastDCT<Torig, Tdct>::FastDCT(int block_size)
33       : block_size_(block_size),
34         intermediate_scale_factor_(1.0),
35         orth_scale_factor_(2.0 / static_cast<double>(block_size_)),
36         sqrt_orth_scale_factor_(sqrt(orth_scale_factor_)),
37         sqrt_one_half_(sqrt(0.5)),
38         ddct2d_in_out_ptr_array_(new double *[block_size_]),
39         contiguous_double_array_(new double[block_size_ * block_size_]),
40         cos_sin_table_(new double[block_size * 3 / 2 + 1]),
41         ddct_in_out_(new double[block_size_]),
42         t_workspace_(new double[4 * block_size]),
43         bit_reversal_workspace_(new int[2 + static_cast<int>(ceil(sqrt(block_size)))]) {
44     // Check that block_size_ is a power of two
45     int block_size_shifted = block_size_;
46     int num_shifts = 0;
47     while (block_size_shifted >>= 1) num_shifts++;
48 
49     // As per the comments above the definition of ddct and ddct2d in
50     // the external fft2d library, setting bit_reversal_workspace_[0] to
51     // zero will cause it and the cos_sin table to be initialized the
52     // first time one of the low-level C routines is called.
53     bit_reversal_workspace_[0] = 0;
54 
55     // The ddct routine requires that the elements of
56     // ddct2d_in_out_ptr_array_ point to successive block_size_-sized
57     // blocks in a contiguous array.
58     ddct2d_in_out_ptr_array_[0] = contiguous_double_array_.get();
59     for (int i = 1; i < block_size_; i++) {
60         ddct2d_in_out_ptr_array_[i] = ddct2d_in_out_ptr_array_[i - 1] + block_size_;
61     }
62 }
63 
64 // Perform one-dimensional forward DCT.
65 template <typename Torig, typename Tdct>
66 void FastDCT<Torig, Tdct>::ForwardTransform1D(const std::vector<Torig> &v_in,
67                                               std::vector<Tdct> *v_out) const {
68     for (int i = 0; i < block_size_; i++) {
69         ddct_in_out_[i] = static_cast<double>(v_in[i]);
70     }
71 
72     ddct(block_size_, -1, ddct_in_out_.get(), bit_reversal_workspace_.get(), cos_sin_table_.get());
73 
74     for (int i = 0; i < block_size_; i++) {
75         ddct_in_out_[i] *= sqrt_orth_scale_factor_ * intermediate_scale_factor_;
76     }
77     ddct_in_out_[0] *= sqrt_one_half_;
78     for (int i = 0; i < block_size_; i++) {
79         (*v_out)[i] = static_cast<Tdct>(ddct_in_out_[i] + 0.5) - static_cast<Tdct>(0.5);
80     }
81 }
82 
83 // Perform one-dimensional reverse (inverse) DCT.
84 template <typename Torig, typename Tdct>
85 void FastDCT<Torig, Tdct>::ReverseTransform1D(const std::vector<Tdct> &v_in,
86                                               std::vector<Torig> *v_out) const {
87     for (int i = 0; i < block_size_; i++) {
88         ddct_in_out_[i] = static_cast<double>(v_in[i]);
89     }
90     ddct_in_out_[0] *= sqrt_one_half_;
91 
92     ddct(block_size_, 1, ddct_in_out_.get(), bit_reversal_workspace_.get(), cos_sin_table_.get());
93 
94     for (int i = 0; i < block_size_; i++) {
95         ddct_in_out_[i] *= sqrt_orth_scale_factor_ / intermediate_scale_factor_;
96     }
97     for (int i = 0; i < block_size_; i++) {
98         (*v_out)[i] = static_cast<Torig>(ddct_in_out_[i] + 0.5) - static_cast<Torig>(0.5);
99     }
100 }
101 
102 // Perform two-dimensional forward DCT.
103 template <typename Torig, typename Tdct>
104 void FastDCT<Torig, Tdct>::ForwardTransform2D(const Array2D<Torig> &m_in,
105                                               Array2D<Tdct> *m_out) const {
106     for (int row = 0; row < block_size_; row++) {
107         for (int col = 0; col < block_size_; col++) {
108             ddct2d_in_out_ptr_array_[row][col] = static_cast<double>(m_in(row, col));
109         }
110     }
111 
112     ddct2d(block_size_, block_size_, -1, ddct2d_in_out_ptr_array_.get(), t_workspace_.get(),
113            bit_reversal_workspace_.get(), cos_sin_table_.get());
114 
115     double isf_squared = intermediate_scale_factor_ * intermediate_scale_factor_;
116     for (int row = 0; row < block_size_; row++) {
117         for (int col = 0; col < block_size_; col++) {
118             ddct2d_in_out_ptr_array_[row][col] *= orth_scale_factor_ * isf_squared;
119         }
120     }
121     for (int i = 0; i < block_size_; i++) {
122         ddct2d_in_out_ptr_array_[i][0] *= sqrt_one_half_;
123         ddct2d_in_out_ptr_array_[0][i] *= sqrt_one_half_;
124     }
125 
126     for (int row = 0; row < block_size_; row++) {
127         for (int col = 0; col < block_size_; col++) {
128             (*m_out)(row, col) = static_cast<Tdct>(ddct2d_in_out_ptr_array_[row][col] + 0.5) -
129                     static_cast<Tdct>(0.5);
130         }
131     }
132 }
133 
134 // Perform two-dimensional reverse (inverse) DCT.
135 template <typename Torig, typename Tdct>
136 void FastDCT<Torig, Tdct>::ReverseTransform2D(const Array2D<Tdct> &m_in,
137                                               Array2D<Torig> *m_out) const {
138     for (int row = 0; row < block_size_; row++) {
139         for (int col = 0; col < block_size_; col++) {
140             ddct2d_in_out_ptr_array_[row][col] = static_cast<double>(m_in(row, col));
141         }
142     }
143 
144     for (int i = 0; i < block_size_; i++) {
145         ddct2d_in_out_ptr_array_[i][0] *= sqrt_one_half_;
146         ddct2d_in_out_ptr_array_[0][i] *= sqrt_one_half_;
147     }
148 
149     ddct2d(block_size_, block_size_, 1, ddct2d_in_out_ptr_array_.get(), t_workspace_.get(),
150            bit_reversal_workspace_.get(), cos_sin_table_.get());
151 
152     double inv_isf_squared = 1.0 / (intermediate_scale_factor_ * intermediate_scale_factor_);
153     for (int row = 0; row < block_size_; row++) {
154         for (int col = 0; col < block_size_; col++) {
155             ddct2d_in_out_ptr_array_[row][col] *= orth_scale_factor_ * inv_isf_squared;
156         }
157     }
158 
159     for (int row = 0; row < block_size_; row++) {
160         for (int col = 0; col < block_size_; col++) {
161             (*m_out)(row, col) = static_cast<Torig>(ddct2d_in_out_ptr_array_[row][col] + 0.5) -
162                     static_cast<Torig>(0.5);
163         }
164     }
165 }
166 } // namespace android
167