convolve.cpp

基于小波的图像配准

/* This procedure is for transforming array to array by convolving
   source: given two types, one is double,another is unsigned char
   target: double array
   coefficient: parameters for convolving
   substar,subend: the start and end of subscript of coefficient
*/

#include "convolve.h"
#include "grayscale.h"	

void dc_convolve(double* source,unsigned char* target,
			     long width,long height,
			     double* coefficient,int substar,int subend)
{
	long i,j,k,l,row,ij,dij,sub1,sub2,m,n;
	int count;
	int LExtend=0,RExtend=0; 
	double* value;
	value = (double *)malloc(sizeof(double)*((width/2)*height));

	count = subend - substar +1;
	if(substar<0)
		LExtend = abs(substar);
	if(subend>0)
		RExtend = subend;
	
	for(j=0;j<height;j++){
		row = j* width;
		for(i=0;i<width/2;i++){
			
			ij = height*i + j; //result subscript
			dij = row + 2*i; // data subscript
			
			value[ij] = 0;

			//processing the  left boundary that should be extended
			//边界对称扩展
			sub1 = LExtend-2*i;
			if(sub1>0){
				if(sub1>count) m = sub1 + subend + 1;
				else m = sub1;
				for(k=0;k<m;k++){
					n = row+sub1-k;
					//越界处理
					if(n<row)  n = row;
					if(n>=(row+width)) n = row+width-1;

					value[ij] += coefficient[k] * source[n];
				}
				for(l=0;l<(count-m);l++){
					n = row+l;
					//越界处理
					if(n<row) n = row;
					if(n>=(row+width)) n = row+width-1;

					value[ij] += coefficient[l+k] * source[n];
				}
			}

			//processing the  right boundary that should be extended
			//边界对称扩展
			sub2 = RExtend - (width-2*i-1);
			if(sub2>0){
				if(sub2>count) m = sub2- substar + 1;
				else m = sub2;
				for(k=0;k<(count-m);k++){
					n = row+width-(count-m)+k;
					//越界处理
					if(n<row) n = row;
					if(n>=(row+width)) n = row+width-1;

					value[ij] += coefficient[k] * source[n];
				}
				for(l=0;l<m;l++){
					n = row+width-1-sub2+l;
					//越界处理
					if(n<row) n = row;
					if(n>=(row+width)) n = row+width-1;

					value[ij] += coefficient[count-1-l] * source[n];
				}
			}

			//processing the left pixels with low filter
			if(sub1<=0&&sub2<=0)
			for(k=0;k<count;k++){
				value[ij] += coefficient[k] * source[dij+substar+k];
			}			 
		}
	}
	grayscale(width/2,height,value,target);
	free(value);
	return;
}

void dd_convolve(double* source,double* target,
			     long width,long height,
			     double* coefficient,int substar,int subend)
{
	long i,j,k,l,row,ij,dij,sub1,sub2,m,n;
	int count;
	int LExtend=0,RExtend=0; 

	count = subend - substar +1;
	if(substar<0)
		LExtend = abs(substar);
	if(subend>0)
		RExtend = subend;
	
	for(j=0;j<height;j++){
		row = j* width;
		for(i=0;i<width/2;i++){
			
			ij = height*i + j; //result subscript
			dij = row + 2*i; // data subscript
			
			target[ij] = 0;

			//processing the  left boundary that should be extended
			//边界对称扩展
			sub1 = LExtend-2*i;
			if(sub1>0){
				if(sub1>count) m = sub1 + subend + 1;
				else m = sub1;
				for(k=0;k<m;k++){
					n = row+sub1-k;
					//越界处理
					if(n<row)  n = row;
					if(n>=(row+width)) n = row+width-1;

					target[ij] = target[ij]+
						coefficient[k] * source[n];
				}
				for(l=0;l<(count-m);l++){
					n = row+l;
					//越界处理
					if(n<row) n = row;
					if(n>=(row+width)) n = row+width-1;

					target[ij] = target[ij]+
						coefficient[l+k] * source[n];
				}
			}

			//processing the  right boundary that should be extended
			//边界对称扩展
			sub2 = RExtend - (width-2*i-1);
			if(sub2>0){
				if(sub2>count) m = sub2- substar + 1;
				else m = sub2;
				for(k=0;k<(count-m);k++){
					n = row+width-(count-m)+k;
					//越界处理
					if(n<row) n = row;
					if(n>=(row+width)) n = row+width-1;

					target[ij] = target[ij]+
						coefficient[k] * source[n];
				}
				for(l=0;l<m;l++){
					n = row+width-1-sub2+l;
					//越界处理
					if(n<row) n = row;
					if(n>=(row+width)) n = row+width-1;

					target[ij] = target[ij]+
						coefficient[count-1-l] * source[n];
				}
			}

			//processing the left pixels with low filter
			if(sub1<=0&&sub2<=0)
			for(k=0;k<count;k++){
				target[ij] = target[ij]+
					coefficient[k] * source[dij+substar+k];
			}
		}
	}
	return;
}

void cd_convolve(unsigned char* source,double* target,
			     long width,long height,
			     double* coefficient,int substar,int subend)
{
	long i,j,k,l,row,ij,dij,sub1,sub2,m,n;
	int count;
	int LExtend=0,RExtend=0; 

	count = subend - substar +1;
	if(substar<0)
		LExtend = abs(substar);
	if(subend>0)
		RExtend = subend;
	
	for(j=0;j<height;j++){
		row = j* width;
		for(i=0;i<width/2;i++){
			
			ij = height*i + j; //result subscript
			dij = row + 2*i; // data subscript
			
			target[ij] = 0;

			//processing the  left boundary that should be extended
			//边界对称扩展
			sub1 = LExtend-2*i;
			if(sub1>0){
				if(sub1>count) m = sub1 + subend + 1;
				else m = sub1;
				for(k=0;k<m;k++){
					n = row+sub1-k;
					//越界处理
					if(n<row)  n = row;
					if(n>=(row+width)) n = row+width-1;

					target[ij] = target[ij]+
						coefficient[k] * source[n];
				}
				for(l=0;l<(count-m);l++){
					n = row+l;
					//越界处理
					if(n<row) n = row;
					if(n>=(row+width)) n = row+width-1;

					target[ij] = target[ij]+
						coefficient[l+k] * source[n];
				}
			}

			//processing the  right boundary that should be extended
			//边界对称扩展
			sub2 = RExtend - (width-2*i-1);
			if(sub2>0){
				if(sub2>count) m = sub2- substar + 1;
				else m = sub2;
				for(k=0;k<(count-m);k++){
					n = row+width-(count-m)+k;
					//越界处理
					if(n<row) n = row;
					if(n>=(row+width)) n = row+width-1;

					target[ij] = target[ij]+
						coefficient[k] * source[n];
				}
				for(l=0;l<m;l++){
					n = row+width-1-sub2+l;
					//越界处理
					if(n<row) n = row;
					if(n>=(row+width)) n = row+width-1;

					target[ij] = target[ij]+
						coefficient[count-1-l] * source[n];
				}
			}

			//processing the left pixels with low filter
			if(sub1<=0&&sub2<=0)
			for(k=0;k<count;k++){
				target[ij] = target[ij]+
					coefficient[k] * source[dij+substar+k];
			}
		}
	}
	return;
}