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