DÉTECTION OPTIMALE DES CONTOURS: CANNY, DÉRICHE, SHENCASTAN

//*****************************************************************************
//MORARD Vincent					
//04 mars 2008
//vincent.morard@cpe.fr
//http://pistol.petesampras.free.fr
//*****************************************************************************



//******************************************************************************
//Shen_Castan.cpp
//Ce fichier regroupe les fonctions permettant de trouver les contours selon les 
//critères de Shen-Castan. 
//OUTLINE corespond au nombre de pixels que l'on ne prendra pas en compte lors 
//de la convolution.
//Aussi, pour pouvoir effectuer la détection de contour à tous les pixels, y 
//compris ceux du bord de l'image, on augmentera la taille de l'image de 
//2*OUTLINE. C'est le rôle des fonctions embded et debed
//WINDOW_SIZE est un paramètre qui determine la taille de la fenêtre pour 
//calculer le gradient adaptatif.
//
//plus b est faible plus il y a de détail (plus de bruit aussi)
//******************************************************************************
#include "CImage.h"
#include "AdvancEdge.h"

#define WINDOW_SIZE		7	
#define OUTLINE 		25
double b;


//*************************************************************************
//ShenCastan: Detection des contours optimaux:
//Procedure: On augmente d'abord la taille de l'image source en on convertit 
//l'image initiale en float. On envoie cette 
//nouvelle image à la fonction Shen avec ses nouvelles dimensions
//Une fois la fin du traitement, on remet les dimensions initiales et l'on place 
//l'image obtenue dans le buffer de sortie.
//*************************************************************************
bool CImage::ShenCastan(CImage *ImgDest,double B)
{
	float **Img=0;
	int Largeur,Hauteur;

	b=B;
	if(hBmp==0){
		MessageBox(NULL,"Shen-Castan : L'image source est vide",
			NULL,MB_OK|MB_ICONWARNING);
		return 0;
	}

	if(ImgDest != 0 && ImgDest != this)
		ImgDest->Copy(this);

	

	if(ImgDest == 0)
		ImgDest=this;

	GetBitmapBits(hBmp,Width*Height*4,ucBits);

	//Ajout d'une bordure de OUTLINE pixel autours de l'image
	Img=embed(OUTLINE,&Largeur,&Hauteur);

	Shen(Img,Largeur,Hauteur);
	
	//Retrait de la bordure
	ImgDest->debed(Img,OUTLINE,Largeur,Hauteur);

	SetBitmapBits(ImgDest->hBmp,Width*Height*4,ImgDest->ucBits);
	ImgDest->ImgType=GRAY;

	DesAllocT_float(Img,Largeur);
	return 1;
}


//******************************************************************************
//Shen: 
//C'est dans cette fonction que tout le traitement est effectué
//On reçoit en paramètre d'entrée l'image directement accessible en pixel 
//ainsi que ses dimensions.
//L'image de sortie sera dans la variable Img
//******************************************************************************
void Shen(float **Img,int Largeur,int Hauteur)
{
	float **BufFiltrer=0;
	int   **ImgBli=0;

	BufFiltrer=AllocT_float(Largeur,Hauteur);

	//On filtre le bruit en appliquant l'algo ISEF
	ComputeISEF(Img,BufFiltrer,Largeur,Hauteur);
	
	//On trouve les pixels de l'image qui on un Laplacian positif 
	ImgBli=ComputeBli(BufFiltrer,Img,Largeur,Hauteur);

	//Detection des coutours à partir de l'ImgBli et de l'image filtrer
	LocateZeroCrossings(Img,BufFiltrer,ImgBli,Largeur,Hauteur);

	MaxContraste(Img,Largeur,Hauteur);

	//Desallocation de la memoire allouer
	DesAllocT_float(BufFiltrer,Largeur);
	DesAllocT_int(ImgBli,Largeur);


	
}

//******************************************************************************
//ISEF  (Infinite Symetrical Exponential Filter)
//Filtrage de l'image horizontalement et verticalement. L'image filtrer sera 
//placée dans la variable y
//******************************************************************************
void ComputeISEF(float **x,float **y,int Largeur,int Hauteur)
{
	float **Causal=0,**AntiCausal=0;

	Causal=AllocT_float(Largeur,Hauteur);
	AntiCausal=AllocT_float(Largeur,Hauteur);

	//On applique d'abord le filtre dans la direction verticale
	ApplyISEF_Vertical(x,y,Causal,AntiCausal,Largeur,Hauteur);
	ApplyISEF_Horizontal(y,y,Causal,AntiCausal,Largeur,Hauteur);
	
	

	//Libération de la mémoire
	DesAllocT_float(Causal,Largeur);
	DesAllocT_float(AntiCausal,Largeur);
	

}


