chamfer_sdt3.c

The salience distance transform incorporates edge strength information into the distance transform.

/*
 * computes salience distance transform
 * can also generate Voronoi diagram
 * distance weighted by:
 *        edge magnitude
 *        pixel list length
 *        clutterness (averaged over each pixel list)
 *
 * (serial) chamfer algorithm used, either 3-4 or 5-7-11
 *
 * points not in pixel lists (e.g. junctions) are removed from
 * gradient image
 *
 * Paul Rosin
 * March 1995
 */

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

#ifndef FALSE
# define FALSE 0
# define TRUE (!FALSE)
#endif

#define SIMPLE TRUE /* quick hack to get XY images output without saliency stuff*/

#define K 1

#define MAX_SIZE 1500

#include "pgmio.h"

#define CHAMFER34   0
#define CHAMFER5711 1

#define ALGORITHM CHAMFER34

#if ALGORITHM == CHAMFER34
# define MASK_SIZE 5
  int x_offset[MASK_SIZE] = { -1,  0,  1, -1,  0 };
  int y_offset[MASK_SIZE] = { -1, -1, -1,  0,  0 };
  int i_offset[MASK_SIZE] = {  4,  3,  4,  3,  0 };
#elif ALGORITHM == CHAMFER5711
# define MASK_SIZE 9
  int x_offset[MASK_SIZE] = { -1,  1, -2, -1,  0,  1,  2, -1,  0 };
  int y_offset[MASK_SIZE] = { -2, -2, -1, -1, -1, -1, -1,  0,  0 };
  int i_offset[MASK_SIZE] = { 11, 11, 11,  7,  5,  7, 11,  5,  0 };
#endif

#define MAX_PIXELS 50000
int x[MAX_PIXELS],y[MAX_PIXELS];
int no_pixels;

unsigned char tmp[MAX_SIZE][MAX_SIZE];

int image_dist[MAX_SIZE][MAX_SIZE];
unsigned char image_mag[MAX_SIZE][MAX_SIZE];
double image_measure[MAX_SIZE][MAX_SIZE];
int image_length[MAX_SIZE][MAX_SIZE];

int X[MAX_SIZE][MAX_SIZE];
int Y[MAX_SIZE][MAX_SIZE];

double clutter_area[MAX_SIZE][MAX_SIZE];
double tmp2[MAX_SIZE][MAX_SIZE];

int curve_labels[MAX_SIZE][MAX_SIZE];

int height,width,depth;

int max_iter = 10;    /* usually terminates after 4 or 5 iterations anyway */

int cancel_length = FALSE;

char fnx[] = "xxx";
char fny[] = "yyy";

FILE *fp_in;

