Support DNG 1.4 files with Floating Point data, Issue 2369, thanks to Javier Celaya

rtengine/dcraw.patch

@@ -1,5 +1,5 @@
---- dcraw.c	2014-02-19 17:25:45.051457734 +0100
-+++ dcraw.cc	2014-03-25 10:57:44.977344962 +0100
+--- dcraw.c	2014-05-11 10:30:58 +0000
++++ dcraw.cc	2014-05-11 12:30:50 +0000
 @@ -1,3 +1,15 @@
 +/*RT*/#include <glib.h>
 +/*RT*/#include <glib/gstdio.h>
@@ -609,7 +609,8 @@
 -  int code[16][16][32], size=16, *ip, sum[4];
 -  int f, c, i, x, y, row, col, shift, color;
 -  ushort *pix;
--
++/* RT: delete interpolation functions */
+ 
 -  if (verbose) fprintf (stderr,_("Bilinear interpolation...\n"));
 -  if (filters == 9) size = 6;
 -  border_interpolate(1);
@@ -795,8 +796,7 @@
 -
 -  border_interpolate(3);
 -  if (verbose) fprintf (stderr,_("PPG interpolation...\n"));
-+/* RT: delete interpolation functions */
- 
+-
 -/*  Fill in the green layer with gradients and pattern recognition: */
 -  for (row=3; row < height-3; row++)
 -    for (col=3+(FC(row,3) & 1), c=FC(row,col); col < width-3; col+=2) {
@@ -857,7 +857,7 @@
 -  ushort (*rgb)[TS][TS][3], (*rix)[3], (*pix)[4];
 -   short (*lab)[TS][TS][3], (*lix)[3];
 -   char (*homo)[TS][TS], *buffer;
--
+ 
 -  if (verbose) fprintf (stderr,_("AHD interpolation...\n"));
 -
 -  cielab (0,0);
@@ -870,7 +870,7 @@
 -
 -  for (top=2; top < height-5; top += TS-6)
 -    for (left=2; left < width-5; left += TS-6) {
- 
+-
 -/*  Interpolate green horizontally and vertically:		*/
 -      for (row=top; row < top+TS && row < height-2; row++) {
 -	col = left + (FC(row,left) & 1);
@@ -1011,7 +1011,37 @@
  
    if (tiff_nifds >= sizeof tiff_ifd / sizeof tiff_ifd[0])
      return 1;
-@@ -5660,10 +5509,21 @@
+@@ -5435,6 +5284,9 @@
+       case 315:				/* Artist */
+ 	fread (artist, 64, 1, ifp);
+ 	break;
++      case 317:				/* Predictor */
++	tiff_ifd[ifd].predictor = getint(type);
++	break;
+       case 322:				/* TileWidth */
+ 	tiff_ifd[ifd].tile_width = getint(type);
+ 	break;
+@@ -5448,6 +5300,9 @@
+ 	  is_raw = 5;
+ 	}
+ 	break;
++      case 325:                         /* TileByteCounts */
++        tiff_ifd[ifd].bytes = len > 1 ? ftell(ifp) : get4();
++        break;
+       case 330:				/* SubIFDs */
+ 	if (!strcmp(model,"DSLR-A100") && tiff_ifd[ifd].width == 3872) {
+ 	  load_raw = &CLASS sony_arw_load_raw;
+@@ -5461,6 +5316,9 @@
+ 	  fseek (ifp, i+4, SEEK_SET);
+ 	}
+ 	break;
++      case 339:
++        tiff_ifd[ifd].sample_format = getint(type);
++        break;
+       case 400:
+ 	strcpy (make, "Sarnoff");
+ 	maximum = 0xfff;
+@@ -5660,10 +5518,21 @@
        case 61450:
  	cblack[4] = cblack[5] = MIN(sqrt(len),64);
        case 50714:			/* BlackLevel */
@@ -1037,7 +1067,7 @@
        case 50715:			/* BlackLevelDeltaH */
        case 50716:			/* BlackLevelDeltaV */
  	for (num=i=0; i < len; i++)
-@@ -5741,12 +5601,15 @@
+@@ -5741,12 +5610,15 @@
      fread (buf, sony_length, 1, ifp);
      sony_decrypt (buf, sony_length/4, 1, sony_key);
      sfp = ifp;
@@ -1059,7 +1089,7 @@
      ifp = sfp;
      free (buf);
    }
-@@ -5770,6 +5633,7 @@
+@@ -5770,6 +5642,7 @@
  int CLASS parse_tiff (int base)
  {
    int doff;
@@ -1067,7 +1097,7 @@
  
    fseek (ifp, base, SEEK_SET);
    order = get2();
-@@ -5847,7 +5711,7 @@
+@@ -5847,7 +5720,7 @@
  	  case  8: load_raw = &CLASS eight_bit_load_raw;	break;
  	  case 12: if (tiff_ifd[raw].phint == 2)
  		     load_flags = 6;
@@ -1076,7 +1106,15 @@
  	  case 14: load_flags = 0;
  	  case 16: load_raw = &CLASS unpacked_load_raw;
  		   if (!strncmp(make,"OLYMPUS",7) &&
-@@ -5963,7 +5827,7 @@
+@@ -5878,6 +5751,7 @@
+ 	  case 32803: load_raw = &CLASS kodak_65000_load_raw;
+ 	}
+       case 32867: case 34892: break;
++      case 8: break;
+       default: is_raw = 0;
+     }
+   if (!dng_version)
+@@ -5963,7 +5837,7 @@
  {
    const char *file, *ext;
    char *jname, *jfile, *jext;
@@ -1085,7 +1123,7 @@
  
    ext  = strrchr (ifname, '.');
    file = strrchr (ifname, '/');
-@@ -5985,13 +5849,14 @@
+@@ -5985,13 +5859,14 @@
    } else
      while (isdigit(*--jext)) {
        if (*jext != '9') {
@@ -1102,7 +1140,7 @@
        if (verbose)
  	fprintf (stderr,_("Reading metadata from %s ...\n"), jname);
        parse_tiff (12);
-@@ -6256,6 +6121,8 @@
+@@ -6256,6 +6131,8 @@
        case 0x21d:  ph1.black     = data;		break;
        case 0x222:  ph1.split_col = data;		break;
        case 0x223:  ph1.black_off = data+base;		break;
@@ -1111,7 +1149,7 @@
        case 0x301:
  	model[63] = 0;
  	fread (model, 1, 63, ifp);
-@@ -6334,7 +6201,11 @@
+@@ -6334,7 +6211,11 @@
      order = get2();
      hlen  = get4();
      if (get4() == 0x48454150)		/* "HEAP" */
@@ -1123,7 +1161,7 @@
      if (parse_tiff (save+6)) apply_tiff();
      fseek (ifp, save+len, SEEK_SET);
    }
-@@ -6586,7 +6457,8 @@
+@@ -6586,7 +6467,8 @@
  {
    static const struct {
      const char *prefix;
@@ -1133,7 +1171,7 @@
    } table[] = {
      { "AgfaPhoto DC-833m", 0, 0,	/* DJC */
  	{ 11438,-3762,-1115,-2409,9914,2497,-1227,2295,5300 } },
-@@ -7457,6 +7329,27 @@
+@@ -7457,6 +7339,27 @@
        }
        break;
      }
@@ -1161,7 +1199,7 @@
  }
  
  void CLASS simple_coeff (int index)
-@@ -7764,13 +7657,20 @@
+@@ -7764,13 +7667,20 @@
    fread (head, 1, 32, ifp);
    fseek (ifp, 0, SEEK_END);
    flen = fsize = ftell(ifp);
@@ -1184,7 +1222,7 @@
        parse_ciff (hlen, flen-hlen, 0);
        load_raw = &CLASS canon_load_raw;
      } else if (parse_tiff(0)) apply_tiff();
