Add calculation of capture sharpen radius

1. Add a button beside the radius slider. Pressing that leads to a calculation of a very good guess of the
   radius to be used.
   Visual feedback via control log plus setting of the slider via a ui pipe finished signal.
2. You may also set the radius to zero, in that case we do the radius calculation at runtime

Author: jenshannoschwalm

src/iop/demosaic.c

@@ -148,6 +148,8 @@ typedef struct dt_iop_demosaic_gui_data_t
   GtkWidget *cs_strength;
   gboolean cs_mask;
   gboolean dual_mask;
+  gboolean autoradius;
+  float newradius;
 } dt_iop_demosaic_gui_data_t;
 
 typedef struct dt_iop_demosaic_global_data_t
@@ -1436,6 +1438,9 @@ void gui_update(dt_iop_module_t *self)
 {
   gui_changed(self, NULL, NULL);
   gtk_stack_set_visible_child_name(GTK_STACK(self->widget), self->default_enabled ? "raw" : "non_raw");
+  dt_iop_demosaic_gui_data_t *g = self->gui_data;
+  g->autoradius = FALSE;
+  g->newradius = 0.0f;
 }
 
 static void _dual_quad_callback(GtkWidget *quad, dt_iop_module_t *self)
@@ -1461,6 +1466,28 @@ static void _cs_quad_callback(GtkWidget *quad, dt_iop_module_t *self)
   dt_dev_reprocess_center(self->dev);
 }
 
+static void _cs_autoradius_callback(GtkWidget *quad, dt_iop_module_t *self)
+{
+  if(darktable.gui->reset) return;
+  dt_iop_demosaic_gui_data_t *g = self->gui_data;
+  g->autoradius = TRUE;
+  g->newradius = 2.0f;
+  dt_dev_reprocess_center(self->dev);
+}
+
+static void _check_autoradius(gpointer instance, dt_iop_module_t *self)
+{
+  dt_iop_demosaic_gui_data_t *g = self->gui_data;
+  if(g && g->autoradius)
+  {
+    dt_iop_demosaic_params_t *p = self->params;
+    g->autoradius = FALSE;
+    p->cs_radius = g->newradius;
+    dt_bauhaus_slider_set_val(g->cs_radius, g->newradius);
+    dt_dev_add_history_item(darktable.develop, self, TRUE);
+  }
+}
+
 void gui_focus(dt_iop_module_t *self, gboolean in)
 {
   dt_iop_demosaic_gui_data_t *g = self->gui_data;
@@ -1526,6 +1553,10 @@ void gui_init(dt_iop_module_t *self)
   dt_bauhaus_slider_set_step(g->cs_radius, 0.02f);
   dt_bauhaus_slider_set_format(g->cs_radius, _(_(" px")));
   gtk_widget_set_tooltip_text(g->cs_radius, _("capture sharpen radius. increasing this too far will lead to halos"));
+  dt_bauhaus_widget_set_quad(g->cs_radius, self, dtgtk_cairo_paint_reset, FALSE, _cs_autoradius_callback,
+    _("approximate a safe radius from sensor data"));
+  g->autoradius = FALSE;
+  g->newradius = 0.0f;
 
   g->cs_thrs = dt_bauhaus_slider_from_params(self, "cs_thrs");
   dt_bauhaus_slider_set_format(g->cs_thrs, "%");
@@ -1549,6 +1580,7 @@ void gui_init(dt_iop_module_t *self)
 
   gtk_stack_add_named(GTK_STACK(self->widget), label_non_raw, "non_raw");
   gtk_stack_add_named(GTK_STACK(self->widget), box_raw, "raw");
+  DT_CONTROL_SIGNAL_HANDLE(DT_SIGNAL_DEVELOP_UI_PIPE_FINISHED, _check_autoradius);
 }
 
 // clang-format off

src/iop/demosaicing/capture.c

@@ -335,6 +335,232 @@ static void _build_blend_mask(float *clipmask,
   }
 }
 
