/* * 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 50000int 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);