color_picker.c

Go to the documentation of this file.

/*
    This file is part of darktable,
    Copyright (C) 2016 Roman Lebedev.
    Copyright (C) 2019 Andreas Schneider.
    Copyright (C) 2019 Edgardo Hoszowski.
    Copyright (C) 2019-2021 Pascal Obry.
    Copyright (C) 2020-2021 Harold le Clément de Saint-Marcq.
    Copyright (C) 2020 Heiko Bauke.
    Copyright (C) 2020 Hubert Kowalski.
    Copyright (C) 2020-2021 Ralf Brown.
    Copyright (C) 2021 Hanno Schwalm.
    Copyright (C) 2022, 2025-2026 Aurélien PIERRE.
    Copyright (C) 2022 Martin Bařinka.
    Copyright (C) 2022 Philipp Lutz.
    
    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/>.
*/
 
#include "common/color_picker.h"
#include "common/colorspaces_inline_conversions.h"
#include "pixel/bspline.h"
#include "system/macros.h"
#include "system/openmp.h"
#include "system/mem_alloc.h"
#include "system/simd.h"
#include "common/logging.h"
#include "common/times.h"
#include "caches/pixelpipe_cache_alloc.h"
#include "colorprofiles/iop_profile.h"
#include "pixel/format.h"
#include "develop/imageop_math.h"
 
static inline size_t _box_size(const int *const box)
{
  return (size_t)((box[3] - box[1]) * (box[2] - box[0]));
}
 
static inline void rgb_to_JzCzhz(const dt_aligned_pixel_t rgb, dt_aligned_pixel_t JzCzhz,
                                 const dt_iop_order_iccprofile_info_t *const profile);
 
static inline gboolean _picker_cst_wants_lab(const dt_iop_colorspace_type_t picker_cst)
{
  return picker_cst == IOP_CS_LAB || picker_cst == IOP_CS_LCH;
}
 
static void _color_picker_convert_buffer(const float *const restrict input, float *const restrict output,
                                         const size_t pixels, const dt_iop_colorspace_type_t image_cst,
                                         const dt_iop_colorspace_type_t picker_cst,
                                         const dt_iop_order_iccprofile_info_t *const profile)
{
  /* Two steps, so every pairing of the blending tabs is served: first bring the pixel into the
   * family the picker space derives from (Lab for Lab/LCh, RGB for RGB/HSL/JzCzhz) -- the only
   * step that needs the profile -- then derive the picker space from it. A Lab module blended in
   * RGB (scene) asks Lab -> JzCzhz, an RGB module blended in Lab asks RGB -> LCh. */
  const gboolean image_lab = (image_cst == IOP_CS_LAB);
  const gboolean want_lab = _picker_cst_wants_lab(picker_cst);
 
  __OMP_PARALLEL_FOR__()
  for(size_t k = 0; k < pixels; k++)
  {
    const size_t offset = 4 * k;
    dt_aligned_pixel_t base = { input[offset], input[offset + 1], input[offset + 2], 0.0f };
 
    if(image_lab && !want_lab)
      dt_ioppr_lab_to_rgb_matrix(input + offset, base, profile->matrix_out_transposed, profile->lut_out,
                                 profile->unbounded_coeffs_out, profile->lutsize, profile->nonlinearlut);
    else if(!image_lab && want_lab)
      dt_ioppr_rgb_matrix_to_lab(input + offset, base, profile->matrix_in_transposed, profile->lut_in,
                                 profile->unbounded_coeffs_in, profile->lutsize, profile->nonlinearlut);
 
    switch(picker_cst)
    {
      case IOP_CS_LCH:
        dt_Lab_2_LCH(base, output + offset);
        break;
      case IOP_CS_HSL:
        dt_RGB_2_HSL(base, output + offset);
        break;
      case IOP_CS_JZCZHZ:
        rgb_to_JzCzhz(base, output + offset, profile);
        break;
      default:
        for(int c = 0; c < 3; c++) output[offset + c] = base[c];
        break;
    }
    output[offset + 3] = input[offset + 3];
  }
}
 
