detail.c

Go to the documentation of this file.

/*
    This file is part of darktable,
    Copyright (C) 2021 Hanno Schwalm.
    Copyright (C) 2021 luzpaz.
    Copyright (C) 2021 Pascal Obry.
    Copyright (C) 2021 Ralf Brown.
    Copyright (C) 2021 Roman Lebedev.
    Copyright (C) 2022 Martin Bařinka.
    Copyright (C) 2023, 2025 Aurélien PIERRE.
    
    darktable is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.
    
    darktable is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.
    
    You should have received a copy of the GNU General Public License
    along with darktable.  If not, see <http://www.gnu.org/licenses/>.
*/
 
/* Declares what this file implements. This file is textually #included by masks.c rather than
 * compiled on its own, so the include has to live HERE and not in the host: clang-tidy attributes
 * a symbol's use to the file that names it, and an include placed in masks.c for detail.c's
 * benefit reads as unused there. */
#include "develop/masks_detail.h"
 
/* How are "detail masks" implemented?
 
  The detail masks (DM) are used by the dual demosaicer and as a further refinement step for
  shape / parametric masks.
  They contain threshold weighed values of pixel-wise local signal changes so they can be
  understood as "areas with or without local detail".
 
  As the DM using algorithms (like dual demosaicing, sharpening ...) are all pixel peeping we
  want the "original data" from the sensor to calculate it.
  (Calculating the mask from the modules roi might not detect such regions at all because of
  scaling / rotating artifacts, some blurring earlier in the pipeline, color changes ...)
 
  In all cases the user interface is pretty simple, we just pass a threshold value, which
  is in the range of -1.0 to 1.0 by an additional slider in the masks refinement section.
  Positive values will select regions with lots of local detail, negatives select for flat areas.
  (The dual demosaicer only wants positives as we always look for high frequency content.)
  A threshold value of 0.0 means bypassing.
 
  So the first important point is:
  We make sure taking the input data for the DM from a dedicated hidden pipeline stage placed
  right after demosaic. This means some additional housekeeping for the pixelpipe.
  If any mask in any module selects a threshold of != 0.0 we leave a flag in the pipe struct
  telling that we want a DM from that dedicated stage. If such a flag has not been previously
  set we will force a pipeline reprocessing.
 
  The hidden `detailmask` module writes a preliminary mask holding signal-change values for every pixel
  in its CPU and OpenCL `process()` callbacks.
  These mask values are calculated as
  a) get Y0 for every pixel
  b) apply a scharr operator on it
 
  This raw detail mask (RM) is not scaled but only cropped to the roi of the writing module.
  The pipe gets roi copy of the writing module so we can later scale/distort the LM.
 
  Calculating the RM is done for performance and lower mem pressure reasons, so we don't have to
  pass full data to the module. Also the RM can be used by other modules.
 
  If a mask uses the details refinement step it takes the raw details mask RM and calculates an
  intermediate mask (IM) which is still not scaled but has the roi of the writing module.
 
  For every pixel we calculate the IM value via a sigmoid function with the threshold and RM as parameters.
 
  At last the IM is slightly blurred to avoid hard transitions, as there still is no scaling we can use
  a constant sigma. As the blur_9x9 is pretty fast both in openmp/cl code paths - much faster than dt
  gaussians - it is used here.
  Now we have an unscaled detail mask which requires to be transformed through the pipeline using
 
  float *dt_dev_distort_detail_mask(const dt_dev_pixelpipe_t *pipe, float *src, const dt_iop_module_t *target_module)
 
  returning a pointer to a distorted mask (DT) with same size as used in the module wanting the refinement.
  This DM is finally used to refine the original mask.
 
  All other refinements and parametric parameters are untouched.
 
  Some additional comments:
  1. the detail mask is authored from the RGBA float pipeline stage immediately after demosaic,
     which keeps the source buffer format stable across RAW and non-RAW inputs.
  2. In the gui the slider is above the rest of the refinemt sliders to emphasize that blurring & feathering use the
     mask corrected by detail refinemnt.
  3. Of course credit goes to Ingo @heckflosse from rt team for the original idea. (in the rt world this is knowb
     as details mask)
  4. Thanks to rawfiner for pointing out how to use Y0 and scharr for better maths.
 
  hanno@schwalm-brmouseemen.de 21/04/29
*/
 
