[qcms] Keep the output of the TRC between 0 and 1

third_party/qcms/src/transform_util.c

 //  qcms
 //  Copyright (C) 2009 Mozilla Foundation
 //
 // Permission is hereby granted, free of charge, to any person obtaining
 // a copy of this software and associated documentation files (the "Software"),
 // to deal in the Software without restriction, including without limitation
 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
 // and/or sell copies of the Software, and to permit persons to whom the Software
 // is furnished to do so, subject to the following conditions:
 //
@@ skipping 18 matching lines @@
 #include "matrix.h"
 
 #if !defined(INFINITY)
 #define INFINITY HUGE_VAL
 #endif
 
 #define PARAMETRIC_CURVE_TYPE 0x70617261 //'para'
 
 /* value must be a value between 0 and 1 */
 //XXX: is the above a good restriction to have?
-float lut_interp_linear(double value, uint16_t *table, size_t length)
+// the output range of this function is 0..1
+float lut_interp_linear(double input_value, uint16_t *table, size_t length)
 {
         int upper, lower;
-        value = value * (length - 1); // scale to length of the array
-        upper = ceil(value);
-        lower = floor(value);
+        float value;
+        input_value = input_value * (length - 1); // scale to length of the array
+        upper = ceil(input_value);
+        lower = floor(input_value);
         //XXX: can we be more performant here?
-        value = table[upper]*(1. - (upper - value)) + table[lower]*(upper - value);
+        value = table[upper]*(1. - (upper - input_value)) + table[lower]*(upper - input_value);
         /* scale the value */
-        return value * (1./65535.);
+        return value * (1.f/65535.f);
 }
 
 /* same as above but takes and returns a uint16_t value representing a range from 0..1 */
 uint16_t lut_interp_linear16(uint16_t input_value, uint16_t *table, size_t length)
 {
         /* Start scaling input_value to the length of the array: 65535*(length-1).
          * We'll divide out the 65535 next */
         uintptr_t value = (input_value * (length - 1));
         uint32_t upper = (value + 65534) / 65535; /* equivalent to ceil(value/65535) */
         uint32_t lower = value / 65535;           /* equivalent to floor(value/65535) */
@@ skipping 54 matching lines @@
         uint32_t upper = (value + 4095) / 4096; /* equivalent to ceil(value/4096) */
         uint32_t lower = value / 4096;           /* equivalent to floor(value/4096) */
         uint32_t interp = value % 4096;
 
         value = (table[upper]*(interp) + table[lower]*(4096 - interp))/4096; // 0..4096*4096
 
         return value;
 }
 #endif
 
-void compute_curve_gamma_table_type1(float gamma_table[256], double gamma)
+void compute_curve_gamma_table_type1(float gamma_table[256], uint16_t gamma)
+{
+        unsigned int i;
+        float gamma_float = u8Fixed8Number_to_float(gamma);
+        for (i = 0; i < 256; i++) {
+                // 0..1^(0..255 + 255/256) will always be between 0 and 1
+                gamma_table[i] = pow(i/255., gamma_float);

[robert.bradford, 2015/07/01 12:10:50]

I don't think we need to do the calculation in double precision?

so what about:

gamma_table[i] = powf(i/255.f, gamma_float);

[Noel Gordon, 2015/07/02 05:51:49]

Resolved by your subsequent measurements: use pow.  If there's anything
perf-wise we should be expending our time on, it'd be a faster floor.

+        }
+}
+
+void compute_curve_gamma_table_type2(float gamma_table[256], uint16_t *table, size_t length)
 {
         unsigned int i;
         for (i = 0; i < 256; i++) {
-                gamma_table[i] = pow(i/255., gamma);
-        }
-}
-
-void compute_curve_gamma_table_type2(float gamma_table[256], uint16_t *table, int length)
-{
-        unsigned int i;
-        for (i = 0; i < 256; i++) {
                 gamma_table[i] = lut_interp_linear(i/255., table, length);
         }
 }
 
 void compute_curve_gamma_table_type_parametric(float gamma_table[256], float parameter[7], int count)
 {
         size_t X;
         float interval;
         float a, b, c, e, f;
         float y = parameter[0];
@@ skipping 35 matching lines @@
         } else {
                 assert(0 && "invalid parametric function type.");
                 a = 1;
                 b = 0;
                 c = 0;
                 e = 0;
                 f = 0;
                 interval = -INFINITY;
         }       
         for (X = 0; X < 256; X++) {
                 if (X >= interval) {

[Noel Gordon, 2015/07/07 08:00:55]

Oh and spot the bug: interval is a float on 0.0 .. 1.0 and X is a size_t on
[0..255].  I mean, what could possibly go wrong?

                         // XXX The equations are not exactly as definied in the spec but are
                         //     algebraic equivilent.
                         // TODO Should division by 255 be for the whole expression.
-                        gamma_table[X] = pow(a * X / 255. + b, y) + c + e;
+                        gamma_table[X] = clamp_float(pow(a * X / 255. + b, y) + c + e);
                 } else {
-                        gamma_table[X] = c * X / 255. + f;
+                        gamma_table[X] = clamp_float(c * X / 255. + f);
                 }
         }
 }
 
 void compute_curve_gamma_table_type0(float gamma_table[256])
 {
         unsigned int i;
         for (i = 0; i < 256; i++) {
                 gamma_table[i] = i/255.;
         }
 }
 
-
 float clamp_float(float a)
 {
-        if (a > 1.)
-                return 1.;
-        else if (a < 0)
-                return 0;
-        else
-                return a;
+        /* One would naturally write this function as the following:
+        if (a > 1.)
+                return 1.;
+        else if (a < 0)
+                return 0;
+        else
+                return a;
+
+        However, that version will let NaNs pass through which is undesirable
+        for most consumers.
+        */
+
+        if (a > 1.)
+                return 1.;
+        else if (a >= 0)
+                return a;
+        else // a < 0 or a is NaN
+                return 0;
 }
 
 unsigned char clamp_u8(float v)
 {
         if (v > 255.)
                 return 255;
         else if (v < 0)
                 return 0;
         else
                 return floor(v+.5);
@@ skipping 28 matching lines @@
 
         if (!TRC) return NULL;
         gamma_table = malloc(sizeof(float)*256);
         if (gamma_table) {
                 if (TRC->type == PARAMETRIC_CURVE_TYPE) {
                         compute_curve_gamma_table_type_parametric(gamma_table, TRC->parameter, TRC->count);
                 } else {
                         if (TRC->count == 0) {
                                 compute_curve_gamma_table_type0(gamma_table);
                         } else if (TRC->count == 1) {
-                                compute_curve_gamma_table_type1(gamma_table, u8Fixed8Number_to_float(TRC->data[0]));
+                                compute_curve_gamma_table_type1(gamma_table, TRC->data[0]);
                         } else {
                                 compute_curve_gamma_table_type2(gamma_table, TRC->data, TRC->count);
                         }
                 }
         }
 
         validate_gamma_table(gamma_table);
 
         return gamma_table;
 }
@@ skipping 320 matching lines @@
                         //     measurement or data, however it is what lcms uses.
                         *output_gamma_lut_length = trc->count;
                         if (*output_gamma_lut_length < 256)
                                 *output_gamma_lut_length = 256;
 
                         *output_gamma_lut = invert_lut(trc->data, trc->count, *output_gamma_lut_length);
                 }
         }
 
 }