vpoint.c

image processing including fourier,wavelet,segmentation etc.

#define TMDLMeaning(T,ie,x,y) (T[ie].mdlmeaning[T[ie].nx*(y)+(x)])
#define TVP(T,ie,x,y)         (T[ie].vp[T[ie].nx*(y)+(x)])
#define TMDLVP(T,ie,x,y)      (T[ie].mdlvp[T[ie].nx*(y)+(x)])
#define TIsValid(T,ie,x,y)    (T[ie].isvalid[T[ie].nx*(y)+(x)])
#define TIsValid2(T,ie,x,y)   ( (0<=(x)) && ((x)<T[ie].nx) && \
			        (0<=(y)) && ((y)<T[ie].ny) && \
			        T[ie].isvalid[T[ie].nx*(y)+(x)] )
#define TSegs(T,ie,x,y)       (T[ie].segs[T[ie].nx*(y)+(x)])
#define TMDLSegs(T,ie,x,y)    (T[ie].mdlsegs[T[ie].nx*(y)+(x)])


/* Identification of pairs of opposite tiles in the second
   semi-circle of the last ring */
#define TOpposite(T,x,y) (((x) >= T[EXTERIOR_TILING].nx/2) ? ((x)-T[EXTERIOR_TILING].nx/2) : (x)+T[EXTERIOR_TILING].nx/2)

/* Alternative to TMeaning that deals with out of bounds indices.
   Useful to compare a tile to its neighbours
 */
double TMeaning2(T,ie,ix,iy)
     Tiling* T;
     int ie,ix,iy;
{
  int nx, ny, x, y;
  nx = T[ie].nx;
  ny = T[ie].ny;

  if (ie==INTERIOR_TILING)
    if ((ix<0) || (iy<0) || (ix>=nx) || (iy>=ny))
      return -Inf;
    else
      return TMeaning(T,ie,ix,iy);
  else {    /* extend x (theta) circularly */
    x = ix;
    while (x<0) x += nx;
    while (x>=nx) x-= nx;
    y = iy;
    if ((y<0) || (y>=2*ny-IdentifyOppositeTiles))
      return -Inf;
    else if (y>=ny)
      return TMeaning(T,ie,TOpposite(T,x,2*ny-1-y-IdentifyOppositeTiles),
		           2*ny-1-y-IdentifyOppositeTiles);
    else
      return TMeaning(T,ie,x,y);
  }

}



/* Create and initialize a new Tiling of the plane into vanishing regions
   for a given angular precision level (ntheta orientations) */
