1 /* -*- C++ -*-
2 * Copyright 2019-2021 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
pre_interpolate()21 void LibRaw::pre_interpolate()
22 {
23 ushort(*img)[4];
24 int row, col, c;
25 RUN_CALLBACK(LIBRAW_PROGRESS_PRE_INTERPOLATE, 0, 2);
26 if (shrink)
27 {
28 if (half_size)
29 {
30 height = iheight;
31 width = iwidth;
32 if (filters == 9)
33 {
34 for (row = 0; row < 3; row++)
35 for (col = 1; col < 4; col++)
36 if (!(image[row * width + col][0] | image[row * width + col][2]))
37 goto break2;
38 break2:
39 for (; row < height; row += 3)
40 for (col = (col - 1) % 3 + 1; col < width - 1; col += 3)
41 {
42 img = image + row * width + col;
43 for (c = 0; c < 3; c += 2)
44 img[0][c] = (img[-1][c] + img[1][c]) >> 1;
45 }
46 }
47 }
48 else
49 {
50 img = (ushort(*)[4])calloc(height, width * sizeof *img);
51 merror(img, "pre_interpolate()");
52 for (row = 0; row < height; row++)
53 for (col = 0; col < width; col++)
54 {
55 c = fcol(row, col);
56 img[row * width + col][c] =
57 image[(row >> 1) * iwidth + (col >> 1)][c];
58 }
59 free(image);
60 image = img;
61 shrink = 0;
62 }
63 }
64 if (filters > 1000 && colors == 3)
65 {
66 mix_green = four_color_rgb ^ half_size;
67 if (four_color_rgb | half_size)
68 colors++;
69 else
70 {
71 for (row = FC(1, 0) >> 1; row < height; row += 2)
72 for (col = FC(row, 1) & 1; col < width; col += 2)
73 image[row * width + col][1] = image[row * width + col][3];
74 filters &= ~((filters & 0x55555555U) << 1);
75 }
76 }
77 if (half_size)
78 filters = 0;
79 RUN_CALLBACK(LIBRAW_PROGRESS_PRE_INTERPOLATE, 1, 2);
80 }
81
border_interpolate(int border)82 void LibRaw::border_interpolate(int border)
83 {
84 unsigned row, col, y, x, f, c, sum[8];
85
86 for (row = 0; row < height; row++)
87 for (col = 0; col < width; col++)
88 {
89 if (col == (unsigned)border && row >= (unsigned)border && row < (unsigned)(height - border))
90 col = width - border;
91 memset(sum, 0, sizeof sum);
92 for (y = row - 1; y != row + 2; y++)
93 for (x = col - 1; x != col + 2; x++)
94 if (y < height && x < width)
95 {
96 f = fcol(y, x);
97 sum[f] += image[y * width + x][f];
98 sum[f + 4]++;
99 }
100 f = fcol(row, col);
101 FORC(unsigned(colors)) if (c != f && sum[c + 4]) image[row * width + col][c] =
102 sum[c] / sum[c + 4];
103 }
104 }
105
lin_interpolate_loop(int * code,int size)106 void LibRaw::lin_interpolate_loop(int *code, int size)
107 {
108 int row;
109 for (row = 1; row < height - 1; row++)
110 {
111 int col, *ip;
112 ushort *pix;
113 for (col = 1; col < width - 1; col++)
114 {
115 int i;
116 int sum[4];
117 pix = image[row * width + col];
118 ip = code + ((((row % size) * 16) + (col % size)) * 32);
119 memset(sum, 0, sizeof sum);
120 for (i = *ip++; i--; ip += 3)
121 sum[ip[2]] += pix[ip[0]] << ip[1];
122 for (i = colors; --i; ip += 2)
123 pix[ip[0]] = sum[ip[0]] * ip[1] >> 8;
124 }
125 }
126 }
127
lin_interpolate()128 void LibRaw::lin_interpolate()
129 {
130 std::vector<int> code_buffer(16 * 16 * 32);
131 int* code = &code_buffer[0], size = 16, *ip, sum[4];
132 int f, c, x, y, row, col, shift, color;
133
134 RUN_CALLBACK(LIBRAW_PROGRESS_INTERPOLATE, 0, 3);
135
136 if (filters == 9)
137 size = 6;
138 border_interpolate(1);
139 for (row = 0; row < size; row++)
140 for (col = 0; col < size; col++)
141 {
142 ip = code + (((row * 16) + col) * 32) + 1;
143 f = fcol(row, col);
144 memset(sum, 0, sizeof sum);
145 for (y = -1; y <= 1; y++)
146 for (x = -1; x <= 1; x++)
147 {
148 shift = (y == 0) + (x == 0);
149 color = fcol(row + y + 48, col + x + 48);
150 if (color == f)
151 continue;
152 *ip++ = (width * y + x) * 4 + color;
153 *ip++ = shift;
154 *ip++ = color;
155 sum[color] += 1 << shift;
156 }
157 code[(row * 16 + col) * 32] = (ip - (code + ((row * 16) + col) * 32)) / 3;
158 FORCC
159 if (c != f)
160 {
161 *ip++ = c;
162 *ip++ = sum[c] > 0 ? 256 / sum[c] : 0;
163 }
164 }
165 RUN_CALLBACK(LIBRAW_PROGRESS_INTERPOLATE, 1, 3);
166 lin_interpolate_loop(code, size);
167 RUN_CALLBACK(LIBRAW_PROGRESS_INTERPOLATE, 2, 3);
168 }
169
170 /*
171 This algorithm is officially called:
172
173 "Interpolation using a Threshold-based variable number of gradients"
174
175 described in
176 http://scien.stanford.edu/pages/labsite/1999/psych221/projects/99/tingchen/algodep/vargra.html
177
178 I've extended the basic idea to work with non-Bayer filter arrays.
179 Gradients are numbered clockwise from NW=0 to W=7.
180 */
vng_interpolate()181 void LibRaw::vng_interpolate()
182 {
183 static const signed char *cp,
184 terms[] =
185 {-2, -2, +0, -1, 0, 0x01, -2, -2, +0, +0, 1, 0x01, -2, -1, -1,
186 +0, 0, 0x01, -2, -1, +0, -1, 0, 0x02, -2, -1, +0, +0, 0, 0x03,
187 -2, -1, +0, +1, 1, 0x01, -2, +0, +0, -1, 0, 0x06, -2, +0, +0,
188 +0, 1, 0x02, -2, +0, +0, +1, 0, 0x03, -2, +1, -1, +0, 0, 0x04,
189 -2, +1, +0, -1, 1, 0x04, -2, +1, +0, +0, 0, 0x06, -2, +1, +0,
190 +1, 0, 0x02, -2, +2, +0, +0, 1, 0x04, -2, +2, +0, +1, 0, 0x04,
191 -1, -2, -1, +0, 0, -128, -1, -2, +0, -1, 0, 0x01, -1, -2, +1,
192 -1, 0, 0x01, -1, -2, +1, +0, 1, 0x01, -1, -1, -1, +1, 0, -120,
193 -1, -1, +1, -2, 0, 0x40, -1, -1, +1, -1, 0, 0x22, -1, -1, +1,
194 +0, 0, 0x33, -1, -1, +1, +1, 1, 0x11, -1, +0, -1, +2, 0, 0x08,
195 -1, +0, +0, -1, 0, 0x44, -1, +0, +0, +1, 0, 0x11, -1, +0, +1,
196 -2, 1, 0x40, -1, +0, +1, -1, 0, 0x66, -1, +0, +1, +0, 1, 0x22,
197 -1, +0, +1, +1, 0, 0x33, -1, +0, +1, +2, 1, 0x10, -1, +1, +1,
198 -1, 1, 0x44, -1, +1, +1, +0, 0, 0x66, -1, +1, +1, +1, 0, 0x22,
199 -1, +1, +1, +2, 0, 0x10, -1, +2, +0, +1, 0, 0x04, -1, +2, +1,
200 +0, 1, 0x04, -1, +2, +1, +1, 0, 0x04, +0, -2, +0, +0, 1, -128,
201 +0, -1, +0, +1, 1, -120, +0, -1, +1, -2, 0, 0x40, +0, -1, +1,
202 +0, 0, 0x11, +0, -1, +2, -2, 0, 0x40, +0, -1, +2, -1, 0, 0x20,
203 +0, -1, +2, +0, 0, 0x30, +0, -1, +2, +1, 1, 0x10, +0, +0, +0,
204 +2, 1, 0x08, +0, +0, +2, -2, 1, 0x40, +0, +0, +2, -1, 0, 0x60,
205 +0, +0, +2, +0, 1, 0x20, +0, +0, +2, +1, 0, 0x30, +0, +0, +2,
206 +2, 1, 0x10, +0, +1, +1, +0, 0, 0x44, +0, +1, +1, +2, 0, 0x10,
207 +0, +1, +2, -1, 1, 0x40, +0, +1, +2, +0, 0, 0x60, +0, +1, +2,
208 +1, 0, 0x20, +0, +1, +2, +2, 0, 0x10, +1, -2, +1, +0, 0, -128,
209 +1, -1, +1, +1, 0, -120, +1, +0, +1, +2, 0, 0x08, +1, +0, +2,
210 -1, 0, 0x40, +1, +0, +2, +1, 0, 0x10},
211 chood[] = {-1, -1, -1, 0, -1, +1, 0, +1, +1, +1, +1, 0, +1, -1, 0, -1};
212 ushort(*brow[5])[4], *pix;
213 int prow = 8, pcol = 2, *ip, *code[16][16], gval[8], gmin, gmax, sum[4];
214 int row, col, x, y, x1, x2, y1, y2, t, weight, grads, color, diag;
215 int g, diff, thold, num, c;
216
217 lin_interpolate();
218
219 if (filters == 1)
220 prow = pcol = 16;
221 if (filters == 9)
222 prow = pcol = 6;
223 ip = (int *)calloc(prow * pcol, 1280);
224 merror(ip, "vng_interpolate()");
225 for (row = 0; row < prow; row++) /* Precalculate for VNG */
226 for (col = 0; col < pcol; col++)
227 {
228 code[row][col] = ip;
229 for (cp = terms, t = 0; t < 64; t++)
230 {
231 y1 = *cp++;
232 x1 = *cp++;
233 y2 = *cp++;
234 x2 = *cp++;
235 weight = *cp++;
236 grads = *cp++;
237 color = fcol(row + y1 + 144, col + x1 + 144);
238 if (fcol(row + y2 + 144, col + x2 + 144) != color)
239 continue;
240 diag = (fcol(row, col + 1) == color && fcol(row + 1, col) == color) ? 2
241 : 1;
242 if (abs(y1 - y2) == diag && abs(x1 - x2) == diag)
243 continue;
244 *ip++ = (y1 * width + x1) * 4 + color;
245 *ip++ = (y2 * width + x2) * 4 + color;
246 *ip++ = weight;
247 for (g = 0; g < 8; g++)
248 if (grads & 1 << g)
249 *ip++ = g;
250 *ip++ = -1;
251 }
252 *ip++ = INT_MAX;
253 for (cp = chood, g = 0; g < 8; g++)
254 {
255 y = *cp++;
256 x = *cp++;
257 *ip++ = (y * width + x) * 4;
258 color = fcol(row, col);
259 if (fcol(row + y + 144, col + x + 144) != color &&
260 fcol(row + y * 2 + 144, col + x * 2 + 144) == color)
261 *ip++ = (y * width + x) * 8 + color;
262 else
263 *ip++ = 0;
264 }
265 }
266 brow[4] = (ushort(*)[4])calloc(width * 3, sizeof **brow);
267 merror(brow[4], "vng_interpolate()");
268 for (row = 0; row < 3; row++)
269 brow[row] = brow[4] + row * width;
270 for (row = 2; row < height - 2; row++)
271 { /* Do VNG interpolation */
272 if (!((row - 2) % 256))
273 RUN_CALLBACK(LIBRAW_PROGRESS_INTERPOLATE, (row - 2) / 256 + 1,
274 ((height - 3) / 256) + 1);
275 for (col = 2; col < width - 2; col++)
276 {
277 pix = image[row * width + col];
278 ip = code[row % prow][col % pcol];
279 memset(gval, 0, sizeof gval);
280 while ((g = ip[0]) != INT_MAX)
281 { /* Calculate gradients */
282 diff = ABS(pix[g] - pix[ip[1]]) << ip[2];
283 gval[ip[3]] += diff;
284 ip += 5;
285 if ((g = ip[-1]) == -1)
286 continue;
287 gval[g] += diff;
288 while ((g = *ip++) != -1)
289 gval[g] += diff;
290 }
291 ip++;
292 gmin = gmax = gval[0]; /* Choose a threshold */
293 for (g = 1; g < 8; g++)
294 {
295 if (gmin > gval[g])
296 gmin = gval[g];
297 if (gmax < gval[g])
298 gmax = gval[g];
299 }
300 if (gmax == 0)
301 {
302 memcpy(brow[2][col], pix, sizeof *image);
303 continue;
304 }
305 thold = gmin + (gmax >> 1);
306 memset(sum, 0, sizeof sum);
307 color = fcol(row, col);
308 for (num = g = 0; g < 8; g++, ip += 2)
309 { /* Average the neighbors */
310 if (gval[g] <= thold)
311 {
312 FORCC
313 if (c == color && ip[1])
314 sum[c] += (pix[c] + pix[ip[1]]) >> 1;
315 else
316 sum[c] += pix[ip[0] + c];
317 num++;
318 }
319 }
320 FORCC
321 { /* Save to buffer */
322 t = pix[color];
323 if (c != color)
324 t += (sum[c] - sum[color]) / num;
325 brow[2][col][c] = CLIP(t);
326 }
327 }
328 if (row > 3) /* Write buffer to image */
329 memcpy(image[(row - 2) * width + 2], brow[0] + 2,
330 (width - 4) * sizeof *image);
331 for (g = 0; g < 4; g++)
332 brow[(g - 1) & 3] = brow[g];
333 }
334 memcpy(image[(row - 2) * width + 2], brow[0] + 2,
335 (width - 4) * sizeof *image);
336 memcpy(image[(row - 1) * width + 2], brow[1] + 2,
337 (width - 4) * sizeof *image);
338 free(brow[4]);
339 free(code[0][0]);
340 }
341
342 /*
343 Patterned Pixel Grouping Interpolation by Alain Desbiolles
344 */
ppg_interpolate()345 void LibRaw::ppg_interpolate()
346 {
347 int dir[5] = {1, width, -1, -width, 1};
348 int row, col, diff[2], guess[2], c, d, i;
349 ushort(*pix)[4];
350
351 border_interpolate(3);
352
353 /* Fill in the green layer with gradients and pattern recognition: */
354 RUN_CALLBACK(LIBRAW_PROGRESS_INTERPOLATE, 0, 3);
355 #ifdef LIBRAW_USE_OPENMP
356 #pragma omp parallel for default(shared) private(guess, diff, row, col, d, c, \
357 i, pix) schedule(static)
358 #endif
359 for (row = 3; row < height - 3; row++)
360 for (col = 3 + (FC(row, 3) & 1), c = FC(row, col); col < width - 3;
361 col += 2)
362 {
363 pix = image + row * width + col;
364 for (i = 0; i < 2; i++)
365 {
366 d = dir[i];
367 guess[i] = (pix[-d][1] + pix[0][c] + pix[d][1]) * 2 - pix[-2 * d][c] -
368 pix[2 * d][c];
369 diff[i] =
370 (ABS(pix[-2 * d][c] - pix[0][c]) + ABS(pix[2 * d][c] - pix[0][c]) +
371 ABS(pix[-d][1] - pix[d][1])) *
372 3 +
373 (ABS(pix[3 * d][1] - pix[d][1]) +
374 ABS(pix[-3 * d][1] - pix[-d][1])) *
375 2;
376 }
377 d = dir[i = diff[0] > diff[1]];
378 pix[0][1] = ULIM(guess[i] >> 2, pix[d][1], pix[-d][1]);
379 }
380 /* Calculate red and blue for each green pixel: */
381 RUN_CALLBACK(LIBRAW_PROGRESS_INTERPOLATE, 1, 3);
382 #ifdef LIBRAW_USE_OPENMP
383 #pragma omp parallel for default(shared) private(guess, diff, row, col, d, c, \
384 i, pix) schedule(static)
385 #endif
386 for (row = 1; row < height - 1; row++)
387 for (col = 1 + (FC(row, 2) & 1), c = FC(row, col + 1); col < width - 1;
388 col += 2)
389 {
390 pix = image + row * width + col;
391 for (i = 0; i < 2; c = 2 - c, i++)
392 {
393 d = dir[i];
394 pix[0][c] = CLIP(
395 (pix[-d][c] + pix[d][c] + 2 * pix[0][1] - pix[-d][1] - pix[d][1]) >>
396 1);
397 }
398 }
399 /* Calculate blue for red pixels and vice versa: */
400 RUN_CALLBACK(LIBRAW_PROGRESS_INTERPOLATE, 2, 3);
401 #ifdef LIBRAW_USE_OPENMP
402 #pragma omp parallel for default(shared) private(guess, diff, row, col, d, c, \
403 i, pix) schedule(static)
404 #endif
405 for (row = 1; row < height - 1; row++)
406 for (col = 1 + (FC(row, 1) & 1), c = 2 - FC(row, col); col < width - 1;
407 col += 2)
408 {
409 pix = image + row * width + col;
410 for (i = 0; i < 2; i++)
411 {
412 d = dir[i] + dir[i+1];
413 diff[i] = ABS(pix[-d][c] - pix[d][c]) + ABS(pix[-d][1] - pix[0][1]) +
414 ABS(pix[d][1] - pix[0][1]);
415 guess[i] =
416 pix[-d][c] + pix[d][c] + 2 * pix[0][1] - pix[-d][1] - pix[d][1];
417 }
418 if (diff[0] != diff[1])
419 pix[0][c] = CLIP(guess[diff[0] > diff[1]] >> 1);
420 else
421 pix[0][c] = CLIP((guess[0] + guess[1]) >> 2);
422 }
423 }
424