1 /* -*- C++ -*-
2 * Copyright 2019-2020 LibRaw LLC (info@libraw.org)
3 *
4 LibRaw uses code from dcraw.c -- Dave Coffin's raw photo decoder,
5 dcraw.c is copyright 1997-2018 by Dave Coffin, dcoffin a cybercom o net.
6 LibRaw do not use RESTRICTED code from dcraw.c
7
8 LibRaw is free software; you can redistribute it and/or modify
9 it under the terms of the one of two licenses as you choose:
10
11 1. GNU LESSER GENERAL PUBLIC LICENSE version 2.1
12 (See file LICENSE.LGPL provided in LibRaw distribution archive for details).
13
14 2. COMMON DEVELOPMENT AND DISTRIBUTION LICENSE (CDDL) Version 1.0
15 (See file LICENSE.CDDL provided in LibRaw distribution archive for details).
16
17 */
18
19 #include "../../internal/dcraw_defs.h"
20
hat_transform(float * temp,float * base,int st,int size,int sc)21 void LibRaw::hat_transform(float *temp, float *base, int st, int size, int sc)
22 {
23 int i;
24 for (i = 0; i < sc; i++)
25 temp[i] = 2 * base[st * i] + base[st * (sc - i)] + base[st * (i + sc)];
26 for (; i + sc < size; i++)
27 temp[i] = 2 * base[st * i] + base[st * (i - sc)] + base[st * (i + sc)];
28 for (; i < size; i++)
29 temp[i] = 2 * base[st * i] + base[st * (i - sc)] +
30 base[st * (2 * size - 2 - (i + sc))];
31 }
32
33 #if !defined(LIBRAW_USE_OPENMP)
wavelet_denoise()34 void LibRaw::wavelet_denoise()
35 {
36 float *fimg = 0, *temp, thold, mul[2], avg, diff;
37 int scale = 1, size, lev, hpass, lpass, row, col, nc, c, i, wlast, blk[2];
38 ushort *window[4];
39 static const float noise[] = {0.8002, 0.2735, 0.1202, 0.0585,
40 0.0291, 0.0152, 0.0080, 0.0044};
41
42 while (maximum << scale < 0x10000)
43 scale++;
44 maximum <<= --scale;
45 black <<= scale;
46 FORC4 cblack[c] <<= scale;
47 if ((size = iheight * iwidth) < 0x15550000)
48 fimg = (float *)malloc((size * 3 + iheight + iwidth) * sizeof *fimg);
49 merror(fimg, "wavelet_denoise()");
50 temp = fimg + size * 3;
51 if ((nc = colors) == 3 && filters)
52 nc++;
53 FORC(nc)
54 { /* denoise R,G1,B,G3 individually */
55 for (i = 0; i < size; i++)
56 fimg[i] = 256 * sqrt((double)(image[i][c] << scale));
57 for (hpass = lev = 0; lev < 5; lev++)
58 {
59 lpass = size * ((lev & 1) + 1);
60 for (row = 0; row < iheight; row++)
61 {
62 hat_transform(temp, fimg + hpass + row * iwidth, 1, iwidth, 1 << lev);
63 for (col = 0; col < iwidth; col++)
64 fimg[lpass + row * iwidth + col] = temp[col] * 0.25;
65 }
66 for (col = 0; col < iwidth; col++)
67 {
68 hat_transform(temp, fimg + lpass + col, iwidth, iheight, 1 << lev);
69 for (row = 0; row < iheight; row++)
70 fimg[lpass + row * iwidth + col] = temp[row] * 0.25;
71 }
72 thold = threshold * noise[lev];
73 for (i = 0; i < size; i++)
74 {
75 fimg[hpass + i] -= fimg[lpass + i];
76 if (fimg[hpass + i] < -thold)
77 fimg[hpass + i] += thold;
78 else if (fimg[hpass + i] > thold)
79 fimg[hpass + i] -= thold;
80 else
81 fimg[hpass + i] = 0;
82 if (hpass)
83 fimg[i] += fimg[hpass + i];
84 }
85 hpass = lpass;
86 }
87 for (i = 0; i < size; i++)
88 image[i][c] = CLIP(SQR(fimg[i] + fimg[lpass + i]) / 0x10000);
89 }
90 if (filters && colors == 3)
91 { /* pull G1 and G3 closer together */
92 for (row = 0; row < 2; row++)
93 {
94 mul[row] = 0.125 * pre_mul[FC(row + 1, 0) | 1] / pre_mul[FC(row, 0) | 1];
95 blk[row] = cblack[FC(row, 0) | 1];
96 }
97 for (i = 0; i < 4; i++)
98 window[i] = (ushort *)fimg + width * i;
99 for (wlast = -1, row = 1; row < height - 1; row++)
100 {
101 while (wlast < row + 1)
102 {
103 for (wlast++, i = 0; i < 4; i++)
104 window[(i + 3) & 3] = window[i];
105 for (col = FC(wlast, 1) & 1; col < width; col += 2)
106 window[2][col] = BAYER(wlast, col);
107 }
108 thold = threshold / 512;
109 for (col = (FC(row, 0) & 1) + 1; col < width - 1; col += 2)
110 {
111 avg = (window[0][col - 1] + window[0][col + 1] + window[2][col - 1] +
112 window[2][col + 1] - blk[~row & 1] * 4) *
113 mul[row & 1] +
114 (window[1][col] + blk[row & 1]) * 0.5;
115 avg = avg < 0 ? 0 : sqrt(avg);
116 diff = sqrt((double)BAYER(row, col)) - avg;
117 if (diff < -thold)
118 diff += thold;
119 else if (diff > thold)
120 diff -= thold;
121 else
122 diff = 0;
123 BAYER(row, col) = CLIP(SQR(avg + diff) + 0.5);
124 }
125 }
126 }
127 free(fimg);
128 }
129 #else /* LIBRAW_USE_OPENMP */
wavelet_denoise()130 void LibRaw::wavelet_denoise()
131 {
132 float *fimg = 0, *temp, thold, mul[2], avg, diff;
133 int scale = 1, size, lev, hpass, lpass, row, col, nc, c, i, wlast, blk[2];
134 ushort *window[4];
135 static const float noise[] = {0.8002, 0.2735, 0.1202, 0.0585,
136 0.0291, 0.0152, 0.0080, 0.0044};
137
138 while (maximum << scale < 0x10000)
139 scale++;
140 maximum <<= --scale;
141 black <<= scale;
142 FORC4 cblack[c] <<= scale;
143 if ((size = iheight * iwidth) < 0x15550000)
144 fimg = (float *)malloc((size * 3 + iheight + iwidth) * sizeof *fimg);
145 merror(fimg, "wavelet_denoise()");
146 temp = fimg + size * 3;
147 if ((nc = colors) == 3 && filters)
148 nc++;
149 #pragma omp parallel default(shared) private( \
150 i, col, row, thold, lev, lpass, hpass, temp, c) firstprivate(scale, size)
151 {
152 #pragma omp critical /* LibRaw's malloc is not local thread-safe */
153 temp = (float *)malloc((iheight + iwidth) * sizeof *fimg);
154 FORC(nc)
155 { /* denoise R,G1,B,G3 individually */
156 #pragma omp for
157 for (i = 0; i < size; i++)
158 fimg[i] = 256 * sqrt((double)(image[i][c] << scale));
159 for (hpass = lev = 0; lev < 5; lev++)
160 {
161 lpass = size * ((lev & 1) + 1);
162 #pragma omp for
163 for (row = 0; row < iheight; row++)
164 {
165 hat_transform(temp, fimg + hpass + row * iwidth, 1, iwidth, 1 << lev);
166 for (col = 0; col < iwidth; col++)
167 fimg[lpass + row * iwidth + col] = temp[col] * 0.25;
168 }
169 #pragma omp for
170 for (col = 0; col < iwidth; col++)
171 {
172 hat_transform(temp, fimg + lpass + col, iwidth, iheight, 1 << lev);
173 for (row = 0; row < iheight; row++)
174 fimg[lpass + row * iwidth + col] = temp[row] * 0.25;
175 }
176 thold = threshold * noise[lev];
177 #pragma omp for
178 for (i = 0; i < size; i++)
179 {
180 fimg[hpass + i] -= fimg[lpass + i];
181 if (fimg[hpass + i] < -thold)
182 fimg[hpass + i] += thold;
183 else if (fimg[hpass + i] > thold)
184 fimg[hpass + i] -= thold;
185 else
186 fimg[hpass + i] = 0;
187 if (hpass)
188 fimg[i] += fimg[hpass + i];
189 }
190 hpass = lpass;
191 }
192 #pragma omp for
193 for (i = 0; i < size; i++)
194 image[i][c] = CLIP(SQR(fimg[i] + fimg[lpass + i]) / 0x10000);
195 }
196 #pragma omp critical
197 free(temp);
198 } /* end omp parallel */
199 /* the following loops are hard to parallize, no idea yes,
200 * problem is wlast which is carrying dependency
201 * second part should be easyer, but did not yet get it right.
202 */
203 if (filters && colors == 3)
204 { /* pull G1 and G3 closer together */
205 for (row = 0; row < 2; row++)
206 {
207 mul[row] = 0.125 * pre_mul[FC(row + 1, 0) | 1] / pre_mul[FC(row, 0) | 1];
208 blk[row] = cblack[FC(row, 0) | 1];
209 }
210 for (i = 0; i < 4; i++)
211 window[i] = (ushort *)fimg + width * i;
212 for (wlast = -1, row = 1; row < height - 1; row++)
213 {
214 while (wlast < row + 1)
215 {
216 for (wlast++, i = 0; i < 4; i++)
217 window[(i + 3) & 3] = window[i];
218 for (col = FC(wlast, 1) & 1; col < width; col += 2)
219 window[2][col] = BAYER(wlast, col);
220 }
221 thold = threshold / 512;
222 for (col = (FC(row, 0) & 1) + 1; col < width - 1; col += 2)
223 {
224 avg = (window[0][col - 1] + window[0][col + 1] + window[2][col - 1] +
225 window[2][col + 1] - blk[~row & 1] * 4) *
226 mul[row & 1] +
227 (window[1][col] + blk[row & 1]) * 0.5;
228 avg = avg < 0 ? 0 : sqrt(avg);
229 diff = sqrt((double)BAYER(row, col)) - avg;
230 if (diff < -thold)
231 diff += thold;
232 else if (diff > thold)
233 diff -= thold;
234 else
235 diff = 0;
236 BAYER(row, col) = CLIP(SQR(avg + diff) + 0.5);
237 }
238 }
239 }
240 free(fimg);
241 }
242
243 #endif
median_filter()244 void LibRaw::median_filter()
245 {
246 ushort(*pix)[4];
247 int pass, c, i, j, k, med[9];
248 static const uchar opt[] = /* Optimal 9-element median search */
249 {1, 2, 4, 5, 7, 8, 0, 1, 3, 4, 6, 7, 1, 2, 4, 5, 7, 8, 0,
250 3, 5, 8, 4, 7, 3, 6, 1, 4, 2, 5, 4, 7, 4, 2, 6, 4, 4, 2};
251
252 for (pass = 1; pass <= med_passes; pass++)
253 {
254 RUN_CALLBACK(LIBRAW_PROGRESS_MEDIAN_FILTER, pass - 1, med_passes);
255 for (c = 0; c < 3; c += 2)
256 {
257 for (pix = image; pix < image + width * height; pix++)
258 pix[0][3] = pix[0][c];
259 for (pix = image + width; pix < image + width * (height - 1); pix++)
260 {
261 if ((pix - image + 1) % width < 2)
262 continue;
263 for (k = 0, i = -width; i <= width; i += width)
264 for (j = i - 1; j <= i + 1; j++)
265 med[k++] = pix[j][3] - pix[j][1];
266 for (i = 0; i < int(sizeof opt); i += 2)
267 if (med[opt[i]] > med[opt[i + 1]])
268 SWAP(med[opt[i]], med[opt[i + 1]]);
269 pix[0][c] = CLIP(med[4] + pix[0][1]);
270 }
271 }
272 }
273 }
274
blend_highlights()275 void LibRaw::blend_highlights()
276 {
277 int clip = INT_MAX, row, col, c, i, j;
278 static const float trans[2][4][4] = {
279 {{1, 1, 1}, {1.7320508, -1.7320508, 0}, {-1, -1, 2}},
280 {{1, 1, 1, 1}, {1, -1, 1, -1}, {1, 1, -1, -1}, {1, -1, -1, 1}}};
281 static const float itrans[2][4][4] = {
282 {{1, 0.8660254, -0.5}, {1, -0.8660254, -0.5}, {1, 0, 1}},
283 {{1, 1, 1, 1}, {1, -1, 1, -1}, {1, 1, -1, -1}, {1, -1, -1, 1}}};
284 float cam[2][4], lab[2][4], sum[2], chratio;
285
286 if ((unsigned)(colors - 3) > 1)
287 return;
288 RUN_CALLBACK(LIBRAW_PROGRESS_HIGHLIGHTS, 0, 2);
289 FORCC if (clip > (i = 65535 * pre_mul[c])) clip = i;
290 for (row = 0; row < height; row++)
291 for (col = 0; col < width; col++)
292 {
293 FORCC if (image[row * width + col][c] > clip) break;
294 if (c == colors)
295 continue;
296 FORCC
297 {
298 cam[0][c] = image[row * width + col][c];
299 cam[1][c] = MIN(cam[0][c], clip);
300 }
301 for (i = 0; i < 2; i++)
302 {
303 FORCC for (lab[i][c] = j = 0; j < colors; j++) lab[i][c] +=
304 trans[colors - 3][c][j] * cam[i][j];
305 for (sum[i] = 0, c = 1; c < colors; c++)
306 sum[i] += SQR(lab[i][c]);
307 }
308 chratio = sqrt(sum[1] / sum[0]);
309 for (c = 1; c < colors; c++)
310 lab[0][c] *= chratio;
311 FORCC for (cam[0][c] = j = 0; j < colors; j++) cam[0][c] +=
312 itrans[colors - 3][c][j] * lab[0][j];
313 FORCC image[row * width + col][c] = cam[0][c] / colors;
314 }
315 RUN_CALLBACK(LIBRAW_PROGRESS_HIGHLIGHTS, 1, 2);
316 }
317
318 #define SCALE (4 >> shrink)
recover_highlights()319 void LibRaw::recover_highlights()
320 {
321 float *map, sum, wgt, grow;
322 int hsat[4], count, spread, change, val, i;
323 unsigned high, wide, mrow, mcol, row, col, kc, c, d, y, x;
324 ushort *pixel;
325 static const signed char dir[8][2] = {{-1, -1}, {-1, 0}, {-1, 1}, {0, 1},
326 {1, 1}, {1, 0}, {1, -1}, {0, -1}};
327
328 grow = pow(2.0, 4 - highlight);
329 FORC(unsigned(colors)) hsat[c] = 32000 * pre_mul[c];
330 for (kc = 0, c = 1; c < (unsigned)colors; c++)
331 if (pre_mul[kc] < pre_mul[c])
332 kc = c;
333 high = height / SCALE;
334 wide = width / SCALE;
335 map = (float *)calloc(high, wide * sizeof *map);
336 merror(map, "recover_highlights()");
337 FORC(unsigned(colors)) if (c != kc)
338 {
339 RUN_CALLBACK(LIBRAW_PROGRESS_HIGHLIGHTS, c - 1, colors - 1);
340 memset(map, 0, high * wide * sizeof *map);
341 for (mrow = 0; mrow < high; mrow++)
342 for (mcol = 0; mcol < wide; mcol++)
343 {
344 sum = wgt = count = 0;
345 for (row = mrow * SCALE; row < (mrow + 1) * SCALE; row++)
346 for (col = mcol * SCALE; col < (mcol + 1) * SCALE; col++)
347 {
348 pixel = image[row * width + col];
349 if (pixel[c] / hsat[c] == 1 && pixel[kc] > 24000)
350 {
351 sum += pixel[c];
352 wgt += pixel[kc];
353 count++;
354 }
355 }
356 if (count == SCALE * SCALE)
357 map[mrow * wide + mcol] = sum / wgt;
358 }
359 for (spread = 32 / grow; spread--;)
360 {
361 for (mrow = 0; mrow < high; mrow++)
362 for (mcol = 0; mcol < wide; mcol++)
363 {
364 if (map[mrow * wide + mcol])
365 continue;
366 sum = count = 0;
367 for (d = 0; d < 8; d++)
368 {
369 y = mrow + dir[d][0];
370 x = mcol + dir[d][1];
371 if (y < high && x < wide && map[y * wide + x] > 0)
372 {
373 sum += (1 + (d & 1)) * map[y * wide + x];
374 count += 1 + (d & 1);
375 }
376 }
377 if (count > 3)
378 map[mrow * wide + mcol] = -(sum + grow) / (count + grow);
379 }
380 for (change = i = 0; i < int(high * wide); i++)
381 if (map[i] < 0)
382 {
383 map[i] = -map[i];
384 change = 1;
385 }
386 if (!change)
387 break;
388 }
389 for (i = 0; i < int(high * wide); i++)
390 if (map[i] == 0)
391 map[i] = 1;
392 for (mrow = 0; mrow < high; mrow++)
393 for (mcol = 0; mcol < wide; mcol++)
394 {
395 for (row = mrow * SCALE; row < (mrow + 1) * SCALE; row++)
396 for (col = mcol * SCALE; col < (mcol + 1) * SCALE; col++)
397 {
398 pixel = image[row * width + col];
399 if (pixel[c] / hsat[c] > 1)
400 {
401 val = pixel[kc] * map[mrow * wide + mcol];
402 if (pixel[c] < val)
403 pixel[c] = CLIP(val);
404 }
405 }
406 }
407 }
408 free(map);
409 }
410 #undef SCALE
411