Merge pull request #20811 from jenshannoschwalm/various_fixes_and_maintenance

Various fixes and maintenance

Author: TurboGit

data/kernels/basic.cl

@@ -3264,7 +3264,7 @@ colorzones_v3 (read_only image2d_t in,
     default:
     case DT_IOP_COLORZONES_h:
       select = h;
-      blend = dtcl_pow(1.0f - C/128.0f, 2.0f);
+      blend = fsquare(1.0f - C/128.0f);
       break;
   }
 

data/kernels/capture.cl

@@ -205,7 +205,7 @@ __kernel void prepare_blend(__read_only image2d_t cfa,
   if(row > 1 && col > 1 && (row < height-2) && (col < w -2))
   {
     const int w2 = 2 * w;
-    const int color = (filters == 9u) ? FCxtrans(row, col, xtrans) : FC(row, col, filters);
+    const int color = fcol(row, col, filters, xtrans);
     const float val = Areadsingle(cfa, col, row);
     if(val > whites[color] || Y < CAPTURE_YMIN)
     {

data/kernels/colorreconstruction.cl

@@ -93,7 +93,7 @@ colorreconstruction_splat(
   switch(precedence)
   {
     case COLORRECONSTRUCT_PRECEDENCE_CHROMA:
-      weight = sqrt(pixel.y * pixel.y + pixel.z * pixel.z);
+      weight = dt_fast_hypot(pixel.y, pixel.z);
       break;
 
     case COLORRECONSTRUCT_PRECEDENCE_HUE:

data/kernels/demosaic_rcd.cl

@@ -57,35 +57,62 @@ __kernel void rcd_write_output (__write_only image2d_t out, global float *rgb0,
 #define eps 1e-5f              // Tolerance to avoid dividing by zero
 #define epssq 1e-10f
 
-// Step 1.1: Calculate a squared vertical and horizontal high pass filter on color differences
-__kernel void rcd_step_1_1 (global float *cfa, global float *v_diff, global float *h_diff, const int w, const int height)
+static inline float rcd_vdiff_local(local const float *buf, const int stride)
 {
-  const int col = 3 + get_global_id(0);
-  const int row = 3 + get_global_id(1);
-  if((row > height - 4) || (col > w - 4)) return;
-  const int idx = mad24(row, w, col);
-  const int w2 = 2 * w;
-  const int w3 = 3 * w;
+  return fsquare(buf[-3 * stride] - buf[-stride] - buf[stride] + buf[3 * stride] - 3.0f *(buf[-2 * stride] + buf[2 * stride]) + 6.0f * buf[0]);
+}
 
-  v_diff[idx] = fsquare(cfa[idx - w3] - 3.0f * cfa[idx - w2] - cfa[idx - w] + 6.0f * cfa[idx] - cfa[idx + w] - 3.0f * cfa[idx + w2] + cfa[idx + w3]);
-  h_diff[idx] = fsquare(cfa[idx -  3] - 3.0f * cfa[idx -  2] - cfa[idx - 1] + 6.0f * cfa[idx] - cfa[idx + 1] - 3.0f * cfa[idx +  2] + cfa[idx +  3]);
+static inline float rcd_hdiff_local(local const float *buf)
+{
+  return fsquare(buf[-3] - buf[-1] - buf[1] + buf[3] - 3.0f *(buf[-2] + buf[2]) + 6.0f * buf[0]);
 }
 
-// Step 1.2: Calculate vertical and horizontal local discrimination
-__kernel void rcd_step_1_2 (global float *VH_dir, global float *v_diff, global float *h_diff, const int w, const int height)
+// Step 1.1 + 1.2: preload one CFA tile and derive the directional discrimination locally
+// so we avoid materializing two full-frame high-pass buffers in global memory.
+// helpers and rcd_step_1 from ansel code @aurelienpierre
+__kernel void rcd_step_1(global float *cfa, global float *VH_dir, const int w, const int height, local float *buffer)
 {
+  const int xlsz = get_local_size(0);
+  const int ylsz = get_local_size(1);
+  const int xlid = get_local_id(0);
+  const int ylid = get_local_id(1);
+  const int xgid = get_group_id(0);
+  const int ygid = get_group_id(1);
+  const int l = mad24(ylid, xlsz, xlid);
+  const int lsz = mul24(xlsz, ylsz);
+  const int stride = xlsz + 8;
+  const int maxbuf = mul24(stride, ylsz + 8);
+  const int xul = mul24(xgid, xlsz) - 2;
+  const int yul = mul24(ygid, ylsz) - 2;
+
+  for(int n = 0; n <= maxbuf / lsz; n++)
+  {
+    const int bufidx = mad24(n, lsz, l);
+    if(bufidx >= maxbuf) continue;
+    const int xx = clamp(xul + bufidx % stride, 0, w - 1);
+    const int yy = clamp(yul + bufidx / stride, 0, height - 1);
+    buffer[bufidx] = cfa[mad24(yy, w, xx)];
+  }
+
+  barrier(CLK_LOCAL_MEM_FENCE);
+
   const int col = 2 + get_global_id(0);
   const int row = 2 + get_global_id(1);
   if((row > height - 3) || (col > w - 3)) return;
   const int idx = mad24(row, w, col);
-
-  const float V_Stat = fmax(epssq, v_diff[idx - w] + v_diff[idx] + v_diff[idx + w]);
-  const float H_Stat = fmax(epssq, h_diff[idx - 1] + h_diff[idx] + h_diff[idx + 1]);
+  local const float *buf = buffer + mad24(ylid + 4, stride, xlid + 4);
+
+  const float V_Stat = fmax(epssq, rcd_vdiff_local(buf - stride, stride)
+                                  + rcd_vdiff_local(buf, stride)
+                                  + rcd_vdiff_local(buf + stride, stride));
+  const float H_Stat = fmax(epssq, rcd_hdiff_local(buf - 1)
+                                  + rcd_hdiff_local(buf)
+                                  + rcd_hdiff_local(buf + 1));
   VH_dir[idx] = V_Stat / (V_Stat + H_Stat);
 }
 
-// Step 2.1: Low pass filter incorporating green, red and blue local samples from the raw data
-__kernel void rcd_step_2_1(global float *lpf, global float *cfa, const int w, const int height, const unsigned int filters)
+// Step 2: Low pass filter incorporating green, red and blue local samples from the raw data
+__kernel void rcd_step_2(global float *lpf, global float *cfa, const int w, const int height, const unsigned int filters)
 {
   const int row = 2 + get_global_id(1);
   const int col = 2 + (FC(row, 0, filters) & 1) + 2 *get_global_id(0);
@@ -97,8 +124,8 @@ __kernel void rcd_step_2_1(global float *lpf, global float *cfa, const int w, co
     + 0.25f * (cfa[idx - w - 1] + cfa[idx - w + 1] + cfa[idx + w - 1] + cfa[idx + w + 1]);
 }
 
-// Step 3.1: Populate the green channel at blue and red CFA positions
-__kernel void rcd_step_3_1(global float *lpf, global float *cfa, global float *rgb1, global float *VH_Dir, const int w, const int height, const unsigned int filters)
+// Step 3: Populate the green channel at blue and red CFA positions
+__kernel void rcd_step_3(global float *lpf, global float *cfa, global float *rgb1, global float *VH_Dir, const int w, const int height, const unsigned int filters)
 {
   const int row = 4 + get_global_id(1);
   const int col = 4 + (FC(row, 0, filters) & 1) + 2 * get_global_id(0);
@@ -133,11 +160,11 @@ __kernel void rcd_step_3_1(global float *lpf, global float *cfa, global float *r
   const float H_Est = (W_Grad * E_Est + E_Grad * W_Est) / (E_Grad + W_Grad);
 
   // G@B and G@R interpolation
-  rgb1[idx] = mix(V_Est, H_Est, VH_Disc);
+  rgb1[idx] = mix(V_Est, H_Est, clipf(VH_Disc));
 }
 
 // Step 4.0: Calculate the square of the P/Q diagonals color difference high pass filter
-__kernel void rcd_step_4_1(global float *cfa, global float *p_diff, global float *q_diff, const int w, const int height, const unsigned int filters)
+__kernel void rcd_step_4_0(global float *cfa, global float *p_diff, global float *q_diff, const int w, const int height, const unsigned int filters)
 {
   const int row = 3 + get_global_id(1);
   const int col = 3 + 2 * get_global_id(0);
@@ -152,7 +179,7 @@ __kernel void rcd_step_4_1(global float *cfa, global float *p_diff, global float
 }
 
 // Step 4.1: Calculate P/Q diagonals local discrimination strength
-__kernel void rcd_step_4_2(global float *PQ_dir, global float *p_diff, global float *q_diff, const int w, const int height, const unsigned int filters)
+__kernel void rcd_step_4_1(global float *PQ_dir, global float *p_diff, global float *q_diff, const int w, const int height, const unsigned int filters)
 {
   const int row = 2 + get_global_id(1);
   const int col = 2 + (FC(row, 0, filters) & 1) + 2 *get_global_id(0);
@@ -168,7 +195,7 @@ __kernel void rcd_step_4_2(global float *PQ_dir, global float *p_diff, global fl
 }
 
 // Step 4.2: Populate the red and blue channels at blue and red CFA positions
-__kernel void rcd_step_5_1(global float *PQ_dir, global float *rgb0, global float *rgb1, global float *rgb2, const int w, const int height, const unsigned int filters)
+__kernel void rcd_step_4_2(global float *PQ_dir, global float *rgb0, global float *rgb1, global float *rgb2, const int w, const int height, const unsigned int filters)
 {
   const int row = 4 + get_global_id(1);
   const int col = 4 + (FC(row, 0, filters) & 1) + 2 * get_global_id(0);
@@ -204,11 +231,11 @@ __kernel void rcd_step_5_1(global float *PQ_dir, global float *rgb0, global floa
   const float P_Est = (NW_Grad * SE_Est + SE_Grad * NW_Est) / (NW_Grad + SE_Grad);
   const float Q_Est = (NE_Grad * SW_Est + SW_Grad * NE_Est) / (NE_Grad + SW_Grad);
 
-  rgbc[idx]= rgb1[idx] + mix(P_Est, Q_Est, PQ_Disc);
+  rgbc[idx]= rgb1[idx] + mix(P_Est, Q_Est, clipf(PQ_Disc));
 }
 
 // Step 4.3: Populate the red and blue channels at green CFA positions
-__kernel void rcd_step_5_2(global float *VH_dir, global float *rgb0, global float *rgb1, global float *rgb2, const int w, const int height, const unsigned int filters)
+__kernel void rcd_step_4_3(global float *VH_dir, global float *rgb0, global float *rgb1, global float *rgb2, const int w, const int height, const unsigned int filters)
 {
   const int row = 4 + get_global_id(1);
   const int col = 4 + (FC(row, 1, filters) & 1) + 2 * get_global_id(0);
@@ -259,7 +286,7 @@ __kernel void rcd_step_5_2(global float *VH_dir, global float *rgb0, global floa
     const float H_Est = (E_Grad * W_Est + W_Grad * E_Est) / (E_Grad + W_Grad);
 
     // R@G and B@G interpolation
-    rgbc[idx] = rgb1[idx] + mix(V_Est, H_Est, VH_Disc);
+    rgbc[idx] = rgb1[idx] + mix(V_Est, H_Est, clipf(VH_Disc));
   }
 }
 

data/kernels/demosaic_vng.cl

@@ -177,7 +177,7 @@ vng_interpolate(read_only image2d_t in,
     if(bufidx >= maxbuf) continue;
     const int xx = xul + bufidx % stride;
     const int yy = yul + bufidx / stride;
-    const float4 pixel = fmax(0.0f, readpixel(in, xx, yy));
+    const float4 pixel = readpixel(in, xx, yy);
     vstore4(pixel, bufidx, buffer);
   }
 

data/kernels/extended.cl

@@ -263,9 +263,9 @@ vibrance (read_only image2d_t in, write_only image2d_t out, const int width, con
 
   if(x >= width || y >= height) return;
 
-  float4 pixel = read_imagef(in, sampleri, (int2)(x, y));
+  float4 pixel = readpixel(in, x, y);
 
-  const float sw = sqrt(pixel.y*pixel.y + pixel.z*pixel.z)/256.0f;
+  const float sw = dt_fast_hypot(pixel.y, pixel.z)/256.0f;
   const float ls = 1.0f - amount * sw * 0.25f;
   const float ss = 1.0f + amount * sw;
 

src/control/jobs/control_jobs.c

@@ -554,7 +554,7 @@ static int _control_merge_hdr_process(dt_imageio_module_data_t *datai,
   const float eap = image.exif_aperture > 0.0f ? image.exif_aperture : 22.0f;
   const float efl = image.exif_focal_length > 0.0f ? image.exif_focal_length : 8.0f;
   const float rad = .5f * efl / eap;
-  const float aperture = M_PI * rad * rad;
+  const float aperture = M_PI_F * rad * rad;
   const float iso = image.exif_iso > 0.0f ? image.exif_iso : 100.0f;
   const float exp = image.exif_exposure > 0.0f ? image.exif_exposure : 1.0f;
   const float cal = 100.0f / (aperture * exp * iso);

src/develop/blend.c

@@ -252,7 +252,7 @@ static inline float _detail_mask_threshold(const float level,
                                            const gboolean detail)
 {
   // this does some range calculation for smoother ui experience
-  return 0.005f * (detail ? powf(level, 2.0f) : 1.0f - powf(fabs(level), 0.5f ));
+  return 0.005f * (detail ? sqrf(level) : 1.0f - sqrtf(fabs(level)));
 }
 
 static void _refine_with_detail_mask(dt_iop_module_t *self,
@@ -1478,7 +1478,7 @@ void tiling_callback_blendop(dt_iop_module_t *self,
       if(devid > DT_DEVICE_CPU)
       {
         /* OpenCL feathering does simple internal tiling for less mem pressure,
-           we still need some mem here for this. 
+           we still need some mem here for this.
         */
         tiling->factor_cl = MAX(tiling->factor, 1.0f);
       }

src/develop/blends/blendif_lab.c

@@ -1347,7 +1347,7 @@ static void _display_channel(const float *const restrict a,
     }
     case DT_DEV_PIXELPIPE_DISPLAY_LCH_C:
     {
-      const float factor = 1.0f / (128.0f * sqrtf(2.0f) * exp2f(boost_factors[DEVELOP_BLENDIF_C_in]));
+      const float factor = 1.0f / (128.0f * M_SQRT2_F * exp2f(boost_factors[DEVELOP_BLENDIF_C_in]));
       for(size_t i = 0, j = 0; i < stride; i++, j += DT_BLENDIF_LAB_CH)
       {
         dt_aligned_pixel_t LCH;

src/develop/masks/brush.c

@@ -248,7 +248,7 @@ static void _brush_border_get_XY(const float p0x,
     *yb = DT_INVALID_COORDINATE;
     return;
   }
-  const float l = 1.0f / sqrtf(dx * dx + dy * dy);
+  const float l = 1.0f / dt_fast_hypotf(dx, dy);
   *xb = (*xc) + rad * dy * l;
   *yb = (*yc) - rad * dx * l;
 }
@@ -430,18 +430,16 @@ static void _brush_points_recurs_border_gaps(float *cmax,
   // we have to be sure that we turn in the correct direction
   if(a2 < a1 && clockwise)
   {
-    a2 += 2.0f * M_PI;
+    a2 += DT_2PI_F;
   }
   if(a2 > a1 && !clockwise)
   {
-    a1 += 2.0f * M_PI;
+    a1 += DT_2PI_F;
   }
 
   // we determine start and end radius too
-  float r1 = sqrtf((bmin[1] - cmax[1]) * (bmin[1] - cmax[1])
-                   + (bmin[0] - cmax[0]) * (bmin[0] - cmax[0]));
-  float r2 = sqrtf((bmax[1] - cmax[1]) * (bmax[1] - cmax[1])
-                   + (bmax[0] - cmax[0]) * (bmax[0] - cmax[0]));
+  float r1 = dt_fast_hypotf(bmin[1] - cmax[1], bmin[0] - cmax[0]);
+  float r2 = dt_fast_hypotf(bmax[1] - cmax[1], bmax[0] - cmax[0]);
 
   // and the max length of the circle arc
   int l;
@@ -489,21 +487,19 @@ static void _brush_points_recurs_border_small_gaps(float *cmax,
 {
   // we want to find the start and end angles
   const float a1 = fmodf(atan2f(bmin[1] - cmax[1], bmin[0] - cmax[0])
-                         + 2.0f * M_PI, 2.0f * M_PI);
+                         + DT_2PI_F, DT_2PI_F);
   const float a2 = fmodf(atan2f(bmax[1] - cmax[1], bmax[0] - cmax[0])
-                         + 2.0f * M_PI, 2.0f * M_PI);
+                         + DT_2PI_F, DT_2PI_F);
 
   if(a1 == a2) return;
 
   // we determine start and end radius too
-  const float r1 = sqrtf((bmin[1] - cmax[1]) * (bmin[1] - cmax[1])
-                         + (bmin[0] - cmax[0]) * (bmin[0] - cmax[0]));
-  const float r2 = sqrtf((bmax[1] - cmax[1]) * (bmax[1] - cmax[1])
-                         + (bmax[0] - cmax[0]) * (bmax[0] - cmax[0]));
+  const float r1 = dt_fast_hypotf(bmin[1] - cmax[1], bmin[0] - cmax[0]);
+  const float r2 = dt_fast_hypotf(bmax[1] - cmax[1], bmax[0] - cmax[0]);
 
   // we close the gap in the shortest direction
   float delta = a2 - a1;
-  if(fabsf(delta) > M_PI) delta = delta - copysignf(2.0f * M_PI, delta);
+  if(fabsf(delta) > M_PI_F) delta = delta - copysignf(DT_2PI_F, delta);
 
   // get the max length of the circle arc
   const int l = fabsf(delta) * fmaxf(r1, r2);
@@ -547,15 +543,14 @@ static void _brush_points_stamp(float *cmax,
   const float a1 = atan2f(bmin[1] - cmax[1], bmin[0] - cmax[0]);
 
   // we determine the radius too
-  const float rad = sqrtf((bmin[1] - cmax[1]) * (bmin[1] - cmax[1])
-                          + (bmin[0] - cmax[0]) * (bmin[0] - cmax[0]));
+  const float rad = dt_fast_hypotf(bmin[1] - cmax[1], bmin[0] - cmax[0]);
 
   // determine the max length of the circle arc
   const int l = 2.0f * M_PI * rad;
   if(l < 2) return;
 
   // and now we add the points
-  const float incra = 2.0f * M_PI / l;
+  const float incra = DT_2PI_F / l;
   float aa = a1 + incra;
   // allocate entries in the dynbufs
   float *dpoints_ptr = dt_masks_dynbuf_reserve_n(dpoints, 2*(l-1));
@@ -2273,7 +2268,7 @@ static int _brush_events_mouse_moved(struct dt_iop_module_t *module,
     dt_masks_point_brush_t *point = g_list_nth_data(form->points, k);
     const float nx = point->corner[0] * iwidth;
     const float ny = point->corner[1] * iheight;
-    const float nr = sqrtf((pts[0] - nx) * (pts[0] - nx) + (pts[1] - ny) * (pts[1] - ny));
+    const float nr = dt_fast_hypotf(pts[0] - nx, pts[1] - ny);
     const float bdr = nr / fminf(iwidth, iheight);
 
     point->border[0] = point->border[1] = bdr;