chamfer_sdt2.c

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

/*
 * computes salience distance transform
 * distance weighted by:
 *        edge magnitude
 *        pixel list length
 *
 * (serial) chamfer algorithm used, either 3-4 or 5-7-11
 *
 * Paul Rosin
 * August 1994
 */

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

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

#define K 1

#define MAX_SIZE 1024

#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 height,width,depth;

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

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 measure,min_measure,prev_measure;
    char *binary_file,*outfile,*edge_file,*pixel_file;
    double max_val;
    int min_loc;
    int xn,yn;
    FILE *fp;
    char file_type[50];
    int endoffile;

    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 'b':
                    i++;
                    binary_file = argv[i];
                    break;
                case 'e':
                    i++;
                    edge_file = argv[i];
                    break;
                case 'o':
                    i++;
                    outfile = argv[i];
                    break;
                case 'p':
                    i++;
                    pixel_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 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();
    }

    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

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

    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] = tmp[x][y];
    }
    else {
        /* threshold magnitude image */
        for (y = 0; y < height; y++)
            for (x = 0; x < width; x++)
                if ((unsigned int)image_mag[x][y] > 0)
                    image_dist[x][y] = 0;
                else
                    image_dist[x][y] = 255;
    }

    /* propagate distance etc */
    j = 0;
    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 int)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];

                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 int)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];

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

    /* calculate measure values */
    for (y = 1; y < height-1; y++) {
        for (x = 1; x < width-1; x++) {
            image_measure[x][y] = (double) (image_dist[x][y]+K) /
                                 ((double) ((unsigned int)image_mag[x][y]+K) *
                                            (image_length[x][y]+K));
        }
    }

    /* rescale */
    max_val = image_measure[1][1];
    for (y = 1; y < height-1; y++)
        for (x = 1; x < width-1; x++)
            if (image_measure[x][y] > max_val)
                max_val = image_measure[x][y];

    for (y = 1; y < height-1; y++)
        for (x = 1; x < width-1; x++)
            tmp[x][y] = (unsigned int)(image_measure[x][y] * 255 / max_val);

    /* reset borders */
    for (y = 0; y < height; y++)
        tmp[0][y] = tmp[width-1][y] = 255;
    for (x = 0; x < width; x++)
        tmp[x][0] = tmp[x][height-1] = 255;

    write_pgm(tmp,outfile,width,height);
}

store_length()
{
    int i;

    for (i = 0; i < no_pixels; i++)
        image_length[x[i]][y[i]] = no_pixels;
}

read_link_data(fp,endoffile)
FILE *fp;
int *endoffile;
{
    char dumstring[50];
    int j;

    fscanf(fp,"%s %d\n",dumstring,&j);
    j = -1;
    do{
       j++;
       fscanf(fp,"%d %d\n",&x[j],&y[j]);
    } while(x[j] != -1);
    *endoffile = (y[j] == -1);
    no_pixels = j;
}

options(progname)
char *progname;
{
    printf("usage: %s [options]\n",progname);
    printf("     -b file    binary edge map (optional)\n");
    printf("     -e file    edge magnitude map\n");
    printf("     -o file    output distance map\n");
    printf("     -p file    pixel lists\n");
    printf("     -n int     maximum number of iterations (default: %d)\n",max_iter);
    exit(-1);
}