canny_d.c

canny 边缘检测算法

/* ************************************************************************ */
/* *************** SIMPLE SYMMETRIC CANNY EDGE OPERATOR ******************* */
/* ************************************************************************ */
/*
 *   PROGRAM NAME : canny_d.c  (modified from canny.c - GAWW - Nov 1992)
 *
 *   AUTHOR       : Andrew Zisserman (adapted from a Richard Tobin program)
 *
 *   MODIFIED BY  : Robert Lloyd (Heriot-Watt Uni, Dept. of Computer Sci)
 *
 *   MODIFIED ON  : 16th August 1986
 *
 *   MODIFIED BY  : G A W West, Beatrice Brillault - 6th March 1989
 *   MODIFIED BY  : G A W West 26th. Sept 1990, Perth, WA for use with
 *         command line arguments or defaults. Files must be
 *         inputted with arguments i and o on the command line
 *         otherwise the prog askes for them. Other parameters
 *         not on the command line are taken as the defaults.
 *   MODIFIED BY : G A W West 20-July-1992
 *         Works with any size images
 *
 *   MODIFIED BY : G A W West November 1992
 *         Uses Deriche's recursive filtering for first order stuff
 *
 *   MODIFIED BY : Paul Rosin November 1994
 *         Edge stability estimate (see Colchester etc.)
 *             but it has little effect!
 *         Can dump separate nms X & Y magnitude images
 *             useful for relaxation labelling using "crack edges"
 *
 *
 *   DESCRIPTION  : The program applies a simple canny edge operator to an
 *                  input file (image) up to a maximum of 512 pixels square.
 *                  Included are Simpson's rule and a correct treatment of
 *                  boundary.
 *                     Information is handled as a byte, with the direction
 *                  of the edge stored as a signed byte. For the edge
 *                  magnitude a zero entry represents no edge present at that
 *                  position.
 *
 *
 *   DEFAULTS     : Sigma of the Gaussian = 1
 *                  width = 512, height = 512
 */

#include <stdio.h>
#include <math.h>
#include "malloc_image.h"
#include "recursive_filter.h"

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

#define ABS(x)      (((x)<0.0) ? -(x): (x))

/* perform stability estimate using Laplacian */
#define LAPLACIAN FALSE
/* dump both X & Y edge magnitudes separately after nms */

#define PI 3.14159265

#define WIDTH 512
#define HEIGHT 512

unsigned char **in;        /* input grey image */
unsigned char **out;        /* output magnitude image */
unsigned char **dir;        /* output direction image */
float **x;
float **y;
float **rtemp;
float **g;

#if LAPLACIAN
float **gaussian;
float **laplacian;
float **stability;
#endif

#define FILTER_SIZE 20
static float m[FILTER_SIZE], m1[FILTER_SIZE];

double s = 1.0;            /* sigma for Gaussian */
int verbose = TRUE;
int both_output_files = FALSE;
char file_in[255],file_out_mag[255],file_out_dir[255],*file_out_pre;
int width = WIDTH, height = HEIGHT;
int border_width = 1;
int border_factor = 3;
char fname[255];
int dump_pre_nms = FALSE;

