extract_blobs.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 <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include "extract_blobs.h"
#include "pnmio.h"
#include "sbinfilt.h"
#include "region_image.h"
#include "merge_blobs.h"
#include "file_and_time.h"

/*
 *   Calculate number of scales
 *   required to build an octave pyramid
 *   of bf_image.
 */
int number_of_scales(buffer *bf_image) {
  int maxdim;
  maxdim=bf_image->rows;
  if(maxdim<bf_image->cols)
    maxdim=bf_image->cols;
  return((int)ceil(log(maxdim)/log(2))+1);
}

/*
 *   Create a pyramid and insert
 *   input image at scale 0.
 */
buffer **imk_pyramid_new(buffer *bf_image) {
  int nsc,k,rows,cols,ndim;
  buffer **bf_list;

  /* Calculate number of scales */
  nsc=number_of_scales(bf_image);
#ifdef DEBUG
  printf("Scales %d\n",nsc);
#endif
  /* Allocate array */
  bf_list=(buffer **)calloc(nsc,sizeof(buffer *));
  bf_list[0]=bf_image;    /* Insert image at pos 0 */
  rows=bf_image->rows;
  cols=bf_image->cols;
  ndim=bf_image->ndim;
  for(k=1;k<nsc;k++) {
    rows=(rows+1)/2;  /* Half size output */
    cols=(cols+1)/2;   
    bf_list[k]=buffer_new(rows,cols,ndim);
#ifdef DEBUG
    printf("scale %d ",k);
    buffer_pdims(bf_list[k]);   
#endif
  }
  return(bf_list);
}

/*
 *   Create a certainty pyramid
 *   and insert input array at scale 0
 */
ibuffer **imc_pyramid_new(ibuffer *bf_imc) {
int nsc,k,rows,cols;
ibuffer **bf_list;
 
 rows=bf_imc->rows; 
 cols=bf_imc->cols; 
 /* Calculate number of scales */
 nsc=number_of_scales((buffer *)bf_imc);
#ifdef DEBUG
 printf("Scales %d\n",nsc);
#endif
 /* Allocate array */
 bf_list=(ibuffer **)calloc(nsc,sizeof(ibuffer *));
 bf_list[0]=bf_imc;    /* Insert image at pos 0 */
 for(k=1;k<nsc;k++) {
   rows=(rows+1)/2;  /* Half size output */
   cols=(cols+1)/2;   
   bf_list[k]=ibuffer_new(rows,cols,1);
#ifdef DEBUG
   printf("scale %d ",k);
   ibuffer_pdims(bf_list[k]);
#endif
 }
 return(bf_list);
}

/*
 *  Set an ibuffer to all ones
 *  The buffer is assumed to be
 *  of size MxNx1
 */
void set_to_ones(ibuffer *bf_image) {
int rows,cols,k,l;

 rows=bf_image->rows;
 cols=bf_image->cols;
 printf("setting %dx%d ibuffer to ones.\n",rows,cols);
 for(k=0;k<rows;k++) {
   for(l=0;l<cols;l++) {
     bf_image->data[k+l*rows]=1;
   }
 }
}

/*
 *  Filter bl_imk[0] and bl_imc[0] to generate
 *  successive scales in the pyramid.
 * 
 */
void sbinfilt_pyramid(buffer **bl_imk,ibuffer **bl_imc,int nsc,
		  fpnum dmax,fpnum cmin,int roi_side,int miter) {
int k;
 if(roi_side>0) {
   for(k=0;k<nsc-1;k++) {
#ifdef DEBUG
     printf("scale %d->%d\n",k,k+1);
#endif
     sbinfilt2d(bl_imk[k],bl_imc[k],bl_imk[k+1],bl_imc[k+1],dmax*dmax,cmin,roi_side,miter);
   }
 } else {
   for(k=0;k<nsc-1;k++) {
#ifdef DEBUG
     printf("scale(d12) %d->%d\n",k,k+1);
#endif
     sbinfilt2d12(bl_imk[k],bl_imc[k],bl_imk[k+1],bl_imc[k+1],dmax*dmax,cmin,miter);
   }
 }
}
/*
 *  Make label image
 */