__OMP_DECLARE_SIMD__(aligned(rgb, JzCzhz: 16) uniform(profile))
static inline void rgb_to_JzCzhz(const dt_aligned_pixel_t rgb, dt_aligned_pixel_t JzCzhz,
                                 const dt_iop_order_iccprofile_info_t *const profile)
{
  dt_aligned_pixel_t XYZ_D65 = { 0.0f, 0.0f, 0.0f };
  dt_aligned_pixel_t JzAzBz = { 0.0f, 0.0f, 0.0f };
 
  if(profile)
  {
    dt_aligned_pixel_t XYZ_D50 = { 0.0f, 0.0f, 0.0f };
    dt_ioppr_rgb_matrix_to_xyz(rgb, XYZ_D50, profile->matrix_in_transposed, profile->lut_in, profile->unbounded_coeffs_in,
                               profile->lutsize, profile->nonlinearlut);
    dt_XYZ_D50_2_XYZ_D65(XYZ_D50, XYZ_D65);
  }
  else
  {
    // This should not happen (we don't know what RGB is), but use this when profile is not defined
    dt_XYZ_D50_2_XYZ_D65(rgb, XYZ_D65);
  }
 
  dt_XYZ_2_JzAzBz(XYZ_D65, JzAzBz);
  dt_JzAzBz_2_JzCzhz(JzAzBz, JzCzhz);
}
static inline void _color_picker_rgb_or_lab(dt_aligned_pixel_t avg, dt_aligned_pixel_t min, dt_aligned_pixel_t max,
                                            const float *const pixels, const float w, const size_t width)
{
  for(size_t i = 0; i < width; i += 4)
  {
    dt_aligned_pixel_t pick = { pixels[i], pixels[i + 1], pixels[i + 2], 0.0f };
    for(size_t k = 0; k < 4; k++)
    {
      avg[k] += w * pick[k];
      min[k] = fminf(min[k], pick[k]);
      max[k] = fmaxf(max[k], pick[k]);
    }
  }
}
static inline void _color_picker_direct_hsl(dt_aligned_pixel_t avg, dt_aligned_pixel_t min, dt_aligned_pixel_t max,
                                            const float *const pixels, const float w, const size_t width)
{
  for(size_t i = 0; i < width; i += 4)
  {
    dt_aligned_pixel_t pick = { pixels[i], pixels[i + 1], pixels[i + 2], 0.0f };
    pick[3] = pick[0] < 0.5f ? pick[0] + 0.5f : pick[0] - 0.5f;
    for(size_t k = 0; k < 4; k++)
    {
      avg[k] += w * pick[k];
      min[k] = fminf(min[k], pick[k]);
      max[k] = fmaxf(max[k], pick[k]);
    }
  }
}
static inline void _color_picker_direct_lch_or_jzczhz(dt_aligned_pixel_t avg, dt_aligned_pixel_t min,
                                                       dt_aligned_pixel_t max, const float *const pixels,
                                                       const float w, const size_t width)
{
  for(size_t i = 0; i < width; i += 4)
  {
    dt_aligned_pixel_t pick = { pixels[i], pixels[i + 1], pixels[i + 2], 0.0f };
    pick[3] = pick[2] < 0.5f ? pick[2] + 0.5f : pick[2] - 0.5f;
    for(size_t k = 0; k < 4; k++)
    {
      avg[k] += w * pick[k];
      min[k] = fminf(min[k], pick[k]);
      max[k] = fmaxf(max[k], pick[k]);
    }
  }
}
static inline void _color_picker_lch(dt_aligned_pixel_t avg, dt_aligned_pixel_t min, dt_aligned_pixel_t max,
                                     const float *const pixels, const float w, const size_t width)
{
  for(size_t i = 0; i < width; i += 4)
  {
    dt_aligned_pixel_t pick;
    dt_Lab_2_LCH(pixels + i, pick);
    pick[3] = pick[2] < 0.5f ? pick[2] + 0.5f : pick[2] - 0.5f;
    for(size_t k = 0; k < 4; k++)
    {
      avg[k] += w * pick[k];
      min[k] = fminf(min[k], pick[k]);
      max[k] = fmaxf(max[k], pick[k]);
    }
  }
}
static inline void _color_picker_hsl(dt_aligned_pixel_t avg, dt_aligned_pixel_t min, dt_aligned_pixel_t max,
                                     const float *const pixels, const float w, const size_t width)
{
  for(size_t i = 0; i < width; i += 4)
  {
    dt_aligned_pixel_t pick;
    dt_RGB_2_HSL(pixels + i, pick);
    pick[3] = pick[0] < 0.5f ? pick[0] + 0.5f : pick[0] - 0.5f;
    for(size_t k = 0; k < 4; k++)
    {
      avg[k] += w * pick[k];
      min[k] = fminf(min[k], pick[k]);
      max[k] = fmaxf(max[k], pick[k]);
    }
  }
}
static inline void _color_picker_jzczhz(dt_aligned_pixel_t avg, dt_aligned_pixel_t min, dt_aligned_pixel_t max,
                                        const float *const pixels, const float w, const size_t width,
                                        const dt_iop_order_iccprofile_info_t *const profile)
{
  for(size_t i = 0; i < width; i += 4)
  {
    dt_aligned_pixel_t pick;
    rgb_to_JzCzhz(pixels + i, pick, profile);
    pick[3] = pick[2] < 0.5f ? pick[2] + 0.5f : pick[2] - 0.5f;
    for(size_t k = 0; k < 4; k++)
    {
      avg[k] += w * pick[k];
      min[k] = fminf(min[k], pick[k]);
      max[k] = fmaxf(max[k], pick[k]);
    }
  }
}
 
