vpoint.c

image processing including fourier,wavelet,segmentation etc.

/*--------------------------- Commande MegaWave -----------------------------*/
/* mwcommand
  name = {vpoint};
  version = {"1.0"};
  author = {"Andres Almansa"};
  function = {"Detect (MDL) maximal meaningful vanishing regions"};
  usage = {  
 'a'->all           "to get all meaningful v.points (before MDL)",
 'm'->masked        "to get also masked v.points",
 'v'->verbose       "display verbose messages during computations",
 'e':[eps=0.]->eps  "-log10(max. number of false alarms), default 0",
 's':csegs<-segs    "store segments indexes associated to each vanishing region (output Flists)",
  imagein->imagein  "input image (only used to extract image domain)",
  allsegs->allsegs  "input segments (5xn)-Fimage (x1 y1 x2 y2 -log10(nfa))",
  output<-output    "detected vanishing regions: 9-Flist (x1 y1 x2 y2 x3 y3 x4 y4 -log10(nfa))",
  NVP<-vpoint       "number of detected vanishing regions",
  NVP_masked<-maskedVPs  "number of detected masked vanishing regions"
          };
*/
/*
V 1.01 : (JF) backslash removed (lexical error on SunOS 5.x) 
*/

#include <stdio.h>
#include <math.h>
#include "mw.h"

/*
   IdentifyOppositeTiles -- Global setting
                            defines how infinite opposite tiles are treated.

   IdentifyOppositeTiles = 1
     With this setting opposite infinite tiles are considered as a single tile.
     This is more symmetric, but we have to use an upper bound for the
     probability of these tiles, not the real probability.
   IdentifyOppositeTiles = 0
     With this setting opposite infinite tiles are considered as different.
     Here we use the real probability for each one, and opposite infinite
     tiles are considered as neighbours when computing local maxima of
     meaningfulness.
     The disadvantage of this setting is that exactly parallel lines
     fall exactly between both opposite infinite tiles, and noise will
     decide which one wins.
*/ 
static int IdentifyOppositeTiles=0;

/*
   SegmentsAsFimage -- Global setting
   if defined, input segments are read from an Nx5 or Nx6 Fimage,
   as in the old version of align.
   otherwise, input segments are read from a 5-Flist or a 6-Flist
   with N records, as produced by the new version of align
 */
#define SegmentsAsFimage

/* ****************************************************************
   Ref-counted Set implementation using MegaWave Flist
   **************************************************************** */
typedef  Flist SegList;
SegList  newSegList(n)
     int n;
{
  SegList l = mw_change_flist(NULL,n,0,1);
  l->data = malloc(sizeof(int)); /* used as ref_count */
  l->data_size = 1;
  (*(int*)l->data) = 1;
  return l;
}
int  deleteSegList(L)
     SegList L;
{
  if (--(*(int*)L->data) == 0) {
    free(L->data);
    L->data = NULL;
    L->data_size = 0;
    mw_delete_flist(L);
  }
  return 0;
}
SegList  newSegListRef(L)
     SegList L;
{
  (*(int*)L->data)++;
  return L;
}
#define  SegListInsert(L,i)   mwerror(FATAL,1,"SegListInsert not implemented in non-HasSTL mode.\n Use FlistAdd instead.\n")
#define  SegListSize(L)       (L->size)
typedef  float* SegListIterator;
#define  SegListBegin(L)      (L->values)
#define  SegListEnd(L,it)     (it>=(L->values+(L->size)))
#define  SegListNext(it)      (it++)
#define  SegListValue(it)     ((int)*it)
Flist    SegListCopy(L)
     SegList L;
{
  SegListIterator p;
  Flist l;

  l = mw_change_flist(NULL,SegListSize(L),0,1);
  
  for(p=SegListBegin(L); !SegListEnd(L,p); SegListNext(p))
    l->values[l->size++] = (float) SegListValue(p);

  l->data = malloc(sizeof(int)); /* used as ref_count */
  l->data_size = 1;
  (*(int*)l->data) = 1;

  return l;
}

/* ****************************************************************
   End of: Ref-counted Set implementation using MegaWave Flist
   **************************************************************** */



#define EPS 1.0e-16 /* Machine double precision */
#define Inf (-log(0.0))

#define max(x,y) (x>y)?x:y
#define _(I,x,y) ((I)->gray[(y)*((I)->ncol)+(x)]) /* 0 <= x < I->ncol
						     0 <= y < I->nrow */



