Matrix Operations for Image Processing

Paul Haeberli

Nov 1993

Introduction

Four by four matrices are commonly used to transform geometry for 3D rendering. These matrices may also be used to transform RGB colors, to scale RGB colors, and to control hue, saturation and contrast. The most important advantage of using matrices is that any number of color transformations can be composed using standard matrix multiplication.

Please note that for these operations to be correct, we really must operate on linear brightness values. If the input image is in a non-linear brightness space RGB colors must be transformed into a linear space before these matrix operations are used.

Color Transformation

RGB colors are transformed by a four by four matrix as shown here:

    xformrgb(mat,r,g,b,tr,tg,tb)
    float mat[4][4];
    float r,g,b;
    float *tr,*tg,*tb;
    {
        *tr = r*mat[0][0] + g*mat[1][0] +
		    b*mat[2][0] + mat[3][0];
        *tg = r*mat[0][1] + g*mat[1][1] +
		    b*mat[2][1] + mat[3][1];
        *tb = r*mat[0][2] + g*mat[1][2] +
		    b*mat[2][2] + mat[3][2];
    }

The Identity

This is the identity matrix:

    float mat[4][4] = {
        1.0,    0.0,    0.0,    0.0,
        0.0,    1.0,    0.0,    0.0,
        0.0,    0.0,    1.0,    0.0,
        0.0,    0.0,    0.0,    1.0,
    };

Transforming colors by the identity matrix will leave them unchanged.

Changing Brightness

To scale RGB colors a matrix like this is used:

    float mat[4][4] = {
        rscale, 0.0,    0.0,    0.0,
        0.0,    gscale, 0.0,    0.0,
        0.0,    0.0,    bscale, 0.0,
        0.0,    0.0,    0.0,    1.0,
    };

Where rscale, gscale, and bscale specify how much to scale the r, g, and b components of colors. This can be used to alter the color balance of an image.

In effect, this calculates:

	tr = r*rscale;
	tg = g*gscale;
	tb = b*bscale;

Modifying Saturation

Converting to Luminance

To convert a color image into a black and white image, this matrix is used:

    float mat[4][4] = {
        rwgt,   rwgt,   rwgt,   0.0,
        gwgt,   gwgt,   gwgt,   0.0,
        bwgt,   bwgt,   bwgt,   0.0,
        0.0,    0.0,    0.0,    1.0,
    };

Where rwgt is 0.3086, gwgt is 0.6094, and bwgt is 0.0820. This is the luminance vector. Notice here that we do not use the standard NTSC weights of 0.299, 0.587, and 0.114. The NTSC weights are only applicable to RGB colors in a gamma 2.2 color space. For linear RGB colors the values above are better.

In effect, this calculates:

	tr = r*rwgt + g*gwgt + b*bwgt;
	tg = r*rwgt + g*gwgt + b*bwgt;
	tb = r*rwgt + g*gwgt + b*bwgt;

Modifying Saturation

To saturate RGB colors, this matrix is used:

     float mat[4][4] = {
        a,      b,      c,      0.0,
        d,      e,      f,      0.0,
        g,      h,      i,      0.0,
        0.0,    0.0,    0.0,    1.0,
    };

Where the constants are derived from the saturation value s as shown below:

    a = (1.0-s)*rwgt + s;
    b = (1.0-s)*rwgt;
    c = (1.0-s)*rwgt;
    d = (1.0-s)*gwgt;
    e = (1.0-s)*gwgt + s;
    f = (1.0-s)*gwgt;
    g = (1.0-s)*bwgt;
    h = (1.0-s)*bwgt;
    i = (1.0-s)*bwgt + s;

One nice property of this saturation matrix is that the luminance of input RGB colors is maintained. This matrix can also be used to complement the colors in an image by specifying a saturation value of -1.0.

Notice that when s is set to 0.0, the matrix is exactly the "convert to luminance" matrix described above. When s is set to 1.0 the matrix becomes the identity. All saturation matrices can be derived by interpolating between or extrapolating beyond these two matrices.

This is discussed in more detail in the note on Image Processing By Interpolation and Extrapolation.

Applying Offsets to Color Components

To offset the r, g, and b components of colors in an image this matrix is used:

    float mat[4][4] = {
        1.0,    0.0,    0.0,    0.0,
        0.0,    1.0,    0.0,    0.0,
        0.0,    0.0,    1.0,    0.0,
        roffset,goffset,boffset,1.0,
    };

This can be used along with color scaling to alter the contrast of RGB images.

Simple Hue Rotation

To rotate the hue, we perform a 3D rotation of RGB colors about the diagonal vector [1.0 1.0 1.0]. The transformation matrix is derived as shown here:

