vpoint.c

image processing including fourier,wavelet,segmentation etc.

/*
  Append tile (i,ie,ix,iy) to the end of L, only if:
  a) this tile is valid and
  b) it is not yet in L, and
  c) no tile equivalent to (i,ie,ix,iy) has been added to L.
  These three tests are done via l=TSegs(T,ie,ix,iy), with T = Tilings[i],
  as follows (if HasSTL is undefined)
  a)   l is non-null
  b,c) there is no l2 = TSegs(Tilings[i2],ie2,ix2,iy2) such that l2==l and
       the tuple (i2,ie2,ix2,iy2) is in L.
  (This is used to check that a segment is not counted twice.
   The search is not too expensive because we build one of these lists
   per segment, precision level. The size of ls is at most 2*nx+4*ny).
   Observe that this implementation ignores the last argument j.
*/
void STilesAddUnique(L,Tilings,i,ie,ix,iy,j)
     Flist L;   /* 4-Flist containing the tiles met by the j-th segment */
     Tiling **Tilings;
     int i,ie,ix,iy; /* 4-tuple representing the tile to be added */
     int j;          /* current segment */
{
  int k,ismember;
  Flist l2;
  SegList l;
  int i2,ie2,ix2,iy2;
  Tiling *T;

  T = Tilings[i];

  /* We test equality via l = TSegs(T,ie,ix,iy)
     since in case IdentifyOppositeTiles==1
     two different 4-tuples (i,ie,ix,iy) might refer to the same tile
  */
  l = TSegs(T,ie,ix,iy);
  if (l!=NULL) {
#ifdef HasSTL
    ismember = SegListInsert(l,j);
#else
    for (ismember=0, k=0; (!ismember) && (k<L->size); k++) {
      i2  = (int) (L->values[4*k+0]);
      ie2 = (int) (L->values[4*k+1]);
      ix2 = (int) (L->values[4*k+2]);
      iy2 = (int) (L->values[4*k+3]);
      l2 = TSegs(Tilings[i2],ie2,ix2,iy2);
      if (l2==l) {
	ismember=1;
	break;
      }
    }
#endif
    if (!ismember)
      if ((L->size == L->max_size) && (!mw_enlarge_flist(L)))
	mwerror(FATAL,1,"Not enough memory to enlarge Flist");
      else {
	L->values[4*L->size + 0] = (float) i;
	L->values[4*L->size + 1] = (float) ie;
	L->values[4*L->size + 2] = (float) ix;
	L->values[4*L->size + 3] = (float) iy;
	L->size++;
      }
  }
}