Tiling* newTiling(ntheta,p,p_inf,M)
   int ntheta;   /* Number of orientations in the exterior tiling */
   double *p;    /* (output) probability of a line meeting any of the tiles */
   double *p_inf;/* (output) probability of a line meeting infinite tiles
		             (p_inf <= p) */
   int *M;       /* (input) total number of segments, and */
                 /* (output) total number of valid tiles */
{
  double dtheta, pint, *theta, *q, dxy, *xx, *yy;
  int t,nq,nx,ny;
  int dx[4],dy[4],isvalid,max_size;
  int ie,ix,iy,i;
  Tiling* T;

  /* Used to iterate on 4-neighbours */
  /* dx={0,1,0,1}; */ dx[0]=0;dx[1]=1;dx[2]=0;dx[3]=1;
  /* dy={0,0,1,1}; */ dy[0]=0;dy[1]=0;dy[2]=1;dy[3]=1;

  /* Create the tiling for this precision level */
  T =  (Tiling*) calloc(2,sizeof(Tiling));

  /* Angular precision for all regions at this precision level */
  dtheta  = M_PI/(double)ntheta;
  /* Probability of internal tiles */
  pint = 4.0*sin(dtheta)/M_PI;

  /* Compute boundaries of external tiles in
     normalized polar coordinates (theta,q) */
  theta = malloc((ntheta+1)*sizeof(double));
  for(t=0;t<=ntheta;t++) theta[t] = 2.0*dtheta*(double)t;
  q = qtile(pint,dtheta,&nq,p_inf);


  /* Compute boundaries of internal tiles in
     normalized pixel coordinates x,y. */
  dxy = 2*sin(dtheta);
  nx = ny = (int) ceil(2.0/dxy);
  xx = malloc((nx+1)*sizeof(double));
  for (t=0;t<=nx;t++)
    xx[t] = -1.0+dxy*t;
  yy = malloc((ny+1)*sizeof(double));
  for (t=0;t<=ny;t++)
    yy[t] = -1.0+dxy*t;


  /* Create data structure for exterior tiles */
  T[EXTERIOR_TILING].nx = ntheta;
  T[EXTERIOR_TILING].ny = nq;
  T[EXTERIOR_TILING].x = theta;
  T[EXTERIOR_TILING].y = q;
  T[EXTERIOR_TILING].meaning = (double*) calloc(ntheta*nq,sizeof(double));
  T[EXTERIOR_TILING].vp      = (int*)    calloc(ntheta*nq,sizeof(int));
  T[EXTERIOR_TILING].isvalid = (int*)    calloc(ntheta*nq,sizeof(int));
  T[EXTERIOR_TILING].segs    = (SegList*)  calloc(ntheta*nq,sizeof(SegList));
  T[EXTERIOR_TILING].mdlmeaning = (double*) calloc(ntheta*nq,sizeof(double));
  T[EXTERIOR_TILING].mdlvp      = (int*)    calloc(ntheta*nq,sizeof(int));
  T[EXTERIOR_TILING].mdlsegs    = (SegList*)  calloc(ntheta*nq,sizeof(SegList));

  /*  Create data structure for interior tiles */
  T[INTERIOR_TILING].nx = nx;
  T[INTERIOR_TILING].ny = ny;
  T[INTERIOR_TILING].x = xx;
  T[INTERIOR_TILING].y = yy;
  T[INTERIOR_TILING].meaning = (double*) calloc(nx*ny,sizeof(double));
  T[INTERIOR_TILING].vp      = (int*)    calloc(nx*ny,sizeof(int));
  T[INTERIOR_TILING].isvalid = (int*)    calloc(nx*ny,sizeof(int));
  T[INTERIOR_TILING].segs    = (SegList*)  calloc(nx*ny,sizeof(SegList));
  T[INTERIOR_TILING].mdlmeaning = (double*) calloc(nx*ny,sizeof(double));
  T[INTERIOR_TILING].mdlvp      = (int*)    calloc(nx*ny,sizeof(int));
  T[INTERIOR_TILING].mdlsegs    = (SegList*)  calloc(nx*ny,sizeof(SegList));


  /* Initial size for segs Flists = expected number of segments at each tile */
  max_size = (int)rint(max(((double)*M)/sqrt((double)(nx*ny+ntheta*nq)),1.0));

  /* Initialize all tiles */
  *M=0;
  for (ie=0;ie<=1;ie++)
    for (iy=0;iy<T[ie].ny;iy++)
      for (ix=0;ix<T[ie].nx;ix++) {
	TMeaning(T,ie,ix,iy) = TMDLMeaning(T,ie,ix,iy) = -Inf;
	TVP(T,ie,ix,iy)      = TMDLVP(T,ie,ix,iy)      = 0;
	if (ie==INTERIOR_TILING) {
	  /* test if one of the corners is inside unit circle
	     (i.e. inside the image domain) */
	  isvalid=0;
	  for(i=0;i<4;i++)
	    if (hypot(T[ie].x[ix+dx[i]],T[ie].y[iy+dy[i]])<=1.0) {
	      isvalid = 1;
	      break;
	    }
	  TIsValid(T,ie,ix,iy) = isvalid;
	} else { /* ie == EXTERIOR_TILING */
	  if ((IdentifyOppositeTiles==1) && (iy==(T[ie].ny-1)))
	    /* Tile is invalid if it belongs to the second semi-circle
	       of the last (infinite) ring */
	    isvalid = (T[ie].x[ix]<M_PI)?1:0 ;
	  else
	    isvalid = 1;
	  TIsValid(T,ie,ix,iy) = isvalid;
	}
	if (isvalid) {
	  (*M)++;
	  TSegs(T,ie,ix,iy)    = newSegList(max_size);
	  TMDLSegs(T,ie,ix,iy) = newSegList(max_size);
	} else if ((ie == EXTERIOR_TILING) && (IdentifyOppositeTiles == 1)) {
	  /* In the last (infinite) ring of tiles we identify tiles
	     in the second semi-circle with their opposite tile
	     in the first semi-circle */
	  TSegs(T,ie,ix,iy)    = newSegListRef(TSegs(T,ie,TOpposite(T,ix,iy),iy));
	  TMDLSegs(T,ie,ix,iy) = newSegListRef(TMDLSegs(T,ie,TOpposite(T,ix,iy),iy));
	} else
	  TSegs(T,ie,ix,iy) = TMDLSegs(T,ie,ix,iy) = NULL;
      }

  *p = pint;
  return T;
}

