recursive_filter.c

canny 边缘检测算法

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


recursive_filter_horizontal(image1, image2, width, height, sigma, option)
float **image1, **image2;
int width,height;
float sigma;
int option;
{
	/* filter image using 1d filter in horizontal direction type of
	   filter depends on option: 1:	gaussian 2:	1st derivative 3:
	   2nd derivative coeeficients are computed depending on the option */
	int x, y;
	float a0, a1, b0, b1, c0, c1, w0, w1;
	float w0n, w1n, cos0, cos1, sin0, sin1, b0n, b1n;
	float n00, n11, n22, n33;
	float n11M, n22M, n33M, n44M;
	float d11, d22, d33, d44;

	float *yP, *yN;

	yN = (float *) malloc(sizeof(float) * width);
	if(yN == NULL){
		printf("cant allocate yN in recursive_filter_horizontal\n");
		exit(-1);
	}
	yP = (float *) malloc(sizeof(float) * width);
	if(yP == NULL){
		printf("cant allocate yP in recursive_filter_horizontal\n");
		exit(-1);
	}


	/* coefficients for Gaussian filter */
	if (option == GAUSS) {
		a0 = 1.68;
		a1 = 3.735;
		b0 = 1.783;
		b1 = 1.723;
		c0 = -0.6803;
		c1 = -0.2598;
		w0 = 0.6318;
		w1 = 1.997;
		printf("computing gaussian horizontally\n");
	}
	else if (option == FIRST_DERIV) {
		a0 = 0.6472;
		a1 = 4.531;
		b0 = 1.527;
		b1 = 1.516;
		c0 = -0.6494;
		c1 = -0.9557;
		w0 = 0.6719;
		w1 = 2.072;
		printf("computing first derivative horizontally\n");
	}
	else if (option == SECOND_DERIV) {
		a0 = 1.331;
		a1 = -3.661;
		b0 = 1.24;
		b1 = 1.314;
		c0 = -0.3225;
		c1 = 1.738;
		w0 = 0.748;
		w1 = 2.166;
		printf("computing second derivative of the gaussian horizontally\n");
	}
	else {
		printf("error - unknown option - aborting\n");
		exit(-1);
	}

	/* coefficients for filter */
	w0n = w0 / sigma;
	w1n = w1 / sigma;
	cos0 = cos(w0n);
	cos1 = cos(w1n);
	sin0 = sin(w0n);
	sin1 = sin(w1n);
	b0n = b0 / sigma;
	b1n = b1 / sigma;

	n33 = exp(-b1n - 2 * b0n) * (c1 * sin1 - cos1 * c0) + exp(-b0n - 2 * b1n) * (a1 * sin0 - cos0 * a0);
	n22 = 2 * exp(-b0n - b1n) * ((a0 + c0) * cos1 * cos0 - cos1 * a1 * sin0 - cos0 * c1 * sin1) + c0 * exp(-2 * b0n) + a0 * exp(-2 * b1n);
	n11 = exp(-b1n) * (c1 * sin1 - (c0 + 2 * a0) * cos1) + exp(-b0n) * (a1 * sin0 - (2 * c0 + a0) * cos0);
	n00 = a0 + c0;

	d44 = exp(-2 * b0n - 2 * b1n);
	d33 = -2 * cos0 * exp(-b0n - 2 * b1n) - 2 * cos1 * exp(-b1n - 2 * b0n);
	d22 = 4 * cos1 * cos0 * exp(-b0n - b1n) + exp(-2 * b1n) + exp(-2 * b0n);
	d11 = -2 * exp(-b1n) * cos1 - 2 * exp(-b0n) * cos0;

	n44M = -d44 * n00;
	n33M = n33 - d33 * n00;
	n22M = n22 - d22 * n00;
	n11M = n11 - d11 * n00;

	if(option == FIRST_DERIV){
		n11M = -n11M;
		n22M = -n22M;
		n33M = -n33M;
		n44M = -n44M;
	}

	/* clear temp arrays */
	for (x = 0; x < width; x++)
			yP[x] = yN[x] = 0;

	/* filter in positive direction - left to right */
	for (y = 0; y < height; y++){
		/* filter in positive direction - left to right */
		for (x = 4; x < width; x++) {
			yP[x] = n00 * (float) image1[y][x] +
				n11 * (float) image1[y][x - 1] +
				n22 * (float) image1[y][x - 2] +
				n33 * (float) image1[y][x - 3] -
				d11 * yP[x - 1] -
				d22 * yP[x - 2] -
				d33 * yP[x - 3] -
				d44 * yP[x - 4];
		}
		/* filter in negative direction - right to left */
		for (x = width - 5; x >= 0; x--) {
			yN[x] = n11M * (float) image1[y][x + 1] +
				n22M * (float) image1[y][x + 2] +
				n33M * (float) image1[y][x + 3] +
				n44M * (float) image1[y][x + 4] -
				d11 * yN[x + 1] -
				d22 * yN[x + 2] -
				d33 * yN[x + 3] -
				d44 * yN[x + 4];
		}
		for (x = 0; x < width; x++)
			image2[y][x] = yP[x] + yN[x];
	}
	free(yP);
	free(yN);
}


