LinuxQuestions.org
Review your favorite Linux distribution.
Home Forums Tutorials Articles Register
Go Back   LinuxQuestions.org > Forums > Non-*NIX Forums > Programming
Reload this Page C code for median filter on color images
User Name
Password
Programming This forum is for all programming questions.
The question does not have to be directly related to Linux and any language is fair game.

Notices


Reply
Page 7 of 9 « First 56 7 89
  Search this Thread
Old 03-20-2012, 05:52 AM   #91
jayshri
Member
 
Registered: Feb 2012
Posts: 71

Rep: Reputation: Disabled

hi,
Code:
#include<stdint.h>
#include<stdio.h>
#include<stdlib.h>
#include<float.h>
#include<math.h>
#include<string.h>
#include<limits.h>
#define MAX_BRIGHTNESS 255
#define M_PI 3.14159265358979323846264338327
//#define OUTPUT_FILE "out.bmp"
struct bmpfile_magic{
unsigned char magic[2];
};
 struct bmpfile_header{
uint32_t filesz;
uint16_t creator1;
uint16_t creator2;
uint32_t  bmp_offset;
};
typedef struct{
uint32_t header_sz;
int32_t width;
int32_t height;
uint16_t nplanes;
uint16_t bitspp;
uint32_t compress_type;
uint32_t bmp_bytesz;
int32_t hres;
int32_t vres;
uint32_t ncolors;
uint32_t nimpcolors;
}BITMAPINFOHEADER;
typedef struct{
uint8_t r;
uint8_t g;
uint8_t b;
uint8_t null;
}rgb_t;
BITMAPINFOHEADER ih;


typedef int pixel_t;
pixel_t *load_bmp( char *filename,BITMAPINFOHEADER *bitmapInfoHeader)
{
  FILE *filePtr;
struct bmpfile_magic mag;
struct bmpfile_header bitmapFileHeader;
pixel_t *bitmapImage;
filePtr=fopen(filename,"r");
if(filePtr==NULL)
{
perror("fopen()");
exit(1);
}

fread(&mag,sizeof(struct bmpfile_magic),1,filePtr);
//verify that this is a bmp file by check bitmap id
if(*((uint16_t*) mag.magic)!=0x4D42)
{
fprintf(stderr,"Not a bmp file:magic=%c%c\n",mag.magic[0],mag.magic[1]);
fclose(filePtr);
return NULL;
}
//read the bitmap header
fread(&bitmapFileHeader,sizeof(struct bmpfile_header),1,filePtr);

//read the bitmap info header
fread(bitmapInfoHeader,sizeof(BITMAPINFOHEADER),1,filePtr);

if(bitmapInfoHeader->compress_type!=0)
fprintf(stderr,"Warning,compression is not supported.\n");
//move file point to the beginning of bitmap data
fseek(filePtr,bitmapFileHeader.bmp_offset,SEEK_SET);

//allocate enough memory for bitmao image data
bitmapImage=(pixel_t*)malloc(bitmapInfoHeader->bmp_bytesz*(sizeof(pixel_t)));
//verify memory allocation
if(!bitmapImage)
{
free(bitmapImage);
fclose(filePtr);
return NULL;
}

//read in the bitmap image data
int i;
unsigned char c;
for( i=0;i<bitmapInfoHeader->bmp_bytesz;i++){
fread(&c,sizeof(unsigned char),1,filePtr);
bitmapImage[i]=(int)c;
}
fclose(filePtr);
return bitmapImage;
}

//return:nonzero on error.
int save_bmp(char *filename,BITMAPINFOHEADER *bmp_ih, pixel_t *data)
{
unsigned int offset=sizeof(struct bmpfile_magic)
	+sizeof(struct bmpfile_header)
	+sizeof(BITMAPINFOHEADER)
	+(1<<bmp_ih->bitspp)*4;
struct bmpfile_header  bmp_fh={
	.filesz=offset+bmp_ih->bmp_bytesz,
	.creator1=0,
	.creator2=0,
	.bmp_offset=offset
};
FILE* fp=fopen(filename,"w");
if(fp==NULL)return 1;
struct bmpfile_magic mag={{0x42,0x4d}};
fwrite(&mag,1,sizeof(struct bmpfile_magic),fp);
fwrite(&bmp_fh,1,sizeof(struct bmpfile_header),fp);
fwrite(bmp_ih,1,sizeof(BITMAPINFOHEADER),fp);
// palette
rgb_t color={0,0,0,0};
int i;
for( i=0;i<(1<<bmp_ih->bitspp);i++)
{
	color.r=i;
	color.g=i;
	color.b=i;
fwrite(&color,1,sizeof(rgb_t),fp);
}
unsigned char c;
for(i=0;i<bmp_ih->bmp_bytesz;i++){
c=(unsigned char)data[i];
fwrite(&c,sizeof(unsigned char),1,fp);
}

fclose(fp);
return 0;
}
void convolution( pixel_t *in,pixel_t *out,float *kernel, int nx, int  ny, int kn,int norm)
{
int i,j,m,n,c;
 int khalf=floor(kn/2.);
float pixel, min=DBL_MAX,max=DBL_MIN;
if(norm)
for( m=khalf;m<nx-khalf;m++)
for( n=khalf;n<ny-khalf;n++)
{
pixel=0;
c=0;
for( j=-khalf;j<=khalf;j++)
for( i=-khalf;i<=khalf;i++)
pixel+=in[(n-j)*nx+m-i]*kernel[c++];

if(pixel<min)min=pixel;
if(pixel>max)max=pixel;
}
for( m=khalf;m<nx-khalf;m++)
for( n=khalf;n<ny-khalf;n++)
{
 pixel=0;
 c=0;
for( j=-khalf;j<=khalf;j++)
for( i=-khalf;i<=khalf;i++)
pixel+=in[(n-j)*nx+m-i]*kernel[c++];
if (norm)pixel=MAX_BRIGHTNESS*(pixel-min)/(max-min);
out[n*nx+m]=(pixel_t)pixel;
}
}
void gaussian_filter(pixel_t *in,pixel_t *out, int nx, int ny, float sigma)
{
int i,j,c=0;
const int n=2*(int)(2*sigma)+3;
float mean=(float)floor(n/2.);
float kernel[n*n];
fprintf(stderr,"gaussian_filter:kernel size %d,sigma=%g\n",n,sigma);
for( i=0;i<n;i++)
for( j=0;j<n;j++)
kernel[c++]=exp(-0.5*(pow((i-mean)/sigma,2.0)+pow((j-mean)/sigma,2.0)))/(2*M_PI*sigma*sigma);

convolution(in,out,kernel,nx,ny,n,1);
}
void  canny_edge_detection( pixel_t *in,pixel_t *out,int nx,int ny,  int tmin,const int tmax, float sigma)
{
int i,j,c,Gmax;
int counter=0;
float dir;
float Gx[]={
	-1,0,1,
	-2,0,2,
	-1,0,1};
float Gy[]={
	1,2,1,
	0,0,0,
	-1,-2,-1};
pixel_t *G =calloc(nx*ny*sizeof(pixel_t),1),
	 *after_Gx=calloc(nx*ny*sizeof(pixel_t),1),
	*after_Gy=calloc(nx*ny*sizeof(pixel_t),1),
	*nms=calloc(nx*ny*sizeof(pixel_t),1);

gaussian_filter(in,out,nx,ny,sigma);
convolution(out,after_Gx,Gx,nx,ny,3,0);
convolution(out,after_Gy,Gy,nx,ny,3,0);
for(i=1;i<nx-1;i++)
for( j=1;j<ny-1;j++)
{

 c=i+nx*j;
G[c]=hypot(after_Gx[c],after_Gy[c]);
if(G[c]>Gmax)Gmax=G[c];
}
//if(G[c]>Gmax) Gmax=G[c];
//}
int nn,ss,ww,ee,nw,ne,sw,se;
//non-maximum suppresion
for( i=1;i<nx-1;i++)
for( j=1;j<ny-1;j++)
{
 c=i+nx*j;
 nn=c-nx;
 ss=c+nx;
 ww=c+1;
 ee=c-1;
 nw=nn+1;
 ne=nn-1;
 sw=ss+1;
 se=ss-1;
dir=atan2(after_Gy[c],after_Gx[c])+8.0/M_PI;
if(((fabs(dir)<=1||fabs(fabs(dir)-8)<=1) && G[c]>G[ee]&&G[c]>G[ww])
||
((fabs(dir-2)<=1||fabs(dir+6)<=1)&& G[c]>G[nw]&&G[c]>G[se])
||
((fabs(fabs(fabs(dir)-4))<=1)&& G[c]>G[nn]&&G[c]>G[ss])
||
((fabs(dir-6)<=1||fabs(dir+2)<=1)&&G[c]>G[ne]&&G[c]>G[sw])
)

nms[c]=G[c];
else
nms[c]=0;
}
//Reuse array
pixel_t *edges=after_Gy;
memset(out,0,sizeof(pixel_t)*nx*ny);
memset(edges,0,sizeof(pixel_t)*nx*ny);
int nedges,k,t;
int nbs[8];
counter=0;

//Tracing edges with hysteresis.Non-recrsive implementation.

for(c=1,j=1;j<ny-1;j++)
for(i=1;i<nx-1;i++)
{
if(nms[c]>=tmax&&out[c]==0)
{
out[c]=MAX_BRIGHTNESS;
nedges=1;
edges[0]=c;
do{
nedges--;
t=edges[nedges];
nbs[0]=t-nx;
nbs[1]=t+nx;
nbs[2]=t+1;
nbs[3]=t-1;
nbs[4]=nbs[0]+1;
nbs[5]=nbs[0]-1;
nbs[6]=nbs[1]+1;
nbs[7]=nbs[1]-1;
for( k=0;k<8;k++)
if(nms[nbs[k]]>=tmin&&out[nbs[k]]==0)
{
out[nbs[k]]=MAX_BRIGHTNESS;
edges[nedges++]=nbs[k];
}
}while(nedges>0);
}
c++;
}
free(after_Gx);
free(after_Gy);
free(G);
free(nms);
}
void cross(pixel_t *in,BITMAPINFOHEADER *ih,int x,int y,int r)
{
int i, Nx=r,Ny=r;
int xmin=(x-Nx>0)? x-Nx:0,
xmax=(x+Nx<ih->width) ? x+Nx : ih->width-1,
ymin=(y-Ny>0) ? y-Ny:0,
ymax=(y+Ny<ih->height) ? y+Ny: ih->height-1;
for(i=xmin;i<=xmax;i++) in[i+ih->width*y]=255;
for(i=ymin;i<=ymax;i++) in[x+ih->width*i]=255;
}

