canny_d.c

canny 边缘检测算法

                g[j][i] *= stability[j][i];
            }
        /*
        rescale_image_3(g, out, border_width+1);
        write_image(out, file_out_mag, width, height);
        exit(-1);
        */
#endif

        if (both_output_files)
            calc_direction(n);

        printf("performing non maximal suppression\n");
        non_maximum_suppression(n);

        if (dump_x_and_y) {
            for (j = 0; j < height; j++)
                for (i = 0; i < width; i++) {
                    if (rtemp[j][i] == 0)
                        x[j][i] = y[j][i] = 0;
                    else {
                        x[j][i] = ABS(x[j][i]);
                        y[j][i] = ABS(y[j][i]);
                    }
                }
            rescale_image_1(x, out);
            sprintf(fname,"%s.X",file_out_mag);
            write_image(out, fname, width, height);
            rescale_image_1(y, out);
            sprintf(fname,"%s.Y",file_out_mag);
            write_image(out, fname, width, height);
        }

        printf("rescaling magnitude image\n");
        rescale_image_1(rtemp, out);
        /*** not necessary
        if (both_output_files) {
            printf("rescaling direction image\n");
            rescale_image_2(dir);
        }
        ***/

        write_image(out, file_out_mag, width, height);
        if (both_output_files) {
            write_image(dir, file_out_dir, width, height);
        }

    }
    else {
        printf("usage:  %s -i file_in -o file_mag -s st.dev\n",argv[0]);
        printf("options:\n");
        printf(" -S dump X & Y edge magnitudes in separate files\n");
        printf(" -s val for sigma of gaussian (real)\n");
        printf(" -a output magnitude and direction images\n");
        printf(" -i input image name\n");
        printf(" -o output magnitude image name\n");
        printf(" -d output direction image name\n");
        printf(" -p output magnitude image name (pre nms)\n");
        printf(" -w val for image width\n");
        printf(" -h val for image height\n");
        exit(-1);
    }
}



/* ************************************************************************ */
/* * CALC_FILTER_COEFFS *                                                   */
/* **********************                                                   */
/* The gaussian is a two-dimensional smoothing filter applied to the input  */
/* image to reduce the effects of noise. The weightings for each discrete   */
/* position along the curve are calculated and stored here.                 */
/* ************************************************************************ */
calc_filter_coeffs()
{
    int n;
    double a, b, c, gauss_const, gauss_expon();

    gauss_const = 1.0 / (6.0 * (s * s) * (2 * PI));    /* 6 for simpsons rule */

    for (n = 0; n < FILTER_SIZE; ++n)
        if (((a = gauss_expon((double) n, s)) > 0.001) || (n < 2)) {
            b = gauss_expon((n - 0.5), s);
            c = gauss_expon((n + 0.5), s);
            m[n] = (gauss_const * (b + (4 * a) + c));    /* simpsons rule */
            m1[n] = (c - b);
        }
        else {
            --n;
            break;
        };

    return (n);
}

/* Calculate edge orientation in range (-Pi, Pi) and store in [0,255] */
calc_direction(n)
int n;
{
    int i, j, npl1;
    double angle;
    double ux, uy;

    npl1 = (n + 1);

    for (i = npl1; i < (height - npl1); ++i) {
        for (j = npl1; j < (width - npl1); ++j) { 
            ux = x[i][j];
            uy = y[i][j];

            /* ********************************************************** */
            /* pi*40 ~= 128  - direction is (thought of as) a signed byte */
            /* ********************************************************** */
            /* dir[i][j] = (atan2( uy, ux) * 40); changed to sobel */
            angle = atan2(uy, ux); /* -Pi -> Pi */
            angle += PI; 
            angle = angle * 360.0 / (2*PI);    /* radians to degrees */
            dir[i][j] = (unsigned char) (angle * 255.0 / 360.0);
        }
    }
}

/* ************************************************************************ */
/* * NON_MAXIMUM_SUPPRESSION *                                              */
/* ***************************                                              */
/* To fix the position of the maximum edge strength and direction, maximum  */
/* suppression is applied. Each pixel in turn forms the centre of a nine    */
/* pixel neighbourhood. By interpolation, the strengths are calculated at   */
/* the boundary in both directions perpendicular to the centre pixel. The   */
/* centre pixel is marked as the maximum strength if it is greater than     */
/* the two values at the neighbourhood boundary.                            */
/* ************************************************************************ */
non_maximum_suppression(n)
int n;
{
    int i, j, npl1;
    int continue_flag = FALSE;
    double ux, uy, g0;

    /* zero temp image */
    for (j=0;j<height;j++)
        for (i=0;i<width;i++)
            rtemp[j][i] = 0;
    
    npl1 = (n + 1);

    for (i = npl1; i < (height - npl1); ++i) {
        for (j = npl1; j < (width - npl1); ++j) {    /* n-m s */
            ux = x[i][j];
            uy = y[i][j];
            if ((ux * uy) > 0)
                continue_flag = suppression_1(i, j, ux, uy, &g0);
            else
                continue_flag = suppression_2(i, j, ux, uy, &g0);

            if (!continue_flag) {
                rtemp[i][j] = g[i][j];
            }
        }
    }
}