//*******************************************************************************
//ApplyISEF_Vertical
//Filtrage vertical : Calcul des composantes causales et anticausales
//*******************************************************************************
void ApplyISEF_Vertical(float **x,float **y,float **Causal,float **AntiCausal,
						int Largeur,int Hauteur)
{
	int i,j;
	float b1,b2;
	b1 = (float)((1.0-b)/(1.0+b));
	b2 = (float)(b*b1);

	//Methode recurssive donc on calcule les pixels des bords de l'image
	for(j=0;j<Hauteur;j++)
	{
		Causal[0][j]=b1*x[0][j];
		AntiCausal[Largeur-1][j]=b2*x[Largeur-1][j];
	}

	//Calcul des composantes causales
	for(i=1;i<Largeur;i++)
		for(j=0;j<Hauteur;j++)
			Causal[i][j]=(float)(b1*x[i][j]+b*Causal[i-1][j]);

	//Calcul des composantes anti-causales
	for(i=Largeur-2;i>=0;i--)
		for(j=0;j<Hauteur;j++)
			AntiCausal[i][j]=(float)(b2*x[i][j]+b*AntiCausal[i+1][j]); 

	//on calcule les pixels des bords de l'image de sortie
	for(j=0;j<Hauteur-1;j++)
		y[Largeur-1][j]=Causal[Largeur-1][j];

	
	//On calcule l'image de sortie du premier filtre que l'on place dans la variable y
	//Correspond a la somme des composantes causal et anti-causales
	for(i=0;i<Largeur-2;i++)				
		for(j=0;j<Hauteur-1;j++)
			y[i][j]=Causal[i][j]+AntiCausal[i+1][j];

}

//*******************************************************************************
//ApplyISEF_Horizontal
//Filtrage vertical : Calcul des composantes causales et anticausales
//*******************************************************************************
void ApplyISEF_Horizontal(float **x,float **y,float **Causal,float **AntiCausal,
						  int Largeur,int Hauteur)
{
	int i,j;
	float b1,b2;
	b1 = (float)((1.0-b)/(1.0+b));
	b2 = (float)(b*b1);

	//on calcule les pixels des bords de l'image
	for(i=0;i<Largeur;i++)
	{
		Causal[i][0]=b1*x[i][0];
		AntiCausal[i][Hauteur-2]=b2*x[i][Hauteur-2];
	}

	//Calcul des composantes causales
	for(j=1;j<Hauteur;j++)
		for(i=0;i<Largeur;i++)
			Causal[i][j]=(float)(b1*x[i][j]+b*Causal[i][j-1]);

	//Compute anti causal component
	for(j=Hauteur-3;j>=0;j--)
		for(i=0;i<Largeur;i++)
			AntiCausal[i][j]=(float)(b2*x[i][j]+b*AntiCausal[i][j+1]); 

	//Calcul des composantes anti-causales
	for(i=0;i<Largeur-1;i++)
		y[i][Hauteur-1]=Causal[i][Hauteur-1];

	//On calcule l'image de sortie du premier filtre que l'on place dans la variable y
	//Correspond a la somme des composantes causal et anti-causales
	for(i=0;i<Largeur-1;i++)
		for(j=0;j<Hauteur-2;j++)
			y[i][j]=Causal[i][j]+AntiCausal[i][j+1];

}

//*******************************************************************************
//ComputeBli:
//On fait la différence des deux images et on compare le resultat à 0
//ImgBli est donc une image composé uniquement de 0 et de 1
//*******************************************************************************
int **ComputeBli(float **ImgFiltrer,float **ImgBuf,int nRows,int nCols)
{
	int Row,Col;
	int **ImgBli=0;

	ImgBli=AllocT_int(nRows,nCols);

	//On prend la difference entre l'image lisse et l'image originale.
	//On calcule l'ImgBli en mettant a 1 tous les pixels ou le Laplacian est positif. 0 sinon
	for(Row=0;Row<nRows;Row++)
	{
		for(Col=0;Col<nCols;Col++)
		{
			ImgBli[Row][Col]=0;
			if(Row<OUTLINE || Row>=nRows-OUTLINE || Col<OUTLINE || Col>=nCols-OUTLINE)
				continue;
			ImgBli[Row][Col]=((ImgFiltrer[Row][Col]-ImgBuf[Row][Col])>0.0);

		}
	}
	return ImgBli;

}


//*************************************************************************************
//LocateZeroCrossings
//Cette fonction permettra de déterminer pour tous les pixels de l'image s'il est un pixel
//appartenant à un contour ou non.
//**************************************************************************************
void LocateZeroCrossings(float **Orig,float **BufFiltrer,int **ImgBli,int nRows,int nCols)
{
	int Row,Col;
	
	for(Row=0;Row<nRows;Row++)
		for(Col=0;Col<nCols;Col++)
		{
			//On ignore les pixels que l'on a ajoute pour le calcule
			if(Row<OUTLINE || Row>=nRows-OUTLINE || Col<OUTLINE || Col>=nCols-OUTLINE)
				Orig[Row][Col]=0.0;

			//On verifie si ce pixel est un "zero crossing" pour le Laplacian
			else if(IsCandidateEdge(ImgBli,BufFiltrer,Row,Col))
			{
				//On Calcule le gradian adaptatif
				Orig[Row][Col]=ComputeAdaptativeGradient(ImgBli,BufFiltrer,Row,Col);
				
			}
			else
				Orig[Row][Col]=0.0;
		}
}

