1 /*
2 * UFRaw - Unidentified Flying Raw converter for digital camera images
3 *
4 * dcraw_indi.c - DCRaw functions made independent
5 * Copyright 2004-2016 by Udi Fuchs
6 *
7 * based on dcraw by Dave Coffin
8 * http://www.cybercom.net/~dcoffin/
9 *
10 * This program is free software; you can redistribute it and/or modify
11 * it under the terms of the GNU General Public License as published by
12 * the Free Software Foundation; either version 2 of the License, or
13 * (at your option) any later version.
14 */
15
16 #ifdef HAVE_CONFIG_H
17 #include "config.h"
18 #endif
19
20 #include <math.h>
21 #include <stdio.h>
22 #include <stdlib.h>
23 #include <string.h>
24 #include <time.h>
25 #include <glib.h>
26 #include <glib/gi18n.h> /*For _(String) definition - NKBJ*/
27 #include "dcraw_api.h"
28 #include "uf_progress.h"
29
30 #ifdef _OPENMP
31 #include <omp.h>
32 #define uf_omp_get_thread_num() omp_get_thread_num()
33 #define uf_omp_get_num_threads() omp_get_num_threads()
34 #else
35 #define uf_omp_get_thread_num() 0
36 #define uf_omp_get_num_threads() 1
37 #endif
38
39 #if !defined(ushort)
40 #define ushort unsigned short
41 #endif
42
43 extern const double xyz_rgb[3][3];
44 extern const float d65_white[3];
45
46 #define CLASS
47
48 #define FORC(cnt) for (c=0; c < cnt; c++)
49 #define FORC3 FORC(3)
50 #define FORC4 FORC(4)
51 #define FORCC FORC(colors)
52
53 #define SQR(x) ((x)*(x))
54 #define LIM(x,min,max) MAX(min,MIN(x,max))
55 #define ULIM(x,y,z) ((y) < (z) ? LIM(x,y,z) : LIM(x,z,y))
56 #define CLIP(x) LIM((int)(x),0,65535)
57 #define SWAP(a,b) { a ^= b; a ^= (b ^= a); }
58
59 /*
60 In order to inline this calculation, I make the risky
61 assumption that all filter patterns can be described
62 by a repeating pattern of eight rows and two columns
63
64 Return values are either 0/1/2/3 = G/M/C/Y or 0/1/2/3 = R/G1/B/G2
65 */
66 #define FC(row,col) \
67 (int)(filters >> ((((row) << 1 & 14) + ((col) & 1)) << 1) & 3)
68
69 #define BAYER(row,col) \
70 image[((row) >> shrink)*iwidth + ((col) >> shrink)][FC(row,col)]
71
fcol_INDI(const unsigned filters,const int row,const int col,const int top_margin,const int left_margin,char xtrans[6][6])72 int CLASS fcol_INDI(const unsigned filters, const int row, const int col,
73 const int top_margin, const int left_margin,
74 /*const*/ char xtrans[6][6])
75 {
76 static const char filter[16][16] = {
77 { 2, 1, 1, 3, 2, 3, 2, 0, 3, 2, 3, 0, 1, 2, 1, 0 },
78 { 0, 3, 0, 2, 0, 1, 3, 1, 0, 1, 1, 2, 0, 3, 3, 2 },
79 { 2, 3, 3, 2, 3, 1, 1, 3, 3, 1, 2, 1, 2, 0, 0, 3 },
80 { 0, 1, 0, 1, 0, 2, 0, 2, 2, 0, 3, 0, 1, 3, 2, 1 },
81 { 3, 1, 1, 2, 0, 1, 0, 2, 1, 3, 1, 3, 0, 1, 3, 0 },
82 { 2, 0, 0, 3, 3, 2, 3, 1, 2, 0, 2, 0, 3, 2, 2, 1 },
83 { 2, 3, 3, 1, 2, 1, 2, 1, 2, 1, 1, 2, 3, 0, 0, 1 },
84 { 1, 0, 0, 2, 3, 0, 0, 3, 0, 3, 0, 3, 2, 1, 2, 3 },
85 { 2, 3, 3, 1, 1, 2, 1, 0, 3, 2, 3, 0, 2, 3, 1, 3 },
86 { 1, 0, 2, 0, 3, 0, 3, 2, 0, 1, 1, 2, 0, 1, 0, 2 },
87 { 0, 1, 1, 3, 3, 2, 2, 1, 1, 3, 3, 0, 2, 1, 3, 2 },
88 { 2, 3, 2, 0, 0, 1, 3, 0, 2, 0, 1, 2, 3, 0, 1, 0 },
89 { 1, 3, 1, 2, 3, 2, 3, 2, 0, 2, 0, 1, 1, 0, 3, 0 },
90 { 0, 2, 0, 3, 1, 0, 0, 1, 1, 3, 3, 2, 3, 2, 2, 1 },
91 { 2, 1, 3, 2, 3, 1, 2, 1, 0, 3, 0, 2, 0, 2, 0, 2 },
92 { 0, 3, 1, 0, 0, 2, 0, 3, 2, 1, 3, 1, 1, 3, 1, 3 }
93 };
94
95 if (filters == 1) return filter[(row + top_margin) & 15][(col + left_margin) & 15];
96 if (filters == 9) return xtrans[(row + 6) % 6][(col + 6) % 6];
97 return FC(row, col);
98 }
99
merror(void * ptr,char * where)100 static void CLASS merror(void *ptr, char *where)
101 {
102 if (ptr) return;
103 g_error("Out of memory in %s\n", where);
104 }
105
hat_transform(float * temp,float * base,int st,int size,int sc)106 static void CLASS hat_transform(float *temp, float *base, int st, int size, int sc)
107 {
108 int i;
109 for (i = 0; i < sc; i++)
110 temp[i] = 2 * base[st * i] + base[st * (sc - i)] + base[st * (i + sc)];
111 for (; i + sc < size; i++)
112 temp[i] = 2 * base[st * i] + base[st * (i - sc)] + base[st * (i + sc)];
113 for (; i < size; i++)
114 temp[i] = 2 * base[st * i] + base[st * (i - sc)] + base[st * (2 * size - 2 - (i + sc))];
115 }
116
wavelet_denoise_INDI(ushort (* image)[4],const int black,const int iheight,const int iwidth,const int height,const int width,const int colors,const int shrink,const float pre_mul[4],const float threshold,const unsigned filters)117 void CLASS wavelet_denoise_INDI(ushort(*image)[4], const int black,
118 const int iheight, const int iwidth,
119 const int height, const int width,
120 const int colors, const int shrink,
121 const float pre_mul[4], const float threshold,
122 const unsigned filters)
123 {
124 float *fimg = 0, thold, mul[2], avg, diff;
125 int size, lev, hpass, lpass, row, col, nc, c, i, wlast;
126 ushort *window[4];
127 static const float noise[] =
128 { 0.8002, 0.2735, 0.1202, 0.0585, 0.0291, 0.0152, 0.0080, 0.0044 };
129
130 // dcraw_message (dcraw, DCRAW_VERBOSE,_("Wavelet denoising...\n")); /*UF*/
131
132 /* Scaling is done somewhere else - NKBJ*/
133 size = iheight * iwidth;
134 float temp[iheight + iwidth];
135 if ((nc = colors) == 3 && filters) nc++;
136 progress(PROGRESS_WAVELET_DENOISE, -nc * 5);
137 #ifdef _OPENMP
138 #if defined(__sun) && !defined(__GNUC__) /* Fix bug #3205673 - NKBJ */
139 #pragma omp parallel for \
140 shared(nc,image,size,noise) \
141 private(c,i,hpass,lev,lpass,row,col,thold,fimg,temp)
142 #else
143 #pragma omp parallel for \
144 shared(nc,image,size) \
145 private(c,i,hpass,lev,lpass,row,col,thold,fimg,temp)
146 #endif
147 #endif
148 FORC(nc) { /* denoise R,G1,B,G3 individually */
149 fimg = (float *) malloc(size * 3 * sizeof * fimg);
150 for (i = 0; i < size; i++)
151 fimg[i] = 256 * sqrt(image[i][c] /*<< scale*/);
152 for (hpass = lev = 0; lev < 5; lev++) {
153 progress(PROGRESS_WAVELET_DENOISE, 1);
154 lpass = size * ((lev & 1) + 1);
155 for (row = 0; row < iheight; row++) {
156 hat_transform(temp, fimg + hpass + row * iwidth, 1, iwidth, 1 << lev);
157 for (col = 0; col < iwidth; col++)
158 fimg[lpass + row * iwidth + col] = temp[col] * 0.25;
159 }
160 for (col = 0; col < iwidth; col++) {
161 hat_transform(temp, fimg + lpass + col, iwidth, iheight, 1 << lev);
162 for (row = 0; row < iheight; row++)
163 fimg[lpass + row * iwidth + col] = temp[row] * 0.25;
164 }
165 thold = threshold * noise[lev];
166 for (i = 0; i < size; i++) {
167 fimg[hpass + i] -= fimg[lpass + i];
168 if (fimg[hpass + i] < -thold) fimg[hpass + i] += thold;
169 else if (fimg[hpass + i] > thold) fimg[hpass + i] -= thold;
170 else fimg[hpass + i] = 0;
171 if (hpass) fimg[i] += fimg[hpass + i];
172 }
173 hpass = lpass;
174 }
175 for (i = 0; i < size; i++)
176 image[i][c] = CLIP(SQR(fimg[i] + fimg[lpass + i]) / 0x10000);
177 free(fimg);
178 }
179 if (filters && colors == 3) { /* pull G1 and G3 closer together */
180 for (row = 0; row < 2; row++)
181 mul[row] = 0.125 * pre_mul[FC(row + 1, 0) | 1] / pre_mul[FC(row, 0) | 1];
182 ushort window_mem[4][width];
183 for (i = 0; i < 4; i++)
184 window[i] = window_mem[i]; /*(ushort *) fimg + width*i;*/
185 for (wlast = -1, row = 1; row < height - 1; row++) {
186 while (wlast < row + 1) {
187 for (wlast++, i = 0; i < 4; i++)
188 window[(i + 3) & 3] = window[i];
189 for (col = FC(wlast, 1) & 1; col < width; col += 2)
190 window[2][col] = BAYER(wlast, col);
191 }
192 thold = threshold / 512;
193 for (col = (FC(row, 0) & 1) + 1; col < width - 1; col += 2) {
194 avg = (window[0][col - 1] + window[0][col + 1] +
195 window[2][col - 1] + window[2][col + 1] - black * 4)
196 * mul[row & 1] + (window[1][col] - black) * 0.5 + black;
197 avg = avg < 0 ? 0 : sqrt(avg);
198 diff = sqrt(BAYER(row, col)) - avg;
199 if (diff < -thold) diff += thold;
200 else if (diff > thold) diff -= thold;
201 else diff = 0;
202 BAYER(row, col) = CLIP(SQR(avg + diff) + 0.5);
203 }
204 }
205 }
206 }
207
scale_colors_INDI(const int maximum,const int black,const int use_camera_wb,const float cam_mul[4],const int colors,float pre_mul[4],const unsigned filters,ushort white[8][8],const char * ifname_display,void * dcraw)208 void CLASS scale_colors_INDI(const int maximum, const int black,
209 const int use_camera_wb, const float cam_mul[4], const int colors,
210 float pre_mul[4], const unsigned filters, /*const*/ ushort white[8][8],
211 const char *ifname_display, void *dcraw)
212 {
213 unsigned row, col, c, sum[8];
214 int val;
215 double dmin, dmax;
216
217 if (use_camera_wb && cam_mul[0] != -1) {
218 memset(sum, 0, sizeof sum);
219 for (row = 0; row < 8; row++)
220 for (col = 0; col < 8; col++) {
221 c = FC(row, col);
222 if ((val = white[row][col] - black) > 0)
223 sum[c] += val;
224 sum[c + 4]++;
225 }
226 if (sum[0] && sum[1] && sum[2] && sum[3])
227 FORC4 pre_mul[c] = (float) sum[c + 4] / sum[c];
228 else if (cam_mul[0] && cam_mul[2])
229 /* 'sizeof pre_mul' does not work because pre_mul is an argument (UF)*/
230 memcpy(pre_mul, cam_mul, 4 * sizeof(float));
231 else
232 dcraw_message(dcraw, DCRAW_NO_CAMERA_WB,
233 _("%s: Cannot use camera white balance.\n"), ifname_display);
234 } else {
235 dcraw_message(dcraw, DCRAW_NO_CAMERA_WB,
236 _("%s: Cannot use camera white balance.\n"), ifname_display);
237 }
238 if (pre_mul[1] == 0) pre_mul[1] = 1;
239 if (pre_mul[3] == 0) pre_mul[3] = colors < 4 ? pre_mul[1] : 1;
240 for (dmin = DBL_MAX, dmax = c = 0; c < 4; c++) {
241 if (dmin > pre_mul[c])
242 dmin = pre_mul[c];
243 if (dmax < pre_mul[c])
244 dmax = pre_mul[c];
245 }
246 FORC4 pre_mul[c] /= dmax;
247 dcraw_message(dcraw, DCRAW_VERBOSE,
248 _("Scaling with darkness %d, saturation %d, and\nmultipliers"),
249 black, maximum);
250 FORC4 dcraw_message(dcraw, DCRAW_VERBOSE, " %f", pre_mul[c]);
251 dcraw_message(dcraw, DCRAW_VERBOSE, "\n");
252
253 /* The rest of the scaling is done somewhere else UF*/
254 }
255
border_interpolate_INDI(const int height,const int width,ushort (* image)[4],const unsigned filters,int colors,int border,dcraw_data * h)256 void CLASS border_interpolate_INDI(const int height, const int width,
257 ushort(*image)[4], const unsigned filters, int colors, int border, dcraw_data *h)
258 {
259 int row, col, y, x, f, c, sum[8];
260
261 for (row = 0; row < height; row++)
262 for (col = 0; col < width; col++) {
263 if (col == border && row >= border && row < height - border)
264 col = width - border;
265 memset(sum, 0, sizeof sum);
266 for (y = row - 1; y != row + 2; y++)
267 for (x = col - 1; x != col + 2; x++)
268 if (y >= 0 && y < height && x >= 0 && x < width) {
269 f = fcol_INDI(filters, y, x, h->top_margin, h->left_margin, h->xtrans);
270 sum[f] += image[y * width + x][f];
271 sum[f + 4]++;
272 }
273 f = fcol_INDI(filters, row, col, h->top_margin, h->left_margin, h->xtrans);
274 FORCC if (c != f && sum[c + 4])
275 image[row * width + col][c] = sum[c] / sum[c + 4];
276 }
277 }
278
lin_interpolate_INDI(ushort (* image)[4],const unsigned filters,const int width,const int height,const int colors,void * dcraw,dcraw_data * h)279 void CLASS lin_interpolate_INDI(ushort(*image)[4], const unsigned filters,
280 const int width, const int height, const int colors, void *dcraw, dcraw_data *h) /*UF*/
281 {
282 int code[16][16][32], size = 16, *ip, sum[4];
283 int f, c, i, x, y, row, col, shift, color;
284 ushort *pix;
285
286 dcraw_message(dcraw, DCRAW_VERBOSE, _("Bilinear interpolation...\n")); /*UF*/
287 if (filters == 9) size = 6;
288 border_interpolate_INDI(height, width, image, filters, colors, 1, h);
289 for (row = 0; row < size; row++) {
290 for (col = 0; col < size; col++) {
291 ip = code[row][col] + 1;
292 f = fcol_INDI(filters, row, col, h->top_margin, h->left_margin, h->xtrans);
293 memset(sum, 0, sizeof sum);
294 for (y = -1; y <= 1; y++)
295 for (x = -1; x <= 1; x++) {
296 shift = (y == 0) + (x == 0);
297 color = fcol_INDI(filters, row + y, col + x, h->top_margin, h->left_margin, h->xtrans);
298 if (color == f) continue;
299 *ip++ = (width * y + x) * 4 + color;
300 *ip++ = shift;
301 *ip++ = color;
302 sum[color] += 1 << shift;
303 }
304 code[row][col][0] = (ip - code[row][col]) / 3;
305 FORCC
306 if (c != f) {
307 *ip++ = c;
308 *ip++ = 256 / sum[c];
309 }
310 }
311 }
312 #ifdef _OPENMP
313 #pragma omp parallel for default(shared) private(row,col,pix,ip,sum,i)
314 #endif
315 for (row = 1; row < height - 1; row++) {
316 for (col = 1; col < width - 1; col++) {
317 pix = image[row * width + col];
318 ip = code[row % size][col % size];
319 memset(sum, 0, sizeof sum);
320 for (i = *ip++; i--; ip += 3)
321 sum[ip[2]] += pix[ip[0]] << ip[1];
322 for (i = colors; --i; ip += 2)
323 pix[ip[0]] = sum[ip[0]] * ip[1] >> 8;
324 }
325 }
326 }
327
328 /*
329 This algorithm is officially called:
330
331 "Interpolation using a Threshold-based variable number of gradients"
332
333 described in http://scien.stanford.edu/class/psych221/projects/99/tingchen/algodep/vargra.html
334
335 I've extended the basic idea to work with non-Bayer filter arrays.
336 Gradients are numbered clockwise from NW=0 to W=7.
337 */
vng_interpolate_INDI(ushort (* image)[4],const unsigned filters,const int width,const int height,const int colors,void * dcraw,dcraw_data * h)338 void CLASS vng_interpolate_INDI(ushort(*image)[4], const unsigned filters,
339 const int width, const int height, const int colors, void *dcraw, dcraw_data *h) /*UF*/
340 {
341 static const signed char *cp, terms[] = {
342 -2, -2, +0, -1, 0, 0x01, -2, -2, +0, +0, 1, 0x01, -2, -1, -1, +0, 0, 0x01,
343 -2, -1, +0, -1, 0, 0x02, -2, -1, +0, +0, 0, 0x03, -2, -1, +0, +1, 1, 0x01,
344 -2, +0, +0, -1, 0, 0x06, -2, +0, +0, +0, 1, 0x02, -2, +0, +0, +1, 0, 0x03,
345 -2, +1, -1, +0, 0, 0x04, -2, +1, +0, -1, 1, 0x04, -2, +1, +0, +0, 0, 0x06,
346 -2, +1, +0, +1, 0, 0x02, -2, +2, +0, +0, 1, 0x04, -2, +2, +0, +1, 0, 0x04,
347 -1, -2, -1, +0, 0, 0x80, -1, -2, +0, -1, 0, 0x01, -1, -2, +1, -1, 0, 0x01,
348 -1, -2, +1, +0, 1, 0x01, -1, -1, -1, +1, 0, 0x88, -1, -1, +1, -2, 0, 0x40,
349 -1, -1, +1, -1, 0, 0x22, -1, -1, +1, +0, 0, 0x33, -1, -1, +1, +1, 1, 0x11,
350 -1, +0, -1, +2, 0, 0x08, -1, +0, +0, -1, 0, 0x44, -1, +0, +0, +1, 0, 0x11,
351 -1, +0, +1, -2, 1, 0x40, -1, +0, +1, -1, 0, 0x66, -1, +0, +1, +0, 1, 0x22,
352 -1, +0, +1, +1, 0, 0x33, -1, +0, +1, +2, 1, 0x10, -1, +1, +1, -1, 1, 0x44,
353 -1, +1, +1, +0, 0, 0x66, -1, +1, +1, +1, 0, 0x22, -1, +1, +1, +2, 0, 0x10,
354 -1, +2, +0, +1, 0, 0x04, -1, +2, +1, +0, 1, 0x04, -1, +2, +1, +1, 0, 0x04,
355 +0, -2, +0, +0, 1, 0x80, +0, -1, +0, +1, 1, 0x88, +0, -1, +1, -2, 0, 0x40,
356 +0, -1, +1, +0, 0, 0x11, +0, -1, +2, -2, 0, 0x40, +0, -1, +2, -1, 0, 0x20,
357 +0, -1, +2, +0, 0, 0x30, +0, -1, +2, +1, 1, 0x10, +0, +0, +0, +2, 1, 0x08,
358 +0, +0, +2, -2, 1, 0x40, +0, +0, +2, -1, 0, 0x60, +0, +0, +2, +0, 1, 0x20,
359 +0, +0, +2, +1, 0, 0x30, +0, +0, +2, +2, 1, 0x10, +0, +1, +1, +0, 0, 0x44,
360 +0, +1, +1, +2, 0, 0x10, +0, +1, +2, -1, 1, 0x40, +0, +1, +2, +0, 0, 0x60,
361 +0, +1, +2, +1, 0, 0x20, +0, +1, +2, +2, 0, 0x10, +1, -2, +1, +0, 0, 0x80,
362 +1, -1, +1, +1, 0, 0x88, +1, +0, +1, +2, 0, 0x08, +1, +0, +2, -1, 0, 0x40,
363 +1, +0, +2, +1, 0, 0x10
364 }, chood[] = { -1, -1, -1, 0, -1, +1, 0, +1, +1, +1, +1, 0, +1, -1, 0, -1 };
365 ushort(*brow[4])[4], *pix;
366 int prow = 8, pcol = 2, *ip, *code[16][16], gval[8], gmin, gmax, sum[4];
367 int row, col, x, y, x1, x2, y1, y2, t, weight, grads, color, diag;
368 int g, diff, thold, num, c;
369 ushort rowtmp[4][width * 4];
370
371 lin_interpolate_INDI(image, filters, width, height, colors, dcraw, h); /*UF*/
372 dcraw_message(dcraw, DCRAW_VERBOSE, _("VNG interpolation...\n")); /*UF*/
373
374 if (filters == 1) prow = pcol = 16;
375 if (filters == 9) prow = pcol = 6;
376 int *ipalloc = ip = (int *) calloc(prow * pcol, 1280);
377 merror(ip, "vng_interpolate()");
378 for (row = 0; row < prow; row++) /* Precalculate for VNG */
379 for (col = 0; col < pcol; col++) {
380 code[row][col] = ip;
381 for (cp = terms, t = 0; t < 64; t++) {
382 y1 = *cp++;
383 x1 = *cp++;
384 y2 = *cp++;
385 x2 = *cp++;
386 weight = *cp++;
387 grads = *cp++;
388 color = fcol_INDI(filters, row + y1, col + x1, h->top_margin, h->left_margin, h->xtrans);
389 if (fcol_INDI(filters, row + y2, col + x2, h->top_margin, h->left_margin, h->xtrans) != color) continue;
390 diag = (fcol_INDI(filters, row, col + 1, h->top_margin, h->left_margin, h->xtrans) == color && fcol_INDI(filters, row + 1, col, h->top_margin, h->left_margin, h->xtrans) == color) ? 2 : 1;
391 if (abs(y1 - y2) == diag && abs(x1 - x2) == diag) continue;
392 *ip++ = (y1 * width + x1) * 4 + color;
393 *ip++ = (y2 * width + x2) * 4 + color;
394 *ip++ = weight;
395 for (g = 0; g < 8; g++)
396 if (grads & 1 << g) *ip++ = g;
397 *ip++ = -1;
398 }
399 *ip++ = INT_MAX;
400 for (cp = chood, g = 0; g < 8; g++) {
401 y = *cp++;
402 x = *cp++;
403 *ip++ = (y * width + x) * 4;
404 color = fcol_INDI(filters, row, col, h->top_margin, h->left_margin, h->xtrans);
405 if (fcol_INDI(filters, row + y, col + x, h->top_margin, h->left_margin, h->xtrans) != color && fcol_INDI(filters, row + y * 2, col + x * 2, h->top_margin, h->left_margin, h->xtrans) == color)
406 *ip++ = (y * width + x) * 8 + color;
407 else
408 *ip++ = 0;
409 }
410 }
411 progress(PROGRESS_INTERPOLATE, -height);
412 #ifdef _OPENMP
413 #pragma omp parallel \
414 shared(image,code,prow,pcol,h) \
415 private(row,col,g,brow,rowtmp,pix,ip,gval,diff,gmin,gmax,thold,sum,color,num,c,t)
416 #endif
417 {
418 int slice = (height - 4) / uf_omp_get_num_threads();
419 int start_row = 2 + slice * uf_omp_get_thread_num();
420 int end_row = MIN(start_row + slice, height - 2);
421 for (row = start_row; row < end_row; row++) { /* Do VNG interpolation */
422 progress(PROGRESS_INTERPOLATE, 1);
423 for (g = 0; g < 4; g++)
424 brow[g] = &rowtmp[(row + g - 2) % 4];
425 for (col = 2; col < width - 2; col++) {
426 pix = image[row * width + col];
427 ip = code[row % prow][col % pcol];
428 memset(gval, 0, sizeof gval);
429 while ((g = ip[0]) != INT_MAX) { /* Calculate gradients */
430 diff = ABS(pix[g] - pix[ip[1]]) << ip[2];
431 gval[ip[3]] += diff;
432 ip += 5;
433 if ((g = ip[-1]) == -1) continue;
434 gval[g] += diff;
435 while ((g = *ip++) != -1)
436 gval[g] += diff;
437 }
438 ip++;
439 gmin = gmax = gval[0]; /* Choose a threshold */
440 for (g = 1; g < 8; g++) {
441 if (gmin > gval[g]) gmin = gval[g];
442 if (gmax < gval[g]) gmax = gval[g];
443 }
444 if (gmax == 0) {
445 memcpy(brow[2][col], pix, sizeof * image);
446 continue;
447 }
448 thold = gmin + (gmax >> 1);
449 memset(sum, 0, sizeof sum);
450 color = fcol_INDI(filters, row, col, h->top_margin, h->left_margin, h->xtrans);
451 for (num = g = 0; g < 8; g++, ip += 2) { /* Average the neighbors */
452 if (gval[g] <= thold) {
453 FORCC
454 if (c == color && ip[1])
455 sum[c] += (pix[c] + pix[ip[1]]) >> 1;
456 else
457 sum[c] += pix[ip[0] + c];
458 num++;
459 }
460 }
461 FORCC { /* Save to buffer */
462 t = pix[color];
463 if (c != color)
464 t += (sum[c] - sum[color]) / num;
465 brow[2][col][c] = CLIP(t);
466 }
467 }
468 /* Write buffer to image */
469 if ((row > start_row + 1) || (row == height - 2))
470 memcpy(image[(row - 2)*width + 2], brow[0] + 2, (width - 4)*sizeof * image);
471 if (row == height - 2) {
472 memcpy(image[(row - 1)*width + 2], brow[1] + 2, (width - 4)*sizeof * image);
473 break;
474 }
475 }
476 }
477 free(ipalloc);
478 }
479
480 /*
481 Patterned Pixel Grouping Interpolation by Alain Desbiolles
482 */
ppg_interpolate_INDI(ushort (* image)[4],const unsigned filters,const int width,const int height,const int colors,void * dcraw,dcraw_data * h)483 void CLASS ppg_interpolate_INDI(ushort(*image)[4], const unsigned filters,
484 const int width, const int height,
485 const int colors, void *dcraw, dcraw_data *h)
486 {
487 int dir[5] = { 1, width, -1, -width, 1 };
488 int row, col, diff[2] = { 0, 0 }, guess[2], c, d, i;
489 ushort(*pix)[4];
490
491 border_interpolate_INDI(height, width, image, filters, colors, 3, h);
492 dcraw_message(dcraw, DCRAW_VERBOSE, _("PPG interpolation...\n")); /*UF*/
493
494 #ifdef _OPENMP
495 #pragma omp parallel \
496 shared(image,dir,diff) \
497 private(row,col,i,d,c,pix,guess)
498 #endif
499 {
500 /* Fill in the green layer with gradients and pattern recognition: */
501 #ifdef _OPENMP
502 #pragma omp for
503 #endif
504 for (row = 3; row < height - 3; row++) {
505 for (col = 3 + (FC(row, 3) & 1), c = FC(row, col); col < width - 3; col += 2) {
506 pix = image + row * width + col;
507 for (i = 0; (d = dir[i]) > 0; i++) {
508 guess[i] = (pix[-d][1] + pix[0][c] + pix[d][1]) * 2
509 - pix[-2 * d][c] - pix[2 * d][c];
510 diff[i] = (ABS(pix[-2 * d][c] - pix[ 0][c]) +
511 ABS(pix[ 2 * d][c] - pix[ 0][c]) +
512 ABS(pix[ -d][1] - pix[ d][1])) * 3 +
513 (ABS(pix[ 3 * d][1] - pix[ d][1]) +
514 ABS(pix[-3 * d][1] - pix[-d][1])) * 2;
515 }
516 d = dir[i = diff[0] > diff[1]];
517 pix[0][1] = ULIM(guess[i] >> 2, pix[d][1], pix[-d][1]);
518 }
519 }
520 /* Calculate red and blue for each green pixel: */
521 #ifdef _OPENMP
522 #pragma omp for
523 #endif
524 for (row = 1; row < height - 1; row++) {
525 for (col = 1 + (FC(row, 2) & 1), c = FC(row, col + 1); col < width - 1; col += 2) {
526 pix = image + row * width + col;
527 for (i = 0; (d = dir[i]) > 0; c = 2 - c, i++)
528 pix[0][c] = CLIP((pix[-d][c] + pix[d][c] + 2 * pix[0][1]
529 - pix[-d][1] - pix[d][1]) >> 1);
530 }
531 }
532 /* Calculate blue for red pixels and vice versa: */
533 #ifdef _OPENMP
534 #pragma omp for
535 #endif
536 for (row = 1; row < height - 1; row++) {
537 for (col = 1 + (FC(row, 1) & 1), c = 2 - FC(row, col); col < width - 1; col += 2) {
538 pix = image + row * width + col;
539 for (i = 0; (d = dir[i] + dir[i + 1]) > 0; i++) {
540 diff[i] = ABS(pix[-d][c] - pix[d][c]) +
541 ABS(pix[-d][1] - pix[0][1]) +
542 ABS(pix[ d][1] - pix[0][1]);
543 guess[i] = pix[-d][c] + pix[d][c] + 2 * pix[0][1]
544 - pix[-d][1] - pix[d][1];
545 }
546 if (diff[0] != diff[1])
547 pix[0][c] = CLIP(guess[diff[0] > diff[1]] >> 1);
548 else
549 pix[0][c] = CLIP((guess[0] + guess[1]) >> 2);
550 }
551 }
552 }
553 }
554
cielab_INDI(ushort rgb[3],short lab[3],const int colors,const float rgb_cam[3][4])555 void CLASS cielab_INDI(ushort rgb[3], short lab[3], const int colors,
556 const float rgb_cam[3][4])
557 {
558 int c, i, j, k;
559 float r, xyz[3];
560 static float cbrt[0x10000], xyz_cam[3][4];
561
562 if (!rgb) {
563 for (i = 0; i < 0x10000; i++) {
564 r = i / 65535.0;
565 cbrt[i] = r > 0.008856 ? pow(r, (float)(1 / 3.0)) : 7.787 * r + 16 / 116.0;
566 }
567 for (i = 0; i < 3; i++)
568 for (j = 0; j < colors; j++)
569 for (xyz_cam[i][j] = k = 0; k < 3; k++)
570 xyz_cam[i][j] += xyz_rgb[i][k] * rgb_cam[k][j] / d65_white[i];
571 return;
572 }
573 xyz[0] = xyz[1] = xyz[2] = 0.5;
574 FORCC {
575 xyz[0] += xyz_cam[0][c] * rgb[c];
576 xyz[1] += xyz_cam[1][c] * rgb[c];
577 xyz[2] += xyz_cam[2][c] * rgb[c];
578 }
579 xyz[0] = cbrt[CLIP((int) xyz[0])];
580 xyz[1] = cbrt[CLIP((int) xyz[1])];
581 xyz[2] = cbrt[CLIP((int) xyz[2])];
582 lab[0] = 64 * (116 * xyz[1] - 16);
583 lab[1] = 64 * 500 * (xyz[0] - xyz[1]);
584 lab[2] = 64 * 200 * (xyz[1] - xyz[2]);
585 }
586
587 #define TS 512 /* Tile Size */
588 /*
589 Frank Markesteijn's algorithm for Fuji X-Trans sensors
590 */
xtrans_interpolate_INDI(ushort (* image)[4],const unsigned filters,const int width,const int height,const int colors,const float rgb_cam[3][4],void * dcraw,dcraw_data * hh,const int passes)591 void CLASS xtrans_interpolate_INDI(ushort(*image)[4], const unsigned filters,
592 const int width, const int height,
593 const int colors, const float rgb_cam[3][4],
594 void *dcraw, dcraw_data *hh, const int passes)
595 {
596 int c, d, f, g, h, i, v, ng, row, col, top, left, mrow, mcol;
597 int val, ndir, pass, hm[8], avg[4], color[3][8];
598 static const short orth[12] = { 1, 0, 0, 1, -1, 0, 0, -1, 1, 0, 0, 1 },
599 patt[2][16] = { { 0, 1, 0, -1, 2, 0, -1, 0, 1, 1, 1, -1, 0, 0, 0, 0 },
600 { 0, 1, 0, -2, 1, 0, -2, 0, 1, 1, -2, -2, 1, -1, -1, 1 }
601 },
602 dir[4] = { 1, TS, TS + 1, TS - 1 };
603 short allhex[3][3][2][8], *hex;
604 ushort min, max, sgrow = 0, sgcol = 0;
605 ushort(*rgb)[TS][TS][3], (*rix)[3], (*pix)[4];
606 short(*lab) [TS][3], (*lix)[3];
607 float(*drv)[TS][TS], diff[6], tr;
608 char(*homo)[TS][TS], *buffer;
609
610 dcraw_message(dcraw, DCRAW_VERBOSE, _("%d-pass X-Trans interpolation...\n"), passes); /*NKBJ*/
611
612 cielab_INDI(0, 0, colors, rgb_cam);
613 ndir = 4 << (passes > 1);
614
615 /* Map a green hexagon around each non-green pixel and vice versa: */
616 for (row = 0; row < 3; row++)
617 for (col = 0; col < 3; col++)
618 for (ng = d = 0; d < 10; d += 2) {
619 g = fcol_INDI(filters, row, col, hh->top_margin, hh->left_margin, hh->xtrans) == 1;
620 if (fcol_INDI(filters, row + orth[d], col + orth[d + 2], hh->top_margin, hh->left_margin, hh->xtrans) == 1) ng = 0;
621 else ng++;
622 if (ng == 4) {
623 sgrow = row;
624 sgcol = col;
625 }
626 if (ng == g + 1) FORC(8) {
627 v = orth[d ] * patt[g][c * 2] + orth[d + 1] * patt[g][c * 2 + 1];
628 h = orth[d + 2] * patt[g][c * 2] + orth[d + 3] * patt[g][c * 2 + 1];
629 allhex[row][col][0][c ^ (g * 2 & d)] = h + v * width;
630 allhex[row][col][1][c ^ (g * 2 & d)] = h + v * TS;
631 }
632 }
633
634 /* Set green1 and green3 to the minimum and maximum allowed values: */
635 for (row = 2; row < height - 2; row++)
636 for (min = ~(max = 0), col = 2; col < width - 2; col++) {
637 if (fcol_INDI(filters, row, col, hh->top_margin, hh->left_margin, hh->xtrans) == 1 && (min = ~(max = 0))) continue;
638 pix = image + row * width + col;
639 hex = allhex[row % 3][col % 3][0];
640 if (!max) FORC(6) {
641 val = pix[hex[c]][1];
642 if (min > val) min = val;
643 if (max < val) max = val;
644 }
645 pix[0][1] = min;
646 pix[0][3] = max;
647 switch ((row - sgrow) % 3) {
648 case 1:
649 if (row < height - 3) {
650 row++;
651 col--;
652 }
653 break;
654 case 2:
655 if ((min = ~(max = 0)) && (col += 2) < width - 3 && row > 2) row--;
656 }
657 }
658
659
660 #ifdef _OPENMP
661 #pragma omp parallel \
662 default(shared) \
663 private(top, left, row, col, pix, mrow, mcol, hex, color, c, pass, rix, val, d, f, g, h, i, diff, lix, tr, avg, v, buffer, rgb, lab, drv, homo, hm, max)
664 #endif
665 {
666 buffer = (char *) malloc(TS * TS * (ndir * 11 + 6));
667 merror(buffer, "xtrans_interpolate()");
668 rgb = (ushort(*)[TS][TS][3]) buffer;
669 lab = (short(*) [TS][3])(buffer + TS * TS * (ndir * 6));
670 drv = (float(*)[TS][TS])(buffer + TS * TS * (ndir * 6 + 6));
671 homo = (char(*)[TS][TS])(buffer + TS * TS * (ndir * 10 + 6));
672
673 progress(PROGRESS_INTERPOLATE, -height);
674
675 #ifdef _OPENMP
676 #pragma omp for
677 #endif
678
679 for (top = 3; top < height - 19; top += TS - 16) {
680 progress(PROGRESS_INTERPOLATE, TS - 16);
681 for (left = 3; left < width - 19; left += TS - 16) {
682 mrow = MIN(top + TS, height - 3);
683 mcol = MIN(left + TS, width - 3);
684 for (row = top; row < mrow; row++)
685 for (col = left; col < mcol; col++)
686 memcpy(rgb[0][row - top][col - left], image[row * width + col], 6);
687 FORC3 memcpy(rgb[c + 1], rgb[0], sizeof * rgb);
688
689 /* Interpolate green horizontally, vertically, and along both diagonals: */
690 for (row = top; row < mrow; row++)
691 for (col = left; col < mcol; col++) {
692 if ((f = fcol_INDI(filters, row, col, hh->top_margin, hh->left_margin, hh->xtrans)) == 1) continue;
693 pix = image + row * width + col;
694 hex = allhex[row % 3][col % 3][0];
695 color[1][0] = 174 * (pix[ hex[1]][1] + pix[ hex[0]][1]) -
696 46 * (pix[2 * hex[1]][1] + pix[2 * hex[0]][1]);
697 color[1][1] = 223 * pix[ hex[3]][1] + pix[ hex[2]][1] * 33 +
698 92 * (pix[ 0 ][f] - pix[ -hex[2]][f]);
699 FORC(2) color[1][2 + c] =
700 164 * pix[hex[4 + c]][1] + 92 * pix[-2 * hex[4 + c]][1] + 33 *
701 (2 * pix[0][f] - pix[3 * hex[4 + c]][f] - pix[-3 * hex[4 + c]][f]);
702 FORC4 rgb[c ^ !((row - sgrow) % 3)][row - top][col - left][1] =
703 LIM(color[1][c] >> 8, pix[0][1], pix[0][3]);
704 }
705
706 for (pass = 0; pass < passes; pass++) {
707 if (pass == 1)
708 memcpy(rgb += 4, buffer, 4 * sizeof * rgb);
709
710 /* Recalculate green from interpolated values of closer pixels: */
711 if (pass) {
712 for (row = top + 2; row < mrow - 2; row++)
713 for (col = left + 2; col < mcol - 2; col++) {
714 if ((f = fcol_INDI(filters, row, col, hh->top_margin, hh->left_margin, hh->xtrans)) == 1) continue;
715 pix = image + row * width + col;
716 hex = allhex[row % 3][col % 3][1];
717 for (d = 3; d < 6; d++) {
718 rix = &rgb[(d - 2) ^ !((row - sgrow) % 3)][row - top][col - left];
719 val = rix[-2 * hex[d]][1] + 2 * rix[hex[d]][1]
720 - rix[-2 * hex[d]][f] - 2 * rix[hex[d]][f] + 3 * rix[0][f];
721 rix[0][1] = LIM(val / 3, pix[0][1], pix[0][3]);
722 }
723 }
724 }
725
726 /* Interpolate red and blue values for solitary green pixels: */
727 for (row = (top - sgrow + 4) / 3 * 3 + sgrow; row < mrow - 2; row += 3)
728 for (col = (left - sgcol + 4) / 3 * 3 + sgcol; col < mcol - 2; col += 3) {
729 rix = &rgb[0][row - top][col - left];
730 h = fcol_INDI(filters, row, col + 1, hh->top_margin, hh->left_margin, hh->xtrans);
731 memset(diff, 0, sizeof diff);
732 for (i = 1, d = 0; d < 6; d++, i ^= TS ^ 1, h ^= 2) {
733 for (c = 0; c < 2; c++, h ^= 2) {
734 g = 2 * rix[0][1] - rix[i << c][1] - rix[-i << c][1];
735 color[h][d] = g + rix[i << c][h] + rix[-i << c][h];
736 if (d > 1)
737 diff[d] += SQR(rix[i << c][1] - rix[-i << c][1]
738 - rix[i << c][h] + rix[-i << c][h]) + SQR(g);
739 }
740 if (d > 1 && (d & 1))
741 if (diff[d - 1] < diff[d])
742 FORC(2) color[c * 2][d] = color[c * 2][d - 1];
743 if (d < 2 || (d & 1)) {
744 FORC(2) rix[0][c * 2] = CLIP(color[c * 2][d] / 2);
745 rix += TS * TS;
746 }
747 }
748 }
749
750 /* Interpolate red for blue pixels and vice versa: */
751 for (row = top + 3; row < mrow - 3; row++)
752 for (col = left + 3; col < mcol - 3; col++) {
753 if ((f = 2 - fcol_INDI(filters, row, col, hh->top_margin, hh->left_margin, hh->xtrans)) == 1) continue;
754 rix = &rgb[0][row - top][col - left];
755 c = (row - sgrow) % 3 ? TS : 1;
756 h = 3 * (c ^ TS ^ 1);
757 for (d = 0; d < 4; d++, rix += TS * TS) {
758 i = d > 1 || ((d ^ c) & 1) ||
759 ((ABS(rix[0][1] - rix[c][1]) + ABS(rix[0][1] - rix[-c][1])) <
760 2 * (ABS(rix[0][1] - rix[h][1]) + ABS(rix[0][1] - rix[-h][1]))) ? c : h;
761 rix[0][f] = CLIP((rix[i][f] + rix[-i][f] +
762 2 * rix[0][1] - rix[i][1] - rix[-i][1]) / 2);
763 }
764 }
765
766 /* Fill in red and blue for 2x2 blocks of green: */
767 for (row = top + 2; row < mrow - 2; row++) if ((row - sgrow) % 3)
768 for (col = left + 2; col < mcol - 2; col++) if ((col - sgcol) % 3) {
769 rix = &rgb[0][row - top][col - left];
770 hex = allhex[row % 3][col % 3][1];
771 for (d = 0; d < ndir; d += 2, rix += TS * TS)
772 if (hex[d] + hex[d + 1]) {
773 g = 3 * rix[0][1] - 2 * rix[hex[d]][1] - rix[hex[d + 1]][1];
774 for (c = 0; c < 4; c += 2) rix[0][c] =
775 CLIP((g + 2 * rix[hex[d]][c] + rix[hex[d + 1]][c]) / 3);
776 } else {
777 g = 2 * rix[0][1] - rix[hex[d]][1] - rix[hex[d + 1]][1];
778 for (c = 0; c < 4; c += 2) rix[0][c] =
779 CLIP((g + rix[hex[d]][c] + rix[hex[d + 1]][c]) / 2);
780 }
781 }
782 }
783 rgb = (ushort(*)[TS][TS][3]) buffer;
784 mrow -= top;
785 mcol -= left;
786
787 /* Convert to CIELab and differentiate in all directions: */
788 for (d = 0; d < ndir; d++) {
789 for (row = 2; row < mrow - 2; row++)
790 for (col = 2; col < mcol - 2; col++)
791 cielab_INDI(rgb[d][row][col], lab[row][col], colors, rgb_cam);
792 for (f = dir[d & 3], row = 3; row < mrow - 3; row++)
793 for (col = 3; col < mcol - 3; col++) {
794 lix = &lab[row][col];
795 g = 2 * lix[0][0] - lix[f][0] - lix[-f][0];
796 drv[d][row][col] = SQR(g)
797 + SQR((2 * lix[0][1] - lix[f][1] - lix[-f][1] + g * 500 / 232))
798 + SQR((2 * lix[0][2] - lix[f][2] - lix[-f][2] - g * 500 / 580));
799 }
800 }
801
802 /* Build homogeneity maps from the derivatives: */
803 memset(homo, 0, ndir * TS * TS);
804 for (row = 4; row < mrow - 4; row++)
805 for (col = 4; col < mcol - 4; col++) {
806 for (tr = FLT_MAX, d = 0; d < ndir; d++)
807 if (tr > drv[d][row][col])
808 tr = drv[d][row][col];
809 tr *= 8;
810 for (d = 0; d < ndir; d++)
811 for (v = -1; v <= 1; v++)
812 for (h = -1; h <= 1; h++)
813 if (drv[d][row + v][col + h] <= tr)
814 homo[d][row][col]++;
815 }
816
817 /* Average the most homogenous pixels for the final result: */
818 if (height - top < TS + 4) mrow = height - top + 2;
819 if (width - left < TS + 4) mcol = width - left + 2;
820 for (row = MIN(top, 8); row < mrow - 8; row++)
821 for (col = MIN(left, 8); col < mcol - 8; col++) {
822 for (d = 0; d < ndir; d++)
823 for (hm[d] = 0, v = -2; v <= 2; v++)
824 for (h = -2; h <= 2; h++)
825 hm[d] += homo[d][row + v][col + h];
826 for (d = 0; d < ndir - 4; d++)
827 if (hm[d] < hm[d + 4]) hm[d ] = 0;
828 else if (hm[d] > hm[d + 4]) hm[d + 4] = 0;
829 for (max = hm[0], d = 1; d < ndir; d++)
830 if (max < hm[d]) max = hm[d];
831 max -= max >> 3;
832 memset(avg, 0, sizeof avg);
833 for (d = 0; d < ndir; d++)
834 if (hm[d] >= max) {
835 FORC3 avg[c] += rgb[d][row][col][c];
836 avg[3]++;
837 }
838 FORC3 image[(row + top)*width + col + left][c] = avg[c] / avg[3];
839 }
840 }
841 }
842 free(buffer);
843 } /* _OPENMP */
844 border_interpolate_INDI(height, width, image, filters, colors, 8, hh);
845 }
846
847 /*
848 Adaptive Homogeneity-Directed interpolation is based on
849 the work of Keigo Hirakawa, Thomas Parks, and Paul Lee.
850 */
ahd_interpolate_INDI(ushort (* image)[4],const unsigned filters,const int width,const int height,const int colors,const float rgb_cam[3][4],void * dcraw,dcraw_data * h)851 void CLASS ahd_interpolate_INDI(ushort(*image)[4], const unsigned filters,
852 const int width, const int height,
853 const int colors, const float rgb_cam[3][4],
854 void *dcraw, dcraw_data *h)
855 {
856 int i, j, top, left, row, col, tr, tc, c, d, val, hm[2];
857 static const int dir[4] = { -1, 1, -TS, TS };
858 unsigned ldiff[2][4], abdiff[2][4], leps, abeps;
859 ushort(*rgb)[TS][TS][3], (*rix)[3], (*pix)[4];
860 short(*lab)[TS][TS][3], (*lix)[3];
861 char(*homo)[TS][TS], *buffer;
862
863 dcraw_message(dcraw, DCRAW_VERBOSE, _("AHD interpolation...\n")); /*UF*/
864
865 #ifdef _OPENMP
866 #pragma omp parallel \
867 default(shared) \
868 private(top, left, row, col, pix, rix, lix, c, val, d, tc, tr, i, j, ldiff, abdiff, leps, abeps, hm, buffer, rgb, lab, homo)
869 #endif
870 {
871 cielab_INDI(0, 0, colors, rgb_cam);
872 border_interpolate_INDI(height, width, image, filters, colors, 5, h);
873 buffer = (char *) malloc(26 * TS * TS);
874 merror(buffer, "ahd_interpolate()");
875 rgb = (ushort(*)[TS][TS][3]) buffer;
876 lab = (short(*)[TS][TS][3])(buffer + 12 * TS * TS);
877 homo = (char(*)[TS][TS])(buffer + 24 * TS * TS);
878
879 progress(PROGRESS_INTERPOLATE, -height);
880 #ifdef _OPENMP
881 #pragma omp for
882 #endif
883 for (top = 2; top < height - 5; top += TS - 6) {
884 progress(PROGRESS_INTERPOLATE, TS - 6);
885 for (left = 2; left < width - 5; left += TS - 6) {
886
887 /* Interpolate green horizontally and vertically: */
888 for (row = top; row < top + TS && row < height - 2; row++) {
889 col = left + (FC(row, left) & 1);
890 for (c = FC(row, col); col < left + TS && col < width - 2; col += 2) {
891 pix = image + row * width + col;
892 val = ((pix[-1][1] + pix[0][c] + pix[1][1]) * 2
893 - pix[-2][c] - pix[2][c]) >> 2;
894 rgb[0][row - top][col - left][1] = ULIM(val, pix[-1][1], pix[1][1]);
895 val = ((pix[-width][1] + pix[0][c] + pix[width][1]) * 2
896 - pix[-2 * width][c] - pix[2 * width][c]) >> 2;
897 rgb[1][row - top][col - left][1] = ULIM(val, pix[-width][1], pix[width][1]);
898 }
899 }
900 /* Interpolate red and blue, and convert to CIELab: */
901 for (d = 0; d < 2; d++)
902 for (row = top + 1; row < top + TS - 1 && row < height - 3; row++)
903 for (col = left + 1; col < left + TS - 1 && col < width - 3; col++) {
904 pix = image + row * width + col;
905 rix = &rgb[d][row - top][col - left];
906 lix = &lab[d][row - top][col - left];
907 if ((c = 2 - FC(row, col)) == 1) {
908 c = FC(row + 1, col);
909 val = pix[0][1] + ((pix[-1][2 - c] + pix[1][2 - c]
910 - rix[-1][1] - rix[1][1]) >> 1);
911 rix[0][2 - c] = CLIP(val);
912 val = pix[0][1] + ((pix[-width][c] + pix[width][c]
913 - rix[-TS][1] - rix[TS][1]) >> 1);
914 } else
915 val = rix[0][1] + ((pix[-width - 1][c] + pix[-width + 1][c]
916 + pix[+width - 1][c] + pix[+width + 1][c]
917 - rix[-TS - 1][1] - rix[-TS + 1][1]
918 - rix[+TS - 1][1] - rix[+TS + 1][1] + 1) >> 2);
919 rix[0][c] = CLIP(val);
920 c = FC(row, col);
921 rix[0][c] = pix[0][c];
922 cielab_INDI(rix[0], lix[0], colors, rgb_cam);
923 }
924 /* Build homogeneity maps from the CIELab images: */
925 memset(homo, 0, 2 * TS * TS);
926 for (row = top + 2; row < top + TS - 2 && row < height - 4; row++) {
927 tr = row - top;
928 for (col = left + 2; col < left + TS - 2 && col < width - 4; col++) {
929 tc = col - left;
930 for (d = 0; d < 2; d++) {
931 lix = &lab[d][tr][tc];
932 for (i = 0; i < 4; i++) {
933 ldiff[d][i] = ABS(lix[0][0] - lix[dir[i]][0]);
934 abdiff[d][i] = SQR(lix[0][1] - lix[dir[i]][1])
935 + SQR(lix[0][2] - lix[dir[i]][2]);
936 }
937 }
938 leps = MIN(MAX(ldiff[0][0], ldiff[0][1]),
939 MAX(ldiff[1][2], ldiff[1][3]));
940 abeps = MIN(MAX(abdiff[0][0], abdiff[0][1]),
941 MAX(abdiff[1][2], abdiff[1][3]));
942 for (d = 0; d < 2; d++)
943 for (i = 0; i < 4; i++)
944 if (ldiff[d][i] <= leps && abdiff[d][i] <= abeps)
945 homo[d][tr][tc]++;
946 }
947 }
948 /* Combine the most homogenous pixels for the final result: */
949 for (row = top + 3; row < top + TS - 3 && row < height - 5; row++) {
950 tr = row - top;
951 for (col = left + 3; col < left + TS - 3 && col < width - 5; col++) {
952 tc = col - left;
953 for (d = 0; d < 2; d++)
954 for (hm[d] = 0, i = tr - 1; i <= tr + 1; i++)
955 for (j = tc - 1; j <= tc + 1; j++)
956 hm[d] += homo[d][i][j];
957 if (hm[0] != hm[1])
958 FORC3 image[row * width + col][c] = rgb[hm[1] > hm[0]][tr][tc][c];
959 else
960 FORC3 image[row * width + col][c] =
961 (rgb[0][tr][tc][c] + rgb[1][tr][tc][c]) >> 1;
962 }
963 }
964 }
965 }
966 free(buffer);
967 } /* _OPENMP */
968 }
969 #undef TS
970
971
972 #define DTOP(x) ((x>65535) ? (unsigned short)65535 : (x<0) ? (unsigned short)0 : (unsigned short) x)
973
974 /*
975 * MG - This comment applies to the 3x3 optimized median function
976 *
977 * The following routines have been built from knowledge gathered
978 * around the Web. I am not aware of any copyright problem with
979 * them, so use it as you want.
980 * N. Devillard - 1998
981 */
982
983 #define PIX_SORT(a,b) { if ((a)>(b)) PIX_SWAP((a),(b)); }
984 #define PIX_SWAP(a,b) { int temp=(a);(a)=(b);(b)=temp; }
985
median9(int * p)986 static inline int median9(int *p)
987 {
988 PIX_SORT(p[1], p[2]) ;
989 PIX_SORT(p[4], p[5]) ;
990 PIX_SORT(p[7], p[8]) ;
991 PIX_SORT(p[0], p[1]) ;
992 PIX_SORT(p[3], p[4]) ;
993 PIX_SORT(p[6], p[7]) ;
994 PIX_SORT(p[1], p[2]) ;
995 PIX_SORT(p[4], p[5]) ;
996 PIX_SORT(p[7], p[8]) ;
997 PIX_SORT(p[0], p[3]) ;
998 PIX_SORT(p[5], p[8]) ;
999 PIX_SORT(p[4], p[7]) ;
1000 PIX_SORT(p[3], p[6]) ;
1001 PIX_SORT(p[1], p[4]) ;
1002 PIX_SORT(p[2], p[5]) ;
1003 PIX_SORT(p[4], p[7]) ;
1004 PIX_SORT(p[4], p[2]) ;
1005 PIX_SORT(p[6], p[4]) ;
1006 PIX_SORT(p[4], p[2]) ;
1007 return (p[4]) ;
1008 }
1009
1010 #undef PIX_SWAP
1011 #undef PIX_SORT
1012
1013 // Just making this function inline speeds up ahd_interpolate_INDI() up to 15%
eahd_median(int row,int col,int color,ushort (* image)[4],const int width)1014 static inline ushort eahd_median(int row, int col, int color,
1015 ushort(*image)[4], const int width)
1016 {
1017 //declare the pixel array
1018 int pArray[9];
1019 int result;
1020
1021 //perform the median filter (this only works for red or blue)
1022 // result = median(R-G)+G or median(B-G)+G
1023 //
1024 // to perform the filter on green, it needs to be the average
1025 // results = (median(G-R)+median(G-B)+R+B)/2
1026
1027 //no checks are done here to speed up the inlining
1028 pArray[0] = image[width * (row) + col + 1][color] - image[width * (row) + col + 1][1];
1029 pArray[1] = image[width * (row - 1) + col + 1][color] - image[width * (row - 1) + col + 1][1];
1030 pArray[2] = image[width * (row - 1) + col ][color] - image[width * (row - 1) + col ][1];
1031 pArray[3] = image[width * (row - 1) + col - 1][color] - image[width * (row - 1) + col - 1][1];
1032 pArray[4] = image[width * (row) + col - 1][color] - image[width * (row) + col - 1][1];
1033 pArray[5] = image[width * (row + 1) + col - 1][color] - image[width * (row + 1) + col - 1][1];
1034 pArray[6] = image[width * (row + 1) + col ][color] - image[width * (row + 1) + col ][1];
1035 pArray[7] = image[width * (row + 1) + col + 1][color] - image[width * (row + 1) + col + 1][1];
1036 pArray[8] = image[width * (row) + col ][color] - image[width * (row) + col ][1];
1037
1038 median9(pArray);
1039 result = pArray[4] + image[width * (row) + col ][1];
1040 return DTOP(result);
1041
1042 }
1043
1044 // Add the color smoothing from Kimmel as suggested in the AHD paper
1045 // Algorithm updated by Michael Goertz
color_smooth(ushort (* image)[4],const int width,const int height,const int passes)1046 void CLASS color_smooth(ushort(*image)[4], const int width, const int height,
1047 const int passes)
1048 {
1049 int row, col;
1050 int row_start = 2;
1051 int row_stop = height - 2;
1052 int col_start = 2;
1053 int col_stop = width - 2;
1054 //interate through all the colors
1055 int count;
1056
1057 ushort *mpix;
1058
1059 for (count = 0; count < passes; count++) {
1060 //perform 3 iterations - seems to be a commonly settled upon number of iterations
1061 #ifdef _OPENMP
1062 #pragma omp parallel for default(shared) private(row,col,mpix)
1063 #endif
1064 for (row = row_start; row < row_stop; row++) {
1065 for (col = col_start; col < col_stop; col++) {
1066 //calculate the median only over the red and blue
1067 //calculating over green seems to offer very little additional quality
1068 mpix = image[row * width + col];
1069 mpix[0] = eahd_median(row, col, 0, image, width);
1070 mpix[2] = eahd_median(row, col, 2, image, width);
1071 }
1072 }
1073 }
1074 }
1075
fuji_rotate_INDI(ushort (** image_p)[4],int * height_p,int * width_p,int * fuji_width_p,const int colors,const double step,void * dcraw)1076 void CLASS fuji_rotate_INDI(ushort(**image_p)[4], int *height_p,
1077 int *width_p, int *fuji_width_p, const int colors,
1078 const double step, void *dcraw)
1079 {
1080 int height = *height_p, width = *width_p, fuji_width = *fuji_width_p; /*UF*/
1081 ushort(*image)[4] = *image_p; /*UF*/
1082 int i, row, col;
1083 float r, c, fr, fc;
1084 int ur, uc;
1085 ushort wide, high, (*img)[4], (*pix)[4];
1086
1087 if (!fuji_width) return;
1088 dcraw_message(dcraw, DCRAW_VERBOSE, _("Rotating image 45 degrees...\n"));
1089 fuji_width = (fuji_width - 1/* + shrink*/)/* >> shrink*/;
1090 wide = fuji_width / step;
1091 high = (height - fuji_width) / step;
1092 img = (ushort(*)[4]) calloc(wide * high, sizeof * img);
1093 merror(img, "fuji_rotate()");
1094
1095 #ifdef _OPENMP
1096 #pragma omp parallel for default(shared) private(row,col,ur,uc,r,c,fr,fc,pix,i)
1097 #endif
1098 for (row = 0; row < high; row++) {
1099 for (col = 0; col < wide; col++) {
1100 ur = r = fuji_width + (row - col) * step;
1101 uc = c = (row + col) * step;
1102 if (ur > height - 2 || uc > width - 2) continue;
1103 fr = r - ur;
1104 fc = c - uc;
1105 pix = image + ur * width + uc;
1106 for (i = 0; i < colors; i++)
1107 img[row * wide + col][i] =
1108 (pix[ 0][i] * (1 - fc) + pix[ 1][i] * fc) * (1 - fr) +
1109 (pix[width][i] * (1 - fc) + pix[width + 1][i] * fc) * fr;
1110 }
1111 }
1112 free(image);
1113 width = wide;
1114 height = high;
1115 image = img;
1116 fuji_width = 0;
1117 *height_p = height; /* INDI - UF*/
1118 *width_p = width;
1119 *fuji_width_p = fuji_width;
1120 *image_p = image;
1121 }
1122
flip_image_INDI(ushort (* image)[4],int * height_p,int * width_p,int flip)1123 void CLASS flip_image_INDI(ushort(*image)[4], int *height_p, int *width_p,
1124 /*const*/ int flip) /*UF*/
1125 {
1126 unsigned *flag;
1127 int size, base, dest, next, row, col;
1128 gint64 *img, hold;
1129 int height = *height_p, width = *width_p;/* INDI - UF*/
1130
1131 // Message is suppressed because error handling is not enabled here.
1132 // dcraw_message (dcraw, DCRAW_VERBOSE,_("Flipping image %c:%c:%c...\n"),
1133 // flip & 1 ? 'H':'0', flip & 2 ? 'V':'0', flip & 4 ? 'T':'0'); /*UF*/
1134
1135 img = (gint64 *) image;
1136 size = height * width;
1137 flag = calloc((size + 31) >> 5, sizeof * flag);
1138 merror(flag, "flip_image()");
1139 for (base = 0; base < size; base++) {
1140 if (flag[base >> 5] & (1 << (base & 31)))
1141 continue;
1142 dest = base;
1143 hold = img[base];
1144 while (1) {
1145 if (flip & 4) {
1146 row = dest % height;
1147 col = dest / height;
1148 } else {
1149 row = dest / width;
1150 col = dest % width;
1151 }
1152 if (flip & 2)
1153 row = height - 1 - row;
1154 if (flip & 1)
1155 col = width - 1 - col;
1156 next = row * width + col;
1157 if (next == base) break;
1158 flag[next >> 5] |= 1 << (next & 31);
1159 img[dest] = img[next];
1160 dest = next;
1161 }
1162 img[dest] = hold;
1163 }
1164 free(flag);
1165 if (flip & 4) SWAP(height, width);
1166 *height_p = height; /* INDI - UF*/
1167 *width_p = width;
1168 }
1169