1 /*
2 * This file is part of RawTherapee.
3 *
4 * Copyright (c) 2018 Ingo Weyrich (heckflosse67@gmx.de)
5 *
6 * RawTherapee is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation, either version 3 of the License, or
9 * (at your option) any later version.
10 *
11 * RawTherapee is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with RawTherapee. If not, see <http://www.gnu.org/licenses/>.
18 */
19
20 //
21 // Adaptive Homogeneity-Directed interpolation is based on
22 // the work of Keigo Hirakawa, Thomas Parks, and Paul Lee.
23 // Optimized for speed and reduced memory usage 2018 Ingo Weyrich
24 //
25
26 #include <cmath>
27 #include <climits>
28 #include "bayerhelper.h"
29 #include "LUT.h"
30 #include "librtprocess.h"
31 #include "opthelper.h"
32 #include "rt_math.h"
33 #include "median.h"
34 #include "StopWatch.h"
35
36 #define TS 144
37
38 using namespace librtprocess;
39
ahd_demosaic(int width,int height,const float * const * rawData,float ** red,float ** green,float ** blue,const unsigned cfarray[2][2],const float rgb_cam[3][4],const std::function<bool (double)> & setProgCancel)40 rpError ahd_demosaic(int width, int height, const float * const *rawData, float **red, float **green, float **blue, const unsigned cfarray[2][2], const float rgb_cam[3][4], const std::function<bool(double)> &setProgCancel)
41 {
42 BENCHFUN
43
44 if (!validateBayerCfa(3, cfarray)) {
45 return RP_WRONG_CFA;
46 }
47
48 rpError rc = RP_NO_ERROR;
49
50 constexpr int dir[4] = { -1, 1, -TS, TS };
51 float xyz_cam[3][3];
52 LUTf cbrt(65536);
53
54 constexpr float xyz_rgb[3][3] = { /* XYZ from RGB */
55 { 0.412453, 0.357580, 0.180423 },
56 { 0.212671, 0.715160, 0.072169 },
57 { 0.019334, 0.119193, 0.950227 }
58 };
59
60 constexpr float d65_white[3] = { 0.950456, 1, 1.088754 };
61
62 double progress = 0.0;
63 setProgCancel(progress);
64
65 for (int i = 0; i < 65536; i++) {
66 const double r = i / 65535.0;
67 cbrt[i] = r > 0.008856 ? std::cbrt(r) : 7.787 * r + 16 / 116.0;
68 }
69
70 for (int i = 0; i < 3; i++)
71 for (unsigned int j = 0; j < 3; j++) {
72 xyz_cam[i][j] = 0;
73 for (int k = 0; k < 3; k++) {
74 xyz_cam[i][j] += xyz_rgb[i][k] * rgb_cam[k][j] / d65_white[i];
75 }
76 }
77
78 rc = bayerborder_demosaic(width, height, 5, rawData, red, green, blue, cfarray);
79
80 #ifdef _OPENMP
81 #pragma omp parallel
82 #endif
83 {
84 int progresscounter = 0;
85 float *buffer = new (std::nothrow) float[13 * TS * TS]; /* 1053 kB per core */
86 #ifdef _OPENMP
87 #pragma omp critical
88 #endif
89 {
90 if (!buffer) {
91 rc = RP_MEMORY_ERROR;
92 }
93 }
94 #ifdef _OPENMP
95 #pragma omp barrier
96 #endif
97 if (!rc) {
98 auto rgb = (float(*)[TS][TS][3]) buffer;
99 auto lab = (float(*)[TS][TS][3])(buffer + 6 * TS * TS);
100 auto homo = (uint16_t(*)[TS][TS])(buffer + 12 * TS * TS);
101
102 #ifdef _OPENMP
103 #pragma omp for collapse(2) schedule(dynamic) nowait
104 #endif
105 for (int top = 2; top < height - 5; top += TS - 6) {
106 for (int left = 2; left < width - 5; left += TS - 6) {
107 // Interpolate green horizontally and vertically:
108 for (int row = top; row < top + TS && row < height - 2; row++) {
109 for (int col = left + (fc(cfarray, row, left) & 1); col < std::min(left + TS, width - 2); col += 2) {
110 auto pix = &rawData[row][col];
111 float val0 = 0.25f * ((pix[-1] + pix[0] + pix[1]) * 2
112 - pix[-2] - pix[2]) ;
113 rgb[0][row - top][col - left][1] = median(val0, pix[-1], pix[1]);
114 float val1 = 0.25f * ((pix[-width] + pix[0] + pix[width]) * 2
115 - pix[-2 * width] - pix[2 * width]) ;
116 rgb[1][row - top][col - left][1] = median(val1, pix[-width], pix[width]);
117 }
118 }
119
120 // Interpolate red and blue, and convert to CIELab:
121 for (int d = 0; d < 2; d++)
122 for (int row = top + 1; row < top + TS - 1 && row < height - 3; row++) {
123 int cng = fc(cfarray, row + 1, fc(cfarray, row + 1, 0) & 1);
124 for (int col = left + 1; col < std::min(left + TS - 1, width - 3); col++) {
125 auto pix = &rawData[row][col];
126 auto rix = &rgb[d][row - top][col - left];
127 auto lix = lab[d][row - top][col - left];
128 if (fc(cfarray, row, col) == 1) {
129 rix[0][2 - cng] = CLIP(pix[0] + (0.5f * (pix[-1] + pix[1]
130 - rix[-1][1] - rix[1][1] ) ));
131 rix[0][cng] = CLIP(pix[0] + (0.5f * (pix[-width] + pix[width]
132 - rix[-TS][1] - rix[TS][1])));
133 rix[0][1] = pix[0];
134 } else {
135 rix[0][cng] = CLIP(rix[0][1] + (0.25f * (pix[-width - 1] + pix[-width + 1]
136 + pix[+width - 1] + pix[+width + 1]
137 - rix[-TS - 1][1] - rix[-TS + 1][1]
138 - rix[+TS - 1][1] - rix[+TS + 1][1])));
139 rix[0][2 - cng] = pix[0];
140 }
141 float xyz[3] = {};
142
143 for(unsigned int c = 0; c < 3; ++c) {
144 xyz[0] += xyz_cam[0][c] * rix[0][c];
145 xyz[1] += xyz_cam[1][c] * rix[0][c];
146 xyz[2] += xyz_cam[2][c] * rix[0][c];
147 }
148
149 xyz[0] = cbrt[xyz[0]];
150 xyz[1] = cbrt[xyz[1]];
151 xyz[2] = cbrt[xyz[2]];
152
153 lix[0] = 116.f * xyz[1] - 16.f;
154 lix[1] = 500.f * (xyz[0] - xyz[1]);
155 lix[2] = 200.f * (xyz[1] - xyz[2]);
156 }
157 }
158
159 // Build homogeneity maps from the CIELab images:
160
161 for (int row = top + 2; row < top + TS - 2 && row < height - 4; row++) {
162 int tr = row - top;
163 float ldiff[2][4], abdiff[2][4];
164
165 for (int col = left + 2, tc = 2; col < left + TS - 2 && col < width - 4; col++, tc++) {
166 for (int d = 0; d < 2; d++) {
167 auto lix = &lab[d][tr][tc];
168
169 for (int i = 0; i < 4; i++) {
170 ldiff[d][i] = std::fabs(lix[0][0] - lix[dir[i]][0]);
171 abdiff[d][i] = SQR(lix[0][1] - lix[dir[i]][1])
172 + SQR(lix[0][2] - lix[dir[i]][2]);
173 }
174 }
175
176 float leps = std::min(std::max(ldiff[0][0], ldiff[0][1]),
177 std::max(ldiff[1][2], ldiff[1][3]));
178 float abeps = std::min(std::max(abdiff[0][0], abdiff[0][1]),
179 std::max(abdiff[1][2], abdiff[1][3]));
180
181 for (int d = 0; d < 2; d++) {
182 homo[d][tr][tc] = 0;
183 for (int i = 0; i < 4; i++) {
184 homo[d][tr][tc] += (ldiff[d][i] <= leps) * (abdiff[d][i] <= abeps);
185 }
186 }
187 }
188 }
189
190 // Combine the most homogenous pixels for the final result:
191 for (int row = top + 3; row < top + TS - 3 && row < height - 5; row++) {
192 int tr = row - top;
193
194 for (int col = left + 3, tc = 3; col < std::min(left + TS - 3, width - 5); col++, tc++) {
195 uint16_t hm0 = 0, hm1 = 0;
196 for (int i = tr - 1; i <= tr + 1; i++)
197 for (int j = tc - 1; j <= tc + 1; j++) {
198 hm0 += homo[0][i][j];
199 hm1 += homo[1][i][j];
200 }
201
202 if (hm0 != hm1) {
203 int ldir = hm1 > hm0;
204 red[row][col] = rgb[ldir][tr][tc][0];
205 green[row][col] = rgb[ldir][tr][tc][1];
206 blue[row][col] = rgb[ldir][tr][tc][2];
207 } else {
208 red[row][col] = 0.5f * (rgb[0][tr][tc][0] + rgb[1][tr][tc][0]);
209 green[row][col] = 0.5f * (rgb[0][tr][tc][1] + rgb[1][tr][tc][1]);
210 blue[row][col] = 0.5f * (rgb[0][tr][tc][2] + rgb[1][tr][tc][2]);
211 }
212 }
213 }
214
215
216 progresscounter++;
217 if(progresscounter % 32 == 0) {
218 #ifdef _OPENMP
219 #pragma omp critical (ahdprogress)
220 #endif
221 {
222 progress += 32.0 * SQR(TS - 6) / (height * width);
223 progress = std::min(progress, 1.0);
224 setProgCancel(progress);
225 }
226 }
227
228 }
229 }
230 }
231 delete [] buffer;
232 }
233
234 setProgCancel(1.0);
235
236 return rc;
237 }
238 #undef TS
239
240
241
242
243