+#define RAWEPS 0.0001f
+float calcRadiusBayer(const float *in,
+                      const int width,
+                      const int height,
+                      const float lowerLimit,
+                      const float upperLimit,
+                      const uint32_t filters)
+{
+  const unsigned int fc[2] = {FC(0, 0, filters), FC(1, 0, filters)};
+  float maxRatio = 1.f;
+  DT_OMP_FOR(reduction(max: maxRatio))
+  for(int row = 4; row < height - 4; ++row)
+  {
+    for(int col = 5 + (fc[row & 1] & 1); col < width - 4; col += 2)
+    {
+      const float *cfa = in + row*width + col;
+      const float val00 = cfa[0];
+      if(val00 > RAWEPS)
+      {
+        const float val1m1 = cfa[width-1];
+        const float val1p1 = cfa[width+1];
+        const float maxVal0 = MAX(val00, val1m1);
+        if(val1m1 > RAWEPS && maxVal0 > lowerLimit)
+        {
+          const float minVal = MIN(val00, val1m1);
+          if(maxVal0 > maxRatio * minVal)
+          {
+            gboolean clipped = FALSE;
+            if(maxVal0 == val00)
+            { // check for influence by clipped green in neighborhood
+              if(MAX(MAX(cfa[-width-1], cfa[-width+1]), val1p1) >= upperLimit)
+              {
+                clipped = TRUE;
+              }
+            }
+            else
+            { // check for influence by clipped green in neighborhood
+              if(MAX(MAX(MAX(cfa[-2], val00), cfa[2*width-2]), cfa[2*width]) >= upperLimit)
+              {
+                clipped = TRUE;
+              }
+            }
+            if(!clipped)
+            {
+              maxRatio = maxVal0 / minVal;
+            }
+          }
+        }
+
+        const float maxVal1 = MAX(val00, val1p1);
+        if(val1p1 > RAWEPS && maxVal1 > lowerLimit)
+        {
+          const float minVal = MIN(val00, val1p1);
+          if(maxVal1 > maxRatio * minVal)
+          {
+            if(maxVal1 == val00)
+            { // check for influence by clipped green in neighborhood
+              if(MAX(MAX(cfa[-width-1], cfa[-width+1]), val1p1) >= upperLimit)
+              {
+                continue;
+              }
+            }
+            else
+            { // check for influence by clipped green in neighborhood
+              if(MAX(MAX(MAX(val00, cfa[2]), cfa[2*width]), cfa[2*width+2]) >= upperLimit)
+              {
+                continue;
+              }
+            }
+            maxRatio = maxVal1 / minVal;
+          }
+        }
+      }
+    }
+  }
+  return sqrtf(1.0f / logf(maxRatio));
+}
+
+float calcRadiusXtrans(const float *in,
+                       const float lowerLimit,
+                       const float upperLimit,
+                       const dt_iop_roi_t *const roi,
+                       const uint8_t(*const xtrans)[6])
+{
+  const int width = roi->width;
+  const int height = roi->height;
+
+  int startx, starty;
+  gboolean found = FALSE;
+  for(starty = 6; starty < 12 && !found; starty++)
+  {
+    for(startx = 6; startx < 12 && !found; startx++)
+    {
+      if(FCxtrans(starty, startx, roi, xtrans) == 1)
+      {
+        if(FCxtrans(starty, startx - 1, roi, xtrans) != FCxtrans(starty, startx + 1, roi, xtrans))
+        {
+          if(FCxtrans(starty -1, startx, roi, xtrans) != 1)
+          {
+            if(FCxtrans(starty, startx -1, roi, xtrans) != 1)
+            {
+              found = TRUE;
+              break;
+            }
+          }
+        }
+      }
+    }
+  }
+
+  float maxRatio = 1.0f;
+  DT_OMP_FOR(reduction(max: maxRatio))
+  for(int row = starty + 2; row < height - 4; row += 3)
+  {
+    for(int col = startx + 2; col < width - 4; col += 3)
+    {
+      const float *cfa = in + row*width + col;
+      const float valp1p1 = cfa[width+1];
+      const gboolean squareClipped = MAX(MAX(MAX(valp1p1, cfa[width+2]), cfa[2*width+1]), cfa[2*width+2]) >= upperLimit;
+      const float greenSolitary = cfa[0];
+      if(greenSolitary > RAWEPS && MAX(cfa[-width-1], cfa[-width+1]) < upperLimit)
+      {
+        if(greenSolitary < upperLimit)
+        {
+          const float valp1m1 = cfa[width-1];
+          if(valp1m1 > RAWEPS && MAX(MAX(MAX(cfa[width-2], valp1m1), cfa[2*width-2]), cfa[width-1]) < upperLimit)
+          {
+            const float maxVal = MAX(greenSolitary, valp1m1);
+            if(maxVal > lowerLimit)
+            {
+              const float minVal = MIN(greenSolitary, valp1m1);
+              if(maxVal > maxRatio * minVal)
+              {
+                maxRatio = maxVal / minVal;
+              }
+            }
+          }
+          if(valp1p1 > RAWEPS && !squareClipped)
+          {
+            const float maxVal = MAX(greenSolitary, valp1p1);
+            if(maxVal > lowerLimit)
+            {
+              const float minVal = MIN(greenSolitary, valp1p1);
+              if(maxVal > maxRatio * minVal)
+              {
+                maxRatio = maxVal / minVal;
+              }
+            }
+          }
+        }
+      }
+
+      if(!squareClipped)
+      {
+        const float valp2p2 = cfa[2*width+2];
+        if(valp2p2 > RAWEPS)
+        {
+          if(valp1p1 > RAWEPS)
+          {
+            const float maxVal = MAX(valp1p1, valp2p2);
+            if(maxVal > lowerLimit)
+            {
+              const float minVal = MIN(valp1p1, valp2p2);
+              if(maxVal > maxRatio * minVal)
+              {
+                maxRatio = maxVal / minVal;
+              }
+            }
+          }
+          const float greenSolitaryRight = cfa[3*width+3];
+          if(MAX(MAX(greenSolitaryRight, cfa[4*width+2]), cfa[4*width+4]) < upperLimit)
+          {
+            if(greenSolitaryRight > RAWEPS)
+            {
+              const float maxVal = MAX(greenSolitaryRight, valp2p2);
+              if(maxVal > lowerLimit)
+              {
+                const float minVal = MIN(greenSolitaryRight, valp2p2);
+                if(maxVal > maxRatio * minVal)
+                {
+                  maxRatio = maxVal / minVal;
+                }
+              }
+            }
+          }
+        }
+        const float valp1p2 = cfa[width+2];
+        const float valp2p1 = cfa[2*width+1];
+        if(valp2p1 > RAWEPS)
+        {
+          if(valp1p2 > RAWEPS)
+          {
+            const float maxVal = MAX(valp1p2, valp2p1);
+            if(maxVal > lowerLimit)
+            {
+              const float minVal = MIN(valp1p2, valp2p1);
+              if(maxVal > maxRatio * minVal)
+              {
+                maxRatio = maxVal / minVal;
+              }
+            }
+          }
+          const float greenSolitaryLeft = cfa[3*width];
+          if(MAX(MAX(greenSolitaryLeft, cfa[4*width-1]), cfa[4*width+1]) < upperLimit)
+          {
+            if(greenSolitaryLeft > RAWEPS)
+            {
+              const float maxVal = MAX(greenSolitaryLeft, valp2p1);
+              if(maxVal > lowerLimit)
+              {
+                const float minVal = MIN(greenSolitaryLeft, valp2p1);
+                if(maxVal > maxRatio * minVal)
+                {
+                  maxRatio = maxVal / minVal;
+                }
+              }
+            }
+          }
+        }
+      }
+    }
+  }
+  return sqrtf(1.0f / logf(maxRatio));
+}
+#undef RAWEPS
+
 void _capture_sharpen(dt_iop_module_t *self,
                       dt_dev_pixelpipe_iop_t *piece,
                       float *in,
@@ -350,6 +576,7 @@ void _capture_sharpen(dt_iop_module_t *self,
   const size_t pixels = width * height;
   const dt_iop_demosaic_data_t *d = piece->data;
   const dt_iop_demosaic_global_data_t *gd = self->global_data;
+  dt_iop_demosaic_gui_data_t *g = self->gui_data;
 
   const uint8_t(*const xtrans)[6] = (const uint8_t(*const)[6])piece->pipe->dsc.xtrans;
   const uint32_t filters = piece->pipe->dsc.filters;
@@ -359,6 +586,26 @@ void _capture_sharpen(dt_iop_module_t *self,
                                        wbon ? dsc->temperature.coeffs[1] : 1.0f,
                                        wbon ? dsc->temperature.coeffs[2] : 1.0f,
                                        0.0f };
+  const gboolean autoradius = (pipe->type & DT_DEV_PIXELPIPE_FULL)
+                                && g && g->autoradius;
+  float radius = d->cs_radius;
+  const gboolean noradius = radius < 0.01f;
+  if(autoradius || noradius)
+  {
+    radius = filters != 9u
+              ? calcRadiusBayer(in, width, height, 0.01f, 1.0f, filters)
+              : calcRadiusXtrans(in, 0.01f, 1.0f, roi, xtrans);
+
+    dt_print_pipe(DT_DEBUG_PIPE, filters != 9u ? "bayer autoradius" : "xtrans autoradius",
+      pipe, self, DT_DEVICE_CPU, roi, NULL, "autoradius=%.2f", radius);
+    if(autoradius)
+    {
+      dt_control_log(_("calculated radius: %.2f"), radius);
+      dt_iop_gui_enter_critical_section(self);
+      g->newradius = MIN(g->newradius, radius);
+      dt_iop_gui_leave_critical_section(self);
+    }
+  }
 
   char *gauss_idx = NULL;
   gboolean error = TRUE;
@@ -385,7 +632,7 @@ void _capture_sharpen(dt_iop_module_t *self,
     goto finalize;
   }
 
-  gauss_idx = _cs_precalc_gauss_idx(self, roi, d->cs_radius, d->cs_boost);
+  gauss_idx = _cs_precalc_gauss_idx(self, roi, radius, d->cs_boost);
   if(!gauss_idx) goto finalize;
 
   for(int iter = 0; iter < iterations; iter++)
@@ -434,6 +681,7 @@ int _capture_sharpen_cl(dt_iop_module_t *self,
   const int devid = piece->pipe->devid;
   const dt_iop_demosaic_data_t *d = piece->data;
   dt_iop_demosaic_global_data_t *const gd = self->global_data;
+  dt_iop_demosaic_gui_data_t *g = self->gui_data;
 
   const float threshold = 0.09f * d->cs_thrs;
   const uint32_t filters = piece->pipe->dsc.filters;
@@ -443,6 +691,33 @@ int _capture_sharpen_cl(dt_iop_module_t *self,
                                  wbon ? dsc->temperature.coeffs[1] : 1.0f,
                                  wbon ? dsc->temperature.coeffs[2] : 1.0f,
                                  0.0f };
+  const gboolean autoradius = (pipe->type & DT_DEV_PIXELPIPE_FULL)
+                                && g && g->autoradius;
+  float radius = d->cs_radius;
+  const gboolean noradius = radius < 0.01f;
+  if(autoradius || noradius)
+  {
+    float *in = dt_alloc_align_float((size_t)width * height);
+    if(in)
+    {
+      if(dt_opencl_read_host_from_device(devid, in, dev_in, width, height, sizeof(float)) == CL_SUCCESS)
+      {
+        radius = filters != 9u
+                ? calcRadiusBayer(in, width, height, 0.01f, 1.0f, filters)
+                : calcRadiusXtrans(in, 0.01f, 1.0f, roi, (const uint8_t(*const)[6])piece->pipe->dsc.xtrans);
+        dt_print_pipe(DT_DEBUG_PIPE, filters != 9u ? "bayer autoradius" : "xtrans autoradius",
+            pipe, self, devid, roi, NULL, "autoradius=%.2f", radius);
+        if(autoradius)
+        {
+          dt_control_log(_("calculated radius: %.2f"), radius);
+          dt_iop_gui_enter_critical_section(self);
+          g->newradius = MIN(g->newradius, radius);
+          dt_iop_gui_leave_critical_section(self);
+        }
+      }
+      dt_free_align(in);
+    }
+  }
 
   cl_mem gcoeffs = NULL;
   cl_mem gauss_idx = NULL;
@@ -481,7 +756,7 @@ int _capture_sharpen_cl(dt_iop_module_t *self,
     goto finish;
   }
 
-  char *f_gauss_idx = _cs_precalc_gauss_idx(self, roi, d->cs_radius, d->cs_boost);
+  char *f_gauss_idx = _cs_precalc_gauss_idx(self, roi, radius, d->cs_boost);
   if(f_gauss_idx)
   {
     gcoeffs = dt_opencl_copy_host_to_device_constant(devid, sizeof(float) * CHAR_MAX * CAPTURE_KERNEL_ALIGN, gd->gauss_coeffs);