void circular_hough_transform(pixel_t *in,BITMAPINFOHEADER *ih,int T1,int T2)
{
int nx=ih->width;
int ny=ih->height;
unsigned long long  i,j,k,l,x0,y0,
r2,
nx0=nx,
ny0=ny,
N=nx0*ny0,
nr2=50,
r2max=pow(fmax(nx,ny),2);
fprintf(stderr,"circular Hough space size %llu\n;r2max=%lld\n",nx*ny*nr2,r2max);
unsigned short *out=calloc(nx0*ny0*nr2,sizeof(unsigned short));
if(!out){
perror("circular_hough_transform()");
exit(1);
}
for(j=0;j<ny;j++)
for(i=0;i<nx;i++){
if(in[i+nx*j]<T1) continue;
for(y0=0;y0<ny0;y0++){
for(x0=0;x0<nx0;x0++){
r2=(i-x0)*(i-x0)+(j-y0)*(j-y0);
if(r2>r2max||r2<1) continue;
r2=floor(r2*nr2/r2max);
k=x0+nx0*y0+N*r2;
if(out[k]<USHRT_MAX)
out[k]++;
}
}
}
int counter=0;
//find maximum in the out array.
for(y0=0;y0<ny0;y0++)
for(x0=0;x0<nx0;x0++)
for(r2=0;r2<nr2;r2++)
if(out[x0+nx0*y0+N*r2]>T2)
{
printf("%d:\tx0=%d\ty0=%d\tr=%g\n",counter++,(int)x0,(int)y0,sqrt(r2*r2max/(double)nr2));
cross(in,ih,x0,y0,(int)sqrt(r2*r2max/(double)nr2));
}
free(out);
}
int main(int argc,char **argv)
{
if(argc<2){
printf("Usage:%s image.bmp\n",argv[0]);
exit( 1);
}

 pixel_t *bitmap_data,*temp_image;
bitmap_data=load_bmp(argv[1],&ih);
temp_image=(pixel_t *)malloc(ih.bmp_bytesz*sizeof(pixel_t));
printf("Info:%dx%dx%d\n",ih.width,ih.height,ih.bitspp);
canny_edge_detection(bitmap_data,temp_image,ih.width,ih.height,45,50,1);
circular_hough_transform(temp_image,&ih,200,220);
save_bmp("out.bmp",&ih,temp_image);
free(bitmap_data);
free(temp_image);
return 0;
}
...
i have made the changes as u have said .now out.bmp is created.but itz giving a full white blank image for 8 bit pic of 512x512 pixel.
is there any fault my program????
Last edited by jayshri; 03-20-2012 at 06:26 AM.
 
Old 03-20-2012, 08:19 AM   #92
firstfire
Member
 
Registered: Mar 2006
Location: Ekaterinburg, Russia
Distribution: Debian, Ubuntu
Posts: 709

Rep: Reputation: 428Reputation: 428Reputation: 428Reputation: 428Reputation: 428
Hough threshold T2 is too low.
 
Old 03-21-2012, 02:29 AM   #93
jayshri
Member
 
Registered: Feb 2012
Posts: 71

Rep: Reputation: Disabled
hi
first fire i have changed the treshold to 350 and 600..then also itz giving a white blank image with lena512.

actulay i am doing my project with a TEM image of nano particles of 8 bit and 512X512..

But when i use this image the loop doest start after a long time also..

in case of lena it start within 5 min but in Case of TEM it doen'st start.

here i am giving the link of a sample picture of tem. itz of dimension 600X600 and a jpg mage. but my image is of 512X512 and 8bit and a bmp image.
can u plz help me figure this problem out.
if the program run on the image i am attaching then it will run on my project image 2..
here is the link
thank u...
 
Old 03-21-2012, 02:41 AM   #94
jayshri
Member
 
Registered: Feb 2012
Posts: 71

Rep: Reputation: Disabled
and also will this program run with image 512x512 ???? because i have merged the program with canny edge ..so it accpets only 512x512 images..
when i have reduced it to 200x200 then itz giving segmentation fault..
here is my program
Code:
#include<stdint.h>
#include<stdio.h>
#include<stdlib.h>
#include<float.h>
#include<math.h>
#include<string.h>
#include<limits.h>
#define MAX_BRIGHTNESS 255
#define M_PI 3.14159265358979323846264338327
//#define OUTPUT_FILE "out.bmp"
struct bmpfile_magic{
unsigned char magic[2];
};
 struct bmpfile_header{
uint32_t filesz;
uint16_t creator1;
uint16_t creator2;
uint32_t  bmp_offset;
};
typedef struct{
uint32_t header_sz;
int32_t width;
int32_t height;
uint16_t nplanes;
uint16_t bitspp;
uint32_t compress_type;
uint32_t bmp_bytesz;
int32_t hres;
int32_t vres;
uint32_t ncolors;
uint32_t nimpcolors;
}BITMAPINFOHEADER;
typedef struct{
uint8_t r;
uint8_t g;
uint8_t b;
uint8_t null;
}rgb_t;
BITMAPINFOHEADER ih;


