/* * link_dynamic.c * * extracts lists of edge pixels from an image * automatically determines a threshold value to reject poor edge lists * * two threshold methods: *     1/ linear decision rule *     2/ 3*sigma/length decision rule *  * usage: run the program twice *        1st time: edge list statistics are generated *        then calculate median information, using: *             1/ median4 or *             2/ estimate_noise *        2nd time: edge lists above threshold are output * it's much easier to use the shell scripts provided! * * original program (based on "linknew.c" written by Geoff West): * * partly described in: * Svetha Venkatesh & Paul L. Rosin, * Dynamic threshold determination by local and global edge evaluation, * Graphical Models & Image Processing, * Vol. 75, No. 2, pp. 146-160, 1995. * * Paul Rosin & Svetha Venkatesh * Curtin University * Perth, Western Australia * February 1993 * * ------------------------------------------------------------------ * * updated with non-linear decision rule: * Paul Rosin * Joint Research Centre * Ispra, Italy * December 1994 * paul.rosin@jrc.it */#include <stdio.h>#include <stdlib.h>#include <math.h>#include "malloc_image.h"#ifndef FALSE# define FALSE 0# define TRUE (!FALSE)#endif#define MAX_PIX    10000#define BLACK      0#define WHITE      255#define WIDTH      512#define HEIGHT     512#define OLD_INFO   0    /* output all edge list information */#define OLD_THRESH 1    /* output edge lists above threshold */#define NEW_INFO   2    /* output all edge list information - new format */#define NEW_THRESH 3    /* output edge lists above threshold - new version */int   list_no;float xpix[MAX_PIX],ypix[MAX_PIX];int   Index;long  weight;int   flag;int   no_lists_written;float aspect_ratio;    /* defined as y/x and divide y value by this */int   closed_only;int   floating_point;FILE  *fp_out;unsigned char  **image;int   height,width;char  *file_in=NULL,*file_out=NULL;int   round();FILE *fp_in;/* additional variables added for automatic thresholding */float xx1,yy1,xx2,yy2;float slope,intercept;/* new additional variables added for automatic thresholding */float med_mag,dev;float fact = 2;main(argc,argv)int  argc;char *argv[];{    int i;    int option = OLD_INFO;    closed_only = FALSE;    aspect_ratio = 1.0;    floating_point = FALSE;    for (i = 1; i < argc; i++) {        if (argv[i][0] == '-') {            switch(argv[i][1]) {                case 'i':                    i++;                    file_in = argv[i];                    break;                case 'o':                    i++;                    file_out = argv[i];                    break;                case 'n':                    option = NEW_INFO;                    break;                case 'c':                    closed_only = TRUE;                    break;                case 'a':                    i++;                    aspect_ratio = atof(argv[i]);                    printf("aspect ratio is: %f\n",aspect_ratio);                    break;                case 'd':                    i++;                    xx1 = atof(argv[i]);                    i++;                    yy1 = atof(argv[i]);                    i++;                    xx2 = atof(argv[i]);                    i++;                    yy2 = atof(argv[i]);                    option = OLD_THRESH;                    slope = (float)(yy2 - yy1) / (float)(xx2 -xx1);                    intercept = yy1;                    break;                case 'D':                    i++;                    med_mag = atof(argv[i]);                    i++;                    dev = atof(argv[i]);                    option = NEW_THRESH;                    break;                case 'F':                    i++;                    fact = atof(argv[i]);                    break;                case 'f':                    floating_point = TRUE;                    break;                default:                    fprintf(stderr,"unknown option %s\n",argv[i]);                    options(argv[0]);            }        }        else {            fprintf(stderr,"unknown option %s\n",argv[i]);            options(argv[0]);        }    }    if ((file_in == NULL) || (file_out == NULL)) {        fprintf(stderr,"ERROR: need filenames\n");        options(argv[0]);    }    link_and_save(option);}link_and_save(option)int option;{    flag = FALSE;        /* at start only - no lists */    no_lists_written = 0;    read_pgm_header(file_in,&width,&height);    image = malloc_char_image(width,height);    read_image_body(image,width,height);    /* clean up image, make all object 8 connected */    clean();    /* remove isolated points - of no interest to us here */    remove_isolated();    if ((fp_out = fopen(file_out, "w")) == NULL) {        printf("file %s cannot be created - aborting\n", file_out);        exit(-1);    }    /* put magic name at top of file */    if (option == OLD_INFO) {        fprintf(fp_out,"super\nlist:\n");    }    else if ((option == OLD_THRESH) || (option == NEW_THRESH)) {        if (floating_point == TRUE)            fprintf(fp_out, "pixel_float\n");        else            fprintf(fp_out, "pixel\n");    }    /* link open edges */    list_no = 0;    link_open(option);    /* remove isolated points - of no interest to us here */    remove_isolated();    /* link closed edges */    link_closed(option);    /* remove isolated points - of no interest to us here */    remove_isolated();    if (option == OLD_INFO)        fprintf(fp_out,"endl:\nendf:\n");    else if ((option == OLD_THRESH) || (option == NEW_THRESH))        fprintf(fp_out, "-1 -1\n");    fclose(fp_out);}clean(){    int             loop1, loop2;    unsigned char   i1, i2, i3, i4;    /* clear border */    for (loop1 = 0; loop1 < height; loop1++)        image[loop1][0] = image[loop1][width - 1] = BLACK;    for (loop1 = 0; loop1 < width; loop1++)        image[0][loop1] = image[height - 1][loop1] = BLACK;    for (loop1 = 1; loop1 < height - 1; loop1++)        for (loop2 = 1; loop2 < width - 1; loop2++)            if (image[loop1][loop2] != BLACK) {                i1 = image[loop1 - 1][loop2];                i2 = image[loop1][loop2 - 1];                i3 = image[loop1 + 1][loop2];                i4 = image[loop1][loop2 + 1];                if ((i1 != BLACK) && (i2 != BLACK))                    image[loop1][loop2] = BLACK;                else if ((i2 != BLACK) && (i3 != BLACK))                    image[loop1][loop2] = BLACK;                else if ((i3 != BLACK) && (i4 != BLACK))                    image[loop1][loop2] = BLACK;                else if ((i4 != BLACK) && (i1 != BLACK))                    image[loop1][loop2] = BLACK;            }}remove_isolated(){    int loop1,loop2;    unsigned char i1,i2,i3,i4,i6,i7,i8,i9;    for (loop1 = 1; loop1 < height - 1; loop1++)        for (loop2 = 1; loop2 < width - 1; loop2++)            if (image[loop1][loop2] != BLACK) {                i1 = image[loop1 - 1][loop2 - 1];                i2 = image[loop1][loop2 - 1];                i3 = image[loop1 + 1][loop2 - 1];                i4 = image[loop1 - 1][loop2];                i6 = image[loop1 + 1][loop2];                i7 = image[loop1 - 1][loop2 + 1];                i8 = image[loop1][loop2 + 1];                i9 = image[loop1 + 1][loop2 + 1];                if ((i1 + i2 + i3 + i4 + i6 + i7 + i8 + i9) == (8 * BLACK)) {                    image[loop1][loop2] = BLACK;                }            }}link_open(option)int option;{    int             loop1,loop2,loop3;    unsigned char   i1,i2,i3,i4,i6,i7,i8,i9;    int             xp,yp;    int             end_of_line;    int ok;    for (loop1 = 0; loop1 < height; loop1++)         /* for each row y */        for (loop2 = 0; loop2 < width; loop2++) {    /* for each column x */            /* find pixel at end of line */            if (image[loop1][loop2] != BLACK) {                i1 = i2 = i3 = i4 = i6 = i7 = i8 = i9 = 0;                if (image[loop1 - 1][loop2 - 1] != BLACK)                    i1 = 1;                if (image[loop1][loop2 - 1] != BLACK)                    i2 = 1;                if (image[loop1 + 1][loop2 - 1] != BLACK)                    i3 = 1;                if (image[loop1 - 1][loop2] != BLACK)                    i4 = 1;                if (image[loop1 + 1][loop2] != BLACK)