region_image.c

it runs for blob detection with opencv

/*
**  Source file for blobdetect. (c) 2004 Per-Erik Forssen
**
**  This program is free software; you can redistribute it and/or modify
**  it under the terms of the GNU General Public License as published by
**  the Free Software Foundation; either version 2 of the License, or
**  (at your option) any later version.
**  
**  See the file COPYING for details.
**
*/
#include "region_image.h"
#include <stdio.h>
#include <math.h>
/*
** Region propagation without counting
**
**  bf_labelim1  Input label image (YxX)
**  bf_labelim2  Output label image (YxX)
**  bf_imk       Property image (YxXxD)
**  bf_imc       Confidence image (YxX)
**  bf_pvec      Prototype list (DxN)
**  dmax2        Squared max property distance
** Returns number of new seeds (for later allocation by find_new_seeds)
*/

int propagate_regions(ibuffer *bf_labelim1,ibuffer *bf_labelim2,
		      buffer *bf_imk,ibuffer *bf_imc,buffer *bf_pvec,
		      fpnum dmax2)
{
int x,y,xx,yy,xx2,yy2,k,new_seeds;
int r_out,c_out,r_in,c_in,ndim,l1,l2,l3;
int *in_labels,*out_labels,*imcdata;
fpnum *imkdata,*pvecdata;
fpnum diff,cp,d1,d2,d3;
int new_label;

 in_labels=bf_labelim1->data;
 r_in=bf_labelim1->rows;
 c_in=bf_labelim1->cols;

 out_labels=bf_labelim2->data;

 imkdata = bf_imk->data;
 r_out   = bf_imk->rows;
 c_out   = bf_imk->cols;
 ndim    = bf_imk->ndim;
 imcdata = bf_imc->data;

 pvecdata= bf_pvec->data;

 new_seeds=0;

 for(x=0;x<c_out;x++) {
   for(y=0;y<r_out;y++) {
     if(imcdata[y+x*r_out]>0) {
       /* This is a homogeneous pixel */

       /* Compute coordinates at input scale */
       xx=(int)floor(x/2);
       yy=(int)floor(y/2);
       if((x-2*xx)==0) {
	 xx2=xx-1;if(xx2<0) xx2=0;
       } else {
	 xx2=xx+1;if(xx2>=c_in) xx2=c_in-1;
       }
       if((y-2*yy)==0) {
	 yy2=yy-1;if(yy2<0) yy2=0;
       } else {
	 yy2=yy+1;if(yy2>=r_in) yy2=r_in-1;
       }
       /* Read labels */
       l1=in_labels[yy+xx*r_in];
       l2=in_labels[yy2+xx*r_in];
       l3=in_labels[yy+xx2*r_in];

       /* Compute distances */
       d1=0;d2=0;d3=0;
       for(k=0;k<ndim;k++) {
	 cp=imkdata[y+x*r_out+k*r_out*c_out];
	 if(l1>0) {diff=cp-pvecdata[k+(l1-1)*ndim];d1+=diff*diff;}
	 if(l2>0) {diff=cp-pvecdata[k+(l2-1)*ndim];d2+=diff*diff;}
	 if(l3>0) {diff=cp-pvecdata[k+(l3-1)*ndim];d3+=diff*diff;}
       }

       /* Decide label */
       new_label=0;
       if((d1<dmax2)&&(l1>0)&&((d1<=d2)||(l2==0))&&((d1<=d3)||(l3==0))) {
	 new_label=in_labels[yy+xx*r_in];
       }
       if((d2<dmax2)&&(l2>0)&&((d2<=d1)||(l1==0))&&((d2<=d3)||(l3==0))) {
	 new_label=in_labels[yy2+xx*r_in];
       }
       if((d3<dmax2)&&(l3>0)&&((d3<=d1)||(l1==0))&&((d3<=d2)||(l2==0))) {
	 new_label=in_labels[yy+xx2*r_in];
       }
       out_labels[y+x*r_out]=new_label;
       if(new_label==0) new_seeds++;
     }
   }
 }
 return(new_seeds);
}
/*
** Assign new seeds and fill in prototypes in bf_pvec2 array.
**
**  bf_labelim  Input label image (YxX)
**  bf_imk      Property image (YxXxD)
**  bf_imc      Confidence image (YxX)
**  bf_pvec1    Input prototype list (Dx<regions>)
**  bf_pvec2    Output prototype list (Dx<regions>+<new_seeds>)
**  regions     Length of pvec1
**  new_seeds   Number of new seeds (as found by propagate_regions)
**
*/
void find_new_seeds(ibuffer *bf_labelim,buffer *bf_imk,
		    ibuffer *bf_imc,buffer *bf_pvec1,
		    buffer *bf_pvec2,int regions,int new_seeds)
{
  int rows,cols,ndim,x,y,d,cind,nlabel;
  int *out_labelim,*imcdata;
  fpnum *imkdata,*pvec1data,*pvec2data;

  out_labelim=bf_labelim->data;
  imkdata = bf_imk->data;
  rows    = bf_imk->rows;
  cols    = bf_imk->cols;
  ndim    = bf_imk->ndim;
  imcdata = bf_imc->data;
  pvec1data = bf_pvec1->data;
  pvec2data = bf_pvec2->data;

  /* First copy old pvec */
  for(x=0;x<regions;x++) {
    for(d=0;d<ndim;d++) {
      pvec2data[x*ndim+d]=pvec1data[x*ndim+d];
    }
  }
  /* Now assign new labels and read prototype values */
  nlabel=regions+1; /* next label */
  for(x=0;x<cols;x++) {
    for(y=0;y<rows;y++) {
      cind=y+x*rows;
      if((imcdata[cind]>0)&&(out_labelim[cind]==0)) {
	/* This is a new seed */
	out_labelim[cind]=nlabel;
	/* Copy prototype value */
	for(d=0;d<ndim;d++) {
	  pvec2data[(nlabel-1)*ndim+d]=imkdata[cind+d*rows*cols];
	}
	nlabel++;
	if(nlabel>(regions+new_seeds+1))
	  printf("OOPS!\n");
      }      
    }
  }
}
/*
** Region propagation with counting
**
**  bf_labelim1  Input label image ([Y/2]x[X/2])
**  bf_labelim2  Output label image (YxX)
**  bf_imk       Property image (YxXxD)
**  bf_imc       Confidence image (YxX)
**  bf_pvec      Prototype list (DxN)
**  bf_cntl      boundary count list  
**  dmax2        Squared max property distance
**
*/