__DT_CLONE_TARGETS__
void dt_masks_extend_border(float *const restrict mask, const int width, const int height, const int border)
{
  if(border <= 0) return;
  const int max_col = width - border - 1;
  __OMP_PARALLEL_FOR_SIMD__(aligned(mask : 64) if((size_t)width * height > 10000))
  for(int row = border; row < height - border; row++)
  {
    float *const rowptr = mask + (size_t)(row * width);
    for(int i = 0; i < border; i++)
    {
      rowptr[i] = rowptr[border];
      rowptr[width - i - 1] = rowptr[max_col];
    }
  }
  const float *const top_row = mask + (size_t)(border * width);
  const float *const bot_row = mask + (size_t)(height - border - 1) * width;
  __OMP_FOR_SIMD__(aligned(mask : 64) if((size_t)width * height > 10000))
  for(int col = 0; col < width; col++)
  {
    const int c = MIN(max_col, MAX(col, border));
    const float top = top_row[c];
    const float bot = bot_row[c];
    for(int i = 0; i < border; i++)
    {
      mask[col + i * width] = top;
      mask[col + (height - i - 1) * width] = bot;
    }
  }
}
 
void _masks_blur_5x5_coeff(float *c, const float sigma)
{
  float kernel[5][5];
  const float temp = -2.0f * sqf(sigma);
  const float range = sqf(3.0f * 0.84f);
  float sum = 0.0f;
  for(int k = -2; k <= 2; k++)
  {
    for(int j = -2; j <= 2; j++)
    {
      if((sqf(k) + sqf(j)) <= range)
      {
        kernel[k + 2][j + 2] = expf((sqf(k) + sqf(j)) / temp);
        sum += kernel[k + 2][j + 2];
      }
      else
        kernel[k + 2][j + 2] = 0.0f;
    }
  }
  for(int i = 0; i < 5; i++)
    for(int j = 0; j < 5; j++)
      kernel[i][j] /= sum;
 
  // FIXME: are you for real ? managing arrays with loops and index shift much ?
  /* c21 */ c[0]  = kernel[0][1];
  /* c20 */ c[1]  = kernel[0][2];
  /* c11 */ c[2]  = kernel[1][1];
  /* c10 */ c[3]  = kernel[1][2];
  /* c00 */ c[4]  = kernel[2][2];
}
#define FAST_BLUR_5 ( \
  blurmat[0] * ((src[i - w2 - 1] + src[i - w2 + 1]) + (src[i - w1 - 2] + src[i - w1 + 2]) + (src[i + w1 - 2] + src[i + w1 + 2]) + (src[i + w2 - 1] + src[i + w2 + 1])) + \
  blurmat[1] * (src[i - w2] + src[i - 2] + src[i + 2] + src[i + w2]) + \
  blurmat[2] * (src[i - w1 - 1] + src[i - w1 + 1] + src[i + w1 - 1] + src[i + w1 + 1]) + \
  blurmat[3] * (src[i - w1] + src[i - 1] + src[i + 1] + src[i + w1]) + \
  blurmat[4] * src[i] )
 
void dt_masks_blur_9x9_coeff(float *c, const float sigma)
{
  float kernel[9][9];
  const float temp = -2.0f * sqf(sigma);
  const float range = sqf(3.0f * 1.5f);
  float sum = 0.0f;
  for(int k = -4; k <= 4; k++)
  {
    for(int j = -4; j <= 4; j++)
    {
      if((sqf(k) + sqf(j)) <= range)
      {
        kernel[k + 4][j + 4] = expf((sqf(k) + sqf(j)) / temp);
        sum += kernel[k + 4][j + 4];
      }
      else
        kernel[k + 4][j + 4] = 0.0f;
    }
  }
  for(int i = 0; i < 9; i++)
    for(int j = 0; j < 9; j++)
      kernel[i][j] /= sum;
 
  // FIXME: are you for real ? managing arrays with loops and index shift much ?
  /* c00 */ c[0]  = kernel[4][4];
  /* c10 */ c[1]  = kernel[3][4];
  /* c11 */ c[2]  = kernel[3][3];
  /* c20 */ c[3]  = kernel[2][4];
  /* c21 */ c[4]  = kernel[2][3];
  /* c22 */ c[5]  = kernel[2][2];
  /* c30 */ c[6]  = kernel[1][4];
  /* c31 */ c[7]  = kernel[1][3];
  /* c32 */ c[8]  = kernel[1][2];
  /* c33 */ c[9]  = kernel[1][1];
  /* c40 */ c[10] = kernel[0][4];
  /* c41 */ c[11] = kernel[0][3];
  /* c42 */ c[12] = kernel[0][2];
}
 
