Add Panasonic RW2 embedded lens distortion correction

src/common/exif.cc

@@ -1140,6 +1140,64 @@ static gboolean _check_lens_correction_data(Exiv2::ExifData &exifData,
     }
   }
 
+  /*
+   * Panasonic distortion correction (RW2/RWL)
+   *
+   * The DistortionInfo tag (0x0119 in IFD0 of the RW2) is a 32-byte blob of
+   * 16 signed 16-bit little-endian entries. Exiv2 exposes it as an Undefined
+   * byte array under the numeric key; we parse the fields ourselves.
+   *
+   * Layout (from ExifTool PanasonicRaw.pm + trou/panasonic-rw2 notes):
+   *   [0,1]  checksums (ignored)
+   *   [2]    DistortionParam02 (unknown; leave for future investigation)
+   *   [3]    usually 0
+   *   [4]    DistortionParam04 -> b (r^5)
+   *   [5]    DistortionScale
+   *   [6]    unused
+   *   [7]    DistortionCorrection flag: low nibble 0=off, 1=on
+   *   [8]    DistortionParam08 -> a (r^3)
+   *   [9]    DistortionParam09 (unknown)
+   *   [10]   unused
+   *   [11]   DistortionParam11 -> c (r^7)
+   *   [12]   DistortionN (constant 2500; disables correction if changed)
+   *   [13..] unused / trailing
+   *
+   * Formula (Rigo, https://github.com/trou/panasonic-rw2/blob/master/notes.txt):
+   *   Ru = Rd + scale * (a*Rd^3 + b*Rd^5 + c*Rd^7)  (Rd, Ru normalised to half-diag)
+   *   scale = 1 / (1 + DistortionScale/32768)
+   *   a = DistortionParam08 / 32768
+   *   b = DistortionParam04 / 32768
+   *   c = DistortionParam11 / 32768
+   * we store the raw a/b/c and scale here; _init_coeffs_md_v2 applies scale
+   * when evaluating the polynomial so the two factors stay separable
+   */
+  if((_exif_read_exif_tag(exifData, &pos, "Exif.PanasonicRaw.0x0119")
+      // TIFF/DNG round-trip: dt_exif_read_blob copies the blob here
+      || _exif_read_exif_tag(exifData, &pos, "Exif.Image.0xf119"))
+     && pos->size() == 32)
+  {
+    uint8_t buf[32];
+    pos->copy(buf, Exiv2::littleEndian);
+    int16_t v[16];
+    memcpy(v, buf, sizeof(v));
+
+    const int enabled = (v[7] & 0x0f) == 1;
+    if(enabled)
+    {
+      img->exif_correction_type = CORRECTION_TYPE_PANASONIC;
+      img->exif_correction_data.panasonic.scale = 1.0f / (1.0f + (float)v[5] / 32768.0f);
+      img->exif_correction_data.panasonic.a = (float)v[8]  / 32768.0f;
+      img->exif_correction_data.panasonic.b = (float)v[4]  / 32768.0f;
+      img->exif_correction_data.panasonic.c = (float)v[11] / 32768.0f;
+      dt_print(DT_DEBUG_IMAGEIO,
+               "[exif] Panasonic distortion: scale=%.6f a=%.6f b=%.6f c=%.6f",
+               img->exif_correction_data.panasonic.scale,
+               img->exif_correction_data.panasonic.a,
+               img->exif_correction_data.panasonic.b,
+               img->exif_correction_data.panasonic.c);
+    }
+  }
+
   return img->exif_correction_type != CORRECTION_TYPE_NONE;
 }
 
@@ -2799,6 +2857,23 @@ int dt_exif_read_blob(uint8_t **buf,
       }
     }
 