static void color_picker_helper_4ch_seq(const dt_iop_buffer_dsc_t *const dsc, const float *const pixel,
                                        const dt_iop_roi_t *const roi, const int *const box,
                                        dt_aligned_pixel_t picked_color, dt_aligned_pixel_t picked_color_min,
                                        dt_aligned_pixel_t picked_color_max, const dt_iop_colorspace_type_t cst_to,
                                        const dt_iop_order_iccprofile_info_t *const profile)
{
  const int width = roi->width;
 
  const size_t size = _box_size(box);
  const size_t stride = 4 * (size_t)(box[2] - box[0]);
  const size_t off_mul = 4 * width;
  const size_t off_add = 4 * box[0];
 
  const float w = 1.0f / (float)size;
 
  // code path for small region, especially for color picker point mode
  if(cst_to == IOP_CS_LCH)
  {
    for(size_t j = box[1]; j < box[3]; j++)
    {
      const size_t offset = j * off_mul + off_add;
      _color_picker_lch(picked_color, picked_color_min, picked_color_max, pixel + offset, w, stride);
    }
  }
  else if(cst_to == IOP_CS_HSL)
  {
    for(size_t j = box[1]; j < box[3]; j++)
    {
      const size_t offset = j * off_mul + off_add;
      _color_picker_hsl(picked_color, picked_color_min, picked_color_max, pixel + offset, w, stride);
    }
  }
  else if(cst_to == IOP_CS_JZCZHZ)
  {
    for(size_t j = box[1]; j < box[3]; j++)
    {
      const size_t offset = j * off_mul + off_add;
      _color_picker_jzczhz(picked_color, picked_color_min, picked_color_max, pixel + offset, w, stride, profile);
    }
  }
  else
  {
    for(size_t j = box[1]; j < box[3]; j++)
    {
      const size_t offset = j * off_mul + off_add;
      _color_picker_rgb_or_lab(picked_color, picked_color_min, picked_color_max, pixel + offset, w, stride);
    }
  }
}
 