typedef int pixel_t;
pixel_t *load_bmp( char *filename,BITMAPINFOHEADER *bitmapInfoHeader)
{
  FILE *filePtr;
struct bmpfile_magic mag;
struct bmpfile_header bitmapFileHeader;
pixel_t *bitmapImage;
filePtr=fopen(filename,"r");
if(filePtr==NULL)
{
perror("fopen()");
exit(1);
}

fread(&mag,sizeof(struct bmpfile_magic),1,filePtr);
//verify that this is a bmp file by check bitmap id
if(*((uint16_t*) mag.magic)!=0x4D42)
{
fprintf(stderr,"Not a bmp file:magic=%c%c\n",mag.magic[0],mag.magic[1]);
fclose(filePtr);
return NULL;
}
//read the bitmap header
fread(&bitmapFileHeader,sizeof(struct bmpfile_header),1,filePtr);

//read the bitmap info header
fread(bitmapInfoHeader,sizeof(BITMAPINFOHEADER),1,filePtr);

if(bitmapInfoHeader->compress_type!=0)
fprintf(stderr,"Warning,compression is not supported.\n");
//move file point to the beginning of bitmap data
fseek(filePtr,bitmapFileHeader.bmp_offset,SEEK_SET);

//allocate enough memory for bitmao image data
bitmapImage=(pixel_t*)malloc(bitmapInfoHeader->bmp_bytesz*(sizeof(pixel_t)));
//verify memory allocation
if(!bitmapImage)
{
free(bitmapImage);
fclose(filePtr);
return NULL;
}

//read in the bitmap image data
int i;
unsigned char c;
for( i=0;i<bitmapInfoHeader->bmp_bytesz;i++){
fread(&c,sizeof(unsigned char),1,filePtr);
bitmapImage[i]=(int)c;
}
fclose(filePtr);
return bitmapImage;
}

//return:nonzero on error.
int save_bmp(char *filename,BITMAPINFOHEADER *bmp_ih, pixel_t *data)
{
unsigned int offset=sizeof(struct bmpfile_magic)
	+sizeof(struct bmpfile_header)
	+sizeof(BITMAPINFOHEADER)
	+(1<<bmp_ih->bitspp)*4;
struct bmpfile_header  bmp_fh={
	.filesz=offset+bmp_ih->bmp_bytesz,
	.creator1=0,
	.creator2=0,
	.bmp_offset=offset
};
FILE* fp=fopen(filename,"w");
if(fp==NULL)return 1;
struct bmpfile_magic mag={{0x42,0x4d}};
fwrite(&mag,1,sizeof(struct bmpfile_magic),fp);
fwrite(&bmp_fh,1,sizeof(struct bmpfile_header),fp);
fwrite(bmp_ih,1,sizeof(BITMAPINFOHEADER),fp);
// palette
rgb_t color={0,0,0,0};
int i;
for( i=0;i<(1<<bmp_ih->bitspp);i++)
{
	color.r=i;
	color.g=i;
	color.b=i;
fwrite(&color,1,sizeof(rgb_t),fp);
}
unsigned char c;
for(i=0;i<bmp_ih->bmp_bytesz;i++){
c=(unsigned char)data[i];
fwrite(&c,sizeof(unsigned char),1,fp);
}

fclose(fp);
return 0;
}
void convolution( pixel_t *in,pixel_t *out,float *kernel, int nx, int  ny, int kn,int norm)
{
int i,j,m,n,c;
 int khalf=floor(kn/2.);
float pixel, min=DBL_MAX,max=DBL_MIN;
if(norm)
for( m=khalf;m<nx-khalf;m++)
for( n=khalf;n<ny-khalf;n++)
{
pixel=0;
c=0;
for( j=-khalf;j<=khalf;j++)
for( i=-khalf;i<=khalf;i++)
pixel+=in[(n-j)*nx+m-i]*kernel[c++];

if(pixel<min)min=pixel;
if(pixel>max)max=pixel;
}
for( m=khalf;m<nx-khalf;m++)
for( n=khalf;n<ny-khalf;n++)
{
 pixel=0;
 c=0;
for( j=-khalf;j<=khalf;j++)
for( i=-khalf;i<=khalf;i++)
pixel+=in[(n-j)*nx+m-i]*kernel[c++];
if (norm)pixel=MAX_BRIGHTNESS*(pixel-min)/(max-min);
out[n*nx+m]=(pixel_t)pixel;
}
}
void gaussian_filter(pixel_t *in,pixel_t *out, int nx, int ny, float sigma)
{
int i,j,c=0;
const int n=2*(int)(2*sigma)+3;
float mean=(float)floor(n/2.);
float kernel[n*n];
fprintf(stderr,"gaussian_filter:kernel size %d,sigma=%g\n",n,sigma);
for( i=0;i<n;i++)
for( j=0;j<n;j++)
kernel[c++]=exp(-0.5*(pow((i-mean)/sigma,2.0)+pow((j-mean)/sigma,2.0)))/(2*M_PI*sigma*sigma);

convolution(in,out,kernel,nx,ny,n,1);
}
void  canny_edge_detection( pixel_t *in,pixel_t *out,int nx,int ny,  int tmin,const int tmax, float sigma)
{
int i,j,c,Gmax;
int counter=0;
float dir;
float Gx[]={
	-1,0,1,
	-2,0,2,
	-1,0,1};
float Gy[]={
	1,2,1,
	0,0,0,
	-1,-2,-1};
pixel_t *G =calloc(nx*ny*sizeof(pixel_t),1),
	 *after_Gx=calloc(nx*ny*sizeof(pixel_t),1),
	*after_Gy=calloc(nx*ny*sizeof(pixel_t),1),
	*nms=calloc(nx*ny*sizeof(pixel_t),1);

gaussian_filter(in,out,nx,ny,sigma);
convolution(out,after_Gx,Gx,nx,ny,3,0);
convolution(out,after_Gy,Gy,nx,ny,3,0);
for(i=1;i<nx-1;i++)
for( j=1;j<ny-1;j++)
{

 c=i+nx*j;
G[c]=hypot(after_Gx[c],after_Gy[c]);
if(G[c]>Gmax)Gmax=G[c];
}
//if(G[c]>Gmax) Gmax=G[c];
//}
int nn,ss,ww,ee,nw,ne,sw,se;
//non-maximum suppresion
for( i=1;i<nx-1;i++)
for( j=1;j<ny-1;j++)
{
 c=i+nx*j;
 nn=c-nx;
 ss=c+nx;
 ww=c+1;
 ee=c-1;
 nw=nn+1;
 ne=nn-1;
 sw=ss+1;
 se=ss-1;
dir=atan2(after_Gy[c],after_Gx[c])+8.0/M_PI;
if(((fabs(dir)<=1||fabs(fabs(dir)-8)<=1) && G[c]>G[ee]&&G[c]>G[ww])
||
((fabs(dir-2)<=1||fabs(dir+6)<=1)&& G[c]>G[nw]&&G[c]>G[se])
||
((fabs(fabs(fabs(dir)-4))<=1)&& G[c]>G[nn]&&G[c]>G[ss])
||
((fabs(dir-6)<=1||fabs(dir+2)<=1)&&G[c]>G[ne]&&G[c]>G[sw])
)

nms[c]=G[c];
else
nms[c]=0;
}
//Reuse array
pixel_t *edges=after_Gy;
memset(out,0,sizeof(pixel_t)*nx*ny);
memset(edges,0,sizeof(pixel_t)*nx*ny);
int nedges,k,t;
int nbs[8];
counter=0;

//Tracing edges with hysteresis.Non-recrsive implementation.

for(c=1,j=1;j<ny-1;j++)
for(i=1;i<nx-1;i++)
{
if(nms[c]>=tmax&&out[c]==0)
{
out[c]=MAX_BRIGHTNESS;
nedges=1;
edges[0]=c;
do{
nedges--;
t=edges[nedges];
nbs[0]=t-nx;
nbs[1]=t+nx;
nbs[2]=t+1;
nbs[3]=t-1;
nbs[4]=nbs[0]+1;
nbs[5]=nbs[0]-1;
nbs[6]=nbs[1]+1;
nbs[7]=nbs[1]-1;
for( k=0;k<8;k++)
if(nms[nbs[k]]>=tmin&&out[nbs[k]]==0)
{
out[nbs[k]]=MAX_BRIGHTNESS;
edges[nedges++]=nbs[k];
}
}while(nedges>0);
}
c++;
}
free(after_Gx);
free(after_Gy);
free(G);
free(nms);
}
void cross(pixel_t *in,BITMAPINFOHEADER *ih,int x,int y,int r)
{
int i, Nx=r,Ny=r;
int xmin=(x-Nx>0)? x-Nx:0,
xmax=(x+Nx<ih->width) ? x+Nx : ih->width-1,
ymin=(y-Ny>0) ? y-Ny:0,
ymax=(y+Ny<ih->height) ? y+Ny: ih->height-1;
for(i=xmin;i<=xmax;i++) in[i+ih->width*y]=255;
for(i=ymin;i<=ymax;i++) in[x+ih->width*i]=255;
}

