Search code examples
cparallel-processingopenmp

Parallelization of a reverse-lookup with OpenMP

I'm trying to find the closest correlated color temperature to any arbitrary chromaticity. That is, for any (x, y) point in the following graph, I want the closest point that belong to the Planckian locus, and from that point, I want the correlated black body temperature : enter image description here

The parametric function of the black body curve is a polynomial :


#ifdef _OPENMP
#pragma omp declare simd
#endif
static inline void CCT_to_xy_blackbody(const float t, float *x, float *y)
{
  // Take correlated color temperature in K and find the closest blackbody illuminant in 1667 K - 250000 K
  float x_temp = 0.f;
  float y_temp = 0.f;

  if(t >= 1667.f && t <= 4000.f)
    x_temp = -0.2661239e9f / (t * t * t) - 0.2343589e6f / (t * t)  + 0.8776956e3f / t + 0.179910f;
  else if(t > 4000.f && t <= 25000.f)
    x_temp = -3.0258469e9f / (t * t * t) + 2.1070379e6f / (t * t)  + 0.2226347e3f / t + 0.240390f;

  if(t >= 1667.f && t <= 2222.f)
    y_temp = -1.1063814f * x_temp * x_temp * x_temp - 1.34811020f * x_temp * x_temp + 2.18555832f * x_temp - 0.20219683f;
  else if(t > 2222.f && t <= 4000.f)
    y_temp = -0.9549476f * x_temp * x_temp * x_temp - 1.37418593f * x_temp * x_temp + 2.09137015f * x_temp - 0.16748867f;
  else if(t > 4000.f && t <= 25000.f)
    y_temp =  3.0817580f * x_temp * x_temp * x_temp - 5.87338670f * x_temp * x_temp + 3.75112997f * x_temp - 0.37001483f;

  *x = x_temp;
  *y = y_temp;
}

So the idea is to build a 2D LUT : T -> x_bb, y_bb, measure the distance between (x, y) and every set of (x_bb, y_bb), find the minimum of that distance, and the corresponding index in the LUT will yield the temperature.

Here is the function where I build and search the LUT at the same time:

static inline float CCT_reverse_lookup(const float x, const float y)
{
  // Find out the closest correlated color temperature (closest point over the planckian locus)
  // for any arbitrary x, y chromaticity, by brute-force reverse-lookup.
  // Note that the LUT computation could be defered somewhere else, and computed once

  static const float T_min = 1700.f;
  static const float T_max = 25000.f;
  static const float T_range = T_max - T_min;
  static const size_t LUT_samples = 1<<16;

  // just init radius with something big.
  float radius = 2.f;
  float temperature = 0.f;

#ifdef _OPENMP
#pragma omp parallel for simd default(none) \
  firstprivate(x, y) shared(radius, temperature)\
  schedule(simd:static)
#endif
  for(size_t i = 0; i < LUT_samples; i++)
  {
    // we need more values for the low temperatures, so we scale the step with a power
    const float step = powf((float)i / (float)(LUT_samples - 1), 4.0f);

    // Current temperature in the lookup range
    const float T = T_min +  step * T_range;

    // Current x, y chromaticity
    float x_bb, y_bb;
    CCT_to_xy_blackbody(T, &x_bb, &y_bb);

    // Compute distance between current planckian chromaticity and input
    float radius_tmp = hypotf((x_bb - x), (y_bb - y));

    // If we found a smaller radius, save it
    const int match = (radius_tmp < radius);
    radius = (match) ? radius_tmp : radius;
    temperature = (match) ? T : temperature;
  }

  return temperature;
}

So, here, I need to share radius and temperature between threads and it is slower than what I would like.

I know I can use reduction(min:radius) if I'm interested in the minimum, so I would like to use a similar reduction here in order to make radius and temperature private in each thread, then at the end, return the temperature correlated to the minimum radius of all threads.

Is that possible ?

Solution

  • Your current code has a nasty race condition in

    // If we found a smaller radius, save it
const int match = (radius_tmp < radius);
radius = (match) ? radius_tmp : radius;
temperature = (match) ? T : temperature;

    Multiple threads can be executing those lines at the same time, causing the values of radius and temperature to get out of sync. It should be:

    #ifdef _OPENMP
#pragma omp critical
#endif
if {radius_tmp < radius) {
  radius = radius_tmp;
  temperature = T;
}

    instead.

    Anyways, OpenMP 4.0 added user-defined reduction operations, so if your compiler supports at least that version, you can try that. Here's an example that uses a struct to wrap multiple values:

    #include <stdio.h>
#include <float.h>

struct pair {
  float radius;
  float temperature;
};

struct pair pair_min(struct pair r, struct pair n) {
  /* r is the current min value, n in the value to compare against it */
  if (n.radius < r.radius) {
    return n;
  } else {
    return r;
  }
}
#ifdef _OPENMP
// Define a new reduction operation
#pragma omp declare reduction \
  (pairmin:struct pair:omp_out=pair_min(omp_out,omp_in))    \
  initializer(omp_priv = { FLT_MAX, 0.0f })
#endif

int main(void) {
  struct pair min_radius = { FLT_MAX, 0.0f };
  struct pair values[4] = {
    {1.0f, 0.1f},
    {2.0f, 0.2f},
    {4.0f, 0.4f},
    {3.0f, 0.3f}
  };

#ifdef _OPENMP
#pragma omp parallel for reduction(pairmin:min_radius)
#endif
  for (int i = 0; i < 4; i++) {
    min_radius = pair_min(min_radius, values[i]);
  }

  printf("{%f, %f}\n", min_radius.radius, min_radius.temperature);
  return 0;
}

    For more information on user-defined reductions, see section 2.19.5.7 of the OpenMP 5.0 specification (or the equivalent in the spec for the version your compiler uses).