static void color_picker_helper_4ch_parallel(const dt_iop_buffer_dsc_t *const dsc, const float *const pixel,
                                             const dt_iop_roi_t *const roi, const int *const box,
                                             dt_aligned_pixel_t picked_color, dt_aligned_pixel_t picked_color_min,
                                             dt_aligned_pixel_t picked_color_max, const dt_iop_colorspace_type_t cst_to,
                                             const dt_iop_order_iccprofile_info_t *const profile)
{
  const int width = roi->width;
 
  const size_t size = _box_size(box);
  const size_t stride = 4 * (size_t)(box[2] - box[0]);
  const size_t off_mul = 4 * width;
  const size_t off_add = 4 * box[0];
 
  const float w = 1.0f / (float)size;
 
  const size_t numthreads = MAX(1, (size_t)omp_get_max_threads());
 
  size_t allocsize;
  float *const restrict mean = dt_pixelpipe_cache_alloc_perthread_float(4, &allocsize);
  float *const restrict mmin = dt_pixelpipe_cache_alloc_perthread_float(4, &allocsize);
  float *const restrict mmax = dt_pixelpipe_cache_alloc_perthread_float(4, &allocsize);
 
  if(IS_NULL_PTR(mean) || IS_NULL_PTR(mmax) || IS_NULL_PTR(mmin))
    goto error;
 
  for(int n = 0; n < allocsize * numthreads; n++)
  {
    mean[n] = 0.0f;
    mmin[n] = INFINITY;
    mmax[n] = -INFINITY;
  }
 
  if(cst_to == IOP_CS_LCH)
  {
    __OMP_PARALLEL__()
    {
      float *const restrict tmean = dt_get_perthread(mean,allocsize);
      float *const restrict tmmin = dt_get_perthread(mmin,allocsize);
      float *const restrict tmmax = dt_get_perthread(mmax,allocsize);
      __OMP_FOR__()
      for(size_t j = box[1]; j < box[3]; j++)
      {
        const size_t offset = j * off_mul + off_add;
        _color_picker_lch(tmean, tmmin, tmmax, pixel + offset, w, stride);
      }
    }
  }
  else if(cst_to == IOP_CS_HSL)
  {
    __OMP_PARALLEL__()
    {
      float *const restrict tmean = dt_get_perthread(mean,allocsize);
      float *const restrict tmmin = dt_get_perthread(mmin,allocsize);
      float *const restrict tmmax = dt_get_perthread(mmax,allocsize);
      __OMP_FOR__()
      for(size_t j = box[1]; j < box[3]; j++)
      {
        const size_t offset = j * off_mul + off_add;
        _color_picker_hsl(tmean, tmmin, tmmax, pixel + offset, w, stride);
      }
    }
  }
  else if(cst_to == IOP_CS_JZCZHZ)
  {
    __OMP_PARALLEL__()
    {
      float *const restrict tmean = dt_get_perthread(mean,allocsize);
      float *const restrict tmmin = dt_get_perthread(mmin,allocsize);
      float *const restrict tmmax = dt_get_perthread(mmax,allocsize);
      __OMP_FOR__()
      for(size_t j = box[1]; j < box[3]; j++)
      {
        const size_t offset = j * off_mul + off_add;
        _color_picker_jzczhz(tmean, tmmin, tmmax, pixel + offset, w, stride, profile);
      }
    }
  }
  else
  {
    __OMP_PARALLEL__()
    {
      float *const restrict tmean = dt_get_perthread(mean,allocsize);
      float *const restrict tmmin = dt_get_perthread(mmin,allocsize);
      float *const restrict tmmax = dt_get_perthread(mmax,allocsize);
      __OMP_FOR__()
      for(size_t j = box[1]; j < box[3]; j++)
      {
        const size_t offset = j * off_mul + off_add;
        _color_picker_rgb_or_lab(tmean, tmmin, tmmax, pixel + offset, w, stride);
      }
    }
  }
 
  for(int n = 0; n < numthreads; n++)
  {
    for(int k = 0; k < 4; k++)
    {
      picked_color[k] += mean[allocsize * n + k];
      picked_color_min[k] = fminf(picked_color_min[k], mmin[allocsize * n + k]);
      picked_color_max[k] = fmaxf(picked_color_max[k], mmax[allocsize * n + k]);
    }
  }
 
error:;
  dt_pixelpipe_cache_free_align(mmax);
  dt_pixelpipe_cache_free_align(mmin);
  dt_pixelpipe_cache_free_align(mean);
}
 
static void color_picker_helper_4ch(const dt_iop_buffer_dsc_t *dsc, const float *const pixel,
                                    const dt_iop_roi_t *roi, const int *const box, dt_aligned_pixel_t picked_color,
                                    dt_aligned_pixel_t picked_color_min, dt_aligned_pixel_t picked_color_max,
                                    const dt_iop_colorspace_type_t cst_to,
                                    const dt_iop_order_iccprofile_info_t *const profile)
{
  const size_t size = _box_size(box);
 
  if(size > 100) // avoid inefficient multi-threading in case of small region size (arbitrary limit)
    return color_picker_helper_4ch_parallel(dsc, pixel, roi, box, picked_color, picked_color_min,
                                            picked_color_max, cst_to, profile);
  else
    return color_picker_helper_4ch_seq(dsc, pixel, roi, box, picked_color, picked_color_min, picked_color_max,
                                       cst_to, profile);
}
 
static void color_picker_helper_4ch_converted_seq(const float *const pixel, const dt_iop_roi_t *const roi,
                                                  const int *const box, dt_aligned_pixel_t picked_color,
                                                  dt_aligned_pixel_t picked_color_min,
                                                  dt_aligned_pixel_t picked_color_max,
                                                  const dt_iop_colorspace_type_t picker_cst)
{
  const int width = roi->width;
  const size_t size = _box_size(box);
  const size_t stride = 4 * (size_t)(box[2] - box[0]);
  const size_t off_mul = 4 * width;
  const size_t off_add = 4 * box[0];
  const float w = 1.0f / (float)size;
 
  for(size_t j = box[1]; j < box[3]; j++)
  {
    const size_t offset = j * off_mul + off_add;
    if(picker_cst == IOP_CS_HSL)
      _color_picker_direct_hsl(picked_color, picked_color_min, picked_color_max, pixel + offset, w, stride);
    else if(picker_cst == IOP_CS_LCH || picker_cst == IOP_CS_JZCZHZ)
      _color_picker_direct_lch_or_jzczhz(picked_color, picked_color_min, picked_color_max, pixel + offset, w, stride);
    else
      _color_picker_rgb_or_lab(picked_color, picked_color_min, picked_color_max, pixel + offset, w, stride);
  }
}
 