void circular_hough_transform(pixel_t *in,BITMAPINFOHEADER *ih,int T1,int T2)
{
int nx=ih->width;
int ny=ih->height;
unsigned long long  i,j,k,l,x0,y0,
r2,
nx0=nx,
ny0=ny,
N=nx0*ny0,
nr2=50,
r2max=pow(fmax(nx,ny),2);
fprintf(stderr,"circular Hough space size %llu\n;r2max=%lld\n",nx*ny*nr2,r2max);
unsigned short *out=calloc(nx0*ny0*nr2,sizeof(unsigned short));
if(!out){
perror("circular_hough_transform()");
exit(1);
}
for(j=0;j<ny;j++)
for(i=0;i<nx;i++){
if(in[i+nx*j]<T1) continue;
for(y0=0;y0<ny0;y0++){
for(x0=0;x0<nx0;x0++){
r2=(i-x0)*(i-x0)+(j-y0)*(j-y0);
if(r2>r2max||r2<1) continue;
r2=floor(r2*nr2/r2max);
k=x0+nx0*y0+N*r2;
if(out[k]<USHRT_MAX)
out[k]++;
}
}
}
int counter=0;
//find maximum in the out array.
for(y0=0;y0<ny0;y0++)
for(x0=0;x0<nx0;x0++)
for(r2=0;r2<nr2;r2++)
if(out[x0+nx0*y0+N*r2]>T2)
{
printf("%d:\tx0=%d\ty0=%d\tr=%g\n",counter++,(int)x0,(int)y0,sqrt(r2*r2max/(double)nr2));
cross(in,ih,x0,y0,(int)sqrt(r2*r2max/(double)nr2));
}
free(out);
}
int main(int argc,char **argv)
{
if(argc<2){
printf("Usage:%s image.bmp\n",argv[0]);
exit( 1);
}

 pixel_t *bitmap_data,*temp_image;
bitmap_data=load_bmp(argv[1],&ih);
temp_image=(pixel_t *)malloc(ih.bmp_bytesz*sizeof(pixel_t));
printf("Info:%dx%dx%d\n",ih.width,ih.height,ih.bitspp);
canny_edge_detection(bitmap_data,temp_image,ih.width,ih.height,45,50,1);
circular_hough_transform(temp_image,&ih,200,220);
save_bmp("out.bmp",&ih,temp_image);
free(bitmap_data);
free(temp_image);
return 0;
}
.
thank u
 
Old 03-23-2012, 12:45 PM   #95
firstfire
Member
 
Registered: Mar 2006
Location: Ekaterinburg, Russia
Distribution: Debian, Ubuntu
Posts: 709

Rep: Reputation: 428Reputation: 428Reputation: 428Reputation: 428Reputation: 428
Hi.

I implemented an improved Hough algorithm for finding circles described here. It uses edge orientation information and is much faster than previous straightforward implementation. The search for maxima in the 3-dimensional accumulator array is not very simple. I implemented two types of search -- ordinary thresholding and recursive (though without direct use of recursion) search. The latter should reduce multiple matches. Neither is perfect.

Code:
typedef unsigned int stack_t;
#define sign(x) (x>0?1:-1)
void hough_circ_grad(pixel_t *orig, pixel_t *edge, BITMAPINFOHEADER *ih)
{
	int T;
	unsigned int nx = ih->width, ny = ih->height,
	    X, Y, R,	// coordinates in Hough space
	    i, j, k, max_votes = 0;
	unsigned int	mx = nx / 4,	// Hough space discretization
	    		my = ny / 4,
			mr = 40,
			rmax = 20;	// maximum radius of circles to be found, px.

	int Gx, Gy;	// Gradients
	float slope;
	float	hx = mx / (float)nx,
	    	hy = my / (float)ny,
	    	hr = mr / (float)rmax,
		dx, dy, r;
	unsigned int Ntotal = mx * my * mr, Nxy = mx * my;
	unsigned short *out = calloc(Ntotal, sizeof(unsigned short));
	printf("circular Hough space size %ux%ux%u = %u\n", mx, my, mr, Ntotal);
	if(!out){
		perror("circular_hough_transform()");
		exit(1);
	}
	puts("Hough: init.");
	for(j=1; j<ny-1; j++)
		for(i=1; i<nx-1; i++) {
			if(edge[i+nx*j] < 100) continue;
			Gx = orig[i+1 + nx*j] - orig[i-1+nx*j];
			Gy = orig[i + nx*(j+1)] - orig[i+nx*(j-1)];
			slope = tan( atan2(Gy, Gx) - M_PI/2. );
			for(R=0; R<mr; R++) {
				r = R/(float)hr;
				if(fabs(slope) < 1e-12)
					dx = 0;
				else
					dx =  -sign(Gx) * r / sqrt(1 + 1/slope/slope);
				dy = -sign(Gy) * r / sqrt(1 + slope*slope);
				X = (int)(i + dx) * hx;
				Y = (int)(j + dy) * hy;
				k = X + mx * Y + Nxy * R;
				if(out[k] < USHRT_MAX && k<Ntotal)
				{
					out[k]++;
					if(out[k]>max_votes) max_votes = out[k];
				}
			}
		}
	puts("Hough: accumulation done.");

	T = max_votes*0.5;	// search threshold
	printf("T=%d\n", T);

	unsigned int counter = 0;

	for(i=0; i<Ntotal; i++) if(out[i] < T) counter++;
	printf("zeros=%d of %d: %g\%\n", counter, Ntotal, counter/(float)Ntotal*100);

	int u,v,w, test,
	    sr = 1;	// recursive search radius

	counter = 0;

// #if 1 -- recursive maxima search; should reduce multiple matches for simgle circle.
// #if 0 -- thresholding.
#if 1
	unsigned int smaxsize = 1000,  ssize, k2, maxv, maxk, len;
	stack_t *stack = malloc(sizeof(stack_t)*smaxsize);
	stack_t se;

	for(Y=1; Y<my-1; Y++)
		for(X=1; X<mx-1; X++)
			for(R=3; R<mr-1; R++) {
				k = X + mx * Y + Nxy * R;
				if( out[k] < T ) continue;

				len = 0;
				maxv  = 0;
				ssize = 1;
				stack[0] = k;

				do {
					len++;
					ssize--;
					se = stack[ssize];
					if( out[se] > maxv ) {
						maxv = out[se];
						maxk = se;
					}
					out[se] = 0;	// we were here

					// scan all neighbours
					for(u=-sr;u<=sr;u++)
					for(v=-sr;v<=sr;v++)
					for(w=-sr;w<=sr;w++)
					{
						if(u==0 && v==0 && w==0) continue;
						k2 = se + u+v*mx+w*Nxy;
						if( k2 < Ntotal && out[k2] >= T ) {
							stack[ssize] = k2;
							ssize++;
							if(ssize >= smaxsize){
							  puts("NOTE: Incrementing stack size");
							  smaxsize += 1000;
							  stack = realloc(stack, sizeof(stack_t)*smaxsize);
							}
						}
					}
				} while(ssize > 0);

				//if( len < 3 ) continue;
				R = maxk / (float)Nxy;
				Y = (maxk % Nxy) / (float) mx;
				X = (maxk % Nxy) % mx;// maxk - mx*Y - Nxy*R;
				i = (int)(X/(float)hx);
				j = (int)(Y/(float)hy);
				r = (int)(R/(float)hr);
				printf("%d:\tx=%d\ty=%d\tr=%d (R=%d)\tvotes=%u, k=%d, m=%d\n",
						counter++, X, Y, (int)r, R, maxv, maxk, len);
				cross(edge, ih, i, j, r);

			}
	printf("max stack size = %d\n", smaxsize);
	free(stack);

#else
	for(Y=1; Y<my-1; Y++)
		for(X=1; X<mx-1; X++)
			for(R=5; R<mr-1; R++) {
				k = X + mx * Y + Nxy * R;
				if( out[k] < T /*T*(R+1)*/) continue;
				test = 0;
				
				for(u=-1;u<2;u++)
				for(v=-1;v<2;v++)
				for(w=-1;w<2;w++)
					if( !(u==0 && v==0 && w==0) )
						test += (out[k + u+v*mx+w*Nxy] > out[k]);
						
				if( test == 0 );
				{
					i = (int)(X/(float)hx);
					j = (int)(Y/(float)hy);
					r = (int)(R/(float)hr);
					printf("%d:\tx=%d\ty=%d\tr=%d (R=%d)\tvotes=%u\n",
							counter++, i, j, (int)r, R, out[X + mx * Y + Nxy * R]);
					cross(edge, ih, i, j, r);
				}
			}
#endif
	puts("Hough: search done.");
	printf("Hough: max number of votes %d\n", max_votes);
	free(out);
}
Adjustable parameters are in bold font.