+    // copy Panasonic DistortionInfo to IFD0 by numeric ID; the
+    // Exif.PanasonicRaw namespace is dropped by exiv2 on TIFF/DNG write.
+    // 0xF119 is in TIFF's private range (>= 0x8000) so no standard collision
+    try
+    {
+      Exiv2::ExifData::iterator src =
+        exifData.findKey(Exiv2::ExifKey("Exif.PanasonicRaw.0x0119"));
+      if(src != exifData.end())
+        exifData.add(Exiv2::ExifKey("Exif.Image.0xf119"), &src->value());
+    }
+    catch(const Exiv2::AnyError &e)
+    {
+      dt_print(DT_DEBUG_IMAGEIO,
+               "[exiv2 dt_exif_read_blob] failed to copy Panasonic distortion: %s",
+               e.what());
+    }
+
     {
       static const char *keys[] = {
         // Canon color space info

src/common/image.h

@@ -155,7 +155,8 @@ typedef enum dt_image_correction_type_t
   CORRECTION_TYPE_SONY,
   CORRECTION_TYPE_FUJI,
   CORRECTION_TYPE_DNG,
-  CORRECTION_TYPE_OLYMPUS
+  CORRECTION_TYPE_OLYMPUS,
+  CORRECTION_TYPE_PANASONIC
 } dt_image_correction_type_t;
 
 typedef union dt_image_correction_data_t
@@ -184,6 +185,14 @@ typedef union dt_image_correction_data_t
     gboolean has_ca;
     float ca[6];
   } olympus;
+  struct {
+    // raw polynomial coefficients from Exif.PanasonicRaw.0x0119, decoded as
+    // word/32768. applied per Rigo's formulation in _init_coeffs_md_v2:
+    //   Ru = Rd + scale * (a*Rd^3 + b*Rd^5 + c*Rd^7)
+    // (Rd, Ru normalised to half-diagonal)
+    float a, b, c;
+    float scale;
+  } panasonic;
 } dt_image_correction_data_t;
 
 typedef enum dt_image_loader_t

src/iop/lens.cc

@@ -2417,6 +2417,52 @@ static int _init_coeffs_md_v2(const dt_image_t *img,
         vig[i] = 1;
     }
   }
+  else if(img->exif_correction_type == CORRECTION_TYPE_PANASONIC)
+  {
+    // Panasonic distortion polynomial (Rigo, trou/panasonic-rw2):
+    //   Ru = Rd + scale * (a*Rd^3 + b*Rd^5 + c*Rd^7)
+    // where Rd is the source (distorted) radius, Ru is the destination
+    // (undistorted) radius, both normalised to the half-diagonal, and
+    // scale is DistortionScale from Panasonic firmware. The pipeline
+    // needs the multiplier dr = Rd/Ru at each destination radius r;
+    // invert with two fixed-point iterations (Rd_new = Ru / f(Rd_old)),
+    // which reduce error by ~1000x for the small distortions typical
+    // of Panasonic bodies
+    const float a  = cd->panasonic.a;
+    const float b  = cd->panasonic.b;
+    const float c  = cd->panasonic.c;
+    const float sc = cd->panasonic.scale;
+
+    nc = MAXKNOTS;
+
+    for(int i = 0; i < nc; i++)
+    {
+      const float r = (float)i / (float)(nc - 1);
+      knots_dist[i] = knots_vig[i] = r;
+
+      if(cor_rgb && p->modify_flags & DT_IOP_LENS_MODIFY_FLAG_DISTORTION)
+      {
+        // invert Ru -> Rd via two fixed-point iterations
+        float rd = r;
+        for(int k = 0; k < 2; k++)
+        {
+          const float rd2 = rd * rd;
+          const float rd4 = rd2 * rd2;
+          const float rd6 = rd4 * rd2;
+          const float f = 1.0f + sc * (a*rd2 + b*rd4 + c*rd6);
+          rd = (f > 1e-6f) ? r / f : r;
+        }
+        const float dr = (r > 0.0f) ? rd / r : 1.0f;
+        const float fine = p->cor_dist_ft * (dr - 1.0f) + 1.0f;
+        cor_rgb[0][i] = cor_rgb[1][i] = cor_rgb[2][i] = fine;
+      }
+      else if(cor_rgb)
+        cor_rgb[0][i] = cor_rgb[1][i] = cor_rgb[2][i] = 1.0f;
+
+      if(vig)
+        vig[i] = 1.0f;
+    }
+  }
 
 
   // calculate the optimal scaling value to show the maximum