michael@0: #define _ISOC99_SOURCE /* for INFINITY */ michael@0: michael@0: #include michael@0: #include michael@0: #include //memcpy michael@0: #include "qcmsint.h" michael@0: #include "transform_util.h" michael@0: #include "matrix.h" michael@0: michael@0: #if !defined(INFINITY) michael@0: #define INFINITY HUGE_VAL michael@0: #endif michael@0: michael@0: #define PARAMETRIC_CURVE_TYPE 0x70617261 //'para' michael@0: michael@0: /* value must be a value between 0 and 1 */ michael@0: //XXX: is the above a good restriction to have? michael@0: // the output range of this functions is 0..1 michael@0: float lut_interp_linear(double input_value, uint16_t *table, int length) michael@0: { michael@0: int upper, lower; michael@0: float value; michael@0: input_value = input_value * (length - 1); // scale to length of the array michael@0: upper = ceil(input_value); michael@0: lower = floor(input_value); michael@0: //XXX: can we be more performant here? michael@0: value = table[upper]*(1. - (upper - input_value)) + table[lower]*(upper - input_value); michael@0: /* scale the value */ michael@0: return value * (1.f/65535.f); michael@0: } michael@0: michael@0: /* same as above but takes and returns a uint16_t value representing a range from 0..1 */ michael@0: uint16_t lut_interp_linear16(uint16_t input_value, uint16_t *table, int length) michael@0: { michael@0: /* Start scaling input_value to the length of the array: 65535*(length-1). michael@0: * We'll divide out the 65535 next */ michael@0: uint32_t value = (input_value * (length - 1)); michael@0: uint32_t upper = (value + 65534) / 65535; /* equivalent to ceil(value/65535) */ michael@0: uint32_t lower = value / 65535; /* equivalent to floor(value/65535) */ michael@0: /* interp is the distance from upper to value scaled to 0..65535 */ michael@0: uint32_t interp = value % 65535; michael@0: michael@0: value = (table[upper]*(interp) + table[lower]*(65535 - interp))/65535; // 0..65535*65535 michael@0: michael@0: return value; michael@0: } michael@0: michael@0: /* same as above but takes an input_value from 0..PRECACHE_OUTPUT_MAX michael@0: * and returns a uint8_t value representing a range from 0..1 */ michael@0: static michael@0: uint8_t lut_interp_linear_precache_output(uint32_t input_value, uint16_t *table, int length) michael@0: { michael@0: /* Start scaling input_value to the length of the array: PRECACHE_OUTPUT_MAX*(length-1). michael@0: * We'll divide out the PRECACHE_OUTPUT_MAX next */ michael@0: uint32_t value = (input_value * (length - 1)); michael@0: michael@0: /* equivalent to ceil(value/PRECACHE_OUTPUT_MAX) */ michael@0: uint32_t upper = (value + PRECACHE_OUTPUT_MAX-1) / PRECACHE_OUTPUT_MAX; michael@0: /* equivalent to floor(value/PRECACHE_OUTPUT_MAX) */ michael@0: uint32_t lower = value / PRECACHE_OUTPUT_MAX; michael@0: /* interp is the distance from upper to value scaled to 0..PRECACHE_OUTPUT_MAX */ michael@0: uint32_t interp = value % PRECACHE_OUTPUT_MAX; michael@0: michael@0: /* the table values range from 0..65535 */ michael@0: value = (table[upper]*(interp) + table[lower]*(PRECACHE_OUTPUT_MAX - interp)); // 0..(65535*PRECACHE_OUTPUT_MAX) michael@0: michael@0: /* round and scale */ michael@0: value += (PRECACHE_OUTPUT_MAX*65535/255)/2; michael@0: value /= (PRECACHE_OUTPUT_MAX*65535/255); // scale to 0..255 michael@0: return value; michael@0: } michael@0: michael@0: /* value must be a value between 0 and 1 */ michael@0: //XXX: is the above a good restriction to have? michael@0: float lut_interp_linear_float(float value, float *table, int length) michael@0: { michael@0: int upper, lower; michael@0: value = value * (length - 1); michael@0: upper = ceilf(value); michael@0: lower = floorf(value); michael@0: //XXX: can we be more performant here? michael@0: value = table[upper]*(1. - (upper - value)) + table[lower]*(upper - value); michael@0: /* scale the value */ michael@0: return value; michael@0: } michael@0: michael@0: #if 0 michael@0: /* if we use a different representation i.e. one that goes from 0 to 0x1000 we can be more efficient michael@0: * because we can avoid the divisions and use a shifting instead */ michael@0: /* same as above but takes and returns a uint16_t value representing a range from 0..1 */ michael@0: uint16_t lut_interp_linear16(uint16_t input_value, uint16_t *table, int length) michael@0: { michael@0: uint32_t value = (input_value * (length - 1)); michael@0: uint32_t upper = (value + 4095) / 4096; /* equivalent to ceil(value/4096) */ michael@0: uint32_t lower = value / 4096; /* equivalent to floor(value/4096) */ michael@0: uint32_t interp = value % 4096; michael@0: michael@0: value = (table[upper]*(interp) + table[lower]*(4096 - interp))/4096; // 0..4096*4096 michael@0: michael@0: return value; michael@0: } michael@0: #endif michael@0: michael@0: void compute_curve_gamma_table_type1(float gamma_table[256], uint16_t gamma) michael@0: { michael@0: unsigned int i; michael@0: float gamma_float = u8Fixed8Number_to_float(gamma); michael@0: for (i = 0; i < 256; i++) { michael@0: // 0..1^(0..255 + 255/256) will always be between 0 and 1 michael@0: gamma_table[i] = pow(i/255., gamma_float); michael@0: } michael@0: } michael@0: michael@0: void compute_curve_gamma_table_type2(float gamma_table[256], uint16_t *table, int length) michael@0: { michael@0: unsigned int i; michael@0: for (i = 0; i < 256; i++) { michael@0: gamma_table[i] = lut_interp_linear(i/255., table, length); michael@0: } michael@0: } michael@0: michael@0: void compute_curve_gamma_table_type_parametric(float gamma_table[256], float parameter[7], int count) michael@0: { michael@0: size_t X; michael@0: float interval; michael@0: float a, b, c, e, f; michael@0: float y = parameter[0]; michael@0: if (count == 0) { michael@0: a = 1; michael@0: b = 0; michael@0: c = 0; michael@0: e = 0; michael@0: f = 0; michael@0: interval = -INFINITY; michael@0: } else if(count == 1) { michael@0: a = parameter[1]; michael@0: b = parameter[2]; michael@0: c = 0; michael@0: e = 0; michael@0: f = 0; michael@0: interval = -1 * parameter[2] / parameter[1]; michael@0: } else if(count == 2) { michael@0: a = parameter[1]; michael@0: b = parameter[2]; michael@0: c = 0; michael@0: e = parameter[3]; michael@0: f = parameter[3]; michael@0: interval = -1 * parameter[2] / parameter[1]; michael@0: } else if(count == 3) { michael@0: a = parameter[1]; michael@0: b = parameter[2]; michael@0: c = parameter[3]; michael@0: e = -c; michael@0: f = 0; michael@0: interval = parameter[4]; michael@0: } else if(count == 4) { michael@0: a = parameter[1]; michael@0: b = parameter[2]; michael@0: c = parameter[3]; michael@0: e = parameter[5] - c; michael@0: f = parameter[6]; michael@0: interval = parameter[4]; michael@0: } else { michael@0: assert(0 && "invalid parametric function type."); michael@0: a = 1; michael@0: b = 0; michael@0: c = 0; michael@0: e = 0; michael@0: f = 0; michael@0: interval = -INFINITY; michael@0: } michael@0: for (X = 0; X < 256; X++) { michael@0: if (X >= interval) { michael@0: // XXX The equations are not exactly as definied in the spec but are michael@0: // algebraic equivilent. michael@0: // TODO Should division by 255 be for the whole expression. michael@0: gamma_table[X] = clamp_float(pow(a * X / 255. + b, y) + c + e); michael@0: } else { michael@0: gamma_table[X] = clamp_float(c * X / 255. + f); michael@0: } michael@0: } michael@0: } michael@0: michael@0: void compute_curve_gamma_table_type0(float gamma_table[256]) michael@0: { michael@0: unsigned int i; michael@0: for (i = 0; i < 256; i++) { michael@0: gamma_table[i] = i/255.; michael@0: } michael@0: } michael@0: michael@0: float *build_input_gamma_table(struct curveType *TRC) michael@0: { michael@0: float *gamma_table; michael@0: michael@0: if (!TRC) return NULL; michael@0: gamma_table = malloc(sizeof(float)*256); michael@0: if (gamma_table) { michael@0: if (TRC->type == PARAMETRIC_CURVE_TYPE) { michael@0: compute_curve_gamma_table_type_parametric(gamma_table, TRC->parameter, TRC->count); michael@0: } else { michael@0: if (TRC->count == 0) { michael@0: compute_curve_gamma_table_type0(gamma_table); michael@0: } else if (TRC->count == 1) { michael@0: compute_curve_gamma_table_type1(gamma_table, TRC->data[0]); michael@0: } else { michael@0: compute_curve_gamma_table_type2(gamma_table, TRC->data, TRC->count); michael@0: } michael@0: } michael@0: } michael@0: return gamma_table; michael@0: } michael@0: michael@0: struct matrix build_colorant_matrix(qcms_profile *p) michael@0: { michael@0: struct matrix result; michael@0: result.m[0][0] = s15Fixed16Number_to_float(p->redColorant.X); michael@0: result.m[0][1] = s15Fixed16Number_to_float(p->greenColorant.X); michael@0: result.m[0][2] = s15Fixed16Number_to_float(p->blueColorant.X); michael@0: result.m[1][0] = s15Fixed16Number_to_float(p->redColorant.Y); michael@0: result.m[1][1] = s15Fixed16Number_to_float(p->greenColorant.Y); michael@0: result.m[1][2] = s15Fixed16Number_to_float(p->blueColorant.Y); michael@0: result.m[2][0] = s15Fixed16Number_to_float(p->redColorant.Z); michael@0: result.m[2][1] = s15Fixed16Number_to_float(p->greenColorant.Z); michael@0: result.m[2][2] = s15Fixed16Number_to_float(p->blueColorant.Z); michael@0: result.invalid = false; michael@0: return result; michael@0: } michael@0: michael@0: /* The following code is copied nearly directly from lcms. michael@0: * I think it could be much better. For example, Argyll seems to have better code in michael@0: * icmTable_lookup_bwd and icmTable_setup_bwd. However, for now this is a quick way michael@0: * to a working solution and allows for easy comparing with lcms. */ michael@0: uint16_fract_t lut_inverse_interp16(uint16_t Value, uint16_t LutTable[], int length) michael@0: { michael@0: int l = 1; michael@0: int r = 0x10000; michael@0: int x = 0, res; // 'int' Give spacing for negative values michael@0: int NumZeroes, NumPoles; michael@0: int cell0, cell1; michael@0: double val2; michael@0: double y0, y1, x0, x1; michael@0: double a, b, f; michael@0: michael@0: // July/27 2001 - Expanded to handle degenerated curves with an arbitrary michael@0: // number of elements containing 0 at the begining of the table (Zeroes) michael@0: // and another arbitrary number of poles (FFFFh) at the end. michael@0: // First the zero and pole extents are computed, then value is compared. michael@0: michael@0: NumZeroes = 0; michael@0: while (LutTable[NumZeroes] == 0 && NumZeroes < length-1) michael@0: NumZeroes++; michael@0: michael@0: // There are no zeros at the beginning and we are trying to find a zero, so michael@0: // return anything. It seems zero would be the less destructive choice michael@0: /* I'm not sure that this makes sense, but oh well... */ michael@0: if (NumZeroes == 0 && Value == 0) michael@0: return 0; michael@0: michael@0: NumPoles = 0; michael@0: while (LutTable[length-1- NumPoles] == 0xFFFF && NumPoles < length-1) michael@0: NumPoles++; michael@0: michael@0: // Does the curve belong to this case? michael@0: if (NumZeroes > 1 || NumPoles > 1) michael@0: { michael@0: int a, b; michael@0: michael@0: // Identify if value fall downto 0 or FFFF zone michael@0: if (Value == 0) return 0; michael@0: // if (Value == 0xFFFF) return 0xFFFF; michael@0: michael@0: // else restrict to valid zone michael@0: michael@0: a = ((NumZeroes-1) * 0xFFFF) / (length-1); michael@0: b = ((length-1 - NumPoles) * 0xFFFF) / (length-1); michael@0: michael@0: l = a - 1; michael@0: r = b + 1; michael@0: } michael@0: michael@0: michael@0: // Seems not a degenerated case... apply binary search michael@0: michael@0: while (r > l) { michael@0: michael@0: x = (l + r) / 2; michael@0: michael@0: res = (int) lut_interp_linear16((uint16_fract_t) (x-1), LutTable, length); michael@0: michael@0: if (res == Value) { michael@0: michael@0: // Found exact match. michael@0: michael@0: return (uint16_fract_t) (x - 1); michael@0: } michael@0: michael@0: if (res > Value) r = x - 1; michael@0: else l = x + 1; michael@0: } michael@0: michael@0: // Not found, should we interpolate? michael@0: michael@0: michael@0: // Get surrounding nodes michael@0: michael@0: val2 = (length-1) * ((double) (x - 1) / 65535.0); michael@0: michael@0: cell0 = (int) floor(val2); michael@0: cell1 = (int) ceil(val2); michael@0: michael@0: if (cell0 == cell1) return (uint16_fract_t) x; michael@0: michael@0: y0 = LutTable[cell0] ; michael@0: x0 = (65535.0 * cell0) / (length-1); michael@0: michael@0: y1 = LutTable[cell1] ; michael@0: x1 = (65535.0 * cell1) / (length-1); michael@0: michael@0: a = (y1 - y0) / (x1 - x0); michael@0: b = y0 - a * x0; michael@0: michael@0: if (fabs(a) < 0.01) return (uint16_fract_t) x; michael@0: michael@0: f = ((Value - b) / a); michael@0: michael@0: if (f < 0.0) return (uint16_fract_t) 0; michael@0: if (f >= 65535.0) return (uint16_fract_t) 0xFFFF; michael@0: michael@0: return (uint16_fract_t) floor(f + 0.5); michael@0: michael@0: } michael@0: michael@0: /* michael@0: The number of entries needed to invert a lookup table should not michael@0: necessarily be the same as the original number of entries. This is michael@0: especially true of lookup tables that have a small number of entries. michael@0: michael@0: For example: michael@0: Using a table like: michael@0: {0, 3104, 14263, 34802, 65535} michael@0: invert_lut will produce an inverse of: michael@0: {3, 34459, 47529, 56801, 65535} michael@0: which has an maximum error of about 9855 (pixel difference of ~38.346) michael@0: michael@0: For now, we punt the decision of output size to the caller. */ michael@0: static uint16_t *invert_lut(uint16_t *table, int length, int out_length) michael@0: { michael@0: int i; michael@0: /* for now we invert the lut by creating a lut of size out_length michael@0: * and attempting to lookup a value for each entry using lut_inverse_interp16 */ michael@0: uint16_t *output = malloc(sizeof(uint16_t)*out_length); michael@0: if (!output) michael@0: return NULL; michael@0: michael@0: for (i = 0; i < out_length; i++) { michael@0: double x = ((double) i * 65535.) / (double) (out_length - 1); michael@0: uint16_fract_t input = floor(x + .5); michael@0: output[i] = lut_inverse_interp16(input, table, length); michael@0: } michael@0: return output; michael@0: } michael@0: michael@0: static void compute_precache_pow(uint8_t *output, float gamma) michael@0: { michael@0: uint32_t v = 0; michael@0: for (v = 0; v < PRECACHE_OUTPUT_SIZE; v++) { michael@0: //XXX: don't do integer/float conversion... and round? michael@0: output[v] = 255. * pow(v/(double)PRECACHE_OUTPUT_MAX, gamma); michael@0: } michael@0: } michael@0: michael@0: void compute_precache_lut(uint8_t *output, uint16_t *table, int length) michael@0: { michael@0: uint32_t v = 0; michael@0: for (v = 0; v < PRECACHE_OUTPUT_SIZE; v++) { michael@0: output[v] = lut_interp_linear_precache_output(v, table, length); michael@0: } michael@0: } michael@0: michael@0: void compute_precache_linear(uint8_t *output) michael@0: { michael@0: uint32_t v = 0; michael@0: for (v = 0; v < PRECACHE_OUTPUT_SIZE; v++) { michael@0: //XXX: round? michael@0: output[v] = v / (PRECACHE_OUTPUT_SIZE/256); michael@0: } michael@0: } michael@0: michael@0: qcms_bool compute_precache(struct curveType *trc, uint8_t *output) michael@0: { michael@0: michael@0: if (trc->type == PARAMETRIC_CURVE_TYPE) { michael@0: float gamma_table[256]; michael@0: uint16_t gamma_table_uint[256]; michael@0: uint16_t i; michael@0: uint16_t *inverted; michael@0: int inverted_size = 256; michael@0: michael@0: compute_curve_gamma_table_type_parametric(gamma_table, trc->parameter, trc->count); michael@0: for(i = 0; i < 256; i++) { michael@0: gamma_table_uint[i] = (uint16_t)(gamma_table[i] * 65535); michael@0: } michael@0: michael@0: //XXX: the choice of a minimum of 256 here is not backed by any theory, michael@0: // measurement or data, howeve r it is what lcms uses. michael@0: // the maximum number we would need is 65535 because that's the michael@0: // accuracy used for computing the pre cache table michael@0: if (inverted_size < 256) michael@0: inverted_size = 256; michael@0: michael@0: inverted = invert_lut(gamma_table_uint, 256, inverted_size); michael@0: if (!inverted) michael@0: return false; michael@0: compute_precache_lut(output, inverted, inverted_size); michael@0: free(inverted); michael@0: } else { michael@0: if (trc->count == 0) { michael@0: compute_precache_linear(output); michael@0: } else if (trc->count == 1) { michael@0: compute_precache_pow(output, 1./u8Fixed8Number_to_float(trc->data[0])); michael@0: } else { michael@0: uint16_t *inverted; michael@0: int inverted_size = trc->count; michael@0: //XXX: the choice of a minimum of 256 here is not backed by any theory, michael@0: // measurement or data, howeve r it is what lcms uses. michael@0: // the maximum number we would need is 65535 because that's the michael@0: // accuracy used for computing the pre cache table michael@0: if (inverted_size < 256) michael@0: inverted_size = 256; michael@0: michael@0: inverted = invert_lut(trc->data, trc->count, inverted_size); michael@0: if (!inverted) michael@0: return false; michael@0: compute_precache_lut(output, inverted, inverted_size); michael@0: free(inverted); michael@0: } michael@0: } michael@0: return true; michael@0: } michael@0: michael@0: michael@0: static uint16_t *build_linear_table(int length) michael@0: { michael@0: int i; michael@0: uint16_t *output = malloc(sizeof(uint16_t)*length); michael@0: if (!output) michael@0: return NULL; michael@0: michael@0: for (i = 0; i < length; i++) { michael@0: double x = ((double) i * 65535.) / (double) (length - 1); michael@0: uint16_fract_t input = floor(x + .5); michael@0: output[i] = input; michael@0: } michael@0: return output; michael@0: } michael@0: michael@0: static uint16_t *build_pow_table(float gamma, int length) michael@0: { michael@0: int i; michael@0: uint16_t *output = malloc(sizeof(uint16_t)*length); michael@0: if (!output) michael@0: return NULL; michael@0: michael@0: for (i = 0; i < length; i++) { michael@0: uint16_fract_t result; michael@0: double x = ((double) i) / (double) (length - 1); michael@0: x = pow(x, gamma); //XXX turn this conversion into a function michael@0: result = floor(x*65535. + .5); michael@0: output[i] = result; michael@0: } michael@0: return output; michael@0: } michael@0: michael@0: void build_output_lut(struct curveType *trc, michael@0: uint16_t **output_gamma_lut, size_t *output_gamma_lut_length) michael@0: { michael@0: if (trc->type == PARAMETRIC_CURVE_TYPE) { michael@0: float gamma_table[256]; michael@0: uint16_t i; michael@0: uint16_t *output = malloc(sizeof(uint16_t)*256); michael@0: michael@0: if (!output) { michael@0: *output_gamma_lut = NULL; michael@0: return; michael@0: } michael@0: michael@0: compute_curve_gamma_table_type_parametric(gamma_table, trc->parameter, trc->count); michael@0: *output_gamma_lut_length = 256; michael@0: for(i = 0; i < 256; i++) { michael@0: output[i] = (uint16_t)(gamma_table[i] * 65535); michael@0: } michael@0: *output_gamma_lut = output; michael@0: } else { michael@0: if (trc->count == 0) { michael@0: *output_gamma_lut = build_linear_table(4096); michael@0: *output_gamma_lut_length = 4096; michael@0: } else if (trc->count == 1) { michael@0: float gamma = 1./u8Fixed8Number_to_float(trc->data[0]); michael@0: *output_gamma_lut = build_pow_table(gamma, 4096); michael@0: *output_gamma_lut_length = 4096; michael@0: } else { michael@0: //XXX: the choice of a minimum of 256 here is not backed by any theory, michael@0: // measurement or data, however it is what lcms uses. michael@0: *output_gamma_lut_length = trc->count; michael@0: if (*output_gamma_lut_length < 256) michael@0: *output_gamma_lut_length = 256; michael@0: michael@0: *output_gamma_lut = invert_lut(trc->data, trc->count, *output_gamma_lut_length); michael@0: } michael@0: } michael@0: michael@0: } michael@0: