Misc Notes‎ > ‎

Sum Tables/Integral Images for Convolution

Sum tables are an old trick for computing sums over areas in images quickly.  They were well-known in the statistics literature and in the 1980's started to be applied to convolutions of textures in the computer graphics community.   They were described in the Graphics Gems series of books.  More recently, they have become very popular in the computer vision community under the name "integral images".

I implemented them as part of the VIS library, in use at the MIT AI lab in the 1990's.  Here are two compact, self-contained implementations that I sent to the Gwydion project for benchmarking.

Sum tables turn out to be less useful than they seem at first glance; other methods are just as fast but give better results.

I found this code again by accident doing a search on Google Code.  The link is here:

http://www.cs.cmu.edu/afs/cs.cmu.edu/project/gwydion/hackers/clisp-hackers/ram/work/

C Version of Sum Table Code

/* Convolution using sum tables.
Thomas M. Breuel, 29 Jan 1992 */

#include <math.h>
#include <stdio.h>

void *malloc();

float **fialloc(w,h) {
float *base = malloc(w*h*sizeof (float));
float **result = malloc(w*sizeof (float*));
int i;
for(i=0;i<w;i++) result[i]=base+i*h;
return result;
}

fifree(p) float **p; {
free(p[0]);
free(p);
}

unsigned char **ucialloc(w,h) {
unsigned char *base = malloc(w*h*sizeof (unsigned char));
unsigned char **result = malloc(w*sizeof (unsigned char*));
int i;
for(i=0;i<w;i++) result[i]=base+i*h;
return result;
}

ucifree(p) unsigned char **p; {
free(p[0]);
free(p);
}

unsigned char **ucifread_raw(stream,w,h) FILE *stream; {
int i,j;
unsigned char **result=ucialloc(w,h);
for(i=0;i<w;i++) for(j=0;j<h;j++) result[i][j]=getc(stream);
return result;
}

ucifwrite_raw(stream,p,w,h) FILE *stream; unsigned char **p; {
int i,j;
for(i=0;i<w;i++) for(j=0;j<h;j++) putc(p[i][j],stream);
}

float **uci_to_fi(uci,w,h) unsigned char **uci; {
int i,j;
float **result=fialloc(w,h);
for(i=0;i<w;i++) for(j=0;j<h;j++) result[i][j]=uci[i][j];
return result;
}

unsigned char **fi_to_uci(fi,w,h) float **fi; {
int i,j;
unsigned char **result=ucialloc(w,h);
float min=1e38,max=-1e38;
float scale;
for(i=0;i<w;i++) for(j=0;j<h;j++) if(fi[i][j]<min) min=fi[i][j];
for(i=0;i<w;i++) for(j=0;j<h;j++) if(fi[i][j]>max) max=fi[i][j];
if(max<=min) scale=1.0; else scale=255.0/(max-min);
for(i=0;i<w;i++) for(j=0;j<h;j++) result[i][j]=(int)(scale*(fi[i][j]-min));
return result;
}

float **sumtable(fi,w,h) float **fi; {
int i,j; float sum;
float **result=fialloc(w,h);
for(i=0;i<w;i++) for(j=0,sum=0.0;j<h;j++) {
float t=fi[i][j];
result[i][j]=sum;
sum+=t;
}
for(j=0;j<h;j++) for(i=0,sum=0.0;i<w;i++) {
float t=result[i][j];
result[i][j]=sum;
sum+=t;
}
return result;
}

float sub(fi,w,h,x,y) float **fi; {
if(x<0) x=0; else if(x>=w) x=w-1;
if(y<0) y=0; else if(y>=h) y=h-1;
return fi[x][y];
}

float **locsum(fi,w,h,dx,dy) float **fi; {
int i,j;
float **sum=sumtable(fi,w,h);
float **result=fialloc(w,h);
#define S(x,y) sub(sum,w,h,x,y)
for(i=0;i<w;i++) for(j=0;j<h;j++)
result[i][j]=S(i-dx,j-dy)+S(i+dx,j+dy)-S(i-dx,j+dy)-S(i+dx,j-dy);
fifree(sum);
return result;
}

main() {
unsigned char **input,**output;
float **finput,**foutput;
fprintf(stderr,"reading ");
input = ucifread_raw(stdin,640,480);
fprintf(stderr,"converting ");
finput = uci_to_fi(input,640,480);
fprintf(stderr,"summing ");
foutput = locsum(finput,640,480,5,5);
fprintf(stderr,"converting ");
output = fi_to_uci(foutput,640,480);
fprintf(stderr,"writing ");
printf("VISF='n");
printf("ETYPE=uchar'n");
printf("DIMS=640 480'n");
printf("'n");
ucifwrite_raw(stdout,output,640,480);
fprintf(stderr,"done'n");
}

C++ Version of Sum Table Code

//
// Convolution using sum tables
//
// Thomas M. Breuel, 29 Jan 1992
//

extern "C" {
#include <math.h>
#include <stdio.h>
int abort();
}

class ByteImage;
class FloatImage;