/*---------------------------------------------------------------------------
   Compute a partition of the plane in regions of equal probability p,
   that a random line visible in the image meets the region.
  ---------------------------------------------------------------------------
   For the iterior of the image domain we choose square regions of side
      N*2*sin(theta)
   where
      theta is the minimal angular precision of the line segments, and
      N is the size of the image.
   Thus the probability of an image line meeting such regions is:
      pint = 8*N*sin(theta)
   For the exterior of the image we choose portions of sectors of a circle
   centered at the image center. These sectors have an angle 2*theta, and
   the portions of sectors are limited between two cords at distance
   d1 and d2 from the image center. Instead of di it will be more convenient
   to work with
      betai = acos(cos(theta)*R/di),
   where R=sqrt(2)*N is the radius of the circle circumscribed to the image.
   For such an exterior region the probability of an image line meeting the
   region is
      pext(beta1,beta2,theta)
   as computed by the corresponding function below. We also provide a
   function to compute the derivative of pext with respect to beta2:
      pext_prime(beta1,theta)
   When d1 remains fixed and d2 tends to infinity (beta2 tends to pi/2)
   the probability of an image line meeting the region becomes
      pinf(beta1,theta)
   as computed by the corresponding function below.
   fzero_convex finds a zero within the given interval of the function
      f(x) = pext(beta1,x,theta)-pint
   which is convex, with known derivative, so we can use a modified Newton
   method, which ensures a given precision both on x and f(x).
   Finally qtile uses fzero_convex to find the good values of beta which are
   needed to obtain a tiling of the plane into regions with equal
   probability. It proceeds as follows: first it fixes
      d_1 = R*cos(theta) or equivalently beta_1=0
   then for all i>1 it uses fzero_convex to find beta_i or d_i such that
      pext(beta_i-1,beta_i,theta) = pint
   The iteration stops when
      pinf(beta_i,theta) <= pint
   so the last (infinite) region has a probability between 0 and pint.
  ---------------------------------------------------------------------------*/

double pext(beta1,beta2,theta)
     double beta1,beta2,theta;
{
  double f1,f2;
  if (cos(beta1)<EPS*1.0e5)
    f1 = M_PI_2;
  else
    f1 = beta1+1.0/cos(beta1)-tan(beta1);

  if (cos(beta2)<EPS*1.0e5)
    f2 = M_PI_2;
  else
    f2 = beta2+1.0/cos(beta2)-tan(beta2);

  return ( 2.0*theta + f2-f1 ) / (double) M_PI;
}

double pext_prime(beta2,theta)
     double beta2,theta;
{
  double fp2;
  if (cos(beta2)<EPS*1.0e5)
    fp2 = -0.5;
  else
    fp2 = (sin(beta2)-1.0)/cos(beta2)/cos(beta2);

  return ( 1.0 + fp2 ) / (double) M_PI;
}

double pinf(beta1,theta)
     double beta1,theta;
{
  return pext(beta1,M_PI_2,theta);
}

double fzero_convex(a0,b0,tolx,tolf,beta1,dtheta,p)
     double a0,b0,tolx,tolf,beta1,dtheta,p;
     /* Finds the zero x0 of a convex function myfun(x) on an interval [a0,b0].
        We assume that:
	- its derivative myfunprime(x) is also known and
        - since myfun is convex,
          myfunprime is nonnegative and increasing on [a,b]
	The root x0 is given with precision tolx on the x axis
        and tolf on the y axis, i.e.  abs(myfun(x0))<tolf.

     */
{
  double a,b,fa,fb,fpa,fpb;
  a = a0; b=b0;
  fa = pext(beta1,a,dtheta)-p; /* == myfun(a); */
  fb = pext(beta1,b,dtheta)-p; /* == myfun(b); */

  mwdebug("Iterating fzero_convex...\n");
  while (fabs(fb-fa)>tolf && fabs(b-a)>tolx) {
    fpb = pext_prime(b,dtheta); /* == myfunprime(b); */
    a = a-fa/fpb; fa =  pext(beta1,a,dtheta)-p; /* == myfun(a); */
    b = b-fb/fpb; fb =  pext(beta1,b,dtheta)-p; /* == myfun(b); */
    mwdebug(".");
 }

  return (a+b)/2.0;

}

double* qtile(p,dtheta,nq,p_inf)
     double p;
     double dtheta;
     int* nq;
     double* p_inf;
{
  double tol,b;
  double *beta;
  int i;
  tol = 1.0e-8;
  
  *nq  = 1;
  beta = malloc(sizeof(double));
  beta[0] = b = 0.0;

  while ((*p_inf=pinf(b,dtheta))>p) {
    b = fzero_convex(b,M_PI_2,tol,tol,b,dtheta,p);
    beta = realloc(beta,((*nq)+1)*sizeof(double));
    beta[(*nq)++] = b;
  }

  for (i=1; i<*nq; i++) {
    mwdebug("beta=%8.3f, q=%8.3f (p-pint)/pint=%f\n",beta[i],1.0/cos(beta[i]),
	    (pext(beta[i-1],beta[i],dtheta)-p)/p);
  }
    mwdebug("beta=%8.3f, q=%8.3f (p-pint)/pint=%f\n",M_PI_2,HUGE_VAL,
	    (pinf(beta[*nq-1],dtheta)-p)/p);
    mwdebug("\n");

  /* convert beta to normalized distances q */
  for (i=0; i<*nq; i++) {
    beta[i] = 1.0/cos(beta[i]);
  }

  return beta;
}