-@@ -7816,6 +7716,7 @@
+@@ -7816,6 +7726,7 @@
      fseek (ifp, 100+28*(shot_select > 0), SEEK_SET);
      parse_tiff (data_offset = get4());
      parse_tiff (thumb_offset+12);
@@ -1192,7 +1230,7 @@
      apply_tiff();
    } else if (!memcmp (head,"RIFF",4)) {
      fseek (ifp, 0, SEEK_SET);
-@@ -7925,15 +7826,18 @@
+@@ -7925,15 +7836,18 @@
    if (make[0] == 0) parse_smal (0, flen);
    if (make[0] == 0) {
      parse_jpeg(0);
@@ -1220,7 +1258,7 @@
    }
  
    for (i=0; i < sizeof corp / sizeof *corp; i++)
-@@ -7966,7 +7870,7 @@
+@@ -7966,7 +7880,7 @@
    if (height == 3136 && width == 4864)  /* Pentax K20D and Samsung GX20 */
      { height  = 3124;   width  = 4688; filters = 0x16161616; }
    if (width == 4352 && (!strcmp(model,"K-r") || !strcmp(model,"K-x")))
@@ -1229,7 +1267,15 @@
    if (width >= 4960 && !strncmp(model,"K-5",3))
      { left_margin = 10; width  = 4950; filters = 0x16161616; }
    if (width == 4736 && !strcmp(model,"K-7"))
-@@ -8112,7 +8016,7 @@
+@@ -7985,6 +7899,7 @@
+     switch (tiff_compress) {
+       case 1:     load_raw = &CLASS   packed_dng_load_raw;  break;
+       case 7:     load_raw = &CLASS lossless_dng_load_raw;  break;
++      case 8:     load_raw = &CLASS  deflate_dng_load_raw;  break;
+       case 34892: load_raw = &CLASS    lossy_dng_load_raw;  break;
+       default:    load_raw = 0;
+     }
+@@ -8112,7 +8027,7 @@
      width -= 44;
    } else if (!strcmp(model,"D3200") ||
  	     !strcmp(model,"D600")  ||
@@ -1238,7 +1284,7 @@
      width -= 46;
    } else if (!strcmp(model,"D4") ||
  	     !strcmp(model,"Df")) {
-@@ -8394,6 +8298,7 @@
+@@ -8394,6 +8309,7 @@
        filters = 0x16161616;
      }
    } else if (!strcmp(make,"Leica") || !strcmp(make,"Panasonic")) {
@@ -1246,7 +1292,7 @@
      if ((flen - data_offset) / (raw_width*8/7) == raw_height)
        load_raw = &CLASS panasonic_load_raw;
      if (!load_raw) {
-@@ -8411,6 +8316,7 @@
+@@ -8411,6 +8327,7 @@
        }
      filters = 0x01010101 * (uchar) "\x94\x61\x49\x16"
  	[((filters-1) ^ (left_margin & 1) ^ (top_margin << 1)) & 3];
@@ -1254,7 +1300,7 @@
    } else if (!strcmp(model,"C770UZ")) {
      height = 1718;
      width  = 2304;
-@@ -8630,6 +8536,10 @@
+@@ -8630,6 +8547,10 @@
      memcpy (rgb_cam, cmatrix, sizeof cmatrix);
      raw_color = 0;
    }