Call this code as follows:
Code:
hough_circ_grad(bitmap_data, temp_image, &ih);
where `temp_image' is the result of Canny edge detection.

For tests I used the picture you provided. Here is the result
Click image for larger version Name: out.png Views: 34 Size: 35.2 KB ID: 9297

Hope that helps.
Last edited by firstfire; 03-23-2012 at 12:50 PM.
 
Old 03-23-2012, 08:41 PM   #96
jayshri
Member
 
Registered: Feb 2012
Posts: 71

Rep: Reputation: Disabled
hi Firstfire i have tried the code but while compilation its giving a warning : Implicit declaration of cross and undefined reference to cross.

how will i rectify it?

Thank You
 
Old 03-23-2012, 10:51 PM   #97
firstfire
Member
 
Registered: Mar 2006
Location: Ekaterinburg, Russia
Distribution: Debian, Ubuntu
Posts: 709

Rep: Reputation: 428Reputation: 428Reputation: 428Reputation: 428Reputation: 428
Hi.

Function cross was defined few posts back here.
 
Old 03-24-2012, 01:22 AM   #98
firstfire
Member
 
Registered: Mar 2006
Location: Ekaterinburg, Russia
Distribution: Debian, Ubuntu
Posts: 709

Rep: Reputation: 428Reputation: 428Reputation: 428Reputation: 428Reputation: 428
Hi.

Previously I relied on the direction of gradient when computing center of the circle, which is wrong. Only the gradient orientation should be used. Changes are in bold. Function `cross' draws gray pixels to not interfere with Hough algorithm.
Code:
#define GRAY_PIXEL	200
void cross(pixel_t *in, BITMAPINFOHEADER *ih, int x, int y, int r)
{
	int i, Nx, Ny;
	Nx = Ny = r;
	int xmin = (x-Nx>0) ? x-Nx : 0,
	    xmax = (x+Nx<ih->width) ? x+Nx : ih->width-1,
	    ymin = (y-Ny>0) ? y-Ny : 0,
	    ymax = (y+Ny<ih->height) ? y+Ny : ih->height-1;
	for(i=xmin; i<=xmax; i++) in[i+ih->width*y] = GRAY_PIXEL;
	for(i=ymin; i<=ymax; i++) in[x+ih->width*i] = GRAY_PIXEL;
}

typedef unsigned int stack_t;
#define sign(x) (x>0?1:-1)
void hough_circ_grad(pixel_t *orig, pixel_t *edge, BITMAPINFOHEADER *ih)
{
	int T;
	unsigned int nx = ih->width, ny = ih->height,
	    X, Y, R,	// coordinates in Hough space
	    i, j, k, max_votes = 0;
	unsigned int	mx = nx / 4,	// Hough space discretization
	    		my = ny / 4,
			mr = 80,
			rmax = 80;	// maximum radius of circles to be found, px.

	int Gx, Gy;	// Gradients
	float slope;
	float	hx = mx / (float)nx,
	    	hy = my / (float)ny,
	    	hr = mr / (float)rmax,
		dx, dy, r;
	unsigned int Ntotal = mx * my * mr, Nxy = mx * my;
	unsigned short *out = calloc(Ntotal, sizeof(unsigned short));
	printf("circular Hough space size %ux%ux%u = %u\n", mx, my, mr, Ntotal);
	if(!out){
		perror("circular_hough_transform()");
		exit(1);
	}
	puts("Hough: init.");
	for(j=1; j<ny-1; j++)
		for(i=1; i<nx-1; i++) {
			if(edge[i+nx*j] <= GRAY_PIXEL) continue;
			Gx = orig[i+1 + nx*j] - orig[i-1+nx*j];
			Gy = orig[i + nx*(j+1)] - orig[i+nx*(j-1)];
			slope = tan( atan2(Gy, Gx) - M_PI/2. );
			//printf("i=%d, j=%d, slope = %g\n", i, j, slope);
			//if(fabs(slope)<1e-3 ) cross(edge, ih, i, j, 2);
			int c=0;
			for(R=0; R<mr; R++) {
				r = R/(float)hr;
				if(fabs(slope) < 1e-12)
					dx = 0;
				else
					dx =  -sign(Gx) * r / sqrt(1 + 1/slope/slope);
				dy = -sign(Gy) * r / sqrt(1 + slope*slope);
				X = (int)(i + dx) * hx;
				Y = (int)(j + dy) * hy;
				
				k = X + mx * Y + Nxy * R;
				if(out[k] < USHRT_MAX && k<Ntotal)
				{
					out[k]++;
					X = (int)(i - dx) * hx;
				Y = (int)(j - dy) * hy;
				
				k = X + mx * Y + Nxy * R;
				if(out[k] < USHRT_MAX && k<Ntotal)
				{
					out[k]++;
					if(out[k]>max_votes) max_votes = out[k];
				}if(out[k]>max_votes) max_votes = out[k];
				}
                                
                                // we should not rely on gradient direction
				X = (int)(i - dx) * hx;
				Y = (int)(j - dy) * hy;
				
				k = X + mx * Y + Nxy * R;
				if(out[k] < USHRT_MAX && k<Ntotal)
				{
					out[k]++;
					if(out[k]>max_votes) max_votes = out[k];
				}
			}
		}
	puts("Hough: accumulation done.");

	T = max_votes*0.6;	// search threshold
	printf("T=%d\n", T);

	unsigned int counter = 0;

	for(i=0; i<Ntotal; i++) if(out[i] < T) counter++;
	printf("zeros=%d of %d: %g\%\n", counter, Ntotal, counter/(float)Ntotal*100);

	int u,v,w, test, T2,
	    sr = 2;	// recursive search radius

	counter = 0;

// #if 1 -- recursive maxima search; should reduce multiple matches for simgle circle.
// #if 0 -- thresholding.
#if 1
	unsigned int smaxsize = 1000,  ssize, k2, maxv, maxk, len;
	stack_t *stack = malloc(sizeof(stack_t)*smaxsize);
	stack_t se;

	for(Y=1; Y<my-1; Y++)
		for(X=1; X<mx-1; X++)
			for(R=3; R<mr-1; R++) {
				k = X + mx * Y + Nxy * R;
				if( out[k] < T ) continue;

				len = 0;
				maxv  = 0;
				ssize = 1;
				stack[0] = k;

				do {
					len++;
					ssize--;
					se = stack[ssize];
					if( out[se] > maxv ) {
						maxv = out[se];
						maxk = se;
					}
					out[se] = 0;	// we were here

					// scan all neighbours
					for(u=-sr;u<=sr;u++)
					for(v=-sr;v<=sr;v++)
					for(w=-sr;w<=sr;w++)
					{
						if(u==0 && v==0 && w==0) continue;
						k2 = se + u+v*mx+w*Nxy;
						if( k2 < Ntotal && out[k2] >= T ) {
							stack[ssize] = k2;
							ssize++;
							if(ssize >= smaxsize){
							  puts("NOTE: Incrementing stack size");
							  smaxsize += 1000;
							  stack = realloc(stack, sizeof(stack_t)*smaxsize);
							}
						}
					}
				} while(ssize > 0);

				//if( len < 3 ) continue;
				R = maxk / (float)Nxy;
				Y = (maxk % Nxy) / (float) mx;
				X = (maxk % Nxy) % mx;// maxk - mx*Y - Nxy*R;
				i = (int)(X/(float)hx);
				j = (int)(Y/(float)hy);
				r = (int)(R/(float)hr);
				printf("%d:\tx=%d\ty=%d\tr=%d (R=%d)\tvotes=%u, k=%d, m=%d\n",
						counter++, X, Y, (int)r, R, maxv, maxk, len);
				cross(edge, ih, i, j, r);

			}
	printf("max stack size = %d\n", smaxsize);
	free(stack);

#else
	for(Y=1; Y<my-1; Y++)
		for(X=1; X<mx-1; X++)
			for(R=5; R<mr-1; R++) {
				k = X + mx * Y + Nxy * R;
				if( out[k] < T) continue;
				test = 0;
				
				for(u=-1;u<2;u++)
				for(v=-1;v<2;v++)
				for(w=-1;w<2;w++)
					if( !(u==0 && v==0 && w==0) )
						test += (out[k + u+v*mx+w*Nxy] > out[k]);
						
				if( test == 0 );
				{
					i = (int)(X/(float)hx);
					j = (int)(Y/(float)hy);
					r = (int)(R/(float)hr);
					printf("%d:\tx=%d\ty=%d\tr=%d (R=%d)\tvotes=%u\n",
							counter++, i, j, (int)r, R, out[X + mx * Y + Nxy * R]);
					cross(edge, ih, i, j, r);
				}
			}
#endif
	puts("Hough: search done.");
	printf("Hough: max number of votes %d\n", max_votes);
	free(out);
}
 