If we have functions:

• identmat(mat)
  • that creates an identity matrix.
• xrotatemat(mat,rsin,rcos)
  • that multiplies a matrix that rotates about the x (red) axis.
• yrotatemat(mat,rsin,rcos)
  • that multiplies a matrix that rotates about the y (green) axis.
• zrotatemat(mat,rsin,rcos)
  • that multiplies a matrix that rotates about the z (blue) axis.

Then a matrix that rotates about the 1.0,1.0,1.0 diagonal can be constructed like this:

First we make an identity matrix

    identmat(mat);

Rotate the grey vector into positive Z

    mag = sqrt(2.0);
    xrs = 1.0/mag;
    xrc = 1.0/mag;
    xrotatemat(mat,xrs,xrc);

    mag = sqrt(3.0);
    yrs = -1.0/mag;
    yrc = sqrt(2.0)/mag;
    yrotatemat(mat,yrs,yrc);

Rotate the hue

    zrs = sin(rot*PI/180.0);
    zrc = cos(rot*PI/180.0);
    zrotatemat(mat,zrs,zrc);

Rotate the grey vector back into place

    yrotatemat(mat,-yrs,yrc);
    xrotatemat(mat,-xrs,xrc);

The resulting matrix will rotate the hue of the input RGB colors. A rotation of 120.0 degrees will exactly map Red into Green, Green into Blue and Blue into Red. This transformation has one problem, however, the luminance of the input colors is not preserved. This can be fixed with the following refinement:

Hue Rotation While Preserving Luminance

We make an identity matrix

   identmat(mmat);

Rotate the grey vector into positive Z

    mag = sqrt(2.0);
    xrs = 1.0/mag;
    xrc = 1.0/mag;
    xrotatemat(mmat,xrs,xrc);
    mag = sqrt(3.0);
    yrs = -1.0/mag;
    yrc = sqrt(2.0)/mag;
    yrotatemat(mmat,yrs,yrc);
    matrixmult(mmat,mat,mat);

Shear the space to make the luminance plane horizontal

    xformrgb(mmat,rwgt,gwgt,bwgt,&lx;,&ly;,&lz;);
    zsx = lx/lz;
    zsy = ly/lz;
    zshearmat(mat,zsx,zsy);

Rotate the hue

    zrs = sin(rot*PI/180.0);
    zrc = cos(rot*PI/180.0);
    zrotatemat(mat,zrs,zrc);

Unshear the space to put the luminance plane back

    zshearmat(mat,-zsx,-zsy);

Rotate the grey vector back into place

    yrotatemat(mat,-yrs,yrc);
    xrotatemat(mat,-xrs,xrc);

Conclusion

I've presented several matrix transformations that may be applied to RGB colors. Each color transformation is represented by a 4 by 4 matrix, similar to matrices commonly used to transform 3D geometry.

Example C code that demonstrates these concepts is provided for your enjoyment (below).

These transformations allow us to adjust image contrast, brightness, hue and saturation individually. In addition, color matrix transformations concatenate in a way similar to geometric transformations. Any sequence of operations can be combined into a single matrix using matrix multiplication.

/*
 *	matrix - 
 *		Use 4x4 matricies to process color images.
 *
 *	To compile:
	    cc matrix.c -o matrix -lgutil -limage -lgl -lm
 *
 *				Paul Haeberli - 1993
 */
#include "math.h"
#include "stdio.h"

#define RLUM    (0.3086)
#define GLUM    (0.6094)
#define BLUM    (0.0820)

#define OFFSET_R        3
#define OFFSET_G        2
#define OFFSET_B        1
#define OFFSET_A        0

/* 
 *	printmat -	
 *		print a 4 by 4 matrix
 */
printmat(mat)
float mat[4][4];
{
    int x, y;

    fprintf(stderr,"\n");
    for(y=0; y<4; y++) {
	for(x=0; x<4; x++) 
	   fprintf(stderr,"%f ",mat[y][x]);
       	fprintf(stderr,"\n");
    }
    fprintf(stderr,"\n");
}

/* 
 *	applymatrix -	
 *		use a matrix to transform colors.
 */