@@ -1265,7 +1311,16 @@
    if (raw_color) adobe_coeff (make, model);
    if (load_raw == &CLASS kodak_radc_load_raw)
      if (raw_color) adobe_coeff ("Apple","Quicktake");
-@@ -8725,194 +8635,7 @@
+@@ -8646,7 +8567,7 @@
+   if (!tiff_bps) tiff_bps = 12;
+   if (!maximum) maximum = (1 << tiff_bps) - 1;
+   if (!load_raw || height < 22 || width < 22 ||
+-	tiff_bps > 16 || tiff_samples > 4 || colors > 4)
++	tiff_samples > 4 || colors > 4)
+     is_raw = 0;
+ #ifdef NO_JASPER
+   if (load_raw == &CLASS redcine_load_raw) {
+@@ -8725,195 +8646,269 @@
  }
  #endif
  
@@ -1354,7 +1409,8 @@
 -  if (verbose)
 -    fprintf (stderr, raw_color ? _("Building histograms...\n") :
 -	_("Converting to %s colorspace...\n"), name[output_color-1]);
--
++/* RT: DNG Float */
+ 
 -  memset (histogram, 0, sizeof histogram);
 -  for (img=image[0], row=0; row < height; row++)
 -    for (col=0; col < width; col++, img+=4) {
@@ -1382,7 +1438,9 @@
 -  float r, c, fr, fc;
 -  unsigned ur, uc;
 -  ushort wide, high, (*img)[4], (*pix)[4];
--
++#include <zlib.h>
++#include <stdint.h>
+ 
 -  if (!fuji_width) return;
 -  if (verbose)
 -    fprintf (stderr,_("Rotating image 45 degrees...\n"));
@@ -1405,14 +1463,74 @@
 -	img[row*wide+col][i] =
 -	  (pix[    0][i]*(1-fc) + pix[      1][i]*fc) * (1-fr) +
 -	  (pix[width][i]*(1-fc) + pix[width+1][i]*fc) * fr;
--    }
++static void decodeFPDeltaRow(Bytef * src, Bytef * dst, size_t tileWidth, size_t realTileWidth, int bytesps, int factor) {
++  // DecodeDeltaBytes
++  for (size_t col = factor; col < realTileWidth*bytesps; ++col) {
++    src[col] += src[col - factor];
++  }
++  // Reorder bytes into the image
++  // 16 and 32-bit versions depend on local architecture, 24-bit does not
++  if (bytesps == 3) {
++    for (size_t col = 0; col < tileWidth; ++col) {
++      dst[col*3] = src[col];
++      dst[col*3 + 1] = src[col + realTileWidth];
++      dst[col*3 + 2] = src[col + realTileWidth*2];
+     }
 -  free (image);
 -  width  = wide;
 -  height = high;
 -  image  = img;
 -  fuji_width = 0;
--}
--
++  } else {
++    if (((union { uint32_t x; uint8_t c; }){1}).c) {
++		for (size_t col = 0; col < tileWidth; ++col) {
++			for (size_t byte = 0; byte < bytesps; ++byte)
++				dst[col*bytesps + byte] = src[col + realTileWidth*(bytesps-byte-1)];  // Little endian
++		} 
++    } else {
++		for (size_t col = 0; col < tileWidth; ++col) {
++			for (size_t byte = 0; byte < bytesps; ++byte)
++				dst[col*bytesps + byte] = src[col + realTileWidth*byte];
++        }
++      }
++    }
++  
++}
++
++// From DNG SDK dng_utils.h
++static inline uint32_t DNG_HalfToFloat(uint16_t halfValue) {
++  int32_t sign     = (halfValue >> 15) & 0x00000001;
++  int32_t exponent = (halfValue >> 10) & 0x0000001f;
++  int32_t mantissa =  halfValue        & 0x000003ff;
++  if (exponent == 0) {
++    if (mantissa == 0) {
++      // Plus or minus zero
++      return (uint32_t) (sign << 31);
++    } else {
++      // Denormalized number -- renormalize it
++      while (!(mantissa & 0x00000400)) {
++        mantissa <<= 1;
++        exponent -=  1;
++      }
++      exponent += 1;
++      mantissa &= ~0x00000400;
++    }
++  } else if (exponent == 31) {
++    if (mantissa == 0) {
++      // Positive or negative infinity, convert to maximum (16 bit) values.
++      return (uint32_t) ((sign << 31) | ((0x1eL + 127 - 15) << 23) |  (0x3ffL << 13));
++    } else {
++      // Nan -- Just set to zero.
++      return 0;
++    }
++  }
++  // Normalized number
++  exponent += (127 - 15);
++  mantissa <<= 13;
++  // Assemble sign, exponent and mantissa.
++  return (uint32_t) ((sign << 31) | (exponent << 23) | mantissa);
+ }
+ 
 -void CLASS stretch()
 -{
 -  ushort newdim, (*img)[4], *pix0, *pix1;
@@ -1431,9 +1549,63 @@
 -      if (c+1 < height) pix1 += width*4;
 -      for (col=0; col < width; col++, pix0+=4, pix1+=4)
 -	FORCC img[row*width+col][c] = pix0[c]*(1-frac) + pix1[c]*frac + 0.5;
--    }
++static inline uint32_t DNG_FP24ToFloat(const uint8_t * input) {
++  int32_t sign     = (input [0] >> 7) & 0x01;
++  int32_t exponent = (input [0]     ) & 0x7F;
++  int32_t mantissa = (((int32_t) input [1]) << 8) | input[2];
++  if (exponent == 0) {
++    if (mantissa == 0) {
++      // Plus or minus zero
++      return (uint32_t) (sign << 31);
++    } else {
++      // Denormalized number -- renormalize it
++      while (!(mantissa & 0x00010000)) {
++        mantissa <<= 1;
++        exponent -=  1;
++      }
++      exponent += 1;
++      mantissa &= ~0x00010000;
++    }
++  } else if (exponent == 127) {
++    if (mantissa == 0) {
++      // Positive or negative infinity, convert to maximum (24 bit) values.
++      return (uint32_t) ((sign << 31) | ((0x7eL + 128 - 64) << 23) |  (0xffffL << 7));
++    } else {
++      // Nan -- Just set to zero.
++      return 0;
++    }
++  }
++  // Normalized number
++  exponent += (128 - 64);
++  mantissa <<= 7;
++  // Assemble sign, exponent and mantissa.
++  return (uint32_t) ((sign << 31) | (exponent << 23) | mantissa);
++}
++
++static void expandFloats(Bytef * dst, int tileWidth, int bytesps, float &vmax) {
++  float * dstfl = (float *) dst;
++
++  if (bytesps == 2) {
++    uint16_t * dst16 = (uint16_t *) dst;
++    uint32_t * dst32 = (uint32_t *) dst;
++	
++    for (int index = tileWidth - 1; index >= 0; --index) {
++      dst32[index] = DNG_HalfToFloat(dst16[index]);
++      float v = dstfl[index];
++      if(v>vmax)
++		vmax = v;
++    }
++  } else if (bytesps == 3) {
++    uint8_t  * dst8  = ((uint8_t *) dst) + (tileWidth - 1) * 3;
++    uint32_t * dst32 = (uint32_t *) dst;
++    for (int index = tileWidth - 1; index >= 0; --index, dst8 -= 3) {
++      dst32[index] = DNG_FP24ToFloat(dst8);
++      float v = dstfl[index];
++      if(v>vmax)
++		vmax = v;
+     }
 -    height = newdim;
--  } else {
+   } else {
 -    newdim = width * pixel_aspect + 0.5;
 -    img = (ushort (*)[4]) calloc (height, newdim*sizeof *img);
 -    merror (img, "stretch()");
@@ -1443,25 +1615,157 @@
 -      if (c+1 < width) pix1 += 4;
 -      for (row=0; row < height; row++, pix0+=width*4, pix1+=width*4)
 -	FORCC img[row*newdim+col][c] = pix0[c]*(1-frac) + pix1[c]*frac + 0.5;
--    }
++    for (int index = 0; index < tileWidth; ++index) {
++      float v = dstfl[index];
++      if(v>vmax)
++		vmax = v;
+     }
 -    width = newdim;
--  }
+   }
 -  free (image);
 -  image = img;
