home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
vsiftp.vmssoftware.com
/
VSIPUBLIC@vsiftp.vmssoftware.com.tar
/
FREEWARE
/
FREEWARE40.ZIP
/
xpaint-244
/
src
/
iprocess.c
< prev
next >
Wrap
C/C++ Source or Header
|
1996-07-28
|
36KB
|
1,524 lines
/* +-------------------------------------------------------------------+ */
/* | Copyright 1993, David Koblas (koblas@netcom.com) | */
/* | Copyright 1995, 1996 Torsten Martinsen (bullestock@dk-online.dk) | */
/* | | */
/* | Permission to use, copy, modify, and to distribute this software | */
/* | and its documentation for any purpose is hereby granted without | */
/* | fee, provided that the above copyright notice appear in all | */
/* | copies and that both that copyright notice and this permission | */
/* | notice appear in supporting documentation. There is no | */
/* | representations about the suitability of this software for | */
/* | any purpose. this software is provided "as is" without express | */
/* | or implied warranty. | */
/* | | */
/* +-------------------------------------------------------------------+ */
/* $Id: iprocess.c,v 1.11 1996/05/09 07:12:23 torsten Exp $ */
#include <X11/Intrinsic.h>
#include <X11/StringDefs.h>
#include <X11/Xatom.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "xpaint.h"
#include "image.h"
#include "misc.h"
#include "protocol.h"
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
/*
** Some real image processing
*/
extern struct imageprocessinfo ImgProcessInfo;
#define CLAMP(low, value, high) \
if (value < low) value = low; else if (value > high) value = high
#define DEG2RAD (M_PI/180.0)
#define ConvMatrixSize 3
typedef float ConvMatrix[ConvMatrixSize * ConvMatrixSize];
static int *histogram(Image * input);
static Image *
convolve(Image * input, float *mat, int n,
unsigned char *basePixel, Boolean absFlag)
{
int x, y, xx, yy, xv, yv;
float sum;
unsigned char *p;
unsigned char *op;
float r, g, b;
int ir, ig, ib;
Image *output;
output = ImageNew(input->width, input->height);
op = output->data;
for (y = 0; y < input->height; y++) {
for (x = 0; x < input->width; x++) {
r = g = b = 0;
sum = 0;
for (yy = 0; yy < n; yy++) {
for (xx = 0; xx < n; xx++) {
xv = x + xx - 1;
yv = y + yy - 1;
if (xv < 0 || yv < 0)
continue;
if (xv >= input->width || yv >= input->height)
continue;
p = ImagePixel(input, xv, yv);
r += (float) *p++ * mat[xx + n * yy];
g += (float) *p++ * mat[xx + n * yy];
b += (float) *p * mat[xx + n * yy];
sum += mat[xx + n * yy];
}
}
if (sum <= 0)
sum = 1;
if (absFlag) {
if (r < 0)
r = -r;
if (g < 0)
g = -g;
if (b < 0)
b = -b;
}
ir = r / sum;
ig = g / sum;
ib = b / sum;
if (basePixel) {
ir += basePixel[0];
ig += basePixel[1];
ib += basePixel[2];
}
CLAMP(0, ir, 255);
CLAMP(0, ig, 255);
CLAMP(0, ib, 255);
*op++ = ir;
*op++ = ig;
*op++ = ib;
}
if (y % 16 == 0)
StateTimeStep();
}
return output;
}
/*
** rescale values from 0..255
*/
static void
normalize(Image * image)
{
int i, count;
unsigned char *sp;
unsigned char *ip;
int maxval = 0;
if (image->cmapSize != 0) {
sp = image->cmapData;
count = image->cmapSize * 3;
} else {
sp = image->data;
count = image->width * image->height * 3;
}
for (ip = sp, i = 0; i < count; i++, ip++)
if (*ip > maxval)
maxval = *ip;
if (maxval == 0)
return;
for (ip = sp, i = 0; i < count; i++, ip++)
*ip = ((int) *ip * 255) / maxval;
}
/*
** Convert image into a monochrome image
*/
static void
monochrome(Image * image)
{
int i, count;
unsigned char *sp;
unsigned char *ip;
int v;
if (image->cmapSize != 0) {
sp = image->cmapData;
count = image->cmapSize;
} else {
sp = image->data;
count = image->width * image->height;
}
for (ip = sp, i = 0; i < count; i++) {
v = (ip[0] * 11 + ip[1] * 16 + ip[2] * 5) >> 5; /* pp = .33R+.5G+.17B */
*ip++ = v;
*ip++ = v;
*ip++ = v;
}
}
Image *
ImageSmooth(Image * input)
{
float *mat;
int n, i;
Image *output;
/*
* Build n x n convolution matrix (all 1's)
*/
n = ImgProcessInfo.smoothMaskSize;
mat = malloc(n * n * sizeof(float));
for (i = 0; i < n * n; ++i)
mat[i] = 1.0;
output = convolve(input, mat, n, NULL, False);
free(mat);
return output;
}
Image *
ImageSharpen(Image * input)
{
static ConvMatrix mat =
{
-1, -2, -1,
-2, 20, -2,
-1, -2, -1
};
return convolve(input, mat, ConvMatrixSize, NULL, False);
}
Image *
ImageEdge(Image * input)
{
static ConvMatrix mat =
{
-1, -2, 0,
-2, 0, 2,
0, 2, 1
};
Image *image = convolve(input, mat, ConvMatrixSize, NULL, True);
normalize(image);
return image;
}
Image *
ImageEmbose(Image * input)
{
static ConvMatrix mat =
{
-1, -2, 0,
-2, 0, 2,
0, 2, 1
};
static unsigned char base[3] =
{128, 128, 128};
Image *image = convolve(input, mat, ConvMatrixSize, base, False);
monochrome(image);
normalize(image);
return image;
}
Image *
ImageInvert(Image * input)
{
Image *output = ImageNewCmap(input->width, input->height, input->cmapSize);
int i;
unsigned char *ip, *op;
int count;
/*
** If the input has a colormap, just invert that.
*/
if (input->cmapSize != 0) {
ip = input->cmapData;
op = output->cmapData;
count = input->cmapSize;
memcpy(output->data, input->data,
sizeof(char) * input->scale * input->width * input->height);
} else {
ip = input->data;
op = output->data;
count = input->width * input->height;
}
for (i = 0; i < count; i++) {
*op++ = 255 - *ip++;
*op++ = 255 - *ip++;
*op++ = 255 - *ip++;
}
return output;
}
Image *
ImageOilPaint(Image * input)
{
Image *output = ImageNew(input->width, input->height);
unsigned char *op = output->data;
int x, y, xx, yy, i;
int rVal, gVal, bVal;
int rCnt, gCnt, bCnt;
int rHist[256], gHist[256], bHist[256];
int oilArea_2;
oilArea_2 = ImgProcessInfo.oilArea / 2;
for (y = 0; y < input->height; y++) {
for (x = 0; x < input->width; x++) {
/*
** compute histogram of (on-screen hunk of) n*n
** region centered plane
*/
rCnt = gCnt = bCnt = 0;
rVal = gVal = bVal = 0;
for (i = 0; i < XtNumber(rHist); i++)
rHist[i] = gHist[i] = bHist[i] = 0;
for (yy = y - oilArea_2; yy < y + oilArea_2; yy++) {
if (yy < 0 || yy >= input->height)
continue;
for (xx = x - oilArea_2; xx < x + oilArea_2; xx++) {
int c, p;
unsigned char *rgb;
if (xx < 0 || xx >= input->width)
continue;
rgb = ImagePixel(input, xx, yy);
if ((c = ++rHist[(p = rgb[0]) / 4]) > rCnt) {
rVal = p;
rCnt = c;
}
if ((c = ++gHist[(p = rgb[1]) / 4]) > gCnt) {
gVal = p;
gCnt = c;
}
if ((c = ++bHist[(p = rgb[2]) / 4]) > bCnt) {
bVal = p;
bCnt = c;
}
}
}
*op++ = rVal;
*op++ = gVal;
*op++ = bVal;
}
if (y % 16 == 0)
StateTimeStep();
}
return output;
}
/*
* Add random noise to an image.
*/
Image *
ImageAddNoise(Image * input)
{
unsigned char *p;
int x, y, r, g, b, ddelta;
Image *output = ImageNew(input->width, input->height);
unsigned char *op = output->data;
ddelta = ImgProcessInfo.noiseDelta * 2;
for (y = 0; y < input->height; y++) {
for (x = 0; x < input->width; x++) {
p = ImagePixel(input, x, y);
r = p[0] + ddelta * gauss();
CLAMP(0, r, 255);
g = p[1] + ddelta * gauss();
CLAMP(0, g, 255);
b = p[2] + ddelta * gauss();
CLAMP(0, b, 255);
*op++ = r;
*op++ = g;
*op++ = b;
}
}
return output;
}
/*
* This function works in-place.
* Because it swaps pixels in the input image, which may be colour mapped
* or not, is has knowledge about the Image format that should probably
* have stayed in image.h. Too bad.
*/
Image *
ImageSpread(Image * input)
{
int w = input->width, h = input->height;
unsigned char *p, *np;
int x, y, minx, miny, maxx, maxy, xn, yn, dist;
dist = ImgProcessInfo.spreadDistance;
for (y = 0; y < h; y++) {
for (x = 0; x < w; x++) {
p = ImagePixel(input, x, y);
/* find random neighbour pixel within +- dist */
minx = x - dist;
if (minx < 0)
minx = 0;
maxx = x + dist;
if (maxx >= w)
maxx = w - 1;
xn = RANDOMI2(minx, maxx);
miny = y - dist;
if (miny < 0)
miny = 0;
maxy = y + dist;
if (maxy >= h)
maxy = h - 1;
yn = RANDOMI2(miny, maxy);
/* swap pixel with neighbour */
np = ImagePixel(input, xn, yn);
if (input->cmapSize == 0) {
unsigned char rn, gn, bn;
rn = np[0];
gn = np[1];
bn = np[2];
np[0] = p[0];
np[1] = p[1];
np[2] = p[2];
p[0] = rn;
p[1] = gn;
p[2] = bn;
} else if (input->cmapSize > 256) {
unsigned short i, *ps, *nps;
ps = &((unsigned short *) input->data)[y * w + x];
nps = &((unsigned short *) input->data)[yn * w + xn];
i = *nps;
*nps = *ps;
*ps = i;
} else { /* colour map of 256 or less */
unsigned char i;
p = &input->data[y * w + x];
np = &input->data[yn * w + xn];
i = *np;
*np = *p;
*p = i;
}
}
}
return input;
}
Image *
ImageBlend(Image * input)
{
int w = input->width, h = input->height;
Image *output = ImageNew(w, h);
unsigned char *op = output->data;
unsigned char *rgb;
int x, y, n, ox, oy, px, py, dx, dy;
int rSum, gSum, bSum, r, g, b;
float diagonal, gradient, ap;
ox = w / 2;
oy = h / 2;
diagonal = ((double) h) / w;
/*
* Compute average of pixels on edge of region. While we are
* at it, we copy the edge to the output as well.
*/
rSum = gSum = bSum = 0;
for (x = 0; x < w; ++x) {
rgb = ImagePixel(input, x, 0);
r = rgb[0];
g = rgb[1];
b = rgb[2];
rSum += r;
gSum += g;
bSum += b;
*op++ = r;
*op++ = g;
*op++ = b;
}
op = ImagePixel(output, 0, h - 1);
for (x = 0; x < w; ++x) {
rgb = ImagePixel(input, x, h - 1);
r = rgb[0];
g = rgb[1];
b = rgb[2];
rSum += r;
gSum += g;
bSum += b;
*op++ = r;
*op++ = g;
*op++ = b;
}
for (y = 1; y < h - 1; ++y) {
rgb = ImagePixel(input, 0, y);
r = rgb[0];
g = rgb[1];
b = rgb[2];
rSum += r;
gSum += g;
bSum += b;
op = ImagePixel(output, 0, y);
*op++ = r;
*op++ = g;
*op++ = b;
rgb = ImagePixel(input, w - 1, y);
r = rgb[0];
g = rgb[1];
b = rgb[2];
rSum += r;
gSum += g;
bSum += b;
op = ImagePixel(output, w - 1, y);
*op++ = r;
*op++ = g;
*op++ = b;
}
n = 2 * (w + h - 2); /* total # of pixels on edge */
rSum /= n;
gSum /= n;
bSum /= n;
op = ImagePixel(output, 1, 1);
/* compute colour of each point in the interior */
for (y = 1; y < h - 1; y++) {
for (x = 1; x < w - 1; x++) {
dx = x - ox;
dy = y - oy;
if (dx == 0 && dy == 0) { /* this is the centre point */
r = rSum;
g = gSum;
b = bSum;
} else {
if (dx == 0) { /* special case 1 */
px = ox;
py = (dy > 0) ? h - 1 : 0;
} else if (dy == 0) { /* special case 2 */
py = oy;
px = (dx > 0) ? w - 1 : 0;
} else { /* general case */
gradient = ((double) dy) / dx;
if (fabs(gradient) < fabs(diagonal)) {
px = (dx > 0) ? w - 1 : 0;
py = oy + ((px - ox) * dy) / dx;
} else {
py = (dy > 0) ? h - 1 : 0;
px = ox + ((py - oy) * dx) / dy;
}
}
/*
* given O(ox,oy), P(px,py), and A(x,y), compute
* |AO|/|PO|
*/
ap = sqrt((double) (dx * dx) + (double) (dy * dy)) /
sqrt((double) ((px - ox) * (px - ox)) +
(double) ((py - oy) * (py - oy)));
rgb = ImagePixel(input, px, py);
r = (1 - ap) * rSum + ap * rgb[0];
g = (1 - ap) * gSum + ap * rgb[1];
b = (1 - ap) * bSum + ap * rgb[2];
}
*op++ = r;
*op++ = g;
*op++ = b;
}
op += 6; /* skip last on this line and first on next */
}
return output;
}
Image *
ImagePixelize(Image * input)
{
int width = input->width, height = input->height;
Image *output = ImageNew(width, height);
unsigned char *op, *rgb;
int x, y, xx, yy, n;
int rSum, gSum, bSum;
int xsize, ysize;
xsize = ImgProcessInfo.pixelizeXSize;
ysize = ImgProcessInfo.pixelizeYSize;
if (xsize > width)
xsize = width;
if (ysize > height)
ysize = height;
n = xsize * ysize;
for (y = 0; y < height; y += ysize) {
for (x = 0; x < width; x += xsize) {
/*
* compute average of pixels inside megapixel
*/
rSum = gSum = bSum = 0;
for (yy = y; yy < y + ysize; yy++) {
if (yy >= height)
continue;
for (xx = x; xx < x + xsize; xx++) {
if (xx >= width)
continue;
rgb = ImagePixel(input, xx, yy);
rSum += rgb[0];
gSum += rgb[1];
bSum += rgb[2];
}
}
rSum /= n;
gSum /= n;
bSum /= n;
/*
* replace each pixel in megapixel with average
*/
for (yy = y; yy < y + ysize; yy++) {
if (yy >= height)
continue;
for (xx = x; xx < x + xsize; xx++) {
if (xx >= width)
continue;
op = ImagePixel(output, xx, yy);
*op++ = rSum;
*op++ = gSum;
*op = bSum;
}
}
}
if (y % 16 == 0)
StateTimeStep();
}
/* if any incomplete megapixels are left, copy them to the output */
for (x = (width / xsize) * xsize; x < width; ++x)
for (y = 0; y < height; ++y) {
rgb = ImagePixel(input, x, y);
op = ImagePixel(output, x, y);
*op++ = *rgb++;
*op++ = *rgb++;
*op = *rgb;
}
for (y = (height / ysize) * ysize; y < height; ++y)
for (x = 0; x < width; ++x) {
rgb = ImagePixel(input, x, y);
op = ImagePixel(output, x, y);
*op++ = *rgb++;
*op++ = *rgb++;
*op = *rgb;
}
return output;
}
#define SWAP(a, b) { int t = (a); (a) = (b); (b) = t; }
Image *
ImageDespeckle(Image * input)
{
int width = input->width, height = input->height;
Image *output = ImageNew(width, height);
unsigned char *op, *rgb;
int x, y, xx, yy, mask, mask2, i, j, k, l;
int *ra, *ga, *ba;
mask = ImgProcessInfo.despeckleMask;
if (mask > width)
mask = width;
if (mask > height)
mask = height;
mask2 = mask / 2;
/* arrays for storing pixels inside mask */
ra = xmalloc(mask * mask * sizeof(int));
ga = xmalloc(mask * mask * sizeof(int));
ba = xmalloc(mask * mask * sizeof(int));
op = output->data;
for (y = 0; y < height; ++y) {
for (x = 0; x < width; ++x) {
i = 0;
for (yy = MAX(0, y - mask2); yy < MIN(height, y + mask2); ++yy)
for (xx = MAX(0, x - mask2); xx < MIN(width, x + mask2); ++xx) {
rgb = ImagePixel(input, xx, yy);
ra[i] = *rgb++;
ga[i] = *rgb++;
ba[i] = *rgb;
++i;
}
/*
* now find median by (shell-)sorting the arrays and
* picking the center value
*/
for (j = i / 2; j > 0; j = j / 2)
for (k = j; k < i; k++) {
for (l = k - j; l >= 0 && ra[l] > ra[l + j]; l -= j)
SWAP(ra[l], ra[l + j]);
for (l = k - j; l >= 0 && ga[l] > ga[l + j]; l -= j)
SWAP(ga[l], ga[l + j]);
for (l = k - j; l >= 0 && ba[l] > ba[l + j]; l -= j)
SWAP(ba[l], ba[l + j]);
}
if (i & 1) { /* uneven number of data points */
*op++ = ra[i / 2];
*op++ = ga[i / 2];
*op++ = ba[i / 2];
} else { /* even, take average */
*op++ = (ra[i / 2 - 1] + ra[i / 2]) / 2;
*op++ = (ga[i / 2 - 1] + ga[i / 2]) / 2;
*op++ = (ba[i / 2 - 1] + ba[i / 2]) / 2;
}
}
if (y % 16 == 0)
StateTimeStep();
}
free(ra);
free(ga);
free(ba);
return output;
}
/*
* Normalize the contrast of an image.
*/
Image *
ImageNormContrast(Image * input)
{
unsigned char *p;
int *hist;
int x, y, w, h, i, max = 0, cumsum, limit;
Image *output = ImageNew(input->width, input->height);
unsigned char *op = output->data;
int blackpcnt, blackval, whitepcnt, whiteval;
blackpcnt = ImgProcessInfo.contrastB;
whitepcnt = 100 - ImgProcessInfo.contrastW;
hist = histogram(input);
w = input->width;
h = input->height;
/* Find lower knee point in histogram to determine 'black' threshold. */
cumsum = 0;
limit = w * h * blackpcnt / 100;
for (i = 0; i < 256; ++i)
if ((cumsum += hist[i]) > limit)
break;
blackval = i;
/* Likewise for 'white' threshold. */
cumsum = 0;
limit = w * h * whitepcnt / 100;
for (i = 256; i > 0; --i) {
if ((max == 0) && hist[i])
max = i + 1; /* max is index of highest non-zero entry */
if ((cumsum += hist[i]) > limit)
break;
}
whiteval = i;
/*
* Now create a histogram with a linear ramp
* from (blackval, 0) to (whiteval, max)
*/
for (i = 0; i < blackval; ++i)
hist[i] = 0;
for (; i <= whiteval; ++i)
hist[i] = max * (i - blackval) / (whiteval - blackval);
for (; i < 256; ++i)
hist[i] = max;
for (y = 0; y < h; y++) {
for (x = 0; x < w; x++) {
p = ImagePixel(input, x, y);
*op++ = hist[*p++];
*op++ = hist[*p++];
*op++ = hist[*p];
}
}
free(hist);
return output;
}
/*
* Build histogram data for the image.
*/
static int *
histogram(Image * input)
{
unsigned char *p;
int x, y, w, h, i, r, g, b, *hist;
w = input->width;
h = input->height;
/* Make a table of 256 zeros plus one sentinel */
hist = xmalloc(257 * sizeof(int));
for (i = 0; i <= 256; ++i)
hist[i] = 0;
/*
* Build histogram of intensity values.
*/
for (y = 0; y < h; y++) {
for (x = 0; x < w; x++) {
p = ImagePixel(input, x, y);
r = *p++;
g = *p++;
b = *p;
i = (r * 169 + g * 256 + b * 87) / 512; /* i = .33R+.5G+.17B */
++hist[i];
}
}
return hist;
}
#if 0
#define HIST_H 120 /* Height of histogram */
#define HIST_W 256 /* Width of histogram */
#define HIST_B 2 /* Horizontal border width for histogram */
#define HIST_R 12 /* Height of ruler below histogram */
#define HIST_T 9 /* Height of tick marks on ruler */
#define LIGHTGREY 200 /* Greyscale level used for bars */
#define MIDGREY 255 /* Greyscale level used for ruler */
#define DARKGREY 120 /* Greyscale level used for background */
/*
* Create an image (width 256, height HIST_H) depicting the histogram data.
*/
Image *
HistogramImage(char *hist)
{
unsigned char *p;
int x, y, i;
int max, e:
unsigned char *ihist;
Image *output = ImageNewGrey(HIST_W + 2 * HIST_B, HIST_H + HIST_R);;
max = 0;
ihist = xmalloc(128 * sizeof(int));
for (i = 0; i < 128; ++i) {
e = ihist[i] = hist[i * 2] + hist[i * 2 + 1];
if (e > max)
max = e;
}
/*
* Now create histogram image.
* Only display 128 bins to make the histogram easier to read.
* Leave a border of a few pixels in each side.
*/
memset(output->data, DARKGREY, (HIST_W + 2 * HIST_B) * (HIST_H + HIST_R));
for (x = 0; x < 128; ++x) {
i = ihist[x] * (HIST_H - 1) / max;
for (y = HIST_H - 1 - i; y < HIST_H - 1; y++) {
p = output->data + y * (HIST_W + 2 * HIST_B) + x * 2 + HIST_B;
*p++ = LIGHTGREY;
*p = LIGHTGREY;
}
}
/* Ruler */
memset(output->data + HIST_H * (HIST_W + 2 * HIST_B) + HIST_B,
MIDGREY, HIST_W);
for (x = 0; x < 5; ++x) {
int tick[] =
{0, HIST_W / 4, HIST_W / 2, HIST_W * 3 / 4, HIST_W - 1};
for (y = 0; y < HIST_T; y++) {
p = output->data + (y + HIST_H + 1) * (HIST_W + 2 * HIST_B)
+ tick[x] + HIST_B;
*p = MIDGREY;
}
}
free(ihist);
return output;
}
#endif
#ifdef FEATURE_TILT
/*
* Tilt an image.
*/
/*
* Compute the interpolated RGB values of a 1 x 1 pixel centered
* at (x,y) and store these values in op[0] trough op[2].
* The interpolation is based on weighting the RGB values proportionally
* to the overlapping areas.
*/
static void
interpol(float x, float y, Image * input, unsigned char *op)
{
int xl, xh, yl, yh;
float a, b, c, d, A, B, C, D;
unsigned char *plh, *phh, *pll, *phl;
xl = floor(x);
xh = ceil(x);
yl = floor(y);
yh = ceil(y);
pll = ImagePixel(input, xl, yl);
if (xh == xl) {
if (yh == yl) {
/* pixel coincides with one pixel */
*op++ = *pll++;
*op++ = *pll++;
*op++ = *pll++;
} else {
/* pixel overlaps with two pixels on top of each other */
plh = ImagePixel(input, xl, yh);
A = y - yl;
B = yh - y;
*op++ = A * *plh++ + B * *pll++;
*op++ = A * *plh++ + B * *pll++;
*op++ = A * *plh++ + B * *pll++;
}
} else if (yh == yl) {
/* pixel overlaps with two side-by-side pixels */
phl = ImagePixel(input, xh, yl);
A = x - xl;
B = xh - x;
*op++ = A * *phl++ + B * *pll++;
*op++ = A * *phl++ + B * *pll++;
*op++ = A * *phl++ + B * *pll++;
} else {
/* general case: pixel overlaps all four neighbour pixels */
a = xh - x;
b = y - yl;
c = x - xl;
d = yh - y;
A = a * b;
B = b * c;
C = a * d;
D = c * d;
plh = ImagePixel(input, xl, yh);
phl = ImagePixel(input, xh, yl);
phh = ImagePixel(input, xh, yh);
*op++ = A * *plh++ + B * *phh++ + C * *pll++ + D * *phl++;
*op++ = A * *plh++ + B * *phh++ + C * *pll++ + D * *phl++;
*op++ = A * *plh++ + B * *phh++ + C * *pll++ + D * *phl++;
}
}
Image *
ImageTilt(Image * input)
{
int x, y, br, bg, bb;
int width = input->width, height = input->height;
Image *output;
unsigned char *op;
float xt, yt;
float X1, Y1, X2, Y2;
X1 = width * ImgProcessInfo.tiltX1 / 100;
Y1 = height * ImgProcessInfo.tiltY1 / 100;
X2 = width * ImgProcessInfo.tiltX2 / 100;
Y2 = height * ImgProcessInfo.tiltY2 / 100;
if (((Y1 >= 0) && (Y1 < height)) || ((X2 >= 0) && (X2 < width)))
return input;
output = ImageNew(width, height);
/* Get RGB values of background colour */
br = ImgProcessInfo.background->red / 256;
bg = ImgProcessInfo.background->green / 256;
bb = ImgProcessInfo.background->blue / 256;
op = output->data;
for (x = 0; x < width * height; ++x) {
*op++ = br;
*op++ = bg;
*op++ = bb;
}
for (y = 0; y < height; y++) {
for (x = 0; x < width; x++) {
op = ImagePixel(output, x, y);
/* find the input coords corresponding to output coords (x,y) */
xt = (X1 * y - Y1 * x) / (y - Y1);
yt = (-X2 * y + x * Y2) / (x - X2);
if ((xt >= 0) && (xt < width) && (yt >= 0) && (yt < height)) {
#if 0
unsigned char *p;
/* rounding */
p = ImagePixel(input, (int) (xt + 0.5), (int) (yt + 0.5));
*op++ = *p++;
*op++ = *p++;
*op = *p;
#else
interpol(xt, yt, input, op);
#endif
}
}
}
return output;
}
#endif
/*
* Produce the 'solarization' effect seen when exposing a
* photographic film to light during the development process.
* Done here by inverting all pixels above threshold level.
*/
Image *
ImageSolarize(Image * input)
{
Image *output = ImageNewCmap(input->width, input->height, input->cmapSize);
int i;
unsigned char *ip, *op;
int count, limit;
limit = ImgProcessInfo.solarizeThreshold * 255 / 100;
/*
** If the input has a colormap, just tweak that.
*/
if (input->cmapSize != 0) {
ip = input->cmapData;
op = output->cmapData;
count = input->cmapSize;
memcpy(output->data, input->data,
sizeof(char) * input->scale * input->width * input->height);
} else {
ip = input->data;
op = output->data;
count = input->width * input->height;
}
for (i = 0; i < count; i++) {
*op++ = (*ip > limit) ? 255 - *ip++ : *ip++;
*op++ = (*ip > limit) ? 255 - *ip++ : *ip++;
*op++ = (*ip > limit) ? 255 - *ip++ : *ip++;
}
return output;
}
/*
* Colourmap quantization. Hacked from ppmqvga.c, which is part of Netpbm and
* was originally written by Lyle Rains (lrains@netcom.com) and Bill Davidsen
* (davidsen@crd.ge.com).
* You are probably not supposed to understand this.
*/
#define RED_BITS 5
#define GREEN_BITS 6
#define BLUE_BITS 5
#define MAX_RED (1 << RED_BITS)
#define MAX_GREEN (1 << GREEN_BITS)
#define MAX_BLUE (1 << BLUE_BITS)
#define MAXWEIGHT 128
#define STDWEIGHT_DIV (2 << 8)
#define STDWEIGHT_MUL (2 << 10)
#define GAIN 4
static int r, g, b, clutx, rep_threshold, rep_weight, dr, dg, db;
static int *color_cube, ncolors;
static unsigned char *clut;
#define CUBEINDEX(r,g,b) (r)*(MAX_GREEN*MAX_BLUE) + (g)*MAX_BLUE + (b)
static void
diffuse(void)
{
int _7_32nds, _3_32nds, _1_16th;
if (clutx < ncolors) {
if (color_cube[CUBEINDEX(r, g, b)] > rep_threshold) {
clut[clutx * 4 + 0] = ((2 * r + 1) * 256) / (2 * MAX_RED);
clut[clutx * 4 + 1] = ((2 * g + 1) * 256) / (2 * MAX_GREEN);
clut[clutx * 4 + 2] = ((2 * b + 1) * 256) / (2 * MAX_BLUE);
++clutx;
color_cube[CUBEINDEX(r, g, b)] -= rep_weight;
}
_7_32nds = (7 * color_cube[CUBEINDEX(r, g, b)]) / 32;
_3_32nds = (3 * color_cube[CUBEINDEX(r, g, b)]) / 32;
_1_16th = color_cube[CUBEINDEX(r, g, b)] - 3 * (_7_32nds + _3_32nds);
color_cube[CUBEINDEX(r, g, b)] = 0;
/* spread error evenly in color space. */
color_cube[CUBEINDEX(r, g, b + db)] += _7_32nds;
color_cube[CUBEINDEX(r, g + dg, b)] += _7_32nds;
color_cube[CUBEINDEX(r + dr, g, b)] += _7_32nds;
color_cube[CUBEINDEX(r, g + dg, b + db)] += _3_32nds;
color_cube[CUBEINDEX(r + dr, g, b + db)] += _3_32nds;
color_cube[CUBEINDEX(r + dr, g + dg, b)] += _3_32nds;
color_cube[CUBEINDEX(r + dr, g + dg, b + db)] += _1_16th;
/*
* Conserve the error at edges if possible
* (which it is, except the last pixel)
*/
if (color_cube[CUBEINDEX(r, g, b)] != 0) {
if (dg != 0)
color_cube[CUBEINDEX(r, g + dg, b)] +=
color_cube[CUBEINDEX(r, g, b)];
else if (dr != 0)
color_cube[CUBEINDEX(r + dr, g, b)] +=
color_cube[CUBEINDEX(r, g, b)];
else if (db != 0)
color_cube[CUBEINDEX(r, g, b) + db] +=
color_cube[CUBEINDEX(r, g, b)];
else
fprintf(stderr, "lost error term\n");
}
}
color_cube[CUBEINDEX(r, g, b)] = -1;
}
/*
** Find representative color nearest to requested color. Check color cube
** for a cached color index. If not cached, compute nearest and cache result.
*/
static int
nearest_color(unsigned char *p)
{
register unsigned char *test;
register unsigned i;
unsigned long min_dist_sqd, dist_sqd;
int nearest = 0;
int *cache;
int r, g, b;
r = *p++;
g = *p++;
b = *p;
cache = &(color_cube[CUBEINDEX((r << RED_BITS) / 256,
(g << GREEN_BITS) / 256,
(b << BLUE_BITS) / 256)]);
if (*cache >= 0)
return *cache;
min_dist_sqd = ~0;
for (i = 0; i < ncolors; ++i) {
test = &clut[i * 4];
dist_sqd =
3 * (r - test[0]) * (r - test[0]) +
4 * (g - test[1]) * (g - test[1]) +
2 * (b - test[2]) * (b - test[2]);
if (dist_sqd < min_dist_sqd) {
nearest = i;
min_dist_sqd = dist_sqd;
}
}
return (*cache = nearest);
}
Image *
ImageQuantize(Image * input)
{
int w = input->width, h = input->height;
Image *output = ImageNew(w, h);
unsigned char *op = output->data;
int *weight_convert;
int k, x, y, i, j, nearest;
int total_weight, cum_weight[MAX_GREEN];
int *erropt;
int *errP;
unsigned char *p, *clutP;
ncolors = ImgProcessInfo.quantizeColors;
/* Make a color cube, initialized to zeros */
color_cube = calloc(MAX_RED * MAX_GREEN * MAX_BLUE, sizeof(int));
clut = calloc(ncolors * 4, sizeof(int));
erropt = calloc(ncolors * 4, sizeof(int));
clutx = 0;
/* Count all occurrances of each color */
for (y = 0; y < h; ++y)
for (x = 0; x < w; ++x) {
p = ImagePixel(input, x, y);
r = p[0] / (256 / MAX_RED);
g = p[1] / (256 / MAX_GREEN);
b = p[2] / (256 / MAX_BLUE);
++color_cube[CUBEINDEX(r, g, b)];
}
/* Initialize logarithmic weighting table */
weight_convert = malloc(MAXWEIGHT * sizeof(int));
weight_convert[0] = 0;
for (i = 1; i < MAXWEIGHT; ++i) {
weight_convert[i] = (int) (100.0 * log((double) i));
}
k = w * h;
if ((k /= STDWEIGHT_DIV) == 0)
k = 1;
total_weight = i = 0;
for (g = 0; g < MAX_GREEN; ++g) {
for (r = 0; r < MAX_RED; ++r) {
for (b = 0; b < MAX_BLUE; ++b) {
register int weight;
/* Normalize the weights, independent of picture size. */
weight = color_cube[CUBEINDEX(r, g, b)] * STDWEIGHT_MUL;
weight /= k;
if (weight)
++i;
if (weight >= MAXWEIGHT)
weight = MAXWEIGHT - 1;
total_weight += (color_cube[CUBEINDEX(r, g, b)]
= weight_convert[weight]);
}
}
cum_weight[g] = total_weight;
}
rep_weight = total_weight / ncolors;
/* Magic foo-foo dust here. What IS the correct way to select threshold? */
rep_threshold = total_weight * (28 + 110000 / i) / 95000;
/*
* Do a 3-D error diffusion dither on the data in the color cube
* to select the representative colors. Do the dither back and forth in
* such a manner that all the error is conserved (none lost at the edges).
*/
dg = 1;
for (g = 0; g < MAX_GREEN; ++g) {
dr = 1;
for (r = 0; r < MAX_RED; ++r) {
db = 1;
for (b = 0; b < MAX_BLUE - 1; ++b)
diffuse();
db = 0;
diffuse();
++b;
if (++r == MAX_RED - 1)
dr = 0;
db = -1;
while (--b > 0)
diffuse();
db = 0;
diffuse();
}
/* Modify threshold to keep rep points proportionally distributed */
if ((j = clutx - (ncolors * cum_weight[g]) / total_weight) != 0)
rep_threshold += j * GAIN;
if (++g == MAX_GREEN - 1)
dg = 0;
dr = -1;
while (r-- > 0) {
db = 1;
for (b = 0; b < MAX_BLUE - 1; ++b)
diffuse();
db = 0;
diffuse();
++b;
if (--r == 0)
dr = 0;
db = -1;
while (--b > 0)
diffuse();
db = 0;
diffuse();
}
/* Modify threshold to keep rep points proportionally distributed */
if ((j = clutx - (ncolors * cum_weight[g]) / total_weight) != 0)
rep_threshold += j * GAIN;
}
/*
* Check the error associated with the use of each color, and
* change the value of the color to minimize the error.
*/
for (y = 0; y < h; ++y) {
for (x = 0; x < w; ++x) {
p = ImagePixel(input, x, y);
nearest = nearest_color(p);
errP = &erropt[nearest * 4];
clutP = &clut[nearest * 4];
errP[0] += *p++ - clutP[0];
errP[1] += *p++ - clutP[1];
errP[2] += *p - clutP[2];
++errP[3];
}
}
for (i = 0; i < ncolors; ++i) {
clutP = &clut[i * 4];
errP = &erropt[i * 4];
j = errP[3];
if (j > 0)
j *= 4;
else if (j == 0)
j = 1;
clutP[0] += (errP[0] / j) * 4;
clutP[1] += (errP[1] / j) * 4;
clutP[2] += (errP[2] / j) * 4;
}
/* Reset the color cache. */
for (i = 0; i < MAX_RED * MAX_GREEN * MAX_BLUE; ++i)
color_cube[i] = -1;
/*
* Map the colors in the image to their closest match in the new colormap.
*/
for (y = 0; y < h; ++y) {
for (x = 0; x < w; ++x) {
p = ImagePixel(input, x, y);
nearest = nearest_color(p);
clutP = &clut[nearest * 4];
*op++ = *clutP++;
*op++ = *clutP++;
*op++ = *clutP;
}
}
free(clut);
free(erropt);
free(weight_convert);
return output;
}
/*
* Convert an image to grey scale.
*/
Image *
ImageGrey(Image * input)
{
Image *output = ImageNew(input->width, input->height);
unsigned char *p, *op = output->data;
int x, y, v;
for (y = 0; y < input->height; y++) {
for (x = 0; x < input->width; x++) {
p = ImagePixel(input, x, y);
v = (p[0] * 11 + p[1] * 16 + p[2] * 5) >> 5; /* .33R + .5G + .17B */
*op++ = v;
*op++ = v;
*op++ = v;
}
}
return output;
}
#define GRAY(p) (((p)[0]*169 + (p)[1]*256 + (p)[2]*87)/512)
#define SQ(x) ((x)*(x))
/*
* Directional filter, according to the algorithm described on p. 60 of:
* _Algorithms_for_Graphics_and_Image_Processing_. Theo Pavlidis.
* Computer Science Press, 1982.
* For each pixel, detect the most prominent edge and apply a filter that
* does not degrade that edge.
*/
Image *
ImageDirectionalFilter(Image * input)
{
unsigned char *p00, *p10, *p_0, *p11, *p__, *p01, *p0_, *p_1, *p1_;
int g00, g10, g_0, g11, g__, g01, g0_, g_1, g1_;
int v0, v45, v90, v135, vmin, theta;
int x, y;
Image *output = ImageNew(input->width, input->height);
unsigned char *op = output->data;
/*
* We don't process the border of the image tto avoid the hassle
* of having do deal with boundary conditions. Hopefully no one
* will notice.
*/
/* copy first row unchanged */
for (x = 0; x < input->width; x++) {
p00 = ImagePixel(input, x, 0);
*op++ = *p00++;
*op++ = *p00++;
*op++ = *p00++;
}
for (y = 1; y < input->height - 1; y++) {
/* copy first column unchanged */
p00 = ImagePixel(input, 0, y);
*op++ = *p00++;
*op++ = *p00++;
*op++ = *p00++;
for (x = 1; x < input->width - 1; x++) {
/* find values of pixel and all neighbours */
p00 = ImagePixel(input, x, y);
p10 = ImagePixel(input, x + 1, y);
p_0 = ImagePixel(input, x - 1, y);
p11 = ImagePixel(input, x + 1, y + 1);
p__ = ImagePixel(input, x - 1, y - 1);
p01 = ImagePixel(input, x, y + 1);
p0_ = ImagePixel(input, x, y - 1);
p_1 = ImagePixel(input, x - 1, y + 1);
p1_ = ImagePixel(input, x + 1, y - 1);
/* get grayscale values */
g00 = GRAY(p00);
g01 = GRAY(p01);
g10 = GRAY(p10);
g11 = GRAY(p11);
g0_ = GRAY(p0_);
g_0 = GRAY(p_0);
g__ = GRAY(p__);
g_1 = GRAY(p_1);
g1_ = GRAY(p1_);
/* estimate direction of edge, if any */
v0 = SQ(g00 - g10) + SQ(g00 - g_0);
v45 = SQ(g00 - g11) + SQ(g00 - g__);
v90 = SQ(g00 - g01) + SQ(g00 - g0_);
v135 = SQ(g00 - g_1) + SQ(g00 - g1_);
vmin = MIN(MIN(v0, v45), MIN(v90, v135));
theta = 0;
if (vmin == v45)
theta = 1;
else if (vmin == v90)
theta = 2;
else if (vmin == v135)
theta = 3;
/* apply filtering according to direction of edge */
switch (theta) {
case 0: /* 0 degrees */
*op++ = (*p_0++ + *p00++ + *p10++) / 3;
*op++ = (*p_0++ + *p00++ + *p10++) / 3;
*op++ = (*p_0++ + *p00++ + *p10++) / 3;
break;
case 1: /* 45 degrees */
*op++ = (*p__++ + *p00++ + *p11++) / 3;
*op++ = (*p__++ + *p00++ + *p11++) / 3;
*op++ = (*p__++ + *p00++ + *p11++) / 3;
break;
case 2: /* 90 degrees */
*op++ = (*p0_++ + *p00++ + *p01++) / 3;
*op++ = (*p0_++ + *p00++ + *p01++) / 3;
*op++ = (*p0_++ + *p00++ + *p01++) / 3;
break;
case 3: /* 135 degrees */
*op++ = (*p1_++ + *p00++ + *p_1++) / 3;
*op++ = (*p1_++ + *p00++ + *p_1++) / 3;
*op++ = (*p1_++ + *p00++ + *p_1++) / 3;
break;
}
}
/* copy last column unchanged */
p00 = ImagePixel(input, x, y);
*op++ = *p00++;
*op++ = *p00++;
*op++ = *p00++;
}
/* copy last row unchanged */
for (x = 0; x < input->width; x++) {
p00 = ImagePixel(input, x, y);
*op++ = *p00++;
*op++ = *p00++;
*op++ = *p00++;
}
return output;
}
#if 0
#define ISFOREGND(p) ((p) ? \
(deltaR-deltaRV <= p[0]) && (p[0] <= deltaR+deltaRV) && \
(deltaG-deltaGV <= p[1]) && (p[1] <= deltaG+deltaGV) && \
(deltaB-deltaBV <= p[2]) && (p[2] <= deltaB+deltaBV) : 1)
/*
* Thicken an image.
*/
Image *
ImageThicken(Image * input)
{
unsigned char *p00, *p10, *p_0, *p01, *p0_;
int x, y, width, height, br, bg, bb;
int deltaR, deltaRV, deltaG, deltaGV, deltaB, deltaBV;
Image *output;
unsigned char *op;
width = input->width;
height = input->height;
output = ImageNew(width, height);
op = output->data;
/* Get RGB values of background colour */
br = ImgProcessInfo.background->red / 256;
bg = ImgProcessInfo.background->green / 256;
bb = ImgProcessInfo.background->blue / 256;
for (y = 0; y < height; y++)
for (x = 0; x < width; x++) {
/* find values of pixel and all d-neighbours */
p00 = ImagePixel(input, x, y);
p10 = x < width - 1 ? ImagePixel(input, x + 1, y) : NULL;
p_0 = x > 0 ? ImagePixel(input, x - 1, y) : NULL;
p01 = y < height - 1 ? ImagePixel(input, x, y + 1) : NULL;
p0_ = y > 0 ? ImagePixel(input, x, y - 1) : NULL;
if (ISFOREGND(p00)) {
*op++ = *p00++;
*op++ = *p00++;
*op++ = *p00++;
} else {
*op++ = br;
*op++ = bg;
*op++ = bb;
}
}
return output;
}
#endif