--- dcraw.c 2007-03-26 01:26:49.000000000 +0200 +++ dcraw_new.c 2007-04-14 07:41:16.000000000 +0200 @@ -21,6 +21,13 @@ $Revision: 1.377 $ $Date: 2007/03/25 22:56:19 $ */ +/* Modified by Cyril Guyot to blend an image with clipped highlights + * and one with unclipped highlights in LCh space in order to recover + * highlight detail. + * For a saturated pixel, the final LCh coordinates are: L from the + * unclipped one, C from the clipped one and h from the unclipped one. + * To activate LCh blending, use -H 10. + * 2007/04/14 */ #define VERSION "8.69" @@ -81,6 +88,8 @@ typedef unsigned long long UINT64; #define LONG_BIT (8 * sizeof (long)) #endif +#define HIGHLIGHT_BLEND 10 + #define ushort UshORt typedef unsigned char uchar; typedef unsigned short ushort; @@ -113,6 +122,7 @@ int verbose=0, use_auto_wb=0, use_camera int output_color=1, output_bps=8, output_tiff=0; int fuji_layout, fuji_secondary, shot_select=0; float cam_mul[4], pre_mul[4], rgb_cam[3][4]; /* RGB from camera color */ +float scale_correction = 1; const double xyz_rgb[3][3] = { /* XYZ from RGB */ { 0.412453, 0.357580, 0.180423 }, { 0.212671, 0.715160, 0.072169 }, @@ -3438,6 +3448,8 @@ skip_block: dmax = pre_mul[c]; } if (!highlight) dmax = dmin; + if (highlight==HIGHLIGHT_BLEND) + scale_correction = dmax / dmin; FORC4 scale_mul[c] = (pre_mul[c] /= dmax) * 65535.0 / maximum; if (verbose) { fprintf (stderr,_("Scaling with black %d, multipliers"), dblack); @@ -3688,6 +3700,80 @@ void CLASS vng_interpolate() free (code[0][0]); } +void CLASS cam_to_cielab (ushort cam[4], float lab[3]) +{ + int c, i, j, k; + float r, xyz[3]; + static float cbrt[0x10000], xyz_cam[3][4]; + + if (cam == NULL) { + for (i=0; i < 0x10000; i++) { + r = i / 65535.0; + cbrt[i] = r > 0.008856 ? pow(r,1/3.0) : 7.787*r + 16/116.0; + } + for (i=0; i < 3; i++) + for (j=0; j < colors; j++) + for (xyz_cam[i][j] = k=0; k < 3; k++) + xyz_cam[i][j] += xyz_rgb[i][k] * rgb_cam[k][j] / d65_white[i]; + } else { + xyz[0] = xyz[1] = xyz[2] = 0.5; + FORCC { + xyz[0] += xyz_cam[0][c] * cam[c]; + xyz[1] += xyz_cam[1][c] * cam[c]; + xyz[2] += xyz_cam[2][c] * cam[c]; + } + xyz[0] = cbrt[CLIP((int) xyz[0])]; + xyz[1] = cbrt[CLIP((int) xyz[1])]; + xyz[2] = cbrt[CLIP((int) xyz[2])]; + lab[0] = 116 * xyz[1] - 16; + lab[1] = 500 * (xyz[0] - xyz[1]); + lab[2] = 200 * (xyz[1] - xyz[2]); + } +} + +void CLASS cielab_to_cam (ushort cam[4], float lab[3]) +{ + int c, i, j, k; + float xyz[3], camf[3], fx, fy, fz, xr, yr, zr, kappa, epsilon; + double xyz_cam[3][4], xyz_cam_tricol[3][3]; + static double cam_xyz[3][3]; + + if (cam == NULL) { + for (i=0; i < 3; i++) + for (j=0; j < colors; j++) + { + for (xyz_cam[i][j] = k=0; k < 3; k++) + xyz_cam[i][j] += xyz_rgb[i][k] * rgb_cam[k][j] / d65_white[i]; + if(j < 3) + xyz_cam_tricol[i][j] = xyz_cam[i][j]; + else + xyz_cam_tricol[i][1] += xyz_cam[i][3]; + } + pseudoinverse(xyz_cam_tricol,cam_xyz, 3); + } else { + epsilon = 0.008856; kappa = 903.3; + yr = (lab[0]<=kappa*epsilon) ? (lab[0]/kappa) : (pow((lab[0]+16.0)/116.0, 3.0)); + fy = (yr<=epsilon) ? ((kappa*yr+16.0)/116.0) : ((lab[0]+16.0)/116.0); + fz = fy - lab[2]/200.0; + fx = lab[1]/500.0 + fy; + zr = (pow(fz, 3.0)<=epsilon) ? ((116.0*fz-16.0)/kappa) : (pow(fz, 3.0)); + xr = (pow(fx, 3.0)<=epsilon) ? ((116.0*fx-16.0)/kappa) : (pow(fx, 3.0)); + + xyz[0] = xr*65535.0 - 0.5; + xyz[1] = yr*65535.0 - 0.5; + xyz[2] = zr*65535.0 - 0.5; + + camf[0]=camf[1]=camf[2]=0; + FORC3 { + camf[0] += cam_xyz[c][0] * xyz[c]; + camf[1] += cam_xyz[c][1] * xyz[c]; + camf[2] += cam_xyz[c][2] * xyz[c]; + } + FORC3 cam[c]=(ushort)camf[c]; + cam[3]=(ushort)camf[1]; + } +} + /* Adaptive Homogeneity-Directed interpolation is based on the work of Keigo Hirakawa, Thomas Parks, and Paul Lee. @@ -3819,6 +3905,63 @@ void CLASS ahd_interpolate() } #undef TS +void CLASS cielab_to_cielch (float lab[3], float lch[3]) +{ + lch[0] = lab[0]; + lch[1] = sqrt (lab[1]*lab[1] + lab[2]*lab[2]); + lch[2] = atan2 (lab[2], lab[1]); +} + +void CLASS cielch_to_cielab (float lab[3], float lch[3]) +{ + lab[0] = lch[0]; + lab[1] = lch[1] * cos(lch[2]); + lab[2] = lch[1] * sin(lch[2]); +} + +void CLASS recover_highlights_blend() +{ + float lab[3], lch[3], lch_clip[3]; + ushort pix_clipped[4]; + int val, c, pix_diff, is_sat; + unsigned row, col; + + if (verbose) fprintf (stderr, "Highlight recovery: LCh blending...\n"); + + cielab_to_cam (NULL, NULL); + + for(row=0; row2 || pix_diff<-2) is_sat=1; + } + + if (is_sat) + { + cam_to_cielab (image[row*width+col], lab); + cielab_to_cielch (lab, lch); + + cam_to_cielab (pix_clipped, lab); + cielab_to_cielch (lab, lch_clip); + + if(lch[1]>lch_clip[1]) + lch[1]=lch_clip[1]; + if(lch[0]> shrink) void CLASS recover_highlights() { @@ -7218,7 +7361,7 @@ int CLASS main (int argc, char **argv) puts(_("-n Set threshold for wavelet denoising")); puts(_("-k Set black point")); puts(_("-K Subtract dark frame (16-bit raw PGM)")); - puts(_("-H [0-9] Highlight mode (0=clip, 1=no clip, 2+=recover)")); + puts(_("-H [0-10] Highlight mode (0=clip, 1=no clip, 2-9=recover, 10=blend)")); puts(_("-t [0-7] Flip image (0=none, 3=180, 5=90CCW, 6=90CW)")); puts(_("-o [0-5] Output colorspace (raw,sRGB,Adobe,Wide,ProPhoto,XYZ)")); #ifndef NO_LCMS @@ -7447,6 +7590,7 @@ next: if (is_foveon && !document_mode) foveon_interpolate(); if (!is_foveon && document_mode < 2) scale_colors(); pre_interpolate(); + cam_to_cielab (NULL,NULL); if (filters && !document_mode) { if (quality == 0) lin_interpolate(); @@ -7454,7 +7598,13 @@ next: vng_interpolate(); else ahd_interpolate(); } - if (!is_foveon && highlight > 1) recover_highlights(); + if (!is_foveon && highlight > 1) + { + if (highlight == HIGHLIGHT_BLEND) + recover_highlights_blend(); + else + recover_highlights(); + } if (use_fuji_rotate) fuji_rotate(); if (mix_green && (colors = 3)) for (i=0; i < height*width; i++)