struct ByteImage {
unsigned char *data;
int dims[2];
ByteImage(int w,int h) {
dims[0]=w; dims[1]=h;
data = new unsigned char[dims[0]*dims[1]];
}
ByteImage(int *dims) {
this->dims[0]=dims[0];
this->dims[1]=dims[1];
data = new unsigned char[dims[0]*dims[1]];
}
~ByteImage() {
delete data;
}
unsigned char &operator()(int i,int j) {
#ifndef UNSAFE
if(unsigned(i)>=dims[0]||unsigned(j)>=dims[1]) abort();
#endif
return data[i+j*dims[0]];
}
int dim(int n) {return dims[n];}
void operator=(ByteImage &other) {
if(data) delete data;
dims[0]=other.dims[0];
dims[1]=other.dims[1];
int n=dims[0]*dims[1];
data=new unsigned char[n];
for(int i=0;i<n;i++) data[i]=other.data[i];
}
ByteImage(ByteImage &other) {
fprintf(stderr,"[ByteImage copy]");
data=0;
*this=other;
}
operator FloatImage();
};

ByteImage read_raw(char *file,int w,int h) {
FILE *stream=fopen(file,"r");
ByteImage result(w,h);
for(int i=0;i<w;i++) for(int j=0;j<h;j++) {
char c = getc(stream);
result(i,j) = *(unsigned char *)&c;
}
fclose(stream);
return result;
}

void write_vis(char *file,ByteImage &image) {
FILE *stream=fopen(file,"w");
fprintf(stream,"VISF='n");
fprintf(stream,"DIMS=%d %d'n",image.dim(0),image.dim(1));
fprintf(stream,"ETYPE=uchar'n");
fprintf(stream,"'n");
for(int i=0;i<image.dim(0);i++) for(int j=0;j<image.dim(1);j++) {
unsigned char c=image(i,j);
putc(c,stream);
}
fclose(stream);
}

struct FloatImage {
float *data;
int dims[2];
FloatImage(int w,int h) {
dims[0]=w; dims[1]=h;
data = new float[dims[0]*dims[1]];
}
FloatImage(int *dims) {
this->dims[0]=dims[0];
this->dims[1]=dims[1];
data = new float[dims[0]*dims[1]];
}
~FloatImage() {
delete data;
}
float &operator()(int i,int j) {
#ifndef UNSAFE
if(unsigned(i)>=dims[0]||unsigned(j)>=dims[1]) abort();
#endif
return data[i+j*dims[0]];
}
int dim(int n) {return dims[n];}
void operator=(FloatImage &other) {
if(data) delete data;
dims[0]=other.dims[0];
dims[1]=other.dims[1];
int n=dims[0]*dims[1];
data=new float[n];
for(int i=0;i<n;i++) data[i]=other.data[i];
}
FloatImage(FloatImage &other) {
fprintf(stderr,"[FloatImage copy]");
data=0;
*this=other;
}
operator ByteImage();
};

ByteImage::operator FloatImage() {
int i,j;
FloatImage result(dims);
for(i=0;i<dim(0);i++) for(j=0;j<dim(1);j++) result(i,j) = float(operator()(i,j));
return result;
}

FloatImage::operator ByteImage() {
int i,j;
ByteImage result(dims);
float min=operator()(0,0);
for(i=0;i<dim(0);i++) for(j=0;j<dim(1);j++) if(operator()(i,j)<min) min=operator()(i,j);
float max=operator()(0,0);
for(i=0;i<dim(0);i++) for(j=0;j<dim(1);j++) if(operator()(i,j)>max) max=operator()(i,j);
float scale=255.0/(max-min);
for(i=0;i<dim(0);i++) for(j=0;j<dim(1);j++) result(i,j) = int(scale*(operator()(i,j)-min));
return result;
}

struct Sumtable:FloatImage {
Sumtable(int w,int h):FloatImage(w,h) {}
Sumtable(int *dims):FloatImage(dims) {}
float &operator()(int i,int j) {
if(i<0) i=0; else if(i>=dims[0]) i=dims[0]-1;
if(j<0) j=0; else if(j>=dims[1]) j=dims[1]-1;
return data[i+j*dim(0)];
}
void operator=(Sumtable &other) {
this->FloatImage::operator=(other);
}
Sumtable(Sumtable &other):FloatImage(other) {}
};

Sumtable sumtable(FloatImage &image) {
Sumtable result(image.dims);
int w=image.dim(0),h=image.dim(1);
for(int i=0;i<w;i++) {
float sum=0.0;
for(int j=0;j<h;j++) {
float t=image(i,j);
result(i,j)=sum;
sum+=t;
}
}
for(int j=0;j<h;j++) {
float sum=0.0;
for(int i=0;i<w;i++) {
float t=result(i,j);
result(i,j)=sum;
sum+=t;
}
}
return result;
}

FloatImage locsum(FloatImage &image,int dx,int dy) {
Sumtable S = sumtable(image);
FloatImage result(image.dims);
int i,j;
int w=image.dim(0),h=image.dim(1);
for(i=0;i<w;i++) for(j=0;j<h;j++)
result(i,j)=S(i-dx,j-dy)+S(i+dx,j+dy)-S(i-dx,j+dy)-S(i+dx,j-dy);
return result;
}

main() {
fprintf(stderr,"reading ");
ByteImage input = read_raw("/vmunix",640,480);
fprintf(stderr,"converting ");
FloatImage finput = input;
fprintf(stderr,"summing ");
FloatImage foutput = locsum(finput,5,5);
fprintf(stderr,"converting ");
ByteImage output = foutput;
fprintf(stderr,"writing ");
write_vis("/tmp/1",output);
fprintf(stderr,"done'n");
}

Comments