applymatrix(lptr,mat,n)
unsigned long *lptr;
float mat[4][4];
int n;
{
    int ir, ig, ib, r, g, b;
    unsigned char *cptr;

    cptr = (unsigned char *)lptr;
    while(n--) {
	ir = cptr[OFFSET_R];
	ig = cptr[OFFSET_G];
	ib = cptr[OFFSET_B];
	r = ir*mat[0][0] + ig*mat[1][0] + ib*mat[2][0] + mat[3][0];
	g = ir*mat[0][1] + ig*mat[1][1] + ib*mat[2][1] + mat[3][1];
	b = ir*mat[0][2] + ig*mat[1][2] + ib*mat[2][2] + mat[3][2];
	if(r<0) r = 0;
	if(r>255) r = 255;
	if(g<0) g = 0;
	if(g>255) g = 255;
	if(b<0) b = 0;
	if(b>255) b = 255;
	cptr[OFFSET_R] = r;
	cptr[OFFSET_G] = g;
	cptr[OFFSET_B] = b;
	cptr += 4;
    }
}

/* 
 *	matrixmult -	
 *		multiply two matricies
 */
matrixmult(a,b,c)
float a[4][4], b[4][4], c[4][4];
{
    int x, y;
    float temp[4][4];

    for(y=0; y<4 ; y++)
        for(x=0 ; x<4 ; x++) {
            temp[y][x] = b[y][0] * a[0][x]
                       + b[y][1] * a[1][x]
                       + b[y][2] * a[2][x]
                       + b[y][3] * a[3][x];
        }
    for(y=0; y<4; y++)
        for(x=0; x<4; x++)
            c[y][x] = temp[y][x];
}

/* 
 *	identmat -	
 *		make an identity matrix
 */
identmat(matrix)
float *matrix;
{
    *matrix++ = 1.0;    /* row 1        */
    *matrix++ = 0.0;
    *matrix++ = 0.0;
    *matrix++ = 0.0;
    *matrix++ = 0.0;    /* row 2        */
    *matrix++ = 1.0;
    *matrix++ = 0.0;
    *matrix++ = 0.0;
    *matrix++ = 0.0;    /* row 3        */
    *matrix++ = 0.0;
    *matrix++ = 1.0;
    *matrix++ = 0.0;
    *matrix++ = 0.0;    /* row 4        */
    *matrix++ = 0.0;
    *matrix++ = 0.0;
    *matrix++ = 1.0;
}

/* 
 *	xformpnt -	
 *		transform a 3D point using a matrix
 */
xformpnt(matrix,x,y,z,tx,ty,tz)
float matrix[4][4];
float x,y,z;
float *tx,*ty,*tz;
{
    *tx = x*matrix[0][0] + y*matrix[1][0] + z*matrix[2][0] + matrix[3][0];
    *ty = x*matrix[0][1] + y*matrix[1][1] + z*matrix[2][1] + matrix[3][1];
    *tz = x*matrix[0][2] + y*matrix[1][2] + z*matrix[2][2] + matrix[3][2];
}

/* 
 *	cscalemat -	
 *		make a color scale marix
 */
cscalemat(mat,rscale,gscale,bscale)
float mat[4][4];
float rscale, gscale, bscale;
{
    float mmat[4][4];

    mmat[0][0] = rscale;
    mmat[0][1] = 0.0;
    mmat[0][2] = 0.0;
    mmat[0][3] = 0.0;

    mmat[1][0] = 0.0;
    mmat[1][1] = gscale;
    mmat[1][2] = 0.0;
    mmat[1][3] = 0.0;


    mmat[2][0] = 0.0;
    mmat[2][1] = 0.0;
    mmat[2][2] = bscale;
    mmat[2][3] = 0.0;

    mmat[3][0] = 0.0;
    mmat[3][1] = 0.0;
    mmat[3][2] = 0.0;
    mmat[3][3] = 1.0;
    matrixmult(mmat,mat,mat);
}

/* 
 *	lummat -	
 *		make a luminance marix
 */
lummat(mat)
float mat[4][4];
{
    float mmat[4][4];
    float rwgt, gwgt, bwgt;

    rwgt = RLUM;
    gwgt = GLUM;
    bwgt = BLUM;
    mmat[0][0] = rwgt;
    mmat[0][1] = rwgt;
    mmat[0][2] = rwgt;
    mmat[0][3] = 0.0;

    mmat[1][0] = gwgt;
    mmat[1][1] = gwgt;
    mmat[1][2] = gwgt;
    mmat[1][3] = 0.0;

    mmat[2][0] = bwgt;
    mmat[2][1] = bwgt;
    mmat[2][2] = bwgt;
    mmat[2][3] = 0.0;

    mmat[3][0] = 0.0;
    mmat[3][1] = 0.0;
    mmat[3][2] = 0.0;
    mmat[3][3] = 1.0;
    matrixmult(mmat,mat,mat);
}

/* 
 *	saturatemat -	
 *		make a saturation marix
 */