recursive_filter_vertical(image1, image2, width,height,sigma,option)
float **image1, **image2;
int width,height;
float sigma;
int option;
{
	/* filter image using 1d filter in horozontal direction type of
	   filter depends on option: 1:	gaussian 2:	1st derivative 3:
	   2nd derivative coeeficients are computed depending on the option */
	int x, y;
	float a0, a1, b0, b1, c0, c1, w0, w1;
	float w0n, w1n, cos0, cos1, sin0, sin1, b0n, b1n;
	float n00, n11, n22, n33;
	float n11M, n22M, n33M, n44M;
	float d11, d22, d33, d44;
	
	float *yP, *yN;

	yN = (float *) malloc(sizeof(float) * width);
	if(yN == NULL){
		printf("cant allocate yN in recursive_filter_horizontal\n");
		exit(-1);
	}
	yP = (float *) malloc(sizeof(float) * width);
	if(yP == NULL){
		printf("cant allocate yP in recursive_filter_horizontal\n");
		exit(-1);
	}

	/* coefficients for Gaussian filter */
	if (option == GAUSS) {
		a0 = 1.68;
		a1 = 3.735;
		b0 = 1.783;
		b1 = 1.723;
		c0 = -0.6803;
		c1 = -0.2598;
		w0 = 0.6318;
		w1 = 1.997;
		printf("computing gaussian vertically\n");
	}
	else if (option == FIRST_DERIV) {
		a0 = 0.6472;
		a1 = 4.531;
		b0 = 1.527;
		b1 = 1.516;
		c0 = -0.6494;
		c1 = -0.9557;
		w0 = 0.6719;
		w1 = 2.072;
		printf("computing first derivative vertically\n");
	}
	else if (option == SECOND_DERIV) {
		a0 = 1.331;
		a1 = -3.661;
		b0 = 1.24;
		b1 = 1.314;
		c0 = -0.3225;
		c1 = 1.738;
		w0 = 0.748;
		w1 = 2.166;
		printf("computing second derivative of the gaussian vertically\n");
	}
	else {
		printf("error - unknown option - aborting\n");
		exit(-1);
	}

	/* coefficients for filter */
	w0n = w0 / sigma;
	w1n = w1 / sigma;
	cos0 = cos(w0n);
	cos1 = cos(w1n);
	sin0 = sin(w0n);
	sin1 = sin(w1n);
	b0n = b0 / sigma;
	b1n = b1 / sigma;

	n33 = exp(-b1n - 2 * b0n) * (c1 * sin1 - cos1 * c0) + exp(-b0n - 2 * b1n) * (a1 * sin0 - cos0 * a0);
	n22 = 2 * exp(-b0n - b1n) * ((a0 + c0) * cos1 * cos0 - cos1 * a1 * sin0 - cos0 * c1 * sin1) + c0 * exp(-2 * b0n) + a0 * exp(-2 * b1n);
	n11 = exp(-b1n) * (c1 * sin1 - (c0 + 2 * a0) * cos1) + exp(-b0n) * (a1 * sin0 - (2 * c0 + a0) * cos0);
	n00 = a0 + c0;

	d44 = exp(-2 * b0n - 2 * b1n);
	d33 = -2 * cos0 * exp(-b0n - 2 * b1n) - 2 * cos1 * exp(-b1n - 2 * b0n);
	d22 = 4 * cos1 * cos0 * exp(-b0n - b1n) + exp(-2 * b1n) + exp(-2 * b0n);
	d11 = -2 * exp(-b1n) * cos1 - 2 * exp(-b0n) * cos0;

	n44M = -d44 * n00;
	n33M = n33 - d33 * n00;
	n22M = n22 - d22 * n00;
	n11M = n11 - d11 * n00;
	
	if(option == FIRST_DERIV){
		n11M = -n11M;
		n22M = -n22M;
		n33M = -n33M;
		n44M = -n44M;
	}

	/* clear temp images */
	for (y = 0; y < height; y++)
			yP[y] = yN[y] = 0;

	/* filter vertically downwards */
	for (x = 0; x < width; x++){
		for (y = 4; y < height; y++) {
			yP[y] = n00 * image1[y][x] +
				n11 * image1[y - 1][x] +
				n22 * image1[y - 2][x] +
				n33 * image1[y - 3][x] -
				d11 * yP[y - 1] -
				d22 * yP[y - 2] -
				d33 * yP[y - 3] -
				d44 * yP[y - 4];
		}
		/* filter vertically upwards */
		for (y = height - 5; y > 0; y--) {
			yN[y] = n11M * image1[y + 1][x] +
				n22M * image1[y + 2][x] +
				n33M * image1[y + 3][x] +
				n44M * image1[y + 4][x] -
				d11 * yN[y + 1] -
				d22 * yN[y + 2] -
				d33 * yN[y + 3] -
				d44 * yN[y + 4];
		}

		/* sum forward and backward parts */
		for (y = 0; y < height; y++)
			image2[y][x] = yP[y] + yN[y];
	}
	free(yN);
	free(yP);
}