chamfer_sdt3.c

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

        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();
        combine_clutter();
    } while (!endoffile);
    fclose(fp);

    /* ++++++++++++++++++ SECOND PASS ++++++++++++++++++ */

    /* 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) *
                                          (clutter_area[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];
                    c[i] = clutter_area[xn][yn];
                }

                min_loc = 0;
                min_measure = (double) (d[0]+K) /
                              (double) ((m[0]+K) * (l[0]+K) * (c[0]+K));
                for (i = 1; i < MASK_SIZE; i++) {
                    measure = (double) (d[i]+K) /
                              (double) ((m[i]+K) * (l[i]+K) * (c[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];
                clutter_area[x][y] = c[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) *
                                          (clutter_area[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];
                    c[i] = clutter_area[xn][yn];
                }

                min_loc = 0;
                min_measure = (double) (d[0]+K) /
                              (double) ((m[0]+K) * (l[0]+K) * (c[0]+K));
                for (i = 1; i < MASK_SIZE; i++) {
                    measure = (double) (d[i]+K) /
                              (double) ((m[i]+K) * (l[i]+K) * (c[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];
                clutter_area[x][y] = c[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) *
                                            (clutter_area[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 char)(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;

    if (cancel_length)
        /* arbitrary constant value */
        for (i = 0; i < no_pixels; i++)
            image_length[x[i]][y[i]] = 100;
    else
        for (i = 0; i < no_pixels; i++)
            image_length[x[i]][y[i]] = no_pixels;
}

store_curve_labels(label)
{
    int i;

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

combine_clutter()
{
    int i;
    float total = 0;

    for (i = 0; i < no_pixels; i++)
        total += clutter_area[x[i]][y[i]];

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

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;
}

/* save coordinates of nearest edge point to each point in text file */
save_ascii(filename)
char filename[];
{
    int x,y;
    FILE *fp;
    double dx,dy;

    if ((fp=fopen(filename,"w")) == NULL) {
        printf("cant open file: %s\n",filename);
        exit(-1);
    }

    for (y = 1; y < height-1; y++)
        for (x = 1; x < width-1; x++) {
            dx = x - X[x][y];
            dy = y - Y[x][y];
            fprintf(fp,"%d %d ",x,y);
            fprintf(fp,"%d %d\n",X[x][y],Y[x][y]);
        }

    fclose(fp);
}

/* save coordinates of nearest edge point to each point in 2 images */
save_xy(fnx,fny)
char fnx[],fny[];
{
    int x,y;

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

    for (y = 1; y < height-1; y++)
        for (x = 1; x < width-1; x++)
            tmp[x][y] = X[x][y];
    write_pgm(tmp,fnx,width,height);

    for (y = 1; y < height-1; y++)
        for (x = 1; x < width-1; x++)
            tmp[x][y] = Y[x][y];
    write_pgm(tmp,fny,width,height);
}

/* save image of Voronoi cell labels from each point */
save_voronoi1(filename)
char filename[];
{
    static int label[MAX_SIZE][MAX_SIZE];
    static int hist[MAX_SIZE*MAX_SIZE][2];
    int x,y;
    int i;
    int count;

    /* create unique label for each Voronoi cell */
    for (y = 1; y < height-1; y++)
        for (x = 1; x < width-1; x++)
            label[x][y] = X[x][y] + Y[x][y] * MAX_SIZE;

    for (i = 0; i < MAX_SIZE*MAX_SIZE; i++)
        hist[i][0] = hist[i][1] = 0;

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

    count = 0;
    for (i = 0; i < MAX_SIZE*MAX_SIZE; i++) {
        if (hist[i][0] > 0) {
            hist[i][1] = count;
            count++;
        }
    }

    printf("number of Voronoi cells: %d\n",count);

    for (x = 0; x < width; x++)
        for (y = 0; y < height; y++)
            tmp[x][y] = (unsigned char) hist[label[x][y]][1];

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

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

/* save image of Voronoi cell labels from each curve */
save_voronoi2(pixel_file,voronoi_file)
char *pixel_file,*voronoi_file;
{
    int x,y;
    int i,j;
    int count;
    FILE *fp;
    char file_type[50];
    int endoffile;

    /* 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++)
            curve_labels[j][i] = 0;

    count = 1;
    do {
        read_link_data(fp,&endoffile);
        store_curve_labels(count);
        count++;
    } while (!endoffile);
    fclose(fp);

    if (count > 255) {
        fprintf(stderr,"ERROR: too many pixel lists (%d) to store as Voronoi labels\n",count);
        return;
    }

    /* store curve label for each Voronoi cell */
    for (y = 1; y < height-1; y++)
        for (x = 1; x < width-1; x++)
            tmp[x][y] = (unsigned char) curve_labels[X[x][y]][Y[x][y]];

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

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

/* save image of clutter values at each curve */
save_clutter(pixel_file,clutter_file)
char *pixel_file,*clutter_file;
{
    int xx,yy;
    int i,j;
    FILE *fp;
    char file_type[50];
    int endoffile;
    double avg,max_val;

    /* 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 (yy = 0; yy < height; yy++)
       for (xx = 0; xx < width; xx++)
          tmp2[xx][yy] = 0;

    max_val = 0;
    do {
        read_link_data(fp,&endoffile);

        avg = 0;
        for (i = 0; i < no_pixels; i++)
            avg += clutter_area[x[i]][y[i]];
        avg /= (double)no_pixels;

        for (i = 0; i < no_pixels; i++)
            tmp2[x[i]][y[i]] = avg;

        if (avg >  max_val)
            max_val = avg;
    } while (!endoffile);
    fclose(fp);

    printf("rescaling from maximum clutter value: %f\n",max_val);
    for (yy = 0; yy < height; yy++)
       for (xx = 0; xx < width; xx++)
          tmp[xx][yy] = (unsigned char)(tmp2[xx][yy] * 255.0 / max_val);

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

options(progname)
char *progname;
{
    printf("usage: %s [options]\n",progname);
    printf("     -a file    save coords of nearest edge as ASCII file\n");
    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);
    printf("     -v file    save Voronoi diagram from 1st pass\n");
    printf("     -c file    save edge clutter map from 1st pass\n");
    printf("     -L         do NOT use length attribute\n");
    exit(-1);
}