From: ben@ph.tn.tudelft.nl (Ben Verwer)
Subject: Re: image processing loops: a c programming question
Organization: Delft University of Technology
Distribution: sci.image.processing
Date: Tue, 18 Aug 1992 17:03:57 GMT

The basic question is how to program for images with variable sizes. The 
solution presented by the author here (malloc a linear array of xsize*ysize) is
imo the best. Those numerical recipies in C readers who propose the famous 
array of pointers solution are thinking rather 2d. What if you go to 3d: an 
image of pointers? Which is of course implemented as an array of pointers to 
pointers. Good luck!
As for the author who proposed the function call in the innerloop: you must 
drink lots of coffee at work :)

The loop I most often use is:

pixel *p = beginofdata;
int i;  /* make sure this is 32 bits! */

for (i=xsize*ysize; --i>=0;) *p++ = ...;

This is the fastest on 68000Us! And also on sparcstations. Try it. (thanx to 
Gert van Antwerpen)
Make sure the optimizer is on. Else write register pixel *p and register int i;

What if you want to do something with *p++?

Of course you will only in the most exceptional of cases want to write a loop 
like in the original posting. Most filters can be written as a sequence of 2 
1-dimensional filters. One in x and one in y (which also scales nicely to 3d).
Such as maximum, minimum, uniform, gaussian filters.
But if you have to you should keep the innerloop simple. So the 
straightforward implementation:

for ysize
   for xsize
      for yfiltersize
          for xfiltersize


is hopeless. The innerloop is only a few iterations. And the work to be done 
there not really simple.

The following is a lot better:

for yfiltersize
    for xfiltersize
        for ysize
            for xsize


In that case the innerloop is the largest and simplest. Also handy if you have
a vector processing machine. Of course you need a little imaganation to 
construct the innermost part of the loop, with may include a temporary image or
two to hold intermediate results. Whereas in the most often used construct you
only need one or two ints. But memory is cheap, and my time is not. On 68000Us
and the like the speedup can be a factor of 3! On sparcUs however the speedup 
is negligable, because the cache is not well used. So the following is best:

for ysize
   for yfiltersize
      for xfiltersize
          for xsize


Ben Verwer. Below see a real program (coauthored by Henri Vrooman).



/* The rolling ball filter smooths an image by rolling a ball over the surface
of the image, (imagine grey-value as height) and tracing the contour of the 
ball as it goes */

static int ball[63*63]; /* rolling ball shape, maximum diameter 63 */

qball(inp,out,xsize,ysize,fsize)
short *inp,*out;
int xsize, ysize, fsize;
        {

        register int *q1, *q2;
        register int x,xend,v;

        int *input, *output;
        int i,j,y,b;
        int hf,len;

        char *malloc();
        double dtmp, sqrt();

        input = (int *)malloc(xsize*ysize*sizeof(int));
        output = (int *)malloc(xsize*ysize*sizeof(int));

        hf = fsize/2;

        len = xsize*ysize;
        for (j=0; j<len; j++) output[j] = input[j] = (int)inp[j];

        b = 0;
        for (j=-hf; j<=hf; j++) {
                for (i=-hf; i<=hf; i++) {

                        dtmp = .25*fsize*fsize-i*i-j*j;
                        if (dtmp < 0.0) ball[b++] = -1000; 
                        else ball[b++] = (int) sqrt(dtmp);
                }
        }

        xend = xsize-hf;

        for (y=hf; y<ysize-hf; y++) {

                b = 0;
                for (j=-hf; j<=hf; j++) {
                        for (i=-hf; i<=hf; i++) {

                                v = ball[b++];

                                q1 = input+y*xsize+hf;
                                q2 = output+y*xsize+hf+j*xsize+i;

                                for (x=0; x<xend; x++, q1++, q2++){
                                        if (*q2 < *q1+v) *q2 = *q1+v;
                                }
                        }
                }
        }


        len = xsize*ysize;
        for (j=0; j<len; j++) out[j] = output[j];

        free(output);
        free(input);

        return(0);
}

Newsgroups: sci.image.processing,rec.photo