void make_label_image(buffer **bl_result,buffer **bl_imk,ibuffer **bl_imc,
		      int nsc,int lowsc,fpnum dmax) {
int k,l,regions,new_seeds,rows,cols,ndim;
buffer *bf_pvec0,*bf_pvec1;
ibuffer *bf_labelim0,*bf_labelim1;
ibuffer *bf_cntl,*bf_csc0,*bf_csc1; 

  bf_labelim0=ibuffer_new(1,1,1);        /* Empty labelim */
  ndim=bl_imk[0]->ndim;
  regions=0;
  bf_pvec0=buffer_new(ndim,regions,1);  /* Empty pvec */
  bf_csc0=ibuffer_new(1,regions,1);     /* Empty csc */

  for(k=nsc-1;k>lowsc;k--) {
#ifdef DEBUG
    printf("scale %d->%d regions=%d\n",k,k-1,regions);
#endif
    rows=bl_imk[k-1]->rows;
    cols=bl_imk[k-1]->cols;
    bf_labelim1=ibuffer_new(rows,cols,1);
    new_seeds=propagate_regions(bf_labelim0,bf_labelim1,
				bl_imk[k-1],bl_imc[k-1],bf_pvec0,dmax*dmax);
    
    bf_pvec1=buffer_new(ndim,regions+new_seeds,1);

    /* Set detection scale for blobs */
    bf_csc1=ibuffer_new(1,regions+new_seeds,1);
    for(l=0;l<regions;l++) {
      bf_csc1->data[l]=bf_csc0->data[l];
    }
    for(l=regions;l<regions+new_seeds;l++) {
      bf_csc1->data[l]=k;
    }
    find_new_seeds(bf_labelim1,bl_imk[k-1],bl_imc[k-1],bf_pvec0,
		   bf_pvec1,regions,new_seeds);
    regions+=new_seeds;
    /* Free old buffers */
    buffer_free(bf_pvec0);bf_pvec0=bf_pvec1;
    ibuffer_free(bf_labelim0);bf_labelim0=bf_labelim1;
    ibuffer_free(bf_csc0);bf_csc0=bf_csc1;
  }
  for(k=lowsc;k>1;k--) {
#ifdef DEBUG
    printf("scale %d->%d no new regions\n",k,k-1);
#endif
    rows=bl_imk[k-1]->rows;
    cols=bl_imk[k-1]->cols;
    bf_labelim1=ibuffer_new(rows,cols,1);
    new_seeds=propagate_regions(bf_labelim0,bf_labelim1,
				bl_imk[k-1],bl_imc[k-1],bf_pvec0,dmax*dmax);
    /* Free old buffers */
    ibuffer_free(bf_labelim0);bf_labelim0=bf_labelim1;    
  }
  /* Last scale */
#ifdef DEBUG
  printf("scale 1->0 counting\n");
#endif
  bf_labelim1=ibuffer_new(bl_imk[0]->rows,bl_imk[0]->cols,1);
  bf_cntl=ibuffer_new(regions*(regions-1)/2,1,1);
  propagate_regions_cnt(bf_labelim0,bf_labelim1,bl_imk[0],
			bl_imc[0],bf_pvec0,bf_cntl,dmax*dmax);
    
  /* Free old buffers */
  ibuffer_free(bf_labelim0);bf_labelim0=bf_labelim1;

  /* Refine outlines of blobs */
  bf_labelim1=ibuffer_new(bl_imk[0]->rows,bl_imk[0]->cols,1);
  expand_regions_cnt(bf_labelim0,bf_labelim1,bl_imk[0],bl_imc[0],
		     bf_pvec0,bf_cntl,dmax*dmax);
  ibuffer_free(bf_labelim0);bf_labelim0=bf_labelim1;

  /* Store result in out array */
  bl_result[0]=(buffer *)bf_labelim0;
  bl_result[1]=(buffer *)bf_cntl;
  bl_result[2]=bf_pvec0;
  bl_result[3]=(buffer *)bf_csc0;
}
/*
 *
 */
void merge_and_cleanup(buffer **bl_result,buffer *bf_mvec1,buffer *bf_pvec1,
		       ibuffer *bf_csc1,ibuffer *bf_cntl,fpnum minc,
		       int amin,fpnum dmax) {
int ndim,blobs,mblobs;
buffer *bf_mvec2,*bf_pvec2,*bf_mvec3,*bf_pvec3;
ibuffer *bf_csc3,*bf_out_ind;

  ndim  = bf_pvec1->rows;
  blobs = bf_pvec1->cols;
  
  bf_mvec2=buffer_new(6,blobs,1);
  bf_pvec2=buffer_new(ndim,blobs,1);
  bf_out_ind=ibuffer_new(1,blobs,1);
  
  /*mblobs=bloblist_merge_cnt(bf_mvec1,bf_pvec1,bf_mvec2,bf_pvec2,bf_out_ind,bf_cntl,minc);*/
  mblobs=bloblist_merge_cnt2(bf_mvec1,bf_pvec1,bf_mvec2,bf_pvec2,bf_out_ind,bf_cntl,minc,dmax*dmax); 
  mblobs=bloblist_mark_invalid(bf_mvec2,bf_out_ind,amin);
  
  bf_mvec3=buffer_new(6,mblobs,1);
  bf_pvec3=buffer_new(ndim,mblobs,1);
  bf_csc3=ibuffer_new(1,mblobs,1);
  
  bloblist_compact(bf_mvec2,bf_pvec2,bf_csc1,bf_mvec3,bf_pvec3,bf_csc3,bf_out_ind);
  buffer_free(bf_mvec2);
  buffer_free(bf_pvec2);
  /* Store result in out array */
  bl_result[0]=bf_mvec3;
  bl_result[1]=bf_pvec3;
  bl_result[2]=(buffer *)bf_csc3;
  bl_result[3]=(buffer *)bf_out_ind;
}
/*
 *
 */