void deleteTiling(T)
     Tiling* T;
{
  int ie, ix, iy;
  SegList S;

  for (ie=0;ie<=1;ie++) {
    /**/
    for (iy=0;iy<T[ie].ny;iy++)
      for (ix=0;ix<T[ie].nx;ix++) {
	S = (TSegs(T,ie,ix,iy));    if (S!=NULL) {deleteSegList(S); TSegs(T,ie,ix,iy)=NULL;};
	S = (TMDLSegs(T,ie,ix,iy)); if (S!=NULL) {deleteSegList(S); TMDLSegs(T,ie,ix,iy)=NULL;};
      }
    /**/
    free(T[ie].x);
    free(T[ie].y);
    free(T[ie].isvalid);
    free(T[ie].meaning);
    free(T[ie].vp);
    free(T[ie].segs);
    free(T[ie].mdlmeaning);
    free(T[ie].mdlvp);
    free(T[ie].mdlsegs);
  }
  free(T);

}

/* -------------------------------------------------------------------------
   Data structure to trace back the tiles met by a given segment.
   Segments is an Flists.
   Each element s of Segments is a 4-Flist which contains:
   s->data = { isvalid, x0, y0, x1, y1, nfa [, precision] }
             (float array with 6-7 elements identifying the segment)
   s->values = a list of tiles met by the segment s.
               Each 4-tuple (i,ie,ix,iy) in the list represents
	       one of these tiles and can be used e.g. as
	       TMeaning(Tilings[i],ie,iy,ix) to find its meaningfulness.
   ------------------------------------------------------------------------- */
#ifdef SegmentsAsFimage
Flists newSegments(allsegs)
     Fimage allsegs;
{
  int N, dim, i, j;
  Flists Segments;
  Flist s;

  N = allsegs->nrow;
  dim = allsegs->ncol;
  if ((dim!=6) && (dim!=5))
    mwerror(ERROR,1,"Segments Fimage should be 5xN or 6xN");

  /* convert list of segments to Flists */
  Segments = mw_change_flists(NULL,N,N);
  for(i=0;i<N;i++) {
    s = Segments->list[i] = mw_change_flist(NULL,0,0,4);
    s->data = calloc(dim+1,sizeof(float));
    s->data_size = (dim+1)*sizeof(float);
    ((float*) (s->data))[0] = 1.0; /* isvalid field */
    for(j=0;j<dim;j++)
      ((float*) (s->data))[j+1]=allsegs->gray[i*dim+j];
  }

  return Segments;
}

void deleteSegments(Segments)
     Flists Segments;
{
  int i;
  Flist s;

  for(i=0;i<Segments->size;i++) {
    s = Segments->list[i];
    free(s->data);
    s->data = NULL;
    s->data_size = 0;
  }

  mw_delete_flists(Segments);
  
}
#else
#endif

/* Convenience macros to access the attributes of a segment */
#define SegIsValid 0
#define SegX0 1
#define SegY0 2
#define SegX1 3
#define SegY1 4
#define SegMeaning 5
#define SegPrecision 6
#define STiles(S,j)     ((S)->list[j]) /* 4-Flist of tiles
					  met by j-th segment */
#define SSeg(S,j)       ((float*) ((S)->list[j]->data))
#define SIsValid(S,j)   (SSeg(S,j)[SegIsValid])
#define SX0(S,j)        (SSeg(S,j)[SegX0])
#define SY0(S,j)        (SSeg(S,j)[SegY0])
#define SX1(S,j)        (SSeg(S,j)[SegX1])
#define SY1(S,j)        (SSeg(S,j)[SegY1])
#define SMeaning(S,j)   (SSeg(S,j)[SegMeaning])
#define SPrecision(S,j) (((S)->list[j]->data_size >= 6*sizeof(float))? \
                          SSeg(S,j)[SegPrecision] : \
                          NAN)


/* -------------------------------------------------------------------------
   Auxiliary functions for TAddSegment
   ------------------------------------------------------------------------- */

/* Cross-product in R3 or join or meet in P2 */
void cross_prod(a,b,c)
     double *a,*b,*c;
{
  c[2] =  a[0]*b[1]-a[1]*b[0];
  c[0] =  a[1]*b[2]-a[2]*b[1];
  c[1] =  a[2]*b[0]-a[0]*b[2];
}


