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