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 "rawimagesource.h"
22 #include "rawimagesource_i.h"
23 #include "jaggedarray.h"
24 #include "rawimage.h"
25 #include "rt_math.h"
26 #include "../rtgui/multilangmgr.h"
27 #include "opthelper.h"
28 //#define BENCHMARK
29 #include "StopWatch.h"
30 #ifdef _OPENMP
31 #include <omp.h>
32 #endif
33
34 using namespace std;
35 using namespace rtengine;
36
37 namespace {
38
hphd_vertical(const array2D<float> & rawData,float ** hpmap,int col_from,int col_to,int H)39 void hphd_vertical(const array2D<float> &rawData, float** hpmap, int col_from, int col_to, int H)
40 {
41
42 // process 'numCols' columns for better usage of L1 cpu cache (especially faster for large values of H)
43 constexpr int numCols = 8;
44 JaggedArray<float> temp(numCols, H, true);
45 JaggedArray<float> avg(numCols, H, true);
46 JaggedArray<float> dev(numCols, H, true);
47
48 int k = col_from;
49 #ifdef __SSE2__
50 const vfloat ninev = F2V(9.f);
51 const vfloat epsv = F2V(0.001f);
52 #endif
53 for (; k < col_to - 7; k += numCols) {
54 for (int i = 5; i < H - 5; i++) {
55 #ifdef _OPENMP
56 #pragma omp simd
57 #endif
58 for(int h = 0; h < numCols; ++h) {
59 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]));
60 }
61 }
62
63 for (int j = 4; j < H - 4; j++) {
64 #ifdef __SSE2__
65 // faster than #pragma omp simd...
66 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;
67 STVFU(avg[j][0], avgL1);
68 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)));
69 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;
70 STVFU(avg[j][4], avgL2);
71 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)));
72 #else
73 #ifdef _OPENMP
74 #pragma omp simd
75 #endif
76 for(int h = 0; h < numCols; ++h) {
77 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;
78 avg[j][h] = avgL;
79 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));
80 }
81 #endif
82 }
83
84 for (int j = 5; j < H - 5; j++) {
85 #ifdef _OPENMP
86 #pragma omp simd
87 #endif
88 for(int h = 0; h < numCols; ++h) {
89 const float avgL = avg[j - 1][h];
90 const float avgR = avg[j + 1][h];
91 const float devL = dev[j - 1][h];
92 const float devR = dev[j + 1][h];
93 hpmap[j][k + h] = avgL + (avgR - avgL) * devL / (devL + devR);
94 }
95 }
96 }
97 for (; k < col_to; k++) {
98 for (int i = 5; i < H - 5; i++) {
99 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]));
100 }
101
102 for (int j = 4; j < H - 4; j++) {
103 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;
104 avg[j][0] = avgL;
105 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));
106 }
107
108 for (int j = 5; j < H - 5; j++) {
109 const float avgL = avg[j - 1][0];
110 const float avgR = avg[j + 1][0];
111 const float devL = dev[j - 1][0];
112 const float devR = dev[j + 1][0];
113 hpmap[j][k] = avgL + (avgR - avgL) * devL / (devL + devR);
114 }
115 }
116 }
117
hphd_horizontal(const array2D<float> & rawData,float ** hpmap,int row_from,int row_to,int W)118 void hphd_horizontal(const array2D<float> &rawData, float** hpmap, int row_from, int row_to, int W)
119 {
120
121 float* temp = new float[W];
122 float* avg = new float[W];
123 float* dev = new float[W];
124
125 memset(temp, 0, W * sizeof(float));
126 memset(avg, 0, W * sizeof(float));
127 memset(dev, 0, W * sizeof(float));
128
129 #ifdef __SSE2__
130 const vfloat onev = F2V(1.f);
131 const vfloat twov = F2V(2.f);
132 const vfloat zd8v = F2V(0.8f);
133 #endif
134 for (int i = row_from; i < row_to; i++) {
135 #ifdef _OPENMP
136 #pragma omp simd
137 #endif
138 for (int j = 5; j < W - 5; j++) {
139 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]));
140 }
141
142 #ifdef _OPENMP
143 #pragma omp simd
144 #endif
145 for (int j = 4; j < W - 4; j++) {
146 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;
147 avg[j] = avgL;
148 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));
149 }
150
151 int j = 5;
152 #ifdef __SSE2__
153 // faster than #pragma omp simd
154 for (; j < W - 8; j+=4) {
155 const vfloat avgL = LVFU(avg[j - 1]);
156 const vfloat avgR = LVFU(avg[j + 1]);
157 const vfloat devL = LVFU(dev[j - 1]);
158 const vfloat devR = LVFU(dev[j + 1]);
159 const vfloat hpv = avgL + (avgR - avgL) * devL / (devL + devR);
160
161 const vfloat hpmapoldv = LVFU(hpmap[i][j]);
162 const vfloat hpmapv = vselfzero(vmaskf_lt(hpmapoldv, zd8v * hpv), twov);
163 STVFU(hpmap[i][j], vself(vmaskf_lt(hpv, zd8v * hpmapoldv), onev, hpmapv));
164 }
165 #endif
166 for (; j < W - 5; j++) {
167 const float avgL = avg[j - 1];
168 const float avgR = avg[j + 1];
169 const float devL = dev[j - 1];
170 const float devR = dev[j + 1];
171 const float hpv = avgL + (avgR - avgL) * devL / (devL + devR);
172
173 if (hpmap[i][j] < 0.8f * hpv) {
174 hpmap[i][j] = 2;
175 } else if (hpv < 0.8f * hpmap[i][j]) {
176 hpmap[i][j] = 1;
177 } else {
178 hpmap[i][j] = 0;
179 }
180 }
181 }
182
183 delete [] temp;
184 delete [] avg;
185 delete [] dev;
186 }
187
hphd_green(const RawImage * ri,const array2D<float> & rawData,float ** hpmap,int W,int H,array2D<float> & green)188 void hphd_green(const RawImage *ri, const array2D<float> &rawData, float** hpmap, int W, int H, array2D<float> &green)
189 {
190
191 constexpr float eps = 0.001f;
192 #ifdef _OPENMP
193 #pragma omp parallel for schedule(dynamic, 16)
194 #endif
195
196 for (int i = 3; i < H - 3; i++) {
197 for (int j = 3; j < W - 3; j++) {
198 if (ri->ISGREEN(i, j)) {
199 green[i][j] = rawData[i][j];
200 } else {
201 if (hpmap[i][j] == 1) {
202 const float g2 = rawData[i][j + 1] - rawData[i][j + 2] * 0.5f;
203 const float g4 = rawData[i][j - 1] - rawData[i][j - 2] * 0.5f;
204
205 const float dx = eps + std::fabs(rawData[i][j + 1] - rawData[i][j - 1]);
206 float d1 = rawData[i][j + 3] - rawData[i][j + 1];
207 float d2 = rawData[i][j + 2] - rawData[i][j];
208 float d3 = rawData[i - 1][j + 2] - rawData[i - 1][j];
209 float d4 = rawData[i + 1][j + 2] - rawData[i + 1][j];
210
211 const float e2 = 1.f / (dx + (std::fabs(d1) + std::fabs(d2)) + (std::fabs(d3) + std::fabs(d4)) * 0.5f);
212
213 d1 = rawData[i][j - 3] - rawData[i][j - 1];
214 d2 = rawData[i][j - 2] - rawData[i][j];
215 d3 = rawData[i - 1][j - 2] - rawData[i - 1][j];
216 d4 = rawData[i + 1][j - 2] - rawData[i + 1][j];
217
218 const float e4 = 1.f / (dx + (std::fabs(d1) + std::fabs(d2)) + (std::fabs(d3) + std::fabs(d4)) * 0.5f);
219
220 green[i][j] = std::max(0.f, rawData[i][j] * 0.5f + (e2 * g2 + e4 * g4) / (e2 + e4));
221 } else if (hpmap[i][j] == 2) {
222 const float g1 = rawData[i - 1][j] - rawData[i - 2][j] * 0.5f;
223 const float g3 = rawData[i + 1][j] - rawData[i + 2][j] * 0.5f;
224
225 const float dy = eps + std::fabs(rawData[i + 1][j] - rawData[i - 1][j]);
226 float d1 = rawData[i - 1][j] - rawData[i - 3][j];
227 float d2 = rawData[i][j] - rawData[i - 2][j];
228 float d3 = rawData[i][j - 1] - rawData[i - 2][j - 1];
229 float d4 = rawData[i][j + 1] - rawData[i - 2][j + 1];
230
231 const float e1 = 1.f / (dy + (std::fabs(d1) + std::fabs(d2)) + (std::fabs(d3) + std::fabs(d4)) * 0.5f);
232
233 d1 = rawData[i + 1][j] - rawData[i + 3][j];
234 d2 = rawData[i][j] - rawData[i + 2][j];
235 d3 = rawData[i][j - 1] - rawData[i + 2][j - 1];
236 d4 = rawData[i][j + 1] - rawData[i + 2][j + 1];
237
238 const float e3 = 1.f / (dy + (std::fabs(d1) + std::fabs(d2)) + (std::fabs(d3) + std::fabs(d4)) * 0.5f);
239
240 green[i][j] = std::max(0.f, rawData[i][j] * 0.5f + (e1 * g1 + e3 * g3) / (e1 + e3));
241 } else {
242 const float g1 = rawData[i - 1][j] - rawData[i - 2][j] * 0.5f;
243 const float g2 = rawData[i][j + 1] - rawData[i][j + 2] * 0.5f;
244 const float g3 = rawData[i + 1][j] - rawData[i + 2][j] * 0.5f;
245 const float g4 = rawData[i][j - 1] - rawData[i][j - 2] * 0.5f;
246
247 const float dx = eps + std::fabs(rawData[i][j + 1] - rawData[i][j - 1]);
248 const float dy = eps + std::fabs(rawData[i + 1][j] - rawData[i - 1][j]);
249
250 float d1 = rawData[i - 1][j] - rawData[i - 3][j];
251 float d2 = rawData[i][j] - rawData[i - 2][j];
252 float d3 = rawData[i][j - 1] - rawData[i - 2][j - 1];
253 float d4 = rawData[i][j + 1] - rawData[i - 2][j + 1];
254
255 const float e1 = 1.f / (dy + (std::fabs(d1) + std::fabs(d2)) + (std::fabs(d3) + std::fabs(d4)) * 0.5f);
256
257 d1 = rawData[i][j + 3] - rawData[i][j + 1];
258 d2 = rawData[i][j + 2] - rawData[i][j];
259 d3 = rawData[i - 1][j + 2] - rawData[i - 1][j];
260 d4 = rawData[i + 1][j + 2] - rawData[i + 1][j];
261
262 const float e2 = 1.f / (dx + (std::fabs(d1) + std::fabs(d2)) + (std::fabs(d3) + std::fabs(d4)) * 0.5f);
263
264 d1 = rawData[i + 1][j] - rawData[i + 3][j];
265 d2 = rawData[i][j] - rawData[i + 2][j];
266 d3 = rawData[i][j - 1] - rawData[i + 2][j - 1];
267 d4 = rawData[i][j + 1] - rawData[i + 2][j + 1];
268
269 const float e3 = 1.f / (dy + (std::fabs(d1) + std::fabs(d2)) + (std::fabs(d3) + std::fabs(d4)) * 0.5f);
270
271 d1 = rawData[i][j - 3] - rawData[i][j - 1];
272 d2 = rawData[i][j - 2] - rawData[i][j];
273 d3 = rawData[i - 1][j - 2] - rawData[i - 1][j];
274 d4 = rawData[i + 1][j - 2] - rawData[i + 1][j];
275
276 const float e4 = 1.f / (dx + (std::fabs(d1) + std::fabs(d2)) + (std::fabs(d3) + std::fabs(d4)) * 0.5f);
277
278 green[i][j] = std::max(0.f, rawData[i][j] * 0.5f + ((e1 * g1 + e2 * g2) + (e3 * g3 + e4 * g4)) / (e1 + e2 + e3 + e4));
279 }
280 }
281 }
282 }
283 }
284
285 }
286
287 namespace rtengine
288 {
289
hphd_demosaic()290 void RawImageSource::hphd_demosaic ()
291 {
292 BENCHFUN
293 if (plistener) {
294 plistener->setProgressStr(Glib::ustring::compose(M("TP_RAW_DMETHOD_PROGRESSBAR"), RAWParams::BayerSensor::getMethodString(RAWParams::BayerSensor::Method::HPHD)));
295 plistener->setProgress(0.0);
296 }
297
298 JaggedArray<float> hpmap(W, H, true);
299
300 #ifdef _OPENMP
301 #pragma omp parallel
302 {
303 int tid = omp_get_thread_num();
304 int nthreads = omp_get_num_threads();
305 int blk = W / nthreads;
306
307 if (tid < nthreads - 1) {
308 hphd_vertical(rawData, hpmap, tid * blk, (tid + 1)*blk, H);
309 } else {
310 hphd_vertical(rawData, hpmap, tid * blk, W, H);
311 }
312 }
313 #else
314 hphd_vertical(rawData, hpmap, 0, W, H);
315 #endif
316
317 if (plistener) {
318 plistener->setProgress(0.35);
319 }
320
321 #ifdef _OPENMP
322 #pragma omp parallel
323 {
324 int tid = omp_get_thread_num();
325 int nthreads = omp_get_num_threads();
326 int blk = H / nthreads;
327
328 if (tid < nthreads - 1) {
329 hphd_horizontal(rawData, hpmap, tid * blk, (tid + 1)*blk, W);
330 } else {
331 hphd_horizontal(rawData, hpmap, tid * blk, H, W);
332 }
333 }
334 #else
335 hphd_horizontal(rawData, hpmap, 0, H, W);
336 #endif
337
338 if (plistener) {
339 plistener->setProgress(0.43);
340 }
341
342 hphd_green(ri, rawData, hpmap, W, H, green);
343
344 if (plistener) {
345 plistener->setProgress(0.65);
346 }
347
348 #ifdef _OPENMP
349 #pragma omp parallel for
350 #endif
351 for (int i = 4; i < H - 4; i++) {
352 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);
353 }
354
355 border_interpolate2(W, H, 4, rawData, red, green, blue);
356
357 if (plistener) {
358 plistener->setProgress(1.0);
359 }
360 }
361
362
363
364 } /* namespace */
365