Old 03-24-2012, 06:38 AM   #99
jayshri
Member
 
Registered: Feb 2012
Posts: 71

Rep: Reputation: Disabled
image processing.
hi firstfire,
just now i saw your recent posted code.. i wil modify my code. but i got a problem with previous code.. when i run my code loop doesn't start. at first i thought the code is a bit slow but then i wait for another 15 min but the out.bmp is not created..
i havent modiefied the code yet with ur recently posted code.. can u please tel me what's going wrong in my code that out.bmp is nt created???..because previously with ur program your out.bmp ws created.. but in my case it just showing threshold value ,hough value nd the display of radius doesn't start.
Code:
#include<stdint.h>
#include<stdio.h>`
#include<stdlib.h>
#include<float.h>
#include<math.h>
#include<string.h>
#include<limits.h>
#define MAX_BRIGHTNESS 255
#define M_PI 3.14159265358979323846264338327
//#define sign(x)(x>0?1:-1)
//#define OUTPUT_FILE "out.bmp"
struct bmpfile_magic{
unsigned char magic[2];
};
 struct bmpfile_header{
uint32_t filesz;
uint16_t creator1;
uint16_t creator2;
uint32_t  bmp_offset;
};
typedef struct{
uint32_t header_sz;
int32_t width;
int32_t height;
uint16_t nplanes;
uint16_t bitspp;
uint32_t compress_type;
uint32_t bmp_bytesz;
int32_t hres;
int32_t vres;
uint32_t ncolors;
uint32_t nimpcolors;
}BITMAPINFOHEADER;
typedef struct{
uint8_t r;
uint8_t g;
uint8_t b;
uint8_t null;
}rgb_t;
BITMAPINFOHEADER ih;


typedef int pixel_t;
pixel_t *load_bmp( char *filename,BITMAPINFOHEADER *bitmapInfoHeader)
{
  FILE *filePtr;
struct bmpfile_magic mag;
struct bmpfile_header bitmapFileHeader;
pixel_t *bitmapImage;
filePtr=fopen(filename,"r");
if(filePtr==NULL)
{
perror("fopen()");
exit(1);
}

fread(&mag,sizeof(struct bmpfile_magic),1,filePtr);
//verify that this is a bmp file by check bitmap id
if(*((uint16_t*) mag.magic)!=0x4D42)
{
fprintf(stderr,"Not a bmp file:magic=%c%c\n",mag.magic[0],mag.magic[1]);
fclose(filePtr);
return NULL;
}
//read the bitmap header
fread(&bitmapFileHeader,sizeof(struct bmpfile_header),1,filePtr);

//read the bitmap info header
fread(bitmapInfoHeader,sizeof(BITMAPINFOHEADER),1,filePtr);

if(bitmapInfoHeader->compress_type!=0)
fprintf(stderr,"Warning,compression is not supported.\n");
//move file point to the beginning of bitmap data
fseek(filePtr,bitmapFileHeader.bmp_offset,SEEK_SET);

//allocate enough memory for bitmao image data
bitmapImage=(pixel_t*)malloc(bitmapInfoHeader->bmp_bytesz*(sizeof(pixel_t)));
//verify memory allocation
if(!bitmapImage)
{
free(bitmapImage);
fclose(filePtr);
return NULL;
}

//read in the bitmap image data
int i;
unsigned char c;
for( i=0;i<bitmapInfoHeader->bmp_bytesz;i++){
fread(&c,sizeof(unsigned char),1,filePtr);
bitmapImage[i]=(int)c;
}
fclose(filePtr);
return bitmapImage;
}