main(argc, argv)
int argc;
char *argv[];
{
    int n;
    char *temp, ch;
    int count, set_file_in, set_file_out_dir, set_file_out_mag;
    FILE *fp_in, *fp_out;
    int i, j;
    int dump_x_and_y = FALSE;

    if (argc > 1) {
        /* use command line options */
        /* -i file_in(string) -o file_out(string) -a (output both
           images) -s sigma (real) -w image_width (integer) -h
           image_height (integer) */
        count = 0;
        do {
            count++;
            temp = argv[count];
            if (*argv[count] == '-') {
                ch = *(++temp);
                switch (ch) {
                case 'o':
                    count++;
                    strcpy(file_out_mag, argv[count]);
                    set_file_out_mag = TRUE;
                    printf("read argument -o %s\n", file_out_mag);
                    break;
                case 'i':
                    count++;
                    strcpy(file_in, argv[count]);
                    set_file_in = TRUE;
                    printf("read argument -i %s\n", file_in);
                    break;
                case 'S':
                    dump_x_and_y = TRUE;
                    printf("dumping X & Y magnitudes in separate files\n");
                    break;
                case 's':
                    count++;
                    s = atof(argv[count]);
                    printf("read argument -s %f\n", s);
                    break;
                case 'w':
                    count++;
                    width = atoi(argv[count]);
                    break;
                case 'h':
                    count++;
                    height = atoi(argv[count]);
                    break;
                case 'a':
                    both_output_files = TRUE;
                    break;
                case 'p':
                    dump_pre_nms = TRUE;
                    count++;
                    file_out_pre = argv[count];
                    break;
                case 'd':
                    count++;
                    set_file_out_dir = TRUE;
                    strcpy(file_out_dir, argv[count]);
                    printf("read argument -d %s\n", file_out_dir);
                    break;
                default:
                    printf("error on command line (1)\n");
                }
            }
            else {
                printf("error on command line (2)\n");
                exit(-1);
            }
        } while (count < argc - 1);
        if (!set_file_out_mag) {
            printf("need output file name\n");
            exit(-1);
        }
        if (!set_file_in) {
            printf("need input file name\n");
            exit(-1);
        }
        if (both_output_files) {
            if ((!set_file_out_dir) || (!set_file_out_mag)) {
                printf("need both output files as used option a on command line\n");
                exit(-1);
            }
        }
        else {
            if ((set_file_out_dir) && (set_file_out_mag)) {
                printf("only need magnitude image as option a not specified\n");
                printf("ignoring dir image\n");
            }
        }
        printf("read arguments ok\n");
        printf("standard deviation: %f\n", s);
        printf("image size: %d wide %d high\n", width, height);
        printf("reading from file %s\n", file_in);
        printf("writing to file %s\n", file_out_mag);
        if (both_output_files)
            printf("writing direction image to %s\n", file_out_dir);

        /* check to see if input file available */ ;
        if ((fp_in = fopen(file_in, "r")) == NULL) {
            fprintf(stderr,"file %s not found - aborting\n", file_in);
            exit(-1);
        }
        else
            fclose(fp_in);

        /* check to see if can make the output magnitude file */ ;
        if ((fp_out = fopen(file_out_mag, "w")) == NULL) {
            fprintf(stderr,"file %s cannot be created - aborting\n", file_out_mag);
            exit(-1);
        }
        else
            fclose(fp_out);
        /* check to see if can make the output direction file */ ;
        if (both_output_files) {
            if ((fp_out = fopen(file_out_dir, "w")) == NULL) {
                fprintf(stderr,"file %s cannot be created - aborting\n", file_out_dir);
                exit(-1);
            }
            else
                fclose(fp_out);
        }

        /* allocate arrays */
        in = malloc_char_image(width, height);
        out = malloc_char_image(width, height);
        dir = malloc_char_image(width, height);
        x = malloc_float_image(width, height);
        y = malloc_float_image(width, height);
        rtemp = malloc_float_image(width, height);
        g = malloc_float_image(width, height);
#if LAPLACIAN
        gaussian = malloc_float_image(width, height);
        laplacian = malloc_float_image(width, height);
        stability = malloc_float_image(width, height);
#endif

        /* just used for non-maximal suppression here */
        n = calc_filter_coeffs();
        printf("n: %d\n", n);

        border_width = 3*s;
        printf("border width: %d\n", border_width);

        printf("reading image file\n");
        read_image(in, file_in, width, height);

        /* convert from char to float */
        for (j = 0; j < height; j++)
            for (i = 0; i < width; i++)
                rtemp[j][i] = (float) in[j][i];

        if (both_output_files)
            for (j = 0; j < height; j++)
                for (i = 0; i < width; i++)
                    dir[j][i] = 0;

        recursive_filter_horizontal(rtemp, x, width, height, s, GAUSS);
        recursive_filter_vertical(rtemp, y, width, height, s, GAUSS);

#if LAPLACIAN
        recursive_filter_vertical(x, gaussian, width, height, s, GAUSS);
        /* sum two outputs to get Gaussian smoothed image */
        for (j = 0; j < height; j++)
            for (i = 0; i < width; i++)
                gaussian[j][i] = hypot((double) x[j][i], (double) gaussian[j][i]);
        /* perform Laplacian of Gaussian */
        for (j = 0; j < height; j++)
            for (i = 0; i < width; i++)
                laplacian[j][i] = 0;

        for (j = 1; j < height-1; j++)
            for (i = 1; i < width-1; i++) {
                int ii, jj;

                for (jj = -1; jj <= 1; jj++)
                    for (ii = -1; ii <= 1; ii++)
                        laplacian[j][i] -= gaussian[j+jj][i+ii];
                laplacian[j][i] += gaussian[j][i] * 9;

                laplacian[j][i] /= 3.0;

                if (laplacian[j][i] < 0)
                    laplacian[j][i] = -laplacian[j][i];
            }

        /*
        rescale_image_3(laplacian, out, border_width+1);
        write_image(out, file_out_mag, width, height);
        exit(-1);
        */
#endif

        recursive_filter_vertical(x, x, width, height, s, FIRST_DERIV);
        recursive_filter_horizontal(y, y, width, height, s, FIRST_DERIV);

        /* sum two outputs to get g */
        for (j = 0; j < height; j++)
            for (i = 0; i < width; i++)
                g[j][i] = hypot((double) x[j][i], (double) y[j][i]);

        if (dump_pre_nms) {
            rescale_image_3(g, out, border_width+1);
            write_image(out, file_out_pre, width, height);
        }

#if LAPLACIAN
        /* calculate & incorporate stability */
        for (j = 0; j < height; j++)
            for (i = 0; i < width; i++) {
                stability[j][i] =
                    atan2((double) g[j][i], (double) laplacian[j][i]);
                /*
                */