void extract_blobs(buffer *bf_image,ibuffer *bf_cert,buffer **bl_lout,
		   fpnum dmax,fpnum cmin,int roi_side,int miter,
		   int lowsc,fpnum minc,int amin) {
  
  buffer **bl_imk;
  ibuffer **bl_imc;
  ibuffer *bf_labelim;
  buffer *bf_pvec0,*bf_mvec1,*bf_pvec1;
  buffer *bf_mvec2,*bf_pvec2;
  ibuffer *bf_csc0,*bf_csc2,*bf_cntl,*bf_out_ind;
  int ndim,blobs,mblobs,nsc,k;
  clock_t t0,t1,t2;

#ifdef DEBUG
  printf("Image ");
  buffer_pdims(bf_image);
#endif
  t0=clock();  /* TIC */

  bl_imk=imk_pyramid_new(bf_image);
  bl_imc=imc_pyramid_new(bf_cert);

  nsc=number_of_scales(bf_image);
  
  printf("sbinfilt_pyramid: nsc=%d  dmax=%g cmin=%g roi_side=%d miter=%d\n",
	 nsc,dmax,cmin,roi_side,miter);

  t1=clock();  /* TIC */

  sbinfilt_pyramid(bl_imk,bl_imc,nsc,dmax,cmin,roi_side,miter);

  t2=clock();  /* TOC */
  write_time_diff("sbinfilt_pyramid took",t1,t2);

  printf("make_label_image: lowsc=%d dmax=%g\n",lowsc,dmax);

  t1=clock();  /* TIC */

  make_label_image(bl_lout,bl_imk,bl_imc,nsc,lowsc,dmax);
  bf_labelim = (ibuffer *)bl_lout[0];
  bf_cntl    = (ibuffer *)bl_lout[1];
  bf_pvec0   = bl_lout[2]; 
  bf_csc0    = (ibuffer *)bl_lout[3]; 

  t2=clock();  /* TOC */
  write_time_diff("make_label_image took",t1,t2);

  ndim =bf_pvec0->rows;
  blobs=bf_pvec0->cols;
  bf_pvec1=buffer_new(ndim,blobs,1);
  bf_mvec1=buffer_new(6,blobs,1);
  printf("blobs=%d\n",blobs);
 
  printf("compute_moments\n");

  t1=clock();  /* TIC */

  compute_moments(bf_labelim,bf_image,bf_pvec1,bf_mvec1);

  t2=clock();  /* TOC */
  write_time_diff("compute_moments took",t1,t2);

  printf("merge_and_cleanup: minc=%g amin=%d\n",minc,amin);

  t1=clock();  /* TIC */

  merge_and_cleanup(bl_lout,bf_mvec1,bf_pvec1,
		    bf_csc0,bf_cntl,minc,amin,dmax);
  bf_mvec2   = bl_lout[0];
  bf_pvec2   = bl_lout[1];
  bf_csc2    = (ibuffer *)bl_lout[2];
  bf_out_ind = (ibuffer *)bl_lout[3];
  
  labelim_compact(bf_labelim,bf_out_ind);

  t2=clock();  /* TOC */
  write_time_diff("merge_and_cleanup took",t1,t2);
  write_time_diff("Total time: ",t0,t2);

  mblobs=bf_pvec2->cols;
  printf("blobs=%d\n",mblobs);
 
  /* Free all buffers except at scale 0 */

  for(k=1;k<nsc;k++) {
    if(bl_imk[k]!=NULL) {
#ifdef DEBUG
      printf("Freeing imk buffer %d\n",k);
#endif
      buffer_free(bl_imk[k]);     
    }
    if(bl_imc[k]!=NULL) {
#ifdef DEBUG
      printf("Freeing imc buffer %d\n",k);
#endif
      ibuffer_free(bl_imc[k]);     
    }
 }
  free(bl_imk);
  free(bl_imc);
  ibuffer_free(bf_cntl);
  buffer_free(bf_pvec0);
  ibuffer_free(bf_csc0);
  buffer_free(bf_mvec1);
  buffer_free(bf_pvec1);
  ibuffer_free(bf_out_ind);
  bl_lout[3]=(buffer *)bf_labelim;  /* Return label image */
}