1 ////////////////////////////////////////////////////////////////
2 //
3 // Fast demosaicing algorithm
4 //
5 // copyright (c) 2008-2010 Emil Martinec <ejmartin@uchicago.edu>
6 //
7 //
8 // code dated: August 26, 2010
9 //
10 // fast_demo.cc is free software: you can redistribute it and/or modify
11 // it under the terms of the GNU General Public License as published by
12 // the Free Software Foundation, either version 3 of the License, or
13 // (at your option) any later version.
14 //
15 // This program is distributed in the hope that it will be useful,
16 // but WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 // GNU General Public License for more details.
19 //
20 // You should have received a copy of the GNU General Public License
21 // along with this program. If not, see <http://www.gnu.org/licenses/>.
22 //
23 ////////////////////////////////////////////////////////////////
24
25 #include <cmath>
26 #include "rawimagesource.h"
27 #include "../rtgui/multilangmgr.h"
28 #include "procparams.h"
29 #include "opthelper.h"
30
31 using namespace std;
32 using namespace rtengine;
33
34 #define TS 224
35 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
36
37 /*
38 LUTf RawImageSource::initInvGrad()
39 {
40 LUTf invGrad (0x10000);
41
42 //set up directional weight function
43 for (int i=0; i<0x10000; i++)
44 invGrad[i] = 1.0/SQR(1.0+i);
45
46 return invGrad;
47 }
48 */
49 #define INVGRAD(i) (16.0f/SQR(4.0f+i))
50 #ifdef __SSE2__
51 #define INVGRADV(i) (c16v*_mm_rcp_ps(SQRV(fourv+i)))
52 #endif
53 //LUTf RawImageSource::invGrad = RawImageSource::initInvGrad();
54
fast_demosaic()55 void RawImageSource::fast_demosaic()
56 {
57
58 double progress = 0.0;
59 const bool plistenerActive = plistener;
60
61 //int winx=0, winy=0;
62 //int winw=W, winh=H;
63
64 if (plistener) {
65 plistener->setProgressStr (Glib::ustring::compose(M("TP_RAW_DMETHOD_PROGRESSBAR"), RAWParams::BayerSensor::getMethodString(RAWParams::BayerSensor::Method::FAST)));
66 plistener->setProgress (progress);
67 }
68
69
70 const int bord = 5;
71
72 float clip_pt = 4 * 65535 * initialGain;
73
74 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
75 #ifdef _OPENMP
76 #pragma omp parallel
77 #endif
78 {
79
80 char *buffer;
81 float *greentile;
82 float *redtile;
83 float *bluetile;
84 #define CLF 1
85 // assign working space
86 buffer = (char *) calloc(3 * sizeof(float) * TS * TS + 3 * CLF * 64 + 63, 1);
87 char *data;
88 data = (char*)( ( uintptr_t(buffer) + uintptr_t(63)) / 64 * 64);
89
90 greentile = (float (*)) data; //pointers to array
91 redtile = (float (*)) ((char*)greentile + sizeof(float) * TS * TS + CLF * 64);
92 bluetile = (float (*)) ((char*)redtile + sizeof(float) * TS * TS + CLF * 64);
93
94 #ifdef _OPENMP
95 #pragma omp sections
96 #endif
97 {
98 #ifdef _OPENMP
99 #pragma omp section
100 #endif
101 {
102
103 //first, interpolate borders using bilinear
104 for (int i = 0; i < H; i++)
105 {
106
107 float sum[6];
108 int imin = max(0, i - 1);
109 int imax = min(i + 2, H);
110
111 for (int j = 0; j < bord; j++) { //first few columns
112 for (int c = 0; c < 6; c++) {
113 sum[c] = 0;
114 }
115
116 int jmin = max(0, j - 1);
117
118 for (int i1 = imin; i1 < imax; i1++)
119 for (int j1 = jmin; j1 < j + 2; j1++) {
120 int c = FC(i1, j1);
121 sum[c] += rawData[i1][j1];
122 sum[c + 3]++;
123 }
124
125 int c = FC(i, j);
126
127 if (c == 1) {
128 red[i][j] = sum[0] / sum[3];
129 green[i][j] = rawData[i][j];
130 blue[i][j] = sum[2] / sum[5];
131 } else {
132 green[i][j] = sum[1] / sum[4];
133
134 if (c == 0) {
135 red[i][j] = rawData[i][j];
136 blue[i][j] = sum[2] / sum[5];
137 } else {
138 red[i][j] = sum[0] / sum[3];
139 blue[i][j] = rawData[i][j];
140 }
141 }
142 }//j
143
144 for (int j = W - bord; j < W; j++) { //last few columns
145 for (int c = 0; c < 6; c++) {
146 sum[c] = 0;
147 }
148
149 int jmax = min(j + 2, W);
150
151 for (int i1 = imin; i1 < imax; i1++)
152 for (int j1 = j - 1; j1 < jmax; j1++) {
153 int c = FC(i1, j1);
154 sum[c] += rawData[i1][j1];
155 sum[c + 3]++;
156 }
157
158 int c = FC(i, j);
159
160 if (c == 1) {
161 red[i][j] = sum[0] / sum[3];
162 green[i][j] = rawData[i][j];
163 blue[i][j] = sum[2] / sum[5];
164 } else {
165 green[i][j] = sum[1] / sum[4];
166
167 if (c == 0) {
168 red[i][j] = rawData[i][j];
169 blue[i][j] = sum[2] / sum[5];
170 } else {
171 red[i][j] = sum[0] / sum[3];
172 blue[i][j] = rawData[i][j];
173 }
174 }
175 }//j
176 }//i
177
178 }
179 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
180 #ifdef _OPENMP
181 #pragma omp section
182 #endif
183 {
184
185 for (int j = bord; j < W - bord; j++)
186 {
187 float sum[6];
188
189 for (int i = 0; i < bord; i++) { //first few rows
190 for (int c = 0; c < 6; c++) {
191 sum[c] = 0;
192 }
193
194 for (int i1 = max(0, i - 1); i1 < i + 2; i1++)
195 for (int j1 = j - 1; j1 < j + 2; j1++) {
196 int c = FC(i1, j1);
197 sum[c] += rawData[i1][j1];
198 sum[c + 3]++;
199 }
200
201 int c = FC(i, j);
202
203 if (c == 1) {
204 red[i][j] = sum[0] / sum[3];
205 green[i][j] = rawData[i][j];
206 blue[i][j] = sum[2] / sum[5];
207 } else {
208 green[i][j] = sum[1] / sum[4];
209
210 if (c == 0) {
211 red[i][j] = rawData[i][j];
212 blue[i][j] = sum[2] / sum[5];
213 } else {
214 red[i][j] = sum[0] / sum[3];
215 blue[i][j] = rawData[i][j];
216 }
217 }
218 }//i
219
220 for (int i = H - bord; i < H; i++) { //last few rows
221 for (int c = 0; c < 6; c++) {
222 sum[c] = 0;
223 }
224
225 for (int i1 = i - 1; i1 < min(i + 2, H); i1++)
226 for (int j1 = j - 1; j1 < j + 2; j1++) {
227 int c = FC(i1, j1);
228 sum[c] += rawData[i1][j1];
229 sum[c + 3]++;
230 }
231
232 int c = FC(i, j);
233
234 if (c == 1) {
235 red[i][j] = sum[0] / sum[3];
236 green[i][j] = rawData[i][j];
237 blue[i][j] = sum[2] / sum[5];
238 } else {
239 green[i][j] = sum[1] / sum[4];
240
241 if (c == 0) {
242 red[i][j] = rawData[i][j];
243 blue[i][j] = sum[2] / sum[5];
244 } else {
245 red[i][j] = sum[0] / sum[3];
246 blue[i][j] = rawData[i][j];
247 }
248 }
249 }//i
250 }//j
251
252 }
253 }
254 #ifdef _OPENMP
255 #pragma omp single
256 #endif
257 {
258 if(plistenerActive) {
259 progress += 0.1;
260 plistener->setProgress(progress);
261 }
262 }
263 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
264 int progressCounter = 0;
265 const double progressInc = 16.0 * (1.0 - progress) / ((H * W) / ((TS - 4) * (TS - 4)));
266
267 #ifdef _OPENMP
268 #pragma omp for nowait
269 #endif
270
271 for (int top = bord - 2; top < H - bord + 2; top += TS - (4))
272 for (int left = bord - 2; left < W - bord + 2; left += TS - (4)) {
273 int bottom = min(top + TS, H - bord + 2);
274 int right = min(left + TS, W - bord + 2);
275
276 #ifdef __SSE2__
277 int j, cc;
278 __m128 wtuv, wtdv, wtlv, wtrv;
279 __m128 greenv, tempv, absv, abs2v;
280 __m128 c16v = _mm_set1_ps( 16.0f );
281 __m128 fourv = _mm_set1_ps( 4.0f );
282 vmask selmask;
283 vmask andmask = _mm_set_epi32( 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff );
284
285 if(FC(top, left) == 1) {
286 selmask = _mm_set_epi32( 0, 0xffffffff, 0, 0xffffffff );
287 } else {
288 selmask = _mm_set_epi32( 0xffffffff, 0, 0xffffffff, 0 );
289 }
290
291 #endif
292
293 // interpolate G using gradient weights
294 for (int i = top, rr = 0; i < bottom; i++, rr++) {
295 float wtu, wtd, wtl, wtr;
296 #ifdef __SSE2__
297 selmask = (vmask)_mm_andnot_ps( (__m128)selmask, (__m128)andmask);
298
299 for (j = left, cc = 0; j < right - 3; j += 4, cc += 4) {
300 tempv = LVFU(rawData[i][j]);
301 absv = vabsf(LVFU(rawData[i - 1][j]) - LVFU(rawData[i + 1][j]));
302 wtuv = INVGRADV(absv + vabsf(tempv - LVFU(rawData[i - 2][j])) + vabsf(LVFU(rawData[i - 1][j]) - LVFU(rawData[i - 3][j])));
303 wtdv = INVGRADV(absv + vabsf(tempv - LVFU(rawData[i + 2][j])) + vabsf(LVFU(rawData[i + 1][j]) - LVFU(rawData[i + 3][j])));
304 abs2v = vabsf(LVFU(rawData[i][j - 1]) - LVFU(rawData[i][j + 1]));
305 wtlv = INVGRADV(abs2v + vabsf(tempv - LVFU(rawData[i][j - 2])) + vabsf(LVFU(rawData[i][j - 1]) - LVFU(rawData[i][j - 3])));
306 wtrv = INVGRADV(abs2v + vabsf(tempv - LVFU(rawData[i][j + 2])) + vabsf(LVFU(rawData[i][j + 1]) - LVFU(rawData[i][j + 3])));
307 greenv = (wtuv * LVFU(rawData[i - 1][j]) + wtdv * LVFU(rawData[i + 1][j]) + wtlv * LVFU(rawData[i][j - 1]) + wtrv * LVFU(rawData[i][j + 1])) / (wtuv + wtdv + wtlv + wtrv);
308 _mm_store_ps(&greentile[rr * TS + cc], vself(selmask, greenv, tempv));
309 _mm_store_ps(&redtile[rr * TS + cc], tempv);
310 _mm_store_ps(&bluetile[rr * TS + cc], tempv);
311 }
312
313 for (; j < right; j++, cc++) {
314
315 if (FC(i, j) == 1) {
316 greentile[rr * TS + cc] = rawData[i][j];
317
318 } else {
319 //compute directional weights using image gradients
320 wtu = INVGRAD((abs(rawData[i + 1][j] - rawData[i - 1][j]) + abs(rawData[i][j] - rawData[i - 2][j]) + abs(rawData[i - 1][j] - rawData[i - 3][j])));
321 wtd = INVGRAD((abs(rawData[i - 1][j] - rawData[i + 1][j]) + abs(rawData[i][j] - rawData[i + 2][j]) + abs(rawData[i + 1][j] - rawData[i + 3][j])));
322 wtl = INVGRAD((abs(rawData[i][j + 1] - rawData[i][j - 1]) + abs(rawData[i][j] - rawData[i][j - 2]) + abs(rawData[i][j - 1] - rawData[i][j - 3])));
323 wtr = INVGRAD((abs(rawData[i][j - 1] - rawData[i][j + 1]) + abs(rawData[i][j] - rawData[i][j + 2]) + abs(rawData[i][j + 1] - rawData[i][j + 3])));
324
325 //store in rgb array the interpolated G value at R/B grid points using directional weighted average
326 greentile[rr * TS + cc] = (wtu * rawData[i - 1][j] + wtd * rawData[i + 1][j] + wtl * rawData[i][j - 1] + wtr * rawData[i][j + 1]) / (wtu + wtd + wtl + wtr);
327 }
328
329 redtile[rr * TS + cc] = rawData[i][j];
330 bluetile[rr * TS + cc] = rawData[i][j];
331 }
332
333 #else
334
335 for (int j = left, cc = 0; j < right; j++, cc++) {
336 if (FC(i, j) == 1) {
337 greentile[rr * TS + cc] = rawData[i][j];
338 } else {
339 //compute directional weights using image gradients
340 wtu = INVGRAD((abs(rawData[i + 1][j] - rawData[i - 1][j]) + abs(rawData[i][j] - rawData[i - 2][j]) + abs(rawData[i - 1][j] - rawData[i - 3][j])));
341 wtd = INVGRAD((abs(rawData[i - 1][j] - rawData[i + 1][j]) + abs(rawData[i][j] - rawData[i + 2][j]) + abs(rawData[i + 1][j] - rawData[i + 3][j])));
342 wtl = INVGRAD((abs(rawData[i][j + 1] - rawData[i][j - 1]) + abs(rawData[i][j] - rawData[i][j - 2]) + abs(rawData[i][j - 1] - rawData[i][j - 3])));
343 wtr = INVGRAD((abs(rawData[i][j - 1] - rawData[i][j + 1]) + abs(rawData[i][j] - rawData[i][j + 2]) + abs(rawData[i][j + 1] - rawData[i][j + 3])));
344
345 //store in rgb array the interpolated G value at R/B grid points using directional weighted average
346 greentile[rr * TS + cc] = (wtu * rawData[i - 1][j] + wtd * rawData[i + 1][j] + wtl * rawData[i][j - 1] + wtr * rawData[i][j + 1]) / (wtu + wtd + wtl + wtr);
347 }
348
349 redtile[rr * TS + cc] = rawData[i][j];
350 bluetile[rr * TS + cc] = rawData[i][j];
351 }
352
353 #endif
354 }
355
356 #ifdef __SSE2__
357 __m128 zd25v = _mm_set1_ps(0.25f);
358 __m128 clip_ptv = _mm_set1_ps( clip_pt );
359 #endif
360
361 for (int i = top + 1, rr = 1; i < bottom - 1; i++, rr++) {
362 if (FC(i, left + (FC(i, 2) & 1) + 1) == 0)
363 #ifdef __SSE2__
364 for (int j = left + 1, cc = 1; j < right - 1; j += 4, cc += 4) {
365 //interpolate B/R colors at R/B sites
366 _mm_storeu_ps(&bluetile[rr * TS + cc], LVFU(greentile[rr * TS + cc]) - zd25v * ((LVFU(greentile[(rr - 1)*TS + (cc - 1)]) + LVFU(greentile[(rr - 1)*TS + (cc + 1)]) + LVFU(greentile[(rr + 1)*TS + cc + 1]) + LVFU(greentile[(rr + 1)*TS + cc - 1])) -
367 vminf(LVFU(rawData[i - 1][j - 1]) + LVFU(rawData[i - 1][j + 1]) + LVFU(rawData[i + 1][j + 1]) + LVFU(rawData[i + 1][j - 1]), clip_ptv)));
368 }
369
370 #else
371
372 for (int cc = (FC(i, 2) & 1) + 1, j = left + cc; j < right - 1; j += 2, cc += 2) {
373 //interpolate B/R colors at R/B sites
374 bluetile[rr * TS + cc] = greentile[rr * TS + cc] - 0.25f * ((greentile[(rr - 1) * TS + (cc - 1)] + greentile[(rr - 1) * TS + (cc + 1)] + greentile[(rr + 1) * TS + cc + 1] + greentile[(rr + 1) * TS + cc - 1]) -
375 min(clip_pt, rawData[i - 1][j - 1] + rawData[i - 1][j + 1] + rawData[i + 1][j + 1] + rawData[i + 1][j - 1]));
376 }
377
378 #endif
379 else
380 #ifdef __SSE2__
381 for (int j = left + 1, cc = 1; j < right - 1; j += 4, cc += 4) {
382 //interpolate B/R colors at R/B sites
383 _mm_storeu_ps(&redtile[rr * TS + cc], LVFU(greentile[rr * TS + cc]) - zd25v * ((LVFU(greentile[(rr - 1)*TS + cc - 1]) + LVFU(greentile[(rr - 1)*TS + cc + 1]) + LVFU(greentile[(rr + 1)*TS + cc + 1]) + LVFU(greentile[(rr + 1)*TS + cc - 1])) -
384 vminf(LVFU(rawData[i - 1][j - 1]) + LVFU(rawData[i - 1][j + 1]) + LVFU(rawData[i + 1][j + 1]) + LVFU(rawData[i + 1][j - 1]), clip_ptv)));
385 }
386
387 #else
388
389 for (int cc = (FC(i, 2) & 1) + 1, j = left + cc; j < right - 1; j += 2, cc += 2) {
390 //interpolate B/R colors at R/B sites
391 redtile[rr * TS + cc] = greentile[rr * TS + cc] - 0.25f * ((greentile[(rr - 1) * TS + cc - 1] + greentile[(rr - 1) * TS + cc + 1] + greentile[(rr + 1) * TS + cc + 1] + greentile[(rr + 1) * TS + cc - 1]) -
392 min(clip_pt, rawData[i - 1][j - 1] + rawData[i - 1][j + 1] + rawData[i + 1][j + 1] + rawData[i + 1][j - 1]));
393 }
394
395 #endif
396 }
397
398
399 #ifdef __SSE2__
400 __m128 temp1v, temp2v, greensumv;
401 selmask = _mm_set_epi32( 0xffffffff, 0, 0xffffffff, 0 );
402 #endif
403
404 // interpolate R/B using color differences
405 for (int i = top + 2, rr = 2; i < bottom - 2; i++, rr++) {
406 #ifdef __SSE2__
407
408 for (int cc = 2 + (FC(i, 2) & 1), j = left + cc; j < right - 2; j += 4, cc += 4) {
409 // no need to take care about the borders of the tile. There's enough free space.
410 //interpolate R and B colors at G sites
411 greenv = LVFU(greentile[rr * TS + cc]);
412 greensumv = LVFU(greentile[(rr - 1) * TS + cc]) + LVFU(greentile[(rr + 1) * TS + cc]) + LVFU(greentile[rr * TS + cc - 1]) + LVFU(greentile[rr * TS + cc + 1]);
413
414 temp1v = LVFU(redtile[rr * TS + cc]);
415 temp2v = greenv - zd25v * (greensumv - LVFU(redtile[(rr - 1) * TS + cc]) - LVFU(redtile[(rr + 1) * TS + cc]) - LVFU(redtile[rr * TS + cc - 1]) - LVFU(redtile[rr * TS + cc + 1]));
416
417 // temp2v = greenv - zd25v*((LVFU(greentile[(rr-1)*TS+cc])-LVFU(redtile[(rr-1)*TS+cc]))+(LVFU(greentile[(rr+1)*TS+cc])-LVFU(redtile[(rr+1)*TS+cc]))+
418 // (LVFU(greentile[rr*TS+cc-1])-LVFU(redtile[rr*TS+cc-1]))+(LVFU(greentile[rr*TS+cc+1])-LVFU(redtile[rr*TS+cc+1])));
419 _mm_storeu_ps( &redtile[rr * TS + cc], vself(selmask, temp1v, temp2v));
420
421 temp1v = LVFU(bluetile[rr * TS + cc]);
422
423 temp2v = greenv - zd25v * (greensumv - LVFU(bluetile[(rr - 1) * TS + cc]) - LVFU(bluetile[(rr + 1) * TS + cc]) - LVFU(bluetile[rr * TS + cc - 1]) - LVFU(bluetile[rr * TS + cc + 1]));
424
425 // temp2v = greenv - zd25v*((LVFU(greentile[(rr-1)*TS+cc])-LVFU(bluetile[(rr-1)*TS+cc]))+(LVFU(greentile[(rr+1)*TS+cc])-LVFU(bluetile[(rr+1)*TS+cc]))+
426 // (LVFU(greentile[rr*TS+cc-1])-LVFU(bluetile[rr*TS+cc-1]))+(LVFU(greentile[rr*TS+cc+1])-LVFU(bluetile[rr*TS+cc+1])));
427 _mm_storeu_ps( &bluetile[rr * TS + cc], vself(selmask, temp1v, temp2v));
428 }
429
430 #else
431
432 for (int cc = 2 + (FC(i, 2) & 1), j = left + cc; j < right - 2; j += 2, cc += 2) {
433 //interpolate R and B colors at G sites
434 redtile[rr * TS + cc] = greentile[rr * TS + cc] - 0.25f * ((greentile[(rr - 1) * TS + cc] - redtile[(rr - 1) * TS + cc]) + (greentile[(rr + 1) * TS + cc] - redtile[(rr + 1) * TS + cc]) +
435 (greentile[rr * TS + cc - 1] - redtile[rr * TS + cc - 1]) + (greentile[rr * TS + cc + 1] - redtile[rr * TS + cc + 1]));
436 bluetile[rr * TS + cc] = greentile[rr * TS + cc] - 0.25f * ((greentile[(rr - 1) * TS + cc] - bluetile[(rr - 1) * TS + cc]) + (greentile[(rr + 1) * TS + cc] - bluetile[(rr + 1) * TS + cc]) +
437 (greentile[rr * TS + cc - 1] - bluetile[rr * TS + cc - 1]) + (greentile[rr * TS + cc + 1] - bluetile[rr * TS + cc + 1]));
438 }
439
440 #endif
441 }
442
443
444 for (int i = top + 2, rr = 2; i < bottom - 2; i++, rr++) {
445 #ifdef __SSE2__
446
447 for (j = left + 2, cc = 2; j < right - 5; j += 4, cc += 4) {
448 _mm_storeu_ps(&red[i][j], vmaxf(LVFU(redtile[rr * TS + cc]), ZEROV));
449 _mm_storeu_ps(&green[i][j], vmaxf(LVFU(greentile[rr * TS + cc]), ZEROV));
450 _mm_storeu_ps(&blue[i][j], vmaxf(LVFU(bluetile[rr * TS + cc]), ZEROV));
451 }
452
453 for (; j < right - 2; j++, cc++) {
454 red[i][j] = std::max(0.f, redtile[rr * TS + cc]);
455 green[i][j] = std::max(0.f, greentile[rr * TS + cc]);
456 blue[i][j] = std::max(0.f, bluetile[rr * TS + cc]);
457 }
458
459 #else
460
461 for (int j = left + 2, cc = 2; j < right - 2; j++, cc++) {
462 red[i][j] = std::max(0.f, redtile[rr * TS + cc]);
463 green[i][j] = std::max(0.f, greentile[rr * TS + cc]);
464 blue[i][j] = std::max(0.f, bluetile[rr * TS + cc]);
465 }
466
467 #endif
468
469
470 }
471
472 if(plistenerActive && ((++progressCounter) % 16 == 0)) {
473 #ifdef _OPENMP
474 #pragma omp critical (updateprogress)
475 #endif
476 {
477 progress += progressInc;
478 progress = min(1.0, progress);
479 plistener->setProgress (progress);
480 }
481 }
482
483 }
484
485 free(buffer);
486 } // End of parallelization
487
488 if(plistenerActive) {
489 plistener->setProgress(1.00);
490 }
491
492
493
494 }
495 #undef TS
496 #undef CLF
497