main(argc,argv)
int argc;
char *argv[];
{
    int i,j,x,y;
    int change;
    int d[MASK_SIZE],m[MASK_SIZE],l[MASK_SIZE];
    double c[MASK_SIZE];
    double measure,min_measure,prev_measure;
    char *binary_file,*outfile,*edge_file,*pixel_file;
    char *voronoi_file,*clutter_file,*ascii_file;
    double max_val;
    int min_loc;
    int xn,yn;
    FILE *fp;
    char file_type[50];
    int endoffile;

    ascii_file = clutter_file = voronoi_file = NULL;
    pixel_file = binary_file = outfile = edge_file = NULL;

    /* parse command line */
    for (i = 1; i < argc; i++) {
        if (argv[i][0] == '-') {
            switch(argv[i][1]) {
                case 'a':
                    i++;
                    ascii_file = argv[i];
                    break;
                case 'b':
                    i++;
                    binary_file = argv[i];
                    break;
                case 'e':
                    i++;
                    edge_file = argv[i];
                    break;
                case 'o':
                    i++;
                    outfile = argv[i];
                    break;
                case 'L':
                    cancel_length = TRUE;
                    break;
                case 'p':
                    i++;
                    pixel_file = argv[i];
                    break;
                case 'c':
                    i++;
                    clutter_file = argv[i];
                    break;
                case 'v':
                    i++;
                    voronoi_file = argv[i];
                    break;
                case 'n':
                    i++;
                    max_iter = atoi(argv[i]);
                    printf("number of iterations: %d\n",max_iter);
                    break;
                default:
                    printf("unknown option %s\n",argv[i]);
                    options(argv[0]);
            }
        }
        else {
            printf("unknown option %s\n",argv[i]);
            options(argv[0]);
        }
    }

    if (outfile == NULL || edge_file == NULL || pixel_file == NULL)
        options(argv[0]);

    read_pgm(image_mag,edge_file,&width,&height,&depth);

    if ((width > MAX_SIZE) || (height > MAX_SIZE)) {
        fprintf(stderr,"ERROR: Maximum image size is %d x %d\n",
            MAX_SIZE,MAX_SIZE);
        exit(-1);
    }

    if (cancel_length)
        printf("NOT using length attribute\n");

    /* read and pixel lists and store lengths */
    if ((fp=fopen(pixel_file,"r")) == NULL) {
        printf("cant open pixel file: %s\n",pixel_file);
        exit(-1);
    }
    /* read magic word for format of file */
    fscanf(fp,"%s\n",file_type);
    j = strcmp(file_type,"pixel");
    if (j != 0){
        printf("not link data file - aborting\n");
        exit(-1);
    }

    for (j=0;j<height;j++)
        for (i=0;i<width;i++)
            image_length[j][i] = 0;

    do {
        read_link_data(fp,&endoffile);
        store_length();
    } while (!endoffile);
    fclose(fp);

#if (ALGORITHM == CHAMFER34)
    printf("performing chamfer3-4\n");
#elif (ALGORITHM == CHAMFER5711)
    printf("performing chamfer5-7-11\n");
#endif

    /* remove junctions in gradient image not in pixel lists */
    for (y = 0; y < height; y++)
        for (x = 0; x < width; x++)
            if (image_length[x][y] == 0)
                image_mag[x][y] = 0;

    if (binary_file != NULL) {
        read_pgm(tmp,binary_file,&width,&height,&depth);
        /* copy into integer array */
        for (y = 0; y < height; y++)
            for (x = 0; x < width; x++)
                image_dist[x][y] = (int)(tmp[x][y]) * 10;
    }
    else {
        /* threshold magnitude image */
        for (y = 0; y < height; y++)
            for (x = 0; x < width; x++)
                if ((unsigned char)image_mag[x][y] != 0)
                    image_dist[x][y] = 0;
                else
                    image_dist[x][y] = 255 * 10;
    }

    /* +++++++++++++ FIRST PASS - calculate Voronoi diagram +++++++++++++ */

    printf("calculating Voronoi diagram\n");

    /* initialise measure arrays */
    for (y = 0; y < height; y++)
        for (x = 0; x < width; x++) {
            X[x][y] = x;
            Y[x][y] = y;
        }

    /* propagate distance etc */
    j = 0;
#if SIMPLE
	{
        /* forward pass */
        for (y = 1; y < height-1; y++) {
            for (x = 1; x < width-1; x++) {
                for (i = 0; i < MASK_SIZE; i++) {
                    xn = x + x_offset[i];
                    yn = y + y_offset[i];
                    d[i] = image_dist[xn][yn] + i_offset[i];
                }

                min_loc = 0;
                min_measure = d[0];
                for (i = 1; i < MASK_SIZE; i++) {
                    measure = d[i];
                    if (measure < min_measure) {
                        min_measure = measure;
                        min_loc = i;
                    }
                }
                image_dist[x][y] = d[min_loc];
                xn = x + x_offset[min_loc];
                yn = y + y_offset[min_loc];
                X[x][y] = X[xn][yn];
                Y[x][y] = Y[xn][yn];
            }
        }

        /* backward pass */
        for (y = height-2; y >= 1; y--) {
            for (x = width-2; x >= 1; x--) {
                for (i = 0; i < MASK_SIZE; i++) {
                    xn = x - x_offset[i];
                    yn = y - y_offset[i];
                    d[i] = image_dist[xn][yn] + i_offset[i];
                }

                min_loc = 0;
                min_measure = d[0];
                for (i = 1; i < MASK_SIZE; i++) {
                    measure = d[i];
                    if (measure < min_measure) {
                        min_measure = measure;
                        min_loc = i;
                    }
                }
                image_dist[x][y] = d[min_loc];
                xn = x - x_offset[min_loc];
                yn = y - y_offset[min_loc];
                X[x][y] = X[xn][yn];
                Y[x][y] = Y[xn][yn];
            }
        }
	}
    save_xy(fnx,fny);
    fprintf(stderr,"Generated XY images; aborting now\n");
    exit(-1);
#else
    do {
        j++;
        printf("iteration %d\n",j);
        change = FALSE;

        /* forward pass */
        for (y = 1; y < height-1; y++) {
            for (x = 1; x < width-1; x++) {
                prev_measure = (double) (image_dist[x][y]+K) /
                               (double) (((unsigned int)image_mag[x][y]+K) *
                                          (image_length[x][y]+K));
                for (i = 0; i < MASK_SIZE; i++) {
                    xn = x + x_offset[i];
                    yn = y + y_offset[i];
                    d[i] = image_dist[xn][yn] + i_offset[i];
                    m[i] = (unsigned char)image_mag[xn][yn];
                    l[i] = image_length[xn][yn];
                }

                min_loc = 0;
                min_measure = (double) (d[0]+K) /
                              (double) ((m[0]+K) * (l[0]+K));
                for (i = 1; i < MASK_SIZE; i++) {
                    measure = (double) (d[i]+K) / (double) ((m[i]+K) * (l[i]+K));
                    if (measure < min_measure) {
                        min_measure = measure;
                        min_loc = i;
                    }
                }
                image_dist[x][y] = d[min_loc];
                image_mag[x][y] = m[min_loc];
                image_length[x][y] = l[min_loc];
                xn = x + x_offset[min_loc];
                yn = y + y_offset[min_loc];
                X[x][y] = X[xn][yn];
                Y[x][y] = Y[xn][yn];

                if (prev_measure != min_measure)
                    change = TRUE;
            }
        }

        /* backward pass */
        for (y = height-2; y >= 1; y--) {
            for (x = width-2; x >= 1; x--) {
                prev_measure = (double) (image_dist[x][y]+K) /
                               (double) (((unsigned int)image_mag[x][y]+K) *
                                          (image_length[x][y]+K));
                for (i = 0; i < MASK_SIZE; i++) {
                    xn = x - x_offset[i];
                    yn = y - y_offset[i];
                    d[i] = image_dist[xn][yn] + i_offset[i];
                    m[i] = (unsigned char)image_mag[xn][yn];
                    l[i] = image_length[xn][yn];
                }

                min_loc = 0;
                min_measure = (double) (d[0]+K) /
                              (double) ((m[0]+K) * (l[0]+K));
                for (i = 1; i < MASK_SIZE; i++) {
                    measure = (double) (d[i]+K) / (double) ((m[i]+K) * (l[i]+K));
                    if (measure < min_measure) {
                        min_measure = measure;
                        min_loc = i;
                    }
                }
                image_dist[x][y] = d[min_loc];
                image_mag[x][y] = m[min_loc];
                image_length[x][y] = l[min_loc];
                xn = x - x_offset[min_loc];
                yn = y - y_offset[min_loc];
                X[x][y] = X[xn][yn];
                Y[x][y] = Y[xn][yn];

                if (prev_measure != min_measure)
                    change = TRUE;
            }
        }
    } while (change && j < max_iter);
#endif

    if (voronoi_file != NULL) {
        save_voronoi2(pixel_file,voronoi_file);
    }

    if (ascii_file != NULL)
        save_ascii(ascii_file);

    /* ++++++++++++++++++ INTERMEDIATE PROCESSING ++++++++++++++++++ */

    for (y = 0; y < height; y++)
        for (x = 0; x < width; x++)
            clutter_area[x][y] = 0;

    for (y = 1; y < height-1; y++)
        for (x = 1; x < width-1; x++) {
            int xx,yy;

            xx = X[x][y];
            yy = Y[x][y];
            clutter_area[xx][yy]++;
        }

    if (clutter_file != NULL) {
        save_clutter(pixel_file,clutter_file);
    }

    /* read and pixel lists and store lengths */
    if ((fp=fopen(pixel_file,"r")) == NULL) {
        printf("cant open %s\n",pixel_file);