static void color_picker_helper_4ch_converted_parallel(const float *const pixel, const dt_iop_roi_t *const roi,
                                                       const int *const box, dt_aligned_pixel_t picked_color,
                                                       dt_aligned_pixel_t picked_color_min,
                                                       dt_aligned_pixel_t picked_color_max,
                                                       const dt_iop_colorspace_type_t picker_cst)
{
  const int width = roi->width;
  const size_t size = _box_size(box);
  const size_t stride = 4 * (size_t)(box[2] - box[0]);
  const size_t off_mul = 4 * width;
  const size_t off_add = 4 * box[0];
  const float w = 1.0f / (float)size;
  const size_t numthreads = dt_get_num_openmp_threads();
 
  size_t allocsize;
  float *const restrict mean = dt_pixelpipe_cache_alloc_perthread_float(4, &allocsize);
  float *const restrict mmin = dt_pixelpipe_cache_alloc_perthread_float(4, &allocsize);
  float *const restrict mmax = dt_pixelpipe_cache_alloc_perthread_float(4, &allocsize);
 
  if(IS_NULL_PTR(mean) || IS_NULL_PTR(mmax) || IS_NULL_PTR(mmin))
    goto error;
 
  for(int n = 0; n < allocsize * numthreads; n++)
  {
    mean[n] = 0.0f;
    mmin[n] = INFINITY;
    mmax[n] = -INFINITY;
  }
  __OMP_PARALLEL__()
  {
    float *const restrict tmean = dt_get_perthread(mean,allocsize);
    float *const restrict tmmin = dt_get_perthread(mmin,allocsize);
    float *const restrict tmmax = dt_get_perthread(mmax,allocsize);
    __OMP_FOR__()
    for(size_t j = box[1]; j < box[3]; j++)
    {
      const size_t offset = j * off_mul + off_add;
      if(picker_cst == IOP_CS_HSL)
        _color_picker_direct_hsl(tmean, tmmin, tmmax, pixel + offset, w, stride);
      else if(picker_cst == IOP_CS_LCH || picker_cst == IOP_CS_JZCZHZ)
        _color_picker_direct_lch_or_jzczhz(tmean, tmmin, tmmax, pixel + offset, w, stride);
      else
        _color_picker_rgb_or_lab(tmean, tmmin, tmmax, pixel + offset, w, stride);
    }
  }
 
  for(int n = 0; n < numthreads; n++)
  {
    float *const restrict tmean = dt_get_bythread(mean,allocsize,n);
    float *const restrict tmmin = dt_get_bythread(mmin,allocsize,n);
    float *const restrict tmmax = dt_get_bythread(mmax,allocsize,n);
 
    for_four_channels(k)
    {
      picked_color[k] += tmean[k];
      picked_color_min[k] = fminf(picked_color_min[k], tmmin[k]);
      picked_color_max[k] = fmaxf(picked_color_max[k], tmmax[k]);
    }
  }
 
error:
  dt_pixelpipe_cache_free_align(mmax);
  dt_pixelpipe_cache_free_align(mmin);
  dt_pixelpipe_cache_free_align(mean);
}
 
static void color_picker_helper_4ch_converted(const float *const pixel, const dt_iop_roi_t *const roi,
                                              const int *const box, dt_aligned_pixel_t picked_color,
                                              dt_aligned_pixel_t picked_color_min,
                                              dt_aligned_pixel_t picked_color_max,
                                              const dt_iop_colorspace_type_t picker_cst)
{
  if(_box_size(box) > 10000)
    color_picker_helper_4ch_converted_parallel(pixel, roi, box, picked_color, picked_color_min,
                                               picked_color_max, picker_cst);
  else
    color_picker_helper_4ch_converted_seq(pixel, roi, box, picked_color, picked_color_min,
                                          picked_color_max, picker_cst);
}
 
static void color_picker_helper_bayer_seq(const dt_iop_buffer_dsc_t *const dsc, const float *const pixel,
                                          const dt_iop_roi_t *const roi, const int *const box,
                                          dt_aligned_pixel_t picked_color, dt_aligned_pixel_t picked_color_min,
                                          dt_aligned_pixel_t picked_color_max)
{
  const int width = roi->width;
  const uint32_t filters = dsc->filters;
 
  uint32_t weights[4] = { 0u, 0u, 0u, 0u };
 
  // code path for small region, especially for color picker point mode
  for(size_t j = box[1]; j < box[3]; j++)
  {
    for(size_t i = box[0]; i < box[2]; i++)
    {
      const int c = FC(j + roi->y, i + roi->x, filters);
      const size_t k = width * j + i;
 
      const float v = pixel[k];
 
      picked_color[c] += v;
      picked_color_min[c] = fminf(picked_color_min[c], v);
      picked_color_max[c] = fmaxf(picked_color_max[c], v);
      weights[c]++;
    }
  }
 
  // and finally normalize data. For bayer, there is twice as much green.
  for(int c = 0; c < 4; c++)
  {
    picked_color[c] = weights[c] ? (picked_color[c] / (float)weights[c]) : 0.0f;
  }
}
 
