Ansel 0.0
A darktable fork - bloat + design vision
Loading...
Searching...
No Matches
markesteijn.c
Go to the documentation of this file.
1/*
2 This file is part of the Ansel project.
3 Copyright (C) 2023, 2025-2026 Aurélien PIERRE.
4
5 Ansel is free software: you can redistribute it and/or modify
6 it under the terms of the GNU General Public License as published by
7 the Free Software Foundation, either version 3 of the License, or
8 (at your option) any later version.
9
10 Ansel is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU General Public License for more details.
14
15 You should have received a copy of the GNU General Public License
16 along with Ansel. If not, see <http://www.gnu.org/licenses/>.
17*/
18
19//
20// x-trans specific demosaicing algorithms
21//
22
23// xtrans_interpolate adapted from dcraw 9.20
24
25#define SQR(x) ((x) * (x))
26// tile size, optimized to keep data in L2 cache
27#define TS 122
28
30static inline const short * hexmap(const int row, const int col, short (*const allhex)[3][8])
31{
32 // Row and column offsets may be negative, but C's modulo function
33 // is not useful here with a negative dividend. To be safe, add a
34 // fairly large multiple of 3. In current code row and col will
35 // never be less than -9 (1-pass) or -14 (3-pass).
36 int irow = row + 600;
37 int icol = col + 600;
38 assert(irow >= 0 && icol >= 0);
39 return allhex[irow % 3][icol % 3];
40}
41
42
43/*
44 Frank Markesteijn's algorithm for Fuji X-Trans sensors
45*/
46__DT_CLONE_TARGETS__
47static void xtrans_markesteijn_interpolate(float *out, const float *const in,
48 const dt_iop_roi_t *const roi_out,
49 const dt_iop_roi_t *const roi_in,
50 const uint8_t (*const xtrans)[6], const int passes)
51{
52 static const short orth[12] = { 1, 0, 0, 1, -1, 0, 0, -1, 1, 0, 0, 1 },
53 patt[2][16] = { { 0, 1, 0, -1, 2, 0, -1, 0, 1, 1, 1, -1, 0, 0, 0, 0 },
54 { 0, 1, 0, -2, 1, 0, -2, 0, 1, 1, -2, -2, 1, -1, -1, 1 } },
55 dir[4] = { 1, TS, TS + 1, TS - 1 };
56
57 short allhex[3][3][8];
58 // sgrow/sgcol is the offset in the sensor matrix of the solitary
59 // green pixels (initialized here only to avoid compiler warning)
60 unsigned short sgrow = 0, sgcol = 0;
61
62 const int width = roi_out->width;
63 const int height = roi_out->height;
64 const int ndir = (passes > 1) ? 8 : 4; // only ever 4 or 8
65
66 const size_t buffer_size = (size_t)TS * TS * (ndir * 4 + 3) * sizeof(float);
67 size_t padded_buffer_size;
68 char *const all_buffers = (char *)dt_pixelpipe_cache_alloc_perthread(buffer_size, sizeof(char), &padded_buffer_size);
69 if(IS_NULL_PTR(all_buffers))
70 {
71 printf("[demosaic] not able to allocate Markesteijn buffers\n");
72 return;
73 }
74
75 /* Map a green hexagon around each non-green pixel and vice versa. allhex[] is later
76 * indexed by hexmap() using tile-local (roi_in-relative) row/col taken mod 3, so it
77 * must be built against that same absolute phase: pass roi_in here (not NULL) so the
78 * color lookups land on the actual sensor position, not xtrans[][]'s own local origin. */
79 for(int row = 0; row < 3; row++)
80 for(int col = 0; col < 3; col++)
81 for(int ng = 0, d = 0; d < 10; d += 2)
82 {
83 const int g = FCxtrans(row, col, roi_in, xtrans) == 1;
84 if(FCxtrans(row + orth[d], col + orth[d + 2], roi_in, xtrans) == 1)
85 ng = 0;
86 else
87 ng++;
88 // if there are four non-green pixels adjacent in cardinal
89 // directions, this is the solitary green pixel
90 if(ng == 4)
91 {
92 sgrow = row;
93 sgcol = col;
94 }
95 if(ng == g + 1)
96 for(int c = 0; c < 8; c++)
97 {
98 const int v = orth[d] * patt[g][c * 2] + orth[d + 1] * patt[g][c * 2 + 1];
99 const int h = orth[d + 2] * patt[g][c * 2] + orth[d + 3] * patt[g][c * 2 + 1];
100 // offset within TSxTS buffer
101 allhex[row][col][c ^ (g * 2 & d)] = h + v * TS;
102 }
103 }
104
105 // extra passes propagates out errors at edges, hence need more padding
106 const int pad_tile = (passes == 1) ? 12 : 17;
107 __OMP_PARALLEL_FOR__()
108 // step through TSxTS cells of image, each tile overlapping the
109 // prior as interpolation needs a substantial border
110 for(int top = -pad_tile; top < height - pad_tile; top += TS - (pad_tile*2))
111 {
112 char *const buffer = dt_get_perthread(all_buffers, padded_buffer_size);
113 // rgb points to ndir TSxTS tiles of 3 channels (R, G, and B)
114 float(*rgb)[TS][TS][3] = (float(*)[TS][TS][3])buffer;
115 // yuv points to 3 channel (Y, u, and v) TSxTS tiles
116 // note that channels come before tiles to allow for a
117 // vectorization optimization when building drv[] from yuv[]
118 float (*const yuv)[TS][TS] = (float(*)[TS][TS])(buffer + TS * TS * (ndir * 3) * sizeof(float));
119 // drv points to ndir TSxTS tiles, each a single channel of derivatives
120 float (*const drv)[TS][TS] = (float(*)[TS][TS])(buffer + TS * TS * (ndir * 3 + 3) * sizeof(float));
121 // gmin and gmax reuse memory which is used later by yuv buffer;
122 // each points to a TSxTS tile of single channel data
123 float (*const gmin)[TS] = (float(*)[TS])(buffer + TS * TS * (ndir * 3) * sizeof(float));
124 float (*const gmax)[TS] = (float(*)[TS])(buffer + TS * TS * (ndir * 3 + 1) * sizeof(float));
125 // homo and homosum reuse memory which is used earlier in the
126 // loop; each points to ndir single-channel TSxTS tiles
127 uint8_t (*const homo)[TS][TS] = (uint8_t(*)[TS][TS])(buffer + TS * TS * (ndir * 3) * sizeof(float));
128 uint8_t (*const homosum)[TS][TS] = (uint8_t(*)[TS][TS])(buffer + TS * TS * (ndir * 3) * sizeof(float)
129 + TS * TS * ndir * sizeof(uint8_t));
130
131 for(int left = -pad_tile; left < width - pad_tile; left += TS - (pad_tile*2))
132 {
133 int mrow = MIN(top + TS, height + pad_tile);
134 int mcol = MIN(left + TS, width + pad_tile);
135
136 // Copy current tile from in to image buffer. If border goes
137 // beyond edges of image, fill with mirrored/interpolated edges.
138 // The extra border avoids discontinuities at image edges.
139 for(int row = top; row < mrow; row++)
140 for(int col = left; col < mcol; col++)
141 {
142 float(*const pix) = rgb[0][row - top][col - left];
143 if((col >= 0) && (row >= 0) && (col < width) && (row < height))
144 {
145 const int f = FCxtrans(row, col, roi_in, xtrans);
146 for(int c = 0; c < 3; c++) pix[c] = (c == f) ? in[roi_in->width * row + col] : 0.f;
147 }
148 else
149 {
150 // mirror a border pixel if beyond image edge
151 const int c = FCxtrans(row, col, roi_in, xtrans);
152 for(int cc = 0; cc < 3; cc++)
153 {
154 if(cc != c)
155 pix[cc] = 0.0f;
156 else
157 {
158#define TRANSLATE(n, size) ((n >= size) ? (2 * size - n - 2) : abs(n))
159 const int cy = TRANSLATE(row, height), cx = TRANSLATE(col, width);
160 if(c == FCxtrans(cy, cx, roi_in, xtrans))
161 pix[c] = in[roi_in->width * cy + cx];
162 else
163 {
164 // interpolate if mirror pixel is a different color
165 float sum = 0.0f;
166 uint8_t count = 0;
167 for(int y = row - 1; y <= row + 1; y++)
168 for(int x = col - 1; x <= col + 1; x++)
169 {
170 const int yy = TRANSLATE(y, height), xx = TRANSLATE(x, width);
171 const int ff = FCxtrans(yy, xx, roi_in, xtrans);
172 if(ff == c)
173 {
174 sum += in[roi_in->width * yy + xx];
175 count++;
176 }
177 }
178 pix[c] = sum / count;
179 }
180 }
181 }
182 }
183 }
184
185 // duplicate rgb[0] to rgb[1], rgb[2], and rgb[3]
186 for(int c = 1; c <= 3; c++) memcpy(rgb[c], rgb[0], sizeof(*rgb));
187
188 // note that successive calculations are inset within the tile
189 // so as to give enough border data, and there needs to be a 6
190 // pixel border initially to allow allhex to find neighboring
191 // pixels
192
193 /* Set green1 and green3 to the minimum and maximum allowed values: */
194 // Run through each red/blue or blue/red pair, setting their g1
195 // and g3 values to the min/max of green pixels surrounding the
196 // pair. Use a 3 pixel border as gmin/gmax is used by
197 // interpolate green which has a 3 pixel border.
198 const int pad_g1_g3 = 3;
199 for(int row = top + pad_g1_g3; row < mrow - pad_g1_g3; row++)
200 {
201 // setting max to 0.0f signifies that this is a new pair, which
202 // requires a new min/max calculation of its neighboring greens
203 float min = FLT_MAX, max = 0.0f;
204 for(int col = left + pad_g1_g3; col < mcol - pad_g1_g3; col++)
205 {
206 // if in row of horizontal red & blue pairs (or processing
207 // vertical red & blue pairs near image bottom), reset min/max
208 // between each pair
209 if(FCxtrans(row, col, roi_in, xtrans) == 1)
210 {
211 min = FLT_MAX, max = 0.0f;
212 continue;
213 }
214 // if at start of red & blue pair, calculate min/max of green
215 // pixels surrounding it; note that while normally using == to
216 // compare floats is suspect, here the check is if 0.0f has
217 // explicitly been assigned to max (which signifies a new
218 // red/blue pair)
219 if(max == 0.0f)
220 {
221 float (*const pix)[3] = &rgb[0][row - top][col - left];
222 const short *const hex = hexmap(row,col,allhex);
223 for(int c = 0; c < 6; c++)
224 {
225 const float val = pix[hex[c]][1];
226 if(min > val) min = val;
227 if(max < val) max = val;
228 }
229 }
230 gmin[row - top][col - left] = min;
231 gmax[row - top][col - left] = max;
232 // handle vertical red/blue pairs
233 switch((row - sgrow) % 3)
234 {
235 // hop down a row to second pixel in vertical pair
236 case 1:
237 if(row < mrow - 4) row++, col--;
238 break;
239 // then if not done with the row hop up and right to next
240 // vertical red/blue pair, resetting min/max
241 case 2:
242 min = FLT_MAX, max = 0.0f;
243 if((col += 2) < mcol - 4 && row > top + 3) row--;
244 }
245 }
246 }
247
248 /* Interpolate green horizontally, vertically, and along both diagonals: */
249 // need a 3 pixel border here as 3*hex[] can have a 3 unit offset
250 const int pad_g_interp = 3;
251 for(int row = top + pad_g_interp; row < mrow - pad_g_interp; row++)
252 for(int col = left + pad_g_interp; col < mcol - pad_g_interp; col++)
253 {
254 float color[8];
255 const int f = FCxtrans(row, col, roi_in, xtrans);
256 if(f == 1) continue;
257 float (*const pix)[3] = &rgb[0][row - top][col - left];
258 const short *const hex = hexmap(row,col,allhex);
259 // TODO: these constants come from integer math constants in
260 // dcraw -- calculate them instead from interpolation math
261 color[0] = 0.6796875f * (pix[hex[1]][1] + pix[hex[0]][1])
262 - 0.1796875f * (pix[2 * hex[1]][1] + pix[2 * hex[0]][1]);
263 color[1] = 0.87109375f * pix[hex[3]][1] + pix[hex[2]][1] * 0.13f
264 + 0.359375f * (pix[0][f] - pix[-hex[2]][f]);
265 for(int c = 0; c < 2; c++)
266 color[2 + c] = 0.640625f * pix[hex[4 + c]][1] + 0.359375f * pix[-2 * hex[4 + c]][1]
267 + 0.12890625f * (2 * pix[0][f] - pix[3 * hex[4 + c]][f] - pix[-3 * hex[4 + c]][f]);
268 for(int c = 0; c < 4; c++)
269 rgb[c ^ !((row - sgrow) % 3)][row - top][col - left][1]
270 = CLAMPS(color[c], gmin[row - top][col - left], gmax[row - top][col - left]);
271 }
272
273 for(int pass = 0; pass < passes; pass++)
274 {
275 if(pass == 1)
276 {
277 // if on second pass, copy rgb[0] to [3] into rgb[4] to [7],
278 // and process that second set of buffers
279 memcpy(rgb + 4, rgb, sizeof(*rgb) * 4);
280 rgb += 4;
281 }
282
283 /* Recalculate green from interpolated values of closer pixels: */
284 if(pass)
285 {
286 const int pad_g_recalc = 6;
287 for(int row = top + pad_g_recalc; row < mrow - pad_g_recalc; row++)
288 for(int col = left + pad_g_recalc; col < mcol - pad_g_recalc; col++)
289 {
290 const int f = FCxtrans(row, col, roi_in, xtrans);
291 if(f == 1) continue;
292 const short *const hex = hexmap(row,col,allhex);
293 for(int d = 3; d < 6; d++)
294 {
295 float(*rfx)[3] = &rgb[(d - 2) ^ !((row - sgrow) % 3)][row - top][col - left];
296 const float val = rfx[-2 * hex[d]][1]
297 + 2 * rfx[hex[d]][1] - rfx[-2 * hex[d]][f]
298 - 2 * rfx[hex[d]][f] + 3 * rfx[0][f];
299 rfx[0][1] = CLAMPS(val / 3.0f, gmin[row - top][col - left], gmax[row - top][col - left]);
300 }
301 }
302 }
303
304 /* Interpolate red and blue values for solitary green pixels: */
305 const int pad_rb_g = (passes == 1) ? 6 : 5;
306 for(int row = (top - sgrow + pad_rb_g + 2) / 3 * 3 + sgrow; row < mrow - pad_rb_g; row += 3)
307 for(int col = (left - sgcol + pad_rb_g + 2) / 3 * 3 + sgcol; col < mcol - pad_rb_g; col += 3)
308 {
309 float(*rfx)[3] = &rgb[0][row - top][col - left];
310 int h = FCxtrans(row, col + 1, roi_in, xtrans);
311 float diff[6] = { 0.0f };
312 // interplated color: first index is red/blue, second is
313 // pass, is double actual result
314 float color[2][6];
315 // Six passes, alternating hori/vert interp (i),
316 // starting with R or B (h) depending on which is closest.
317 // Passes 0,1 to rgb[0], rgb[1] of hori/vert interp. Pass
318 // 3,5 to rgb[2], rgb[3] of best of interp hori/vert
319 // results. Each pass which outputs moves on to the next
320 // rgb[] for input of interp greens.
321 for(int i = 1, d = 0; d < 6; d++, i ^= TS ^ 1, h ^= 2)
322 {
323 // look 1 and 2 pixels distance from solitary green to
324 // red then blue or blue then red
325 for(int c = 0; c < 2; c++, h ^= 2)
326 {
327 // rate of change in greens between current pixel and
328 // interpolated pixels 1 or 2 distant: a quick
329 // derivative which will be divided by two later to be
330 // rate of luminance change for red/blue between known
331 // red/blue neighbors and the current unknown pixel
332 const float g = 2 * rfx[0][1] - rfx[i << c][1] - rfx[-(i << c)][1];
333 // color is halved before being stored in rgb, hence
334 // this becomes green rate of change plus the average
335 // of the near red or blue pixels on current axis
336 color[h != 0][d] = g + rfx[i << c][h] + rfx[-(i << c)][h];
337 // Note that diff will become the slope for both red
338 // and blue differentials in the current direction.
339 // For 2nd and 3rd hori+vert passes, create a sum of
340 // steepness for both cardinal directions.
341 if(d > 1)
342 diff[d] += SQR(rfx[i << c][1] - rfx[-(i << c)][1] - rfx[i << c][h] + rfx[-(i << c)][h])
343 + SQR(g);
344 }
345 if((d < 2) || (d & 1))
346 { // output for passes 0, 1, 3, 5
347 // for 0, 1 just use hori/vert, for 3, 5 use best of x/y dir
348 const int d_out = d - ((d > 1) && (diff[d-1] < diff[d]));
349 rfx[0][0] = color[0][d_out] / 2.f;
350 rfx[0][2] = color[1][d_out] / 2.f;
351 rfx += TS * TS;
352 }
353 }
354 }
355
356 /* Interpolate red for blue pixels and vice versa: */
357 const int pad_rb_br = (passes == 1) ? 6 : 5;
358 for(int row = top + pad_rb_br; row < mrow - pad_rb_br; row++)
359 for(int col = left + pad_rb_br; col < mcol - pad_rb_br; col++)
360 {
361 const int f = 2 - FCxtrans(row, col, roi_in, xtrans);
362 if(f == 1) continue;
363 float(*rfx)[3] = &rgb[0][row - top][col - left];
364 const int c = (row - sgrow) % 3 ? TS : 1;
365 const int h = 3 * (c ^ TS ^ 1);
366 for(int d = 0; d < 4; d++, rfx += TS * TS)
367 {
368 const int i = d > 1 || ((d ^ c) & 1) ||
369 ((fabsf(rfx[0][1]-rfx[c][1]) + fabsf(rfx[0][1]-rfx[-c][1])) <
370 2.f*(fabsf(rfx[0][1]-rfx[h][1]) + fabsf(rfx[0][1]-rfx[-h][1]))) ? c:h;
371 rfx[0][f] = (rfx[i][f] + rfx[-i][f] + 2.f * rfx[0][1] - rfx[i][1] - rfx[-i][1]) / 2.f;
372 }
373 }
374
375 /* Fill in red and blue for 2x2 blocks of green: */
376 const int pad_g22 = (passes == 1) ? 8 : 4;
377 for(int row = top + pad_g22; row < mrow - pad_g22; row++)
378 {
379 if((row - sgrow) % 3)
380 for(int col = left + pad_g22; col < mcol - pad_g22; col++)
381 if((col - sgcol) % 3)
382 {
383 float(*rfx)[3] = &rgb[0][row - top][col - left];
384 const short *const hex = hexmap(row,col,allhex);
385 for(int d = 0; d < ndir; d += 2, rfx += TS * TS)
386 if(hex[d] + hex[d + 1])
387 {
388 const float g = 3.f * rfx[0][1] - 2.f * rfx[hex[d]][1] - rfx[hex[d + 1]][1];
389 for(int c = 0; c < 4; c += 2)
390 rfx[0][c] = (g + 2.f * rfx[hex[d]][c] + rfx[hex[d + 1]][c]) / 3.f;
391 }
392 else
393 {
394 const float g = 2.f * rfx[0][1] - rfx[hex[d]][1] - rfx[hex[d + 1]][1];
395 for(int c = 0; c < 4; c += 2)
396 rfx[0][c] = (g + rfx[hex[d]][c] + rfx[hex[d + 1]][c]) / 2.f;
397 }
398 }
399 }
400 } // end of multipass loop
401
402 // jump back to the first set of rgb buffers (this is a nop
403 // unless on the second pass)
404 rgb = (float(*)[TS][TS][3])buffer;
405 // from here on out, mainly are working within the current tile
406 // rather than in reference to the image, so don't offset
407 // mrow/mcol by top/left of tile
408 mrow -= top;
409 mcol -= left;
410
411 /* Convert to perceptual colorspace and differentiate in all directions: */
412 // Original dcraw algorithm uses CIELab as perceptual space
413 // (presumably coming from original AHD) and converts taking
414 // camera matrix into account. Now use YPbPr which requires much
415 // less code and is nearly indistinguishable. It assumes the
416 // camera RGB is roughly linear.
417 for(int d = 0; d < ndir; d++)
418 {
419 const int pad_yuv = (passes == 1) ? 8 : 13;
420 for(int row = pad_yuv; row < mrow - pad_yuv; row++)
421 for(int col = pad_yuv; col < mcol - pad_yuv; col++)
422 {
423 const float *rx = rgb[d][row][col];
424 // use ITU-R BT.2020 YPbPr, which is great, but could use
425 // a better/simpler choice? note that imageop.h provides
426 // dt_iop_RGB_to_YCbCr which uses Rec. 601 conversion,
427 // which appears less good with specular highlights
428 const float y = 0.2627f * rx[0] + 0.6780f * rx[1] + 0.0593f * rx[2];
429 yuv[0][row][col] = y;
430 yuv[1][row][col] = (rx[2] - y) * 0.56433f;
431 yuv[2][row][col] = (rx[0] - y) * 0.67815f;
432 }
433 // Note that f can offset by a column (-1 or +1) and by a row
434 // (-TS or TS). The row-wise offsets cause the undefined
435 // behavior sanitizer to warn of an out of bounds index, but
436 // as yfx is multi-dimensional and there is sufficient
437 // padding, that is not actually so.
438 const int f = dir[d & 3];
439 const int pad_drv = (passes == 1) ? 9 : 14;
440 for(int row = pad_drv; row < mrow - pad_drv; row++)
441 for(int col = pad_drv; col < mcol - pad_drv; col++)
442 {
443 const float(*yfx)[TS][TS] = (float(*)[TS][TS]) & yuv[0][row][col];
444 drv[d][row][col] = SQR(2 * yfx[0][0][0] - yfx[0][0][f] - yfx[0][0][-f])
445 + SQR(2 * yfx[1][0][0] - yfx[1][0][f] - yfx[1][0][-f])
446 + SQR(2 * yfx[2][0][0] - yfx[2][0][f] - yfx[2][0][-f]);
447 }
448 }
449
450 /* Build homogeneity maps from the derivatives: */
451 // The cast is deliberate. ndir is 4 or 8 by construction, but this runs inside an
452 // OpenMP-outlined function where it arrives through the shared closure, so GCC can no
453 // longer see that range and must assume the int -> size_t conversion sign-extends a
454 // negative -- the astronomically large memset bound it reports. Going through unsigned
455 // states the invariant that holds here.
456 memset_zero(homo, sizeof(uint8_t) * (size_t)(unsigned)ndir * TS * TS);
457 const int pad_homo = (passes == 1) ? 10 : 15;
458 for(int row = pad_homo; row < mrow - pad_homo; row++)
459 for(int col = pad_homo; col < mcol - pad_homo; col++)
460 {
461 float tr = FLT_MAX;
462 for(int d = 0; d < ndir; d++)
463 if(tr > drv[d][row][col]) tr = drv[d][row][col];
464 tr *= 8;
465 for(int d = 0; d < ndir; d++)
466 for(int v = -1; v <= 1; v++)
467 for(int h = -1; h <= 1; h++)
468 homo[d][row][col] += ((drv[d][row + v][col + h] <= tr) ? 1 : 0);
469 }
470
471 /* Build 5x5 sum of homogeneity maps for each pixel & direction */
472 for(int d = 0; d < ndir; d++)
473 for(int row = pad_tile; row < mrow - pad_tile; row++)
474 {
475 // start before first column where homo[d][row][col+2] != 0,
476 // so can know v5sum and homosum[d][row][col] will be 0
477 int col = pad_tile-5;
478 uint8_t v5sum[5] = { 0 };
479 homosum[d][row][col] = 0;
480 // calculate by rolling through column sums
481 for(col++; col < mcol - pad_tile; col++)
482 {
483 uint8_t colsum = 0;
484 for(int v = -2; v <= 2; v++) colsum += homo[d][row + v][col + 2];
485 homosum[d][row][col] = homosum[d][row][col - 1] - v5sum[col % 5] + colsum;
486 v5sum[col % 5] = colsum;
487 }
488 }
489
490 /* Average the most homogeneous pixels for the final result: */
491 for(int row = pad_tile; row < mrow - pad_tile; row++)
492 for(int col = pad_tile; col < mcol - pad_tile; col++)
493 {
494 uint8_t hm[8] = { 0 };
495 uint8_t maxval = 0;
496 for(int d = 0; d < ndir; d++)
497 {
498 hm[d] = homosum[d][row][col];
499 maxval = (maxval < hm[d] ? hm[d] : maxval);
500 }
501 maxval -= maxval >> 3;
502 for(int d = 0; d < ndir - 4; d++)
503 {
504 if(hm[d] < hm[d + 4])
505 hm[d] = 0;
506 else if(hm[d] > hm[d + 4])
507 hm[d + 4] = 0;
508 }
509 dt_aligned_pixel_t avg = { 0.0f };
510 for(int d = 0; d < ndir; d++)
511 {
512 if(hm[d] >= maxval)
513 {
514 for(int c = 0; c < 3; c++) avg[c] += rgb[d][row][col][c];
515 avg[3]++;
516 }
517 }
518 for(int c = 0; c < 3; c++)
519 out[4 * (width * (row + top) + col + left) + c] = avg[c]/avg[3];
520 }
521 }
522 }
523
524
525 dt_pixelpipe_cache_free_align(all_buffers);
526}
527
528#undef TS
529
530#define TS 122
531__DT_CLONE_TARGETS__
532static void xtrans_fdc_interpolate(struct dt_iop_module_t *self, float *out, const float *const in,
533 const dt_iop_roi_t *const roi_out, const dt_iop_roi_t *const roi_in,
534 const uint8_t (*const xtrans)[6])
535{
536
537 static const short orth[12] = { 1, 0, 0, 1, -1, 0, 0, -1, 1, 0, 0, 1 },
538 patt[2][16] = { { 0, 1, 0, -1, 2, 0, -1, 0, 1, 1, 1, -1, 0, 0, 0, 0 },
539 { 0, 1, 0, -2, 1, 0, -2, 0, 1, 1, -2, -2, 1, -1, -1, 1 } },
540 dir[4] = { 1, TS, TS + 1, TS - 1 };
541
542 static const float directionality[8] = { 1.0f, 0.0f, 0.5f, 0.5f, 1.0f, 0.0f, 0.5f, 0.5f };
543
544 short allhex[3][3][8];
545 // sgrow/sgcol is the offset in the sensor matrix of the solitary
546 // green pixels (initialized here only to avoid compiler warning)
547 unsigned short sgrow = 0, sgcol = 0;
548
549 const int width = roi_out->width;
550 const int height = roi_out->height;
551 static const int ndir = 4;
552
553 static const float complex Minv[3][8] = {
554 { 1.000000e+00f, 2.500000e-01f - 4.330127e-01f * _Complex_I, -2.500000e-01f - 4.330127e-01f * _Complex_I,
555 -1.000000e+00f, 7.500000e-01f - 1.299038e+00f * _Complex_I, -2.500000e-01f + 4.330127e-01f * _Complex_I,
556 7.500000e-01f + 1.299038e+00f * _Complex_I, 2.500000e-01f + 4.330127e-01f * _Complex_I },
557 { 1.000000e+00f, -2.000000e-01f + 3.464102e-01f * _Complex_I, 2.000000e-01f + 3.464102e-01f * _Complex_I,
558 8.000000e-01f, 0.0f, 2.000000e-01f - 3.464102e-01f * _Complex_I, 0.0f,
559 -2.000000e-01f - 3.464102e-01f * _Complex_I },
560 { 1.000000e+00f, 2.500000e-01f - 4.330127e-01f * _Complex_I, -2.500000e-01f - 4.330127e-01f * _Complex_I,
561 -1.000000e+00f, -7.500000e-01f + 1.299038e+00f * _Complex_I, -2.500000e-01f + 4.330127e-01f * _Complex_I,
562 -7.500000e-01f - 1.299038e+00f * _Complex_I, 2.500000e-01f + 4.330127e-01f * _Complex_I },
563 };
564
565 static const float complex modarr[6][6][8] = {
566 { { 1.000000e+00f + 0.000000e+00f * _Complex_I, 1.000000e+00f + 0.000000e+00f * _Complex_I,
567 1.000000e+00f + 0.000000e+00f * _Complex_I, 1.000000e+00f + 0.000000e+00f * _Complex_I,
568 1.000000e+00f + 0.000000e+00f * _Complex_I, 1.000000e+00f + 0.000000e+00f * _Complex_I,
569 1.000000e+00f + 0.000000e+00f * _Complex_I, 1.000000e+00f + 0.000000e+00f * _Complex_I },
570 { -1.000000e+00f - 1.224647e-16f * _Complex_I, 5.000000e-01f + 8.660254e-01f * _Complex_I,
571 -1.000000e+00f - 1.224647e-16f * _Complex_I, 5.000000e-01f - 8.660254e-01f * _Complex_I,
572 -5.000000e-01f + 8.660254e-01f * _Complex_I, 1.000000e+00f + 0.000000e+00f * _Complex_I,
573 -5.000000e-01f + 8.660254e-01f * _Complex_I, -5.000000e-01f - 8.660254e-01f * _Complex_I },
574 { 1.000000e+00f + 2.449294e-16f * _Complex_I, -5.000000e-01f + 8.660254e-01f * _Complex_I,
575 1.000000e+00f + 2.449294e-16f * _Complex_I, -5.000000e-01f - 8.660254e-01f * _Complex_I,
576 -5.000000e-01f - 8.660254e-01f * _Complex_I, 1.000000e+00f + 0.000000e+00f * _Complex_I,
577 -5.000000e-01f - 8.660254e-01f * _Complex_I, -5.000000e-01f + 8.660254e-01f * _Complex_I },
578 { -1.000000e+00f - 3.673940e-16f * _Complex_I, -1.000000e+00f + 1.224647e-16f * _Complex_I,
579 -1.000000e+00f - 3.673940e-16f * _Complex_I, -1.000000e+00f - 1.224647e-16f * _Complex_I,
580 1.000000e+00f - 2.449294e-16f * _Complex_I, 1.000000e+00f + 0.000000e+00f * _Complex_I,
581 1.000000e+00f - 2.449294e-16f * _Complex_I, 1.000000e+00f + 2.449294e-16f * _Complex_I },
582 { 1.000000e+00f + 4.898587e-16f * _Complex_I, -5.000000e-01f - 8.660254e-01f * _Complex_I,
583 1.000000e+00f + 4.898587e-16f * _Complex_I, -5.000000e-01f + 8.660254e-01f * _Complex_I,
584 -5.000000e-01f + 8.660254e-01f * _Complex_I, 1.000000e+00f + 0.000000e+00f * _Complex_I,
585 -5.000000e-01f + 8.660254e-01f * _Complex_I, -5.000000e-01f - 8.660254e-01f * _Complex_I },
586 { -1.000000e+00f - 6.123234e-16f * _Complex_I, 5.000000e-01f - 8.660254e-01f * _Complex_I,
587 -1.000000e+00f - 6.123234e-16f * _Complex_I, 5.000000e-01f + 8.660254e-01f * _Complex_I,
588 -5.000000e-01f - 8.660254e-01f * _Complex_I, 1.000000e+00f + 0.000000e+00f * _Complex_I,
589 -5.000000e-01f - 8.660254e-01f * _Complex_I, -5.000000e-01f + 8.660254e-01f * _Complex_I } },
590 { { 5.000000e-01f + 8.660254e-01f * _Complex_I, -1.000000e+00f + 1.224647e-16f * _Complex_I,
591 5.000000e-01f - 8.660254e-01f * _Complex_I, -1.000000e+00f + 1.224647e-16f * _Complex_I,
592 1.000000e+00f + 0.000000e+00f * _Complex_I, -5.000000e-01f + 8.660254e-01f * _Complex_I,
593 -5.000000e-01f + 8.660254e-01f * _Complex_I, -5.000000e-01f + 8.660254e-01f * _Complex_I },
594 { -5.000000e-01f - 8.660254e-01f * _Complex_I, -5.000000e-01f - 8.660254e-01f * _Complex_I,
595 -5.000000e-01f + 8.660254e-01f * _Complex_I, -5.000000e-01f + 8.660254e-01f * _Complex_I,
596 -5.000000e-01f + 8.660254e-01f * _Complex_I, -5.000000e-01f + 8.660254e-01f * _Complex_I,
597 -5.000000e-01f - 8.660254e-01f * _Complex_I, 1.000000e+00f + 0.000000e+00f * _Complex_I },
598 { 5.000000e-01f + 8.660254e-01f * _Complex_I, 5.000000e-01f - 8.660254e-01f * _Complex_I,
599 5.000000e-01f - 8.660254e-01f * _Complex_I, 5.000000e-01f + 8.660254e-01f * _Complex_I,
600 -5.000000e-01f - 8.660254e-01f * _Complex_I, -5.000000e-01f + 8.660254e-01f * _Complex_I,
601 1.000000e+00f - 2.449294e-16f * _Complex_I, -5.000000e-01f - 8.660254e-01f * _Complex_I },
602 { -5.000000e-01f - 8.660254e-01f * _Complex_I, 1.000000e+00f - 2.449294e-16f * _Complex_I,
603 -5.000000e-01f + 8.660254e-01f * _Complex_I, 1.000000e+00f + 0.000000e+00f * _Complex_I,
604 1.000000e+00f - 2.449294e-16f * _Complex_I, -5.000000e-01f + 8.660254e-01f * _Complex_I,
605 -5.000000e-01f + 8.660254e-01f * _Complex_I, -5.000000e-01f + 8.660254e-01f * _Complex_I },
606 { 5.000000e-01f + 8.660254e-01f * _Complex_I, 5.000000e-01f + 8.660254e-01f * _Complex_I,
607 5.000000e-01f - 8.660254e-01f * _Complex_I, 5.000000e-01f - 8.660254e-01f * _Complex_I,
608 -5.000000e-01f + 8.660254e-01f * _Complex_I, -5.000000e-01f + 8.660254e-01f * _Complex_I,
609 -5.000000e-01f - 8.660254e-01f * _Complex_I, 1.000000e+00f + 2.449294e-16f * _Complex_I },
610 { -5.000000e-01f - 8.660254e-01f * _Complex_I, -5.000000e-01f + 8.660254e-01f * _Complex_I,
611 -5.000000e-01f + 8.660254e-01f * _Complex_I, -5.000000e-01f - 8.660254e-01f * _Complex_I,
612 -5.000000e-01f - 8.660254e-01f * _Complex_I, -5.000000e-01f + 8.660254e-01f * _Complex_I,
613 1.000000e+00f - 2.266216e-15f * _Complex_I, -5.000000e-01f - 8.660254e-01f * _Complex_I } },
614 { { -5.000000e-01f + 8.660254e-01f * _Complex_I, 1.000000e+00f - 2.449294e-16f * _Complex_I,
615 -5.000000e-01f - 8.660254e-01f * _Complex_I, 1.000000e+00f - 2.449294e-16f * _Complex_I,
616 1.000000e+00f + 0.000000e+00f * _Complex_I, -5.000000e-01f - 8.660254e-01f * _Complex_I,
617 -5.000000e-01f - 8.660254e-01f * _Complex_I, -5.000000e-01f - 8.660254e-01f * _Complex_I },
618 { 5.000000e-01f - 8.660254e-01f * _Complex_I, 5.000000e-01f + 8.660254e-01f * _Complex_I,
619 5.000000e-01f + 8.660254e-01f * _Complex_I, 5.000000e-01f - 8.660254e-01f * _Complex_I,
620 -5.000000e-01f + 8.660254e-01f * _Complex_I, -5.000000e-01f - 8.660254e-01f * _Complex_I,
621 1.000000e+00f - 2.449294e-16f * _Complex_I, -5.000000e-01f + 8.660254e-01f * _Complex_I },
622 { -5.000000e-01f + 8.660254e-01f * _Complex_I, -5.000000e-01f + 8.660254e-01f * _Complex_I,
623 -5.000000e-01f - 8.660254e-01f * _Complex_I, -5.000000e-01f - 8.660254e-01f * _Complex_I,
624 -5.000000e-01f - 8.660254e-01f * _Complex_I, -5.000000e-01f - 8.660254e-01f * _Complex_I,
625 -5.000000e-01f + 8.660254e-01f * _Complex_I, 1.000000e+00f + 0.000000e+00f * _Complex_I },
626 { 5.000000e-01f - 8.660254e-01f * _Complex_I, -1.000000e+00f + 3.673940e-16f * _Complex_I,
627 5.000000e-01f + 8.660254e-01f * _Complex_I, -1.000000e+00f + 1.224647e-16f * _Complex_I,
628 1.000000e+00f - 2.449294e-16f * _Complex_I, -5.000000e-01f - 8.660254e-01f * _Complex_I,
629 -5.000000e-01f - 8.660254e-01f * _Complex_I, -5.000000e-01f - 8.660254e-01f * _Complex_I },
630 { -5.000000e-01f + 8.660254e-01f * _Complex_I, -5.000000e-01f - 8.660254e-01f * _Complex_I,
631 -5.000000e-01f - 8.660254e-01f * _Complex_I, -5.000000e-01f + 8.660254e-01f * _Complex_I,
632 -5.000000e-01f + 8.660254e-01f * _Complex_I, -5.000000e-01f - 8.660254e-01f * _Complex_I,
633 1.000000e+00f - 4.898587e-16f * _Complex_I, -5.000000e-01f + 8.660254e-01f * _Complex_I },
634 { 5.000000e-01f - 8.660254e-01f * _Complex_I, 5.000000e-01f - 8.660254e-01f * _Complex_I,
635 5.000000e-01f + 8.660254e-01f * _Complex_I, 5.000000e-01f + 8.660254e-01f * _Complex_I,
636 -5.000000e-01f - 8.660254e-01f * _Complex_I, -5.000000e-01f - 8.660254e-01f * _Complex_I,
637 -5.000000e-01f + 8.660254e-01f * _Complex_I, 1.000000e+00f + 1.133108e-15f * _Complex_I } },
638 { { -1.000000e+00f + 1.224647e-16f * _Complex_I, -1.000000e+00f + 3.673940e-16f * _Complex_I,
639 -1.000000e+00f - 1.224647e-16f * _Complex_I, -1.000000e+00f + 3.673940e-16f * _Complex_I,
640 1.000000e+00f + 0.000000e+00f * _Complex_I, 1.000000e+00f - 2.449294e-16f * _Complex_I,
641 1.000000e+00f - 2.449294e-16f * _Complex_I, 1.000000e+00f - 2.449294e-16f * _Complex_I },
642 { 1.000000e+00f + 0.000000e+00f * _Complex_I, -5.000000e-01f - 8.660254e-01f * _Complex_I,
643 1.000000e+00f + 2.449294e-16f * _Complex_I, -5.000000e-01f + 8.660254e-01f * _Complex_I,
644 -5.000000e-01f + 8.660254e-01f * _Complex_I, 1.000000e+00f - 2.449294e-16f * _Complex_I,
645 -5.000000e-01f + 8.660254e-01f * _Complex_I, -5.000000e-01f - 8.660254e-01f * _Complex_I },
646 { -1.000000e+00f - 1.224647e-16f * _Complex_I, 5.000000e-01f - 8.660254e-01f * _Complex_I,
647 -1.000000e+00f - 3.673940e-16f * _Complex_I, 5.000000e-01f + 8.660254e-01f * _Complex_I,
648 -5.000000e-01f - 8.660254e-01f * _Complex_I, 1.000000e+00f - 2.449294e-16f * _Complex_I,
649 -5.000000e-01f - 8.660254e-01f * _Complex_I, -5.000000e-01f + 8.660254e-01f * _Complex_I },
650 { 1.000000e+00f + 2.449294e-16f * _Complex_I, 1.000000e+00f - 4.898587e-16f * _Complex_I,
651 1.000000e+00f + 4.898587e-16f * _Complex_I, 1.000000e+00f - 2.449294e-16f * _Complex_I,
652 1.000000e+00f - 2.449294e-16f * _Complex_I, 1.000000e+00f - 2.449294e-16f * _Complex_I,
653 1.000000e+00f - 4.898587e-16f * _Complex_I, 1.000000e+00f + 0.000000e+00f * _Complex_I },
654 { -1.000000e+00f - 3.673940e-16f * _Complex_I, 5.000000e-01f + 8.660254e-01f * _Complex_I,
655 -1.000000e+00f - 6.123234e-16f * _Complex_I, 5.000000e-01f - 8.660254e-01f * _Complex_I,
656 -5.000000e-01f + 8.660254e-01f * _Complex_I, 1.000000e+00f - 2.449294e-16f * _Complex_I,
657 -5.000000e-01f + 8.660254e-01f * _Complex_I, -5.000000e-01f - 8.660254e-01f * _Complex_I },
658 { 1.000000e+00f + 4.898587e-16f * _Complex_I, -5.000000e-01f + 8.660254e-01f * _Complex_I,
659 1.000000e+00f + 7.347881e-16f * _Complex_I, -5.000000e-01f - 8.660254e-01f * _Complex_I,
660 -5.000000e-01f - 8.660254e-01f * _Complex_I, 1.000000e+00f - 2.449294e-16f * _Complex_I,
661 -5.000000e-01f - 8.660254e-01f * _Complex_I, -5.000000e-01f + 8.660254e-01f * _Complex_I } },
662 { { -5.000000e-01f - 8.660254e-01f * _Complex_I, 1.000000e+00f - 4.898587e-16f * _Complex_I,
663 -5.000000e-01f + 8.660254e-01f * _Complex_I, 1.000000e+00f - 4.898587e-16f * _Complex_I,
664 1.000000e+00f + 0.000000e+00f * _Complex_I, -5.000000e-01f + 8.660254e-01f * _Complex_I,
665 -5.000000e-01f + 8.660254e-01f * _Complex_I, -5.000000e-01f + 8.660254e-01f * _Complex_I },
666 { 5.000000e-01f + 8.660254e-01f * _Complex_I, 5.000000e-01f + 8.660254e-01f * _Complex_I,
667 5.000000e-01f - 8.660254e-01f * _Complex_I, 5.000000e-01f - 8.660254e-01f * _Complex_I,
668 -5.000000e-01f + 8.660254e-01f * _Complex_I, -5.000000e-01f + 8.660254e-01f * _Complex_I,
669 -5.000000e-01f - 8.660254e-01f * _Complex_I, 1.000000e+00f - 2.449294e-16f * _Complex_I },
670 { -5.000000e-01f - 8.660254e-01f * _Complex_I, -5.000000e-01f + 8.660254e-01f * _Complex_I,
671 -5.000000e-01f + 8.660254e-01f * _Complex_I, -5.000000e-01f - 8.660254e-01f * _Complex_I,
672 -5.000000e-01f - 8.660254e-01f * _Complex_I, -5.000000e-01f + 8.660254e-01f * _Complex_I,
673 1.000000e+00f - 4.898587e-16f * _Complex_I, -5.000000e-01f - 8.660254e-01f * _Complex_I },
674 { 5.000000e-01f + 8.660254e-01f * _Complex_I, -1.000000e+00f + 6.123234e-16f * _Complex_I,
675 5.000000e-01f - 8.660254e-01f * _Complex_I, -1.000000e+00f + 3.673940e-16f * _Complex_I,
676 1.000000e+00f - 2.449294e-16f * _Complex_I, -5.000000e-01f + 8.660254e-01f * _Complex_I,
677 -5.000000e-01f + 8.660254e-01f * _Complex_I, -5.000000e-01f + 8.660254e-01f * _Complex_I },
678 { -5.000000e-01f - 8.660254e-01f * _Complex_I, -5.000000e-01f - 8.660254e-01f * _Complex_I,
679 -5.000000e-01f + 8.660254e-01f * _Complex_I, -5.000000e-01f + 8.660254e-01f * _Complex_I,
680 -5.000000e-01f + 8.660254e-01f * _Complex_I, -5.000000e-01f + 8.660254e-01f * _Complex_I,
681 -5.000000e-01f - 8.660254e-01f * _Complex_I, 1.000000e+00f + 0.000000e+00f * _Complex_I },
682 { 5.000000e-01f + 8.660254e-01f * _Complex_I, 5.000000e-01f - 8.660254e-01f * _Complex_I,
683 5.000000e-01f - 8.660254e-01f * _Complex_I, 5.000000e-01f + 8.660254e-01f * _Complex_I,
684 -5.000000e-01f - 8.660254e-01f * _Complex_I, -5.000000e-01f + 8.660254e-01f * _Complex_I,
685 1.000000e+00f - 7.347881e-16f * _Complex_I, -5.000000e-01f - 8.660254e-01f * _Complex_I } },
686 { { 5.000000e-01f - 8.660254e-01f * _Complex_I, -1.000000e+00f + 6.123234e-16f * _Complex_I,
687 5.000000e-01f + 8.660254e-01f * _Complex_I, -1.000000e+00f + 6.123234e-16f * _Complex_I,
688 1.000000e+00f + 0.000000e+00f * _Complex_I, -5.000000e-01f - 8.660254e-01f * _Complex_I,
689 -5.000000e-01f - 8.660254e-01f * _Complex_I, -5.000000e-01f - 8.660254e-01f * _Complex_I },
690 { -5.000000e-01f + 8.660254e-01f * _Complex_I, -5.000000e-01f - 8.660254e-01f * _Complex_I,
691 -5.000000e-01f - 8.660254e-01f * _Complex_I, -5.000000e-01f + 8.660254e-01f * _Complex_I,
692 -5.000000e-01f + 8.660254e-01f * _Complex_I, -5.000000e-01f - 8.660254e-01f * _Complex_I,
693 1.000000e+00f - 2.266216e-15f * _Complex_I, -5.000000e-01f + 8.660254e-01f * _Complex_I },
694 { 5.000000e-01f - 8.660254e-01f * _Complex_I, 5.000000e-01f - 8.660254e-01f * _Complex_I,
695 5.000000e-01f + 8.660254e-01f * _Complex_I, 5.000000e-01f + 8.660254e-01f * _Complex_I,
696 -5.000000e-01f - 8.660254e-01f * _Complex_I, -5.000000e-01f - 8.660254e-01f * _Complex_I,
697 -5.000000e-01f + 8.660254e-01f * _Complex_I, 1.000000e+00f - 1.133108e-15f * _Complex_I },
698 { -5.000000e-01f + 8.660254e-01f * _Complex_I, 1.000000e+00f - 7.347881e-16f * _Complex_I,
699 -5.000000e-01f - 8.660254e-01f * _Complex_I, 1.000000e+00f - 4.898587e-16f * _Complex_I,
700 1.000000e+00f - 2.449294e-16f * _Complex_I, -5.000000e-01f - 8.660254e-01f * _Complex_I,
701 -5.000000e-01f - 8.660254e-01f * _Complex_I, -5.000000e-01f - 8.660254e-01f * _Complex_I },
702 { 5.000000e-01f - 8.660254e-01f * _Complex_I, 5.000000e-01f + 8.660254e-01f * _Complex_I,
703 5.000000e-01f + 8.660254e-01f * _Complex_I, 5.000000e-01f - 8.660254e-01f * _Complex_I,
704 -5.000000e-01f + 8.660254e-01f * _Complex_I, -5.000000e-01f - 8.660254e-01f * _Complex_I,
705 1.000000e+00f - 7.347881e-16f * _Complex_I, -5.000000e-01f + 8.660254e-01f * _Complex_I },
706 { -5.000000e-01f + 8.660254e-01f * _Complex_I, -5.000000e-01f + 8.660254e-01f * _Complex_I,
707 -5.000000e-01f - 8.660254e-01f * _Complex_I, -5.000000e-01f - 8.660254e-01f * _Complex_I,
708 -5.000000e-01f - 8.660254e-01f * _Complex_I, -5.000000e-01f - 8.660254e-01f * _Complex_I,
709 -5.000000e-01f + 8.660254e-01f * _Complex_I, 1.000000e+00f + 0.000000e+00f * _Complex_I } },
710 };
711
712 static const float complex harr[4][13][13]
713 = { { { 1.326343e-03f - 1.299441e-18f * _Complex_I, 7.091837e-04f - 1.228342e-03f * _Complex_I,
714 -6.278557e-04f - 1.087478e-03f * _Complex_I, -1.157216e-03f + 9.920263e-19f * _Complex_I,
715 -4.887166e-04f + 8.464820e-04f * _Complex_I, 5.758687e-04f + 9.974338e-04f * _Complex_I,
716 1.225183e-03f - 9.002496e-19f * _Complex_I, 5.758687e-04f - 9.974338e-04f * _Complex_I,
717 -4.887166e-04f - 8.464820e-04f * _Complex_I, -1.157216e-03f + 7.085902e-19f * _Complex_I,
718 -6.278557e-04f + 1.087478e-03f * _Complex_I, 7.091837e-04f + 1.228342e-03f * _Complex_I,
719 1.326343e-03f - 6.497206e-19f * _Complex_I },
720 { -1.980815e-03f + 1.698059e-18f * _Complex_I, -1.070384e-03f + 1.853959e-03f * _Complex_I,
721 7.924697e-04f + 1.372598e-03f * _Complex_I, 1.876584e-03f - 1.378892e-18f * _Complex_I,
722 1.225866e-03f - 2.123262e-03f * _Complex_I, -1.569320e-03f - 2.718142e-03f * _Complex_I,
723 -3.273971e-03f + 2.004729e-18f * _Complex_I, -1.569320e-03f + 2.718142e-03f * _Complex_I,
724 1.225866e-03f + 2.123262e-03f * _Complex_I, 1.876584e-03f - 9.192611e-19f * _Complex_I,
725 7.924697e-04f - 1.372598e-03f * _Complex_I, -1.070384e-03f - 1.853959e-03f * _Complex_I,
726 -1.980815e-03f + 7.277398e-19f * _Complex_I },
727 { 1.457023e-03f - 1.070603e-18f * _Complex_I, 8.487143e-04f - 1.470016e-03f * _Complex_I,
728 -6.873776e-04f - 1.190573e-03f * _Complex_I, -2.668335e-03f + 1.633884e-18f * _Complex_I,
729 -2.459813e-03f + 4.260521e-03f * _Complex_I, 3.238772e-03f + 5.609717e-03f * _Complex_I,
730 7.074895e-03f - 3.465699e-18f * _Complex_I, 3.238772e-03f - 5.609717e-03f * _Complex_I,
731 -2.459813e-03f - 4.260521e-03f * _Complex_I, -2.668335e-03f + 9.803302e-19f * _Complex_I,
732 -6.873776e-04f + 1.190573e-03f * _Complex_I, 8.487143e-04f + 1.470016e-03f * _Complex_I,
733 1.457023e-03f - 3.568678e-19f * _Complex_I },
734 { -1.017660e-03f + 6.231370e-19f * _Complex_I, -5.415171e-04f + 9.379351e-04f * _Complex_I,
735 7.255109e-04f + 1.256622e-03f * _Complex_I, 3.699792e-03f - 1.812375e-18f * _Complex_I,
736 4.090356e-03f - 7.084704e-03f * _Complex_I, -6.006283e-03f - 1.040319e-02f * _Complex_I,
737 -1.391431e-02f + 5.112034e-18f * _Complex_I, -6.006283e-03f + 1.040319e-02f * _Complex_I,
738 4.090356e-03f + 7.084704e-03f * _Complex_I, 3.699792e-03f - 9.061876e-19f * _Complex_I,
739 7.255109e-04f - 1.256622e-03f * _Complex_I, -5.415171e-04f - 9.379351e-04f * _Complex_I,
740 -1.017660e-03f + 1.246274e-19f * _Complex_I },
741 { 9.198983e-04f - 4.506202e-19f * _Complex_I, 6.815900e-04f - 1.180548e-03f * _Complex_I,
742 -1.287335e-03f - 2.229729e-03f * _Complex_I, -5.023856e-03f + 1.845735e-18f * _Complex_I,
743 -5.499048e-03f + 9.524630e-03f * _Complex_I, 9.797672e-03f + 1.697006e-02f * _Complex_I,
744 2.504795e-02f - 6.134977e-18f * _Complex_I, 9.797672e-03f - 1.697006e-02f * _Complex_I,
745 -5.499048e-03f - 9.524630e-03f * _Complex_I, -5.023856e-03f + 6.152449e-19f * _Complex_I,
746 -1.287335e-03f + 2.229729e-03f * _Complex_I, 6.815900e-04f + 1.180548e-03f * _Complex_I,
747 9.198983e-04f + 0.000000e+00f * _Complex_I },
748 { -7.972663e-04f + 2.929109e-19f * _Complex_I, -1.145605e-03f + 1.984247e-03f * _Complex_I,
749 1.983334e-03f + 3.435235e-03f * _Complex_I, 6.730096e-03f - 1.648398e-18f * _Complex_I,
750 6.782033e-03f - 1.174683e-02f * _Complex_I, -1.392077e-02f - 2.411147e-02f * _Complex_I,
751 -3.906939e-02f + 4.784620e-18f * _Complex_I, -1.392077e-02f + 2.411147e-02f * _Complex_I,
752 6.782033e-03f + 1.174683e-02f * _Complex_I, 6.730096e-03f + 0.000000e+00f * _Complex_I,
753 1.983334e-03f - 3.435235e-03f * _Complex_I, -1.145605e-03f - 1.984247e-03f * _Complex_I,
754 -7.972663e-04f - 9.763696e-20f * _Complex_I },
755 { 8.625458e-04f - 2.112628e-19f * _Complex_I, 1.431113e-03f - 2.478760e-03f * _Complex_I,
756 -2.310309e-03f - 4.001572e-03f * _Complex_I, -7.706486e-03f + 9.437723e-19f * _Complex_I,
757 -7.220186e-03f + 1.250573e-02f * _Complex_I, 1.587118e-02f + 2.748969e-02f * _Complex_I,
758 4.765675e-02f + 0.000000e+00f * _Complex_I, 1.587118e-02f - 2.748969e-02f * _Complex_I,
759 -7.220186e-03f - 1.250573e-02f * _Complex_I, -7.706486e-03f - 9.437723e-19f * _Complex_I,
760 -2.310309e-03f + 4.001572e-03f * _Complex_I, 1.431113e-03f + 2.478760e-03f * _Complex_I,
761 8.625458e-04f + 2.112628e-19f * _Complex_I },
762 { -7.972663e-04f + 9.763696e-20f * _Complex_I, -1.145605e-03f + 1.984247e-03f * _Complex_I,
763 1.983334e-03f + 3.435235e-03f * _Complex_I, 6.730096e-03f + 0.000000e+00f * _Complex_I,
764 6.782033e-03f - 1.174683e-02f * _Complex_I, -1.392077e-02f - 2.411147e-02f * _Complex_I,
765 -3.906939e-02f - 4.784620e-18f * _Complex_I, -1.392077e-02f + 2.411147e-02f * _Complex_I,
766 6.782033e-03f + 1.174683e-02f * _Complex_I, 6.730096e-03f + 1.648398e-18f * _Complex_I,
767 1.983334e-03f - 3.435235e-03f * _Complex_I, -1.145605e-03f - 1.984247e-03f * _Complex_I,
768 -7.972663e-04f - 2.929109e-19f * _Complex_I },
769 { 9.198983e-04f + 0.000000e+00f * _Complex_I, 6.815900e-04f - 1.180548e-03f * _Complex_I,
770 -1.287335e-03f - 2.229729e-03f * _Complex_I, -5.023856e-03f - 6.152449e-19f * _Complex_I,
771 -5.499048e-03f + 9.524630e-03f * _Complex_I, 9.797672e-03f + 1.697006e-02f * _Complex_I,
772 2.504795e-02f + 6.134977e-18f * _Complex_I, 9.797672e-03f - 1.697006e-02f * _Complex_I,
773 -5.499048e-03f - 9.524630e-03f * _Complex_I, -5.023856e-03f - 1.845735e-18f * _Complex_I,
774 -1.287335e-03f + 2.229729e-03f * _Complex_I, 6.815900e-04f + 1.180548e-03f * _Complex_I,
775 9.198983e-04f + 4.506202e-19f * _Complex_I },
776 { -1.017660e-03f - 1.246274e-19f * _Complex_I, -5.415171e-04f + 9.379351e-04f * _Complex_I,
777 7.255109e-04f + 1.256622e-03f * _Complex_I, 3.699792e-03f + 9.061876e-19f * _Complex_I,
778 4.090356e-03f - 7.084704e-03f * _Complex_I, -6.006283e-03f - 1.040319e-02f * _Complex_I,
779 -1.391431e-02f - 5.112034e-18f * _Complex_I, -6.006283e-03f + 1.040319e-02f * _Complex_I,
780 4.090356e-03f + 7.084704e-03f * _Complex_I, 3.699792e-03f + 1.812375e-18f * _Complex_I,
781 7.255109e-04f - 1.256622e-03f * _Complex_I, -5.415171e-04f - 9.379351e-04f * _Complex_I,
782 -1.017660e-03f - 6.231370e-19f * _Complex_I },
783 { 1.457023e-03f + 3.568678e-19f * _Complex_I, 8.487143e-04f - 1.470016e-03f * _Complex_I,
784 -6.873776e-04f - 1.190573e-03f * _Complex_I, -2.668335e-03f - 9.803302e-19f * _Complex_I,
785 -2.459813e-03f + 4.260521e-03f * _Complex_I, 3.238772e-03f + 5.609717e-03f * _Complex_I,
786 7.074895e-03f + 3.465699e-18f * _Complex_I, 3.238772e-03f - 5.609717e-03f * _Complex_I,
787 -2.459813e-03f - 4.260521e-03f * _Complex_I, -2.668335e-03f - 1.633884e-18f * _Complex_I,
788 -6.873776e-04f + 1.190573e-03f * _Complex_I, 8.487143e-04f + 1.470016e-03f * _Complex_I,
789 1.457023e-03f + 1.070603e-18f * _Complex_I },
790 { -1.980815e-03f - 7.277398e-19f * _Complex_I, -1.070384e-03f + 1.853959e-03f * _Complex_I,
791 7.924697e-04f + 1.372598e-03f * _Complex_I, 1.876584e-03f + 9.192611e-19f * _Complex_I,
792 1.225866e-03f - 2.123262e-03f * _Complex_I, -1.569320e-03f - 2.718142e-03f * _Complex_I,
793 -3.273971e-03f - 2.004729e-18f * _Complex_I, -1.569320e-03f + 2.718142e-03f * _Complex_I,
794 1.225866e-03f + 2.123262e-03f * _Complex_I, 1.876584e-03f + 1.378892e-18f * _Complex_I,
795 7.924697e-04f - 1.372598e-03f * _Complex_I, -1.070384e-03f - 1.853959e-03f * _Complex_I,
796 -1.980815e-03f - 1.698059e-18f * _Complex_I },
797 { 1.326343e-03f + 6.497206e-19f * _Complex_I, 7.091837e-04f - 1.228342e-03f * _Complex_I,
798 -6.278557e-04f - 1.087478e-03f * _Complex_I, -1.157216e-03f - 7.085902e-19f * _Complex_I,
799 -4.887166e-04f + 8.464820e-04f * _Complex_I, 5.758687e-04f + 9.974338e-04f * _Complex_I,
800 1.225183e-03f + 9.002496e-19f * _Complex_I, 5.758687e-04f - 9.974338e-04f * _Complex_I,
801 -4.887166e-04f - 8.464820e-04f * _Complex_I, -1.157216e-03f - 9.920263e-19f * _Complex_I,
802 -6.278557e-04f + 1.087478e-03f * _Complex_I, 7.091837e-04f + 1.228342e-03f * _Complex_I,
803 1.326343e-03f + 1.299441e-18f * _Complex_I } },
804 { { 9.129120e-04f - 8.943958e-19f * _Complex_I, -5.925973e-04f - 1.026409e-03f * _Complex_I,
805 -5.989682e-04f + 1.037443e-03f * _Complex_I, 1.158755e-03f - 8.514393e-19f * _Complex_I,
806 -8.992493e-04f - 1.557545e-03f * _Complex_I, -1.283187e-03f + 2.222546e-03f * _Complex_I,
807 2.730635e-03f - 1.337625e-18f * _Complex_I, -1.283187e-03f - 2.222546e-03f * _Complex_I,
808 -8.992493e-04f + 1.557545e-03f * _Complex_I, 1.158755e-03f - 2.838131e-19f * _Complex_I,
809 -5.989682e-04f - 1.037443e-03f * _Complex_I, -5.925973e-04f + 1.026409e-03f * _Complex_I,
810 9.129120e-04f + 0.000000e+00f * _Complex_I },
811 { -5.588854e-04f - 9.680179e-04f * _Complex_I, -6.474856e-04f + 1.121478e-03f * _Complex_I,
812 1.536588e-03f - 1.129066e-18f * _Complex_I, -9.123802e-04f - 1.580289e-03f * _Complex_I,
813 -1.541434e-03f + 2.669842e-03f * _Complex_I, 4.379825e-03f - 9.925627e-18f * _Complex_I,
814 -2.394173e-03f - 4.146830e-03f * _Complex_I, -2.189912e-03f + 3.793039e-03f * _Complex_I,
815 3.082869e-03f - 3.493222e-18f * _Complex_I, -9.123802e-04f - 1.580289e-03f * _Complex_I,
816 -7.682939e-04f + 1.330724e-03f * _Complex_I, 1.294971e-03f + 0.000000e+00f * _Complex_I,
817 -5.588854e-04f - 9.680179e-04f * _Complex_I },
818 { -5.883876e-04f + 1.019117e-03f * _Complex_I, 1.714796e-03f - 1.260012e-18f * _Complex_I,
819 -1.180365e-03f - 2.044451e-03f * _Complex_I, -1.483082e-03f + 2.568774e-03f * _Complex_I,
820 4.933362e-03f - 2.416651e-18f * _Complex_I, -3.296542e-03f - 5.709779e-03f * _Complex_I,
821 -3.546477e-03f + 6.142678e-03f * _Complex_I, 6.593085e-03f - 1.614840e-18f * _Complex_I,
822 -2.466681e-03f - 4.272417e-03f * _Complex_I, -1.483082e-03f + 2.568774e-03f * _Complex_I,
823 2.360729e-03f + 0.000000e+00f * _Complex_I, -8.573982e-04f - 1.485057e-03f * _Complex_I,
824 -5.883876e-04f + 1.019117e-03f * _Complex_I },
825 { 1.483526e-03f - 1.090077e-18f * _Complex_I, -1.074793e-03f - 1.861596e-03f * _Complex_I,
826 -1.447448e-03f + 2.507053e-03f * _Complex_I, 3.952416e-03f - 1.936126e-18f * _Complex_I,
827 -3.496688e-03f - 6.056441e-03f * _Complex_I, -4.898024e-03f + 8.483627e-03f * _Complex_I,
828 1.070518e-02f - 2.622012e-18f * _Complex_I, -4.898024e-03f - 8.483627e-03f * _Complex_I,
829 -3.496688e-03f + 6.056441e-03f * _Complex_I, 3.952416e-03f + 0.000000e+00f * _Complex_I,
830 -1.447448e-03f - 2.507053e-03f * _Complex_I, -1.074793e-03f + 1.861596e-03f * _Complex_I,
831 1.483526e-03f + 3.633590e-19f * _Complex_I },
832 { -9.966429e-04f - 1.726236e-03f * _Complex_I, -1.478281e-03f + 2.560458e-03f * _Complex_I,
833 4.306274e-03f - 2.109466e-18f * _Complex_I, -3.294955e-03f - 5.707029e-03f * _Complex_I,
834 -5.436890e-03f + 9.416970e-03f * _Complex_I, 1.556418e-02f - 3.812124e-18f * _Complex_I,
835 -8.842875e-03f - 1.531631e-02f * _Complex_I, -7.782088e-03f + 1.347897e-02f * _Complex_I,
836 1.087378e-02f + 0.000000e+00f * _Complex_I, -3.294955e-03f - 5.707029e-03f * _Complex_I,
837 -2.153137e-03f + 3.729342e-03f * _Complex_I, 2.956562e-03f + 3.350104e-18f * _Complex_I,
838 -9.966429e-04f - 1.726236e-03f * _Complex_I },
839 { -1.291288e-03f + 2.236576e-03f * _Complex_I, 3.942788e-03f - 8.935208e-18f * _Complex_I,
840 -2.798347e-03f - 4.846880e-03f * _Complex_I, -4.448869e-03f + 7.705666e-03f * _Complex_I,
841 1.522441e-02f - 3.728906e-18f * _Complex_I, -1.175443e-02f - 2.035927e-02f * _Complex_I,
842 -1.417872e-02f + 2.455826e-02f * _Complex_I, 2.350886e-02f + 0.000000e+00f * _Complex_I,
843 -7.612206e-03f - 1.318473e-02f * _Complex_I, -4.448869e-03f + 7.705666e-03f * _Complex_I,
844 5.596695e-03f + 1.370795e-18f * _Complex_I, -1.971394e-03f - 3.414555e-03f * _Complex_I,
845 -1.291288e-03f + 2.236576e-03f * _Complex_I },
846 { 2.779286e-03f - 1.361458e-18f * _Complex_I, -2.194126e-03f - 3.800338e-03f * _Complex_I,
847 -3.057720e-03f + 5.296126e-03f * _Complex_I, 9.725261e-03f - 2.382002e-18f * _Complex_I,
848 -8.649261e-03f - 1.498096e-02f * _Complex_I, -1.417667e-02f + 2.455472e-02f * _Complex_I,
849 3.552610e-02f + 0.000000e+00f * _Complex_I, -1.417667e-02f - 2.455472e-02f * _Complex_I,
850 -8.649261e-03f + 1.498096e-02f * _Complex_I, 9.725261e-03f + 2.382002e-18f * _Complex_I,
851 -3.057720e-03f - 5.296126e-03f * _Complex_I, -2.194126e-03f + 3.800338e-03f * _Complex_I,
852 2.779286e-03f + 1.361458e-18f * _Complex_I },
853 { -1.291288e-03f - 2.236576e-03f * _Complex_I, -1.971394e-03f + 3.414555e-03f * _Complex_I,
854 5.596695e-03f - 1.370795e-18f * _Complex_I, -4.448869e-03f - 7.705666e-03f * _Complex_I,
855 -7.612206e-03f + 1.318473e-02f * _Complex_I, 2.350886e-02f + 0.000000e+00f * _Complex_I,
856 -1.417872e-02f - 2.455826e-02f * _Complex_I, -1.175443e-02f + 2.035927e-02f * _Complex_I,
857 1.522441e-02f + 3.728906e-18f * _Complex_I, -4.448869e-03f - 7.705666e-03f * _Complex_I,
858 -2.798347e-03f + 4.846880e-03f * _Complex_I, 3.942788e-03f + 8.935208e-18f * _Complex_I,
859 -1.291288e-03f - 2.236576e-03f * _Complex_I },
860 { -9.966429e-04f + 1.726236e-03f * _Complex_I, 2.956562e-03f - 3.350104e-18f * _Complex_I,
861 -2.153137e-03f - 3.729342e-03f * _Complex_I, -3.294955e-03f + 5.707029e-03f * _Complex_I,
862 1.087378e-02f + 0.000000e+00f * _Complex_I, -7.782088e-03f - 1.347897e-02f * _Complex_I,
863 -8.842875e-03f + 1.531631e-02f * _Complex_I, 1.556418e-02f + 3.812124e-18f * _Complex_I,
864 -5.436890e-03f - 9.416970e-03f * _Complex_I, -3.294955e-03f + 5.707029e-03f * _Complex_I,
865 4.306274e-03f + 2.109466e-18f * _Complex_I, -1.478281e-03f - 2.560458e-03f * _Complex_I,
866 -9.966429e-04f + 1.726236e-03f * _Complex_I },
867 { 1.483526e-03f - 3.633590e-19f * _Complex_I, -1.074793e-03f - 1.861596e-03f * _Complex_I,
868 -1.447448e-03f + 2.507053e-03f * _Complex_I, 3.952416e-03f + 0.000000e+00f * _Complex_I,
869 -3.496688e-03f - 6.056441e-03f * _Complex_I, -4.898024e-03f + 8.483627e-03f * _Complex_I,
870 1.070518e-02f + 2.622012e-18f * _Complex_I, -4.898024e-03f - 8.483627e-03f * _Complex_I,
871 -3.496688e-03f + 6.056441e-03f * _Complex_I, 3.952416e-03f + 1.936126e-18f * _Complex_I,
872 -1.447448e-03f - 2.507053e-03f * _Complex_I, -1.074793e-03f + 1.861596e-03f * _Complex_I,
873 1.483526e-03f + 1.090077e-18f * _Complex_I },
874 { -5.883876e-04f - 1.019117e-03f * _Complex_I, -8.573982e-04f + 1.485057e-03f * _Complex_I,
875 2.360729e-03f + 0.000000e+00f * _Complex_I, -1.483082e-03f - 2.568774e-03f * _Complex_I,
876 -2.466681e-03f + 4.272417e-03f * _Complex_I, 6.593085e-03f + 1.614840e-18f * _Complex_I,
877 -3.546477e-03f - 6.142678e-03f * _Complex_I, -3.296542e-03f + 5.709779e-03f * _Complex_I,
878 4.933362e-03f + 2.416651e-18f * _Complex_I, -1.483082e-03f - 2.568774e-03f * _Complex_I,
879 -1.180365e-03f + 2.044451e-03f * _Complex_I, 1.714796e-03f + 1.260012e-18f * _Complex_I,
880 -5.883876e-04f - 1.019117e-03f * _Complex_I },
881 { -5.588854e-04f + 9.680179e-04f * _Complex_I, 1.294971e-03f + 0.000000e+00f * _Complex_I,
882 -7.682939e-04f - 1.330724e-03f * _Complex_I, -9.123802e-04f + 1.580289e-03f * _Complex_I,
883 3.082869e-03f + 3.493222e-18f * _Complex_I, -2.189912e-03f - 3.793039e-03f * _Complex_I,
884 -2.394173e-03f + 4.146830e-03f * _Complex_I, 4.379825e-03f + 9.925627e-18f * _Complex_I,
885 -1.541434e-03f - 2.669842e-03f * _Complex_I, -9.123802e-04f + 1.580289e-03f * _Complex_I,
886 1.536588e-03f + 1.129066e-18f * _Complex_I, -6.474856e-04f - 1.121478e-03f * _Complex_I,
887 -5.588854e-04f + 9.680179e-04f * _Complex_I },
888 { 9.129120e-04f + 0.000000e+00f * _Complex_I, -5.925973e-04f - 1.026409e-03f * _Complex_I,
889 -5.989682e-04f + 1.037443e-03f * _Complex_I, 1.158755e-03f + 2.838131e-19f * _Complex_I,
890 -8.992493e-04f - 1.557545e-03f * _Complex_I, -1.283187e-03f + 2.222546e-03f * _Complex_I,
891 2.730635e-03f + 1.337625e-18f * _Complex_I, -1.283187e-03f - 2.222546e-03f * _Complex_I,
892 -8.992493e-04f + 1.557545e-03f * _Complex_I, 1.158755e-03f + 8.514393e-19f * _Complex_I,
893 -5.989682e-04f - 1.037443e-03f * _Complex_I, -5.925973e-04f + 1.026409e-03f * _Complex_I,
894 9.129120e-04f + 8.943958e-19f * _Complex_I } },
895 { { 8.228091e-04f + 0.000000e+00f * _Complex_I, -5.365069e-04f + 9.292572e-04f * _Complex_I,
896 -6.011501e-04f - 1.041223e-03f * _Complex_I, 1.249890e-03f - 3.061346e-19f * _Complex_I,
897 -7.632708e-04f + 1.322024e-03f * _Complex_I, -9.846035e-04f - 1.705383e-03f * _Complex_I,
898 2.080486e-03f - 1.019144e-18f * _Complex_I, -9.846035e-04f + 1.705383e-03f * _Complex_I,
899 -7.632708e-04f - 1.322024e-03f * _Complex_I, 1.249890e-03f - 9.184039e-19f * _Complex_I,
900 -6.011501e-04f + 1.041223e-03f * _Complex_I, -5.365069e-04f - 9.292572e-04f * _Complex_I,
901 8.228091e-04f - 8.061204e-19f * _Complex_I },
902 { -5.616336e-04f - 9.727779e-04f * _Complex_I, 1.382894e-03f + 0.000000e+00f * _Complex_I,
903 -8.694311e-04f + 1.505899e-03f * _Complex_I, -9.721139e-04f - 1.683751e-03f * _Complex_I,
904 2.446785e-03f - 2.772471e-18f * _Complex_I, -1.605471e-03f + 2.780758e-03f * _Complex_I,
905 -1.832781e-03f - 3.174469e-03f * _Complex_I, 3.210942e-03f - 7.276687e-18f * _Complex_I,
906 -1.223392e-03f + 2.118978e-03f * _Complex_I, -9.721139e-04f - 1.683751e-03f * _Complex_I,
907 1.738862e-03f - 1.277695e-18f * _Complex_I, -6.914471e-04f + 1.197621e-03f * _Complex_I,
908 -5.616336e-04f - 9.727779e-04f * _Complex_I },
909 { -5.723872e-04f + 9.914038e-04f * _Complex_I, -8.302721e-04f - 1.438073e-03f * _Complex_I,
910 2.445280e-03f + 0.000000e+00f * _Complex_I, -1.378399e-03f + 2.387458e-03f * _Complex_I,
911 -1.882898e-03f - 3.261274e-03f * _Complex_I, 4.921549e-03f - 1.205432e-18f * _Complex_I,
912 -2.760152e-03f + 4.780723e-03f * _Complex_I, -2.460774e-03f - 4.262186e-03f * _Complex_I,
913 3.765795e-03f - 1.844708e-18f * _Complex_I, -1.378399e-03f + 2.387458e-03f * _Complex_I,
914 -1.222640e-03f - 2.117675e-03f * _Complex_I, 1.660544e-03f - 1.220148e-18f * _Complex_I,
915 -5.723872e-04f + 9.914038e-04f * _Complex_I },
916 { 1.226482e-03f + 3.004015e-19f * _Complex_I, -9.600816e-04f + 1.662910e-03f * _Complex_I,
917 -1.495900e-03f - 2.590974e-03f * _Complex_I, 3.833507e-03f + 0.000000e+00f * _Complex_I,
918 -3.167257e-03f + 5.485850e-03f * _Complex_I, -4.303595e-03f - 7.454046e-03f * _Complex_I,
919 9.412791e-03f - 2.305469e-18f * _Complex_I, -4.303595e-03f + 7.454046e-03f * _Complex_I,
920 -3.167257e-03f - 5.485850e-03f * _Complex_I, 3.833507e-03f - 1.877877e-18f * _Complex_I,
921 -1.495900e-03f + 2.590974e-03f * _Complex_I, -9.600816e-04f - 1.662910e-03f * _Complex_I,
922 1.226482e-03f - 9.012046e-19f * _Complex_I },
923 { -9.898007e-04f - 1.714385e-03f * _Complex_I, 3.215120e-03f + 3.643077e-18f * _Complex_I,
924 -2.507621e-03f + 4.343327e-03f * _Complex_I, -3.557798e-03f - 6.162286e-03f * _Complex_I,
925 1.105198e-02f + 0.000000e+00f * _Complex_I, -7.691179e-03f + 1.332151e-02f * _Complex_I,
926 -8.705793e-03f - 1.507888e-02f * _Complex_I, 1.538236e-02f - 3.767591e-18f * _Complex_I,
927 -5.525988e-03f + 9.571292e-03f * _Complex_I, -3.557798e-03f - 6.162286e-03f * _Complex_I,
928 5.015242e-03f - 2.456760e-18f * _Complex_I, -1.607560e-03f + 2.784375e-03f * _Complex_I,
929 -9.898007e-04f - 1.714385e-03f * _Complex_I },
930 { -1.414655e-03f + 2.450254e-03f * _Complex_I, -2.341263e-03f - 4.055186e-03f * _Complex_I,
931 6.915775e-03f + 1.693876e-18f * _Complex_I, -5.086403e-03f + 8.809908e-03f * _Complex_I,
932 -8.062191e-03f - 1.396412e-02f * _Complex_I, 2.415333e-02f + 0.000000e+00f * _Complex_I,
933 -1.451128e-02f + 2.513428e-02f * _Complex_I, -1.207667e-02f - 2.091740e-02f * _Complex_I,
934 1.612438e-02f - 3.949335e-18f * _Complex_I, -5.086403e-03f + 8.809908e-03f * _Complex_I,
935 -3.457887e-03f - 5.989237e-03f * _Complex_I, 4.682526e-03f - 1.061161e-17f * _Complex_I,
936 -1.414655e-03f + 2.450254e-03f * _Complex_I },
937 { 3.039574e-03f + 1.488962e-18f * _Complex_I, -2.598226e-03f + 4.500260e-03f * _Complex_I,
938 -3.750909e-03f - 6.496765e-03f * _Complex_I, 1.119776e-02f + 2.742661e-18f * _Complex_I,
939 -9.210579e-03f + 1.595319e-02f * _Complex_I, -1.464762e-02f - 2.537042e-02f * _Complex_I,
940 3.672076e-02f + 0.000000e+00f * _Complex_I, -1.464762e-02f + 2.537042e-02f * _Complex_I,
941 -9.210579e-03f - 1.595319e-02f * _Complex_I, 1.119776e-02f - 2.742661e-18f * _Complex_I,
942 -3.750909e-03f + 6.496765e-03f * _Complex_I, -2.598226e-03f - 4.500260e-03f * _Complex_I,
943 3.039574e-03f - 1.488962e-18f * _Complex_I },
944 { -1.414655e-03f - 2.450254e-03f * _Complex_I, 4.682526e-03f + 1.061161e-17f * _Complex_I,
945 -3.457887e-03f + 5.989237e-03f * _Complex_I, -5.086403e-03f - 8.809908e-03f * _Complex_I,
946 1.612438e-02f + 3.949335e-18f * _Complex_I, -1.207667e-02f + 2.091740e-02f * _Complex_I,
947 -1.451128e-02f - 2.513428e-02f * _Complex_I, 2.415333e-02f + 0.000000e+00f * _Complex_I,
948 -8.062191e-03f + 1.396412e-02f * _Complex_I, -5.086403e-03f - 8.809908e-03f * _Complex_I,
949 6.915775e-03f - 1.693876e-18f * _Complex_I, -2.341263e-03f + 4.055186e-03f * _Complex_I,
950 -1.414655e-03f - 2.450254e-03f * _Complex_I },
951 { -9.898007e-04f + 1.714385e-03f * _Complex_I, -1.607560e-03f - 2.784375e-03f * _Complex_I,
952 5.015242e-03f + 2.456760e-18f * _Complex_I, -3.557798e-03f + 6.162286e-03f * _Complex_I,
953 -5.525988e-03f - 9.571292e-03f * _Complex_I, 1.538236e-02f + 3.767591e-18f * _Complex_I,
954 -8.705793e-03f + 1.507888e-02f * _Complex_I, -7.691179e-03f - 1.332151e-02f * _Complex_I,
955 1.105198e-02f + 0.000000e+00f * _Complex_I, -3.557798e-03f + 6.162286e-03f * _Complex_I,
956 -2.507621e-03f - 4.343327e-03f * _Complex_I, 3.215120e-03f - 3.643077e-18f * _Complex_I,
957 -9.898007e-04f + 1.714385e-03f * _Complex_I },
958 { 1.226482e-03f + 9.012046e-19f * _Complex_I, -9.600816e-04f + 1.662910e-03f * _Complex_I,
959 -1.495900e-03f - 2.590974e-03f * _Complex_I, 3.833507e-03f + 1.877877e-18f * _Complex_I,
960 -3.167257e-03f + 5.485850e-03f * _Complex_I, -4.303595e-03f - 7.454046e-03f * _Complex_I,
961 9.412791e-03f + 2.305469e-18f * _Complex_I, -4.303595e-03f + 7.454046e-03f * _Complex_I,
962 -3.167257e-03f - 5.485850e-03f * _Complex_I, 3.833507e-03f + 0.000000e+00f * _Complex_I,
963 -1.495900e-03f + 2.590974e-03f * _Complex_I, -9.600816e-04f - 1.662910e-03f * _Complex_I,
964 1.226482e-03f - 3.004015e-19f * _Complex_I },
965 { -5.723872e-04f - 9.914038e-04f * _Complex_I, 1.660544e-03f + 1.220148e-18f * _Complex_I,
966 -1.222640e-03f + 2.117675e-03f * _Complex_I, -1.378399e-03f - 2.387458e-03f * _Complex_I,
967 3.765795e-03f + 1.844708e-18f * _Complex_I, -2.460774e-03f + 4.262186e-03f * _Complex_I,
968 -2.760152e-03f - 4.780723e-03f * _Complex_I, 4.921549e-03f + 1.205432e-18f * _Complex_I,
969 -1.882898e-03f + 3.261274e-03f * _Complex_I, -1.378399e-03f - 2.387458e-03f * _Complex_I,
970 2.445280e-03f + 0.000000e+00f * _Complex_I, -8.302721e-04f + 1.438073e-03f * _Complex_I,
971 -5.723872e-04f - 9.914038e-04f * _Complex_I },
972 { -5.616336e-04f + 9.727779e-04f * _Complex_I, -6.914471e-04f - 1.197621e-03f * _Complex_I,
973 1.738862e-03f + 1.277695e-18f * _Complex_I, -9.721139e-04f + 1.683751e-03f * _Complex_I,
974 -1.223392e-03f - 2.118978e-03f * _Complex_I, 3.210942e-03f + 7.276687e-18f * _Complex_I,
975 -1.832781e-03f + 3.174469e-03f * _Complex_I, -1.605471e-03f - 2.780758e-03f * _Complex_I,
976 2.446785e-03f + 2.772471e-18f * _Complex_I, -9.721139e-04f + 1.683751e-03f * _Complex_I,
977 -8.694311e-04f - 1.505899e-03f * _Complex_I, 1.382894e-03f + 0.000000e+00f * _Complex_I,
978 -5.616336e-04f + 9.727779e-04f * _Complex_I },
979 { 8.228091e-04f + 8.061204e-19f * _Complex_I, -5.365069e-04f + 9.292572e-04f * _Complex_I,
980 -6.011501e-04f - 1.041223e-03f * _Complex_I, 1.249890e-03f + 9.184039e-19f * _Complex_I,
981 -7.632708e-04f + 1.322024e-03f * _Complex_I, -9.846035e-04f - 1.705383e-03f * _Complex_I,
982 2.080486e-03f + 1.019144e-18f * _Complex_I, -9.846035e-04f + 1.705383e-03f * _Complex_I,
983 -7.632708e-04f - 1.322024e-03f * _Complex_I, 1.249890e-03f + 3.061346e-19f * _Complex_I,
984 -6.011501e-04f + 1.041223e-03f * _Complex_I, -5.365069e-04f - 9.292572e-04f * _Complex_I,
985 8.228091e-04f + 0.000000e+00f * _Complex_I } },
986 { { 1.221201e-03f + 5.982162e-19f * _Complex_I, -1.773498e-03f - 6.515727e-19f * _Complex_I,
987 1.246697e-03f + 3.053526e-19f * _Complex_I, -8.215306e-04f - 1.006085e-19f * _Complex_I,
988 7.609372e-04f + 0.000000e+00f * _Complex_I, -4.863927e-04f + 5.956592e-20f * _Complex_I,
989 4.882100e-04f - 1.195770e-19f * _Complex_I, -4.863927e-04f + 1.786978e-19f * _Complex_I,
990 7.609372e-04f - 3.727517e-19f * _Complex_I, -8.215306e-04f + 5.030424e-19f * _Complex_I,
991 1.246697e-03f - 9.160579e-19f * _Complex_I, -1.773498e-03f + 1.520336e-18f * _Complex_I,
992 1.221201e-03f - 1.196432e-18f * _Complex_I },
993 { 7.406884e-04f - 1.282910e-03f * _Complex_I, -1.025411e-03f + 1.776065e-03f * _Complex_I,
994 7.186273e-04f - 1.244699e-03f * _Complex_I, -4.025606e-04f + 6.972554e-04f * _Complex_I,
995 5.908383e-04f - 1.023362e-03f * _Complex_I, -1.125190e-03f + 1.948886e-03f * _Complex_I,
996 1.432695e-03f - 2.481501e-03f * _Complex_I, -1.125190e-03f + 1.948886e-03f * _Complex_I,
997 5.908383e-04f - 1.023362e-03f * _Complex_I, -4.025606e-04f + 6.972554e-04f * _Complex_I,
998 7.186273e-04f - 1.244699e-03f * _Complex_I, -1.025411e-03f + 1.776065e-03f * _Complex_I,
999 7.406884e-04f - 1.282910e-03f * _Complex_I },
1000 { -7.162255e-04f - 1.240539e-03f * _Complex_I, 8.961176e-04f + 1.552121e-03f * _Complex_I,
1001 -6.705589e-04f - 1.161442e-03f * _Complex_I, 6.187140e-04f + 1.071644e-03f * _Complex_I,
1002 -1.165433e-03f - 2.018589e-03f * _Complex_I, 1.948120e-03f + 3.374242e-03f * _Complex_I,
1003 -2.297663e-03f - 3.979669e-03f * _Complex_I, 1.948120e-03f + 3.374242e-03f * _Complex_I,
1004 -1.165433e-03f - 2.018589e-03f * _Complex_I, 6.187140e-04f + 1.071644e-03f * _Complex_I,
1005 -6.705589e-04f - 1.161442e-03f * _Complex_I, 8.961176e-04f + 1.552121e-03f * _Complex_I,
1006 -7.162255e-04f - 1.240539e-03f * _Complex_I },
1007 { -1.280260e-03f - 7.839331e-19f * _Complex_I, 1.987108e-03f + 9.734024e-19f * _Complex_I,
1008 -2.614019e-03f - 9.603749e-19f * _Complex_I, 3.635167e-03f + 8.903590e-19f * _Complex_I,
1009 -4.954867e-03f - 6.067962e-19f * _Complex_I, 6.653220e-03f + 0.000000e+00f * _Complex_I,
1010 -7.600546e-03f + 9.307984e-19f * _Complex_I, 6.653220e-03f - 1.629569e-18f * _Complex_I,
1011 -4.954867e-03f + 1.820389e-18f * _Complex_I, 3.635167e-03f - 1.780718e-18f * _Complex_I,
1012 -2.614019e-03f + 1.600625e-18f * _Complex_I, 1.987108e-03f - 1.460104e-18f * _Complex_I,
1013 -1.280260e-03f + 1.097506e-18f * _Complex_I },
1014 { -5.756945e-04f + 9.971322e-04f * _Complex_I, 1.268614e-03f - 2.197304e-03f * _Complex_I,
1015 -2.421407e-03f + 4.194000e-03f * _Complex_I, 4.045715e-03f - 7.007384e-03f * _Complex_I,
1016 -5.527367e-03f + 9.573681e-03f * _Complex_I, 6.837207e-03f - 1.184239e-02f * _Complex_I,
1017 -7.288212e-03f + 1.262355e-02f * _Complex_I, 6.837207e-03f - 1.184239e-02f * _Complex_I,
1018 -5.527367e-03f + 9.573681e-03f * _Complex_I, 4.045715e-03f - 7.007384e-03f * _Complex_I,
1019 -2.421407e-03f + 4.194000e-03f * _Complex_I, 1.268614e-03f - 2.197304e-03f * _Complex_I,
1020 -5.756945e-04f + 9.971322e-04f * _Complex_I },
1021 { 7.349896e-04f + 1.273039e-03f * _Complex_I, -1.748057e-03f - 3.027723e-03f * _Complex_I,
1022 3.332671e-03f + 5.772355e-03f * _Complex_I, -6.051736e-03f - 1.048191e-02f * _Complex_I,
1023 9.842376e-03f + 1.704749e-02f * _Complex_I, -1.401169e-02f - 2.426897e-02f * _Complex_I,
1024 1.598601e-02f + 2.768858e-02f * _Complex_I, -1.401169e-02f - 2.426897e-02f * _Complex_I,
1025 9.842376e-03f + 1.704749e-02f * _Complex_I, -6.051736e-03f - 1.048191e-02f * _Complex_I,
1026 3.332671e-03f + 5.772355e-03f * _Complex_I, -1.748057e-03f - 3.027723e-03f * _Complex_I,
1027 7.349896e-04f + 1.273039e-03f * _Complex_I },
1028 { 1.400383e-03f + 1.028985e-18f * _Complex_I, -3.545886e-03f - 2.171229e-18f * _Complex_I,
1029 7.289370e-03f + 3.570761e-18f * _Complex_I, -1.418908e-02f - 5.212982e-18f * _Complex_I,
1030 2.520839e-02f + 6.174275e-18f * _Complex_I, -3.934772e-02f - 4.818706e-18f * _Complex_I,
1031 4.797481e-02f + 0.000000e+00f * _Complex_I, -3.934772e-02f + 4.818706e-18f * _Complex_I,
1032 2.520839e-02f - 6.174275e-18f * _Complex_I, -1.418908e-02f + 5.212982e-18f * _Complex_I,
1033 7.289370e-03f - 3.570761e-18f * _Complex_I, -3.545886e-03f + 2.171229e-18f * _Complex_I,
1034 1.400383e-03f - 1.028985e-18f * _Complex_I },
1035 { 7.349896e-04f - 1.273039e-03f * _Complex_I, -1.748057e-03f + 3.027723e-03f * _Complex_I,
1036 3.332671e-03f - 5.772355e-03f * _Complex_I, -6.051736e-03f + 1.048191e-02f * _Complex_I,
1037 9.842376e-03f - 1.704749e-02f * _Complex_I, -1.401169e-02f + 2.426897e-02f * _Complex_I,
1038 1.598601e-02f - 2.768858e-02f * _Complex_I, -1.401169e-02f + 2.426897e-02f * _Complex_I,
1039 9.842376e-03f - 1.704749e-02f * _Complex_I, -6.051736e-03f + 1.048191e-02f * _Complex_I,
1040 3.332671e-03f - 5.772355e-03f * _Complex_I, -1.748057e-03f + 3.027723e-03f * _Complex_I,
1041 7.349896e-04f - 1.273039e-03f * _Complex_I },
1042 { -5.756945e-04f - 9.971322e-04f * _Complex_I, 1.268614e-03f + 2.197304e-03f * _Complex_I,
1043 -2.421407e-03f - 4.194000e-03f * _Complex_I, 4.045715e-03f + 7.007384e-03f * _Complex_I,
1044 -5.527367e-03f - 9.573681e-03f * _Complex_I, 6.837207e-03f + 1.184239e-02f * _Complex_I,
1045 -7.288212e-03f - 1.262355e-02f * _Complex_I, 6.837207e-03f + 1.184239e-02f * _Complex_I,
1046 -5.527367e-03f - 9.573681e-03f * _Complex_I, 4.045715e-03f + 7.007384e-03f * _Complex_I,
1047 -2.421407e-03f - 4.194000e-03f * _Complex_I, 1.268614e-03f + 2.197304e-03f * _Complex_I,
1048 -5.756945e-04f - 9.971322e-04f * _Complex_I },
1049 { -1.280260e-03f - 1.097506e-18f * _Complex_I, 1.987108e-03f + 1.460104e-18f * _Complex_I,
1050 -2.614019e-03f - 1.600625e-18f * _Complex_I, 3.635167e-03f + 1.780718e-18f * _Complex_I,
1051 -4.954867e-03f - 1.820389e-18f * _Complex_I, 6.653220e-03f + 1.629569e-18f * _Complex_I,
1052 -7.600546e-03f - 9.307984e-19f * _Complex_I, 6.653220e-03f + 0.000000e+00f * _Complex_I,
1053 -4.954867e-03f + 6.067962e-19f * _Complex_I, 3.635167e-03f - 8.903590e-19f * _Complex_I,
1054 -2.614019e-03f + 9.603749e-19f * _Complex_I, 1.987108e-03f - 9.734024e-19f * _Complex_I,
1055 -1.280260e-03f + 7.839331e-19f * _Complex_I },
1056 { -7.162255e-04f + 1.240539e-03f * _Complex_I, 8.961176e-04f - 1.552121e-03f * _Complex_I,
1057 -6.705589e-04f + 1.161442e-03f * _Complex_I, 6.187140e-04f - 1.071644e-03f * _Complex_I,
1058 -1.165433e-03f + 2.018589e-03f * _Complex_I, 1.948120e-03f - 3.374242e-03f * _Complex_I,
1059 -2.297663e-03f + 3.979669e-03f * _Complex_I, 1.948120e-03f - 3.374242e-03f * _Complex_I,
1060 -1.165433e-03f + 2.018589e-03f * _Complex_I, 6.187140e-04f - 1.071644e-03f * _Complex_I,
1061 -6.705589e-04f + 1.161442e-03f * _Complex_I, 8.961176e-04f - 1.552121e-03f * _Complex_I,
1062 -7.162255e-04f + 1.240539e-03f * _Complex_I },
1063 { 7.406884e-04f + 1.282910e-03f * _Complex_I, -1.025411e-03f - 1.776065e-03f * _Complex_I,
1064 7.186273e-04f + 1.244699e-03f * _Complex_I, -4.025606e-04f - 6.972554e-04f * _Complex_I,
1065 5.908383e-04f + 1.023362e-03f * _Complex_I, -1.125190e-03f - 1.948886e-03f * _Complex_I,
1066 1.432695e-03f + 2.481501e-03f * _Complex_I, -1.125190e-03f - 1.948886e-03f * _Complex_I,
1067 5.908383e-04f + 1.023362e-03f * _Complex_I, -4.025606e-04f - 6.972554e-04f * _Complex_I,
1068 7.186273e-04f + 1.244699e-03f * _Complex_I, -1.025411e-03f - 1.776065e-03f * _Complex_I,
1069 7.406884e-04f + 1.282910e-03f * _Complex_I },
1070 { 1.221201e-03f + 1.196432e-18f * _Complex_I, -1.773498e-03f - 1.520336e-18f * _Complex_I,
1071 1.246697e-03f + 9.160579e-19f * _Complex_I, -8.215306e-04f - 5.030424e-19f * _Complex_I,
1072 7.609372e-04f + 3.727517e-19f * _Complex_I, -4.863927e-04f - 1.786978e-19f * _Complex_I,
1073 4.882100e-04f + 1.195770e-19f * _Complex_I, -4.863927e-04f - 5.956592e-20f * _Complex_I,
1074 7.609372e-04f + 0.000000e+00f * _Complex_I, -8.215306e-04f + 1.006085e-19f * _Complex_I,
1075 1.246697e-03f - 3.053526e-19f * _Complex_I, -1.773498e-03f + 6.515727e-19f * _Complex_I,
1076 1.221201e-03f - 5.982162e-19f * _Complex_I } } };
1077
1078 const size_t buffer_size = (size_t)TS * TS * (ndir * 4 + 7) * sizeof(float);
1079 size_t padded_buffer_size;
1080 char *const all_buffers = (char *)dt_pixelpipe_cache_alloc_perthread(buffer_size, sizeof(char), &padded_buffer_size);
1081 if(IS_NULL_PTR(all_buffers))
1082 {
1083 fprintf(stderr, "[demosaic] not able to allocate FDC base buffers\n");
1084 return;
1085 }
1086
1087 /* Map a green hexagon around each non-green pixel and vice versa. allhex[] is later
1088 * indexed by hexmap() using tile-local (roi_in-relative) row/col taken mod 3, so it
1089 * must be built against that same absolute phase: pass roi_in here (not NULL) so the
1090 * color lookups land on the actual sensor position, not xtrans[][]'s own local origin. */
1091 for(int row = 0; row < 3; row++)
1092 for(int col = 0; col < 3; col++)
1093 for(int ng = 0, d = 0; d < 10; d += 2)
1094 {
1095 int g = FCxtrans(row, col, roi_in, xtrans) == 1;
1096 if(FCxtrans(row + orth[d], col + orth[d + 2], roi_in, xtrans) == 1)
1097 ng = 0;
1098 else
1099 ng++;
1100 // if there are four non-green pixels adjacent in cardinal
1101 // directions, this is the solitary green pixel
1102 if(ng == 4)
1103 {
1104 sgrow = row;
1105 sgcol = col;
1106 }
1107 if(ng == g + 1)
1108 for(int c = 0; c < 8; c++)
1109 {
1110 int v = orth[d] * patt[g][c * 2] + orth[d + 1] * patt[g][c * 2 + 1];
1111 int h = orth[d + 2] * patt[g][c * 2] + orth[d + 3] * patt[g][c * 2 + 1];
1112 // offset within TSxTS buffer
1113 allhex[row][col][c ^ (g * 2 & d)] = h + v * TS;
1114 }
1115 }
1116
1117 // extra passes propagates out errors at edges, hence need more padding
1118 const int pad_tile = 13;
1119
1120 // calculate offsets for this roi
1121 int rowoffset = 0;
1122 int coloffset = 0;
1123 for(int row = 0; row < 6; row++)
1124 {
1125 if(!((row - sgrow) % 3))
1126 {
1127 for(int col = 0; col < 6; col++)
1128 {
1129 if(!((col - sgcol) % 3) && (FCxtrans(row, col + 1, roi_in, xtrans) == 0))
1130 {
1131 rowoffset = 37 - row - pad_tile; // 1 plus a generous multiple of 6
1132 coloffset = 37 - col - pad_tile; // to avoid that this value gets negative
1133 break;
1134 }
1135 }
1136 break;
1137 }
1138 }
1139
1140 // depending on the iso, use either a hybrid approach for chroma, or pure fdc
1141 float hybrid_fdc[2] = { 1.0f, 0.0f };
1142 const int xover_iso = dt_conf_get_int("plugins/darkroom/demosaic/fdc_xover_iso");
1143 int iso = self->dev->image_storage.exif_iso;
1144 if(iso > xover_iso)
1145 {
1146 hybrid_fdc[0] = 0.0f;
1147 hybrid_fdc[1] = 1.0f;
1148 }
1149 __OMP_PARALLEL_FOR__()
1150 // step through TSxTS cells of image, each tile overlapping the
1151 // prior as interpolation needs a substantial border
1152 for(int top = -pad_tile; top < height - pad_tile; top += TS - (pad_tile * 2))
1153 {
1154 char *const buffer = dt_get_perthread(all_buffers, padded_buffer_size);
1155 // rgb points to ndir TSxTS tiles of 3 channels (R, G, and B)
1156 float(*rgb)[TS][TS][3] = (float(*)[TS][TS][3])buffer;
1157 // yuv points to 3 channel (Y, u, and v) TSxTS tiles
1158 // note that channels come before tiles to allow for a
1159 // vectorization optimization when building drv[] from yuv[]
1160 float (*const yuv)[TS][TS] = (float(*)[TS][TS])(buffer + TS * TS * (ndir * 3) * sizeof(float));
1161 // drv points to ndir TSxTS tiles, each a single channel of derivatives
1162 float (*const drv)[TS][TS] = (float(*)[TS][TS])(buffer + TS * TS * (ndir * 3 + 3) * sizeof(float));
1163 // gmin and gmax reuse memory which is used later by yuv buffer;
1164 // each points to a TSxTS tile of single channel data
1165 float (*const gmin)[TS] = (float(*)[TS])(buffer + TS * TS * (ndir * 3) * sizeof(float));
1166 float (*const gmax)[TS] = (float(*)[TS])(buffer + TS * TS * (ndir * 3 + 1) * sizeof(float));
1167 // homo and homosum reuse memory which is used earlier in the
1168 // loop; each points to ndir single-channel TSxTS tiles
1169 uint8_t (*const homo)[TS][TS] = (uint8_t(*)[TS][TS])(buffer + TS * TS * (ndir * 3) * sizeof(float));
1170 uint8_t (*const homosum)[TS][TS]
1171 = (uint8_t(*)[TS][TS])(buffer + TS * TS * (ndir * 3) * sizeof(float) + TS * TS * ndir * sizeof(uint8_t));
1172 // append all fdc related buffers
1173 float complex *fdc_buf_start = (float complex *)(buffer + TS * TS * (ndir * 4 + 3) * sizeof(float));
1174 const int fdc_buf_size = TS * TS;
1175 float(*const i_src) = (float *)fdc_buf_start;
1176 float complex(*const o_src) = fdc_buf_start + fdc_buf_size;
1177 // by the time the chroma values are calculated, o_src can be overwritten.
1178 float(*const fdc_chroma) = (float *)o_src;
1179
1180 for(int left = -pad_tile; left < width - pad_tile; left += TS - (pad_tile * 2))
1181 {
1182 int mrow = MIN(top + TS, height + pad_tile);
1183 int mcol = MIN(left + TS, width + pad_tile);
1184
1185 // Copy current tile from in to image buffer. If border goes
1186 // beyond edges of image, fill with mirrored/interpolated edges.
1187 // The extra border avoids discontinuities at image edges.
1188 for(int row = top; row < mrow; row++)
1189 for(int col = left; col < mcol; col++)
1190 {
1191 float(*const pix) = rgb[0][row - top][col - left];
1192 if((col >= 0) && (row >= 0) && (col < width) && (row < height))
1193 {
1194 const int f = FCxtrans(row, col, roi_in, xtrans);
1195 for(int c = 0; c < 3; c++) pix[c] = (c == f) ? in[roi_in->width * row + col] : 0.f;
1196 *(i_src + TS * (row - top) + (col - left)) = in[roi_in->width * row + col];
1197 }
1198 else
1199 {
1200 // mirror a border pixel if beyond image edge
1201 const int c = FCxtrans(row, col, roi_in, xtrans);
1202 for(int cc = 0; cc < 3; cc++)
1203 if(cc != c)
1204 pix[cc] = 0.0f;
1205 else
1206 {
1207#define TRANSLATE(n, size) ((n >= size) ? (2 * size - n - 2) : abs(n))
1208 const int cy = TRANSLATE(row, height), cx = TRANSLATE(col, width);
1209 if(c == FCxtrans(cy, cx, roi_in, xtrans))
1210 {
1211 pix[c] = in[roi_in->width * cy + cx];
1212 *(i_src + TS * (row - top) + (col - left)) = in[roi_in->width * cy + cx];
1213 }
1214 else
1215 {
1216 // interpolate if mirror pixel is a different color
1217 float sum = 0.0f;
1218 uint8_t count = 0;
1219 for(int y = row - 1; y <= row + 1; y++)
1220 for(int x = col - 1; x <= col + 1; x++)
1221 {
1222 const int yy = TRANSLATE(y, height), xx = TRANSLATE(x, width);
1223 const int ff = FCxtrans(yy, xx, roi_in, xtrans);
1224 if(ff == c)
1225 {
1226 sum += in[roi_in->width * yy + xx];
1227 count++;
1228 }
1229 }
1230 pix[c] = sum / count;
1231 *(i_src + TS * (row - top) + (col - left)) = pix[c];
1232 }
1233 }
1234 }
1235 }
1236
1237 // duplicate rgb[0] to rgb[1], rgb[2], and rgb[3]
1238 for(int c = 1; c <= 3; c++) memcpy(rgb[c], rgb[0], sizeof(*rgb));
1239
1240 // note that successive calculations are inset within the tile
1241 // so as to give enough border data, and there needs to be a 6
1242 // pixel border initially to allow allhex to find neighboring
1243 // pixels
1244
1245 /* Set green1 and green3 to the minimum and maximum allowed values: */
1246 // Run through each red/blue or blue/red pair, setting their g1
1247 // and g3 values to the min/max of green pixels surrounding the
1248 // pair. Use a 3 pixel border as gmin/gmax is used by
1249 // interpolate green which has a 3 pixel border.
1250 const int pad_g1_g3 = 3;
1251 for(int row = top + pad_g1_g3; row < mrow - pad_g1_g3; row++)
1252 {
1253 // setting max to 0.0f signifies that this is a new pair, which
1254 // requires a new min/max calculation of its neighboring greens
1255 float min = FLT_MAX, max = 0.0f;
1256 for(int col = left + pad_g1_g3; col < mcol - pad_g1_g3; col++)
1257 {
1258 // if in row of horizontal red & blue pairs (or processing
1259 // vertical red & blue pairs near image bottom), reset min/max
1260 // between each pair
1261 if(FCxtrans(row, col, roi_in, xtrans) == 1)
1262 {
1263 min = FLT_MAX, max = 0.0f;
1264 continue;
1265 }
1266 // if at start of red & blue pair, calculate min/max of green
1267 // pixels surrounding it; note that while normally using == to
1268 // compare floats is suspect, here the check is if 0.0f has
1269 // explicitly been assigned to max (which signifies a new
1270 // red/blue pair)
1271 if(max == 0.0f)
1272 {
1273 float (*const pix)[3] = &rgb[0][row - top][col - left];
1274 const short *const hex = hexmap(row, col, allhex);
1275 for(int c = 0; c < 6; c++)
1276 {
1277 const float val = pix[hex[c]][1];
1278 if(min > val) min = val;
1279 if(max < val) max = val;
1280 }
1281 }
1282 gmin[row - top][col - left] = min;
1283 gmax[row - top][col - left] = max;
1284 // handle vertical red/blue pairs
1285 switch((row - sgrow) % 3)
1286 {
1287 // hop down a row to second pixel in vertical pair
1288 case 1:
1289 if(row < mrow - 4) row++, col--;
1290 break;
1291 // then if not done with the row hop up and right to next
1292 // vertical red/blue pair, resetting min/max
1293 case 2:
1294 min = FLT_MAX, max = 0.0f;
1295 if((col += 2) < mcol - 4 && row > top + 3) row--;
1296 }
1297 }
1298 }
1299
1300 /* Interpolate green horizontally, vertically, and along both diagonals: */
1301 // need a 3 pixel border here as 3*hex[] can have a 3 unit offset
1302 const int pad_g_interp = 3;
1303 for(int row = top + pad_g_interp; row < mrow - pad_g_interp; row++)
1304 for(int col = left + pad_g_interp; col < mcol - pad_g_interp; col++)
1305 {
1306 float color[8];
1307 int f = FCxtrans(row, col, roi_in, xtrans);
1308 if(f == 1) continue;
1309 float (*const pix)[3] = &rgb[0][row - top][col - left];
1310 const short *const hex = hexmap(row, col, allhex);
1311 // TODO: these constants come from integer math constants in
1312 // dcraw -- calculate them instead from interpolation math
1313 color[0] = 0.6796875f * (pix[hex[1]][1] + pix[hex[0]][1])
1314 - 0.1796875f * (pix[2 * hex[1]][1] + pix[2 * hex[0]][1]);
1315 color[1] = 0.87109375f * pix[hex[3]][1] + pix[hex[2]][1] * 0.13f
1316 + 0.359375f * (pix[0][f] - pix[-hex[2]][f]);
1317 for(int c = 0; c < 2; c++)
1318 color[2 + c] = 0.640625f * pix[hex[4 + c]][1] + 0.359375f * pix[-2 * hex[4 + c]][1]
1319 + 0.12890625f * (2 * pix[0][f] - pix[3 * hex[4 + c]][f] - pix[-3 * hex[4 + c]][f]);
1320 for(int c = 0; c < 4; c++)
1321 rgb[c ^ !((row - sgrow) % 3)][row - top][col - left][1]
1322 = CLAMPS(color[c], gmin[row - top][col - left], gmax[row - top][col - left]);
1323 }
1324
1325 /* Interpolate red and blue values for solitary green pixels: */
1326 const int pad_rb_g = 6;
1327 for(int row = (top - sgrow + pad_rb_g + 2) / 3 * 3 + sgrow; row < mrow - pad_rb_g; row += 3)
1328 for(int col = (left - sgcol + pad_rb_g + 2) / 3 * 3 + sgcol; col < mcol - pad_rb_g; col += 3)
1329 {
1330 float(*rfx)[3] = &rgb[0][row - top][col - left];
1331 int h = FCxtrans(row, col + 1, roi_in, xtrans);
1332 float diff[6] = { 0.0f };
1333 float color[3][8];
1334 for(int i = 1, d = 0; d < 6; d++, i ^= TS ^ 1, h ^= 2)
1335 {
1336 for(int c = 0; c < 2; c++, h ^= 2)
1337 {
1338 float g = 2 * rfx[0][1] - rfx[i << c][1] - rfx[-(i << c)][1];
1339 color[h][d] = g + rfx[i << c][h] + rfx[-(i << c)][h];
1340 if(d > 1)
1341 diff[d] += SQR(rfx[i << c][1] - rfx[-(i << c)][1] - rfx[i << c][h] + rfx[-(i << c)][h]) + SQR(g);
1342 }
1343 if(d > 1 && (d & 1))
1344 if(diff[d - 1] < diff[d])
1345 for(int c = 0; c < 2; c++) color[c * 2][d] = color[c * 2][d - 1];
1346 if(d < 2 || (d & 1))
1347 {
1348 for(int c = 0; c < 2; c++) rfx[0][c * 2] = color[c * 2][d] / 2.f;
1349 rfx += TS * TS;
1350 }
1351 }
1352 }
1353
1354 /* Interpolate red for blue pixels and vice versa: */
1355 const int pad_rb_br = 6;
1356 for(int row = top + pad_rb_br; row < mrow - pad_rb_br; row++)
1357 for(int col = left + pad_rb_br; col < mcol - pad_rb_br; col++)
1358 {
1359 int f = 2 - FCxtrans(row, col, roi_in, xtrans);
1360 if(f == 1) continue;
1361 float(*rfx)[3] = &rgb[0][row - top][col - left];
1362 int c = (row - sgrow) % 3 ? TS : 1;
1363 int h = 3 * (c ^ TS ^ 1);
1364 for(int d = 0; d < 4; d++, rfx += TS * TS)
1365 {
1366 int i = d > 1 || ((d ^ c) & 1)
1367 || ((fabsf(rfx[0][1] - rfx[c][1]) + fabsf(rfx[0][1] - rfx[-c][1]))
1368 < 2.f * (fabsf(rfx[0][1] - rfx[h][1]) + fabsf(rfx[0][1] - rfx[-h][1]))) ? c : h;
1369 rfx[0][f] = (rfx[i][f] + rfx[-i][f] + 2.f * rfx[0][1] - rfx[i][1] - rfx[-i][1]) / 2.f;
1370 }
1371 }
1372
1373 /* Fill in red and blue for 2x2 blocks of green: */
1374 const int pad_g22 = 8;
1375 for(int row = top + pad_g22; row < mrow - pad_g22; row++)
1376 if((row - sgrow) % 3)
1377 for(int col = left + pad_g22; col < mcol - pad_g22; col++)
1378 if((col - sgcol) % 3)
1379 {
1380 float redblue[3][3];
1381 float(*rfx)[3] = &rgb[0][row - top][col - left];
1382 const short *const hex = hexmap(row, col, allhex);
1383 for(int d = 0; d < ndir; d += 2, rfx += TS * TS)
1384 if(hex[d] + hex[d + 1])
1385 {
1386 float g = 3.f * rfx[0][1] - 2.f * rfx[hex[d]][1] - rfx[hex[d + 1]][1];
1387 for(int c = 0; c < 4; c += 2)
1388 {
1389 rfx[0][c] = (g + 2.f * rfx[hex[d]][c] + rfx[hex[d + 1]][c]) / 3.f;
1390 redblue[d][c] = rfx[0][c];
1391 }
1392 }
1393 else
1394 {
1395 float g = 2.f * rfx[0][1] - rfx[hex[d]][1] - rfx[hex[d + 1]][1];
1396 for(int c = 0; c < 4; c += 2)
1397 {
1398 rfx[0][c] = (g + rfx[hex[d]][c] + rfx[hex[d + 1]][c]) / 2.f;
1399 redblue[d][c] = rfx[0][c];
1400 }
1401 }
1402 // to fill in red and blue also for diagonal directions
1403 for(int d = 0; d < ndir; d += 2, rfx += TS * TS)
1404 for(int c = 0; c < 4; c += 2) rfx[0][c] = (redblue[0][c] + redblue[2][c]) * 0.5f;
1405 }
1406
1407 // jump back to the first set of rgb buffers (this is a nop
1408 // unless on the second pass)
1409 rgb = (float(*)[TS][TS][3])buffer;
1410 // from here on out, mainly are working within the current tile
1411 // rather than in reference to the image, so don't offset
1412 // mrow/mcol by top/left of tile
1413 mrow -= top;
1414 mcol -= left;
1415
1416 /* Convert to perceptual colorspace and differentiate in all directions: */
1417 // Original dcraw algorithm uses CIELab as perceptual space
1418 // (presumably coming from original AHD) and converts taking
1419 // camera matrix into account. Now use YPbPr which requires much
1420 // less code and is nearly indistinguishable. It assumes the
1421 // camera RGB is roughly linear.
1422 for(int d = 0; d < ndir; d++)
1423 {
1424 const int pad_yuv = 8;
1425 for(int row = pad_yuv; row < mrow - pad_yuv; row++)
1426 for(int col = pad_yuv; col < mcol - pad_yuv; col++)
1427 {
1428 float *rx = rgb[d][row][col];
1429 // use ITU-R BT.2020 YPbPr, which is great, but could use
1430 // a better/simpler choice? note that imageop.h provides
1431 // dt_iop_RGB_to_YCbCr which uses Rec. 601 conversion,
1432 // which appears less good with specular highlights
1433 float y = 0.2627f * rx[0] + 0.6780f * rx[1] + 0.0593f * rx[2];
1434 yuv[0][row][col] = y;
1435 yuv[1][row][col] = (rx[2] - y) * 0.56433f;
1436 yuv[2][row][col] = (rx[0] - y) * 0.67815f;
1437 }
1438 // Note that f can offset by a column (-1 or +1) and by a row
1439 // (-TS or TS). The row-wise offsets cause the undefined
1440 // behavior sanitizer to warn of an out of bounds index, but
1441 // as yfx is multi-dimensional and there is sufficient
1442 // padding, that is not actually so.
1443 const int f = dir[d & 3];
1444 const int pad_drv = 9;
1445 for(int row = pad_drv; row < mrow - pad_drv; row++)
1446 for(int col = pad_drv; col < mcol - pad_drv; col++)
1447 {
1448 float(*yfx)[TS][TS] = (float(*)[TS][TS]) & yuv[0][row][col];
1449 drv[d][row][col] = SQR(2 * yfx[0][0][0] - yfx[0][0][f] - yfx[0][0][-f])
1450 + SQR(2 * yfx[1][0][0] - yfx[1][0][f] - yfx[1][0][-f])
1451 + SQR(2 * yfx[2][0][0] - yfx[2][0][f] - yfx[2][0][-f]);
1452 }
1453 }
1454
1455 /* Build homogeneity maps from the derivatives: */
1456 memset(homo, 0, sizeof(uint8_t) * ndir * TS * TS);
1457 const int pad_homo = 10;
1458 for(int row = pad_homo; row < mrow - pad_homo; row++)
1459 for(int col = pad_homo; col < mcol - pad_homo; col++)
1460 {
1461 float tr = FLT_MAX;
1462 for(int d = 0; d < ndir; d++)
1463 if(tr > drv[d][row][col]) tr = drv[d][row][col];
1464 tr *= 8;
1465 for(int d = 0; d < ndir; d++)
1466 for(int v = -1; v <= 1; v++)
1467 for(int h = -1; h <= 1; h++) homo[d][row][col] += ((drv[d][row + v][col + h] <= tr) ? 1 : 0);
1468 }
1469
1470 /* Build 5x5 sum of homogeneity maps for each pixel & direction */
1471 for(int d = 0; d < ndir; d++)
1472 for(int row = pad_tile; row < mrow - pad_tile; row++)
1473 {
1474 // start before first column where homo[d][row][col+2] != 0,
1475 // so can know v5sum and homosum[d][row][col] will be 0
1476 int col = pad_tile - 5;
1477 uint8_t v5sum[5] = { 0 };
1478 homosum[d][row][col] = 0;
1479 // calculate by rolling through column sums
1480 for(col++; col < mcol - pad_tile; col++)
1481 {
1482 uint8_t colsum = 0;
1483 for(int v = -2; v <= 2; v++) colsum += homo[d][row + v][col + 2];
1484 homosum[d][row][col] = homosum[d][row][col - 1] - v5sum[col % 5] + colsum;
1485 v5sum[col % 5] = colsum;
1486 }
1487 }
1488
1489 /* Calculate chroma values in fdc: */
1490 const int pad_fdc = 6;
1491 for(int row = pad_fdc; row < mrow - pad_fdc; row++)
1492 for(int col = pad_fdc; col < mcol - pad_fdc; col++)
1493 {
1494 int myrow, mycol;
1495 uint8_t hm[8] = { 0 };
1496 uint8_t maxval = 0;
1497 for(int d = 0; d < ndir; d++)
1498 {
1499 hm[d] = homosum[d][row][col];
1500 maxval = (maxval < hm[d] ? hm[d] : maxval);
1501 }
1502 maxval -= maxval >> 3;
1503 float dircount = 0;
1504 float dirsum = 0.f;
1505 for(int d = 0; d < ndir; d++)
1506 if(hm[d] >= maxval)
1507 {
1508 dircount++;
1509 dirsum += directionality[d];
1510 }
1511 float w = dirsum / (float)dircount;
1512 int fdc_row, fdc_col;
1513 float complex C2m, C5m, C7m, C10m;
1514#define CONV_FILT(VAR, FILT) \
1515 VAR = 0.0f + 0.0f * _Complex_I; \
1516 for(fdc_row = 0, myrow = row - 6; fdc_row < 13; fdc_row++, myrow++) \
1517 for(fdc_col = 0, mycol = col - 6; fdc_col < 13; fdc_col++, mycol++) \
1518 VAR += FILT[12 - fdc_row][12 - fdc_col] * *(i_src + TS * myrow + mycol);
1519 CONV_FILT(C2m, harr[0])
1520 CONV_FILT(C5m, harr[1])
1521 CONV_FILT(C7m, harr[2])
1522 CONV_FILT(C10m, harr[3])
1523#undef CONV_FILT
1524 // build the q vector components
1525 myrow = (row + rowoffset) % 6;
1526 mycol = (col + coloffset) % 6;
1527 float complex modulator[8];
1528 for(int c = 0; c < 8; c++) modulator[c] = modarr[myrow][mycol][c];
1529 float complex qmat[8];
1530 qmat[4] = w * C10m * modulator[0] - (1.0f - w) * C2m * modulator[1];
1531 qmat[6] = conjf(qmat[4]);
1532 qmat[1] = C5m * modulator[6];
1533 qmat[2] = conjf(-0.5f * qmat[1]);
1534 qmat[5] = conjf(qmat[2]);
1535 qmat[3] = C7m * modulator[7];
1536 qmat[7] = conjf(qmat[1]);
1537 // get L
1538 C2m = qmat[4] * (conjf(modulator[0]) - conjf(modulator[1]));
1539 float complex C3m = qmat[6] * (modulator[2] - modulator[3]);
1540 float complex C6m = qmat[2] * (conjf(modulator[4]) + conjf(modulator[5]));
1541 float complex C12m = qmat[5] * (modulator[4] + modulator[5]);
1542 float complex C18m = qmat[7] * modulator[6];
1543 qmat[0] = *(i_src + row * TS + col) - C2m - C3m - C5m - C6m - 2.0f * C7m - C12m - C18m;
1544 // get the rgb components from fdc
1545 dt_aligned_pixel_t rgbpix = { 0.f, 0.f, 0.f };
1546 // multiply with the inverse matrix of M
1547 for(int color = 0; color < 3; color++)
1548 for(int c = 0; c < 8; c++)
1549 {
1550 rgbpix[color] += Minv[color][c] * qmat[c];
1551 }
1552 // now separate luma and chroma for
1553 // frequency domain chroma
1554 // and store it in fdc_chroma
1555 float uv[2];
1556 float y = 0.2627f * rgbpix[0] + 0.6780f * rgbpix[1] + 0.0593f * rgbpix[2];
1557 uv[0] = (rgbpix[2] - y) * 0.56433f;
1558 uv[1] = (rgbpix[0] - y) * 0.67815f;
1559 for(int c = 0; c < 2; c++) *(fdc_chroma + c * TS * TS + row * TS + col) = uv[c];
1560 }
1561
1562 /* Average the most homogeneous pixels for the final result: */
1563 for(int row = pad_tile; row < mrow - pad_tile; row++)
1564 for(int col = pad_tile; col < mcol - pad_tile; col++)
1565 {
1566 uint8_t hm[8] = { 0 };
1567 uint8_t maxval = 0;
1568 for(int d = 0; d < ndir; d++)
1569 {
1570 hm[d] = homosum[d][row][col];
1571 maxval = (maxval < hm[d] ? hm[d] : maxval);
1572 }
1573 maxval -= maxval >> 3;
1574 for(int d = 0; d < ndir - 4; d++)
1575 {
1576 if(hm[d] < hm[d + 4])
1577 hm[d] = 0;
1578 else if(hm[d] > hm[d + 4])
1579 hm[d + 4] = 0;
1580 }
1581
1582 dt_aligned_pixel_t avg = { 0.f };
1583 for(int d = 0; d < ndir; d++)
1584 {
1585 if(hm[d] >= maxval)
1586 {
1587 for(int c = 0; c < 3; c++) avg[c] += rgb[d][row][col][c];
1588 avg[3]++;
1589 }
1590 }
1591 dt_aligned_pixel_t rgbpix;
1592 for(int c = 0; c < 3; c++) rgbpix[c] = avg[c] / avg[3];
1593 // preserve all components of Markesteijn for this pixel
1594 const float y = 0.2627f * rgbpix[0] + 0.6780f * rgbpix[1] + 0.0593f * rgbpix[2];
1595 const float um = (rgbpix[2] - y) * 0.56433f;
1596 const float vm = (rgbpix[0] - y) * 0.67815f;
1597 float uvf[2];
1598 // macros for fast meadian filtering
1599#define PIX_SWAP(a, b) \
1600 { \
1601 tempf = (a); \
1602 (a) = (b); \
1603 (b) = tempf; \
1604 }
1605#define PIX_SORT(a, b) \
1606 { \
1607 if((a) > (b)) PIX_SWAP((a), (b)); \
1608 }
1609 // instead of merely reading the values, perform 5 pixel median filter
1610 // one median filter is required to avoid textile artifacts
1611 for(int chrm = 0; chrm < 2; chrm++)
1612 {
1613 float temp[5];
1614 float tempf;
1615 // load the window into temp
1616 memcpy(&temp[0], fdc_chroma + chrm * TS * TS + (row - 1) * TS + (col), sizeof(float) * 1);
1617 memcpy(&temp[1], fdc_chroma + chrm * TS * TS + (row)*TS + (col - 1), sizeof(float) * 3);
1618 memcpy(&temp[4], fdc_chroma + chrm * TS * TS + (row + 1) * TS + (col), sizeof(float) * 1);
1619 PIX_SORT(temp[0], temp[1]);
1620 PIX_SORT(temp[3], temp[4]);
1621 PIX_SORT(temp[0], temp[3]);
1622 PIX_SORT(temp[1], temp[4]);
1623 PIX_SORT(temp[1], temp[2]);
1624 PIX_SORT(temp[2], temp[3]);
1625 PIX_SORT(temp[1], temp[2]);
1626 uvf[chrm] = temp[2];
1627 }
1628 // use hybrid or pure fdc, depending on what was set above.
1629 // in case of hybrid, use the chroma that has the smallest
1630 // absolute value
1631 float uv[2];
1632 uv[0] = (((ABS(uvf[0]) < ABS(um)) & (ABS(uvf[1]) < (1.02f * ABS(vm)))) ? uvf[0] : um) * hybrid_fdc[0] + uvf[0] * hybrid_fdc[1];
1633 uv[1] = (((ABS(uvf[1]) < ABS(vm)) & (ABS(uvf[0]) < (1.02f * ABS(vm)))) ? uvf[1] : vm) * hybrid_fdc[0] + uvf[1] * hybrid_fdc[1];
1634 // combine the luma from Markesteijn with the chroma from above
1635 rgbpix[0] = y + 1.474600014746f * uv[1];
1636 rgbpix[1] = y - 0.15498578286403f * uv[0] - 0.571353132557189f * uv[1];
1637 rgbpix[2] = y + 1.77201282937288f * uv[0];
1638 for(int c = 0; c < 3; c++) out[4 * (width * (row + top) + col + left) + c] = rgbpix[c];
1639 }
1640 }
1641 }
1642
1643 dt_pixelpipe_cache_free_align(all_buffers);
1644}
1645
1646#ifdef HAVE_OPENCL
1647
1648static int process_markesteijn_cl(struct dt_iop_module_t *self, const dt_dev_pixelpipe_t *pipe,
1649 const dt_dev_pixelpipe_iop_t *piece, cl_mem dev_in, cl_mem dev_out,
1650 const dt_iop_roi_t *const roi_in, const dt_iop_roi_t *const roi_out,
1651 const gboolean smooth)
1652{
1653 dt_iop_demosaic_data_t *data = (dt_iop_demosaic_data_t *)piece->data;
1654 dt_iop_demosaic_global_data_t *gd = (dt_iop_demosaic_global_data_t *)self->global_data;
1655
1656 const int devid = pipe->devid;
1657 const uint8_t(*const xtrans)[6] = (const uint8_t(*const)[6])piece->dsc_in.xtrans;
1658
1659 const float processed_maximum[4]
1660 = { piece->dsc_in.processed_maximum[0], piece->dsc_in.processed_maximum[1],
1661 piece->dsc_in.processed_maximum[2], 1.0f };
1662
1663 cl_mem dev_tmp = NULL;
1664 cl_mem dev_tmptmp = NULL;
1665 cl_mem dev_xtrans = NULL;
1666 cl_mem dev_green_eq = NULL;
1667 cl_mem dev_rgbv[8] = { NULL };
1668 cl_mem dev_drv[8] = { NULL };
1669 cl_mem dev_homo[8] = { NULL };
1670 cl_mem dev_homosum[8] = { NULL };
1671 cl_mem dev_gminmax = NULL;
1672 cl_mem dev_allhex = NULL;
1673 cl_mem dev_aux = NULL;
1674 cl_mem dev_edge_in = NULL;
1675 cl_mem dev_edge_out = NULL;
1676 cl_int err = -999;
1677
1678 cl_mem *dev_rgb = dev_rgbv;
1679
1680 dev_xtrans = dt_opencl_copy_host_to_device_constant(devid, sizeof(piece->dsc_in.xtrans), (void *)piece->dsc_in.xtrans);
1681 if(IS_NULL_PTR(dev_xtrans)) goto error;
1682
1683 int width = roi_in->width;
1684 int height = roi_in->height;
1685 const int passes = ((data->demosaicing_method & ~DEMOSAIC_DUAL) == DT_IOP_DEMOSAIC_MARKESTEIJN_3) ? 3 : 1;
1686 const int ndir = (passes > 1) ? 8 : 4; // only ever 4 or 8
1687 const int pad_tile = (passes == 1) ? 12 : 17;
1688
1689 static const short orth[12] = { 1, 0, 0, 1, -1, 0, 0, -1, 1, 0, 0, 1 },
1690 patt[2][16] = { { 0, 1, 0, -1, 2, 0, -1, 0, 1, 1, 1, -1, 0, 0, 0, 0 },
1691 { 0, 1, 0, -2, 1, 0, -2, 0, 1, 1, -2, -2, 1, -1, -1, 1 } };
1692
1693 // allhex contains the offset coordinates (x,y) of a green hexagon around each
1694 // non-green pixel and vice versa
1695 char allhex[3][3][8][2];
1696 // sgreen is the offset in the sensor matrix of the solitary
1697 // green pixels (initialized here only to avoid compiler warning)
1698 char sgreen[2] = { 0 };
1699
1700 // Map a green hexagon around each non-green pixel and vice versa. allhex[] is later
1701 // indexed by the GPU kernels using tile-local (roi_in-relative) row/col taken mod 3, so
1702 // it must be built against that same absolute phase: pass roi_in here (not NULL) so the
1703 // color lookups land on the actual sensor position, not xtrans[][]'s own local origin.
1704 for(int row = 0; row < 3; row++)
1705 for(int col = 0; col < 3; col++)
1706 for(int ng = 0, d = 0; d < 10; d += 2)
1707 {
1708 const int g = FCxtrans(row, col, roi_in, xtrans) == 1;
1709 if(FCxtrans(row + orth[d] + 6, col + orth[d + 2] + 6, roi_in, xtrans) == 1)
1710 ng = 0;
1711 else
1712 ng++;
1713 // if there are four non-green pixels adjacent in cardinal
1714 // directions, this is the solitary green pixel
1715 if(ng == 4)
1716 {
1717 sgreen[0] = col;
1718 sgreen[1] = row;
1719 }
1720 if(ng == g + 1)
1721 for(int c = 0; c < 8; c++)
1722 {
1723 const int v = orth[d] * patt[g][c * 2] + orth[d + 1] * patt[g][c * 2 + 1];
1724 const int h = orth[d + 2] * patt[g][c * 2] + orth[d + 3] * patt[g][c * 2 + 1];
1725
1726 allhex[row][col][c ^ (g * 2 & d)][0] = h;
1727 allhex[row][col][c ^ (g * 2 & d)][1] = v;
1728 }
1729 }
1730
1731 dev_allhex = dt_opencl_copy_host_to_device_constant(devid, sizeof(allhex), allhex);
1732 if(IS_NULL_PTR(dev_allhex)) goto error;
1733
1734 for(int n = 0; n < ndir; n++)
1735 {
1736 dev_rgbv[n] = dt_opencl_alloc_device_buffer(devid, sizeof(float) * 4 * width * height);
1737 if(dev_rgbv[n] == NULL) goto error;
1738 }
1739
1740 dev_gminmax = dt_opencl_alloc_device_buffer(devid, sizeof(float) * 2 * width * height);
1741 if(IS_NULL_PTR(dev_gminmax)) goto error;
1742
1743 dev_aux = dt_opencl_alloc_device_buffer(devid, sizeof(float) * 4 * width * height);
1744 if(IS_NULL_PTR(dev_aux)) goto error;
1745
1746 dev_tmp = dev_out;
1747
1748 {
1749 // copy from dev_in to first rgb image buffer.
1750 size_t sizes[3] = { ROUNDUPDWD(width, devid), ROUNDUPDHT(height, devid), 1 };
1751 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_initial_copy, 0, sizeof(cl_mem), (void *)&dev_in);
1752 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_initial_copy, 1, sizeof(cl_mem), (void *)&dev_rgb[0]);
1753 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_initial_copy, 2, sizeof(int), (void *)&width);
1754 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_initial_copy, 3, sizeof(int), (void *)&height);
1755 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_initial_copy, 4, sizeof(int), (void *)&roi_in->x);
1756 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_initial_copy, 5, sizeof(int), (void *)&roi_in->y);
1757 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_initial_copy, 6, sizeof(cl_mem), (void *)&dev_xtrans);
1758 err = dt_opencl_enqueue_kernel_2d(devid, gd->kernel_markesteijn_initial_copy, sizes);
1759 if(err != CL_SUCCESS) goto error;
1760 }
1761
1762
1763 // duplicate dev_rgb[0] to dev_rgb[1], dev_rgb[2], and dev_rgb[3]
1764 for(int c = 1; c <= 3; c++)
1765 {
1766 err = dt_opencl_enqueue_copy_buffer_to_buffer(devid, dev_rgb[0], dev_rgb[c], 0, 0,
1767 sizeof(float) * 4 * width * height);
1768 if(err != CL_SUCCESS) goto error;
1769 }
1770
1771 // find minimum and maximum allowed green values of red/blue pixel pairs
1772 const int pad_g1_g3 = 3;
1773 dt_opencl_local_buffer_t locopt_g1_g3
1774 = (dt_opencl_local_buffer_t){ .xoffset = 2*3, .xfactor = 1, .yoffset = 2*3, .yfactor = 1,
1775 .cellsize = 1 * sizeof(float), .overhead = 0,
1776 .sizex = 1 << 8, .sizey = 1 << 8 };
1777
1778 if(!dt_opencl_local_buffer_opt(devid, gd->kernel_markesteijn_green_minmax, &locopt_g1_g3))
1779 goto error;
1780
1781 {
1782 size_t sizes[3] = { ROUNDUP(width, locopt_g1_g3.sizex), ROUNDUP(height, locopt_g1_g3.sizey), 1 };
1783 size_t local[3] = { locopt_g1_g3.sizex, locopt_g1_g3.sizey, 1 };
1784 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_green_minmax, 0, sizeof(cl_mem), (void *)&dev_rgb[0]);
1785 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_green_minmax, 1, sizeof(cl_mem), (void *)&dev_gminmax);
1786 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_green_minmax, 2, sizeof(int), (void *)&width);
1787 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_green_minmax, 3, sizeof(int), (void *)&height);
1788 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_green_minmax, 4, sizeof(int), (void *)&pad_g1_g3);
1789 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_green_minmax, 5, sizeof(int), (void *)&roi_in->x);
1790 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_green_minmax, 6, sizeof(int), (void *)&roi_in->y);
1791 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_green_minmax, 7, 2 * sizeof(char), (void *)sgreen);
1792 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_green_minmax, 8, sizeof(cl_mem), (void *)&dev_xtrans);
1793 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_green_minmax, 9, sizeof(cl_mem), (void *)&dev_allhex);
1794 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_green_minmax, 10,
1795 sizeof(float) * (locopt_g1_g3.sizex + 2*3) * (locopt_g1_g3.sizey + 2*3), NULL);
1796 err = dt_opencl_enqueue_kernel_2d_with_local(devid, gd->kernel_markesteijn_green_minmax, sizes, local);
1797 if(err != CL_SUCCESS) goto error;
1798 }
1799
1800 // interpolate green horizontally, vertically, and along both diagonals
1801 const int pad_g_interp = 3;
1802 dt_opencl_local_buffer_t locopt_g_interp
1803 = (dt_opencl_local_buffer_t){ .xoffset = 2*6, .xfactor = 1, .yoffset = 2*6, .yfactor = 1,
1804 .cellsize = 4 * sizeof(float), .overhead = 0,
1805 .sizex = 1 << 8, .sizey = 1 << 8 };
1806
1807 if(!dt_opencl_local_buffer_opt(devid, gd->kernel_markesteijn_interpolate_green, &locopt_g_interp))
1808 goto error;
1809
1810 {
1811 size_t sizes[3] = { ROUNDUP(width, locopt_g_interp.sizex), ROUNDUP(height, locopt_g_interp.sizey), 1 };
1812 size_t local[3] = { locopt_g_interp.sizex, locopt_g_interp.sizey, 1 };
1813 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_interpolate_green, 0, sizeof(cl_mem), (void *)&dev_rgb[0]);
1814 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_interpolate_green, 1, sizeof(cl_mem), (void *)&dev_rgb[1]);
1815 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_interpolate_green, 2, sizeof(cl_mem), (void *)&dev_rgb[2]);
1816 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_interpolate_green, 3, sizeof(cl_mem), (void *)&dev_rgb[3]);
1817 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_interpolate_green, 4, sizeof(cl_mem), (void *)&dev_gminmax);
1818 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_interpolate_green, 5, sizeof(int), (void *)&width);
1819 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_interpolate_green, 6, sizeof(int), (void *)&height);
1820 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_interpolate_green, 7, sizeof(int), (void *)&pad_g_interp);
1821 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_interpolate_green, 8, sizeof(int), (void *)&roi_in->x);
1822 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_interpolate_green, 9, sizeof(int), (void *)&roi_in->y);
1823 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_interpolate_green, 10, 2 * sizeof(char), (void *)sgreen);
1824 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_interpolate_green, 11, sizeof(cl_mem), (void *)&dev_xtrans);
1825 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_interpolate_green, 12, sizeof(cl_mem), (void *)&dev_allhex);
1826 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_interpolate_green, 13,
1827 sizeof(float) * 4 * (locopt_g_interp.sizex + 2*6) * (locopt_g_interp.sizey + 2*6), NULL);
1828 err = dt_opencl_enqueue_kernel_2d_with_local(devid, gd->kernel_markesteijn_interpolate_green, sizes, local);
1829 if(err != CL_SUCCESS) goto error;
1830 }
1831
1832 // multi-pass loop: one pass for Markesteijn-1 and three passes for Markesteijn-3
1833 for(int pass = 0; pass < passes; pass++)
1834 {
1835
1836 // if on second pass, copy rgb[0] to [3] into rgb[4] to [7] ....
1837 if(pass == 1)
1838 {
1839 for(int c = 0; c < 4; c++)
1840 {
1841 err = dt_opencl_enqueue_copy_buffer_to_buffer(devid, dev_rgb[c], dev_rgb[c + 4], 0, 0,
1842 sizeof(float) * 4 * width * height);
1843 if(err != CL_SUCCESS) goto error;
1844 }
1845 // ... and process that second set of buffers
1846 dev_rgb += 4;
1847 }
1848
1849 // second and third pass (only Markesteijn-3)
1850 if(pass)
1851 {
1852 // recalculate green from interpolated values of closer pixels
1853 const int pad_g_recalc = 6;
1854 size_t sizes[3] = { ROUNDUPDWD(width, devid), ROUNDUPDHT(height, devid), 1 };
1855 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_recalculate_green, 0, sizeof(cl_mem), (void *)&dev_rgb[0]);
1856 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_recalculate_green, 1, sizeof(cl_mem), (void *)&dev_rgb[1]);
1857 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_recalculate_green, 2, sizeof(cl_mem), (void *)&dev_rgb[2]);
1858 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_recalculate_green, 3, sizeof(cl_mem), (void *)&dev_rgb[3]);
1859 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_recalculate_green, 4, sizeof(cl_mem), (void *)&dev_gminmax);
1860 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_recalculate_green, 5, sizeof(int), (void *)&width);
1861 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_recalculate_green, 6, sizeof(int), (void *)&height);
1862 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_recalculate_green, 7, sizeof(int), (void *)&pad_g_recalc);
1863 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_recalculate_green, 8, sizeof(int), (void *)&roi_in->x);
1864 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_recalculate_green, 9, sizeof(int), (void *)&roi_in->y);
1865 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_recalculate_green, 10, 2 * sizeof(char), (void *)sgreen);
1866 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_recalculate_green, 11, sizeof(cl_mem), (void *)&dev_xtrans);
1867 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_recalculate_green, 12, sizeof(cl_mem), (void *)&dev_allhex);
1868 err = dt_opencl_enqueue_kernel_2d(devid, gd->kernel_markesteijn_recalculate_green, sizes);
1869 if(err != CL_SUCCESS) goto error;
1870 }
1871
1872 // interpolate red and blue values for solitary green pixels
1873 const int pad_rb_g = (passes == 1) ? 6 : 5;
1874 dt_opencl_local_buffer_t locopt_rb_g
1875 = (dt_opencl_local_buffer_t){ .xoffset = 2*2, .xfactor = 1, .yoffset = 2*2, .yfactor = 1,
1876 .cellsize = 4 * sizeof(float), .overhead = 0,
1877 .sizex = 1 << 8, .sizey = 1 << 8 };
1878
1879 if(!dt_opencl_local_buffer_opt(devid, gd->kernel_markesteijn_solitary_green, &locopt_rb_g))
1880 goto error;
1881
1882 cl_mem *dev_trgb = dev_rgb;
1883 for(int d = 0, i = 1, h = 0; d < 6; d++, i ^= 1, h ^= 2)
1884 {
1885 const char dir[2] = { i, i ^ 1 };
1886
1887 // we use dev_aux to transport intermediate results from one loop run to the next
1888 size_t sizes[3] = { ROUNDUP(width, locopt_rb_g.sizex), ROUNDUP(height, locopt_rb_g.sizey), 1 };
1889 size_t local[3] = { locopt_rb_g.sizex, locopt_rb_g.sizey, 1 };
1890 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_solitary_green, 0, sizeof(cl_mem), (void *)&dev_trgb[0]);
1891 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_solitary_green, 1, sizeof(cl_mem), (void *)&dev_aux);
1892 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_solitary_green, 2, sizeof(int), (void *)&width);
1893 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_solitary_green, 3, sizeof(int), (void *)&height);
1894 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_solitary_green, 4, sizeof(int), (void *)&pad_rb_g);
1895 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_solitary_green, 5, sizeof(int), (void *)&roi_in->x);
1896 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_solitary_green, 6, sizeof(int), (void *)&roi_in->y);
1897 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_solitary_green, 7, sizeof(int), (void *)&d);
1898 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_solitary_green, 8, 2 * sizeof(char), (void *)dir);
1899 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_solitary_green, 9, sizeof(int), (void *)&h);
1900 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_solitary_green, 10, 2 * sizeof(char), (void *)sgreen);
1901 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_solitary_green, 11, sizeof(cl_mem), (void *)&dev_xtrans);
1902 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_solitary_green, 12,
1903 sizeof(float) * 4 * (locopt_rb_g.sizex + 2*2) * (locopt_rb_g.sizey + 2*2), NULL);
1904 err = dt_opencl_enqueue_kernel_2d_with_local(devid, gd->kernel_markesteijn_solitary_green, sizes, local);
1905 if(err != CL_SUCCESS) goto error;
1906
1907 if((d < 2) || (d & 1)) dev_trgb++;
1908 }
1909
1910 // interpolate red for blue pixels and vice versa
1911 const int pad_rb_br = (passes == 1) ? 6 : 5;
1912 dt_opencl_local_buffer_t locopt_rb_br
1913 = (dt_opencl_local_buffer_t){ .xoffset = 2*3, .xfactor = 1, .yoffset = 2*3, .yfactor = 1,
1914 .cellsize = 4 * sizeof(float), .overhead = 0,
1915 .sizex = 1 << 8, .sizey = 1 << 8 };
1916
1917 if(!dt_opencl_local_buffer_opt(devid, gd->kernel_markesteijn_red_and_blue, &locopt_rb_br))
1918 goto error;
1919
1920 for(int d = 0; d < 4; d++)
1921 {
1922 size_t sizes[3] = { ROUNDUP(width, locopt_rb_br.sizex), ROUNDUP(height, locopt_rb_br.sizey), 1 };
1923 size_t local[3] = { locopt_rb_br.sizex, locopt_rb_br.sizey, 1 };
1924 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_red_and_blue, 0, sizeof(cl_mem), (void *)&dev_rgb[d]);
1925 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_red_and_blue, 1, sizeof(int), (void *)&width);
1926 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_red_and_blue, 2, sizeof(int), (void *)&height);
1927 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_red_and_blue, 3, sizeof(int), (void *)&pad_rb_br);
1928 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_red_and_blue, 4, sizeof(int), (void *)&roi_in->x);
1929 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_red_and_blue, 5, sizeof(int), (void *)&roi_in->y);
1930 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_red_and_blue, 6, sizeof(int), (void *)&d);
1931 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_red_and_blue, 7, 2 * sizeof(char), (void *)sgreen);
1932 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_red_and_blue, 8, sizeof(cl_mem), (void *)&dev_xtrans);
1933 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_red_and_blue, 9,
1934 sizeof(float) * 4 * (locopt_rb_br.sizex + 2*3) * (locopt_rb_br.sizey + 2*3), NULL);
1935 err = dt_opencl_enqueue_kernel_2d_with_local(devid, gd->kernel_markesteijn_red_and_blue, sizes, local);
1936 if(err != CL_SUCCESS) goto error;
1937 }
1938
1939 // interpolate red and blue for 2x2 blocks of green
1940 const int pad_g22 = (passes == 1) ? 8 : 4;
1941 dt_opencl_local_buffer_t locopt_g22
1942 = (dt_opencl_local_buffer_t){ .xoffset = 2*2, .xfactor = 1, .yoffset = 2*2, .yfactor = 1,
1943 .cellsize = 4 * sizeof(float), .overhead = 0,
1944 .sizex = 1 << 8, .sizey = 1 << 8 };
1945
1946 if(!dt_opencl_local_buffer_opt(devid, gd->kernel_markesteijn_interpolate_twoxtwo, &locopt_g22))
1947 goto error;
1948
1949 for(int d = 0, n = 0; d < ndir; d += 2, n++)
1950 {
1951 size_t sizes[3] = { ROUNDUP(width, locopt_g22.sizex), ROUNDUP(height, locopt_g22.sizey), 1 };
1952 size_t local[3] = { locopt_g22.sizex, locopt_g22.sizey, 1 };
1953 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_interpolate_twoxtwo, 0, sizeof(cl_mem), (void *)&dev_rgb[n]);
1954 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_interpolate_twoxtwo, 1, sizeof(int), (void *)&width);
1955 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_interpolate_twoxtwo, 2, sizeof(int), (void *)&height);
1956 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_interpolate_twoxtwo, 3, sizeof(int), (void *)&pad_g22);
1957 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_interpolate_twoxtwo, 4, sizeof(int), (void *)&roi_in->x);
1958 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_interpolate_twoxtwo, 5, sizeof(int), (void *)&roi_in->y);
1959 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_interpolate_twoxtwo, 6, sizeof(int), (void *)&d);
1960 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_interpolate_twoxtwo, 7, 2 * sizeof(char), (void *)sgreen);
1961 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_interpolate_twoxtwo, 8, sizeof(cl_mem), (void *)&dev_xtrans);
1962 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_interpolate_twoxtwo, 9, sizeof(cl_mem), (void *)&dev_allhex);
1963 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_interpolate_twoxtwo, 10,
1964 sizeof(float) * 4 * (locopt_g22.sizex + 2*2) * (locopt_g22.sizey + 2*2), NULL);
1965 err = dt_opencl_enqueue_kernel_2d_with_local(devid, gd->kernel_markesteijn_interpolate_twoxtwo, sizes, local);
1966 if(err != CL_SUCCESS) goto error;
1967 }
1968 }
1969 // end of multi pass
1970
1971 // gminmax data no longer needed
1972 dt_opencl_release_mem_object(dev_gminmax);
1973 dev_gminmax = NULL;
1974
1975 // jump back to the first set of rgb buffers (this is a noop for Markesteijn-1)
1976 dev_rgb = dev_rgbv;
1977
1978 // prepare derivatives buffers
1979 for(int n = 0; n < ndir; n++)
1980 {
1981 dev_drv[n] = dt_opencl_alloc_device_buffer(devid, sizeof(float) * width * height);
1982 if(dev_drv[n] == NULL) goto error;
1983 }
1984
1985 // convert to perceptual colorspace and differentiate in all directions
1986 const int pad_yuv = (passes == 1) ? 8 : 13;
1987 dt_opencl_local_buffer_t locopt_diff
1988 = (dt_opencl_local_buffer_t){ .xoffset = 2*1, .xfactor = 1, .yoffset = 2*1, .yfactor = 1,
1989 .cellsize = 4 * sizeof(float), .overhead = 0,
1990 .sizex = 1 << 8, .sizey = 1 << 8 };
1991
1992 if(!dt_opencl_local_buffer_opt(devid, gd->kernel_markesteijn_differentiate, &locopt_diff))
1993 goto error;
1994
1995 for(int d = 0; d < ndir; d++)
1996 {
1997 // convert to perceptual YPbPr colorspace
1998 size_t sizes_yuv[3] = { ROUNDUPDWD(width, devid), ROUNDUPDHT(height, devid), 1 };
1999 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_convert_yuv, 0, sizeof(cl_mem), (void *)&dev_rgb[d]);
2000 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_convert_yuv, 1, sizeof(cl_mem), (void *)&dev_aux);
2001 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_convert_yuv, 2, sizeof(int), (void *)&width);
2002 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_convert_yuv, 3, sizeof(int), (void *)&height);
2003 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_convert_yuv, 4, sizeof(int), (void *)&pad_yuv);
2004 err = dt_opencl_enqueue_kernel_2d(devid, gd->kernel_markesteijn_convert_yuv, sizes_yuv);
2005 if(err != CL_SUCCESS) goto error;
2006
2007
2008 // differentiate in all directions
2009 size_t sizes_diff[3] = { ROUNDUP(width, locopt_diff.sizex), ROUNDUP(height, locopt_diff.sizey), 1 };
2010 size_t local_diff[3] = { locopt_diff.sizex, locopt_diff.sizey, 1 };
2011 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_differentiate, 0, sizeof(cl_mem), (void *)&dev_aux);
2012 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_differentiate, 1, sizeof(cl_mem), (void *)&dev_drv[d]);
2013 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_differentiate, 2, sizeof(int), (void *)&width);
2014 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_differentiate, 3, sizeof(int), (void *)&height);
2015 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_differentiate, 4, sizeof(int), (void *)&pad_yuv);
2016 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_differentiate, 5, sizeof(int), (void *)&d);
2017 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_differentiate, 6,
2018 sizeof(float) * 4 * (locopt_diff.sizex + 2*1) * (locopt_diff.sizey + 2*1), NULL);
2019 err = dt_opencl_enqueue_kernel_2d_with_local(devid, gd->kernel_markesteijn_differentiate, sizes_diff, local_diff);
2020 if(err != CL_SUCCESS) goto error;
2021 }
2022
2023 // reserve buffers for homogeneity maps and sum maps
2024 for(int n = 0; n < ndir; n++)
2025 {
2026 dev_homo[n] = dt_opencl_alloc_device_buffer(devid, sizeof(unsigned char) * width * height);
2027 if(dev_homo[n] == NULL) goto error;
2028
2029 dev_homosum[n] = dt_opencl_alloc_device_buffer(devid, sizeof(unsigned char) * width * height);
2030 if(dev_homosum[n] == NULL) goto error;
2031 }
2032
2033 // get thresholds for homogeneity map (store them in dev_aux)
2034 for(int d = 0; d < ndir; d++)
2035 {
2036 const int pad_homo = (passes == 1) ? 10 : 15;
2037 size_t sizes[3] = { ROUNDUPDWD(width, devid), ROUNDUPDHT(height, devid), 1 };
2038 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_homo_threshold, 0, sizeof(cl_mem), (void *)&dev_drv[d]);
2039 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_homo_threshold, 1, sizeof(cl_mem), (void *)&dev_aux);
2040 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_homo_threshold, 2, sizeof(int), (void *)&width);
2041 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_homo_threshold, 3, sizeof(int), (void *)&height);
2042 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_homo_threshold, 4, sizeof(int), (void *)&pad_homo);
2043 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_homo_threshold, 5, sizeof(int), (void *)&d);
2044 err = dt_opencl_enqueue_kernel_2d(devid, gd->kernel_markesteijn_homo_threshold, sizes);
2045 if(err != CL_SUCCESS) goto error;
2046 }
2047
2048 // set homogeneity maps
2049 const int pad_homo = (passes == 1) ? 10 : 15;
2050 dt_opencl_local_buffer_t locopt_homo
2051 = (dt_opencl_local_buffer_t){ .xoffset = 2*1, .xfactor = 1, .yoffset = 2*1, .yfactor = 1,
2052 .cellsize = 1 * sizeof(float), .overhead = 0,
2053 .sizex = 1 << 8, .sizey = 1 << 8 };
2054
2055 if(!dt_opencl_local_buffer_opt(devid, gd->kernel_markesteijn_homo_set, &locopt_homo))
2056 goto error;
2057
2058 for(int d = 0; d < ndir; d++)
2059 {
2060 size_t sizes[3] = { ROUNDUP(width, locopt_homo.sizex),ROUNDUP(height, locopt_homo.sizey), 1 };
2061 size_t local[3] = { locopt_homo.sizex, locopt_homo.sizey, 1 };
2062 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_homo_set, 0, sizeof(cl_mem), (void *)&dev_drv[d]);
2063 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_homo_set, 1, sizeof(cl_mem), (void *)&dev_aux);
2064 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_homo_set, 2, sizeof(cl_mem), (void *)&dev_homo[d]);
2065 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_homo_set, 3, sizeof(int), (void *)&width);
2066 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_homo_set, 4, sizeof(int), (void *)&height);
2067 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_homo_set, 5, sizeof(int), (void *)&pad_homo);
2068 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_homo_set, 6,
2069 sizeof(float) * (locopt_homo.sizex + 2*1) * (locopt_homo.sizey + 2*1), NULL);
2070 err = dt_opencl_enqueue_kernel_2d_with_local(devid, gd->kernel_markesteijn_homo_set, sizes, local);
2071 if(err != CL_SUCCESS) goto error;
2072 }
2073
2074 // get rid of dev_drv buffers
2075 for(int n = 0; n < 8; n++)
2076 {
2077 dt_opencl_release_mem_object(dev_drv[n]);
2078 dev_drv[n] = NULL;
2079 }
2080
2081 // build 5x5 sum of homogeneity maps for each pixel and direction
2082 dt_opencl_local_buffer_t locopt_homo_sum
2083 = (dt_opencl_local_buffer_t){ .xoffset = 2*2, .xfactor = 1, .yoffset = 2*2, .yfactor = 1,
2084 .cellsize = 1 * sizeof(float), .overhead = 0,
2085 .sizex = 1 << 8, .sizey = 1 << 8 };
2086
2087 if(!dt_opencl_local_buffer_opt(devid, gd->kernel_markesteijn_homo_sum, &locopt_homo_sum))
2088 goto error;
2089
2090 for(int d = 0; d < ndir; d++)
2091 {
2092 size_t sizes[3] = { ROUNDUP(width, locopt_homo_sum.sizex), ROUNDUP(height, locopt_homo_sum.sizey), 1 };
2093 size_t local[3] = { locopt_homo_sum.sizex, locopt_homo_sum.sizey, 1 };
2094 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_homo_sum, 0, sizeof(cl_mem), (void *)&dev_homo[d]);
2095 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_homo_sum, 1, sizeof(cl_mem), (void *)&dev_homosum[d]);
2096 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_homo_sum, 2, sizeof(int), (void *)&width);
2097 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_homo_sum, 3, sizeof(int), (void *)&height);
2098 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_homo_sum, 4, sizeof(int), (void *)&pad_tile);
2099 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_homo_sum, 5,
2100 sizeof(char) * (locopt_homo_sum.sizex + 2*2) * (locopt_homo_sum.sizey + 2*2), NULL);
2101 err = dt_opencl_enqueue_kernel_2d_with_local(devid, gd->kernel_markesteijn_homo_sum, sizes, local);
2102 if(err != CL_SUCCESS) goto error;
2103 }
2104
2105 // get maximum of homogeneity maps (store in dev_aux)
2106 for(int d = 0; d < ndir; d++)
2107 {
2108 size_t sizes[3] = { ROUNDUPDWD(width, devid), ROUNDUPDHT(height, devid), 1 };
2109 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_homo_max, 0, sizeof(cl_mem), (void *)&dev_homosum[d]);
2110 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_homo_max, 1, sizeof(cl_mem), (void *)&dev_aux);
2111 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_homo_max, 2, sizeof(int), (void *)&width);
2112 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_homo_max, 3, sizeof(int), (void *)&height);
2113 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_homo_max, 4, sizeof(int), (void *)&pad_tile);
2114 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_homo_max, 5, sizeof(int), (void *)&d);
2115 err = dt_opencl_enqueue_kernel_2d(devid, gd->kernel_markesteijn_homo_max, sizes);
2116 if(err != CL_SUCCESS) goto error;
2117 }
2118
2119 {
2120 // adjust maximum value
2121 size_t sizes[3] = { ROUNDUPDWD(width, devid), ROUNDUPDHT(height, devid), 1 };
2122 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_homo_max_corr, 0, sizeof(cl_mem), (void *)&dev_aux);
2123 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_homo_max_corr, 1, sizeof(int), (void *)&width);
2124 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_homo_max_corr, 2, sizeof(int), (void *)&height);
2125 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_homo_max_corr, 3, sizeof(int), (void *)&pad_tile);
2126 err = dt_opencl_enqueue_kernel_2d(devid, gd->kernel_markesteijn_homo_max_corr, sizes);
2127 if(err != CL_SUCCESS) goto error;
2128 }
2129
2130 // for Markesteijn-3: use only one of two directions if there is a difference in homogeneity
2131 for(int d = 0; d < ndir - 4; d++)
2132 {
2133 size_t sizes[3] = { ROUNDUPDWD(width, devid), ROUNDUPDHT(height, devid), 1 };
2134 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_homo_quench, 0, sizeof(cl_mem), (void *)&dev_homosum[d]);
2135 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_homo_quench, 1, sizeof(cl_mem), (void *)&dev_homosum[d + 4]);
2136 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_homo_quench, 2, sizeof(int), (void *)&width);
2137 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_homo_quench, 3, sizeof(int), (void *)&height);
2138 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_homo_quench, 4, sizeof(int), (void *)&pad_tile);
2139 err = dt_opencl_enqueue_kernel_2d(devid, gd->kernel_markesteijn_homo_quench, sizes);
2140 if(err != CL_SUCCESS) goto error;
2141 }
2142
2143 {
2144 // initialize output buffer to zero
2145 size_t sizes[3] = { ROUNDUPDWD(width, devid), ROUNDUPDHT(height, devid), 1 };
2146 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_zero, 0, sizeof(cl_mem), (void *)&dev_tmp);
2147 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_zero, 1, sizeof(int), (void *)&width);
2148 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_zero, 2, sizeof(int), (void *)&height);
2149 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_zero, 3, sizeof(int), (void *)&pad_tile);
2150 err = dt_opencl_enqueue_kernel_2d(devid, gd->kernel_markesteijn_zero, sizes);
2151 if(err != CL_SUCCESS) goto error;
2152 }
2153
2154 // need to get another temp buffer for the output image (may use the space of dev_drv[] freed earlier)
2155 dev_tmptmp = dt_opencl_alloc_device(devid, (size_t)width, height, sizeof(float) * 4);
2156 if(IS_NULL_PTR(dev_tmptmp)) goto error;
2157
2158 cl_mem dev_t1 = dev_tmp;
2159 cl_mem dev_t2 = dev_tmptmp;
2160
2161 // accumulate all contributions
2162 for(int d = 0; d < ndir; d++)
2163 {
2164 size_t sizes[3] = { ROUNDUPDWD(width, devid), ROUNDUPDHT(height, devid), 1 };
2165 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_accu, 0, sizeof(cl_mem), (void *)&dev_t1);
2166 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_accu, 1, sizeof(cl_mem), (void *)&dev_t2);
2167 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_accu, 2, sizeof(cl_mem), (void *)&dev_rgbv[d]);
2168 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_accu, 3, sizeof(cl_mem), (void *)&dev_homosum[d]);
2169 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_accu, 4, sizeof(cl_mem), (void *)&dev_aux);
2170 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_accu, 5, sizeof(int), (void *)&width);
2171 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_accu, 6, sizeof(int), (void *)&height);
2172 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_accu, 7, sizeof(int), (void *)&pad_tile);
2173 err = dt_opencl_enqueue_kernel_2d(devid, gd->kernel_markesteijn_accu, sizes);
2174 if(err != CL_SUCCESS) goto error;
2175
2176 // swap buffers
2177 cl_mem dev_t = dev_t2;
2178 dev_t2 = dev_t1;
2179 dev_t1 = dev_t;
2180 }
2181
2182 // copy output to dev_tmptmp (if not already there)
2183 // note: we need to take swap of buffers into account, so current output lies in dev_t1
2184 if(dev_t1 != dev_tmptmp)
2185 {
2186 size_t origin[] = { 0, 0, 0 };
2187 size_t region[] = { width, height, 1 };
2188 err = dt_opencl_enqueue_copy_image(devid, dev_t1, dev_tmptmp, origin, origin, region);
2189 if(err != CL_SUCCESS) goto error;
2190 }
2191
2192 {
2193 // process the final image
2194 size_t sizes[3] = { ROUNDUPDWD(width, devid), ROUNDUPDHT(height, devid), 1 };
2195 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_final, 0, sizeof(cl_mem), (void *)&dev_tmptmp);
2196 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_final, 1, sizeof(cl_mem), (void *)&dev_tmp);
2197 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_final, 2, sizeof(int), (void *)&width);
2198 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_final, 3, sizeof(int), (void *)&height);
2199 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_final, 4, sizeof(int), (void *)&pad_tile);
2200 dt_opencl_set_kernel_arg(devid, gd->kernel_markesteijn_final, 5, 4*sizeof(float), (void *)processed_maximum);
2201 err = dt_opencl_enqueue_kernel_2d(devid, gd->kernel_markesteijn_final, sizes);
2202 if(err != CL_SUCCESS) goto error;
2203 }
2204
2205 // now it's time to get rid of most of the temporary buffers (except of dev_tmp and dev_xtrans)
2206 for(int n = 0; n < 8; n++)
2207 {
2208 dt_opencl_release_mem_object(dev_rgbv[n]);
2209 dev_rgbv[n] = NULL;
2210 }
2211
2212 for(int n = 0; n < 8; n++)
2213 {
2214 dt_opencl_release_mem_object(dev_homo[n]);
2215 dev_homo[n] = NULL;
2216 }
2217
2218 for(int n = 0; n < 8; n++)
2219 {
2220 dt_opencl_release_mem_object(dev_homosum[n]);
2221 dev_homosum[n] = NULL;
2222 }
2223
2224 dt_opencl_release_mem_object(dev_aux);
2225 dev_aux = NULL;
2226
2227 dt_opencl_release_mem_object(dev_xtrans);
2228 dev_xtrans = NULL;
2229
2230 dt_opencl_release_mem_object(dev_allhex);
2231 dev_allhex = NULL;
2232
2233 dt_opencl_release_mem_object(dev_green_eq);
2234 dev_green_eq = NULL;
2235
2236 dt_opencl_release_mem_object(dev_tmptmp);
2237 dev_tmptmp = NULL;
2238
2239 // take care of image borders. the algorithm above leaves an unprocessed border of pad_tile pixels.
2240 // strategy: take the four edges and process them each with process_vng_cl(). as VNG produces
2241 // an image with a border with only linear interpolation we process edges of pad_tile+3px and
2242 // drop 3px on the inner side if possible
2243
2244 // take care of some degenerate cases (which might happen if we are called in a tiling context)
2245 const int wd = (width > pad_tile+3) ? pad_tile+3 : width;
2246 const int ht = (height > pad_tile+3) ? pad_tile+3 : height;
2247 const int wdc = (wd >= pad_tile+3) ? 3 : 0;
2248 const int htc = (ht >= pad_tile+3) ? 3 : 0;
2249
2250 // the data of all four edges:
2251 // total edge: x-offset, y-offset, width, height,
2252 // after dropping: x-offset adjust, y-offset adjust, width adjust, height adjust
2253 const int edges[4][8] = { { 0, 0, wd, height, 0, 0, -wdc, 0 },
2254 { 0, 0, width, ht, 0, 0, 0, -htc },
2255 { width - wd, 0, wd, height, wdc, 0, -wdc, 0 },
2256 { 0, height - ht, width, ht, 0, htc, 0, -htc } };
2257
2258 for(int n = 0; n < 4; n++)
2259 {
2260 dt_iop_roi_t roi = { roi_in->x + edges[n][0], roi_in->y + edges[n][1], edges[n][2], edges[n][3], 1.0f };
2261
2262 size_t iorigin[] = { edges[n][0], edges[n][1], 0 };
2263 size_t oorigin[] = { 0, 0, 0 };
2264 size_t region[] = { edges[n][2], edges[n][3], 1 };
2265
2266 // reserve input buffer for image edge
2267 dev_edge_in = dt_opencl_alloc_device(devid, edges[n][2], edges[n][3], sizeof(float));
2268 if(IS_NULL_PTR(dev_edge_in)) goto error;
2269
2270 // reserve output buffer for VNG processing of edge
2271 dev_edge_out = dt_opencl_alloc_device(devid, edges[n][2], edges[n][3], sizeof(float) * 4);
2272 if(IS_NULL_PTR(dev_edge_out)) goto error;
2273
2274 // copy edge to input buffer
2275 err = dt_opencl_enqueue_copy_image(devid, dev_in, dev_edge_in, iorigin, oorigin, region);
2276 if(err != CL_SUCCESS) goto error;
2277
2278 // VNG processing
2279 if(!process_vng_cl(self, pipe, piece, dev_edge_in, dev_edge_out, &roi, &roi, smooth, FALSE))
2280 goto error;
2281
2282 // adjust for "good" part, dropping linear border where possible
2283 iorigin[0] += edges[n][4];
2284 iorigin[1] += edges[n][5];
2285 oorigin[0] += edges[n][4];
2286 oorigin[1] += edges[n][5];
2287 region[0] += edges[n][6];
2288 region[1] += edges[n][7];
2289
2290 // copy output
2291 err = dt_opencl_enqueue_copy_image(devid, dev_edge_out, dev_tmp, oorigin, iorigin, region);
2292 if(err != CL_SUCCESS) goto error;
2293
2294 // release intermediate buffers
2295 dt_opencl_release_mem_object(dev_edge_in);
2296 dt_opencl_release_mem_object(dev_edge_out);
2297 dev_edge_in = dev_edge_out = NULL;
2298 }
2299 // free remaining temporary buffers
2300 if(dev_tmp != dev_out) dt_opencl_release_mem_object(dev_tmp);
2301 dev_tmp = NULL;
2302
2303 dt_opencl_release_mem_object(dev_xtrans);
2304 dev_xtrans = NULL;
2305
2306
2307 // color smoothing
2308 if(data->color_smoothing)
2309 {
2310 if(!color_smoothing_cl(self, pipe, piece, dev_out, dev_out, roi_out, data->color_smoothing))
2311 goto error;
2312 }
2313
2314 return TRUE;
2315
2316error:
2317 if(dev_tmp != dev_out) dt_opencl_release_mem_object(dev_tmp);
2318
2319 for(int n = 0; n < 8; n++)
2320 dt_opencl_release_mem_object(dev_rgbv[n]);
2321 for(int n = 0; n < 8; n++)
2322 dt_opencl_release_mem_object(dev_drv[n]);
2323 for(int n = 0; n < 8; n++)
2324 dt_opencl_release_mem_object(dev_homo[n]);
2325 for(int n = 0; n < 8; n++)
2326 dt_opencl_release_mem_object(dev_homosum[n]);
2327 dt_opencl_release_mem_object(dev_gminmax);
2328 dt_opencl_release_mem_object(dev_tmptmp);
2329 dt_opencl_release_mem_object(dev_xtrans);
2330 dt_opencl_release_mem_object(dev_allhex);
2331 dt_opencl_release_mem_object(dev_green_eq);
2332 dt_opencl_release_mem_object(dev_aux);
2333 dt_opencl_release_mem_object(dev_edge_in);
2334 dt_opencl_release_mem_object(dev_edge_out);
2335 dt_print(DT_DEBUG_OPENCL, "[opencl_demosaic] couldn't enqueue kernel! %d\n", err);
2336 return FALSE;
2337}
2338
2339#endif // OPENCL
2340
2341#undef PIX_SWAP
2342#undef PIX_SORT
2343#undef CCLIP
2344#undef TS
2345
2346// clang-format off
2347// modelines: These editor modelines have been set for all relevant files by tools/update_modelines.py
2348// vim: shiftwidth=2 expandtab tabstop=2 cindent
2349// kate: tab-indents: off; indent-width 2; replace-tabs on; indent-mode cstyle; remove-trailing-spaces modified;
2350// clang-format on
error
static void error(char *msg)
Definition ashift_lsd.c:202
TRUE
#define TRUE
Definition ashift_lsd.c:162
FALSE
#define FALSE
Definition ashift_lsd.c:158
color_smoothing_cl
static int color_smoothing_cl(struct dt_iop_module_t *self, const dt_dev_pixelpipe_t *pipe, const dt_dev_pixelpipe_iop_t *piece, cl_mem dev_in, cl_mem dev_out, const dt_iop_roi_t *const roi_out, const int passes)
Definition basic.c:334
x
static const float x
Definition colorprofiles/iop_profile.h:560
f
const float f
Definition colorprofiles/iop_profile.h:541
v
const float v
Definition colorprofiles/iop_profile.h:536
i
const float i
Definition colorspaces_inline_conversions.h:441
rgb
static dt_aligned_pixel_t rgb
Definition colorspaces_inline_conversions.h:345
g
const float g
Definition colorspaces_inline_conversions.h:675
d
const float d
Definition colorspaces_inline_conversions.h:681
min
static const float const float const float min
Definition colorspaces_inline_conversions.h:439
max
const float max
Definition colorspaces_inline_conversions.h:491
out
const dt_colormatrix_t dt_aligned_pixel_t out
Definition colorspaces_inline_conversions.h:43
n
const float n
Definition colorspaces_inline_conversions.h:679
top
const float top
Definition colorspaces_inline_conversions.h:444
row
static const int row
Definition colorspaces_inline_conversions.h:36
dt_conf_get_int
int dt_conf_get_int(const char *name)
Integer for name, clamped to the bounds declared in the XML.
Definition common/conf.c:295
FCxtrans
static int FCxtrans(const int row, const int col, global const unsigned char(*const xtrans)[6])
Definition data/kernels/common.h:54
DT_IOP_DEMOSAIC_MARKESTEIJN_3
@ DT_IOP_DEMOSAIC_MARKESTEIJN_3
Definition demosaic.c:136
color
GdkRGBA color[]
Definition geotagging.c:541
DT_DEBUG_OPENCL
@ DT_DEBUG_OPENCL
Definition logging.h:57
dt_print
void dt_print(dt_debug_thread_t thread, const char *msg,...) __attribute__((format(printf
Print to stdout when thread is enabled, prefixed with seconds since startup.
IS_NULL_PTR
#define IS_NULL_PTR(p)
C is way too permissive with !=, == and if(var) checks, which can mean too many things depending on w...
Definition macros.h:96
TRANSLATE
#define TRANSLATE(n, size)
CONV_FILT
#define CONV_FILT(VAR, FILT)
process_markesteijn_cl
static int process_markesteijn_cl(struct dt_iop_module_t *self, const dt_dev_pixelpipe_t *pipe, const dt_dev_pixelpipe_iop_t *piece, cl_mem dev_in, cl_mem dev_out, const dt_iop_roi_t *const roi_in, const dt_iop_roi_t *const roi_out, const gboolean smooth)
Definition markesteijn.c:1648
hexmap
static const short * hexmap(const int row, const int col, short(*const allhex)[3][8])
Definition markesteijn.c:30
SQR
#define SQR(x)
Definition markesteijn.c:25
TS
#define TS
Definition markesteijn.c:27
PIX_SORT
#define PIX_SORT(a, b)
xtrans_markesteijn_interpolate
static __DT_CLONE_TARGETS__ void xtrans_markesteijn_interpolate(float *out, const float *const in, const dt_iop_roi_t *const roi_out, const dt_iop_roi_t *const roi_in, const uint8_t(*const xtrans)[6], const int passes)
Definition markesteijn.c:47
xtrans_fdc_interpolate
static __DT_CLONE_TARGETS__ void xtrans_fdc_interpolate(struct dt_iop_module_t *self, float *out, const float *const in, const dt_iop_roi_t *const roi_out, const dt_iop_roi_t *const roi_in, const uint8_t(*const xtrans)[6])
Definition markesteijn.c:532
CLAMPS
#define CLAMPS(A, L, H)
Definition math.h:78
memset_zero
static void memset_zero(void *const buffer, size_t size)
Set the memory buffer to zero as a pack of unsigned char.
Definition mem_alloc.h:270
width
uint32_t width
Definition mipmap_cache.c:0
height
uint32_t height
Definition mipmap_cache.c:1
dt_opencl_local_buffer_opt
int dt_opencl_local_buffer_opt(const int devid, const int kernel, dt_opencl_local_buffer_t *factors)
Definition opencl.c:3713
dt_opencl_enqueue_kernel_2d
int dt_opencl_enqueue_kernel_2d(const int dev, const int kernel, const size_t *sizes)
Definition opencl.c:2554
dt_opencl_alloc_device_buffer
void * dt_opencl_alloc_device_buffer(const int devid, const size_t size)
Definition opencl.c:2970
dt_opencl_alloc_device
void * dt_opencl_alloc_device(const int devid, const int width, const int height, const int bpp)
Definition opencl.c:2894
dt_opencl_copy_host_to_device_constant
void * dt_opencl_copy_host_to_device_constant(const int devid, const size_t size, void *host)
Definition opencl.c:2750
dt_opencl_enqueue_copy_image
int dt_opencl_enqueue_copy_image(const int devid, cl_mem src, cl_mem dst, size_t *orig_src, size_t *orig_dst, size_t *region)
Definition opencl.c:2679
dt_opencl_set_kernel_arg
int dt_opencl_set_kernel_arg(const int dev, const int kernel, const int num, const size_t size, const void *arg)
Definition opencl.c:2545
dt_opencl_enqueue_kernel_2d_with_local
int dt_opencl_enqueue_kernel_2d_with_local(const int dev, const int kernel, const size_t *sizes, const size_t *local)
Definition opencl.c:2560
dt_opencl_enqueue_copy_buffer_to_buffer
int dt_opencl_enqueue_copy_buffer_to_buffer(const int devid, cl_mem src_buffer, cl_mem dst_buffer, size_t srcoffset, size_t dstoffset, size_t size)
Definition opencl.c:2714
dt_opencl_release_mem_object
void dt_opencl_release_mem_object(cl_mem mem)
Definition opencl.c:2805
ROUNDUP
#define ROUNDUP(a, n)
Definition opencl.h:82
ROUNDUPDHT
#define ROUNDUPDHT(a, b)
Definition opencl.h:86
ROUNDUPDWD
#define ROUNDUPDWD(a, b)
Definition opencl.h:85
__OMP_PARALLEL_FOR__
#define __OMP_PARALLEL_FOR__(...)
Definition openmp.h:95
dt_pixelpipe_cache_alloc_perthread
#define dt_pixelpipe_cache_alloc_perthread(n, objsize, padded_size)
Definition pixelpipe_cache_alloc.h:126
dt_pixelpipe_cache_free_align
#define dt_pixelpipe_cache_free_align(mem)
Definition pixelpipe_cache_alloc.h:106
dt_get_perthread
#define dt_get_perthread(buf, padsize)
Definition pixelpipe_cache_alloc.h:155
ABS
#define ABS(n)
Definition polar_decomposition.h:82
dt_aligned_pixel_t
DT_ALIGNED_PIXEL float dt_aligned_pixel_t[4]
Definition simd.h:53
dt_dev_pixelpipe_iop_t
Definition pixelpipe_hb.h:102
dt_dev_pixelpipe_iop_t::dsc_in
dt_iop_buffer_dsc_t dsc_in
Definition pixelpipe_hb.h:148
dt_dev_pixelpipe_iop_t::data
struct dt_iop_module_t *void * data
Definition pixelpipe_hb.h:103
dt_dev_pixelpipe_t
Definition pixelpipe_hb.h:305
dt_dev_pixelpipe_t::devid
int devid
Definition pixelpipe_hb.h:458
dt_develop_t::image_storage
dt_image_t image_storage
Definition develop.h:225
dt_image_t::exif_iso
float exif_iso
Definition image.h:370
dt_iop_buffer_dsc_t::xtrans
uint8_t xtrans[6][6]
Definition format.h:99
dt_iop_buffer_dsc_t::processed_maximum
dt_aligned_pixel_t processed_maximum
Definition format.h:114
dt_iop_demosaic_data_t
Definition demosaic.c:239
dt_iop_demosaic_data_t::demosaicing_method
uint32_t demosaicing_method
Definition demosaic.c:242
dt_iop_demosaic_data_t::color_smoothing
uint32_t color_smoothing
Definition demosaic.c:241
dt_iop_demosaic_global_data_t
Definition demosaic.c:173
dt_iop_demosaic_global_data_t::kernel_markesteijn_accu
int kernel_markesteijn_accu
Definition demosaic.c:219
dt_iop_demosaic_global_data_t::kernel_markesteijn_solitary_green
int kernel_markesteijn_solitary_green
Definition demosaic.c:206
dt_iop_demosaic_global_data_t::kernel_markesteijn_homo_quench
int kernel_markesteijn_homo_quench
Definition demosaic.c:217
dt_iop_demosaic_global_data_t::kernel_markesteijn_initial_copy
int kernel_markesteijn_initial_copy
Definition demosaic.c:203
dt_iop_demosaic_global_data_t::kernel_markesteijn_final
int kernel_markesteijn_final
Definition demosaic.c:220
dt_iop_demosaic_global_data_t::kernel_markesteijn_homo_max_corr
int kernel_markesteijn_homo_max_corr
Definition demosaic.c:216
dt_iop_demosaic_global_data_t::kernel_markesteijn_interpolate_green
int kernel_markesteijn_interpolate_green
Definition demosaic.c:205
dt_iop_demosaic_global_data_t::kernel_markesteijn_convert_yuv
int kernel_markesteijn_convert_yuv
Definition demosaic.c:210
dt_iop_demosaic_global_data_t::kernel_markesteijn_homo_threshold
int kernel_markesteijn_homo_threshold
Definition demosaic.c:212
dt_iop_demosaic_global_data_t::kernel_markesteijn_homo_sum
int kernel_markesteijn_homo_sum
Definition demosaic.c:214
dt_iop_demosaic_global_data_t::kernel_markesteijn_green_minmax
int kernel_markesteijn_green_minmax
Definition demosaic.c:204
dt_iop_demosaic_global_data_t::kernel_markesteijn_recalculate_green
int kernel_markesteijn_recalculate_green
Definition demosaic.c:207
dt_iop_demosaic_global_data_t::kernel_markesteijn_interpolate_twoxtwo
int kernel_markesteijn_interpolate_twoxtwo
Definition demosaic.c:209
dt_iop_demosaic_global_data_t::kernel_markesteijn_homo_set
int kernel_markesteijn_homo_set
Definition demosaic.c:213
dt_iop_demosaic_global_data_t::kernel_markesteijn_red_and_blue
int kernel_markesteijn_red_and_blue
Definition demosaic.c:208
dt_iop_demosaic_global_data_t::kernel_markesteijn_differentiate
int kernel_markesteijn_differentiate
Definition demosaic.c:211
dt_iop_demosaic_global_data_t::kernel_markesteijn_homo_max
int kernel_markesteijn_homo_max
Definition demosaic.c:215
dt_iop_demosaic_global_data_t::kernel_markesteijn_zero
int kernel_markesteijn_zero
Definition demosaic.c:218
dt_iop_module_t
Definition imageop.h:251
dt_iop_module_t::dev
struct dt_develop_t * dev
Definition imageop.h:311
dt_iop_module_t::global_data
dt_iop_global_data_t * global_data
Definition imageop.h:337
dt_iop_roi_t
Region of interest passed through the pixelpipe.
Definition format.h:49
dt_iop_roi_t::x
int x
Definition format.h:50
dt_iop_roi_t::width
int width
Definition format.h:50
dt_iop_roi_t::height
int height
Definition format.h:50
dt_iop_roi_t::y
int y
Definition format.h:50
dt_opencl_local_buffer_t
Definition opencl.h:325
dt_opencl_local_buffer_t::xoffset
const int xoffset
Definition opencl.h:326
dt_opencl_local_buffer_t::sizey
int sizey
Definition opencl.h:333
dt_opencl_local_buffer_t::sizex
int sizex
Definition opencl.h:332
__DT_CLONE_TARGETS__
#define __DT_CLONE_TARGETS__
Definition target_clones.h:101
MIN
#define MIN(a, b)
Definition thinplate.c:32
process_vng_cl
static int process_vng_cl(struct dt_iop_module_t *self, const dt_dev_pixelpipe_t *pipe, const dt_dev_pixelpipe_iop_t *piece, cl_mem dev_in, cl_mem dev_out, const dt_iop_roi_t *const roi_in, const dt_iop_roi_t *const roi_out, const gboolean smooth, const int only_vng_linear)
Definition vng.c:206