diff options
author | Jeff Muizelaar <jmuizelaar@mozilla.com> | 2010-07-14 14:03:58 -0400 |
---|---|---|
committer | Jeff Muizelaar <jmuizelaar@mozilla.com> | 2010-07-14 14:03:58 -0400 |
commit | 8b943d9ce6451079fd236863ee37613950332182 (patch) | |
tree | 292892cd5482ae7285757408d8633dda0b5d3730 | |
parent | 86771fe9fc6c19856af9316b7269be8cdf4e126e (diff) |
Add missing files
-rw-r--r-- | chain.c | 988 | ||||
-rw-r--r-- | chain.h | 30 | ||||
-rw-r--r-- | matrix.c | 136 | ||||
-rw-r--r-- | matrix.h | 39 | ||||
-rw-r--r-- | transform_util.c | 534 | ||||
-rw-r--r-- | transform_util.h | 59 |
6 files changed, 1786 insertions, 0 deletions
@@ -0,0 +1,988 @@ +/* vim: set ts=8 sw=8 noexpandtab: */ +// qcms +// Copyright (C) 2009 Mozilla Corporation +// Copyright (C) 1998-2007 Marti Maria +// +// 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: +// +// The above copyright notice and this permission notice shall be included in +// all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO +// THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE +// LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION +// OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION +// WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + +#include <stdlib.h> +#include <math.h> +#include <assert.h> +#include <string.h> //memcpy +#include "qcmsint.h" +#include "transform_util.h" +#include "matrix.h" + +static struct matrix build_lut_matrix(struct lutType *lut) +{ + struct matrix result; + if (lut) { + result.m[0][0] = s15Fixed16Number_to_float(lut->e00); + result.m[0][1] = s15Fixed16Number_to_float(lut->e01); + result.m[0][2] = s15Fixed16Number_to_float(lut->e02); + result.m[1][0] = s15Fixed16Number_to_float(lut->e10); + result.m[1][1] = s15Fixed16Number_to_float(lut->e11); + result.m[1][2] = s15Fixed16Number_to_float(lut->e12); + result.m[2][0] = s15Fixed16Number_to_float(lut->e20); + result.m[2][1] = s15Fixed16Number_to_float(lut->e21); + result.m[2][2] = s15Fixed16Number_to_float(lut->e22); + result.invalid = false; + } else { + memset(&result, 0, sizeof(struct matrix)); + result.invalid = true; + } + return result; +} + +static struct matrix build_mAB_matrix(struct lutmABType *lut) +{ + struct matrix result; + if (lut) { + result.m[0][0] = s15Fixed16Number_to_float(lut->e00); + result.m[0][1] = s15Fixed16Number_to_float(lut->e01); + result.m[0][2] = s15Fixed16Number_to_float(lut->e02); + result.m[1][0] = s15Fixed16Number_to_float(lut->e10); + result.m[1][1] = s15Fixed16Number_to_float(lut->e11); + result.m[1][2] = s15Fixed16Number_to_float(lut->e12); + result.m[2][0] = s15Fixed16Number_to_float(lut->e20); + result.m[2][1] = s15Fixed16Number_to_float(lut->e21); + result.m[2][2] = s15Fixed16Number_to_float(lut->e22); + result.invalid = false; + } else { + memset(&result, 0, sizeof(struct matrix)); + result.invalid = true; + } + return result; +} + +//Based on lcms cmsLab2XYZ +#define f(t) (t <= (24.0f/116.0f)*(24.0f/116.0f)*(24.0f/116.0f)) ? ((841.0/108.0) * t + (16.0/116.0)) : pow(t,1.0/3.0) +#define f_1(t) (t <= (24.0f/116.0f)) ? ((108.0/841.0) * (t - (16.0/116.0))) : (t * t * t) +static void qcms_transform_module_LAB_to_XYZ(struct qcms_modular_transform *transform, float *src, float *dest, size_t length) +{ + size_t i; + // lcms: D50 XYZ values + float WhitePointX = 0.9642f; + float WhitePointY = 1.0f; + float WhitePointZ = 0.8249f; + for (i = 0; i < length; i++) { + float device_L = *src++ * 100.0f; + float device_a = *src++ * 255.0f - 128.0f; + float device_b = *src++ * 255.0f - 128.0f; + float y = (device_L + 16.0f) / 116.0f; + + float X = f_1((y + 0.002f * device_a)) * WhitePointX; + float Y = f_1(y) * WhitePointY; + float Z = f_1((y - 0.005f * device_b)) * WhitePointZ; + *dest++ = X / (1.0 + 32767.0/32768.0); + *dest++ = Y / (1.0 + 32767.0/32768.0); + *dest++ = Z / (1.0 + 32767.0/32768.0); + } +} + +//Based on lcms cmsXYZ2Lab +static void qcms_transform_module_XYZ_to_LAB(struct qcms_modular_transform *transform, float *src, float *dest, size_t length) +{ + size_t i; + // lcms: D50 XYZ values + float WhitePointX = 0.9642f; + float WhitePointY = 1.0f; + float WhitePointZ = 0.8249f; + for (i = 0; i < length; i++) { + float device_x = *src++ * (1.0 + 32767.0/32768.0) / WhitePointX; + float device_y = *src++ * (1.0 + 32767.0/32768.0) / WhitePointY; + float device_z = *src++ * (1.0 + 32767.0/32768.0) / WhitePointZ; + + float fx = f(device_x); + float fy = f(device_y); + float fz = f(device_z); + + float L = 116.0f*fy - 16.0f; + float a = 500.0f*(fx - fy); + float b = 200.0f*(fy - fz); + *dest++ = L / 100.0f; + *dest++ = (a+128.0f) / 255.0f; + *dest++ = (b+128.0f) / 255.0f; + } + +} + +static void qcms_transform_module_clut_only(struct qcms_modular_transform *transform, float *src, float *dest, size_t length) +{ + size_t i; + int xy_len = 1; + int x_len = transform->grid_size; + int len = x_len * x_len; + float* r_table = transform->r_clut; + float* g_table = transform->g_clut; + float* b_table = transform->b_clut; + + for (i = 0; i < length; i++) { + float linear_r = *src++; + float linear_g = *src++; + float linear_b = *src++; + + int x = floor(linear_r * (transform->grid_size-1)); + int y = floor(linear_g * (transform->grid_size-1)); + int z = floor(linear_b * (transform->grid_size-1)); + int x_n = ceil(linear_r * (transform->grid_size-1)); + int y_n = ceil(linear_g * (transform->grid_size-1)); + int z_n = ceil(linear_b * (transform->grid_size-1)); + float x_d = linear_r * (transform->grid_size-1) - x; + float y_d = linear_g * (transform->grid_size-1) - y; + float z_d = linear_b * (transform->grid_size-1) - z; + + float r_x1 = lerp(CLU(r_table,x,y,z), CLU(r_table,x_n,y,z), x_d); + float r_x2 = lerp(CLU(r_table,x,y_n,z), CLU(r_table,x_n,y_n,z), x_d); + float r_y1 = lerp(r_x1, r_x2, y_d); + float r_x3 = lerp(CLU(r_table,x,y,z_n), CLU(r_table,x_n,y,z_n), x_d); + float r_x4 = lerp(CLU(r_table,x,y_n,z_n), CLU(r_table,x_n,y_n,z_n), x_d); + float r_y2 = lerp(r_x3, r_x4, y_d); + float clut_r = lerp(r_y1, r_y2, z_d); + + float g_x1 = lerp(CLU(g_table,x,y,z), CLU(g_table,x_n,y,z), x_d); + float g_x2 = lerp(CLU(g_table,x,y_n,z), CLU(g_table,x_n,y_n,z), x_d); + float g_y1 = lerp(g_x1, g_x2, y_d); + float g_x3 = lerp(CLU(g_table,x,y,z_n), CLU(g_table,x_n,y,z_n), x_d); + float g_x4 = lerp(CLU(g_table,x,y_n,z_n), CLU(g_table,x_n,y_n,z_n), x_d); + float g_y2 = lerp(g_x3, g_x4, y_d); + float clut_g = lerp(g_y1, g_y2, z_d); + + float b_x1 = lerp(CLU(b_table,x,y,z), CLU(b_table,x_n,y,z), x_d); + float b_x2 = lerp(CLU(b_table,x,y_n,z), CLU(b_table,x_n,y_n,z), x_d); + float b_y1 = lerp(b_x1, b_x2, y_d); + float b_x3 = lerp(CLU(b_table,x,y,z_n), CLU(b_table,x_n,y,z_n), x_d); + float b_x4 = lerp(CLU(b_table,x,y_n,z_n), CLU(b_table,x_n,y_n,z_n), x_d); + float b_y2 = lerp(b_x3, b_x4, y_d); + float clut_b = lerp(b_y1, b_y2, z_d); + + *dest++ = clamp_float(clut_r); + *dest++ = clamp_float(clut_g); + *dest++ = clamp_float(clut_b); + } +} + +static void qcms_transform_module_clut(struct qcms_modular_transform *transform, float *src, float *dest, size_t length) +{ + size_t i; + int xy_len = 1; + int x_len = transform->grid_size; + int len = x_len * x_len; + float* r_table = transform->r_clut; + float* g_table = transform->g_clut; + float* b_table = transform->b_clut; + for (i = 0; i < length; i++) { + float device_r = *src++; + float device_g = *src++; + float device_b = *src++; + float linear_r = lut_interp_linear_float(device_r, + transform->input_clut_table_r, transform->input_clut_table_length); + float linear_g = lut_interp_linear_float(device_g, + transform->input_clut_table_g, transform->input_clut_table_length); + float linear_b = lut_interp_linear_float(device_b, + transform->input_clut_table_b, transform->input_clut_table_length); + + int x = floor(linear_r * (transform->grid_size-1)); + int y = floor(linear_g * (transform->grid_size-1)); + int z = floor(linear_b * (transform->grid_size-1)); + int x_n = ceil(linear_r * (transform->grid_size-1)); + int y_n = ceil(linear_g * (transform->grid_size-1)); + int z_n = ceil(linear_b * (transform->grid_size-1)); + float x_d = linear_r * (transform->grid_size-1) - x; + float y_d = linear_g * (transform->grid_size-1) - y; + float z_d = linear_b * (transform->grid_size-1) - z; + + float r_x1 = lerp(CLU(r_table,x,y,z), CLU(r_table,x_n,y,z), x_d); + float r_x2 = lerp(CLU(r_table,x,y_n,z), CLU(r_table,x_n,y_n,z), x_d); + float r_y1 = lerp(r_x1, r_x2, y_d); + float r_x3 = lerp(CLU(r_table,x,y,z_n), CLU(r_table,x_n,y,z_n), x_d); + float r_x4 = lerp(CLU(r_table,x,y_n,z_n), CLU(r_table,x_n,y_n,z_n), x_d); + float r_y2 = lerp(r_x3, r_x4, y_d); + float clut_r = lerp(r_y1, r_y2, z_d); + + float g_x1 = lerp(CLU(g_table,x,y,z), CLU(g_table,x_n,y,z), x_d); + float g_x2 = lerp(CLU(g_table,x,y_n,z), CLU(g_table,x_n,y_n,z), x_d); + float g_y1 = lerp(g_x1, g_x2, y_d); + float g_x3 = lerp(CLU(g_table,x,y,z_n), CLU(g_table,x_n,y,z_n), x_d); + float g_x4 = lerp(CLU(g_table,x,y_n,z_n), CLU(g_table,x_n,y_n,z_n), x_d); + float g_y2 = lerp(g_x3, g_x4, y_d); + float clut_g = lerp(g_y1, g_y2, z_d); + + float b_x1 = lerp(CLU(b_table,x,y,z), CLU(b_table,x_n,y,z), x_d); + float b_x2 = lerp(CLU(b_table,x,y_n,z), CLU(b_table,x_n,y_n,z), x_d); + float b_y1 = lerp(b_x1, b_x2, y_d); + float b_x3 = lerp(CLU(b_table,x,y,z_n), CLU(b_table,x_n,y,z_n), x_d); + float b_x4 = lerp(CLU(b_table,x,y_n,z_n), CLU(b_table,x_n,y_n,z_n), x_d); + float b_y2 = lerp(b_x3, b_x4, y_d); + float clut_b = lerp(b_y1, b_y2, z_d); + + float pcs_r = lut_interp_linear_float(clut_r, + transform->output_clut_table_r, transform->output_clut_table_length); + float pcs_g = lut_interp_linear_float(clut_g, + transform->output_clut_table_g, transform->output_clut_table_length); + float pcs_b = lut_interp_linear_float(clut_b, + transform->output_clut_table_b, transform->output_clut_table_length); + + *dest++ = clamp_float(pcs_r); + *dest++ = clamp_float(pcs_g); + *dest++ = clamp_float(pcs_b); + } +} + +/* NOT USED +static void qcms_transform_module_tetra_clut(struct qcms_modular_transform *transform, float *src, float *dest, size_t length) +{ + size_t i; + int xy_len = 1; + int x_len = transform->grid_size; + int len = x_len * x_len; + float* r_table = transform->r_clut; + float* g_table = transform->g_clut; + float* b_table = transform->b_clut; + float c0_r, c1_r, c2_r, c3_r; + float c0_g, c1_g, c2_g, c3_g; + float c0_b, c1_b, c2_b, c3_b; + float clut_r, clut_g, clut_b; + float pcs_r, pcs_g, pcs_b; + for (i = 0; i < length; i++) { + float device_r = *src++; + float device_g = *src++; + float device_b = *src++; + float linear_r = lut_interp_linear_float(device_r, + transform->input_clut_table_r, transform->input_clut_table_length); + float linear_g = lut_interp_linear_float(device_g, + transform->input_clut_table_g, transform->input_clut_table_length); + float linear_b = lut_interp_linear_float(device_b, + transform->input_clut_table_b, transform->input_clut_table_length); + + int x = floor(linear_r * (transform->grid_size-1)); + int y = floor(linear_g * (transform->grid_size-1)); + int z = floor(linear_b * (transform->grid_size-1)); + int x_n = ceil(linear_r * (transform->grid_size-1)); + int y_n = ceil(linear_g * (transform->grid_size-1)); + int z_n = ceil(linear_b * (transform->grid_size-1)); + float rx = linear_r * (transform->grid_size-1) - x; + float ry = linear_g * (transform->grid_size-1) - y; + float rz = linear_b * (transform->grid_size-1) - z; + + c0_r = CLU(r_table, x, y, z); + c0_g = CLU(g_table, x, y, z); + c0_b = CLU(b_table, x, y, z); + if( rx >= ry ) { + if (ry >= rz) { //rx >= ry && ry >= rz + c1_r = CLU(r_table, x_n, y, z) - c0_r; + c2_r = CLU(r_table, x_n, y_n, z) - CLU(r_table, x_n, y, z); + c3_r = CLU(r_table, x_n, y_n, z_n) - CLU(r_table, x_n, y_n, z); + c1_g = CLU(g_table, x_n, y, z) - c0_g; + c2_g = CLU(g_table, x_n, y_n, z) - CLU(g_table, x_n, y, z); + c3_g = CLU(g_table, x_n, y_n, z_n) - CLU(g_table, x_n, y_n, z); + c1_b = CLU(b_table, x_n, y, z) - c0_b; + c2_b = CLU(b_table, x_n, y_n, z) - CLU(b_table, x_n, y, z); + c3_b = CLU(b_table, x_n, y_n, z_n) - CLU(b_table, x_n, y_n, z); + } else { + if (rx >= rz) { //rx >= rz && rz >= ry + c1_r = CLU(r_table, x_n, y, z) - c0_r; + c2_r = CLU(r_table, x_n, y_n, z_n) - CLU(r_table, x_n, y, z_n); + c3_r = CLU(r_table, x_n, y, z_n) - CLU(r_table, x_n, y, z); + c1_g = CLU(g_table, x_n, y, z) - c0_g; + c2_g = CLU(g_table, x_n, y_n, z_n) - CLU(g_table, x_n, y, z_n); + c3_g = CLU(g_table, x_n, y, z_n) - CLU(g_table, x_n, y, z); + c1_b = CLU(b_table, x_n, y, z) - c0_b; + c2_b = CLU(b_table, x_n, y_n, z_n) - CLU(b_table, x_n, y, z_n); + c3_b = CLU(b_table, x_n, y, z_n) - CLU(b_table, x_n, y, z); + } else { //rz > rx && rx >= ry + c1_r = CLU(r_table, x_n, y, z_n) - CLU(r_table, x, y, z_n); + c2_r = CLU(r_table, x_n, y_n, z_n) - CLU(r_table, x_n, y, z_n); + c3_r = CLU(r_table, x, y, z_n) - c0_r; + c1_g = CLU(g_table, x_n, y, z_n) - CLU(g_table, x, y, z_n); + c2_g = CLU(g_table, x_n, y_n, z_n) - CLU(g_table, x_n, y, z_n); + c3_g = CLU(g_table, x, y, z_n) - c0_g; + c1_b = CLU(b_table, x_n, y, z_n) - CLU(b_table, x, y, z_n); + c2_b = CLU(b_table, x_n, y_n, z_n) - CLU(b_table, x_n, y, z_n); + c3_b = CLU(b_table, x, y, z_n) - c0_b; + } + } + } else { + if (rx >= rz) { //ry > rx && rx >= rz + c1_r = CLU(r_table, x_n, y_n, z) - CLU(r_table, x, y_n, z); + c2_r = CLU(r_table, x_n, y_n, z) - c0_r; + c3_r = CLU(r_table, x_n, y_n, z_n) - CLU(r_table, x_n, y_n, z); + c1_g = CLU(g_table, x_n, y_n, z) - CLU(g_table, x, y_n, z); + c2_g = CLU(g_table, x_n, y_n, z) - c0_g; + c3_g = CLU(g_table, x_n, y_n, z_n) - CLU(g_table, x_n, y_n, z); + c1_b = CLU(b_table, x_n, y_n, z) - CLU(b_table, x, y_n, z); + c2_b = CLU(b_table, x_n, y_n, z) - c0_b; + c3_b = CLU(b_table, x_n, y_n, z_n) - CLU(b_table, x_n, y_n, z); + } else { + if (ry >= rz) { //ry >= rz && rz > rx + c1_r = CLU(r_table, x_n, y_n, z_n) - CLU(r_table, x, y_n, z_n); + c2_r = CLU(r_table, x, y_n, z) - c0_r; + c3_r = CLU(r_table, x, y_n, z_n) - CLU(r_table, x, y_n, z); + c1_g = CLU(g_table, x_n, y_n, z_n) - CLU(g_table, x, y_n, z_n); + c2_g = CLU(g_table, x, y_n, z) - c0_g; + c3_g = CLU(g_table, x, y_n, z_n) - CLU(g_table, x, y_n, z); + c1_b = CLU(b_table, x_n, y_n, z_n) - CLU(b_table, x, y_n, z_n); + c2_b = CLU(b_table, x, y_n, z) - c0_b; + c3_b = CLU(b_table, x, y_n, z_n) - CLU(b_table, x, y_n, z); + } else { //rz > ry && ry > rx + c1_r = CLU(r_table, x_n, y_n, z_n) - CLU(r_table, x, y_n, z_n); + c2_r = CLU(r_table, x, y_n, z) - c0_r; + c3_r = CLU(r_table, x_n, y_n, z_n) - CLU(r_table, x_n, y_n, z); + c1_g = CLU(g_table, x_n, y_n, z_n) - CLU(g_table, x, y_n, z_n); + c2_g = CLU(g_table, x, y_n, z) - c0_g; + c3_g = CLU(g_table, x_n, y_n, z_n) - CLU(g_table, x_n, y_n, z); + c1_b = CLU(b_table, x_n, y_n, z_n) - CLU(b_table, x, y_n, z_n); + c2_b = CLU(b_table, x, y_n, z) - c0_b; + c3_b = CLU(b_table, x_n, y_n, z_n) - CLU(b_table, x_n, y_n, z); + } + } + } + + clut_r = c0_r + c1_r*rx + c2_r*ry + c3_r*rz; + clut_g = c0_g + c1_g*rx + c2_g*ry + c3_g*rz; + clut_b = c0_b + c1_b*rx + c2_b*ry + c3_b*rz; + + pcs_r = lut_interp_linear_float(clut_r, + transform->output_clut_table_r, transform->output_clut_table_length); + pcs_g = lut_interp_linear_float(clut_g, + transform->output_clut_table_g, transform->output_clut_table_length); + pcs_b = lut_interp_linear_float(clut_b, + transform->output_clut_table_b, transform->output_clut_table_length); + *dest++ = clamp_float(pcs_r); + *dest++ = clamp_float(pcs_g); + *dest++ = clamp_float(pcs_b); + } +} +*/ + +static void qcms_transform_module_gamma_table(struct qcms_modular_transform *transform, float *src, float *dest, size_t length) +{ + size_t i; + float out_r, out_g, out_b; + for (i = 0; i < length; i++) { + float in_r = *src++; + float in_g = *src++; + float in_b = *src++; + + out_r = lut_interp_linear_float(in_r, transform->input_clut_table_r, 256); + out_g = lut_interp_linear_float(in_g, transform->input_clut_table_g, 256); + out_b = lut_interp_linear_float(in_b, transform->input_clut_table_b, 256); + + *dest++ = clamp_float(out_r); + *dest++ = clamp_float(out_g); + *dest++ = clamp_float(out_b); + } +} + +static void qcms_transform_module_gamma_lut(struct qcms_modular_transform *transform, float *src, float *dest, size_t length) +{ + size_t i; + float out_r, out_g, out_b; + for (i = 0; i < length; i++) { + float in_r = *src++; + float in_g = *src++; + float in_b = *src++; + + out_r = lut_interp_linear(in_r, + transform->output_gamma_lut_r, transform->output_gamma_lut_r_length); + out_g = lut_interp_linear(in_g, + transform->output_gamma_lut_g, transform->output_gamma_lut_g_length); + out_b = lut_interp_linear(in_b, + transform->output_gamma_lut_b, transform->output_gamma_lut_b_length); + + *dest++ = clamp_float(out_r); + *dest++ = clamp_float(out_g); + *dest++ = clamp_float(out_b); + } +} + +static void qcms_transform_module_matrix_translate(struct qcms_modular_transform *transform, float *src, float *dest, size_t length) +{ + size_t i; + struct matrix mat; + + /* store the results in column major mode + * this makes doing the multiplication with sse easier */ + mat.m[0][0] = transform->matrix.m[0][0]; + mat.m[1][0] = transform->matrix.m[0][1]; + mat.m[2][0] = transform->matrix.m[0][2]; + mat.m[0][1] = transform->matrix.m[1][0]; + mat.m[1][1] = transform->matrix.m[1][1]; + mat.m[2][1] = transform->matrix.m[1][2]; + mat.m[0][2] = transform->matrix.m[2][0]; + mat.m[1][2] = transform->matrix.m[2][1]; + mat.m[2][2] = transform->matrix.m[2][2]; + + for (i = 0; i < length; i++) { + float in_r = *src++; + float in_g = *src++; + float in_b = *src++; + + float out_r = mat.m[0][0]*in_r + mat.m[1][0]*in_g + mat.m[2][0]*in_b + transform->tx; + float out_g = mat.m[0][1]*in_r + mat.m[1][1]*in_g + mat.m[2][1]*in_b + transform->ty; + float out_b = mat.m[0][2]*in_r + mat.m[1][2]*in_g + mat.m[2][2]*in_b + transform->tz; + + *dest++ = clamp_float(out_r); + *dest++ = clamp_float(out_g); + *dest++ = clamp_float(out_b); + } +} + +static void qcms_transform_module_matrix(struct qcms_modular_transform *transform, float *src, float *dest, size_t length) +{ + size_t i; + struct matrix mat; + + /* store the results in column major mode + * this makes doing the multiplication with sse easier */ + mat.m[0][0] = transform->matrix.m[0][0]; + mat.m[1][0] = transform->matrix.m[0][1]; + mat.m[2][0] = transform->matrix.m[0][2]; + mat.m[0][1] = transform->matrix.m[1][0]; + mat.m[1][1] = transform->matrix.m[1][1]; + mat.m[2][1] = transform->matrix.m[1][2]; + mat.m[0][2] = transform->matrix.m[2][0]; + mat.m[1][2] = transform->matrix.m[2][1]; + mat.m[2][2] = transform->matrix.m[2][2]; + + for (i = 0; i < length; i++) { + float in_r = *src++; + float in_g = *src++; + float in_b = *src++; + + float out_r = mat.m[0][0]*in_r + mat.m[1][0]*in_g + mat.m[2][0]*in_b; + float out_g = mat.m[0][1]*in_r + mat.m[1][1]*in_g + mat.m[2][1]*in_b; + float out_b = mat.m[0][2]*in_r + mat.m[1][2]*in_g + mat.m[2][2]*in_b; + + *dest++ = clamp_float(out_r); + *dest++ = clamp_float(out_g); + *dest++ = clamp_float(out_b); + } +} + +static struct qcms_modular_transform* qcms_modular_transform_alloc() { + return calloc(1, sizeof(struct qcms_modular_transform)); +} + +static void qcms_modular_transform_release(struct qcms_modular_transform *transform) +{ + struct qcms_modular_transform *next_transform; + while (transform != NULL) { + next_transform = transform->next_transform; + // clut may use a single block of memory. + // Perhaps we should remove this to simply the code. + if (transform->input_clut_table_r + transform->input_clut_table_length == transform->input_clut_table_g && transform->input_clut_table_g + transform->input_clut_table_length == transform->input_clut_table_b) { + if (transform->input_clut_table_r) free(transform->input_clut_table_r); + } else { + if (transform->input_clut_table_r) free(transform->input_clut_table_r); + if (transform->input_clut_table_g) free(transform->input_clut_table_g); + if (transform->input_clut_table_b) free(transform->input_clut_table_b); + } + if (transform->r_clut + 1 == transform->g_clut && transform->g_clut + 1 == transform->b_clut) { + if (transform->r_clut) free(transform->r_clut); + } else { + if (transform->r_clut) free(transform->r_clut); + if (transform->g_clut) free(transform->g_clut); + if (transform->b_clut) free(transform->b_clut); + } + if (transform->output_clut_table_r + transform->output_clut_table_length == transform->output_clut_table_g && transform->output_clut_table_g+ transform->output_clut_table_length == transform->output_clut_table_b) { + if (transform->output_clut_table_r) free(transform->output_clut_table_r); + } else { + if (transform->output_clut_table_r) free(transform->output_clut_table_r); + if (transform->output_clut_table_g) free(transform->output_clut_table_g); + if (transform->output_clut_table_b) free(transform->output_clut_table_b); + } + if (transform->output_gamma_lut_r) free(transform->output_gamma_lut_r); + if (transform->output_gamma_lut_g) free(transform->output_gamma_lut_g); + if (transform->output_gamma_lut_b) free(transform->output_gamma_lut_b); + free(transform); + transform = next_transform; + } +} + +/* Set transform to be the next element in the linked list. */ +static void append_transform(struct qcms_modular_transform *transform, struct qcms_modular_transform ***next_transform) +{ + **next_transform = transform; + while (transform) { + *next_transform = &(transform->next_transform); + transform = transform->next_transform; + } +} + +/* reverse the transformation list (used by mBA) */ +static struct qcms_modular_transform* reverse_transform(struct qcms_modular_transform *transform) +{ + struct qcms_modular_transform *prev_transform = NULL; + while (transform != NULL) { + struct qcms_modular_transform *next_transform = transform->next_transform; + transform->next_transform = prev_transform; + prev_transform = transform; + transform = next_transform; + } + + return prev_transform; +} + +#define EMPTY_TRANSFORM_LIST NULL +static struct qcms_modular_transform* qcms_modular_transform_create_mAB(struct lutmABType *lut) +{ + struct qcms_modular_transform *first_transform = NULL; + struct qcms_modular_transform **next_transform = &first_transform; + struct qcms_modular_transform *transform = NULL; + + if (lut->a_curves[0] != NULL) { + size_t clut_length; + float *clut; + + // If the A curve is present this also implies the + // presence of a CLUT. + if (!lut->clut_table) + goto fail; + + // Prepare A curve. + transform = qcms_modular_transform_alloc(); + if (!transform) + goto fail; + append_transform(transform, &next_transform); + transform->input_clut_table_r = build_input_gamma_table(lut->a_curves[0]); + transform->input_clut_table_g = build_input_gamma_table(lut->a_curves[1]); + transform->input_clut_table_b = build_input_gamma_table(lut->a_curves[2]); + transform->transform_module_fn = qcms_transform_module_gamma_table; + if (lut->num_grid_points[0] != lut->num_grid_points[1] || + lut->num_grid_points[1] != lut->num_grid_points[2] ) { + //XXX: We don't currently support clut that are not squared! + goto fail; + } + + // Prepare CLUT + transform = qcms_modular_transform_alloc(); + if (!transform) + goto fail; + append_transform(transform, &next_transform); + clut_length = sizeof(float)*pow(lut->num_grid_points[0], 3)*3; + clut = malloc(clut_length); + if (!clut) + goto fail; + memcpy(clut, lut->clut_table, clut_length); + transform->r_clut = clut + 0; + transform->g_clut = clut + 1; + transform->b_clut = clut + 2; + transform->grid_size = lut->num_grid_points[0]; + transform->transform_module_fn = qcms_transform_module_clut_only; + } + if (lut->m_curves[0] != NULL) { + // M curve imples the presence of a Matrix + + // Prepare M curve + transform = qcms_modular_transform_alloc(); + if (!transform) + goto fail; + append_transform(transform, &next_transform); + transform->input_clut_table_r = build_input_gamma_table(lut->m_curves[0]); + transform->input_clut_table_g = build_input_gamma_table(lut->m_curves[1]); + transform->input_clut_table_b = build_input_gamma_table(lut->m_curves[2]); + transform->transform_module_fn = qcms_transform_module_gamma_table; + + // Prepare Matrix + transform = qcms_modular_transform_alloc(); + if (!transform) + goto fail; + append_transform(transform, &next_transform); + transform->matrix = build_mAB_matrix(lut); + if (transform->matrix.invalid) + goto fail; + transform->tx = s15Fixed16Number_to_float(lut->e03); + transform->ty = s15Fixed16Number_to_float(lut->e13); + transform->tz = s15Fixed16Number_to_float(lut->e23); + transform->transform_module_fn = qcms_transform_module_matrix_translate; + } + if (lut->b_curves[0] != NULL) { + // Prepare B curve + transform = qcms_modular_transform_alloc(); + if (!transform) + goto fail; + append_transform(transform, &next_transform); + transform->input_clut_table_r = build_input_gamma_table(lut->b_curves[0]); + transform->input_clut_table_g = build_input_gamma_table(lut->b_curves[1]); + transform->input_clut_table_b = build_input_gamma_table(lut->b_curves[2]); + transform->transform_module_fn = qcms_transform_module_gamma_table; + } else { + // B curve is mandatory + goto fail; + } + + if (lut->reversed) { + // mBA are identical to mAB except that the transformation order + // is reversed + first_transform = reverse_transform(first_transform); + } + + return first_transform; +fail: + qcms_modular_transform_release(first_transform); + return NULL; +} + +static struct qcms_modular_transform* qcms_modular_transform_create_lut(struct lutType *lut) +{ + struct qcms_modular_transform *first_transform = NULL; + struct qcms_modular_transform **next_transform = &first_transform; + struct qcms_modular_transform *transform = NULL; + + size_t in_curve_len, clut_length, out_curve_len; + float *in_curves, *clut, *out_curves; + + // Prepare Matrix + transform = qcms_modular_transform_alloc(); + if (!transform) + goto fail; + append_transform(transform, &next_transform); + transform->matrix = build_lut_matrix(lut); + if (transform->matrix.invalid) + goto fail; + transform->transform_module_fn = qcms_transform_module_matrix; + + // Prepare input curves + transform = qcms_modular_transform_alloc(); + if (!transform) + goto fail; + append_transform(transform, &next_transform); + in_curve_len = sizeof(float)*lut->num_input_table_entries * 3; + in_curves = malloc(in_curve_len); + if (!in_curves) + goto fail; + memcpy(in_curves, lut->input_table, in_curve_len); + transform->input_clut_table_r = in_curves + lut->num_input_table_entries * 0; + transform->input_clut_table_g = in_curves + lut->num_input_table_entries * 1; + transform->input_clut_table_b = in_curves + lut->num_input_table_entries * 2; + transform->input_clut_table_length = lut->num_input_table_entries; + + // Prepare table + clut_length = sizeof(float)*pow(lut->num_clut_grid_points, 3)*3; + clut = malloc(clut_length); + if (!clut) + goto fail; + memcpy(clut, lut->clut_table, clut_length); + transform->r_clut = clut + 0; + transform->g_clut = clut + 1; + transform->b_clut = clut + 2; + transform->grid_size = lut->num_clut_grid_points; + + // Prepare output curves + out_curve_len = sizeof(float) * lut->num_output_table_entries * 3; + out_curves = malloc(out_curve_len); + if (!out_curves) + goto fail; + memcpy(out_curves, lut->output_table, out_curve_len); + transform->output_clut_table_r = out_curves + lut->num_output_table_entries * 0; + transform->output_clut_table_g = out_curves + lut->num_output_table_entries * 1; + transform->output_clut_table_b = out_curves + lut->num_output_table_entries * 2; + transform->output_clut_table_length = lut->num_output_table_entries; + transform->transform_module_fn = qcms_transform_module_clut; + + return first_transform; +fail: + qcms_modular_transform_release(first_transform); + return NULL; +} + +struct qcms_modular_transform* qcms_modular_transform_create_input(qcms_profile *in) +{ + struct qcms_modular_transform *first_transform = NULL; + struct qcms_modular_transform **next_transform = &first_transform; + + if (in->A2B0) { + struct qcms_modular_transform *lut_transform; + lut_transform = qcms_modular_transform_create_lut(in->A2B0); + if (!lut_transform) + goto fail; + append_transform(lut_transform, &next_transform); + } else if (in->mAB && in->mAB->num_in_channels == 3 && in->mAB->num_out_channels == 3) { + struct qcms_modular_transform *mAB_transform; + mAB_transform = qcms_modular_transform_create_mAB(in->mAB); + if (!mAB_transform) + goto fail; + append_transform(mAB_transform, &next_transform); + + } else { + struct qcms_modular_transform *transform; + + transform = qcms_modular_transform_alloc(); + if (!transform) + goto fail; + append_transform(transform, &next_transform); + transform->input_clut_table_r = build_input_gamma_table(in->redTRC); + transform->input_clut_table_g = build_input_gamma_table(in->greenTRC); + transform->input_clut_table_b = build_input_gamma_table(in->blueTRC); + transform->transform_module_fn = qcms_transform_module_gamma_table; + if (!transform->input_clut_table_r || !transform->input_clut_table_g || + !transform->input_clut_table_b) { + goto fail; + } + + transform = qcms_modular_transform_alloc(); + if (!transform) + goto fail; + append_transform(transform, &next_transform); + transform->matrix.m[0][0] = 1/1.999969482421875f; + transform->matrix.m[0][1] = 0.f; + transform->matrix.m[0][2] = 0.f; + transform->matrix.m[1][0] = 0.f; + transform->matrix.m[1][1] = 1/1.999969482421875f; + transform->matrix.m[1][2] = 0.f; + transform->matrix.m[2][0] = 0.f; + transform->matrix.m[2][1] = 0.f; + transform->matrix.m[2][2] = 1/1.999969482421875f; + transform->matrix.invalid = false; + transform->transform_module_fn = qcms_transform_module_matrix; + + transform = qcms_modular_transform_alloc(); + if (!transform) + goto fail; + append_transform(transform, &next_transform); + transform->matrix = build_colorant_matrix(in); + transform->transform_module_fn = qcms_transform_module_matrix; + } + + return first_transform; +fail: + qcms_modular_transform_release(first_transform); + return EMPTY_TRANSFORM_LIST; +} +static struct qcms_modular_transform* qcms_modular_transform_create_output(qcms_profile *out) +{ + struct qcms_modular_transform *first_transform = NULL; + struct qcms_modular_transform **next_transform = &first_transform; + + if (out->B2A0) { + struct qcms_modular_transform *lut_transform; + lut_transform = qcms_modular_transform_create_lut(out->B2A0); + if (!lut_transform) + goto fail; + append_transform(lut_transform, &next_transform); + } else if (out->mBA && out->mBA->num_in_channels == 3 && out->mBA->num_out_channels == 3) { + struct qcms_modular_transform *lut_transform; + lut_transform = qcms_modular_transform_create_mAB(out->mBA); + if (!lut_transform) + goto fail; + append_transform(lut_transform, &next_transform); + } else if (out->redTRC && out->greenTRC && out->blueTRC) { + struct qcms_modular_transform *transform; + + transform = qcms_modular_transform_alloc(); + if (!transform) + goto fail; + append_transform(transform, &next_transform); + transform->matrix = matrix_invert(build_colorant_matrix(out)); + transform->transform_module_fn = qcms_transform_module_matrix; + + transform = qcms_modular_transform_alloc(); + if (!transform) + goto fail; + append_transform(transform, &next_transform); + transform->matrix.m[0][0] = 1.999969482421875f; + transform->matrix.m[0][1] = 0.f; + transform->matrix.m[0][2] = 0.f; + transform->matrix.m[1][0] = 0.f; + transform->matrix.m[1][1] = 1.999969482421875f; + transform->matrix.m[1][2] = 0.f; + transform->matrix.m[2][0] = 0.f; + transform->matrix.m[2][1] = 0.f; + transform->matrix.m[2][2] = 1.999969482421875f; + transform->matrix.invalid = false; + transform->transform_module_fn = qcms_transform_module_matrix; + + transform = qcms_modular_transform_alloc(); + if (!transform) + goto fail; + append_transform(transform, &next_transform); + build_output_lut(out->redTRC, &transform->output_gamma_lut_r, + &transform->output_gamma_lut_r_length); + build_output_lut(out->greenTRC, &transform->output_gamma_lut_g, + &transform->output_gamma_lut_g_length); + build_output_lut(out->blueTRC, &transform->output_gamma_lut_b, + &transform->output_gamma_lut_b_length); + transform->transform_module_fn = qcms_transform_module_gamma_lut; + + if (!transform->output_gamma_lut_r || !transform->output_gamma_lut_g || + !transform->output_gamma_lut_b) { + goto fail; + } + } else { + assert(0 && "Unsupported output profile workflow."); + return NULL; + } + + return first_transform; +fail: + qcms_modular_transform_release(first_transform); + return EMPTY_TRANSFORM_LIST; +} + +/* Not Completed +// Simply the transformation chain an equivilent transformation chain +static struct qcms_modular_transform* qcms_modular_transform_reduce(struct qcms_modular_transform *transform) +{ + struct qcms_modular_transform *first_transform = NULL; + struct qcms_modular_transform *curr_trans = transform; + struct qcms_modular_transform *prev_trans = NULL; + while (curr_trans) { + struct qcms_modular_transform *next_trans = curr_trans->next_transform; + if (curr_trans->transform_module_fn == qcms_transform_module_matrix) { + if (next_trans && next_trans->transform_module_fn == qcms_transform_module_matrix) { + curr_trans->matrix = matrix_multiply(curr_trans->matrix, next_trans->matrix); + goto remove_next; + } + } + if (curr_trans->transform_module_fn == qcms_transform_module_gamma_table) { + bool isLinear = true; + uint16_t i; + for (i = 0; isLinear && i < 256; i++) { + isLinear &= (int)(curr_trans->input_clut_table_r[i] * 255) == i; + isLinear &= (int)(curr_trans->input_clut_table_g[i] * 255) == i; + isLinear &= (int)(curr_trans->input_clut_table_b[i] * 255) == i; + } + goto remove_current; + } + +next_transform: + if (!next_trans) break; + prev_trans = curr_trans; + curr_trans = next_trans; + continue; +remove_current: + if (curr_trans == transform) { + //Update head + transform = next_trans; + } else { + prev_trans->next_transform = next_trans; + } + curr_trans->next_transform = NULL; + qcms_modular_transform_release(curr_trans); + //return transform; + return qcms_modular_transform_reduce(transform); +remove_next: + curr_trans->next_transform = next_trans->next_transform; + next_trans->next_transform = NULL; + qcms_modular_transform_release(next_trans); + continue; + } + return transform; +} +*/ + +static struct qcms_modular_transform* qcms_modular_transform_create(qcms_profile *in, qcms_profile *out) +{ + struct qcms_modular_transform *first_transform = NULL; + struct qcms_modular_transform **next_transform = &first_transform; + + if (in->color_space == RGB_SIGNATURE) { + struct qcms_modular_transform* rgb_to_xyz; + rgb_to_xyz = qcms_modular_transform_create_input(in); + if (!rgb_to_xyz) + goto fail; + append_transform(rgb_to_xyz, &next_transform); + } else { + assert(0 && "input color space not supported"); + goto fail; + } + + if (in->pcs == LAB_SIGNATURE && out->pcs == XYZ_SIGNATURE) { + struct qcms_modular_transform* lab_to_xyz; + lab_to_xyz = qcms_modular_transform_alloc(); + if (!lab_to_xyz) + goto fail; + append_transform(lab_to_xyz, &next_transform); + lab_to_xyz->transform_module_fn = qcms_transform_module_LAB_to_XYZ; + } + + //if (in->chromaticAdaption.invalid == false) { + // struct qcms_modular_transform* chromaticAdaption; + // chromaticAdaption = qcms_modular_transform_alloc(); + // if (!chromaticAdaption) + // goto fail; + // append_transform(chromaticAdaption, &next_transform); + // chromaticAdaption->matrix = matrix_invert(in->chromaticAdaption); + // chromaticAdaption->transform_module_fn = qcms_transform_module_matrix; + //} + + if (in->pcs == XYZ_SIGNATURE && out->pcs == LAB_SIGNATURE) { + struct qcms_modular_transform* xyz_to_lab; + xyz_to_lab = qcms_modular_transform_alloc(); + if (!xyz_to_lab) + goto fail; + append_transform(xyz_to_lab, &next_transform); + xyz_to_lab->transform_module_fn = qcms_transform_module_XYZ_to_LAB; + } + + if (out->color_space == RGB_SIGNATURE) { + struct qcms_modular_transform* xyz_to_rgb; + xyz_to_rgb = qcms_modular_transform_create_output(out); + if (!xyz_to_rgb) + goto fail; + append_transform(xyz_to_rgb, &next_transform); + } else { + assert(0 && "output color space not supported"); + goto fail; + } + // Not Completed + //return qcms_modular_transform_reduce(first_transform); + return first_transform; +fail: + qcms_modular_transform_release(first_transform); + return EMPTY_TRANSFORM_LIST; +} + +static float* qcms_modular_transform_data(struct qcms_modular_transform *transform, float *src, float *dest, size_t len) +{ + while (transform != NULL) { + // Keep swaping src/dest when performing a transform to use less memory. + float *new_src = dest; + const void *transform_fn = transform->transform_module_fn; + if (transform_fn != qcms_transform_module_gamma_table && + transform_fn != qcms_transform_module_gamma_lut && + transform_fn != qcms_transform_module_clut && + transform_fn != qcms_transform_module_clut_only && + transform_fn != qcms_transform_module_matrix && + transform_fn != qcms_transform_module_matrix_translate && + transform_fn != qcms_transform_module_LAB_to_XYZ && + transform_fn != qcms_transform_module_XYZ_to_LAB) { + assert(0 && "Unsupported transform module"); + return NULL; + } + transform->transform_module_fn(transform,src,dest,len); + dest = src; + src = new_src; + transform = transform->next_transform; + } + // The results end up in the src buffer because of the switching + return src; +} + +float* qcms_chain_transform(qcms_profile *in, qcms_profile *out, float *src, float *dest, size_t lutSize) +{ + struct qcms_modular_transform *transform_list = qcms_modular_transform_create(in, out); + if (transform_list != NULL) { + float *lut = qcms_modular_transform_data(transform_list, src, dest, lutSize/3); + qcms_modular_transform_release(transform_list); + return lut; + } + return NULL; +} @@ -0,0 +1,30 @@ +/* vim: set ts=8 sw=8 noexpandtab: */ +// qcms +// Copyright (C) 2009 Mozilla Foundation +// Copyright (C) 1998-2007 Marti Maria +// +// 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: +// +// The above copyright notice and this permission notice shall be included in +// all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO +// THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE +// LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION +// OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION +// WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + +#ifndef _QCMS_CHAIN_H +#define _QCMS_CHAIN_H + +// Generates and returns a 3D LUT with lutSize^3 samples using the provided src/dest. +float* qcms_chain_transform(qcms_profile *in, qcms_profile *out, float *src, float *dest, size_t lutSize); + +#endif diff --git a/matrix.c b/matrix.c new file mode 100644 index 0000000..0ce5bd6 --- /dev/null +++ b/matrix.c @@ -0,0 +1,136 @@ +/* vim: set ts=8 sw=8 noexpandtab: */ +// qcms +// Copyright (C) 2009 Mozilla Foundation +// Copyright (C) 1998-2007 Marti Maria +// +// 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: +// +// The above copyright notice and this permission notice shall be included in +// all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO +// THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE +// LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION +// OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION +// WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + +#include <stdlib.h> +#include "qcmsint.h" +#include "matrix.h" + +struct vector matrix_eval(struct matrix mat, struct vector v) +{ + struct vector result; + result.v[0] = mat.m[0][0]*v.v[0] + mat.m[0][1]*v.v[1] + mat.m[0][2]*v.v[2]; + result.v[1] = mat.m[1][0]*v.v[0] + mat.m[1][1]*v.v[1] + mat.m[1][2]*v.v[2]; + result.v[2] = mat.m[2][0]*v.v[0] + mat.m[2][1]*v.v[1] + mat.m[2][2]*v.v[2]; + return result; +} + +//XXX: should probably pass by reference and we could +//probably reuse this computation in matrix_invert +float matrix_det(struct matrix mat) +{ + float det; + det = mat.m[0][0]*mat.m[1][1]*mat.m[2][2] + + mat.m[0][1]*mat.m[1][2]*mat.m[2][0] + + mat.m[0][2]*mat.m[1][0]*mat.m[2][1] - + mat.m[0][0]*mat.m[1][2]*mat.m[2][1] - + mat.m[0][1]*mat.m[1][0]*mat.m[2][2] - + mat.m[0][2]*mat.m[1][1]*mat.m[2][0]; + return det; +} + +/* from pixman and cairo and Mathematics for Game Programmers */ +/* lcms uses gauss-jordan elimination with partial pivoting which is + * less efficient and not as numerically stable. See Mathematics for + * Game Programmers. */ +struct matrix matrix_invert(struct matrix mat) +{ + struct matrix dest_mat; + int i,j; + static int a[3] = { 2, 2, 1 }; + static int b[3] = { 1, 0, 0 }; + + /* inv (A) = 1/det (A) * adj (A) */ + float det = matrix_det(mat); + + if (det == 0) { + dest_mat.invalid = true; + } else { + dest_mat.invalid = false; + } + + det = 1/det; + + for (j = 0; j < 3; j++) { + for (i = 0; i < 3; i++) { + double p; + int ai = a[i]; + int aj = a[j]; + int bi = b[i]; + int bj = b[j]; + + p = mat.m[ai][aj] * mat.m[bi][bj] - + mat.m[ai][bj] * mat.m[bi][aj]; + if (((i + j) & 1) != 0) + p = -p; + + dest_mat.m[j][i] = det * p; + } + } + return dest_mat; +} + +struct matrix matrix_identity(void) +{ + struct matrix i; + i.m[0][0] = 1; + i.m[0][1] = 0; + i.m[0][2] = 0; + i.m[1][0] = 0; + i.m[1][1] = 1; + i.m[1][2] = 0; + i.m[2][0] = 0; + i.m[2][1] = 0; + i.m[2][2] = 1; + i.invalid = false; + return i; +} + +struct matrix matrix_invalid(void) +{ + struct matrix inv = matrix_identity(); + inv.invalid = true; + return inv; +} + + +/* from pixman */ +/* MAT3per... */ +struct matrix matrix_multiply(struct matrix a, struct matrix b) +{ + struct matrix result; + int dx, dy; + int o; + for (dy = 0; dy < 3; dy++) { + for (dx = 0; dx < 3; dx++) { + double v = 0; + for (o = 0; o < 3; o++) { + v += a.m[dy][o] * b.m[o][dx]; + } + result.m[dy][dx] = v; + } + } + result.invalid = a.invalid || b.invalid; + return result; +} + + diff --git a/matrix.h b/matrix.h new file mode 100644 index 0000000..5011988 --- /dev/null +++ b/matrix.h @@ -0,0 +1,39 @@ +/* vim: set ts=8 sw=8 noexpandtab: */ +// qcms +// Copyright (C) 2009 Mozilla Foundation +// Copyright (C) 1998-2007 Marti Maria +// +// 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: +// +// The above copyright notice and this permission notice shall be included in +// all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO +// THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE +// LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION +// OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION +// WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + +#ifndef _QCMS_MATRIX_H +#define _QCMS_MATRIX_H + +struct vector { + float v[3]; +}; + +struct vector matrix_eval(struct matrix mat, struct vector v); +float matrix_det(struct matrix mat); +struct matrix matrix_identity(void); +struct matrix matrix_multiply(struct matrix a, struct matrix b); +struct matrix matrix_invert(struct matrix mat); + +struct matrix matrix_invalid(void); + +#endif diff --git a/transform_util.c b/transform_util.c new file mode 100644 index 0000000..c29a98a --- /dev/null +++ b/transform_util.c @@ -0,0 +1,534 @@ +#include <math.h> +#include <assert.h> +#include <string.h> //memcpy +#include "qcmsint.h" +#include "transform_util.h" +#include "matrix.h" + +#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, int length) +{ + int upper, lower; + value = value * (length - 1); // scale to length of the array + upper = ceil(value); + lower = floor(value); + //XXX: can we be more performant here? + value = table[upper]*(1. - (upper - value)) + table[lower]*(upper - value); + /* scale the value */ + return value * (1./65535.); +} + +/* 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, int length) +{ + /* Start scaling input_value to the length of the array: 65535*(length-1). + * We'll divide out the 65535 next */ + uint32_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) */ + /* interp is the distance from upper to value scaled to 0..65535 */ + uint32_t interp = value % 65535; + + value = (table[upper]*(interp) + table[lower]*(65535 - interp))/65535; // 0..65535*65535 + + return value; +} + +/* same as above but takes an input_value from 0..PRECACHE_OUTPUT_MAX + * and returns a uint8_t value representing a range from 0..1 */ +static +uint8_t lut_interp_linear_precache_output(uint32_t input_value, uint16_t *table, int length) +{ + /* Start scaling input_value to the length of the array: PRECACHE_OUTPUT_MAX*(length-1). + * We'll divide out the PRECACHE_OUTPUT_MAX next */ + uint32_t value = (input_value * (length - 1)); + + /* equivalent to ceil(value/PRECACHE_OUTPUT_MAX) */ + uint32_t upper = (value + PRECACHE_OUTPUT_MAX-1) / PRECACHE_OUTPUT_MAX; + /* equivalent to floor(value/PRECACHE_OUTPUT_MAX) */ + uint32_t lower = value / PRECACHE_OUTPUT_MAX; + /* interp is the distance from upper to value scaled to 0..PRECACHE_OUTPUT_MAX */ + uint32_t interp = value % PRECACHE_OUTPUT_MAX; + + /* the table values range from 0..65535 */ + value = (table[upper]*(interp) + table[lower]*(PRECACHE_OUTPUT_MAX - interp)); // 0..(65535*PRECACHE_OUTPUT_MAX) + + /* round and scale */ + value += (PRECACHE_OUTPUT_MAX*65535/255)/2; + value /= (PRECACHE_OUTPUT_MAX*65535/255); // scale to 0..255 + return value; +} + +/* value must be a value between 0 and 1 */ +//XXX: is the above a good restriction to have? +float lut_interp_linear_float(float value, float *table, int length) +{ + int upper, lower; + value = value * (length - 1); + upper = ceil(value); + lower = floor(value); + //XXX: can we be more performant here? + value = table[upper]*(1. - (upper - value)) + table[lower]*(upper - value); + /* scale the value */ + return value; +} + +#if 0 +/* if we use a different representation i.e. one that goes from 0 to 0x1000 we can be more efficient + * because we can avoid the divisions and use a shifting instead */ +/* 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, int length) +{ + uint32_t value = (input_value * (length - 1)); + 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) +{ + 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]; + if (count == 0) { + a = 1; + b = 0; + c = 0; + e = 0; + f = 0; + interval = -INFINITY; + } else if(count == 1) { + a = parameter[1]; + b = parameter[2]; + c = 0; + e = 0; + f = 0; + interval = -1 * parameter[2] / parameter[1]; + } else if(count == 2) { + a = parameter[1]; + b = parameter[2]; + c = 0; + e = parameter[3]; + f = parameter[3]; + interval = -1 * parameter[2] / parameter[1]; + } else if(count == 3) { + a = parameter[1]; + b = parameter[2]; + c = parameter[3]; + e = -c; + f = 0; + interval = parameter[4]; + } else if(count == 4) { + a = parameter[1]; + b = parameter[2]; + c = parameter[3]; + e = parameter[5] - c; + f = parameter[6]; + interval = parameter[4]; + } 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) { + // 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; + } else { + gamma_table[X] = 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; +} + +unsigned char clamp_u8(float v) +{ + if (v > 255.) + return 255; + else if (v < 0) + return 0; + else + return floor(v+.5); +} + +float u8Fixed8Number_to_float(uint16_t x) +{ + // 0x0000 = 0. + // 0x0100 = 1. + // 0xffff = 255 + 255/256 + return x/256.; +} + +float *build_input_gamma_table(struct curveType *TRC) +{ + float *gamma_table; + + 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])); + } else { + compute_curve_gamma_table_type2(gamma_table, TRC->data, TRC->count); + } + } + } + return gamma_table; +} + +struct matrix build_colorant_matrix(qcms_profile *p) +{ + struct matrix result; + result.m[0][0] = s15Fixed16Number_to_float(p->redColorant.X); + result.m[0][1] = s15Fixed16Number_to_float(p->greenColorant.X); + result.m[0][2] = s15Fixed16Number_to_float(p->blueColorant.X); + result.m[1][0] = s15Fixed16Number_to_float(p->redColorant.Y); + result.m[1][1] = s15Fixed16Number_to_float(p->greenColorant.Y); + result.m[1][2] = s15Fixed16Number_to_float(p->blueColorant.Y); + result.m[2][0] = s15Fixed16Number_to_float(p->redColorant.Z); + result.m[2][1] = s15Fixed16Number_to_float(p->greenColorant.Z); + result.m[2][2] = s15Fixed16Number_to_float(p->blueColorant.Z); + result.invalid = false; + return result; +} + +/* The following code is copied nearly directly from lcms. + * I think it could be much better. For example, Argyll seems to have better code in + * icmTable_lookup_bwd and icmTable_setup_bwd. However, for now this is a quick way + * to a working solution and allows for easy comparing with lcms. */ +uint16_fract_t lut_inverse_interp16(uint16_t Value, uint16_t LutTable[], int length) +{ + int l = 1; + int r = 0x10000; + int x = 0, res; // 'int' Give spacing for negative values + int NumZeroes, NumPoles; + int cell0, cell1; + double val2; + double y0, y1, x0, x1; + double a, b, f; + + // July/27 2001 - Expanded to handle degenerated curves with an arbitrary + // number of elements containing 0 at the begining of the table (Zeroes) + // and another arbitrary number of poles (FFFFh) at the end. + // First the zero and pole extents are computed, then value is compared. + + NumZeroes = 0; + while (LutTable[NumZeroes] == 0 && NumZeroes < length-1) + NumZeroes++; + + // There are no zeros at the beginning and we are trying to find a zero, so + // return anything. It seems zero would be the less destructive choice + /* I'm not sure that this makes sense, but oh well... */ + if (NumZeroes == 0 && Value == 0) + return 0; + + NumPoles = 0; + while (LutTable[length-1- NumPoles] == 0xFFFF && NumPoles < length-1) + NumPoles++; + + // Does the curve belong to this case? + if (NumZeroes > 1 || NumPoles > 1) + { + int a, b; + + // Identify if value fall downto 0 or FFFF zone + if (Value == 0) return 0; + // if (Value == 0xFFFF) return 0xFFFF; + + // else restrict to valid zone + + a = ((NumZeroes-1) * 0xFFFF) / (length-1); + b = ((length-1 - NumPoles) * 0xFFFF) / (length-1); + + l = a - 1; + r = b + 1; + } + + + // Seems not a degenerated case... apply binary search + + while (r > l) { + + x = (l + r) / 2; + + res = (int) lut_interp_linear16((uint16_fract_t) (x-1), LutTable, length); + + if (res == Value) { + + // Found exact match. + + return (uint16_fract_t) (x - 1); + } + + if (res > Value) r = x - 1; + else l = x + 1; + } + + // Not found, should we interpolate? + + + // Get surrounding nodes + + val2 = (length-1) * ((double) (x - 1) / 65535.0); + + cell0 = (int) floor(val2); + cell1 = (int) ceil(val2); + + if (cell0 == cell1) return (uint16_fract_t) x; + + y0 = LutTable[cell0] ; + x0 = (65535.0 * cell0) / (length-1); + + y1 = LutTable[cell1] ; + x1 = (65535.0 * cell1) / (length-1); + + a = (y1 - y0) / (x1 - x0); + b = y0 - a * x0; + + if (fabs(a) < 0.01) return (uint16_fract_t) x; + + f = ((Value - b) / a); + + if (f < 0.0) return (uint16_fract_t) 0; + if (f >= 65535.0) return (uint16_fract_t) 0xFFFF; + + return (uint16_fract_t) floor(f + 0.5); + +} + +/* + The number of entries needed to invert a lookup table should not + necessarily be the same as the original number of entries. This is + especially true of lookup tables that have a small number of entries. + + For example: + Using a table like: + {0, 3104, 14263, 34802, 65535} + invert_lut will produce an inverse of: + {3, 34459, 47529, 56801, 65535} + which has an maximum error of about 9855 (pixel difference of ~38.346) + + For now, we punt the decision of output size to the caller. */ +static uint16_t *invert_lut(uint16_t *table, int length, int out_length) +{ + int i; + /* for now we invert the lut by creating a lut of size out_length + * and attempting to lookup a value for each entry using lut_inverse_interp16 */ + uint16_t *output = malloc(sizeof(uint16_t)*out_length); + if (!output) + return NULL; + + for (i = 0; i < out_length; i++) { + double x = ((double) i * 65535.) / (double) (out_length - 1); + uint16_fract_t input = floor(x + .5); + output[i] = lut_inverse_interp16(input, table, length); + } + return output; +} + +static void compute_precache_pow(uint8_t *output, float gamma) +{ + uint32_t v = 0; + for (v = 0; v <= 0xffff; v++) { + //XXX: don't do integer/float conversion... and round? + output[v] = 255. * pow(v/65535., gamma); + } +} + +void compute_precache_lut(uint8_t *output, uint16_t *table, int length) +{ + uint32_t v = 0; + for (v = 0; v <= 0xffff; v++) { + //XXX: don't do integer/float conversion... round? + output[v] = lut_interp_linear16(v, table, length) >> 8; + } +} + +void compute_precache_linear(uint8_t *output) +{ + uint32_t v = 0; + for (v = 0; v <= 0xffff; v++) { + //XXX: round? + output[v] = v >> 8; + } +} + +qcms_bool compute_precache(struct curveType *trc, uint8_t *output) +{ + + if (trc->type == PARAMETRIC_CURVE_TYPE) { + float gamma_table[256]; + uint16_t gamma_table_uint[256]; + uint16_t i; + uint16_t *inverted; + int inverted_size = 256; + + compute_curve_gamma_table_type_parametric(gamma_table, trc->parameter, trc->count); + for(i = 0; i < 256; i++) { + gamma_table_uint[i] = (uint16_t)(gamma_table[i] * 65535); + } + + //XXX: the choice of a minimum of 256 here is not backed by any theory, + // measurement or data, howeve r it is what lcms uses. + // the maximum number we would need is 65535 because that's the + // accuracy used for computing the pre cache table + if (inverted_size < 256) + inverted_size = 256; + + inverted = invert_lut(gamma_table_uint, 256, inverted_size); + if (!inverted) + return false; + compute_precache_lut(output, inverted, inverted_size); + free(inverted); + } else { + if (trc->count == 0) { + compute_precache_linear(output); + } else if (trc->count == 1) { + compute_precache_pow(output, 1./u8Fixed8Number_to_float(trc->data[0])); + } else { + uint16_t *inverted; + int inverted_size = trc->count; + //XXX: the choice of a minimum of 256 here is not backed by any theory, + // measurement or data, howeve r it is what lcms uses. + // the maximum number we would need is 65535 because that's the + // accuracy used for computing the pre cache table + if (inverted_size < 256) + inverted_size = 256; + + inverted = invert_lut(trc->data, trc->count, inverted_size); + if (!inverted) + return false; + compute_precache_lut(output, inverted, inverted_size); + free(inverted); + } + } + return true; +} + + +static uint16_t *build_linear_table(int length) +{ + int i; + uint16_t *output = malloc(sizeof(uint16_t)*length); + if (!output) + return NULL; + + for (i = 0; i < length; i++) { + double x = ((double) i * 65535.) / (double) (length - 1); + uint16_fract_t input = floor(x + .5); + output[i] = input; + } + return output; +} + +static uint16_t *build_pow_table(float gamma, int length) +{ + int i; + uint16_t *output = malloc(sizeof(uint16_t)*length); + if (!output) + return NULL; + + for (i = 0; i < length; i++) { + uint16_fract_t result; + double x = ((double) i) / (double) (length - 1); + x = pow(x, gamma); //XXX turn this conversion into a function + result = floor(x*65535. + .5); + output[i] = result; + } + return output; +} + +void build_output_lut(struct curveType *trc, + uint16_t **output_gamma_lut, size_t *output_gamma_lut_length) +{ + if (trc->type == PARAMETRIC_CURVE_TYPE) { + float gamma_table[256]; + uint16_t i; + uint16_t *output = malloc(sizeof(uint16_t)*256); + + if (!output) { + *output_gamma_lut = NULL; + return; + } + + compute_curve_gamma_table_type_parametric(gamma_table, trc->parameter, trc->count); + *output_gamma_lut_length = 256; + for(i = 0; i < 256; i++) { + output[i] = (uint16_t)(gamma_table[i] * 65535); + } + *output_gamma_lut = output; + } else { + if (trc->count == 0) { + *output_gamma_lut = build_linear_table(4096); + *output_gamma_lut_length = 4096; + } else if (trc->count == 1) { + float gamma = 1./u8Fixed8Number_to_float(trc->data[0]); + *output_gamma_lut = build_pow_table(gamma, 4096); + *output_gamma_lut_length = 4096; + } else { + //XXX: the choice of a minimum of 256 here is not backed by any theory, + // 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); + } + } + +} + diff --git a/transform_util.h b/transform_util.h new file mode 100644 index 0000000..8f358a8 --- /dev/null +++ b/transform_util.h @@ -0,0 +1,59 @@ +/* vim: set ts=8 sw=8 noexpandtab: */ +// qcms +// Copyright (C) 2009 Mozilla Foundation +// Copyright (C) 1998-2007 Marti Maria +// +// 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: +// +// The above copyright notice and this permission notice shall be included in +// all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO +// THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE +// LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION +// OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION +// WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + +#ifndef _QCMS_TRANSFORM_UTIL_H +#define _QCMS_TRANSFORM_UTIL_H + +#include <stdlib.h> + +#define CLU(table,x,y,z) table[(x*len + y*x_len + z*xy_len)*3] + +//XXX: could use a bettername +typedef uint16_t uint16_fract_t; + +float lut_interp_linear(double value, uint16_t *table, int length); +float lut_interp_linear_float(float value, float *table, int length); +uint16_t lut_interp_linear16(uint16_t input_value, uint16_t *table, int length); + + +static inline float lerp(float a, float b, float t) +{ + return a*(1.f-t) + b*t; +} + +unsigned char clamp_u8(float v); +float clamp_float(float a); + +float u8Fixed8Number_to_float(uint16_t x); + + +float *build_input_gamma_table(struct curveType *TRC); +struct matrix build_colorant_matrix(qcms_profile *p); +void build_output_lut(struct curveType *trc, + uint16_t **output_gamma_lut, size_t *output_gamma_lut_length); + +struct matrix matrix_invert(struct matrix mat); +qcms_bool compute_precache(struct curveType *trc, uint8_t *output); + + +#endif |