--}
--
+ }
+ 
 -int CLASS flip_index (int row, int col)
--{
++static void copyFloatDataToInt(float * src, ushort * dst, size_t size, float max) {
++
++  bool negative = false;
++
++#ifdef _OPENMP
++#pragma omp parallel for
++#endif
++  for (size_t i = 0; i < size; ++i) {
++    // Scale and clip data, because it is used to index some LUTs
++    src[i] *= 65535.0f / max;
++    if (src[i] < 0.0f) {
++      negative = true;
++      src[i] = 0.0f;
++    }
++    // Copy the data to the integer buffer to build the thumbnail
++    dst[i] = (ushort)src[i];
++  }
++  if (negative)
++    fprintf(stderr, "DNG Float: Negative data found in input file\n");
++}
++
++void CLASS deflate_dng_load_raw() {
++  struct tiff_ifd * ifd = &tiff_ifd[0];
++  while (ifd < &tiff_ifd[tiff_nifds] && ifd->offset != data_offset) ++ifd;
++  if (ifd == &tiff_ifd[tiff_nifds]) {
++    fprintf(stderr, "DNG Deflate: Raw image not found???\n");
++    return;
++  }
++
++  int predFactor;
++  switch(ifd->predictor) {
++    case 3: predFactor = 1; break;
++    case 34894: predFactor = 2; break;
++    case 34895: predFactor = 4; break;
++	default: predFactor = 0; break;
++  }
++
++  if (ifd->sample_format == 3) {  // Floating point data
++    float_raw_image = new float[raw_width * raw_height];
++
++#ifdef _OPENMP
++#pragma omp parallel for
++#endif
++    for (size_t i = 0; i < raw_width * raw_height; ++i)
++      float_raw_image[i] = 0.0f;
++    maximum = 65535; black = 0;
++  }
++
++  // NOTE: This reader is based on the official DNG SDK from Adobe.
++  // It assumes tiles without subtiles, but the standard does not say that
++  // subtiles or strips couldn't be used.
++  float maxValue = 0.0f;
++  if (tile_length < INT_MAX) {
++    size_t tilesWide = (raw_width + tile_width - 1) / tile_width;
++    size_t tilesHigh = (raw_height + tile_length - 1) / tile_length;
++    size_t tileCount = tilesWide * tilesHigh;
++    //fprintf(stderr, "%dx%d tiles, %d total\n", tilesWide, tilesHigh, tileCount);
++    size_t tileOffsets[tileCount];
++    for (size_t t = 0; t < tileCount; ++t) {
++      tileOffsets[t] = get4();
++    }
++    size_t tileBytes[tileCount];
++    uLongf maxCompressed = 0;
++    if (tileCount == 1) {
++      tileBytes[0] = maxCompressed = ifd->bytes;
++    } else {
++      fseek(ifp, ifd->bytes, SEEK_SET);
++      for (size_t t = 0; t < tileCount; ++t) {
++        tileBytes[t] = get4();
++        //fprintf(stderr, "Tile %d at %d, size %d\n", t, tileOffsets[t], tileBytes[t]);
++        if (maxCompressed < tileBytes[t])
++          maxCompressed = tileBytes[t];
++      }
++    }
++    uLongf dstLen = tile_width * tile_length * 4;
++
++#ifdef _OPENMP
++#pragma omp parallel
++#endif
++{
++    Bytef cBuffer[maxCompressed];
++    Bytef uBuffer[dstLen];
++    float maxValuethr = 0.0f;
++
++#ifdef _OPENMP
++#pragma omp for collapse(2) nowait
++#endif
++    for (size_t y = 0; y < raw_height; y += tile_length) {
++      for (size_t x = 0; x < raw_width; x += tile_width) {
++		size_t t = (y / tile_length) * tilesWide + (x / tile_width);
++#pragma omp critical
++{
++        fseek(ifp, tileOffsets[t], SEEK_SET);
++        fread(cBuffer, 1, tileBytes[t], ifp);
++}
++        int err = uncompress(uBuffer, &dstLen, cBuffer, tileBytes[t]);
++        if (err != Z_OK) {
++          fprintf(stderr, "DNG Deflate: Failed uncompressing tile %d, with error %d\n", t, err);
++        } else if (ifd->sample_format == 3) {  // Floating point data
++          int bytesps = ifd->bps >> 3;
++          size_t thisTileLength = y + tile_length > raw_height ? raw_height - y : tile_length;
++          size_t thisTileWidth = x + tile_width > raw_width ? raw_width - x : tile_width;
++          for (size_t row = 0; row < thisTileLength; ++row) {
++            Bytef * src = uBuffer + row*tile_width*bytesps;
++            Bytef * dst = (Bytef *)&float_raw_image[(y+row)*raw_width + x];
++            if (predFactor)
++              decodeFPDeltaRow(src, dst, thisTileWidth, tile_width, bytesps, predFactor);
++            expandFloats(dst, thisTileWidth, bytesps, maxValuethr);
++          }
++        } else {  // 32-bit Integer data
++          // TODO
++        }
++      }
++    }
++#ifdef _OPENMP
++#pragma omp critical
++#endif
+ {
 -  if (flip & 4) SWAP(row,col);
 -  if (flip & 2) row = iheight - 1 - row;
 -  if (flip & 1) col = iwidth  - 1 - col;
 -  return row * iwidth + col;
--}
-+/* RT: removed unused functions */
++	if(maxValuethr > maxValue)
++		maxValue = maxValuethr;
++}
++}
++  }
++    
++  if (ifd->sample_format == 3) {  // Floating point data
++    copyFloatDataToInt(float_raw_image, raw_image, raw_width*raw_height, maxValue);
++  }
+ }
  
++/* RT: removed unused functions */
++
  struct tiff_tag {
    ushort tag, type;
-@@ -8935,585 +8658,12 @@
+   int count;
+@@ -8935,585 +8930,12 @@
    unsigned gps[26];
    char desc[512], make[64], model[64], soft[32], date[20], artist[64];
  };