1 /* This code is licensed to the public by its copyright owners under GPL. */
2
3 #define _XOPEN_SOURCE 500 /* get M_PI in math.h */
4
5 #include <assert.h>
6 #include <string.h>
7 #include <float.h>
8 #include <math.h>
9
10 #include "mallocvar.h"
11 #include "pm.h"
12 #include "global_variables.h"
13 #include "decode.h"
14 #include "foveon.h"
15 #include "stdio_nofail.h"
16
17 #if HAVE_INT64
18 typedef int64_t INT64;
19 static bool const have64BitArithmetic = true;
20 #else
21 /* We define INT64 to something that lets the code compile, but we
22 should not execute any INT64 code, because it will get the wrong
23 result.
24 */
25 typedef int INT64;
26 static bool const have64BitArithmetic = false;
27 #endif
28
29
30 #define FORC3 for (c=0; c < 3; c++)
31 #define FORC4 for (c=0; c < colors; c++)
32
33
34
35 static char *
foveon_gets(int const offset,char * const str,int const len)36 foveon_gets(int const offset,
37 char * const str,
38 int const len) {
39
40 unsigned int i;
41 fseek_nofail (ifp, offset, SEEK_SET);
42 for (i=0; i < len-1; ++i) {
43 /* It certains seems wrong that we're reading a 16 bit integer
44 and assigning it to char, but that's what Dcraw does.
45 */
46 unsigned short val;
47 pm_readlittleshortu(ifp, &val);
48 str[i] = val;
49 if (str[i] == 0)
50 break;
51 }
52 str[i] = 0;
53 return str;
54 }
55
56
57
58 void
parse_foveon(FILE * const ifp)59 parse_foveon(FILE * const ifp) {
60 long fliplong;
61 long pos;
62 long magic;
63 long junk;
64 long entries;
65
66 fseek_nofail (ifp, 36, SEEK_SET);
67 pm_readlittlelong(ifp, &fliplong);
68 flip = fliplong;
69 fseek_nofail (ifp, -4, SEEK_END);
70 pm_readlittlelong(ifp, &pos);
71 fseek_nofail (ifp, pos, SEEK_SET);
72 pm_readlittlelong(ifp, &magic);
73 if (magic != 0x64434553)
74 return; /* SECd */
75 pm_readlittlelong(ifp, &junk);
76 pm_readlittlelong(ifp, &entries);
77 while (entries--) {
78 long off, len, tag;
79 long save;
80 long pent;
81 int i, poff[256][2];
82 char name[64];
83 long sec_;
84
85 pm_readlittlelong(ifp, &off);
86 pm_readlittlelong(ifp, &len);
87 pm_readlittlelong(ifp, &tag);
88
89 save = ftell_nofail(ifp);
90 fseek_nofail (ifp, off, SEEK_SET);
91 pm_readlittlelong(ifp, &sec_);
92 if (sec_ != (0x20434553 | (tag << 24))) return;
93 switch (tag) {
94 case 0x47414d49: /* IMAG */
95 if (data_offset)
96 break;
97 data_offset = off + 28;
98 fseek_nofail (ifp, 12, SEEK_CUR);
99 {
100 long wlong, hlong;
101 pm_readlittlelong(ifp, &wlong);
102 pm_readlittlelong(ifp, &hlong);
103 raw_width = wlong;
104 raw_height = hlong;
105 }
106 break;
107 case 0x464d4143: /* CAMF */
108 meta_offset = off + 24;
109 meta_length = len - 28;
110 if (meta_length > 0x20000)
111 meta_length = 0x20000;
112 break;
113 case 0x504f5250: /* PROP */
114 pm_readlittlelong(ifp, &junk);
115 pm_readlittlelong(ifp, &pent);
116 fseek_nofail (ifp, 12, SEEK_CUR);
117 off += pent*8 + 24;
118 if (pent > 256) pent=256;
119 for (i=0; i < pent*2; i++) {
120 long x;
121 pm_readlittlelong(ifp, &x);
122 poff[0][i] = off + x*2;
123 }
124 for (i=0; i < pent; i++) {
125 foveon_gets (poff[i][0], name, 64);
126 if (!strcmp (name, "CAMMANUF"))
127 foveon_gets (poff[i][1], make, 64);
128 if (!strcmp (name, "CAMMODEL"))
129 foveon_gets (poff[i][1], model, 64);
130 if (!strcmp (name, "WB_DESC"))
131 foveon_gets (poff[i][1], model2, 64);
132 if (!strcmp (name, "TIME"))
133 timestamp = atoi (foveon_gets (poff[i][1], name, 64));
134 }
135 }
136 fseek_nofail (ifp, save, SEEK_SET);
137 }
138 is_foveon = 1;
139 }
140
141
142
143 void
foveon_coeff(int * const useCoeffP,float coeff[3][4])144 foveon_coeff(int * const useCoeffP,
145 float coeff[3][4]) {
146
147 static const float foveon[3][3] =
148 { { 1.4032, -0.2231, -0.1016 },
149 { -0.5263, 1.4816, 0.0170 },
150 { -0.0112, 0.0183, 0.9113 } };
151
152 int i, j;
153
154 for (i=0; i < 3; ++i)
155 for (j=0; j < 3; ++j)
156 coeff[i][j] = foveon[i][j];
157 *useCoeffP = 1;
158 }
159
160
161
162 static void
foveon_decoder(unsigned int const huff[1024],unsigned int const code)163 foveon_decoder (unsigned int const huff[1024],
164 unsigned int const code) {
165
166 struct decode *cur;
167 int len;
168 unsigned int code2;
169
170 cur = free_decode++;
171 if (free_decode > first_decode+2048) {
172 pm_error ("decoder table overflow");
173 }
174 if (code) {
175 unsigned int i;
176 for (i=0; i < 1024; i++) {
177 if (huff[i] == code) {
178 cur->leaf = i;
179 return;
180 }
181 }
182 }
183 if ((len = code >> 27) > 26)
184 return;
185 code2 = (len+1) << 27 | (code & 0x3ffffff) << 1;
186
187 cur->branch[0] = free_decode;
188 foveon_decoder (huff, code2);
189 cur->branch[1] = free_decode;
190 foveon_decoder (huff, code2+1);
191 }
192
193
194
195 static void
foveon_load_camf()196 foveon_load_camf() {
197 unsigned int i, val;
198 long key;
199
200 fseek_nofail (ifp, meta_offset, SEEK_SET);
201 pm_readlittlelong(ifp, &key);
202 fread_nofail (meta_data, 1, meta_length, ifp);
203 for (i=0; i < meta_length; i++) {
204 key = (key * 1597 + 51749) % 244944;
205 assert(have64BitArithmetic);
206 val = key * (INT64) 301593171 >> 24;
207 meta_data[i] ^= ((((key << 8) - val) >> 1) + val) >> 17;
208 }
209 }
210
211
212
213 void
foveon_load_raw(Image const image)214 foveon_load_raw(Image const image) {
215
216 struct decode *dindex;
217 short diff[1024], pred[3];
218 unsigned int huff[1024], bitbuf=0;
219 int row, col, bit=-1, c, i;
220
221 assert(have64BitArithmetic);
222
223 for (i=0; i < 1024; ++i)
224 pm_readlittleshort(ifp, &diff[i]);
225
226 for (i=0; i < 1024; ++i) {
227 long x;
228 pm_readlittlelong(ifp, &x);
229 huff[i] = x;
230 }
231 init_decoder();
232 foveon_decoder (huff, 0);
233
234 for (row=0; row < height; row++) {
235 long junk;
236 memset (pred, 0, sizeof pred);
237 if (!bit)
238 pm_readlittlelong(ifp, &junk);
239 for (col=bit=0; col < width; col++) {
240 FORC3 {
241 for (dindex=first_decode; dindex->branch[0]; ) {
242 if ((bit = (bit-1) & 31) == 31)
243 for (i=0; i < 4; i++)
244 bitbuf = (bitbuf << 8) + fgetc_nofail(ifp);
245 dindex = dindex->branch[bitbuf >> bit & 1];
246 }
247 pred[c] += diff[dindex->leaf];
248 }
249 FORC3 image[row*width+col][c] = pred[c];
250 }
251 }
252 foveon_load_camf();
253 maximum = clip_max = 0xffff;
254 }
255
256
257
258 static int
sget4(char const s[])259 sget4(char const s[]) {
260 return s[0] | s[1] << 8 | s[2] << 16 | s[3] << 24;
261 }
262
263
264
265 static char *
foveon_camf_param(const char * const block,const char * const param)266 foveon_camf_param (const char * const block,
267 const char * const param) {
268 unsigned idx, num;
269 char *pos, *cp, *dp;
270
271 for (idx=0; idx < meta_length; idx += sget4(pos+8)) {
272 pos = meta_data + idx;
273 if (strncmp (pos, "CMb", 3)) break;
274 if (pos[3] != 'P') continue;
275 if (strcmp (block, pos+sget4(pos+12))) continue;
276 cp = pos + sget4(pos+16);
277 num = sget4(cp);
278 dp = pos + sget4(cp+4);
279 while (num--) {
280 cp += 8;
281 if (!strcmp (param, dp+sget4(cp)))
282 return dp+sget4(cp+4);
283 }
284 }
285 return NULL;
286 }
287
288
289
290 static void *
foveon_camf_matrix(int dim[3],const char * const name)291 foveon_camf_matrix (int dim[3],
292 const char * const name) {
293
294 unsigned i, idx, type, ndim, size, *mat;
295 char *pos, *cp, *dp;
296
297 for (idx=0; idx < meta_length; idx += sget4(pos+8)) {
298 pos = meta_data + idx;
299 if (strncmp (pos, "CMb", 3)) break;
300 if (pos[3] != 'M') continue;
301 if (strcmp (name, pos+sget4(pos+12))) continue;
302 dim[0] = dim[1] = dim[2] = 1;
303 cp = pos + sget4(pos+16);
304 type = sget4(cp);
305 if ((ndim = sget4(cp+4)) > 3) break;
306 dp = pos + sget4(cp+8);
307 for (i=ndim; i--; ) {
308 cp += 12;
309 dim[i] = sget4(cp);
310 }
311 if ((size = dim[0]*dim[1]*dim[2]) > meta_length/4) break;
312 MALLOCARRAY(mat, size);
313 if (mat == NULL)
314 pm_error("Unable to allocate memory for size=%d", size);
315 for (i=0; i < size; i++)
316 if (type && type != 6)
317 mat[i] = sget4(dp + i*4);
318 else
319 mat[i] = sget4(dp + i*2) & 0xffff;
320 return mat;
321 }
322 pm_message ("'%s' matrix not found!", name);
323 return NULL;
324 }
325
326
327
328 static int
foveon_fixed(void * const ptr,int const size,const char * const name)329 foveon_fixed (void * const ptr,
330 int const size,
331 const char * const name) {
332 void *dp;
333 int dim[3];
334
335 dp = foveon_camf_matrix (dim, name);
336 if (!dp) return 0;
337 memcpy (ptr, dp, size*4);
338 free (dp);
339 return 1;
340 }
341
foveon_avg(unsigned short * pix,int range[2],float cfilt)342 static float foveon_avg (unsigned short *pix, int range[2], float cfilt)
343 {
344 int i;
345 float val, min=FLT_MAX, max=-FLT_MAX, sum=0;
346
347 for (i=range[0]; i <= range[1]; i++) {
348 sum += val =
349 (short)pix[i*4] + ((short)pix[i*4]-(short)pix[(i-1)*4]) * cfilt;
350 if (min > val) min = val;
351 if (max < val) max = val;
352 }
353 return (sum - min - max) / (range[1] - range[0] - 1);
354 }
355
foveon_make_curve(double max,double mul,double filt)356 static short *foveon_make_curve (double max, double mul, double filt)
357 {
358 short *curve;
359 int size;
360 unsigned int i;
361 double x;
362
363 size = 4*M_PI*max / filt;
364 MALLOCARRAY(curve, size+1);
365 if (curve == NULL)
366 pm_error("Out of memory for %d-element curve array", size);
367
368 curve[0] = size;
369 for (i=0; i < size; ++i) {
370 x = i*filt/max/4;
371 curve[i+1] = (cos(x)+1)/2 * tanh(i*filt/mul) * mul + 0.5;
372 }
373 return curve;
374 }
375
foveon_make_curves(short ** curvep,float dq[3],float div[3],float filt)376 static void foveon_make_curves
377 (short **curvep, float dq[3], float div[3], float filt)
378 {
379 double mul[3], max=0;
380 int c;
381
382 FORC3 mul[c] = dq[c]/div[c];
383 FORC3 if (max < mul[c]) max = mul[c];
384 FORC3 curvep[c] = foveon_make_curve (max, mul[c], filt);
385 }
386
foveon_apply_curve(short * curve,int i)387 static int foveon_apply_curve (short *curve, int i)
388 {
389 if (abs(i) >= (unsigned short)curve[0]) return 0;
390 return i < 0 ? -(unsigned short)curve[1-i] : (unsigned short)curve[1+i];
391 }
392
393 void
foveon_interpolate(Image const image,float coeff[3][4])394 foveon_interpolate(Image const image,
395 float coeff[3][4]) {
396
397 static const short hood[] = {
398 -1,-1, -1,0, -1,1, 0,-1, 0,1, 1,-1, 1,0, 1,1 };
399 short *pix, prev[3], *curve[8], (*shrink)[3];
400 float cfilt=0.8, ddft[3][3][2], ppm[3][3][3];
401 float cam_xyz[3][3], correct[3][3], last[3][3], trans[3][3];
402 float chroma_dq[3], color_dq[3], diag[3][3], div[3];
403 float (*black)[3], (*sgain)[3], (*sgrow)[3];
404 float fsum[3], val, frow, num;
405 int row, col, c, i, j, diff, sgx, irow, sum, min, max, limit;
406 int dim[3], dscr[2][2], (*smrow[7])[3], total[4], ipix[3];
407 int work[3][3], smlast, smred, smred_p=0, dev[3];
408 int satlev[3], keep[4], active[4];
409 unsigned *badpix;
410 double dsum=0, trsum[3];
411 char str[128], *cp;
412
413 foveon_fixed (dscr, 4, "DarkShieldColRange");
414 foveon_fixed (ppm[0][0], 27, "PostPolyMatrix");
415 foveon_fixed (ddft[1][0], 12, "DarkDrift");
416 foveon_fixed (&cfilt, 1, "ColumnFilter");
417 foveon_fixed (satlev, 3, "SaturationLevel");
418 foveon_fixed (keep, 4, "KeepImageArea");
419 foveon_fixed (active, 4, "ActiveImageArea");
420 foveon_fixed (chroma_dq, 3, "ChromaDQ");
421 foveon_fixed (color_dq, 3,
422 foveon_camf_param ("IncludeBlocks", "ColorDQ") ?
423 "ColorDQ" : "ColorDQCamRGB");
424
425 if (!(cp = foveon_camf_param ("WhiteBalanceIlluminants", model2)))
426 { pm_message ( "Invalid white balance \"%s\"", model2);
427 return; }
428 foveon_fixed (cam_xyz, 9, cp);
429 foveon_fixed (correct, 9,
430 foveon_camf_param ("WhiteBalanceCorrections", model2));
431 memset (last, 0, sizeof last);
432 for (i=0; i < 3; i++)
433 for (j=0; j < 3; j++)
434 FORC3 last[i][j] += correct[i][c] * cam_xyz[c][j];
435
436 sprintf (str, "%sRGBNeutral", model2);
437 if (foveon_camf_param ("IncludeBlocks", str))
438 foveon_fixed (div, 3, str);
439 else {
440 #define LAST(x,y) last[(i+x)%3][(c+y)%3]
441 for (i=0; i < 3; i++)
442 FORC3 diag[c][i] = LAST(1,1)*LAST(2,2) - LAST(1,2)*LAST(2,1);
443 #undef LAST
444 FORC3 div[c] =
445 diag[c][0]*0.3127 + diag[c][1]*0.329 + diag[c][2]*0.3583;
446 }
447 num = 0;
448 FORC3 if (num < div[c]) num = div[c];
449 FORC3 div[c] /= num;
450
451 memset (trans, 0, sizeof trans);
452 for (i=0; i < 3; i++)
453 for (j=0; j < 3; j++)
454 FORC3 trans[i][j] += coeff[i][c] * last[c][j] * div[j];
455 FORC3 trsum[c] = trans[c][0] + trans[c][1] + trans[c][2];
456 dsum = (6*trsum[0] + 11*trsum[1] + 3*trsum[2]) / 20;
457 for (i=0; i < 3; i++)
458 FORC3 last[i][c] = trans[i][c] * dsum / trsum[i];
459 memset (trans, 0, sizeof trans);
460 for (i=0; i < 3; i++)
461 for (j=0; j < 3; j++)
462 FORC3 trans[i][j] += (i==c ? 32 : -1) * last[c][j] / 30;
463
464 foveon_make_curves (curve, color_dq, div, cfilt);
465 FORC3 chroma_dq[c] /= 3;
466 foveon_make_curves (curve+3, chroma_dq, div, cfilt);
467 FORC3 dsum += chroma_dq[c] / div[c];
468 curve[6] = foveon_make_curve (dsum, dsum, cfilt);
469 curve[7] = foveon_make_curve (dsum*2, dsum*2, cfilt);
470
471 sgain = foveon_camf_matrix (dim, "SpatialGain");
472 if (!sgain) return;
473 sgrow = calloc (dim[1], sizeof *sgrow);
474 sgx = (width + dim[1]-2) / (dim[1]-1);
475
476 black = calloc (height, sizeof(black[0]));
477 for (row=0; row < height; ++row) {
478 unsigned int i;
479 for (i=0; i < 3; ++i) {
480 unsigned int j;
481 for (j = 0; j < 2; ++j)
482 ddft[0][i][j] = ddft[1][i][j] +
483 row / (height-1.0) * (ddft[2][i][j] - ddft[1][i][j]);
484 }
485 FORC3 black[row][c] =
486 ( foveon_avg (image[row*width]+c, dscr[0], cfilt) +
487 foveon_avg (image[row*width]+c, dscr[1], cfilt) * 3
488 - ddft[0][c][0] ) / 4 - ddft[0][c][1];
489 }
490 memcpy (black, black+8, 8 * sizeof(black[0]));
491 memcpy (black+height-11, black+height-22, 11*(sizeof black[0]));
492 memcpy (last, black, sizeof last);
493
494 for (row=1; row < height-1; row++) {
495 FORC3 if (last[1][c] > last[0][c]) {
496 if (last[1][c] > last[2][c])
497 black[row][c] =
498 (last[0][c] > last[2][c]) ? last[0][c]:last[2][c];
499 } else
500 if (last[1][c] < last[2][c])
501 black[row][c] =
502 (last[0][c] < last[2][c]) ? last[0][c]:last[2][c];
503 memmove (last, last+1, 2*sizeof last[0]);
504 memcpy (last[2], black[row+1], sizeof last[2]);
505 }
506 FORC3 black[row][c] = (last[0][c] + last[1][c])/2;
507 FORC3 black[0][c] = (black[1][c] + black[3][c])/2;
508
509 val = 1 - exp(-1/24.0);
510 memcpy (fsum, black, sizeof fsum);
511 for (row=1; row < height; row++)
512 FORC3 fsum[c] += black[row][c] =
513 (black[row][c] - black[row-1][c])*val + black[row-1][c];
514 memcpy (last[0], black[height-1], sizeof last[0]);
515 FORC3 fsum[c] /= height;
516 for (row = height; row--; )
517 FORC3 last[0][c] = black[row][c] =
518 (black[row][c] - fsum[c] - last[0][c])*val + last[0][c];
519
520 memset (total, 0, sizeof total);
521 for (row=2; row < height; row+=4)
522 for (col=2; col < width; col+=4) {
523 FORC3 total[c] += (short) image[row*width+col][c];
524 total[3]++;
525 }
526 for (row=0; row < height; row++)
527 FORC3 black[row][c] += fsum[c]/2 + total[c]/(total[3]*100.0);
528
529 for (row=0; row < height; row++) {
530 unsigned int i;
531 for (i = 0; i < 3; ++i) {
532 unsigned int j;
533 for (j = 0; j < 2; ++j)
534 ddft[0][i][j] = ddft[1][i][j] +
535 row / (height-1.0) * (ddft[2][i][j] - ddft[1][i][j]);
536 }
537 pix = (short*)image[row*width];
538 memcpy (prev, pix, sizeof prev);
539 frow = row / (height-1.0) * (dim[2]-1);
540 if ((irow = frow) == dim[2]-1) irow--;
541 frow -= irow;
542 for (i=0; i < dim[1]; i++)
543 FORC3 sgrow[i][c] = sgain[ irow *dim[1]+i][c] * (1-frow) +
544 sgain[(irow+1)*dim[1]+i][c] * frow;
545 for (col=0; col < width; col++) {
546 FORC3 {
547 diff = pix[c] - prev[c];
548 prev[c] = pix[c];
549 ipix[c] = pix[c] + floor ((diff + (diff*diff >> 14)) * cfilt
550 - ddft[0][c][1]
551 - ddft[0][c][0]
552 * ((float) col/width - 0.5)
553 - black[row][c] );
554 }
555 FORC3 {
556 work[0][c] = ipix[c] * ipix[c] >> 14;
557 work[2][c] = ipix[c] * work[0][c] >> 14;
558 work[1][2-c] = ipix[(c+1) % 3] * ipix[(c+2) % 3] >> 14;
559 }
560 FORC3 {
561 for (val=i=0; i < 3; i++)
562 for ( j=0; j < 3; j++)
563 val += ppm[c][i][j] * work[i][j];
564 ipix[c] = floor ((ipix[c] + floor(val)) *
565 ( sgrow[col/sgx ][c] * (sgx - col%sgx) +
566 sgrow[col/sgx+1][c] * (col%sgx) ) / sgx /
567 div[c]);
568 if (ipix[c] > 32000) ipix[c] = 32000;
569 pix[c] = ipix[c];
570 }
571 pix += 4;
572 }
573 }
574 free (black);
575 free (sgrow);
576 free (sgain);
577
578 if ((badpix = foveon_camf_matrix (dim, "BadPixels"))) {
579 for (i=0; i < dim[0]; i++) {
580 col = (badpix[i] >> 8 & 0xfff) - keep[0];
581 row = (badpix[i] >> 20 ) - keep[1];
582 if ((unsigned)(row-1) > height-3 || (unsigned)(col-1) > width-3)
583 continue;
584 memset (fsum, 0, sizeof fsum);
585 for (sum=j=0; j < 8; j++)
586 if (badpix[i] & (1 << j)) {
587 FORC3 fsum[c] +=
588 image[(row+hood[j*2])*width+col+hood[j*2+1]][c];
589 sum++;
590 }
591 if (sum) FORC3 image[row*width+col][c] = fsum[c]/sum;
592 }
593 free (badpix);
594 }
595
596 /* Array for 5x5 Gaussian averaging of red values */
597 smrow[6] = calloc (width*5, sizeof **smrow);
598 if (smrow[6] == NULL)
599 pm_error("Out of memory");
600 for (i=0; i < 5; i++)
601 smrow[i] = smrow[6] + i*width;
602
603 /* Sharpen the reds against these Gaussian averages */
604 for (smlast=-1, row=2; row < height-2; row++) {
605 while (smlast < row+2) {
606 for (i=0; i < 6; i++)
607 smrow[(i+5) % 6] = smrow[i];
608 pix = (short*)image[++smlast*width+2];
609 for (col=2; col < width-2; col++) {
610 smrow[4][col][0] =
611 (pix[0]*6 + (pix[-4]+pix[4])*4 + pix[-8]+pix[8] + 8) >> 4;
612 pix += 4;
613 }
614 }
615 pix = (short*)image[row*width+2];
616 for (col=2; col < width-2; col++) {
617 smred = ( 6 * smrow[2][col][0]
618 + 4 * (smrow[1][col][0] + smrow[3][col][0])
619 + smrow[0][col][0] + smrow[4][col][0] + 8 ) >> 4;
620 if (col == 2)
621 smred_p = smred;
622 i = pix[0] + ((pix[0] - ((smred*7 + smred_p) >> 3)) >> 3);
623 if (i > 32000) i = 32000;
624 pix[0] = i;
625 smred_p = smred;
626 pix += 4;
627 }
628 }
629
630 /* Adjust the brighter pixels for better linearity */
631 FORC3 {
632 i = satlev[c] / div[c];
633 if (maximum > i) maximum = i;
634 }
635 clip_max = maximum;
636 limit = maximum * 9 >> 4;
637 for (pix=(short*)image[0]; pix < (short *) image[height*width]; pix+=4) {
638 if (pix[0] <= limit || pix[1] <= limit || pix[2] <= limit)
639 continue;
640 min = max = pix[0];
641 for (c=1; c < 3; c++) {
642 if (min > pix[c]) min = pix[c];
643 if (max < pix[c]) max = pix[c];
644 }
645 i = 0x4000 - ((min - limit) << 14) / limit;
646 i = 0x4000 - (i*i >> 14);
647 i = i*i >> 14;
648 FORC3 pix[c] += (max - pix[c]) * i >> 14;
649 }
650 /*
651 Because photons that miss one detector often hit another,
652 the sum R+G+B is much less noisy than the individual colors.
653 So smooth the hues without smoothing the total.
654 */
655 for (smlast=-1, row=2; row < height-2; row++) {
656 while (smlast < row+2) {
657 for (i=0; i < 6; i++)
658 smrow[(i+5) % 6] = smrow[i];
659 pix = (short*)image[++smlast*width+2];
660 for (col=2; col < width-2; col++) {
661 FORC3 smrow[4][col][c] = (pix[c-4]+2*pix[c]+pix[c+4]+2) >> 2;
662 pix += 4;
663 }
664 }
665 pix = (short*)image[row*width+2];
666 for (col=2; col < width-2; col++) {
667 FORC3 dev[c] = -foveon_apply_curve (curve[7], pix[c] -
668 ((smrow[1][col][c] +
669 2*smrow[2][col][c] +
670 smrow[3][col][c]) >> 2));
671 sum = (dev[0] + dev[1] + dev[2]) >> 3;
672 FORC3 pix[c] += dev[c] - sum;
673 pix += 4;
674 }
675 }
676 for (smlast=-1, row=2; row < height-2; row++) {
677 while (smlast < row+2) {
678 for (i=0; i < 6; i++)
679 smrow[(i+5) % 6] = smrow[i];
680 pix = (short*)image[++smlast*width+2];
681 for (col=2; col < width-2; col++) {
682 FORC3 smrow[4][col][c] =
683 (pix[c-8]+pix[c-4]+pix[c]+pix[c+4]+pix[c+8]+2) >> 2;
684 pix += 4;
685 }
686 }
687 pix = (short*)image[row*width+2];
688 for (col=2; col < width-2; col++) {
689 for (total[3]=375, sum=60, c=0; c < 3; c++) {
690 for (total[c]=i=0; i < 5; i++)
691 total[c] += smrow[i][col][c];
692 total[3] += total[c];
693 sum += pix[c];
694 }
695 if (sum < 0) sum = 0;
696 j = total[3] > 375 ? (sum << 16) / total[3] : sum * 174;
697 FORC3 pix[c] += foveon_apply_curve (curve[6],
698 ((j*total[c] + 0x8000) >> 16)
699 - pix[c]);
700 pix += 4;
701 }
702 }
703
704 /* Transform the image to a different colorspace */
705 for (pix=(short*)image[0]; pix < (short *) image[height*width]; pix+=4) {
706 FORC3 pix[c] -= foveon_apply_curve (curve[c], pix[c]);
707 sum = (pix[0]+pix[1]+pix[1]+pix[2]) >> 2;
708 FORC3 pix[c] -= foveon_apply_curve (curve[c], pix[c]-sum);
709 FORC3 {
710 for (dsum=i=0; i < 3; i++)
711 dsum += trans[c][i] * pix[i];
712 if (dsum < 0) dsum = 0;
713 if (dsum > 24000) dsum = 24000;
714 ipix[c] = dsum + 0.5;
715 }
716 FORC3 pix[c] = ipix[c];
717 }
718
719 /* Smooth the image bottom-to-top and save at 1/4 scale */
720 MALLOCARRAY(shrink, (width/4) * (height/4));
721 if (shrink == NULL)
722 pm_error("Out of memory allocating 1/4 scale array");
723
724 for (row = height/4; row > 0; --row) {
725 for (col=0; col < width/4; ++col) {
726 ipix[0] = ipix[1] = ipix[2] = 0;
727 for (i=0; i < 4; i++)
728 for (j=0; j < 4; j++)
729 FORC3 ipix[c] += image[(row*4+i)*width+col*4+j][c];
730 FORC3
731 if (row+2 > height/4)
732 shrink[row*(width/4)+col][c] = ipix[c] >> 4;
733 else
734 shrink[row*(width/4)+col][c] =
735 (shrink[(row+1)*(width/4)+col][c]*1840 +
736 ipix[c]*141 + 2048) >> 12;
737 }
738 }
739 /* From the 1/4-scale image, smooth right-to-left */
740 for (row=0; row < (height & ~3); ++row) {
741 ipix[0] = ipix[1] = ipix[2] = 0;
742 if ((row & 3) == 0)
743 for (col = width & ~3 ; col--; )
744 FORC3 smrow[0][col][c] = ipix[c] =
745 (shrink[(row/4)*(width/4)+col/4][c]*1485 +
746 ipix[c]*6707 + 4096) >> 13;
747
748 /* Then smooth left-to-right */
749 ipix[0] = ipix[1] = ipix[2] = 0;
750 for (col=0; col < (width & ~3); col++)
751 FORC3 smrow[1][col][c] = ipix[c] =
752 (smrow[0][col][c]*1485 + ipix[c]*6707 + 4096) >> 13;
753
754 /* Smooth top-to-bottom */
755 if (row == 0)
756 memcpy (smrow[2], smrow[1], sizeof **smrow * width);
757 else
758 for (col=0; col < (width & ~3); col++)
759 FORC3 smrow[2][col][c] =
760 (smrow[2][col][c]*6707 + smrow[1][col][c]*1485 + 4096)
761 >> 13;
762
763 /* Adjust the chroma toward the smooth values */
764 for (col=0; col < (width & ~3); col++) {
765 for (i=j=30, c=0; c < 3; c++) {
766 i += smrow[2][col][c];
767 j += image[row*width+col][c];
768 }
769 j = (j << 16) / i;
770 for (sum=c=0; c < 3; c++) {
771 ipix[c] = foveon_apply_curve(
772 curve[c+3],
773 ((smrow[2][col][c] * j + 0x8000) >> 16) -
774 image[row*width+col][c]);
775 sum += ipix[c];
776 }
777 sum >>= 3;
778 FORC3 {
779 i = image[row*width+col][c] + ipix[c] - sum;
780 if (i < 0) i = 0;
781 image[row*width+col][c] = i;
782 }
783 }
784 }
785 free (shrink);
786 free (smrow[6]);
787 for (i=0; i < 8; i++)
788 free (curve[i]);
789
790 /* Trim off the black border */
791 active[1] -= keep[1];
792 active[3] -= 2;
793 i = active[2] - active[0];
794 for (row = 0; row < active[3]-active[1]; row++)
795 memcpy (image[row*i], image[(row+active[1])*width+active[0]],
796 i * sizeof *image);
797 width = i;
798 height = row;
799 }
800