static void color_picker_helper_bayer_parallel(const dt_iop_buffer_dsc_t *const dsc, const float *const pixel,
                                               const dt_iop_roi_t *const roi, const int *const box,
                                               dt_aligned_pixel_t picked_color, dt_aligned_pixel_t picked_color_min,
                                               dt_aligned_pixel_t picked_color_max)
{
  const int width = roi->width;
  const uint32_t filters = dsc->filters;
 
  uint32_t weights[4] = { 0u, 0u, 0u, 0u };
 
  const size_t numthreads = dt_get_num_openmp_threads();
 
  //TODO: convert to use dt_pixelpipe_cache_alloc_perthread
  float *const msum = malloc(sizeof(float) * numthreads * 4);
  float *const mmin = malloc(sizeof(float) * numthreads * 4);
  float *const mmax = malloc(sizeof(float) * numthreads * 4);
  uint32_t *const cnt = malloc(sizeof(uint32_t) * numthreads * 4);
 
  if(IS_NULL_PTR(msum) || IS_NULL_PTR(mmin) || IS_NULL_PTR(mmax) || IS_NULL_PTR(cnt))
    goto error;
 
  for(int n = 0; n < 4 * numthreads; n++)
  {
    msum[n] = 0.0f;
    mmin[n] = INFINITY;
    mmax[n] = -INFINITY;
    cnt[n] = 0u;
  }
  __OMP_PARALLEL__(num_threads(numthreads))
  {
    const int tnum = dt_get_thread_num();
 
    float *const tsum = msum + 4 * tnum;
    float *const tmmin = mmin + 4 * tnum;
    float *const tmmax = mmax + 4 * tnum;
    uint32_t *const tcnt = cnt + 4 * tnum;
    __OMP_FOR__(collapse(2))
    for(size_t j = box[1]; j < box[3]; j++)
    {
      for(size_t i = box[0]; i < box[2]; i++)
      {
        const int c = FC(j + roi->y, i + roi->x, filters);
        const size_t k = width * j + i;
 
        const float v = pixel[k];
 
        tsum[c] += v;
        tmmin[c] = fminf(tmmin[c], v);
        tmmax[c] = fmaxf(tmmax[c], v);
        tcnt[c]++;
      }
    }
  }
 
  for(int n = 0; n < numthreads; n++)
  {
    for(int c = 0; c < 4; c++)
    {
      picked_color[c] += msum[4 * n + c];
      picked_color_min[c] = fminf(picked_color_min[c], mmin[4 * n + c]);
      picked_color_max[c] = fmaxf(picked_color_max[c], mmax[4 * n + c]);
      weights[c] += cnt[4 * n + c];
    }
  }
 
  // and finally normalize data. For bayer, there is twice as much green.
  for(int c = 0; c < 4; c++)
  {
    picked_color[c] = weights[c] ? (picked_color[c] / (float)weights[c]) : 0.0f;
  }
 
error:;
  dt_free(cnt);
  dt_free(mmax);
  dt_free(mmin);
  dt_free(msum);
}
 
static void color_picker_helper_bayer(const dt_iop_buffer_dsc_t *dsc, const float *const pixel,
                                      const dt_iop_roi_t *roi, const int *const box, dt_aligned_pixel_t picked_color,
                                      dt_aligned_pixel_t picked_color_min, dt_aligned_pixel_t picked_color_max)
{
  const size_t size = _box_size(box);
 
  if(size > 100) // avoid inefficient multi-threading in case of small region size (arbitrary limit)
    return color_picker_helper_bayer_parallel(dsc, pixel, roi, box, picked_color, picked_color_min,
                                              picked_color_max);
  else
    return color_picker_helper_bayer_seq(dsc, pixel, roi, box, picked_color, picked_color_min, picked_color_max);
}
 