void propagate_regions_cnt(ibuffer *bf_labelim1,ibuffer *bf_labelim2,
			   buffer *bf_imk,ibuffer *bf_imc,buffer *bf_pvec,
			   ibuffer *bf_cntl,fpnum dmax2)
{
int x,y,xx,yy,xx2,yy2,k,ll,hl;
int r_out,c_out,r_in,c_in,ndim,l1,l2,l3;
fpnum diff,cp,dmin,d1,d2,d3,dl12,dl13,dl23;
int new_label;
int *in_labels,*out_labels,*imcdata;
fpnum *imkdata,*pvecdata;
int *cntl;

 in_labels=bf_labelim1->data;
 r_in=bf_labelim1->rows;
 c_in=bf_labelim1->cols;

 out_labels=bf_labelim2->data;

 imkdata = bf_imk->data;
 r_out   = bf_imk->rows;
 c_out   = bf_imk->cols;
 ndim    = bf_imk->ndim;
 imcdata = bf_imc->data;

 pvecdata= bf_pvec->data;
 cntl    = bf_cntl->data;

 for(x=0;x<c_out;x++) {
   for(y=0;y<r_out;y++) {
     if(imcdata[y+x*r_out]>0) {
       /* This is a homogeneous pixel */

       /* Compute coordinates at input scale */
       xx=(int)floor(x/2);
       yy=(int)floor(y/2);
       if((x-2*xx)==0) {
	 xx2=xx-1;if(xx2<0) xx2=0;
       } else {
	 xx2=xx+1;if(xx2>=c_in) xx2=c_in-1;
       }
       if((y-2*yy)==0) {
	 yy2=yy-1;if(yy2<0) yy2=0;
       } else {
	 yy2=yy+1;if(yy2>=r_in) yy2=r_in-1;
       }
       /* Read labels */
       l1=in_labels[yy+xx*r_in];
       l2=in_labels[yy2+xx*r_in];
       l3=in_labels[yy+xx2*r_in];

       /* Compute distances */
       d1=0;d2=0;d3=0;dl12=0;dl13=0;dl23=0;
       for(k=0;k<ndim;k++) {
	 cp=imkdata[y+x*r_out+k*r_out*c_out];
	 if(l1>0) {
	   diff=cp-pvecdata[k+(l1-1)*ndim];d1+=diff*diff;
	   if(l2>0) {
	     diff=pvecdata[k+(l1-1)*ndim]-pvecdata[k+(l2-1)*ndim];
	     dl12+=diff*diff;
	   }
	   if(l3>0) {
	     diff=pvecdata[k+(l1-1)*ndim]-pvecdata[k+(l3-1)*ndim];
	     dl13+=diff*diff;
	   }
	 }
	 if(l2>0) {
	   diff=cp-pvecdata[k+(l2-1)*ndim];d2+=diff*diff;
	   if(l3>0) {
	     diff=pvecdata[k+(l2-1)*ndim]-pvecdata[k+(l3-1)*ndim];
	     dl23+=diff*diff;
	   }
	 }
	 if(l3>0) {diff=cp-pvecdata[k+(l3-1)*ndim];d3+=diff*diff;}
       }
       
       /* Decide label */
       new_label=0;dmin=dmax2;

       if((l1>0)&&(d1<dmin)) {new_label=l1;dmin=d1;}
       if((l2>0)&&(d2<dmin)) {new_label=l2;dmin=d2;}
       if((l3>0)&&(d3<dmin)) {new_label=l3;dmin=d3;}

       out_labels[y+x*r_out]=new_label;

       /* Count boundary pixels */
       if(new_label>0) {
	 if(l1>0) {
	   if((l2>0)&&(l1!=l2)) if(dl12<dmax2) {
	     if(l1<l2) {ll=l1;hl=l2;} else {ll=l2;hl=l1;}
	     cntl[(hl-1)*(hl-2)/2+ll-1]+=1;
	   }	   
	   if((l3>0)&&(l1!=l3)) if(dl13<dmax2) {
	     if(l1<l3) {ll=l1;hl=l3;} else {ll=l3;hl=l1;}
	     cntl[(hl-1)*(hl-2)/2+ll-1]+=1;
	   }	
	 }
	 if(l2>0) if((l3>0)&&(l2!=l3)) if(dl23<dmax2) {
	   if(l2<l3) {ll=l2;hl=l3;} else {ll=l3;hl=l2;}
	   cntl[(hl-1)*(hl-2)/2+ll-1]+=1;
	 }
       }
     }
   }
 }
}
/*
** Region expansion within same scale, with counting
**
**  bf_labelim0  Input label image (YxX)
**  bf_labelim1  Output label image (YxX)
**  bf_imk       Property image (YxXxD)
**  bf_imc       Confidence image (YxX)
**  bf_pvec      Prototype list (DxN)
**  bf_cntl      boundary count list  
**  dmax2        Squared max property distance
**
*/
void expand_regions_cnt(ibuffer *bf_labelim0,ibuffer *bf_labelim1,
			buffer *bf_imk,ibuffer *bf_imc,buffer *bf_pvec,
			ibuffer *bf_cntl,fpnum dmax2)
{
int x,y,k,cind,new_label;
int rows,cols,pixels,ndim,ll,lh;
int *imc,*labelim0,*labelim1,*cntl;
int l1,l2,l3,l4; /* Four labels */
fpnum diff,d1,d2,d3,d4,mindist;
fpnum d12,d13,d14,d23,d24,d34;
fpnum *imk,*pvec;
 rows=bf_labelim0->rows;
 cols=bf_labelim0->cols;
 pixels=rows*cols;
 ndim=bf_imk->ndim;
 labelim0=bf_labelim0->data;
 labelim1=bf_labelim1->data;
 imc=bf_imc->data;
 imk=bf_imk->data;
 pvec=bf_pvec->data;
 cntl=bf_cntl->data;

 for(x=0;x<cols;x++) {
   for(y=0;y<rows;y++) { 
     cind=y+x*rows;
     labelim1[cind]=labelim0[cind];
     if(labelim0[cind]==0) {
       if(imc[cind]>0) {
	 /* now check if a neighbour is within dmax */
	 l1=l2=l3=l4=0;
	 d1=d2=d3=d4=0;
	 if(y>0)      l1=labelim0[cind-1];    	 /* Check above */
	 if(y<rows-1) l2=labelim0[cind+1];	 /* Check below */
	 if(x>0)      l3=labelim0[cind-rows]; 	 /* Check left  */
	 if(x<cols-1) l4=labelim0[cind+rows];	 /* Check right */

	 for(k=0;k<ndim;k++) {
	   if(l1>0) {
	     diff=imk[cind+k*pixels]-pvec[(l1-1)*ndim+k];
	     d1+=diff*diff;	     
	   }
	   if(l2>0) {
	     diff=imk[cind+k*pixels]-pvec[(l2-1)*ndim+k];
	     d2+=diff*diff;	     
	   }
	   if(l3>0) {
	     diff=imk[cind+k*pixels]-pvec[(l3-1)*ndim+k];
	     d3+=diff*diff;	     
	   }
	   if(l4>0) {
	     diff=imk[cind+k*pixels]-pvec[(l4-1)*ndim+k];
	     d4+=diff*diff;	     
	   }
	 }
	 /* Decide label */
	 new_label=0;mindist=dmax2;
	 if((l1>0)&&(d1<mindist)) {new_label=l1;mindist=d1;}
	 if((l2>0)&&(d2<mindist)) {new_label=l2;mindist=d2;}
	 if((l3>0)&&(d3<mindist)) {new_label=l3;mindist=d3;}
	 if((l4>0)&&(d4<mindist)) {new_label=l4;mindist=d4;}

	 labelim1[cind]=new_label;

	 if(new_label>0) {
	   /* Count boundary pixels */
	   d12=d13=d14=d23=d24=d34=0;
	   for(k=0;k<ndim;k++) {
	     if(l1>0) {
	       if(l2>0) {
		 diff=pvec[(l1-1)*ndim+k]-pvec[(l2-1)*ndim+k];
		 d12+=diff*diff;
	       }
	       if(l3>0) {
		 diff=pvec[(l1-1)*ndim+k]-pvec[(l3-1)*ndim+k];
		 d13+=diff*diff;
	       }
	       if(l4>0) {
		 diff=pvec[(l1-1)*ndim+k]-pvec[(l4-1)*ndim+k];
		 d14+=diff*diff;
	       }
	     }
	     if(l2>0) {
	       if(l3>0) {
		 diff=pvec[(l2-1)*ndim+k]-pvec[(l3-1)*ndim+k];
		 d23+=diff*diff;
	       }
	       if(l4>0) {
		 diff=pvec[(l2-1)*ndim+k]-pvec[(l4-1)*ndim+k];
		 d24+=diff*diff;
	       }
	     }
	     if(l3>0) {
	       if(l4>0) {
		 diff=pvec[(l3-1)*ndim+k]-pvec[(l4-1)*ndim+k];
		 d34+=diff*diff;
	       }
	     }
	   }
	   /* Check 1-2 */
	   if((l1>0)&&(l2>0)&&(l1!=l2)&&(d12<dmax2)) {
	     if(l1<l2) {ll=l1;lh=l2;} else {ll=l2;lh=l1;}
	     cntl[(lh-1)*(lh-2)/2+ll-1]+=1;
	   }
	   /* Check 1-3 */
	   if((l1>0)&&(l3>0)&&(l1!=l3)&&(d13<dmax2)) {
	     if(l1<l3) {ll=l1;lh=l3;} else {ll=l3;lh=l1;}
	     cntl[(lh-1)*(lh-2)/2+ll-1]+=1;
	   }
	   /* Check 1-4 */
	   if((l1>0)&&(l4>0)&&(l1!=l4)&&(d14<dmax2)) {
	     if(l1<l4) {ll=l1;lh=l4;} else {ll=l4;lh=l1;}
	     cntl[(lh-1)*(lh-2)/2+ll-1]+=1;
	   }
	   /* Check 2-3 */
	   if((l2>0)&&(l3>0)&&(l2!=l3)&&(d23<dmax2)) {
	     if(l2<l3) {ll=l2;lh=l3;} else {ll=l3;lh=l2;}
	     cntl[(lh-1)*(lh-2)/2+ll-1]+=1;
	   }
	   /* Check 2-4 */
	   if((l2>0)&&(l4>0)&&(l2!=l4)&&(d24<dmax2)) {
	     if(l2<l4) {ll=l2;lh=l4;} else {ll=l4;lh=l2;}
	     cntl[(lh-1)*(lh-2)/2+ll-1]+=1;
	   }
	   /* Check 3-4 */
	   if((l3>0)&&(l4>0)&&(l3!=l4)&&(d34<dmax2)) {
	     if(l3<l4) {ll=l3;lh=l4;} else {ll=l4;lh=l3;}
	     cntl[(lh-1)*(lh-2)/2+ll-1]+=1;
	   }
	 }
       }
     }
   }
 }
 printf("expand_regions\n");
}
/*
** Moment computation on an image
** Upper left pixel is [1,1]
**
**  bf_labelim  Label image (YxX)
**  bf_image    Image (YxXxN)
**  bf_pvec     output property averages
**  bf_mvec     output moments
**
*/
void compute_moments(ibuffer *bf_labelim,buffer *bf_image,
		     buffer *bf_pvec,buffer *bf_mvec)
{
int x,y,k,rind;
int rows,cols,ndim,label,blobs;
fpnum a,mx,my;
int *labelim;
fpnum *image,*out_pvec,*out_mvec;

 labelim = bf_labelim->data; 
 image   = bf_image->data;
 rows    = bf_image->rows;
 cols    = bf_image->cols;
 ndim    = bf_image->ndim;
 
 out_mvec = bf_mvec->data;
 out_pvec = bf_pvec->data;
 blobs    = bf_mvec->cols;
 
 /* Compute raw moments */
 for(x=0;x<cols;x++) {
   for(y=0;y<rows;y++) {
     rind=y+x*rows;
     label=labelim[rind];
     if(label>0) {
       label--;
       out_mvec[label*6]+=1;             /* 0:th moment  */
       out_mvec[label*6+1]+=(x+1);       /* 1:st moments */
       out_mvec[label*6+2]+=(y+1);       /* 1:st moments */
       out_mvec[label*6+3]+=(x+1)*(x+1); /* 2:nd moments */
       out_mvec[label*6+4]+=(x+1)*(y+1); /* 2:nd moments */
       out_mvec[label*6+5]+=(y+1)*(y+1); /* 2:nd moments */
       for(k=0;k<ndim;k++) {
	 out_pvec[label*ndim+k]+=image[rind+k*rows*cols];
       }
     }
   }
 }

 /* Convert raw moments to centroid, inertia etc. */
 for(label=0;label<blobs;label++) {
   a=out_mvec[label*6];
   if(a>0) {
     out_mvec[label*6+1]/=a;mx=out_mvec[label*6+1];
     out_mvec[label*6+2]/=a;my=out_mvec[label*6+2];
     out_mvec[label*6+3]=out_mvec[label*6+3]/a-mx*mx;
     out_mvec[label*6+4]=out_mvec[label*6+4]/a-mx*my;
     out_mvec[label*6+5]=out_mvec[label*6+5]/a-my*my;
     for(k=0;k<ndim;k++) {
       out_pvec[label*ndim+k]/=a;
     }
   }
 }
}
/*
** Loop over label image and replace old labels
** with new that are compatible with mvec and pvec
** lists. This is done according to out_ind list.
** The list out_ind is also simplified.
**
**  bf_labelim  Label image to modify (YxX)
**  bf_out_ind  Compaction list (1xN)
** 
*/
void labelim_compact(ibuffer *bf_labelim,ibuffer *bf_out_ind) {
int x,y,l,rows,cols;
int nblobs,cind,cl,rind;
int *out_ind,*labelim,*out_ind2;
ibuffer *bf_out_ind2;

 rows=bf_labelim->rows;
 cols=bf_labelim->cols;
 labelim=bf_labelim->data;

 nblobs=bf_out_ind->cols;
 out_ind=bf_out_ind->data;

 /* Simplify compaction list */
 for(l=0;l<nblobs;l++) {
   cl=out_ind[l]-1;
   while(cl>=0) {
     if(cl==out_ind[cl]-1) break;
     cl=out_ind[cl]-1;
   }
   out_ind[l]=cl+1;
 }

 /* Make a copy of out_ind with lower numbers */
 bf_out_ind2=ibuffer_new(1,nblobs,1);
 out_ind2=bf_out_ind2->data;
 rind=1;
 for(l=0;l<nblobs;l++) {
   cl=out_ind[l]-1;              /* Get prototype label for this region */
   if(cl>=0) {
     if(out_ind2[cl]>0) {
       out_ind2[l]=out_ind2[cl]; /* Copy new label from prototype */
     } else {
       out_ind2[l]=rind++;       /* This is prototype */
     }
   }
 }

 /* Update labels in label image */
 for(x=0;x<cols;x++) {
   for(y=0;y<rows;y++) {
     cind=x*rows+y;
     cl=labelim[cind];
     if(cl>0) labelim[cind]=out_ind2[cl-1];
   }
 }
 ibuffer_free(bf_out_ind2);
}