// FIMXE: ever heard about loop unrolling ???
#define FAST_BLUR_9 ( \
  blurmat[12] * (src[i - w4 - 2] + src[i - w4 + 2] + src[i - w2 - 4] + src[i - w2 + 4] + src[i + w2 - 4] + src[i + w2 + 4] + src[i + w4 - 2] + src[i + w4 + 2]) + \
  blurmat[11] * (src[i - w4 - 1] + src[i - w4 + 1] + src[i - w1 - 4] + src[i - w1 + 4] + src[i + w1 - 4] + src[i + w1 + 4] + src[i + w4 - 1] + src[i + w4 + 1]) + \
  blurmat[10] * (src[i - w4] + src[i - 4] + src[i + 4] + src[i + w4]) + \
  blurmat[9]  * (src[i - w3 - 3] + src[i - w3 + 3] + src[i + w3 - 3] + src[i + w3 + 3]) + \
  blurmat[8]  * (src[i - w3 - 2] + src[i - w3 + 2] + src[i - w2 - 3] + src[i - w2 + 3] + src[i + w2 - 3] + src[i + w2 + 3] + src[i + w3 - 2] + src[i + w3 + 2]) + \
  blurmat[7]  * (src[i - w3 - 1] + src[i - w3 + 1] + src[i - w1 - 3] + src[i - w1 + 3] + src[i + w1 - 3] + src[i + w1 + 3] + src[i + w3 - 1] + src[i + w3 + 1]) + \
  blurmat[6]  * (src[i - w3] + src[i - 3] + src[i + 3] + src[i + w3]) + \
  blurmat[5]  * (src[i - w2 - 2] + src[i - w2 + 2] + src[i + w2 - 2] + src[i + w2 + 2]) + \
  blurmat[4]  * (src[i - w2 - 1] + src[i - w2 + 1] + src[i - w1 - 2] + src[i - w1 + 2] + src[i + w1 - 2] + src[i + w1 + 2] + src[i + w2 - 1] + src[i + w2 + 1]) + \
  blurmat[3]  * (src[i - w2] + src[i - 2] + src[i + 2] + src[i + w2]) + \
  blurmat[2]  * (src[i - w1 - 1] + src[i - w1 + 1] + src[i + w1 - 1] + src[i + w1 + 1]) + \
  blurmat[1]  * (src[i - w1] + src[i - 1] + src[i + 1] + src[i + w1]) + \
  blurmat[0]  * src[i] )
 
void dt_masks_blur_9x9(float *const restrict src, float *const restrict out, const int width, const int height, const float sigma)
{
  float blurmat[13];
  dt_masks_blur_9x9_coeff(blurmat, sigma);
 
  const int w1 = width;
  const int w2 = 2*width;
  const int w3 = 3*width;
  const int w4 = 4*width;
  __OMP_FOR_SIMD__(aligned(src, out : 64) if((size_t)width * height > 50000))
  for(int row = 4; row < height - 4; row++)
  {
    const int row_off = row * width;
    for(int col = 4; col < width - 4; col++)
    {
      const int i = row_off + col;
      out[i] = fminf(1.0f, fmaxf(0.0f, FAST_BLUR_9));
    }
  }
  dt_masks_extend_border(out, width, height, 4);
}
 