//**************************************************************************************
//IsCandidateEdge
//On regarde le pixel voisin en on regarde s'il y a un franchissement de zero.
//On effectue donc la multiplication des 2 pixels et on compare à 0
//************************************************************************************** 
bool IsCandidateEdge(int **Buff,float **Orig,int Row,int Col)	
{
	//a positive z-c must have a positive 1st derivative,where positive z-c
	//means the second derivative goes from + to - as we cross the edge

	if(Buff[Row][Col]==1 && Buff[Row+1][Col]==0) 
		return (Orig[Row+1][Col]-Orig[Row-1][Col]>0 ? TRUE:FALSE); //positive z-c

	else if(Buff[Row][Col]==1 && Buff[Row][Col+1]==0)
		return (Orig[Row][Col+1]-Orig[Row][Col-1]>0 ? TRUE:FALSE); //positive z-c

	else if(Buff[Row][Col]==1 && Buff[Row-1][Col]==0)
		return (Orig[Row+1][Col]-Orig[Row-1][Col]<0 ? TRUE:FALSE); //negative z-c

	else if(Buff[Row][Col]==1 && Buff[Row][Col-1]==0)
		return (Orig[Row][Col+1]-Orig[Row][Col-1]<0 ? TRUE:FALSE); //negative z-c


	return FALSE;		//not a z-c

}  

//*************************************************************************************
//ComputeAdaptativeGradient
//On calcule un seuil pour chaque pixel. Le seuil sera determiner grace aux pixels voisins
//présent dans la fenêtre WINDOW_SIZE
//*************************************************************************************
float ComputeAdaptativeGradient(int **Bli,float **Orig,int Row,int Col)
{
	int i,j;
	float SumOn,SumOff;
	float AvgOn,AvgOff;
	int NbOn,NbOff;

	SumOn=0.0;
	SumOff=0.0;
	NbOn=0;
	NbOff=0;

	//On regarde par rapport aux pixels voisins, le nombre de pixel a 1 et a 0
	for(i=(-WINDOW_SIZE/2);i<=(WINDOW_SIZE/2);i++)
		for(j=(-WINDOW_SIZE/2);j<=(WINDOW_SIZE/2);j++)
		{
			if(Bli[Row+i][Col+j])
			{
				SumOn+=Orig[Row+i][Col+j];
				NbOn++;
			} 
			else
			{
				SumOff+=Orig[Row+i][Col+j];
				NbOff++;
			}  
		}

		if(SumOff)
			AvgOff=SumOff/(float)NbOff;
		else 
			AvgOff=0.0;

		if(SumOn)
			AvgOn=SumOn/(float)NbOn;
		else
			AvgOn=0.0;
		return (AvgOff-AvgOn);
}



//********************************************************************************
//embed
//Cette fonction retourne un pointeur sur un float qui de l'image initial ou l'on a 
//ajouté des bordure de largeur et de hauteur de Width
//********************************************************************************
float **CImage::embed(int taille,int *L,int *H)
{
	int i,j,I,J;

	float **Img=0;


	taille+=2;
	*L=Width+2*taille;
	*H=Height+2*taille;
	Img=AllocT_float(*L,*H);

	for(i=0;i<*L;i++)
		for(j=0;j<*H;j++)
		{
			I=(i-taille+Width)%Width;
			J=(j-taille+Height)%Height;
			Img[i][j]=(float)GetPixel(I,J,GRAY);		
		}

		return Img;
}

//********************************************************************************
//debed
//Cette fonction retourne inscrit l'image obtenue dans le buffer Dest en enlevant 
//les bordures ajoutées
//********************************************************************************
void CImage::debed(float **Img,int taille,int L,int H)
{
	int i,j;
	taille+=2;
	int Val;
	for(i=taille;i<L-taille;i++)
		for(j=taille;j<H-taille;j++)
		{
			Val=Limit((int)Img[i][j]);

			SetPixel(i-taille,j-taille,Val,Val,Val);		
		}
}


//Etire le contraste au maximum.
//Attention tous les pixels negatifs sont mis a 0
void MaxContraste(float **Img,int nRows,int nCols)
{
	int i,j;
	float X,Scale,Vmin,Vmax;
	Vmin=Img[50][50];
	Vmax=Vmin;
	for(i=0;i<nRows;i++)
	{
		for(j=0;j<nCols;j++)
		{
			if(i<OUTLINE || i>=nRows-OUTLINE || j<OUTLINE || j>=nCols-OUTLINE)
				continue;
			X=Img[i][j];
			if(X<0){Img[i][j]=0;X=0;}
			if(Vmin>X) Vmin=X;
			if(Vmax<X) Vmax=X;
		}
	}

	Scale = (float)(256.0/(Vmax-Vmin+1));

	for(i=0;i<nRows;i++)
		for(j=0;j<nCols;j++)
		{
			if(i<OUTLINE || i>=nRows-OUTLINE || j<OUTLINE || j>=nCols-OUTLINE)
				continue;
			
			Img[i][j]=((Img[i][j]-Vmin)*Scale);

		}



}