//return:nonzero on error.
int save_bmp(char *filename,BITMAPINFOHEADER *bmp_ih, pixel_t *data)
{
unsigned int offset=sizeof(struct bmpfile_magic)
	+sizeof(struct bmpfile_header)
	+sizeof(BITMAPINFOHEADER)
	+(1<<bmp_ih->bitspp)*4;
struct bmpfile_header  bmp_fh={
	.filesz=offset+bmp_ih->bmp_bytesz,
	.creator1=0,
	.creator2=0,
	.bmp_offset=offset
};
FILE* fp=fopen(filename,"w");
if(fp==NULL)return 1;
struct bmpfile_magic mag={{0x42,0x4d}};
fwrite(&mag,1,sizeof(struct bmpfile_magic),fp);
fwrite(&bmp_fh,1,sizeof(struct bmpfile_header),fp);
fwrite(bmp_ih,1,sizeof(BITMAPINFOHEADER),fp);
// palette
rgb_t color={0,0,0,0};
int i;
for( i=0;i<(1<<bmp_ih->bitspp);i++)
{
	color.r=i;
	color.g=i;
	color.b=i;
fwrite(&color,1,sizeof(rgb_t),fp);
}
unsigned char c;
for(i=0;i<bmp_ih->bmp_bytesz;i++){
c=(unsigned char)data[i];
fwrite(&c,sizeof(unsigned char),1,fp);
}

fclose(fp);
return 0;
}
void convolution( pixel_t *in,pixel_t *out,float *kernel, int nx, int  ny, int kn,int norm)
{
int i,j,m,n,c;
 int khalf=floor(kn/2.);
float pixel, min=DBL_MAX,max=DBL_MIN;
if(norm)
for( m=khalf;m<nx-khalf;m++)
for( n=khalf;n<ny-khalf;n++)
{
pixel=0;
c=0;
for( j=-khalf;j<=khalf;j++)
for( i=-khalf;i<=khalf;i++)
pixel+=in[(n-j)*nx+m-i]*kernel[c++];

if(pixel<min)min=pixel;
if(pixel>max)max=pixel;
}
for( m=khalf;m<nx-khalf;m++)
for( n=khalf;n<ny-khalf;n++)
{
 pixel=0;
 c=0;
for( j=-khalf;j<=khalf;j++)
for( i=-khalf;i<=khalf;i++)
pixel+=in[(n-j)*nx+m-i]*kernel[c++];
if (norm)pixel=MAX_BRIGHTNESS*(pixel-min)/(max-min);
out[n*nx+m]=(pixel_t)pixel;
}
}
void gaussian_filter(pixel_t *in,pixel_t *out, int nx, int ny, float sigma)
{
int i,j,c=0;
const int n=2*(int)(2*sigma)+3;
float mean=(float)floor(n/2.);
float kernel[n*n];
fprintf(stderr,"gaussian_filter:kernel size %d,sigma=%g\n",n,sigma);
for( i=0;i<n;i++)
for( j=0;j<n;j++)
kernel[c++]=exp(-0.5*(pow((i-mean)/sigma,2.0)+pow((j-mean)/sigma,2.0)))/(2*M_PI*sigma*sigma);

convolution(in,out,kernel,nx,ny,n,1);
}
void  canny_edge_detection( pixel_t *in,pixel_t *out,int nx,int ny,  int tmin,const int tmax, float sigma)
{
int i,j,c,Gmax;
int counter=0;
float dir;
float Gx[]={
	-1,0,1,
	-2,0,2,
	-1,0,1};
float Gy[]={
	1,2,1,
	0,0,0,
	-1,-2,-1};
pixel_t *G =calloc(nx*ny*sizeof(pixel_t),1),
	 *after_Gx=calloc(nx*ny*sizeof(pixel_t),1),
	*after_Gy=calloc(nx*ny*sizeof(pixel_t),1),
	*nms=calloc(nx*ny*sizeof(pixel_t),1);

gaussian_filter(in,out,nx,ny,sigma);
convolution(out,after_Gx,Gx,nx,ny,3,0);
convolution(out,after_Gy,Gy,nx,ny,3,0);
for(i=1;i<nx-1;i++)
for( j=1;j<ny-1;j++)
{

 c=i+nx*j;
G[c]=hypot(after_Gx[c],after_Gy[c]);
if(G[c]>Gmax)Gmax=G[c];
}
//if(G[c]>Gmax) Gmax=G[c];
//}
int nn,ss,ww,ee,nw,ne,sw,se;
//non-maximum suppresion
for( i=1;i<nx-1;i++)
for( j=1;j<ny-1;j++)
{
 c=i+nx*j;
 nn=c-nx;
 ss=c+nx;
 ww=c+1;
 ee=c-1;
 nw=nn+1;
 ne=nn-1;
 sw=ss+1;
 se=ss-1;
dir=atan2(after_Gy[c],after_Gx[c])+8.0/M_PI;
if(((fabs(dir)<=1||fabs(fabs(dir)-8)<=1) && G[c]>G[ee]&&G[c]>G[ww])
||
((fabs(dir-2)<=1||fabs(dir+6)<=1)&& G[c]>G[nw]&&G[c]>G[se])
||
((fabs(fabs(fabs(dir)-4))<=1)&& G[c]>G[nn]&&G[c]>G[ss])
||
((fabs(dir-6)<=1||fabs(dir+2)<=1)&&G[c]>G[ne]&&G[c]>G[sw])
)

nms[c]=G[c];
else
nms[c]=0;
}
//Reuse array
pixel_t *edges=after_Gy;
memset(out,0,sizeof(pixel_t)*nx*ny);
memset(edges,0,sizeof(pixel_t)*nx*ny);
int nedges,k,t;
int nbs[8];
counter=0;

//Tracing edges with hysteresis.Non-recrsive implementation.

for(c=1,j=1;j<ny-1;j++)
for(i=1;i<nx-1;i++)
{
if(nms[c]>=tmax&&out[c]==0)
{
out[c]=MAX_BRIGHTNESS;
nedges=1;
edges[0]=c;
do{
nedges--;
t=edges[nedges];
nbs[0]=t-nx;
nbs[1]=t+nx;
nbs[2]=t+1;
nbs[3]=t-1;
nbs[4]=nbs[0]+1;
nbs[5]=nbs[0]-1;
nbs[6]=nbs[1]+1;
nbs[7]=nbs[1]-1;
for( k=0;k<8;k++)
if(nms[nbs[k]]>=tmin&&out[nbs[k]]==0)
{
out[nbs[k]]=MAX_BRIGHTNESS;
edges[nedges++]=nbs[k];
}
}while(nedges>0);
}
c++;
}
free(after_Gx);
free(after_Gy);
free(G);
free(nms);
}
typedef unsigned int stack_t;
#define sign(x) (x>0?1:-1)
void cross(pixel_t *in,BITMAPINFOHEADER *ih,int x,int y,int r)
{
int i,Nx=r,Ny=r;
int xmin=(x-Nx>0)? x-Nx:0,
xmax=(x+Nx<ih->width)? x+Nx :ih->width-1,
ymin=(y-Ny>0)? y-Ny:0,
ymax=(y+Ny<ih->height)?y+Ny:ih->height-1;
for(i=xmin;i<=xmax;i++) in[i+ih->width*y]=255;
for(i=ymin;i<=ymax;i++) in[x+ih->width*i]=255;
}
void hough_circ_grad(pixel_t *orig,pixel_t *edge,BITMAPINFOHEADER *ih)
{
int T;
unsigned int nx=ih->width,ny=ih->height,
X,Y,R, //coordinates in Hough Space
i,j,k,max_votes=0;
unsigned int mx=nx/4,
	     my=ny/4,
             mr=40,
             rmax=20; // maximum radius of circles to be found,px.
int Gx,Gy; // Gradients
float slope;
float hx=mx/(float)nx,
      hy=my/(float)ny,
      hr=mr/(float)rmax,
      dx,dy,r;
 unsigned int Ntotal=mx*my*mr,Nxy=mx*my;
 unsigned short *out=calloc(Ntotal,sizeof(unsigned short));
 printf("circular Hough space size %ux%ux%u=%u\n",mx,my,mr,Ntotal);
 if(!out){
   perror("circular_hough_transform()");
   exit(1);
}
puts("Hough:init.");
for(j=1;j<ny-1;j++)
for(i=1;i<nx-1;i++){
if(edge[i+nx*j]<100) continue;
Gx=orig[i+1+nx*j]-orig[i-1+nx*j];
Gy=orig[i+nx*(j+1)]-orig[i+nx*(j-1)];
slope=tan(atan2(Gy,Gx)-M_PI/2.);
for(R=0;R<mr;R++){
r=R/(float)hr;
if(fabs(slope)<1e-12)
dx=0;
else
dx=-sign(Gx) * r/sqrt(1+1/slope/slope);
dy=-sign(Gy) * r/sqrt(1+slope*slope);
X=(int)(i+dx)*hx;
Y=(int)(j+dy)*hy;
k=X+mx*Y+Nxy*R;
if(out[k]<USHRT_MAX && k<Ntotal)
{
out[k]++;
if(out[k]>max_votes)max_votes=out[k];
}
}
}
puts("Hough:accumulation done.");
T=max_votes*0.5; //search thresold
printf("T=%d\n",T);
unsigned int counter=0;
for(i=0;i<Ntotal;i++) if (out[i]<T) counter ++;
printf("zeros=%d of %d:%g\n",counter,Ntotal,counter/(float)Ntotal*100);//In this line if i use \% after %g as u hav written then itz giving a warning.
int u,v,w,test,
sr=1; //recursive search radius
counter=0;
//#if1---- recursive maxima search;should reduce multiple matches for single circle.
//#if 0----thresholding.
#if 1
unsigned int smaxsize=1000,ssize,k2,maxv,maxk,len;
stack_t *stack=malloc(sizeof(stack_t)*smaxsize);
stack_t se;
for(Y=1;Y<my-1;Y++)
for(X=1;X=mx-1;X++)
for(R=3;R<mr-1;R++){
k=X+mx*Y+Nxy*R;
if(out[k]<T)continue;
len=0;
maxv=0;
ssize=1;
stack[0]=k;
do{
len++;
ssize--;
se=stack[ssize];
if(out[se]>maxv){
maxv=out[se];
maxk=se;
}
out[se]=0;
// scan all neighbours
for(u=-sr;u<=sr;u++)
for(v=-sr;v<=sr;v++)
for(w=-sr;w<=sr;w++)
{
if(u==0 && v==0 && w==0)
continue;
k2=se+u+v*mx+w*Nxy;
if(k2<Ntotal && out[k2]>=T){
stack[ssize]=k2;
ssize++;
if(ssize>=smaxsize){
puts("NOTE: Incrementing stack size");
smaxsize+=1000;
stack=realloc(stack,sizeof(stack_t)*smaxsize);
}
}
}
}while(ssize>0);
//if(len<3)continue;
R=maxk/(float)Nxy;
Y=(maxk % Nxy)/(float) mx;
X=(maxk % Nxy) %mx;
i=(int)(X/(float)hx);
j=(int)(Y/(float)hy);
r=(int)(R/(float)hr);
printf("%d:\tx=%d\ty=%d\tr=%d(R=%d)\tvotes=%u,k=%d,m=%d\n",
counter++,X,Y,(int)r,R,maxv,maxk,len);
cross(edge,ih,i,j,r);
}
printf("max stack size=%d\n",smaxsize);
free(stack);
#else
for(Y=1;Y<my-1;Y++)
for(X=1;X=mx-1;X++)
for(R=5;R<mr-1;R++){
k=X+mx*Y+Nxy*R;
if(out[k]<T/*T*(R+1)*/)continue;
test=0;
for(u=-1;u<2;u++)
for(v=-1;v<2;v++)
for(w=-1;w<2;w++)
if(!(u==0 && v==0 && w==0))
test +=(out[k+u+v*mx+w*Nxy]>out[k]);
if(test==0);
{
i=(int)(X/(float)hx);
j=(int)(Y/(float)hy);
r=(int)(R/(float)hr);
printf("%d:\tx=%d\ty=%d\tr=%d (R=%d)\tvotes=%u\n",
counter++,i,j,(int)r,R,out[X+mx*Y+Nxy*R]);
cross(edge,ih,i,j,r);
}
}
#endif
puts("Hough:search done.");
printf("Hough:max number of votes %d\n",max_votes);
free(out);
}