void _masks_blur_13x13_coeff(float *c, const float sigma)
{
  float kernel[13][13];
  const float temp = -2.0f * sqf(sigma);
  const float range = sqf(3.0f * 2.0f);
  float sum = 0.0f;
  for(int k = -6; k <= 6; k++)
  {
    for(int j = -6; j <= 6; j++)
    {
      if((sqf(k) + sqf(j)) <= range)
      {
        kernel[k + 6][j + 6] = expf((sqf(k) + sqf(j)) / temp);
        sum += kernel[k + 6][j + 6];
      }
      else
        kernel[k + 6][j + 6] = 0.0f;
    }
  }
  for(int i = 0; i < 13; i++)
    for(int j = 0; j < 13; j++)
      kernel[i][j] /= sum;
  
  // FIXME: are you for real ? managing arrays with loops and index shift much ?
  /* c60 */ c[0]  = kernel[0][6];
  /* c53 */ c[1]  = kernel[1][3];
  /* c52 */ c[2]  = kernel[1][4];
  /* c51 */ c[3]  = kernel[1][5];
  /* c50 */ c[4]  = kernel[1][6];
  /* c44 */ c[5]  = kernel[2][2];
  /* c42 */ c[6]  = kernel[2][4];
  /* c41 */ c[7]  = kernel[2][5];
  /* c40 */ c[8]  = kernel[2][6];
  /* c33 */ c[9]  = kernel[3][3];
  /* c32 */ c[10] = kernel[3][4];
  /* c31 */ c[11] = kernel[3][5];
  /* c30 */ c[12] = kernel[3][6];
  /* c22 */ c[13] = kernel[4][4];
  /* c21 */ c[14] = kernel[4][5];
  /* c20 */ c[15] = kernel[4][6];
  /* c11 */ c[16] = kernel[5][5];
  /* c10 */ c[17] = kernel[5][6];
  /* c00 */ c[18] = kernel[6][6];
}
 
 
__DT_CLONE_TARGETS__
void dt_masks_calc_rawdetail_mask(float *const restrict src, float *const restrict mask, float *const restrict tmp,
                                  const int width, const int height, const dt_aligned_pixel_t wb)
{
  const int msize = width * height;
  __OMP_FOR_SIMD__(aligned(tmp, src : 64) if(msize > 50000))
  for(int idx =0; idx < msize; idx++)
  {
    const float val = 0.333333333f * (fmaxf(src[4 * idx], 0.0f) / wb[0] + fmaxf(src[4 * idx + 1], 0.0f) / wb[1] + fmaxf(src[4 * idx + 2], 0.0f) / wb[2]);
    tmp[idx] = sqrtf(val); // add a gamma. sqrtf should make noise variance the same for all image
  }
 
  const float scale = 1.0f / 16.0f;
  __OMP_PARALLEL_FOR_SIMD__(aligned(mask, tmp : 64) if((size_t)width * height > 50000))
  for(int row = 1; row < height - 1; row++)
  {
    for(int col = 1, idx = row * width + col; col < width - 1; col++, idx++)
    {
      // scharr operator
      const float gx = 47.0f * (tmp[idx-width-1] - tmp[idx-width+1])
                    + 162.0f * (tmp[idx-1]       - tmp[idx+1])
                     + 47.0f * (tmp[idx+width-1] - tmp[idx+width+1]);
      const float gy = 47.0f * (tmp[idx-width-1] - tmp[idx+width-1])
                    + 162.0f * (tmp[idx-width]   - tmp[idx+width])
                     + 47.0f * (tmp[idx-width+1] - tmp[idx+width+1]);
      const float gradient_magnitude = dt_fast_hypotf(gx / 256.0f, gy / 256.0f);
      mask[idx] = scale * gradient_magnitude;
      // Original code from rt
      // tmp[idx] = scale * sqrtf(sqf(src[idx+1] - src[idx-1]) + sqf(src[idx + width]   - src[idx - width]) +
      //                          sqf(src[idx+2] - src[idx-2]) + sqf(src[idx + 2*width] - src[idx - 2*width]));
    }
  }
  dt_masks_extend_border(mask, width, height, 1);
}
 
static inline float calcBlendFactor(float val, float threshold)
{
    // sigmoid function
    // result is in ]0;1] range
    // inflexion point is at (x, y) (threshold, 0.5)
    return 1.0f / (1.0f + dt_fast_expf(16.0f - (16.0f / threshold) * val));
}
 
void dt_masks_calc_detail_mask(float *const restrict src, float *const restrict out, float *const restrict tmp, const int width, const int height, const float threshold, const gboolean detail)
{
  const int msize = width * height;
  __OMP_FOR_SIMD__(aligned(src, tmp, out : 64) if(msize > 50000))
  for(int idx = 0; idx < msize; idx++)
  {
    const float blend = calcBlendFactor(src[idx], threshold);
    tmp[idx] = detail ? blend : 1.0f - blend;
  }
  dt_masks_blur_9x9(tmp, out, width, height, 2.0f);
}
#undef FAST_BLUR_5
#undef FAST_BLUR_9
 
// clang-format off
// modelines: These editor modelines have been set for all relevant files by tools/update_modelines.py
// vim: shiftwidth=2 expandtab tabstop=2 cindent
// kate: tab-indents: off; indent-width 2; replace-tabs on; indent-mode cstyle; remove-trailing-spaces modified;
// clang-format on