--- dcraw.c	2006-11-04 17:04:42.000000000 +0100
+++ dcraw_lchblend.c	2006-11-20 19:52:15.000000000 +0100
@@ -22,6 +22,15 @@
    $Revision: 1.355 $
    $Date: 2006/11/04 15:55:41 $
  */
+/* 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.
+ * 2006/11/04 */
+
+
 
 #define VERSION "8.42"
 
@@ -76,6 +85,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;
@@ -107,6 +118,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 },
@@ -3243,6 +3255,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, pre_mul[] =", black);
@@ -3501,6 +3515,49 @@ void CLASS cam_to_cielab (ushort cam[4],
   }
 }
 
+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.
@@ -3671,6 +3728,64 @@ void CLASS bilateral_filter()
   free (window);
 }
 
+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; row<height; row++)
+  {
+    for(col=0; col<width; col++)
+    {
+      is_sat=0;
+      FORC4 {
+        val = CLIP(image[row*width+col][c] * scale_correction);
+	pix_clipped[c] = CLIP(val / scale_correction);
+	pix_diff = image[row*width+col][c] - pix_clipped[c];
+	if(pix_diff>2 || 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]<lch_clip[0])
+	  lch[0]=lch_clip[0];
+
+	cielch_to_cielab(lab, lch);
+	cielab_to_cam(image[row*width+col], lab);
+      }
+    }
+  }
+}
+
+
 #define SCALE (4 >> shrink)
 void CLASS recover_highlights()
 {
@@ -6875,7 +6990,7 @@ int CLASS main (int argc, char **argv)
     "\n-r <nums> Set raw white balance (four values required)"
     "\n-b <num>  Adjust brightness (default = 1.0)"
     "\n-k <num>  Set black point"
-    "\n-H [0-9]  Highlight mode (0=clip, 1=no clip, 2+=recover)"
+    "\n-H [0-10] Highlight mode (0=clip, 1=no clip, 2-9=recover, 10=blend)"
     "\n-t [0-7]  Flip image (0=none, 3=180, 5=90CCW, 6=90CW)"
     "\n-o [0-5]  Output colorspace (raw,sRGB,Adobe,Wide,ProPhoto,XYZ)"
 #ifndef NO_LCMS
@@ -7106,7 +7221,13 @@ next:
       else ahd_interpolate();
     }
     if (sigma_d > 0 && sigma_r > 0) bilateral_filter();
-    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();
 #ifndef NO_LCMS
     if (cam_profile) apply_profile (cam_profile, out_profile);