int main(int argc,char **argv)
{
if(argc<2){
printf("Usage:%s image.bmp\n",argv[0]);
exit( 1);
}

 pixel_t *bitmap_data,*temp_image;
bitmap_data=load_bmp(argv[1],&ih);
temp_image=(pixel_t *)malloc(ih.bmp_bytesz*sizeof(pixel_t));
printf("Info:%dx%dx%d\n",ih.width,ih.height,ih.bitspp);
canny_edge_detection(bitmap_data,temp_image,ih.width,ih.height,45,50,1);
hough_circ_grad(bitmap_data,temp_image,&ih);
save_bmp("out.bmp",&ih,temp_image);
free(bitmap_data);
free(temp_image);
return 0;
}
the sample pic i am using is of 8 bit and of 512x512.will 512x512 run smoothly with this code???
Thanx a tonn 4 your help. i am indebt.
regards..
 
Old 03-24-2012, 07:39 AM   #100
firstfire
Member
 
Registered: Mar 2006
Location: Ekaterinburg, Russia
Distribution: Debian, Ubuntu
Posts: 709

Rep: Reputation: 428Reputation: 428Reputation: 428Reputation: 428Reputation: 428
Hi.

I don't know how you do it..
In this code:
Code:
#if 1
	unsigned int smaxsize=1000,ssize,k2,maxv,maxk,len;
	stack_t *stack=malloc(sizeof(stack_t)*smaxsize);
	stack_t se;
	for(Y=1;Y<my-1;Y++)
		for(X=1;X=mx-1;X++)      // BUG is here
			for(R=3;R<mr-1;R++){
				k=X+mx*Y+Nxy*R;
				if(out[k]<T)continue;
				len=0;
should be
Code:
	for(X=1;X<mx-1;X++)
Similarly, later in the code, there is the same error after `#else'.

The code should work on 512x512 image.
Last edited by firstfire; 03-24-2012 at 07:40 AM.
 
Old 03-25-2012, 10:46 AM   #101
jayshri
Member
 
Registered: Feb 2012
Posts: 71

Rep: Reputation: Disabled
hi first fire..
thax alot.. the code is runnig propely now..and it is so well xplained that i can understand it easily...
my prof is asking about a output something like this.
http://http://www.google.co.in/imgre...d65WgDg&zoom=1..
the code u have given me is very easy nd well xpland.. but will it be possible to draw the circles along wid dat code as shown in the pic of the link i have given above???
i mean where the cross is created i need to show that in a circle so that is is clearly visible...
thank u ....
 
Old 03-25-2012, 11:14 PM   #102
firstfire
Member
 
Registered: Mar 2006
Location: Ekaterinburg, Russia
Distribution: Debian, Ubuntu
Posts: 709

Rep: Reputation: 428Reputation: 428Reputation: 428Reputation: 428Reputation: 428
Hi.

Well, IMO crosses are more distinguishable than gray circles.

Anyway here is the code to draw circles, taken from java code on the page you provided:
Code:
void setPixel(pixel_t *in, BITMAPINFOHEADER *ih, int pix, int x, int y)
{
	int k = x + ih->width*y;
	if( k<0 || k>ih->width*ih->height ) return;
	in[x+ih->width*y] = pix;
}
void drawCircle(pixel_t *in, BITMAPINFOHEADER *ih, int xCenter, int yCenter, int r)
{
	cross(in, ih, xCenter, yCenter, r);
	int pix = GRAY_PIXEL;
	int x, y, r2;
	int radius = r;
	r2 = r * r;
	setPixel(in, ih, pix, xCenter, yCenter + radius);
	setPixel(in, ih, pix, xCenter, yCenter - radius);
	setPixel(in, ih, pix, xCenter + radius, yCenter);
	setPixel(in, ih, pix, xCenter - radius, yCenter);

	y = radius;
	x = 1;
	y = (int) (sqrt(r2 - 1) + 0.5);
	while (x < y) {
		setPixel(in, ih, pix, xCenter + x, yCenter + y);
		setPixel(in, ih, pix, xCenter + x, yCenter - y);
		setPixel(in, ih, pix, xCenter - x, yCenter + y);
		setPixel(in, ih, pix, xCenter - x, yCenter - y);
		setPixel(in, ih, pix, xCenter + y, yCenter + x);
		setPixel(in, ih, pix, xCenter + y, yCenter - x);
		setPixel(in, ih, pix, xCenter - y, yCenter + x);
		setPixel(in, ih, pix, xCenter - y, yCenter - x);
		x += 1;
		y = (int) (sqrt(r2 - x*x) + 0.5);
	}
	if (x == y) {
		setPixel(in, ih, pix, xCenter + x, yCenter + y);
		setPixel(in, ih, pix, xCenter + x, yCenter - y);
		setPixel(in, ih, pix, xCenter - x, yCenter + y);
		setPixel(in, ih, pix, xCenter - x, yCenter - y);
	}
}
Call it exactly as the `cross' function.
 
Old 03-26-2012, 06:11 AM   #103
jayshri
Member
 
Registered: Feb 2012
Posts: 71

Rep: Reputation: Disabled
i.m
hi,firstfire.....
i have written your codes and it is running fine.......thanks a lot...
 
Old 03-26-2012, 06:14 AM   #104
jayshri
Member
 
Registered: Feb 2012
Posts: 71

Rep: Reputation: Disabled
i.m
hi firstfire,
i want to work on adaptive weighted median filter..... ?what are the changes necessary to be made in the median filter code????Theoritcaly i got it.. itz a extra work suggested by professor..
thanks alot..
Last edited by jayshri; 03-26-2012 at 09:44 AM.
 
Old 03-28-2012, 05:46 AM   #105
jayshri
Member
 
Registered: Feb 2012
Posts: 71

Rep: Reputation: Disabled
im
hi,
should i call the drawcircle() in main function?
Last edited by jayshri; 03-31-2012 at 05:49 AM.
 
  


Reply
Page 7 of 9 « First 56 7 89



Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts
BB code is On
Smilies are On
[IMG] code is Off
HTML code is Off


Similar Threads
Thread Thread Starter Forum Replies Last Post
Applying median filter to a picture ... tomazN Programming 7 04-03-2006 05:15 AM
Minolta MC 2300dl + foo2zjs driver = color images printed in black? J_Szucs Linux - Hardware 1 02-06-2006 02:30 AM
What are the color code ? cigarstub Linux - Networking 5 10-22-2005 10:51 AM
how to draw color images with xlib? BBB Programming 19 07-15-2005 02:04 PM
filter the code Xiangbuilder Programming 7 10-23-2003 07:28 AM

LinuxQuestions.org > Forums > Non-*NIX Forums > Programming

All times are GMT -5. The time now is 05:48 PM.

Contact Us - Advertising Info - Rules - Privacy - LQ Merchandise - Donations - Contributing Member - LQ Sitemap -
Main Menu
Advertisement
My LQ
Write for LQ
LinuxQuestions.org is looking for people interested in writing Editorials, Articles, Reviews, and more. If you'd like to contribute content, let us know.
Main Menu
Syndicate
RSS1  Latest Threads
RSS1  LQ News
Twitter: @linuxquestions
Open Source Consulting | Domain Registration