/* -------------------------------------------------------------------------
   TAddSegment - Updates Tiling with all intersections of a segment
                 with all tiles
   If the j-th segment in segs doesn't match the current precision level
   do nothing and return 0.
   Otherwise return 1 after doing the following... 
   For all tiles t in T such that the jth segment meets t, update T,
   i.e.:
   - add j to the list of segments meeting t (avoiding duplicates!)
   
   We should also (but this can be done at the end):
   - update meaning to -log10(nfa(n,N))
     (where n is the updated length of the segs list meeting t and
      N is the total number of segments)
   - update vp to meaning>eps
  -------------------------------------------------------------------------
*/
int TAddSegment(Tilings,i,S,j,R,Xcenter,Ycenter)
     Tiling **Tilings;
     int i; /* index to precision level */
     Flists S;
     int j; /* index to current segment */
     double R,Xcenter,Ycenter;
{
  double rho,phi,cos_phi_theta;
  int    nd,ntheta,nx,ny,ie;
  double *d,*theta,dtheta,theta0,*x,*y,dx,dy,xmin,ymin,xmax,ymax;
  int    id,itheta,id1,itheta1,itheta2,ix,iy;
  double            d1, theta1, theta2,xx,yy;
  int k;
  double l[3];
  Flist  L,csegs;
  float tuple;
  int max_tiles;
  Tiling *T;
  double length;

  T = Tilings[i];

  /* compute angular precision of this segment ... */
  length = hypot(SX0(S,j)-SX1(S,j),SY0(S,j)-SY1(S,j));
  ie     = EXTERIOR_TILING;
  dtheta = T[ie].x[1]-T[ie].x[0];
  /* ... if it is too coarse do not add segment */
  if (atan2(1.0,length)>dtheta)
    return 0;

  ntheta = T[EXTERIOR_TILING].nx;
  nd     = T[EXTERIOR_TILING].ny;
  nx     = T[INTERIOR_TILING].nx;
  ny     = T[INTERIOR_TILING].ny;
  max_tiles = 4*nd+2*ntheta + 2*max(nx,ny);

  /*
    we do not start with L = STiles(S,j),
    because it only contains tiles from previous precision levels <i,
    whereas in this call we add tiles at precision level i.
    Thus starting with L = STiles(S,j), would unnecessarily slow down
    STilesAddUnique.
    Instead, we append the list L of tiles at precision i to STiles(S,j)
    at the end of this routine.
  */
  L = mw_change_flist(NULL,max_tiles,0,4);
  
  /*** A) Compute all intersections of j-th segment with exterior tiles ***/

  /* Normalized polar coordinates of the supporting line of j-th segment*/

     polar_coords(SSeg(S,j)                 ,R,Xcenter,Ycenter,&rho,&phi);

  ie     = EXTERIOR_TILING;
  nd     = T[ie].ny;
  d      = T[ie].y;
  ntheta = T[ie].nx;
  theta  = T[ie].x;
  dtheta = T[ie].x[1]-T[ie].x[0];
  theta0 = T[ie].x[0];

  /* Find intersection points of line (rho,phi) with each circle */
  for (id=0;id<nd;id++) {
    theta1 = phi-acos(rho/d[id]);
    theta1 = theta1+((theta1<      0.0)? 2*M_PI : 0.0);
    theta2 = acos(rho/d[id])+phi;
    theta2 = theta2-((theta2>=2.0*M_PI)? 2*M_PI : 0.0);
    itheta1 = (int) floor((theta1-theta0)/dtheta);
    itheta2 = (int) floor((theta2-theta0)/dtheta);
    /*  Add tiles touching both intersection points
	(theta1,d[id]) and (theta2,d[id])
	to the list L, avoiding duplicates */
    STilesAddUnique(L,Tilings,i,ie,itheta1,id  ,j);
    STilesAddUnique(L,Tilings,i,ie,itheta2,id  ,j);
    if (id>0) {
      STilesAddUnique(L,Tilings,i,ie,itheta1,id-1,j);
      STilesAddUnique(L,Tilings,i,ie,itheta2,id-1,j);
    }
  };
  /* Find intersection point of line (rho,phi) with each radius */
  for (itheta=0;itheta<ntheta;itheta++) {
    cos_phi_theta = cos(phi-theta[itheta]); /* when cos==0
					       we force infinite tile */
    d1 = (fabs(cos_phi_theta)<EPS)?Inf:rho/cos(phi-theta[itheta]);
    if (d1>=d[0]) {
      id1=0; while((id1<nd) && (d[id1]<d1)) id1++;
      if (id1>0) id1--;
      /*  Add tiles touching the intersection point (theta[itheta],d1)
	  to the list L, avoiding duplicates */
      STilesAddUnique(L,  Tilings,i,ie,itheta,  id1,j);
      if (itheta>0)
	STilesAddUnique(L,Tilings,i,ie,itheta-1,id1,j);
      else /* itheta==0, its circular neighbour is ntheta-1 */
	STilesAddUnique(L,Tilings,i,ie,ntheta-1,id1,j);
    }
  }

  /*** B) Compute all intersections of j-th segment with interior tiles ***/

     proj_coords(SSeg(S,j)                 ,R,Xcenter,Ycenter,l);

  ie = INTERIOR_TILING;
  nx   = T[ie].nx;
  x    = T[ie].x;
  dx   = T[ie].x[1]-T[ie].x[0];
  xmin = T[ie].x[0];
  xmax = T[ie].x[nx];
  ny   = T[ie].ny;
  y    = T[ie].y;
  dy   = T[ie].y[1]-T[ie].y[0];
  ymin = T[ie].y[0];
  ymax = T[ie].y[ny];

  if ((fabs(phi)         <=M_PI_4) ||
      (fabs(phi-M_PI)    <=M_PI_4) ||
      (fabs(phi-2.0*M_PI)<=M_PI_4)) {
    /* Rather "vertical line" :
       find intersection with horizontal lines of the tiling */
    for (iy=0; iy<ny; iy++) {
      yy = y[iy];
      xx = -(l[1]*yy+l[2])/l[0];
      ix = (int) floor((xx-xmin)/dx);
      if (TIsValid2(T,ie,ix,iy  ))
	STilesAddUnique(L, Tilings,i,ie,ix,iy  ,j);
      if (TIsValid2(T,ie,ix,iy-1))
	STilesAddUnique(L, Tilings,i,ie,ix,iy-1,j);
    }
  } else {
    /* Rather "horizontal line" :
       find intersection with vertical lines of the tiling */
    for (ix=0; ix<ny; ix++) {
      xx = x[ix];
      yy = -(l[0]*xx+l[2])/l[1];
      iy = (int) floor((yy-ymin)/dy);
      if (TIsValid2(T,ie,ix  ,iy))
	STilesAddUnique(L, Tilings,i,ie,ix  ,iy,j);
      if (TIsValid2(T,ie,ix-1,iy))
	STilesAddUnique(L, Tilings,i,ie,ix-1,iy,j);
    }
  }
#ifndef HasSTL
  /* Add segment j to all tiles in L */
  tuple = (float) j;
  for (k=0;k<L->size;k++) {
    i  = (int) (L->values[4*k+0]);
    ie = (int) (L->values[4*k+1]);
    ix = (int) (L->values[4*k+2]);
    iy = (int) (L->values[4*k+3]);
    csegs = TSegs(Tilings[i],ie,ix,iy);
    FlistAdd(csegs,&tuple);
  }
#endif

  STiles(S,j) = FlistCat(STiles(S,j),L);

  mw_delete_flist(L);

  return 1;

}

