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
21 #define fcol(row, col) xtrans[(row + 6) % 6][(col + 6) % 6]
22 /*
23 Frank Markesteijn's algorithm for Fuji X-Trans sensors
24 */
xtrans_interpolate(int passes)25 void LibRaw::xtrans_interpolate(int passes)
26 {
27 int cstat[4] = {0, 0, 0, 0};
28 int ndir;
29 static const short orth[12] = {1, 0, 0, 1, -1, 0, 0, -1, 1, 0, 0, 1},
30 patt[2][16] = {{0, 1, 0, -1, 2, 0, -1, 0, 1, 1, 1, -1, 0,
31 0, 0, 0},
32 {0, 1, 0, -2, 1, 0, -2, 0, 1, 1, -2, -2, 1,
33 -1, -1, 1}},
34 dir[4] = {1, LIBRAW_AHD_TILE, LIBRAW_AHD_TILE + 1,
35 LIBRAW_AHD_TILE - 1};
36 short allhex[3][3][2][8];
37 ushort sgrow = 0, sgcol = 0;
38
39 if (width < LIBRAW_AHD_TILE || height < LIBRAW_AHD_TILE)
40 throw LIBRAW_EXCEPTION_IO_CORRUPT; // too small image
41 /* Check against right pattern */
42 for (int row = 0; row < 6; row++)
43 for (int col = 0; col < 6; col++)
44 cstat[(unsigned)fcol(row, col)]++;
45
46 if (cstat[0] < 6 || cstat[0] > 10 || cstat[1] < 16 || cstat[1] > 24 ||
47 cstat[2] < 6 || cstat[2] > 10 || cstat[3])
48 throw LIBRAW_EXCEPTION_IO_CORRUPT;
49
50 // Init allhex table to unreasonable values
51 for (int i = 0; i < 3; i++)
52 for (int j = 0; j < 3; j++)
53 for (int k = 0; k < 2; k++)
54 for (int l = 0; l < 8; l++)
55 allhex[i][j][k][l] = 32700;
56
57 cielab(0, 0);
58 ndir = 4 << (passes > 1);
59
60 int minv = 0, maxv = 0, minh = 0, maxh = 0;
61 /* Map a green hexagon around each non-green pixel and vice versa: */
62 for (int row = 0; row < 3; row++)
63 for (int col = 0; col < 3; col++)
64 for (int ng = 0, d = 0; d < 10; d += 2)
65 {
66 int g = fcol(row, col) == 1;
67 if (fcol(row + orth[d], col + orth[d + 2]) == 1)
68 ng = 0;
69 else
70 ng++;
71 if (ng == 4)
72 {
73 sgrow = row;
74 sgcol = col;
75 }
76 if (ng == g + 1)
77 {
78 int c;
79 FORC(8)
80 {
81 int v = orth[d] * patt[g][c * 2] + orth[d + 1] * patt[g][c * 2 + 1];
82 int h = orth[d + 2] * patt[g][c * 2] + orth[d + 3] * patt[g][c * 2 + 1];
83 minv = MIN(v, minv);
84 maxv = MAX(v, maxv);
85 minh = MIN(v, minh);
86 maxh = MAX(v, maxh);
87 allhex[row][col][0][c ^ (g * 2 & d)] = h + v * width;
88 allhex[row][col][1][c ^ (g * 2 & d)] = h + v * LIBRAW_AHD_TILE;
89 }
90 }
91 }
92
93 // Check allhex table initialization
94 for (int i = 0; i < 3; i++)
95 for (int j = 0; j < 3; j++)
96 for (int k = 0; k < 2; k++)
97 for (int l = 0; l < 8; l++)
98 if (allhex[i][j][k][l] > maxh + maxv * width + 1 ||
99 allhex[i][j][k][l] < minh + minv * width - 1)
100 throw LIBRAW_EXCEPTION_IO_CORRUPT;
101 int retrycount = 0;
102
103 /* Set green1 and green3 to the minimum and maximum allowed values: */
104 for (int row = 2; row < height - 2; row++)
105 {
106 int col;
107 ushort min, max;
108 for (col = 2, max = 0, min = ~0; col < width - 2; col++)
109 {
110 if (fcol(row, col) == 1 && (min = ~(max = 0)))
111 continue;
112 ushort(*pix)[4];
113 pix = image + row * width + col;
114 short* hex = allhex[row % 3][col % 3][0];
115 if (!max)
116 {
117 int c;
118 FORC(6)
119 {
120 int val = pix[hex[c]][1];
121 if (min > val)
122 min = val;
123 if (max < val)
124 max = val;
125 }
126 }
127 pix[0][1] = min;
128 pix[0][3] = max;
129 switch ((row - sgrow) % 3)
130 {
131 case 1:
132 if (row < height - 3)
133 {
134 row++;
135 col--;
136 }
137 break;
138 case 2:
139 if ((min = ~(max = 0)) && (col += 2) < width - 3 && row > 2)
140 {
141 row--;
142 if (retrycount++ > width * height)
143 throw LIBRAW_EXCEPTION_IO_CORRUPT;
144 }
145 }
146 }
147 }
148
149 for (int row = 3; row < 9 && row < height - 3; row++)
150 for (int col = 3; col < 9 && col < width - 3; col++)
151 {
152 if ((fcol(row, col)) == 1)
153 continue;
154 short* hex = allhex[row % 3][col % 3][0];
155 int c;
156 FORC(2)
157 {
158 int idx3 = 3 * hex[4 + c] + row * width + col;
159 int idx4 = -3 * hex[4 + c] + row * width + col;
160 int maxidx = width * height;
161 if (idx3 < 0 || idx3 >= maxidx)
162 throw LIBRAW_EXCEPTION_IO_CORRUPT;
163 if (idx4 < 0 || idx4 >= maxidx)
164 throw LIBRAW_EXCEPTION_IO_CORRUPT;
165 }
166 }
167
168 #if defined(LIBRAW_USE_OPENMP)
169 int buffer_count = omp_get_max_threads();
170 #else
171 int buffer_count = 1;
172 #endif
173
174 size_t buffer_size = LIBRAW_AHD_TILE * LIBRAW_AHD_TILE * (ndir * 11 + 6);
175 char** buffers = malloc_omp_buffers(buffer_count, buffer_size, "xtrans_interpolate()");
176
177 #if defined(LIBRAW_USE_OPENMP)
178 # pragma omp parallel for schedule(dynamic) default(none) firstprivate(buffers, allhex, passes, sgrow, sgcol, ndir) shared(dir)
179 #endif
180 for (int top = 3; top < height - 19; top += LIBRAW_AHD_TILE - 16)
181 {
182 #if defined(LIBRAW_USE_OPENMP)
183 char* buffer = buffers[omp_get_thread_num()];
184 #else
185 char* buffer = buffers[0];
186 #endif
187
188 ushort(*rgb)[LIBRAW_AHD_TILE][LIBRAW_AHD_TILE][3], (*rix)[3];
189 short(*lab)[LIBRAW_AHD_TILE][3], (*lix)[3];
190 float(*drv)[LIBRAW_AHD_TILE][LIBRAW_AHD_TILE];
191 char(*homo)[LIBRAW_AHD_TILE][LIBRAW_AHD_TILE];
192
193 rgb = (ushort(*)[LIBRAW_AHD_TILE][LIBRAW_AHD_TILE][3])buffer;
194 lab = (short(*)[LIBRAW_AHD_TILE][3])(
195 buffer + LIBRAW_AHD_TILE * LIBRAW_AHD_TILE * (ndir * 6));
196 drv = (float(*)[LIBRAW_AHD_TILE][LIBRAW_AHD_TILE])(
197 buffer + LIBRAW_AHD_TILE * LIBRAW_AHD_TILE * (ndir * 6 + 6));
198 homo = (char(*)[LIBRAW_AHD_TILE][LIBRAW_AHD_TILE])(
199 buffer + LIBRAW_AHD_TILE * LIBRAW_AHD_TILE * (ndir * 10 + 6));
200
201 for (int left = 3; left < width - 19; left += LIBRAW_AHD_TILE - 16)
202 {
203 int mrow = MIN(top + LIBRAW_AHD_TILE, height - 3);
204 int mcol = MIN(left + LIBRAW_AHD_TILE, width - 3);
205 for (int row = top; row < mrow; row++)
206 for (int col = left; col < mcol; col++)
207 memcpy(rgb[0][row - top][col - left], image[row * width + col], 6);
208 int c;
209 FORC3 memcpy(rgb[c + 1], rgb[0], sizeof * rgb);
210
211 /* Interpolate green horizontally, vertically, and along both diagonals:
212 */
213 int color[3][8];
214 for (int row = top; row < mrow; row++)
215 for (int col = left; col < mcol; col++)
216 {
217 int f;
218 if ((f = fcol(row, col)) == 1)
219 continue;
220 ushort (*pix)[4] = image + row * width + col;
221 short* hex = allhex[row % 3][col % 3][0];
222 color[1][0] = 174 * (pix[hex[1]][1] + pix[hex[0]][1]) -
223 46 * (pix[2 * hex[1]][1] + pix[2 * hex[0]][1]);
224 color[1][1] = 223 * pix[hex[3]][1] + pix[hex[2]][1] * 33 +
225 92 * (pix[0][f] - pix[-hex[2]][f]);
226 FORC(2)
227 color[1][2 + c] = 164 * pix[hex[4 + c]][1] +
228 92 * pix[-2 * hex[4 + c]][1] +
229 33 * (2 * pix[0][f] - pix[3 * hex[4 + c]][f] -
230 pix[-3 * hex[4 + c]][f]);
231 FORC4 rgb[c ^ !((row - sgrow) % 3)][row - top][col - left][1] =
232 LIM(color[1][c] >> 8, pix[0][1], pix[0][3]);
233 }
234
235 for (int pass = 0; pass < passes; pass++)
236 {
237 if (pass == 1)
238 memcpy(rgb += 4, buffer, 4 * sizeof * rgb);
239
240 /* Recalculate green from interpolated values of closer pixels: */
241 if (pass)
242 {
243 for (int row = top + 2; row < mrow - 2; row++)
244 for (int col = left + 2; col < mcol - 2; col++)
245 {
246 int f;
247 if ((f = fcol(row, col)) == 1)
248 continue;
249 ushort(*pix)[4] = image + row * width + col;
250 short* hex = allhex[row % 3][col % 3][1];
251 for (int d = 3; d < 6; d++)
252 {
253 rix =
254 &rgb[(d - 2) ^ !((row - sgrow) % 3)][row - top][col - left];
255 int val = rix[-2 * hex[d]][1] + 2 * rix[hex[d]][1] -
256 rix[-2 * hex[d]][f] - 2 * rix[hex[d]][f] + 3 * rix[0][f];
257 rix[0][1] = LIM(val / 3, pix[0][1], pix[0][3]);
258 }
259 }
260 }
261
262 /* Interpolate red and blue values for solitary green pixels: */
263 for (int row = (top - sgrow + 4) / 3 * 3 + sgrow; row < mrow - 2; row += 3)
264 for (int col = (left - sgcol + 4) / 3 * 3 + sgcol; col < mcol - 2; col += 3)
265 {
266 rix = &rgb[0][row - top][col - left];
267 int h = fcol(row, col + 1);
268 float diff[6];
269 memset(diff, 0, sizeof diff);
270 for (int i = 1, d = 0; d < 6; d++, i ^= LIBRAW_AHD_TILE ^ 1, h ^= 2)
271 {
272 for (c = 0; c < 2; c++, h ^= 2)
273 {
274 int g = 2 * rix[0][1] - rix[i << c][1] - rix[-i << c][1];
275 color[h][d] = g + rix[i << c][h] + rix[-i << c][h];
276 if (d > 1)
277 diff[d] += SQR((float)rix[i << c][1] - (float)rix[-i << c][1] -
278 (float)rix[i << c][h] + (float)rix[-i << c][h]) + SQR((float)g);
279 }
280 if (d > 1 && (d & 1))
281 if (diff[d - 1] < diff[d])
282 FORC(2) color[c * 2][d] = color[c * 2][d - 1];
283 if (d < 2 || (d & 1))
284 {
285 FORC(2) rix[0][c * 2] = CLIP(color[c * 2][d] / 2);
286 rix += LIBRAW_AHD_TILE * LIBRAW_AHD_TILE;
287 }
288 }
289 }
290
291 /* Interpolate red for blue pixels and vice versa: */
292 for (int row = top + 3; row < mrow - 3; row++)
293 for (int col = left + 3; col < mcol - 3; col++)
294 {
295 int f;
296 if ((f = 2 - fcol(row, col)) == 1)
297 continue;
298 rix = &rgb[0][row - top][col - left];
299 c = (row - sgrow) % 3 ? LIBRAW_AHD_TILE : 1;
300 int h = 3 * (c ^ LIBRAW_AHD_TILE ^ 1);
301 for (int d = 0; d < 4; d++, rix += LIBRAW_AHD_TILE * LIBRAW_AHD_TILE)
302 {
303 int i = d > 1 || ((d ^ c) & 1) ||
304 ((ABS(rix[0][1] - rix[c][1]) +
305 ABS(rix[0][1] - rix[-c][1])) <
306 2 * (ABS(rix[0][1] - rix[h][1]) +
307 ABS(rix[0][1] - rix[-h][1])))
308 ? c
309 : h;
310 rix[0][f] = CLIP((rix[i][f] + rix[-i][f] + 2 * rix[0][1] -
311 rix[i][1] - rix[-i][1]) /
312 2);
313 }
314 }
315
316 /* Fill in red and blue for 2x2 blocks of green: */
317 for (int row = top + 2; row < mrow - 2; row++)
318 if ((row - sgrow) % 3)
319 for (int col = left + 2; col < mcol - 2; col++)
320 if ((col - sgcol) % 3)
321 {
322 rix = &rgb[0][row - top][col - left];
323 short* hex = allhex[row % 3][col % 3][1];
324 for (int d = 0; d < ndir;
325 d += 2, rix += LIBRAW_AHD_TILE * LIBRAW_AHD_TILE)
326 if (hex[d] + hex[d + 1])
327 {
328 int g = 3 * rix[0][1] - 2 * rix[hex[d]][1] - rix[hex[d + 1]][1];
329 for (c = 0; c < 4; c += 2)
330 rix[0][c] = CLIP(
331 (g + 2 * rix[hex[d]][c] + rix[hex[d + 1]][c]) / 3);
332 }
333 else
334 {
335 int g = 2 * rix[0][1] - rix[hex[d]][1] - rix[hex[d + 1]][1];
336 for (c = 0; c < 4; c += 2)
337 rix[0][c] =
338 CLIP((g + rix[hex[d]][c] + rix[hex[d + 1]][c]) / 2);
339 }
340 }
341 }
342 rgb = (ushort(*)[LIBRAW_AHD_TILE][LIBRAW_AHD_TILE][3])buffer;
343 mrow -= top;
344 mcol -= left;
345
346 /* Convert to CIELab and differentiate in all directions: */
347 // no effect
348 for (int d = 0; d < ndir; d++)
349 {
350 for (int row = 2; row < mrow - 2; row++)
351 for (int col = 2; col < mcol - 2; col++)
352 cielab(rgb[d][row][col], lab[row][col]);
353 for (int f = dir[d & 3], row = 3; row < mrow - 3; row++)
354 for (int col = 3; col < mcol - 3; col++)
355 {
356 lix = &lab[row][col];
357 int g = 2 * lix[0][0] - lix[f][0] - lix[-f][0];
358 drv[d][row][col] =
359 SQR(g) +
360 SQR((2 * lix[0][1] - lix[f][1] - lix[-f][1] + g * 500 / 232)) +
361 SQR((2 * lix[0][2] - lix[f][2] - lix[-f][2] - g * 500 / 580));
362 }
363 }
364
365 /* Build homogeneity maps from the derivatives: */
366 memset(homo, 0, ndir * LIBRAW_AHD_TILE * LIBRAW_AHD_TILE);
367 for (int row = 4; row < mrow - 4; row++)
368 for (int col = 4; col < mcol - 4; col++)
369 {
370 int d;
371 float tr;
372 for (tr = FLT_MAX, d = 0; d < ndir; d++)
373 if (tr > drv[d][row][col])
374 tr = drv[d][row][col];
375 tr *= 8;
376 for (int d = 0; d < ndir; d++)
377 for (int v = -1; v <= 1; v++)
378 for (int h = -1; h <= 1; h++)
379 if (drv[d][row + v][col + h] <= tr)
380 homo[d][row][col]++;
381 }
382
383 /* Average the most homogenous pixels for the final result: */
384 if (height - top < LIBRAW_AHD_TILE + 4)
385 mrow = height - top + 2;
386 if (width - left < LIBRAW_AHD_TILE + 4)
387 mcol = width - left + 2;
388 for (int row = MIN(top, 8); row < mrow - 8; row++)
389 for (int col = MIN(left, 8); col < mcol - 8; col++)
390 {
391 int v;
392 int hm[8];
393 for (int d = 0; d < ndir; d++)
394 for (v = -2, hm[d] = 0; v <= 2; v++)
395 for (int h = -2; h <= 2; h++)
396 hm[d] += homo[d][row + v][col + h];
397 for (int d = 0; d < ndir - 4; d++)
398 if (hm[d] < hm[d + 4])
399 hm[d] = 0;
400 else if (hm[d] > hm[d + 4])
401 hm[d + 4] = 0;
402 ushort max;
403 int d;
404 for (d = 1, max = hm[0]; d < ndir; d++)
405 if (max < hm[d])
406 max = hm[d];
407 max -= max >> 3;
408
409 int avg[4];
410 memset(avg, 0, sizeof avg);
411 for (int d = 0; d < ndir; d++)
412 if (hm[d] >= max)
413 {
414 FORC3 avg[c] += rgb[d][row][col][c];
415 avg[3]++;
416 }
417 FORC3 image[(row + top) * width + col + left][c] = avg[c] / avg[3];
418 }
419 }
420 }
421
422 #ifdef LIBRAW_USE_OPENMP
423 #pragma omp barrier
424 #endif
425
426 free_omp_buffers(buffers, buffer_count);
427
428 border_interpolate(8);
429 }
430 #undef fcol
431