static void color_picker_helper_xtrans_seq(const dt_iop_buffer_dsc_t *const dsc, const float *const pixel,
                                           const dt_iop_roi_t *const roi, const int *const box,
                                           dt_aligned_pixel_t picked_color, dt_aligned_pixel_t picked_color_min,
                                           dt_aligned_pixel_t picked_color_max)
{
  const int width = roi->width;
  const uint8_t(*const xtrans)[6] = (const uint8_t(*const)[6])dsc->xtrans;
 
  uint32_t weights[3] = { 0u, 0u, 0u };
 
  // code path for small region, especially for color picker point mode
  for(size_t j = box[1]; j < box[3]; j++)
  {
    for(size_t i = box[0]; i < box[2]; i++)
    {
      const int c = FCxtrans(j, i, roi, xtrans);
      const size_t k = width * j + i;
 
      const float v = pixel[k];
 
      picked_color[c] += v;
      picked_color_min[c] = fminf(picked_color_min[c], v);
      picked_color_max[c] = fmaxf(picked_color_max[c], v);
      weights[c]++;
    }
  }
 
  // and finally normalize data.
  // X-Trans RGB weighting averages to 2:5:2 for each 3x3 cell
  for(int c = 0; c < 3; c++)
  {
    picked_color[c] /= (float)weights[c];
  }
}
 
static void color_picker_helper_xtrans_parallel(const dt_iop_buffer_dsc_t *const dsc, const float *const pixel,
                                                const dt_iop_roi_t *const roi, const int *const box,
                                                dt_aligned_pixel_t picked_color, dt_aligned_pixel_t picked_color_min,
                                                dt_aligned_pixel_t picked_color_max)
{
  const int width = roi->width;
  const uint8_t(*const xtrans)[6] = (const uint8_t(*const)[6])dsc->xtrans;
 
  uint32_t weights[3] = { 0u, 0u, 0u };
 
  const size_t numthreads = MAX(1, (size_t)omp_get_max_threads());
 
  //TODO: convert to use dt_pixelpipe_cache_alloc_perthread
  float *const mmin = malloc(sizeof(float) * numthreads * 3);
  float *const msum = malloc(sizeof(float) * numthreads * 3);
  float *const mmax = malloc(sizeof(float) * numthreads * 3);
  uint32_t *const cnt = malloc(sizeof(uint32_t) * numthreads * 3);
 
  if(IS_NULL_PTR(mmin) || IS_NULL_PTR(msum) || IS_NULL_PTR(mmax) || IS_NULL_PTR(cnt))
    goto error;
 
 
  for(int n = 0; n < 3 * numthreads; n++)
  {
    msum[n] = 0.0f;
    mmin[n] = INFINITY;
    mmax[n] = -INFINITY;
    cnt[n] = 0u;
  }
  __OMP_PARALLEL__(num_threads(numthreads))
  {
    const int tnum = dt_get_thread_num();
 
    float *const tsum = msum + 3 * tnum;
    float *const tmmin = mmin + 3 * tnum;
    float *const tmmax = mmax + 3 * tnum;
    uint32_t *const tcnt = cnt + 3 * tnum;
    __OMP_FOR__(collapse(2))
    for(size_t j = box[1]; j < box[3]; j++)
    {
      for(size_t i = box[0]; i < box[2]; i++)
      {
        const int c = FCxtrans(j, i, roi, xtrans);
        const size_t k = width * j + i;
 
        const float v = pixel[k];
 
        tsum[c] += v;
        tmmin[c] = fminf(tmmin[c], v);
        tmmax[c] = fmaxf(tmmax[c], v);
        tcnt[c]++;
      }
    }
  }
 
  for(int n = 0; n < numthreads; n++)
  {
    for(int c = 0; c < 3; c++)
    {
      picked_color[c] += msum[3 * n + c];
      picked_color_min[c] = fminf(picked_color_min[c], mmin[3 * n + c]);
      picked_color_max[c] = fmaxf(picked_color_max[c], mmax[3 * n + c]);
      weights[c] += cnt[3 * n + c];
    }
  }
 
  // and finally normalize data.
  // X-Trans RGB weighting averages to 2:5:2 for each 3x3 cell
  for(int c = 0; c < 3; c++)
  {
    picked_color[c] /= (float)weights[c];
  }
 
error:;
  dt_free(cnt);
  dt_free(mmax);
  dt_free(mmin);
  dt_free(msum);
}
 
static void color_picker_helper_xtrans(const dt_iop_buffer_dsc_t *dsc, const float *const pixel,
                                       const dt_iop_roi_t *roi, const int *const box, dt_aligned_pixel_t picked_color,
                                       dt_aligned_pixel_t picked_color_min, dt_aligned_pixel_t picked_color_max)
{
  const size_t size = _box_size(box);
 
  if(size > 100) // avoid inefficient multi-threading in case of small region size (arbitrary limit)
    return color_picker_helper_xtrans_parallel(dsc, pixel, roi, box, picked_color, picked_color_min,
                                               picked_color_max);
  else
    return color_picker_helper_xtrans_seq(dsc, pixel, roi, box, picked_color, picked_color_min, picked_color_max);
}
 
