1 /*
2 * This file is part of RawTherapee.
3 *
4 * Copyright (c) 2004-2019 Gabor Horvath <hgabor@rawtherapee.com>
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 #include <cmath>
20
21 #include "bayerhelper.h"
22 #include "jaggedarray.h"
23 #include "librtprocess.h"
24 #include "rt_math.h"
25 #include "opthelper.h"
26 #include "StopWatch.h"
27 #ifdef _OPENMP
28 #include <omp.h>
29 #endif
30
31 using namespace std;
32 using namespace librtprocess;
33
34 namespace {
35
hphd_vertical(const float * const * rawData,float ** hpmap,int col_from,int col_to,int H)36 rpError hphd_vertical(const float * const *rawData, float** hpmap, int col_from, int col_to, int H)
37 {
38
39 // process 'numCols' columns for better usage of L1 cpu cache (especially faster for large values of H)
40 constexpr int numCols = 8;
41
42 JaggedArray<float> temp(numCols, H, true);
43 JaggedArray<float> avg(numCols, H, true);
44 JaggedArray<float> dev(numCols, H, true);
45
46 if(!(temp && avg && dev)) {
47 return RP_MEMORY_ERROR;
48 }
49
50 int k = col_from;
51 #ifdef __SSE2__
52 const vfloat ninev = F2V(9.f);
53 const vfloat epsv = F2V(0.001f);
54 #endif
55 for (; k < col_to - 7; k += numCols) {
56 for (int i = 5; i < H - 5; i++) {
57 #ifdef _OPENMP
58 #pragma omp simd
59 #endif
60 for(int h = 0; h < numCols; ++h) {
61 temp[i][h] = std::fabs((rawData[i - 5][k + h] - rawData[i + 5][k + h]) - 8 * (rawData[i - 4][k + h] - rawData[i + 4][k + h]) + 27 * (rawData[i - 3][k + h] - rawData[i + 3][k + h]) - 48 * (rawData[i - 2][k + h] - rawData[i + 2][k + h]) + 42 * (rawData[i - 1][k + h] - rawData[i + 1][k + h]));
62 }
63 }
64
65 for (int j = 4; j < H - 4; j++) {
66 #ifdef __SSE2__
67 // faster than #pragma omp simd...
68 const vfloat avgL1 = ((LVFU(temp[j - 4][0]) + LVFU(temp[j - 3][0])) + (LVFU(temp[j - 2][0]) + LVFU(temp[j - 1][0])) + (LVFU(temp[j][0]) + LVFU(temp[j + 1][0])) + (LVFU(temp[j + 2][0]) + LVFU(temp[j + 3][0])) + LVFU(temp[j + 4][0])) / ninev;
69 STVFU(avg[j][0], avgL1);
70 STVFU(dev[j][0], vmaxf(epsv, (SQRV(LVFU(temp[j - 4][0]) - avgL1) + SQRV(LVFU(temp[j - 3][0]) - avgL1)) + (SQRV(LVFU(temp[j - 2][0]) - avgL1) + SQRV(LVFU(temp[j - 1][0]) - avgL1)) + (SQRV(LVFU(temp[j][0]) - avgL1) + SQRV(LVFU(temp[j + 1][0]) - avgL1)) + (SQRV(LVFU(temp[j + 2][0]) - avgL1) + SQRV(LVFU(temp[j + 3][0]) - avgL1)) + SQRV(LVFU(temp[j + 4][0]) - avgL1)));
71 const vfloat avgL2 = ((LVFU(temp[j - 4][4]) + LVFU(temp[j - 3][4])) + (LVFU(temp[j - 2][4]) + LVFU(temp[j - 1][4])) + (LVFU(temp[j][4]) + LVFU(temp[j + 1][4])) + (LVFU(temp[j + 2][4]) + LVFU(temp[j + 3][4])) + LVFU(temp[j + 4][4])) / ninev;
72 STVFU(avg[j][4], avgL2);
73 STVFU(dev[j][4], vmaxf(epsv, (SQRV(LVFU(temp[j - 4][4]) - avgL2) + SQRV(LVFU(temp[j - 3][4]) - avgL2)) + (SQRV(LVFU(temp[j - 2][4]) - avgL2) + SQRV(LVFU(temp[j - 1][4]) - avgL2)) + (SQRV(LVFU(temp[j][4]) - avgL2) + SQRV(LVFU(temp[j + 1][4]) - avgL2)) + (SQRV(LVFU(temp[j + 2][4]) - avgL2) + SQRV(LVFU(temp[j + 3][4]) - avgL2)) + SQRV(LVFU(temp[j + 4][4]) - avgL2)));
74 #else
75 #ifdef _OPENMP
76 #pragma omp simd
77 #endif
78 for(int h = 0; h < numCols; ++h) {
79 const float avgL = ((temp[j - 4][h] + temp[j - 3][h]) + (temp[j - 2][h] + temp[j - 1][h]) + (temp[j][h] + temp[j + 1][h]) + (temp[j + 2][h] + temp[j + 3][h]) + temp[j + 4][h]) / 9.f;
80 avg[j][h] = avgL;
81 dev[j][h] = std::max(0.001f, (SQR(temp[j - 4][h] - avgL) + SQR(temp[j - 3][h] - avgL)) + (SQR(temp[j - 2][h] - avgL) + SQR(temp[j - 1][h] - avgL)) + (SQR(temp[j][h] - avgL) + SQR(temp[j + 1][h] - avgL)) + (SQR(temp[j + 2][h] - avgL) + SQR(temp[j + 3][h] - avgL)) + SQR(temp[j + 4][h] - avgL));
82 }
83 #endif
84 }
85
86 for (int j = 5; j < H - 5; j++) {
87 #ifdef _OPENMP
88 #pragma omp simd
89 #endif
90 for(int h = 0; h < numCols; ++h) {
91 const float avgL = avg[j - 1][h];
92 const float avgR = avg[j + 1][h];
93 const float devL = dev[j - 1][h];
94 const float devR = dev[j + 1][h];
95 hpmap[j][k + h] = avgL + (avgR - avgL) * devL / (devL + devR);
96 }
97 }
98 }
99 for (; k < col_to; k++) {
100 for (int i = 5; i < H - 5; i++) {
101 temp[i][0] = std::fabs((rawData[i - 5][k] - rawData[i + 5][k]) - 8 * (rawData[i - 4][k] - rawData[i + 4][k]) + 27 * (rawData[i - 3][k] - rawData[i + 3][k]) - 48 * (rawData[i - 2][k] - rawData[i + 2][k]) + 42 * (rawData[i - 1][k] -rawData[i + 1][k]));
102 }
103
104 for (int j = 4; j < H - 4; j++) {
105 const float avgL = (temp[j - 4][0] + temp[j - 3][0] + temp[j - 2][0] + temp[j - 1][0] + temp[j][0] + temp[j + 1][0] + temp[j + 2][0] + temp[j + 3][0] + temp[j + 4][0]) / 9.f;
106 avg[j][0] = avgL;
107 dev[j][0] = std::max(0.001f, SQR(temp[j - 4][0] - avgL) + SQR(temp[j - 3][0] - avgL) + SQR(temp[j - 2][0] - avgL) + SQR(temp[j - 1][0] - avgL) + SQR(temp[j][0] - avgL) + SQR(temp[j + 1][0] - avgL) + SQR(temp[j + 2][0] - avgL) + SQR(temp[j + 3][0] - avgL) + SQR(temp[j + 4][0] - avgL));
108 }
109
110 for (int j = 5; j < H - 5; j++) {
111 const float avgL = avg[j - 1][0];
112 const float avgR = avg[j + 1][0];
113 const float devL = dev[j - 1][0];
114 const float devR = dev[j + 1][0];
115 hpmap[j][k] = avgL + (avgR - avgL) * devL / (devL + devR);
116 }
117 }
118
119 return RP_NO_ERROR;
120 }
121
hphd_horizontal(const float * const * rawData,float ** hpmap,int row_from,int row_to,int W)122 rpError hphd_horizontal(const float * const *rawData, float** hpmap, int row_from, int row_to, int W)
123 {
124
125 float* temp = new (std::nothrow) float[W] ();
126 float* avg = new (std::nothrow) float[W] ();
127 float* dev = new (std::nothrow) float[W] ();
128
129 rpError rc = RP_NO_ERROR;
130 if(!(temp && avg && dev)) {
131 rc = RP_MEMORY_ERROR;
132 } else {
133
134 #ifdef __SSE2__
135 const vfloat onev = F2V(1.f);
136 const vfloat twov = F2V(2.f);
137 const vfloat zd8v = F2V(0.8f);
138 #endif
139 for (int i = row_from; i < row_to; i++) {
140 #ifdef _OPENMP
141 #pragma omp simd
142 #endif
143 for (int j = 5; j < W - 5; j++) {
144 temp[j] = std::fabs((rawData[i][j - 5] - rawData[i][j + 5]) - 8 * (rawData[i][j - 4] - rawData[i][j + 4]) + 27 * (rawData[i][j - 3] - rawData[i][j + 3]) - 48 * (rawData[i][j - 2] - rawData[i][j + 2]) + 42 * (rawData[i][j - 1] - rawData[i][j + 1]));
145 }
146
147 #ifdef _OPENMP
148 #pragma omp simd
149 #endif
150 for (int j = 4; j < W - 4; j++) {
151 const float avgL = ((temp[j - 4] + temp[j - 3]) + (temp[j - 2] + temp[j - 1]) + (temp[j] + temp[j + 1]) + (temp[j + 2] + temp[j + 3]) + temp[j + 4]) / 9.f;
152 avg[j] = avgL;
153 dev[j] = std::max(0.001f, (SQR(temp[j - 4] - avgL) + SQR(temp[j - 3] - avgL)) + (SQR(temp[j - 2] - avgL) + SQR(temp[j - 1] - avgL)) + (SQR(temp[j] - avgL) + SQR(temp[j + 1] - avgL)) + (SQR(temp[j + 2] - avgL) + SQR(temp[j + 3] - avgL)) + SQR(temp[j + 4] - avgL));
154 }
155
156 int j = 5;
157 #ifdef __SSE2__
158 // faster than #pragma omp simd
159 for (; j < W - 8; j+=4) {
160 const vfloat avgL = LVFU(avg[j - 1]);
161 const vfloat avgR = LVFU(avg[j + 1]);
162 const vfloat devL = LVFU(dev[j - 1]);
163 const vfloat devR = LVFU(dev[j + 1]);
164 const vfloat hpv = avgL + (avgR - avgL) * devL / (devL + devR);
165
166 const vfloat hpmapoldv = LVFU(hpmap[i][j]);
167 const vfloat hpmapv = vselfzero(vmaskf_lt(hpmapoldv, zd8v * hpv), twov);
168 STVFU(hpmap[i][j], vself(vmaskf_lt(hpv, zd8v * hpmapoldv), onev, hpmapv));
169 }
170 #endif
171 for (; j < W - 5; j++) {
172 const float avgL = avg[j - 1];
173 const float avgR = avg[j + 1];
174 const float devL = dev[j - 1];
175 const float devR = dev[j + 1];
176 const float hpv = avgL + (avgR - avgL) * devL / (devL + devR);
177
178 if (hpmap[i][j] < 0.8f * hpv) {
179 hpmap[i][j] = 2;
180 } else if (hpv < 0.8f * hpmap[i][j]) {
181 hpmap[i][j] = 1;
182 } else {
183 hpmap[i][j] = 0;
184 }
185 }
186 }
187 }
188
189 delete [] temp;
190 delete [] avg;
191 delete [] dev;
192
193 return rc;
194 }
195
hphd_RedGreenBlue(const float * const * rawData,const unsigned cfarray[2][2],const float * const * hpmap,int W,int H,float ** red,float ** green,float ** blue)196 rpError hphd_RedGreenBlue(const float * const *rawData, const unsigned cfarray[2][2], const float * const *hpmap, int W, int H, float **red, float **green, float **blue)
197 {
198
199 rpError rc = RP_NO_ERROR;
200 constexpr float eps = 0.001f;
201 #ifdef _OPENMP
202 #pragma omp parallel
203 #endif
204 {
205 int firstRow = -1;
206 int lastRow = -1;
207 #ifdef _OPENMP
208 // note, static scheduling is important in this implementation
209 #pragma omp for schedule(static)
210 #endif
211 for (int i = 3; i < H - 3; i++) {
212 if (firstRow == -1) {
213 firstRow = i;
214 }
215 lastRow = i;
216 for (int j = 3; j < W - 3; j++) {
217 if (fc(cfarray, i, j) & 1) {
218 green[i][j] = rawData[i][j];
219 } else {
220 if (hpmap[i][j] == 1) {
221 const float g2 = rawData[i][j + 1] - rawData[i][j + 2] * 0.5f;
222 const float g4 = rawData[i][j - 1] - rawData[i][j - 2] * 0.5f;
223
224 const float dx = eps + std::fabs(rawData[i][j + 1] - rawData[i][j - 1]);
225 float d1 = rawData[i][j + 3] - rawData[i][j + 1];
226 float d2 = rawData[i][j + 2] - rawData[i][j];
227 float d3 = rawData[i - 1][j + 2] - rawData[i - 1][j];
228 float d4 = rawData[i + 1][j + 2] - rawData[i + 1][j];
229
230 const float e2 = 1.f / (dx + (std::fabs(d1) + std::fabs(d2)) + (std::fabs(d3) + std::fabs(d4)) * 0.5f);
231
232 d1 = rawData[i][j - 3] - rawData[i][j - 1];
233 d2 = rawData[i][j - 2] - rawData[i][j];
234 d3 = rawData[i - 1][j - 2] - rawData[i - 1][j];
235 d4 = rawData[i + 1][j - 2] - rawData[i + 1][j];
236
237 const float e4 = 1.f / (dx + (std::fabs(d1) + std::fabs(d2)) + (std::fabs(d3) + std::fabs(d4)) * 0.5f);
238
239 green[i][j] = rawData[i][j] * 0.5f + (e2 * g2 + e4 * g4) / (e2 + e4);
240 } else if (hpmap[i][j] == 2) {
241 const float g1 = rawData[i - 1][j] - rawData[i - 2][j] * 0.5f;
242 const float g3 = rawData[i + 1][j] - rawData[i + 2][j] * 0.5f;
243
244 const float dy = eps + std::fabs(rawData[i + 1][j] - rawData[i - 1][j]);
245 float d1 = rawData[i - 1][j] - rawData[i - 3][j];
246 float d2 = rawData[i][j] - rawData[i - 2][j];
247 float d3 = rawData[i][j - 1] - rawData[i - 2][j - 1];
248 float d4 = rawData[i][j + 1] - rawData[i - 2][j + 1];
249
250 const float e1 = 1.f / (dy + (std::fabs(d1) + std::fabs(d2)) + (std::fabs(d3) + std::fabs(d4)) * 0.5f);
251
252 d1 = rawData[i + 1][j] - rawData[i + 3][j];
253 d2 = rawData[i][j] - rawData[i + 2][j];
254 d3 = rawData[i][j - 1] - rawData[i + 2][j - 1];
255 d4 = rawData[i][j + 1] - rawData[i + 2][j + 1];
256
257 const float e3 = 1.f / (dy + (std::fabs(d1) + std::fabs(d2)) + (std::fabs(d3) + std::fabs(d4)) * 0.5f);
258
259 green[i][j] = rawData[i][j] * 0.5f + (e1 * g1 + e3 * g3) / (e1 + e3);
260 } else {
261 const float g1 = rawData[i - 1][j] - rawData[i - 2][j] * 0.5f;
262 const float g2 = rawData[i][j + 1] - rawData[i][j + 2] * 0.5f;
263 const float g3 = rawData[i + 1][j] - rawData[i + 2][j] * 0.5f;
264 const float g4 = rawData[i][j - 1] - rawData[i][j - 2] * 0.5f;
265
266 const float dx = eps + std::fabs(rawData[i][j + 1] - rawData[i][j - 1]);
267 const float dy = eps + std::fabs(rawData[i + 1][j] - rawData[i - 1][j]);
268
269 float d1 = rawData[i - 1][j] - rawData[i - 3][j];
270 float d2 = rawData[i][j] - rawData[i - 2][j];
271 float d3 = rawData[i][j - 1] - rawData[i - 2][j - 1];
272 float d4 = rawData[i][j + 1] - rawData[i - 2][j + 1];
273
274 const float e1 = 1.f / (dy + (std::fabs(d1) + std::fabs(d2)) + (std::fabs(d3) + std::fabs(d4)) * 0.5f);
275
276 d1 = rawData[i][j + 3] - rawData[i][j + 1];
277 d2 = rawData[i][j + 2] - rawData[i][j];
278 d3 = rawData[i - 1][j + 2] - rawData[i - 1][j];
279 d4 = rawData[i + 1][j + 2] - rawData[i + 1][j];
280
281 const float e2 = 1.f / (dx + (std::fabs(d1) + std::fabs(d2)) + (std::fabs(d3) + std::fabs(d4)) * 0.5f);
282
283 d1 = rawData[i + 1][j] - rawData[i + 3][j];
284 d2 = rawData[i][j] - rawData[i + 2][j];
285 d3 = rawData[i][j - 1] - rawData[i + 2][j - 1];
286 d4 = rawData[i][j + 1] - rawData[i + 2][j + 1];
287
288 const float e3 = 1.f / (dy + (std::fabs(d1) + std::fabs(d2)) + (std::fabs(d3) + std::fabs(d4)) * 0.5f);
289
290 d1 = rawData[i][j - 3] - rawData[i][j - 1];
291 d2 = rawData[i][j - 2] - rawData[i][j];
292 d3 = rawData[i - 1][j - 2] - rawData[i - 1][j];
293 d4 = rawData[i + 1][j - 2] - rawData[i + 1][j];
294
295 const float e4 = 1.f / (dx + (std::fabs(d1) + std::fabs(d2)) + (std::fabs(d3) + std::fabs(d4)) * 0.5f);
296
297 green[i][j] = rawData[i][j] * 0.5f + ((e1 * g1 + e2 * g2) + (e3 * g3 + e4 * g4)) / (e1 + e2 + e3 + e4);
298 }
299 }
300 }
301 if (i - 1 > firstRow) {
302 interpolate_row_redblue(rawData, cfarray, red[i - 1], blue[i - 1], green[i - 2], green[i - 1], green[i], i - 1, W);
303 }
304
305 }
306 if (firstRow > 3 && firstRow < H - 4) {
307 interpolate_row_redblue(rawData, cfarray, red[firstRow], blue[firstRow], green[firstRow - 1], green[firstRow], green[firstRow + 1], firstRow, W);
308 }
309
310 if (lastRow > 3 && lastRow < H - 4) {
311 interpolate_row_redblue(rawData, cfarray, red[lastRow], blue[lastRow], green[lastRow - 1], green[lastRow], green[lastRow + 1], lastRow, W);
312 }
313 #ifdef _OPENMP
314 #pragma omp single
315 #endif
316 {
317 // let the first thread, which is out of work, do the border interpolation
318 rc = bayerborder_demosaic(W, H, 4, rawData, red, green, blue, cfarray);
319 }
320 }
321 return rc;
322 }
323
324 }
325
326
hphd_demosaic(int width,int height,const float * const * rawData,float ** red,float ** green,float ** blue,const unsigned cfarray[2][2],const std::function<bool (double)> & setProgCancel)327 rpError hphd_demosaic(int width, int height, const float * const *rawData, float **red, float **green, float **blue, const unsigned cfarray[2][2], const std::function<bool(double)> &setProgCancel)
328 {
329 BENCHFUN
330 if (!validateBayerCfa(3, cfarray)) {
331 return RP_WRONG_CFA;
332 }
333
334 rpError rc = RP_NO_ERROR;
335 const int W = width;
336 const int H = height;
337
338 setProgCancel(0.0);
339
340 JaggedArray<float> hpmap(W, H, true);
341
342 #ifdef _OPENMP
343 #pragma omp parallel
344 {
345 int tid = omp_get_thread_num();
346 int nthreads = omp_get_num_threads();
347 int blk = W / nthreads;
348 rpError errCode;
349 if (tid < nthreads - 1) {
350 errCode = hphd_vertical(rawData, hpmap, tid * blk, (tid + 1) * blk, H);
351 } else {
352 errCode = hphd_vertical(rawData, hpmap, tid * blk, W, H);
353 }
354 #pragma omp critical
355 {
356 if (errCode) {
357 rc = errCode;
358 }
359 }
360 }
361 #else
362 rc = hphd_vertical(rawData, hpmap, 0, W, H);
363 #endif
364
365 if (!rc) {
366 setProgCancel(0.35);
367
368 #ifdef _OPENMP
369 #pragma omp parallel
370 {
371 int tid = omp_get_thread_num();
372 int nthreads = omp_get_num_threads();
373 int blk = H / nthreads;
374 rpError errCode;
375
376 if (tid < nthreads - 1) {
377 errCode = hphd_horizontal(rawData, hpmap, tid * blk, (tid + 1) * blk, W);
378 } else {
379 errCode = hphd_horizontal(rawData, hpmap, tid * blk, H, W);
380 }
381 #pragma omp critical
382 {
383 if (errCode) {
384 rc = errCode;
385 }
386 }
387 }
388 #else
389 rc = hphd_horizontal(rawData, hpmap, 0, H, W);
390 #endif
391 if (!rc) {
392 setProgCancel(0.43);
393 rc = hphd_RedGreenBlue(rawData, cfarray, hpmap, W, H, red, green, blue);
394 }
395 }
396 setProgCancel(1.0);
397 return rc;
398 }
399
400