saturatemat(mat,sat)
float mat[4][4];
float sat;
{
    float mmat[4][4];
    float a, b, c, d, e, f, g, h, i;
    float rwgt, gwgt, bwgt;

    rwgt = RLUM;
    gwgt = GLUM;
    bwgt = BLUM;

    a = (1.0-sat)*rwgt + sat;
    b = (1.0-sat)*rwgt;
    c = (1.0-sat)*rwgt;
    d = (1.0-sat)*gwgt;
    e = (1.0-sat)*gwgt + sat;
    f = (1.0-sat)*gwgt;
    g = (1.0-sat)*bwgt;
    h = (1.0-sat)*bwgt;
    i = (1.0-sat)*bwgt + sat;
    mmat[0][0] = a;
    mmat[0][1] = b;
    mmat[0][2] = c;
    mmat[0][3] = 0.0;

    mmat[1][0] = d;
    mmat[1][1] = e;
    mmat[1][2] = f;
    mmat[1][3] = 0.0;

    mmat[2][0] = g;
    mmat[2][1] = h;
    mmat[2][2] = i;
    mmat[2][3] = 0.0;

    mmat[3][0] = 0.0;
    mmat[3][1] = 0.0;
    mmat[3][2] = 0.0;
    mmat[3][3] = 1.0;
    matrixmult(mmat,mat,mat);
}

/* 
 *	offsetmat -	
 *		offset r, g, and b
 */
offsetmat(mat,roffset,goffset,boffset)
float mat[4][4];
float roffset, goffset, boffset;
{
    float mmat[4][4];

    mmat[0][0] = 1.0;
    mmat[0][1] = 0.0;
    mmat[0][2] = 0.0;
    mmat[0][3] = 0.0;

    mmat[1][0] = 0.0;
    mmat[1][1] = 1.0;
    mmat[1][2] = 0.0;
    mmat[1][3] = 0.0;

    mmat[2][0] = 0.0;
    mmat[2][1] = 0.0;
    mmat[2][2] = 1.0;
    mmat[2][3] = 0.0;

    mmat[3][0] = roffset;
    mmat[3][1] = goffset;
    mmat[3][2] = boffset;
    mmat[3][3] = 1.0;
    matrixmult(mmat,mat,mat);
}

/* 
 *	xrotate -	
 *		rotate about the x (red) axis
 */
xrotatemat(mat,rs,rc)
float mat[4][4];
float rs, rc;
{
    float mmat[4][4];

    mmat[0][0] = 1.0;
    mmat[0][1] = 0.0;
    mmat[0][2] = 0.0;
    mmat[0][3] = 0.0;

    mmat[1][0] = 0.0;
    mmat[1][1] = rc;
    mmat[1][2] = rs;
    mmat[1][3] = 0.0;

    mmat[2][0] = 0.0;
    mmat[2][1] = -rs;
    mmat[2][2] = rc;
    mmat[2][3] = 0.0;

    mmat[3][0] = 0.0;
    mmat[3][1] = 0.0;
    mmat[3][2] = 0.0;
    mmat[3][3] = 1.0;
    matrixmult(mmat,mat,mat);
}

/* 
 *	yrotate -	
 *		rotate about the y (green) axis
 */
yrotatemat(mat,rs,rc)
float mat[4][4];
float rs, rc;
{
    float mmat[4][4];

    mmat[0][0] = rc;
    mmat[0][1] = 0.0;
    mmat[0][2] = -rs;
    mmat[0][3] = 0.0;

    mmat[1][0] = 0.0;
    mmat[1][1] = 1.0;
    mmat[1][2] = 0.0;
    mmat[1][3] = 0.0;

    mmat[2][0] = rs;
    mmat[2][1] = 0.0;
    mmat[2][2] = rc;
    mmat[2][3] = 0.0;

    mmat[3][0] = 0.0;
    mmat[3][1] = 0.0;
    mmat[3][2] = 0.0;
    mmat[3][3] = 1.0;
    matrixmult(mmat,mat,mat);
}

/* 
 *	zrotate -	
 *		rotate about the z (blue) axis
 */
zrotatemat(mat,rs,rc)
float mat[4][4];
float rs, rc;
{
    float mmat[4][4];

    mmat[0][0] = rc;
    mmat[0][1] = rs;
    mmat[0][2] = 0.0;
    mmat[0][3] = 0.0;

    mmat[1][0] = -rs;
    mmat[1][1] = rc;
    mmat[1][2] = 0.0;
    mmat[1][3] = 0.0;

    mmat[2][0] = 0.0;
    mmat[2][1] = 0.0;
    mmat[2][2] = 1.0;
    mmat[2][3] = 0.0;

    mmat[3][0] = 0.0;
    mmat[3][1] = 0.0;
    mmat[3][2] = 0.0;
    mmat[3][3] = 1.0;
    matrixmult(mmat,mat,mat);
}