/* Normalized polar coordinates of the line supporting a segment */
void polar_coords(seg,R,Xcenter,Ycenter,rho,phi)
     float* seg;
     double R,Xcenter,Ycenter;
     double *rho,*phi;
{
  double x0,y0, x1,y1, dx,dy, px,py;
  double X0[3],X1[3],OX[3],P[3],C[3],L1[3],L2[3];

  x0 = seg[SegX0]+0.5; y0 = seg[SegY0]+0.5; /* first end-point */
  x1 = seg[SegX1]+0.5; y1 = seg[SegY1]+0.5; /* second end-point */
  dx = x1-x0;          dy = y1-y0;

  /*
     Projective coordinates of:
     X0 = first endpoint of segment
     X1 = second endpoint of segment
     OX = direction orthogonal to segment
     C  = origin of coordinate system
  */
  X0[0] = x0;          X0[1] = y0;          X0[2] = 1.0;
  X1[0] = x1;          X1[1] = y1;          X1[2] = 1.0;
  OX[0] = -1.0*dy;     OX[1] = dx;          OX[2] = 0.0;
  C[0]  = Xcenter;     C[1]  = Ycenter;     C[2]  = 1.0;

  /* Find the foot P of the perpendicular L2
     from C to supporting line of segment
  */
  cross_prod(X0,X1,L1); /* line L1 through X0,X1 */
  cross_prod(C, OX,L2); /* line L2 from C orthogonal to L1 */
  cross_prod(L1,L2,P ); /* intersection P of L1 and L2 */

  /* Now compute polar coordinates of line:
     rho = (non-signed) distance from C to line, divided by R
     phi = oriented angle in [0,2pi] from e1 to P-C
  */
  px = P[0]/P[2]-Xcenter;  py = P[1]/P[2]-Ycenter;
  *rho = hypot(px,py)/R;
  *phi = atan2(py,px);

  if ((*phi)<0.0)
    *phi = *phi + 2.0*M_PI;

}

/* Normalized projective coordinates of the line supporting a segment */
void proj_coords(seg,R,Xcenter,Ycenter,L)
     float* seg;
     double R,Xcenter,Ycenter;
     double *L; /* should be called with an allocated L[3] */
{
  double x0,y0, x1,y1;
  double X0[3],X1[3];

  x0 = seg[SegX0]+0.5; y0 = seg[SegY0]+0.5; /* first end-point */
  x1 = seg[SegX1]+0.5; y1 = seg[SegY1]+0.5; /* second end-point */

  /*
     Projective coordinates (in normalized coord system):
     X0 = first endpoint of segment
     X1 = second endpoint of segment
  */
  X0[0] = (x0-Xcenter)/R;  X0[1] = (y0-Ycenter)/R;  X0[2] = 1.0;
  X1[0] = (x1-Xcenter)/R;  X1[1] = (y1-Ycenter)/R;  X1[2] = 1.0;

  /* line L through X0,X1 */
  cross_prod(X0,X1,L); 

}

/*
  Append a record to the end of the Flist.
  (This is used to add a segment j to an accumulator TSegs(T,ie,ix,iy))
 */
void FlistAdd(l,tuple)
     Flist l;
     float *tuple;
{
  int i,sz,dim;
  sz = l->size;
  dim = l->dim;

  if ((l->size == l->max_size) && (!mw_enlarge_flist(l)))
    mwerror(FATAL,1,"Not enough memory to enlarge flist");
  for(i=0;i<dim;i++)
    l->values[sz*dim+i] = tuple[i];

  (l->size)++;
}

Flist FlistCat(l1,l2)
     Flist l1,l2;
{
  int i, j; 

  if (l1->dim != l2->dim)
    mwerror(INTERNAL,1,"[FlistCat] Flists should have the same dimension!");

  for (i=0; i<l2->size; i++)
    if ((l1->size == l1->max_size) && (!mw_enlarge_flist(l1)))
      mwerror(ERROR,1,"Not enough memory to enlarge Flist!");
    else {
      for (j=0; j<l1->dim; j++)
	l1->values[l1->size*l1->dim + j] = l2->values[i*l2->dim + j];
      l1->size++;
    }

  return l1;
}

Flist FlistFlush(l1)
     Flist l1;
{
  l1->size=0;
  return l1;
}