rescale_image_1(rimage, image)
float **rimage;
unsigned char **image;
{
    int i, j;
    float max_conv, min_conv;
    float scale_factor;

    /* rescale image to be from 1 to 254 */
    max_conv = rimage[0][0];
    min_conv = rimage[0][0];
    for (i = 0; i < height; i++)
        for (j = 0; j < width; j++)
            if (rimage[i][j] > max_conv) {
                max_conv = rimage[i][j];
            }
            else if (rimage[i][j] < min_conv) {
                min_conv = rimage[i][j];
            }
    printf("image range %f to %f\n", min_conv, max_conv);
    scale_factor = 255.0 / max_conv;
    printf("scale factor: %f\n", scale_factor);
    for (i = 0; i < height; i++) {
        for (j = 0; j < width; j++)
            image[i][j] = (unsigned char) floor(rimage[i][j] * scale_factor + 0.5);
    }
}

rescale_image_2(image)
unsigned char **image;
{
    int i, j;
    int max_conv, min_conv;
    float scale_factor;


    /* rescale image to be from 1 to 254 */
    max_conv = 0;
    min_conv = 256;
    for (i = 0; i < height; i++)
        for (j = 0; j < width; j++)
            if (image[i][j] > max_conv) {
                max_conv = image[i][j];
            }
            else if (image[i][j] < min_conv) {
                min_conv = image[i][j];
            }
    printf("image range %d to %d\n", min_conv, max_conv);
    scale_factor = 254.0 / (float) max_conv;
    printf("scale factor: %f\n", scale_factor);
    for (i = 0; i < height; i++) {
        for (j = 0; j < width; j++)
            image[i][j] = floor((float) image[i][j] * scale_factor + 0.5);
    }
}

/* exclude borders */
rescale_image_3(rimage, image, border)
float **rimage;
unsigned char **image;
int border;
{
    int i, j;
    float max_conv, min_conv;
    float scale_factor;

    /* rescale image to be from 1 to 254 */
    max_conv = rimage[0][0];
    min_conv = rimage[0][0];
    for (i = border; i < height-border; i++)
        for (j = border; j < width-border; j++)
            if (rimage[i][j] > max_conv) {
                max_conv = rimage[i][j];
            }
            else if (rimage[i][j] < min_conv) {
                min_conv = rimage[i][j];
            }
    printf("image range %f to %f\n", min_conv, max_conv);
    scale_factor = 255.0 / max_conv;
    printf("scale factor: %f\n", scale_factor);
    for (i = 0; i < height; i++) {
        for (j = 0; j < width; j++)
            image[i][j] = (unsigned char) floor(rimage[i][j] * scale_factor + 0.5);
    }

    /* zero borders */
    for (i = 0; i < height; i++)
        for (j = 0; j < border; j++)
            image[i][j] = image[i][width-j-1] = 0;
    for (i = 0; i < border; i++)
        for (j = 0; j < width; j++)
            image[i][j] = image[height-i-1][j] = 0;
}

/* ************************************************************************ */
/* * SUPPRESSION_1 *                                                        */
/* *****************                                                        */
suppression_1(i, j, ux, uy, g0_ptr)
int i, j;
double ux, uy, *g0_ptr;
{
    int continue_flag = FALSE;

    if (fabs(ux) < fabs(uy)) {
        if (((*g0_ptr = fabs(uy * g[i][j])) <= fabs((ux * g[i + 1][j + 1]) + ((uy - ux) * g[i][j + 1]))) ||
            (*g0_ptr <= fabs((ux * g[i - 1][j - 1]) + ((uy - ux) * g[i][j - 1]))))
            continue_flag = TRUE;
    }
    else {
        if (((*g0_ptr = fabs(ux * g[i][j])) <= fabs((uy * g[i + 1][j + 1]) + ((ux - uy) * g[i + 1][j]))) ||
            (*g0_ptr <= fabs((uy * g[i - 1][j - 1]) + ((ux - uy) * g[i - 1][j]))))
            continue_flag = TRUE;
    };

    return (continue_flag);
}

/* ************************************************************************ */
/* * SUPPRESSION_2 *                                                        */
/* *****************                                                        */
suppression_2(i, j, ux, uy, g0_ptr)
int i, j;
double ux, uy, *g0_ptr;
{
    int continue_flag = FALSE;

    if (fabs(ux) < fabs(uy)) {
        if (((*g0_ptr = fabs(uy * g[i][j])) <= fabs((ux * g[i + 1][j - 1]) - ((uy + ux) * g[i][j - 1]))) ||
            (*g0_ptr <= fabs((ux * g[i - 1][j + 1]) - ((uy + ux) * g[i][j + 1]))))
            continue_flag = TRUE;
    }
    else {
        if (((*g0_ptr = fabs(ux * g[i][j])) <= fabs((uy * g[i + 1][j - 1]) - ((ux + uy) * g[i + 1][j]))) ||
            (*g0_ptr <= fabs((uy * g[i - 1][j + 1]) - ((ux + uy) * g[i - 1][j]))))
            continue_flag = TRUE;
    };

    return (continue_flag);
}

/* ************************************************************************ */
/* * GAUSS_EXPON *                                                          */
/* ***************                                                          */
double gauss_expon(x, s)
double x, s;
{
    return (exp((-(x * x)) / (2 * s * s)));
}