/* 
 *	zshear -	
 *		shear z using x and y.
 */
zshearmat(mat,dx,dy)
float mat[4][4];
float dx, dy;
{
    float mmat[4][4];

    mmat[0][0] = 1.0;
    mmat[0][1] = 0.0;
    mmat[0][2] = dx;
    mmat[0][3] = 0.0;

    mmat[1][0] = 0.0;
    mmat[1][1] = 1.0;
    mmat[1][2] = dy;
    mmat[1][3] = 0.0;

    mmat[2][0] = 0.0;
    mmat[2][1] = 0.0;
    mmat[2][2] = 1.0;
    mmat[2][3] = 0.0;

    mmat[3][0] = 0.0;
    mmat[3][1] = 0.0;
    mmat[3][2] = 0.0;
    mmat[3][3] = 1.0;
    matrixmult(mmat,mat,mat);
}

/* 
 *	simplehuerotatemat -	
 *		simple hue rotation. This changes luminance 
 */
simplehuerotatemat(mat,rot)
float mat[4][4];
float rot;
{
    float mag;
    float xrs, xrc;
    float yrs, yrc;
    float zrs, zrc;

/* rotate the grey vector into positive Z */
    mag = sqrt(2.0);
    xrs = 1.0/mag;
    xrc = 1.0/mag;
    xrotatemat(mat,xrs,xrc);

    mag = sqrt(3.0);
    yrs = -1.0/mag;
    yrc = sqrt(2.0)/mag;
    yrotatemat(mat,yrs,yrc);

/* rotate the hue */
    zrs = sin(rot*M_PI/180.0);
    zrc = cos(rot*M_PI/180.0);
    zrotatemat(mat,zrs,zrc);

/* rotate the grey vector back into place */
    yrotatemat(mat,-yrs,yrc);
    xrotatemat(mat,-xrs,xrc);
}

/* 
 *	huerotatemat -	
 *		rotate the hue, while maintaining luminance.
 */
huerotatemat(mat,rot)
float mat[4][4];
float rot;
{
    float mmat[4][4];
    float mag;
    float lx, ly, lz;
    float xrs, xrc;
    float yrs, yrc;
    float zrs, zrc;
    float zsx, zsy;

    identmat(mmat);

/* rotate the grey vector into positive Z */
    mag = sqrt(2.0);
    xrs = 1.0/mag;
    xrc = 1.0/mag;
    xrotatemat(mmat,xrs,xrc);
    mag = sqrt(3.0);
    yrs = -1.0/mag;
    yrc = sqrt(2.0)/mag;
    yrotatemat(mmat,yrs,yrc);

/* shear the space to make the luminance plane horizontal */
    xformpnt(mmat,RLUM,GLUM,BLUM,&lx,&ly,&lz);
    zsx = lx/lz;
    zsy = ly/lz;
    zshearmat(mmat,zsx,zsy);

/* rotate the hue */
    zrs = sin(rot*M_PI/180.0);
    zrc = cos(rot*M_PI/180.0);
    zrotatemat(mmat,zrs,zrc);

/* unshear the space to put the luminance plane back */
    zshearmat(mmat,-zsx,-zsy);

/* rotate the grey vector back into place */
    yrotatemat(mmat,-yrs,yrc);
    xrotatemat(mmat,-xrs,xrc);

    matrixmult(mmat,mat,mat);
}

main(argc,argv)
int argc;
char **argv;
{
    unsigned long *lbuf;
    int xsize, ysize;
    float mat[4][4];

    if(argc<2) {
	fprintf(stderr,"usage: matrix in.rgb\n");
	exit(1);
    }
    sizeofimage(argv[1],&xsize,&ysize);
    lbuf = (unsigned long*)longimagedata(argv[1]);
    prefsize(xsize,ysize);
    winopen("matrix color");
    RGBmode();
    gconfig();
    lrectwrite(0,0,xsize-1,ysize-1,lbuf);

    identmat(mat);

    offsetmat(mat,-50.0,-50.0,-50.0);	/* offset color */ 
    cscalemat(mat,1.4,1.5,1.6);		/* scale the colors */
    saturatemat(mat,1.7);		/* saturate by 2.0 */
    huerotatemat(mat,10.0);		/* rotate the hue 10 */
    printmat(mat);
    applymatrix(lbuf,mat,xsize*ysize);

    lrectwrite(0,0,xsize-1,ysize-1,lbuf);
    sleep(10);
    exit(0);
}