void polar2cart(theta,q,R,Xcenter,Ycenter,x,y)
     double theta,q;           /* normalized polar coordinates */
     double R,Xcenter,Ycenter; /* parameters used for normalization */
     float  *x,*y;             /* non-normalized cartesian coordinates */
{
  *x = R*q*cos(theta)+Xcenter;
  *y = R*q*sin(theta)+Ycenter;
}


/*------------------------------------------------------------*/
/*       compute P(k,l) : array out[] of size n+1 * n+1       */
/*   P[k][l] =     sum(i=k..l) binom(l,i) p^i (1-p)^(l-i)     */
/*   P[k][l] = P[ (n+1)*l+k ]                                 */
/*------------------------------------------------------------*/


double *tab(n,p)
int n;
double p;
{
  double *out;
  int adr1,adr2,x,y;
  double lambda,q;

  q = 1.0-p;
  out = (double *)calloc((n+1)*(n+1),sizeof(double));

  /*** compute proba (=x among y) ***/
  out[0] = 1.0;
  for (y=1,adr2=0;y<=n;y++) {
    adr1 = adr2;
    adr2 += n+1;    
    out[adr2] = q*out[adr1];
    for (x=1;x<=y;x++) 
      out[adr2+x] = p*out[adr1+x-1] + q*out[adr1+x];
  }  

  /*** sum to obtain proba (>=k among y) ***/
  for (y=1,adr1=n+1;y<=n;y++,adr1+=n+1) 
    for (x=y-1;x>=0;x--) 
      out[adr1+x] += out[adr1+x+1];

  /*** multiply by m (number of segments) to obtain expectation***
  for (adr1=(n+1)*(n+1);--adr1>=0;)
    out[adr1] *= m;
  */

  return out;
}

double *tab2(n,p)
int n;
double p;
/* alternative implementation: it computes directly the binomial tail,
   without computing the binomial density first */
{
  double *out;
  int adr1,adr2,x,y;
  double lambda,q;

  q = 1.0-p;
  out = (double *)calloc((n+1)*(n+1),sizeof(double));

  /*** compute proba (>=x among y) ***/
  out[0] = 1.0;                                      /* P(0,0) = 1 */
  for (y=1,adr2=0;y<=n;y++) {
    adr1 = adr2; /* column y-1 */
    adr2 += n+1; /* column y   */ 
    out[adr2]   = 1.0;                               /* P(0,y) = 1.0 */
    out[adr2+y] = p*out[adr1+(y-1)];                 /* P(y,y) =
							   p * P(y-1,y-1)
						      */
    for (x=1;x<y;x++) 
      out[adr2+x] = p*out[adr1+x-1] + q*out[adr1+x]; /* P(x,y) =
                                                           p * P(x-1,y-1)
							 + q * P(x,  y-1)
						      */
  }  

  return out;
}

double* binomial_tail(n,p)
int n;
double p;
{
  double *in,*out;
  int x;

  in = tab2(n,p);
  out = calloc(n+1,sizeof(double));

  for (x=0;x<=n;x++)
    out[x] = in[(n+1)*n+x];

  free(in);

  return out;
}

/* -------------------------------------------------------------------------
   Structure to represent a tiling of the image plane.
   Actually a tiling will be
      Tiling T[2]
   where T[INTERIOR_TILING] is a tiling of the image domain
   and   T[EXTERIOR_TILING] is a tiling of its complement
   ------------------------------------------------------------------------- */
#define INTERIOR_TILING 0
#define EXTERIOR_TILING 1
typedef struct tiling {
  int nx,ny;       /* number of tiles = nx*ny */
  double *x;       /* nx+1 array with x-coords of the borders of the tiling */
  double *y;       /* ny+1 array with y-coords the borders of the tiling */
                   /* For exterior tilings (x,y)=(theta,q)
		      in normalized polar coordinates */
  double *meaning; /* nx*ny array with the meaningfulness (-log10(nfa))
                      for each region */
  int *vp;         /* nx*ny array equals 1 if the region is meaningful
                      with locally maximal meaningfulness */
  int *isvalid;    /* nx*ny array equals 1 if the corresponding region
                      is valid. (Used to eliminate interior tiles which are
                      outside of the circular "image domain", 
		      as well as infinite exterior tiles
		      which are split into two.
		      It is also used to remove tiles
		      that have already been found to be MDL-maximal-meaningful). */
  SegList *segs;   /* nx*ny array of Flists. The Flist corresponding to
		      region (i,j) contains pointers to the input segs Flist
		      for all segments whose supporting line meets region (i,j)
		    */
  double *mdlmeaning;
  int *mdlvp;
  SegList *mdlsegs; /* same as meaning, vp and segs, but after MDL */
} Tiling;

/* Convenience macros to access a "pixel" of the tiling */
#define TMeaning(T,ie,x,y)    (T[ie].meaning[T[ie].nx*(y)+(x)])