// picked_color, picked_color_min and picked_color_max should be aligned
void dt_color_picker_helper(const dt_iop_buffer_dsc_t *dsc, const float *const pixel, const dt_iop_roi_t *roi,
                            const int *const box, dt_aligned_pixel_t picked_color, dt_aligned_pixel_t picked_color_min,
                            dt_aligned_pixel_t picked_color_max, const dt_iop_colorspace_type_t image_cst,
                            const dt_iop_colorspace_type_t picker_cst,
                            const dt_iop_order_iccprofile_info_t *const profile)
{
  dt_times_t start_time = { 0 }, end_time = { 0 };
  if(dt_get_debug_flags() & DT_DEBUG_PERF) dt_get_times(&start_time);
 
  if(dsc->channels == 4u)
  {
    // Denoise the image
    size_t padded_size;
    float *const restrict denoised = dt_pixelpipe_cache_alloc_align_float_cache(4 * roi->width * roi->height, 0);
    float *converted = NULL;
    float *const tempbuf = dt_pixelpipe_cache_alloc_perthread_float(4 * roi->width, &padded_size); // TODO: alloc in caller
    if(IS_NULL_PTR(tempbuf) || IS_NULL_PTR(denoised))
      goto error;
 
    // blur without clipping negatives because Lab a and b channels can be legitimately negative
    blur_2D_Bspline(pixel, denoised, tempbuf, roi->width, roi->height, 1, FALSE);
 
    // Lab and RGB buffers convert to every picker space derived from either family; only crossing
    // from one family to the other needs the profile.
    const gboolean image_convertible = (image_cst == IOP_CS_LAB) || dt_iop_colorspace_is_rgb(image_cst);
    const gboolean picker_convertible = _picker_cst_wants_lab(picker_cst) || dt_iop_colorspace_is_rgb(picker_cst)
                                        || (picker_cst == IOP_CS_HSL) || (picker_cst == IOP_CS_JZCZHZ);
    const gboolean crosses_family = (image_cst == IOP_CS_LAB) != _picker_cst_wants_lab(picker_cst);
    const gboolean convertible
        = image_convertible && picker_convertible && (!crosses_family || !IS_NULL_PTR(profile));
 
    if(((image_cst == picker_cst) || (picker_cst == IOP_CS_NONE)))
      color_picker_helper_4ch(dsc, denoised, roi, box, picked_color, picked_color_min, picked_color_max, picker_cst, profile);
    else if(convertible)
    {
      /* The picker samples module input/output buffers after the previous piece has written them.
         When the requested picker colorspace differs from that buffer colorspace, we need a real
         conversion before averaging. Falling back to raw channel statistics would silently report
         Lab values as RGB (or the other way around), which makes module-side picker feedback wrong. */
      converted = dt_pixelpipe_cache_alloc_align_float_cache(4 * roi->width * roi->height, 0);
      if(IS_NULL_PTR(converted))
        goto error;
 
      _color_picker_convert_buffer(denoised, converted, (size_t)roi->width * roi->height, image_cst, picker_cst, profile);
      color_picker_helper_4ch_converted(converted, roi, box, picked_color, picked_color_min, picked_color_max, picker_cst);
    }
    else // This is a fallback, better than crashing as happens with monochromes
      color_picker_helper_4ch(dsc, denoised, roi, box, picked_color, picked_color_min, picked_color_max, picker_cst, profile);
 
  error:;
    dt_pixelpipe_cache_free_align(converted);
    dt_pixelpipe_cache_free_align(denoised);
    dt_pixelpipe_cache_free_align(tempbuf);
  }
  else if(dsc->channels == 1u && dsc->filters != 0u && dsc->filters != 9u)
    color_picker_helper_bayer(dsc, pixel, roi, box, picked_color, picked_color_min, picked_color_max);
  else if(dsc->channels == 1u && dsc->filters == 9u)
    color_picker_helper_xtrans(dsc, pixel, roi, box, picked_color, picked_color_min, picked_color_max);
  else
    dt_unreachable_codepath();
 
  if(dt_get_debug_flags() & DT_DEBUG_PERF)
  {
    dt_get_times(&end_time);
    fprintf(stderr, "colorpicker stats reading took %.3f secs (%.3f CPU)\n",
        end_time.clock - start_time.clock, end_time.user - start_time.user);
  }
}
 
// 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