image-colorspace.c [plain text]
#include "image.h"
#include <math.h>
#ifndef M_PI
# define M_PI 3.14159265358979323846
#endif
#ifndef M_SQRT2
# define M_SQRT2 1.41421356237309504880
#endif
#ifndef M_SQRT1_2
# define M_SQRT1_2 0.70710678118654752440
#endif
typedef int cups_clut_t[3][256];
static int ImageHaveProfile = 0;
static int *ImageDensity;
static cups_clut_t *ImageMatrix;
static cups_cspace_t ImageColorSpace = CUPS_CSPACE_RGB;
static float cielab(float x, float xn);
static void rgb_to_xyz(ib_t *val);
static void rgb_to_lab(ib_t *val);
static void huerotate(float [3][3], float);
static void ident(float [3][3]);
static void mult(float [3][3], float [3][3], float [3][3]);
static void saturate(float [3][3], float);
static void xform(float [3][3], float, float, float, float *, float *, float *);
static void xrotate(float [3][3], float, float);
static void yrotate(float [3][3], float, float);
static void zrotate(float [3][3], float, float);
static void zshear(float [3][3], float, float);
void
ImageSetColorSpace(cups_cspace_t cs)
{
ImageColorSpace = cs;
if (cs >= CUPS_CSPACE_CIEXYZ)
ImageHaveProfile = 0;
}
void
ImageSetProfile(float d,
float g,
float matrix[3][3])
{
int i, j, k;
float m;
int *im;
if (ImageMatrix == NULL)
ImageMatrix = calloc(3, sizeof(cups_clut_t));
if (ImageMatrix == NULL)
return;
if (ImageDensity == NULL)
ImageDensity = calloc(256, sizeof(int));
if (ImageDensity == NULL)
return;
ImageHaveProfile = 1;
for (i = 0, im = ImageMatrix[0][0]; i < 3; i ++)
for (j = 0; j < 3; j ++)
for (k = 0, m = matrix[i][j]; k < 256; k ++)
*im++ = (int)(k * m + 0.5);
for (k = 0, im = ImageDensity; k < 256; k ++)
*im++ = 255.0 * d * pow((float)k / 255.0, g) + 0.5;
}
void
ImageWhiteToWhite(const ib_t *in,
ib_t *out,
int count)
{
if (ImageHaveProfile)
while (count > 0)
{
*out++ = 255 - ImageDensity[255 - *in++];
count --;
}
else if (in != out)
memcpy(out, in, count);
}
void
ImageWhiteToRGB(const ib_t *in,
ib_t *out,
int count)
{
if (ImageHaveProfile)
{
while (count > 0)
{
out[0] = 255 - ImageDensity[255 - *in++];
out[1] = out[0];
out[2] = out[0];
out += 3;
count --;
}
}
else
{
while (count > 0)
{
*out++ = *in;
*out++ = *in;
*out++ = *in++;
if (ImageColorSpace >= CUPS_CSPACE_CIELab)
rgb_to_lab(out - 3);
else if (ImageColorSpace == CUPS_CSPACE_CIEXYZ)
rgb_to_xyz(out - 3);
count --;
}
}
}
void
ImageWhiteToBlack(const ib_t *in,
ib_t *out,
int count)
{
if (ImageHaveProfile)
while (count > 0)
{
*out++ = ImageDensity[255 - *in++];
count --;
}
else
while (count > 0)
{
*out++ = 255 - *in++;
count --;
}
}
void
ImageWhiteToCMY(const ib_t *in,
ib_t *out,
int count)
{
if (ImageHaveProfile)
while (count > 0)
{
out[0] = ImageDensity[255 - *in++];
out[1] = out[0];
out[2] = out[0];
out += 3;
count --;
}
else
while (count > 0)
{
*out++ = 255 - *in;
*out++ = 255 - *in;
*out++ = 255 - *in++;
count --;
}
}
void
ImageWhiteToCMYK(const ib_t *in,
ib_t *out,
int count)
{
if (ImageHaveProfile)
while (count > 0)
{
*out++ = 0;
*out++ = 0;
*out++ = 0;
*out++ = ImageDensity[255 - *in++];
count --;
}
else
while (count > 0)
{
*out++ = 0;
*out++ = 0;
*out++ = 0;
*out++ = 255 - *in++;
count --;
}
}
void
ImageRGBToBlack(const ib_t *in,
ib_t *out,
int count)
{
if (ImageHaveProfile)
while (count > 0)
{
*out++ = ImageDensity[255 - (31 * in[0] + 61 * in[1] + 8 * in[2]) / 100];
in += 3;
count --;
}
else
while (count > 0)
{
*out++ = 255 - (31 * in[0] + 61 * in[1] + 8 * in[2]) / 100;
in += 3;
count --;
}
}
void
ImageRGBToCMY(const ib_t *in,
ib_t *out,
int count)
{
int c, m, y, k;
int cc, cm, cy;
if (ImageHaveProfile)
while (count > 0)
{
c = 255 - *in++;
m = 255 - *in++;
y = 255 - *in++;
k = min(c, min(m, y));
c -= k;
m -= k;
y -= k;
cc = ImageMatrix[0][0][c] +
ImageMatrix[0][1][m] +
ImageMatrix[0][2][y] + k;
cm = ImageMatrix[1][0][c] +
ImageMatrix[1][1][m] +
ImageMatrix[1][2][y] + k;
cy = ImageMatrix[2][0][c] +
ImageMatrix[2][1][m] +
ImageMatrix[2][2][y] + k;
if (cc < 0)
*out++ = 0;
else if (cc > 255)
*out++ = ImageDensity[255];
else
*out++ = ImageDensity[cc];
if (cm < 0)
*out++ = 0;
else if (cm > 255)
*out++ = ImageDensity[255];
else
*out++ = ImageDensity[cm];
if (cy < 0)
*out++ = 0;
else if (cy > 255)
*out++ = ImageDensity[255];
else
*out++ = ImageDensity[cy];
count --;
}
else
while (count > 0)
{
c = 255 - in[0];
m = 255 - in[1];
y = 255 - in[2];
k = min(c, min(m, y));
*out++ = (255 - in[1] / 4) * (c - k) / 255 + k;
*out++ = (255 - in[2] / 4) * (m - k) / 255 + k;
*out++ = (255 - in[0] / 4) * (y - k) / 255 + k;
in += 3;
count --;
}
}
void
ImageRGBToCMYK(const ib_t *in,
ib_t *out,
int count)
{
int c, m, y, k,
km;
int cc, cm, cy;
if (ImageHaveProfile)
while (count > 0)
{
c = 255 - *in++;
m = 255 - *in++;
y = 255 - *in++;
k = min(c, min(m, y));
if ((km = max(c, max(m, y))) > k)
k = k * k * k / (km * km);
c -= k;
m -= k;
y -= k;
cc = (ImageMatrix[0][0][c] +
ImageMatrix[0][1][m] +
ImageMatrix[0][2][y]);
cm = (ImageMatrix[1][0][c] +
ImageMatrix[1][1][m] +
ImageMatrix[1][2][y]);
cy = (ImageMatrix[2][0][c] +
ImageMatrix[2][1][m] +
ImageMatrix[2][2][y]);
if (cc < 0)
*out++ = 0;
else if (cc > 255)
*out++ = ImageDensity[255];
else
*out++ = ImageDensity[cc];
if (cm < 0)
*out++ = 0;
else if (cm > 255)
*out++ = ImageDensity[255];
else
*out++ = ImageDensity[cm];
if (cy < 0)
*out++ = 0;
else if (cy > 255)
*out++ = ImageDensity[255];
else
*out++ = ImageDensity[cy];
*out++ = ImageDensity[k];
count --;
}
else
while (count > 0)
{
c = 255 - *in++;
m = 255 - *in++;
y = 255 - *in++;
k = min(c, min(m, y));
if ((km = max(c, max(m, y))) > k)
k = k * k * k / (km * km);
c -= k;
m -= k;
y -= k;
*out++ = c;
*out++ = m;
*out++ = y;
*out++ = k;
count --;
}
}
void
ImageRGBToWhite(const ib_t *in,
ib_t *out,
int count)
{
if (ImageHaveProfile)
{
while (count > 0)
{
*out++ = 255 - ImageDensity[255 - (31 * in[0] + 61 * in[1] + 8 * in[2]) / 100];
in += 3;
count --;
}
}
else
{
while (count > 0)
{
*out++ = (31 * in[0] + 61 * in[1] + 8 * in[2]) / 100;
in += 3;
count --;
}
}
}
void
ImageRGBToRGB(const ib_t *in,
ib_t *out,
int count)
{
int c, m, y, k;
int cr, cg, cb;
if (ImageHaveProfile)
{
while (count > 0)
{
c = 255 - *in++;
m = 255 - *in++;
y = 255 - *in++;
k = min(c, min(m, y));
c -= k;
m -= k;
y -= k;
cr = ImageMatrix[0][0][c] +
ImageMatrix[0][1][m] +
ImageMatrix[0][2][y] + k;
cg = ImageMatrix[1][0][c] +
ImageMatrix[1][1][m] +
ImageMatrix[1][2][y] + k;
cb = ImageMatrix[2][0][c] +
ImageMatrix[2][1][m] +
ImageMatrix[2][2][y] + k;
if (cr < 0)
*out++ = 255;
else if (cr > 255)
*out++ = 255 - ImageDensity[255];
else
*out++ = 255 - ImageDensity[cr];
if (cg < 0)
*out++ = 255;
else if (cg > 255)
*out++ = 255 - ImageDensity[255];
else
*out++ = 255 - ImageDensity[cg];
if (cb < 0)
*out++ = 255;
else if (cb > 255)
*out++ = 255 - ImageDensity[255];
else
*out++ = 255 - ImageDensity[cb];
count --;
}
}
else
{
if (in != out)
memcpy(out, in, count * 3);
if (ImageColorSpace >= CUPS_CSPACE_CIEXYZ)
{
while (count > 0)
{
if (ImageColorSpace >= CUPS_CSPACE_CIELab)
rgb_to_lab(out);
else
rgb_to_xyz(out);
out += 3;
count --;
}
}
}
}
void
ImageCMYKToBlack(const ib_t *in,
ib_t *out,
int count)
{
int k;
if (ImageHaveProfile)
while (count > 0)
{
k = (31 * in[0] + 61 * in[1] + 8 * in[2]) / 100 + in[3];
if (k < 255)
*out++ = ImageDensity[k];
else
*out++ = ImageDensity[255];
in += 4;
count --;
}
else
while (count > 0)
{
k = (31 * in[0] + 61 * in[1] + 8 * in[2]) / 100 + in[3];
if (k < 255)
*out++ = k;
else
*out++ = 255;
in += 4;
count --;
}
}
void
ImageCMYKToCMY(const ib_t *in,
ib_t *out,
int count)
{
int c, m, y, k;
int cc, cm, cy;
if (ImageHaveProfile)
while (count > 0)
{
c = *in++;
m = *in++;
y = *in++;
k = *in++;
cc = ImageMatrix[0][0][c] +
ImageMatrix[0][1][m] +
ImageMatrix[0][2][y] + k;
cm = ImageMatrix[1][0][c] +
ImageMatrix[1][1][m] +
ImageMatrix[1][2][y] + k;
cy = ImageMatrix[2][0][c] +
ImageMatrix[2][1][m] +
ImageMatrix[2][2][y] + k;
if (cc < 0)
*out++ = 0;
else if (cc > 255)
*out++ = ImageDensity[255];
else
*out++ = ImageDensity[cc];
if (cm < 0)
*out++ = 0;
else if (cm > 255)
*out++ = ImageDensity[255];
else
*out++ = ImageDensity[cm];
if (cy < 0)
*out++ = 0;
else if (cy > 255)
*out++ = ImageDensity[255];
else
*out++ = ImageDensity[cy];
count --;
}
else
while (count > 0)
{
c = *in++;
m = *in++;
y = *in++;
k = *in++;
c += k;
m += k;
y += k;
if (c < 255)
*out++ = c;
else
*out++ = 255;
if (m < 255)
*out++ = y;
else
*out++ = 255;
if (y < 255)
*out++ = y;
else
*out++ = 255;
count --;
}
}
void
ImageCMYKToCMYK(const ib_t *in,
ib_t *out,
int count)
{
int c, m, y, k;
int cc, cm, cy;
if (ImageHaveProfile)
while (count > 0)
{
c = *in++;
m = *in++;
y = *in++;
k = *in++;
cc = (ImageMatrix[0][0][c] +
ImageMatrix[0][1][m] +
ImageMatrix[0][2][y]);
cm = (ImageMatrix[1][0][c] +
ImageMatrix[1][1][m] +
ImageMatrix[1][2][y]);
cy = (ImageMatrix[2][0][c] +
ImageMatrix[2][1][m] +
ImageMatrix[2][2][y]);
if (cc < 0)
*out++ = 0;
else if (cc > 255)
*out++ = ImageDensity[255];
else
*out++ = ImageDensity[cc];
if (cm < 0)
*out++ = 0;
else if (cm > 255)
*out++ = ImageDensity[255];
else
*out++ = ImageDensity[cm];
if (cy < 0)
*out++ = 0;
else if (cy > 255)
*out++ = ImageDensity[255];
else
*out++ = ImageDensity[cy];
*out++ = ImageDensity[k];
count --;
}
else if (in != out)
{
while (count > 0)
{
*out++ = *in++;
*out++ = *in++;
*out++ = *in++;
*out++ = *in++;
count --;
}
}
}
void
ImageCMYKToWhite(const ib_t *in,
ib_t *out,
int count)
{
int w;
if (ImageHaveProfile)
{
while (count > 0)
{
w = 255 - (31 * in[0] + 61 * in[1] + 8 * in[2]) / 100 - in[3];
if (w > 0)
*out++ = ImageDensity[w];
else
*out++ = ImageDensity[0];
in += 4;
count --;
}
}
else
{
while (count > 0)
{
w = 255 - (31 * in[0] + 61 * in[1] + 8 * in[2]) / 100 - in[3];
if (w > 0)
*out++ = w;
else
*out++ = 0;
in += 4;
count --;
}
}
}
void
ImageCMYKToRGB(const ib_t *in,
ib_t *out,
int count)
{
int c, m, y, k;
int cr, cg, cb;
if (ImageHaveProfile)
{
while (count > 0)
{
c = *in++;
m = *in++;
y = *in++;
k = *in++;
cr = ImageMatrix[0][0][c] +
ImageMatrix[0][1][m] +
ImageMatrix[0][2][y] + k;
cg = ImageMatrix[1][0][c] +
ImageMatrix[1][1][m] +
ImageMatrix[1][2][y] + k;
cb = ImageMatrix[2][0][c] +
ImageMatrix[2][1][m] +
ImageMatrix[2][2][y] + k;
if (cr < 0)
*out++ = 255;
else if (cr > 255)
*out++ = 255 - ImageDensity[255];
else
*out++ = 255 - ImageDensity[cr];
if (cg < 0)
*out++ = 255;
else if (cg > 255)
*out++ = 255 - ImageDensity[255];
else
*out++ = 255 - ImageDensity[cg];
if (cb < 0)
*out++ = 255;
else if (cb > 255)
*out++ = 255 - ImageDensity[255];
else
*out++ = 255 - ImageDensity[cb];
count --;
}
}
else
{
while (count > 0)
{
c = 255 - *in++;
m = 255 - *in++;
y = 255 - *in++;
k = *in++;
c -= k;
m -= k;
y -= k;
if (c > 0)
*out++ = c;
else
*out++ = 0;
if (m > 0)
*out++ = m;
else
*out++ = 0;
if (y > 0)
*out++ = y;
else
*out++ = 0;
if (ImageColorSpace >= CUPS_CSPACE_CIELab)
rgb_to_lab(out - 3);
else if (ImageColorSpace == CUPS_CSPACE_CIEXYZ)
rgb_to_xyz(out - 3);
count --;
}
}
}
void
ImageLut(ib_t *pixels,
int count,
const ib_t *lut)
{
while (count > 0)
{
*pixels = lut[*pixels];
pixels ++;
count --;
}
}
void
ImageRGBAdjust(ib_t *pixels,
int count,
int saturation,
int hue)
{
int i, j, k;
float mat[3][3];
static int last_sat = 100,
last_hue = 0;
static cups_clut_t *lut = NULL;
if (saturation != last_sat ||
hue != last_hue)
{
ident(mat);
saturate(mat, saturation * 0.01);
huerotate(mat, (float)hue);
if (lut == NULL)
lut = calloc(3, sizeof(cups_clut_t));
if (lut == NULL)
return;
for (i = 0; i < 3; i ++)
for (j = 0; j < 3; j ++)
for (k = 0; k < 256; k ++)
lut[i][j][k] = mat[i][j] * k + 0.5;
last_sat = saturation;
last_hue = hue;
}
while (count > 0)
{
i = lut[0][0][pixels[0]] +
lut[1][0][pixels[1]] +
lut[2][0][pixels[2]];
if (i < 0)
pixels[0] = 0;
else if (i > 255)
pixels[0] = 255;
else
pixels[0] = i;
i = lut[0][1][pixels[0]] +
lut[1][1][pixels[1]] +
lut[2][1][pixels[2]];
if (i < 0)
pixels[1] = 0;
else if (i > 255)
pixels[1] = 255;
else
pixels[1] = i;
i = lut[0][2][pixels[0]] +
lut[1][2][pixels[1]] +
lut[2][2][pixels[2]];
if (i < 0)
pixels[2] = 0;
else if (i > 255)
pixels[2] = 255;
else
pixels[2] = i;
count --;
pixels += 3;
}
}
#define D65_X (0.412453 + 0.357580 + 0.180423)
#define D65_Y (0.212671 + 0.715160 + 0.072169)
#define D65_Z (0.019334 + 0.119193 + 0.950227)
static float
cielab(float x,
float xn)
{
float x_xn;
x_xn = x / xn;
if (x_xn > 0.008856)
return (cbrt(x_xn));
else
return (7.787 * x_xn + 16.0 / 116.0);
}
static void
rgb_to_xyz(ib_t *val)
{
float r,
g,
b,
ciex,
ciey,
ciez;
r = pow(val[0] / 255.0, 0.58823529412);
g = pow(val[1] / 255.0, 0.58823529412);
b = pow(val[2] / 255.0, 0.58823529412);
ciex = 0.412453 * r + 0.357580 * g + 0.180423 * b;
ciey = 0.212671 * r + 0.715160 * g + 0.072169 * b;
ciez = 0.019334 * r + 0.119193 * g + 0.950227 * b;
if (ciex < 0.0)
val[0] = 0;
else if (ciex < 255.0)
val[0] = (int)ciex;
else
val[0] = 255;
if (ciey < 0.0)
val[1] = 0;
else if (ciey < 255.0)
val[1] = (int)ciey;
else
val[1] = 255;
if (ciez < 0.0)
val[2] = 0;
else if (ciez < 255.0)
val[2] = (int)ciez;
else
val[2] = 255;
}
static void
rgb_to_lab(ib_t *val)
{
float r,
g,
b,
ciex,
ciey,
ciez,
ciey_yn,
ciel,
ciea,
cieb;
r = pow(val[0] / 255.0, 0.58823529412);
g = pow(val[1] / 255.0, 0.58823529412);
b = pow(val[2] / 255.0, 0.58823529412);
ciex = 0.412453 * r + 0.357580 * g + 0.180423 * b;
ciey = 0.212671 * r + 0.715160 * g + 0.072169 * b;
ciez = 0.019334 * r + 0.119193 * g + 0.950227 * b;
ciey_yn = ciey / D65_Y;
if (ciey_yn > 0.008856)
ciel = 116 * cbrt(ciey_yn) - 16;
else
ciel = 903.3 * ciey_yn;
ciel = ciel;
ciea = 500 * (cielab(ciex, D65_X) - cielab(ciey, D65_Y));
cieb = 200 * (cielab(ciey, D65_Y) - cielab(ciez, D65_Z));
ciel *= 2.55;
ciea += 128;
cieb += 128;
if (ciel < 0.0)
val[0] = 0;
else if (ciel < 255.0)
val[0] = (int)ciel;
else
val[0] = 255;
if (ciea < 0.0)
val[1] = 128;
else if (ciea < 255.0)
val[1] = (int)ciea;
else
val[1] = 255;
if (cieb < 0.0)
val[2] = 128;
else if (cieb < 255.0)
val[2] = (int)cieb;
else
val[2] = 255;
}
static void
huerotate(float mat[3][3],
float rot)
{
float hmat[3][3];
float lx, ly, lz;
float xrs, xrc;
float yrs, yrc;
float zrs, zrc;
float zsx, zsy;
ident(hmat);
xrs = M_SQRT1_2;
xrc = M_SQRT1_2;
xrotate(hmat,xrs,xrc);
yrs = -1.0 / sqrt(3.0);
yrc = -M_SQRT2 * yrs;
yrotate(hmat,yrs,yrc);
xform(hmat, 0.3086, 0.6094, 0.0820, &lx, &ly, &lz);
zsx = lx / lz;
zsy = ly / lz;
zshear(hmat, zsx, zsy);
zrs = sin(rot * M_PI / 180.0);
zrc = cos(rot * M_PI / 180.0);
zrotate(hmat, zrs, zrc);
zshear(hmat, -zsx, -zsy);
yrotate(hmat, -yrs, yrc);
xrotate(hmat, -xrs, xrc);
mult(hmat, mat, mat);
}
static void
ident(float mat[3][3])
{
mat[0][0] = 1.0;
mat[0][1] = 0.0;
mat[0][2] = 0.0;
mat[1][0] = 0.0;
mat[1][1] = 1.0;
mat[1][2] = 0.0;
mat[2][0] = 0.0;
mat[2][1] = 0.0;
mat[2][2] = 1.0;
}
static void
mult(float a[3][3],
float b[3][3],
float c[3][3])
{
int x, y;
float temp[3][3];
for (y = 0; y < 3; y ++)
for (x = 0; x < 3; x ++)
temp[y][x] = b[y][0] * a[0][x] +
b[y][1] * a[1][x] +
b[y][2] * a[2][x];
memcpy(c, temp, sizeof(temp));
}
static void
saturate(float mat[3][3],
float sat)
{
float smat[3][3];
smat[0][0] = (1.0 - sat) * 0.3086 + sat;
smat[0][1] = (1.0 - sat) * 0.3086;
smat[0][2] = (1.0 - sat) * 0.3086;
smat[1][0] = (1.0 - sat) * 0.6094;
smat[1][1] = (1.0 - sat) * 0.6094 + sat;
smat[1][2] = (1.0 - sat) * 0.6094;
smat[2][0] = (1.0 - sat) * 0.0820;
smat[2][1] = (1.0 - sat) * 0.0820;
smat[2][2] = (1.0 - sat) * 0.0820 + sat;
mult(smat, mat, mat);
}
static void
xform(float mat[3][3],
float x,
float y,
float z,
float *tx,
float *ty,
float *tz)
{
*tx = x * mat[0][0] + y * mat[1][0] + z * mat[2][0];
*ty = x * mat[0][1] + y * mat[1][1] + z * mat[2][1];
*tz = x * mat[0][2] + y * mat[1][2] + z * mat[2][2];
}
static void
xrotate(float mat[3][3],
float rs,
float rc)
{
float rmat[3][3];
rmat[0][0] = 1.0;
rmat[0][1] = 0.0;
rmat[0][2] = 0.0;
rmat[1][0] = 0.0;
rmat[1][1] = rc;
rmat[1][2] = rs;
rmat[2][0] = 0.0;
rmat[2][1] = -rs;
rmat[2][2] = rc;
mult(rmat, mat, mat);
}
static void
yrotate(float mat[3][3],
float rs,
float rc)
{
float rmat[3][3];
rmat[0][0] = rc;
rmat[0][1] = 0.0;
rmat[0][2] = -rs;
rmat[1][0] = 0.0;
rmat[1][1] = 1.0;
rmat[1][2] = 0.0;
rmat[2][0] = rs;
rmat[2][1] = 0.0;
rmat[2][2] = rc;
mult(rmat,mat,mat);
}
static void
zrotate(float mat[3][3],
float rs,
float rc)
{
float rmat[3][3];
rmat[0][0] = rc;
rmat[0][1] = rs;
rmat[0][2] = 0.0;
rmat[1][0] = -rs;
rmat[1][1] = rc;
rmat[1][2] = 0.0;
rmat[2][0] = 0.0;
rmat[2][1] = 0.0;
rmat[2][2] = 1.0;
mult(rmat,mat,mat);
}
static void
zshear(float mat[3][3],
float dx,
float dy)
{
float smat[3][3];
smat[0][0] = 1.0;
smat[0][1] = 0.0;
smat[0][2] = dx;
smat[1][0] = 0.0;
smat[1][1] = 1.0;
smat[1][2] = dy;
smat[2][0] = 0.0;
smat[2][1] = 0.0;
smat[2][2] = 1.0;
mult(smat, mat, mat);
}