/*
  Given a low-precision tiling TT, and a high-precision tiling T.
  Find all tiles tt(ie2,ix2,iy2) in TT meeting tile t(ie,ix,iy) in T.
  Return the number of intersecting tiles.

  With the current implementation we still miss:
    - intersections between interior and exterior tiles at the same or different levels
   
 */
int FindIntersection(T,ie,ix,iy,TT,ie2,ix2,iy2)
     Tiling *T,*TT;
     int  ie,  ix,   iy;
     int *ie2,*ix2,*iy2;
{
  int k,n,m,new_tile;
  double *x,*y,*y2,*x2;
  double dx2, dy2, xmin2, ymin2;
  int ny2;
  int IE2, IX2, IY2;

  /* Used to iterate on 4-neighbours */
  int DX[4],DY[4];
  DX[0]=0;DX[1]=1;DX[2]=0;DX[3]=1;
  DY[0]=0;DY[1]=0;DY[2]=1;DY[3]=1;

  /* Compute values defining finer grid */
  *ie2 = ie;
  x    = T[ie].x;                   x2    = TT[*ie2].x;
/*dx   = T[ie].x[1]-T[ie].x[0];*/   dx2   = TT[*ie2].x[1]-T[*ie2].x[0];
/*xmin = T[ie].x[0];           */   xmin2 = TT[*ie2].x[0];

/*ny   = T[ie].ny*/                 ny2   = TT[*ie2].ny;
  y    = T[ie].y;                   y2    = TT[*ie2].y;
/*dy   = T[ie].y[1]-T[ie].y[0];*/   dy2   = TT[*ie2].y[1]-TT[*ie2].y[0];
/*ymin = T[ie].y[0];*/              ymin2 = TT[*ie2].y[0];


  for (k=0,n=0; k<4; k++) {
    IE2 = ie2[n] = ie;
    IX2 = ix2[n] = (int) floor((x[ix+DX[k]]-xmin2)/dx2);
    if (ie==INTERIOR_TILING)
      IY2 = iy2[n] = (int) floor((y[iy+DY[k]]-ymin2)/dy2);
    else {
      /* y = "distance to center" is not regularly spaced */
      /* if (y[iy]>=y2[0]) not necessary as in TAddSegment */
      IY2 = iy2[n] = (k==0)?0:iy2[0];
      while( (iy2[n] < ny2) && (y2[iy2[n]]<y[iy+DY[k]]) ) (iy2[n])++;
      if (iy2[n]>0) (iy2[n])--;
      IY2 = iy2[n];
    }
    new_tile = (1==1);
    for (m=0; m<n; m++)
      new_tile = new_tile && ( (ie2[n]!=ie2[m]) || (iy2[n]!=iy2[m]) || (ix2[n]!=ix2[m]) );
    if (new_tile)
      n++;

    if ((fabs(x2[IX2]-x[ix+DX[k]])<dx2/1000.0) && (IX2>0)) {
      ie2[n] = IE2;
      ix2[n] = IX2-1;
      iy2[n] = IY2;
      new_tile = (1==1);
      for (m=0; m<n; m++)
	new_tile = new_tile && ( (ie2[n]!=ie2[m]) || (iy2[n]!=iy2[m]) || (ix2[n]!=ix2[m]) );
      if (new_tile)
	n++;
    }
      
    if ((fabs(y2[IY2]-y[iy+DY[k]])<y2[IY2]/1000.0) && (IY2>0)) {
      ie2[n] = IE2;
      ix2[n] = IX2;
      iy2[n] = IY2-1;
      new_tile = (1==1);
      for (m=0; m<n; m++)
	new_tile = new_tile && ( (ie2[n]!=ie2[m]) || (iy2[n]!=iy2[m]) || (ix2[n]!=ix2[m]) );
      if (new_tile)
	n++;
    }
      
    
  }

  return n;
}


