thanx alot first fire

hey firstfire,
i have got the code from the same page u hav suggested for canny edge detection.
i want to perform circular hough transform of the image after canny
.
will this code work?????
thank u..

Hi.

IFAIK this code does linear Hough transform. Of course you can grasp some ideas from it, but it is not for circular transform. Main difference is that Hough space for linear transform is two(2)-dimentional, while Hough space for circular transform is three(3)-dimentional.

hi,
can u plz suggest any modification with this code so that i can convert i into circular.. i knw very little abt circular hough transform and also didn't able to find any program for it..

Hi.

The ugly code below, when called, tries to find circle in input image and draws a cross in its center:
It works only on small images, say 200x200 pixels, requires to specify sensible threshold values T1, T2, and is very slow. It crashes on lena512. I don't think the circular Hough transform is practical in this form.

I call it in the code like this
Hope that helps.

hi,
if i merge this code with the earlier program of canny edge detection will it work???

hi firstfire,
the program that u have posted for circular hough transform,i have tried it for 200*200 pixel,,,,bt its showing a warning
......hence its showing a segmentation fault at runtym.
reguards.

hi,
also if i run the program with a 8 bit image and of 512x512 dimeantion then itz not showing any segmentation fault.
itz givind the result
bt here itz nt showin any image.. i mean out.bmp is not created,,
thank u.

hi,
will dis code run 4 8 bit image?????

Yes, it will.

This warning has nothing to do with segmentation fault. You can replace `%d' by `%lld' to remove the warning. Or remove that line completely.

Are you sure the program finishes correctly for such dimensions? As I mentioned, lena512.bmp causes core dump (crashes) after a long delay on my machine and no out.bmp created. I did not tried to figure out the reasons.

hi,
ya no output bmpfile created..only showind this dimensions..so on which kind of picture should i test this program????

hi first fire,
heere is the program with canny edge detection and circular hough transform..bt when out.bmp created itz showing a blank image
#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 i, Nx=10,Ny=10;
for(i=x-Nx;i<=x+Nx;i++) in[i+ih->width*y]=255;
for(i=y-Ny;i<=y+Ny;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);
}
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,350);
save_bmp("out.bmp",&ih,temp_image);
free(bitmap_data);
free(temp_image);
return 0;
}




itz a 8 bit image bt when we reduce itz size to 200x200 then the program is giving a segmentation fault...
i am doing with a image of nano particles of 8bit and dimension 512X512
thank u

Hi.
Here is an example image (200x200x8) I use for tests:
circ.bmp.txt (remove the '.txt' suffix)
And here is an output produced by the above program with Hough threshold T2=220 (I modified function cross() so that the size of cross equals actual radius of circle):


I do not observe segmentation fault. If output image is blank, then try to decrease thresholds in Canny and Hough.

hi firstfire'
can u plz xplain me ur logic of modifying the function cross()........where should we change the value of canny and hough???.is it in the main part or any other part?..what are the modifications that r needed while using the image lena512.bmp,,bcoz our pic specification is the same to that of lena512..
reguards...

Hi.

There are two thresholds both in Canny and Hough.
You should first find appropriate pair of thresholds for Canny filter (without Hough to save time) so that output image contains edges you are interested in. Then you should try different values of second threshold in Hough transform (220 in above code snippet) so that circles you are interested in are detected. This number equals to the minimum number of pixels on a circle for the circle to be detected.

As for the segmentation fault, I found the source of all evil -- it is a cross() function, I did not checked for array boundary violation. Here is a new version
Call it as follows: