1 /*
2 * This file is part of RawTherapee.
3 *
4 * Copyright (c) 2004-2018 Gabor Horvath <hgabor@rawtherapee.com> and other RawTherapee contributors
5 * Split out to own compilation unit and made speedup 2018 Ingo Weyrich (heckflosse67@gmx.de)
6 *
7 * RawTherapee is free software: you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation, either version 3 of the License, or
10 * (at your option) any later version.
11 *
12 * RawTherapee is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with RawTherapee. If not, see <http://www.gnu.org/licenses/>.
19 */
20 #include <cmath>
21
22 #include "color.h"
23 #include "rawimagesource.h"
24 #include "rawimagesource_i.h"
25 #include "jaggedarray.h"
26 #include "rawimage.h"
27 #include "iccmatrices.h"
28 #include "rt_math.h"
29 #include "../rtgui/multilangmgr.h"
30 #include "procparams.h"
31 //#define BENCHMARK
32 #include "StopWatch.h"
33
34 using namespace std;
35
36 namespace rtengine
37 {
38
39
interpolate_row_g(float * agh,float * agv,int i)40 inline void RawImageSource::interpolate_row_g (float* agh, float* agv, int i)
41 {
42
43 for (int j = 0; j < W; j++) {
44 if (ri->ISGREEN(i, j)) {
45 agh[j] = rawData[i][j];
46 agv[j] = rawData[i][j];
47 } else {
48 int gh = 0;
49 int gv = 0;
50
51 if (j > 1 && j < W - 2) {
52 gh = (-rawData[i][j - 2] + 2 * rawData[i][j - 1] + 2 * rawData[i][j] + 2 * rawData[i][j + 1] - rawData[i][j + 2]) / 4;
53 int maxgh = max(rawData[i][j - 1], rawData[i][j + 1]);
54 int mingh = min(rawData[i][j - 1], rawData[i][j + 1]);
55
56 if (gh > maxgh) {
57 gh = maxgh;
58 } else if (gh < mingh) {
59 gh = mingh;
60 }
61 } else if (j == 0) {
62 gh = rawData[i][1];
63 } else if (j == 1) {
64 gh = (rawData[i][0] + rawData[i][2]) / 2;
65 } else if (j == W - 1) {
66 gh = rawData[i][W - 2];
67 } else if (j == W - 2) {
68 gh = (rawData[i][W - 1] + rawData[i][W - 3]) / 2;
69 }
70
71 if (i > 1 && i < H - 2) {
72 gv = (-rawData[i - 2][j] + 2 * rawData[i - 1][j] + 2 * rawData[i][j] + 2 * rawData[i + 1][j] - rawData[i + 2][j]) / 4;
73 int maxgv = max(rawData[i - 1][j], rawData[i + 1][j]);
74 int mingv = min(rawData[i - 1][j], rawData[i + 1][j]);
75
76 if (gv > maxgv) {
77 gv = maxgv;
78 } else if (gv < mingv) {
79 gv = mingv;
80 }
81 } else if (i == 0) {
82 gv = rawData[1][j];
83 } else if (i == 1) {
84 gv = (rawData[0][j] + rawData[2][j]) / 2;
85 } else if (i == H - 1) {
86 gv = rawData[H - 2][j];
87 } else if (i == H - 2) {
88 gv = (rawData[H - 1][j] + rawData[H - 3][j]) / 2;
89 }
90
91 agh[j] = gh;
92 agv[j] = gv;
93 }
94 }
95 }
96
interpolate_row_rb(float * ar,float * ab,float * pg,float * cg,float * ng,int i)97 inline void RawImageSource::interpolate_row_rb (float* ar, float* ab, float* pg, float* cg, float* ng, int i)
98 {
99 const auto getPg = [pg](int index) {
100 return
101 pg
102 ? pg[index]
103 : 0.f;
104 };
105
106 const auto getNg = [ng](int index) {
107 return
108 ng
109 ? ng[index]
110 : 0.f;
111 };
112
113 float *nonGreen1 = ar;
114 float *nonGreen2 = ab;
115
116 if ((ri->ISBLUE(i, 0) || ri->ISBLUE(i, 1))) {
117 std::swap(nonGreen1, nonGreen2);
118 }
119 int j = FC(i, 0) & 1;
120 if (j) {
121 // linear R-G interp. horizontally
122 float val1;
123
124 val1 = cg[0] + rawData[i][1] - cg[1];
125
126 nonGreen1[0] = CLIP(val1);
127 // linear B-G interp. vertically
128 float val2;
129
130 if (i == 0) {
131 val2 = getNg(0) + rawData[1][0] - cg[0];
132 } else if (i == H - 1) {
133 val2 = getPg(0) + rawData[H - 2][0] - cg[0];
134 } else {
135 val2 = cg[0] + (rawData[i - 1][0] - getPg(0) + rawData[i + 1][0] - getNg(0)) / 2;
136 }
137
138 nonGreen2[0] = val2;
139 }
140 // RGRGR or GRGRGR line
141 for (; j < W - 1; j += 2) {
142 // nonGreen of row is simple
143 nonGreen1[j] = rawData[i][j];
144 // non green of next row: cross interpolation
145 float nonGreen = 0.f;
146 int n = 0;
147
148 if (i > 0) {
149 if (j > 0) {
150 nonGreen += rawData[i - 1][j - 1] - getPg(j - 1);
151 n++;
152 }
153 nonGreen += rawData[i - 1][j + 1] - getPg(j + 1);
154 n++;
155 }
156
157 if (i < H - 1) {
158 if (j > 0) {
159 nonGreen += rawData[i + 1][j - 1] - getNg(j - 1);
160 n++;
161 }
162 nonGreen += rawData[i + 1][j + 1] - getNg(j + 1);
163 n++;
164 }
165
166 nonGreen2[j] = cg[j] + nonGreen / std::max(n, 1);
167
168 // linear R-G interp. horizontally
169 float val1;
170
171 if (j == W - 2) {
172 val1 = cg[W - 1] + rawData[i][W - 2] - cg[W - 2];
173 } else {
174 val1 = cg[j + 1] + (rawData[i][j] - cg[j] + rawData[i][j + 2] - cg[j + 2]) / 2;
175 }
176
177 nonGreen1[j + 1] = CLIP(val1);
178 // linear B-G interp. vertically
179 float val2;
180
181 if (i == 0) {
182 val2 = getNg(j + 1) + rawData[1][j + 1] - cg[j + 1];
183 } else if (i == H - 1) {
184 val2 = getPg(j + 1) + rawData[H - 2][j + 1] - cg[j + 1];
185 } else {
186 val2 = cg[j + 1] + (rawData[i - 1][j + 1] - getPg(j + 1) + rawData[i + 1][j + 1] - getNg(j + 1)) / 2;
187 }
188
189 nonGreen2[j + 1] = val2;
190 }
191
192 if(j == W - 1) {
193 // nonGreen of row is simple
194 nonGreen1[j] = rawData[i][j];
195 // non green of next row: cross interpolation
196 float nonGreen = 0.f;
197 int n = 0;
198
199 if (i > 0) {
200 nonGreen += rawData[i - 1][j - 1] - getPg(j - 1);
201 n++;
202 }
203
204 if (i < H - 1) {
205 nonGreen += rawData[i + 1][j - 1] - getNg(j - 1);
206 n++;
207 }
208
209 nonGreen2[j] = cg[j] + nonGreen / std::max(n, 1);
210 }
211 }
212
213 #define DIST(a,b) (std::fabs(a-b))
214
eahd_demosaic()215 void RawImageSource::eahd_demosaic ()
216 {
217 BENCHFUN
218 if (plistener) {
219 plistener->setProgressStr (Glib::ustring::compose(M("TP_RAW_DMETHOD_PROGRESSBAR"), RAWParams::BayerSensor::getMethodString(RAWParams::BayerSensor::Method::EAHD)));
220 plistener->setProgress (0.0);
221 }
222
223 // prepare constants for cielab conversion
224 //TODO: revisit after conversion to D50 illuminant
225 const float lc00 = (0.412453 * imatrices.rgb_cam[0][0] + 0.357580 * imatrices.rgb_cam[0][1] + 0.180423 * imatrices.rgb_cam[0][2]) ;// / 0.950456;
226 const float lc01 = (0.412453 * imatrices.rgb_cam[1][0] + 0.357580 * imatrices.rgb_cam[1][1] + 0.180423 * imatrices.rgb_cam[1][2]) ;// / 0.950456;
227 const float lc02 = (0.412453 * imatrices.rgb_cam[2][0] + 0.357580 * imatrices.rgb_cam[2][1] + 0.180423 * imatrices.rgb_cam[2][2]) ;// / 0.950456;
228
229 const float lc10 = 0.212671 * imatrices.rgb_cam[0][0] + 0.715160 * imatrices.rgb_cam[0][1] + 0.072169 * imatrices.rgb_cam[0][2];
230 const float lc11 = 0.212671 * imatrices.rgb_cam[1][0] + 0.715160 * imatrices.rgb_cam[1][1] + 0.072169 * imatrices.rgb_cam[1][2];
231 const float lc12 = 0.212671 * imatrices.rgb_cam[2][0] + 0.715160 * imatrices.rgb_cam[2][1] + 0.072169 * imatrices.rgb_cam[2][2];
232
233 const float lc20 = (0.019334 * imatrices.rgb_cam[0][0] + 0.119193 * imatrices.rgb_cam[0][1] + 0.950227 * imatrices.rgb_cam[0][2]) ;// / 1.088754;
234 const float lc21 = (0.019334 * imatrices.rgb_cam[1][0] + 0.119193 * imatrices.rgb_cam[1][1] + 0.950227 * imatrices.rgb_cam[1][2]) ;// / 1.088754;
235 const float lc22 = (0.019334 * imatrices.rgb_cam[2][0] + 0.119193 * imatrices.rgb_cam[2][1] + 0.950227 * imatrices.rgb_cam[2][2]) ;// / 1.088754;
236
237 const float wp[3][3] = {{lc00, lc01, lc02}, {lc10, lc11, lc12}, {lc20, lc21, lc22}};
238
239 // end of cielab preparation
240
241 JaggedArray<float>
242 rh (W, 3), gh (W, 4), bh (W, 3),
243 rv (W, 3), gv (W, 4), bv (W, 3),
244 lLh (W, 3), lah (W, 3), lbh (W, 3),
245 lLv (W, 3), lav (W, 3), lbv (W, 3);
246 JaggedArray<uint16_t> homh (W, 3), homv (W, 3);
247
248 // interpolate first two lines
249 interpolate_row_g (gh[0], gv[0], 0);
250 interpolate_row_g (gh[1], gv[1], 1);
251 interpolate_row_g (gh[2], gv[2], 2);
252 interpolate_row_rb (rh[0], bh[0], nullptr, gh[0], gh[1], 0);
253 interpolate_row_rb (rv[0], bv[0], nullptr, gv[0], gv[1], 0);
254 interpolate_row_rb (rh[1], bh[1], gh[0], gh[1], gh[2], 1);
255 interpolate_row_rb (rv[1], bv[1], gv[0], gv[1], gv[2], 1);
256
257 Color::RGB2Lab(rh[0], gh[0], bh[0], lLh[0], lah[0], lbh[0], wp, W);
258 Color::RGB2Lab(rv[0], gv[0], bv[0], lLv[0], lav[0], lbv[0], wp, W);
259 Color::RGB2Lab(rh[1], gh[1], bh[1], lLh[1], lah[1], lbh[1], wp, W);
260 Color::RGB2Lab(rv[1], gv[1], bv[1], lLv[1], lav[1], lbv[1], wp, W);
261
262 for (int j = 0; j < W; j++) {
263 homh[0][j] = 0;
264 homv[0][j] = 0;
265 homh[1][j] = 0;
266 homv[1][j] = 0;
267 }
268
269 float dLmaph[9];
270 float dLmapv[9];
271 float dCamaph[9];
272 float dCamapv[9];
273 float dCbmaph[9];
274 float dCbmapv[9];
275
276 for (int i = 1; i < H - 1; i++) {
277 int mod[3] = {(i-1) % 3, i % 3, (i+1) %3};
278 int ix = i % 3;
279 int imx = (i - 1) % 3;
280 int ipx = (i + 1) % 3;
281
282 if (i < H - 2) {
283 interpolate_row_g (gh[(i + 2) % 4], gv[(i + 2) % 4], i + 2);
284 interpolate_row_rb (rh[(i + 1) % 3], bh[(i + 1) % 3], gh[i % 4], gh[(i + 1) % 4], gh[(i + 2) % 4], i + 1);
285 interpolate_row_rb (rv[(i + 1) % 3], bv[(i + 1) % 3], gv[i % 4], gv[(i + 1) % 4], gv[(i + 2) % 4], i + 1);
286 } else {
287 interpolate_row_rb (rh[(i + 1) % 3], bh[(i + 1) % 3], gh[i % 4], gh[(i + 1) % 4], nullptr, i + 1);
288 interpolate_row_rb (rv[(i + 1) % 3], bv[(i + 1) % 3], gv[i % 4], gv[(i + 1) % 4], nullptr, i + 1);
289 }
290
291 Color::RGB2Lab(rh[(i + 1) % 3], gh[(i + 1) % 4], bh[(i + 1) % 3], lLh[(i + 1) % 3], lah[(i + 1) % 3], lbh[(i + 1) % 3], wp, W);
292 Color::RGB2Lab(rv[(i + 1) % 3], gv[(i + 1) % 4], bv[(i + 1) % 3], lLv[(i + 1) % 3], lav[(i + 1) % 3], lbv[(i + 1) % 3], wp, W);
293
294 for (int j = 0; j < W; j++) {
295 homh[ipx][j] = 0;
296 homv[ipx][j] = 0;
297 }
298
299 for (int j = 1; j < W - 1; j++) {
300 int dmi = 0;
301 for (int x = -1; x <= 0; x++) {
302 int idx = mod[x + 1];
303
304 for (int y = -1; y <= 1; y++) {
305 // compute distance in a, b, and L
306 if (dmi < 4) {
307 int sh = homh[idx][j + y];
308 int sv = homv[idx][j + y];
309
310 if (sh > sv) { // fixate horizontal pixel
311 dLmaph[dmi] = DIST(lLh[ix][j], lLh[idx][j + y]);
312 dCamaph[dmi] = DIST(lah[ix][j], lah[idx][j + y]);
313 dCbmaph[dmi] = DIST(lbh[ix][j], lbh[idx][j + y]);
314 dLmapv[dmi] = DIST(lLv[ix][j], lLh[idx][j + y]);
315 dCamapv[dmi] = DIST(lav[ix][j], lah[idx][j + y]);
316 dCbmapv[dmi] = DIST(lbv[ix][j], lbh[idx][j + y]);
317 } else if (sh < sv) {
318 dLmaph[dmi] = DIST(lLh[ix][j], lLv[idx][j + y]);
319 dCamaph[dmi] = DIST(lah[ix][j], lav[idx][j + y]);
320 dCbmaph[dmi] = DIST(lbh[ix][j], lbv[idx][j + y]);
321 dLmapv[dmi] = DIST(lLv[ix][j], lLv[idx][j + y]);
322 dCamapv[dmi] = DIST(lav[ix][j], lav[idx][j + y]);
323 dCbmapv[dmi] = DIST(lbv[ix][j], lbv[idx][j + y]);
324 } else {
325 dLmaph[dmi] = DIST(lLh[ix][j], lLh[idx][j + y]);
326 dCamaph[dmi] = DIST(lah[ix][j], lah[idx][j + y]);
327 dCbmaph[dmi] = DIST(lbh[ix][j], lbh[idx][j + y]);
328 dLmapv[dmi] = DIST(lLv[ix][j], lLv[idx][j + y]);
329 dCamapv[dmi] = DIST(lav[ix][j], lav[idx][j + y]);
330 dCbmapv[dmi] = DIST(lbv[ix][j], lbv[idx][j + y]);
331 }
332 } else {
333 dLmaph[dmi] = DIST(lLh[ix][j], lLh[idx][j + y]);
334 dCamaph[dmi] = DIST(lah[ix][j], lah[idx][j + y]);
335 dCbmaph[dmi] = DIST(lbh[ix][j], lbh[idx][j + y]);
336 dLmapv[dmi] = DIST(lLv[ix][j], lLv[idx][j + y]);
337 dCamapv[dmi] = DIST(lav[ix][j], lav[idx][j + y]);
338 dCbmapv[dmi] = DIST(lbv[ix][j], lbv[idx][j + y]);
339 }
340
341 dmi++;
342 }
343 }
344 int idx = mod[2];
345
346 for (int y = -1; y <= 1; y++) {
347 // compute distance in a, b, and L
348 dLmaph[dmi] = DIST(lLh[ix][j], lLh[idx][j + y]);
349 dCamaph[dmi] = DIST(lah[ix][j], lah[idx][j + y]);
350 dCbmaph[dmi] = DIST(lbh[ix][j], lbh[idx][j + y]);
351 dLmapv[dmi] = DIST(lLv[ix][j], lLv[idx][j + y]);
352 dCamapv[dmi] = DIST(lav[ix][j], lav[idx][j + y]);
353 dCbmapv[dmi] = DIST(lbv[ix][j], lbv[idx][j + y]);
354 dmi++;
355 }
356
357 // compute eL & eC
358 float eL = min(max(dLmaph[3], dLmaph[5]), max(dLmapv[1], dLmapv[7]));
359 float eCa = min(max(dCamaph[3], dCamaph[5]), max(dCamapv[1], dCamapv[7]));
360 float eCb = min(max(dCbmaph[3], dCbmaph[5]), max(dCbmapv[1], dCbmapv[7]));
361
362 int wh = 0;
363
364 for (int dmi = 0; dmi < 9; dmi++) {
365 wh += (dLmaph[dmi] <= eL) * (dCamaph[dmi] <= eCa) * (dCbmaph[dmi] <= eCb);
366 }
367
368 homh[imx][j - 1] += wh;
369 homh[imx][j] += wh;
370 homh[imx][j + 1] += wh;
371 homh[ix][j - 1] += wh;
372 homh[ix][j] += wh;
373 homh[ix][j + 1] += wh;
374 homh[ipx][j - 1] += wh;
375 homh[ipx][j] += wh;
376 homh[ipx][j + 1] += wh;
377
378 int wv = 0;
379
380 for (int dmi = 0; dmi < 9; dmi++) {
381 wv += (dLmapv[dmi] <= eL) * (dCamapv[dmi] <= eCa) * (dCbmapv[dmi] <= eCb);
382 }
383
384 homv[imx][j - 1] += wv;
385 homv[imx][j] += wv;
386 homv[imx][j + 1] += wv;
387 homv[ix][j - 1] += wv;
388 homv[ix][j] += wv;
389 homv[ix][j + 1] += wv;
390 homv[ipx][j - 1] += wv;
391 homv[ipx][j] += wv;
392 homv[ipx][j + 1] += wv;
393 }
394
395 // finalize image
396 for (int j = 0; j < W; j++) {
397 if (ri->ISGREEN(i - 1, j)) {
398 green[i - 1][j] = rawData[i - 1][j];
399 } else {
400 int hc = homh[imx][j];
401 int vc = homv[imx][j];
402
403 if (hc > vc) {
404 green[i - 1][j] = std::max(0.f, gh[(i - 1) % 4][j]);
405 } else if (hc < vc) {
406 green[i - 1][j] = std::max(0.f, gv[(i - 1) % 4][j]);
407 } else {
408 green[i - 1][j] = std::max(0.f, (gh[(i - 1) % 4][j] + gv[(i - 1) % 4][j]) / 2);
409 }
410 }
411 }
412
413 if (!(i % 20) && plistener) {
414 plistener->setProgress ((double)i / (H - 2));
415 }
416 }
417
418 // finish H-2th and H-1th row, homogenity value is still valailable
419 for (int i = H - 1; i < H + 1; i++)
420 for (int j = 0; j < W; j++) {
421 int hc = homh[(i - 1) % 3][j];
422 int vc = homv[(i - 1) % 3][j];
423
424 if (hc > vc) {
425 green[i - 1][j] = std::max(0.f, gh[(i - 1) % 4][j]);
426 } else if (hc < vc) {
427 green[i - 1][j] = std::max(0.f, gv[(i - 1) % 4][j]);
428 } else {
429 green[i - 1][j] = std::max(0.f, (gh[(i - 1) % 4][j] + gv[(i - 1) % 4][j]) / 2);
430 }
431 }
432
433 // Interpolate R and B
434 #ifdef _OPENMP
435 #pragma omp parallel for
436 #endif
437 for (int i = 0; i < H; i++) {
438 if (i == 0) {
439 interpolate_row_rb_mul_pp (rawData, red[i], blue[i], nullptr, green[i], green[i + 1], i, 1.0, 1.0, 1.0, 0, W, 1);
440 } else if (i == H - 1) {
441 interpolate_row_rb_mul_pp (rawData, red[i], blue[i], green[i - 1], green[i], nullptr, i, 1.0, 1.0, 1.0, 0, W, 1);
442 } else {
443 interpolate_row_rb_mul_pp (rawData, red[i], blue[i], green[i - 1], green[i], green[i + 1], i, 1.0, 1.0, 1.0, 0, W, 1);
444 }
445 }
446 }
447
448 }
449