/*------------------------------------------------------------*/
/*                         MAIN MODULE                        */
/*------------------------------------------------------------*/


int vpoint(imagein,allsegs,output,segs,eps,all,masked,verbose,maskedVPs)
Fimage imagein;
Fimage allsegs;
Flist output;
Flists segs;
double *eps;
char   *all;
char   *masked;
char   *verbose;
int    *maskedVPs;
{
  int N,*M,min_pl,max_pl,n_pl,i,j,k,ntheta,itn;
  int    ie, ix, iy;
  int i2,ie2,ix2,iy2;
  int max_i,max_ie,max_ix,max_iy;
  float max_meaning;
  int    IE[9],IX[9],IY[9]; /* used to hold the output of FindIntersection */
  int    n;
  int NVP[2],nvp,nvp1,nvp2,nvp3;
  double R,X0,Y0;
  double **B,**Binf;
  double *p,meaning,*p_inf,threshold;
  Tiling **Tilings, *T, *TT;
  Flists Segments;
  Flist tiles;
  SegListIterator it;
  SegList CSegs;
  float  *csegs;
  int    nsegs,besttile;
  float x1,x2,x3,x4;
  float y1,y2,y3,y4;
  float d1,d2,d3,d4;
  float t1,t2,t3,t4;
  int searching_masked;
  float tuple[4];
  char *msg1, *msg2;
  Flist VPlist[2];


  /* Build data structure to hold pointers from segments to tiles
     (useful for MDL) */
     Segments = newSegments(allsegs);

  /* Number of detected lines */
     N = Segments->size;

  /* Radius and center of circumscribed circle containing image domain */
     R = hypot((double)imagein->nrow,(double)imagein->ncol)/2.0; 
     X0 = (double)imagein->ncol/2.0;
     Y0 = (double)imagein->nrow/2.0;
  /* Establish series of n_pl dyadic angular precision levels
     for the vanishing region tilings */
     min_pl = 4;
     max_pl = 9;
     n_pl = max_pl-min_pl+1;
     M = calloc(n_pl,sizeof(int));
     B = calloc(n_pl,sizeof(double*));
     Binf = calloc(n_pl,sizeof(double*));
     p = calloc(n_pl,sizeof(double));
     p_inf = calloc(n_pl,sizeof(double));
 
  /* Detection threshold for NFA is epsilon
     divided by the number of precision levels.
     Detection threshold for meaningfulness (-log10(NFA)) is then
     -log10(epsilon/n_pl) == -log10(epsilon) + log10(n_pl)
  */
     threshold = *eps + log10(n_pl);

     Tilings = malloc(n_pl*sizeof(Tiling**));


  /*
     Create the Tiling data structure for all precision levels
   */