disocclusion.c

image processing including fourier,wavelet,segmentation etc.

/*--------------------------- MegaWave2 Command -----------------------------*/
/* mwcommand
name = {disocclusion};
version = {"2.0"};
author = {"Simon Masnou"};
function = {"Disocclusion using global minimisation of cost by dynamic recursive programing"};
usage = {
  'e':[energy_type=1]->energy_type    
       "Energy of a level line : 0 = only length, 1 = only angle (default), otherwise = angle+length",
  'a'->angle        "If used then the orientation of each entering level line is computed more accurately on a ball of radius 4",
  input->Input      "Input occluded Cimage",
  holes->Holes      "Input Fimage containing the only occlusions",
  output<-Output    "Output disoccluded Cimage"
  };
*/
/*----------------------------------------------------------------------
 v2.0: revised version, added -e option, no more -m option (S.Masnou)
----------------------------------------------------------------------*/


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



/***************************************************************************/
/** This program performs disocclusion of simply connected sets. Any set with
    a hole will be automatically filled before the disocclusion process

    Crucial convention : x denotes VERTICAL axis and y HORIZONTAL axis !! **/
/***************************************************************************/

/**********************************************************************/
/* IMPORTANT NOTE : at the end of the labeling process, labels of each 
   connected component are sorted increasingly according to their
   order of appearance */
/**********************************************************************/

/* We assume than there are no more than LABELS_NUMBER labels obtained at first pass
 (before updating) */
#define LABELS_NUMBER 100000

int u_compar_i(u,v)  /*  Called by function qsort for sorting decreasingly */
int *u,*v;
  {
    if ((*u)<(*v)) return (1);
    if ((*u)==(*v)) return (0);
    return (-1);
  }

/****************************************************************************/
/*
The algorithm used here is the following : let us consider the case where we want
to compute the connected components of the set of points with value g (the principle
is analogous for the computation of the connected components of the complement).
We create a subsidiary image, one pixel larger in each direction than the input 
image, where all pixels of input frame with value g have value 1, while the others 
are set to 0. Then we scan the interior of this subsidiary image (borders are not
taken into account) from top to bottom and left to right (only one scan is enough). 
Each pixel with value 1 whose 4 upper and left pixels (in case of 8-connectivity) 
are 0 is given a new label (as if it was a new component).
If some of these 4 pixels are nonzero then the current pixel is given the greatest label
among these 4 and we update the transcoding table which tells us which label is 
connected to which one, i.e. which components are part of the same connected components.
Once the image has been totally scanned, we simply have to replace the value (label)
of each nonzero pixel by the smallest label to which it is associated.
*/
/****************************************************************************/

void mise_a_jour_transcode(transcode,a,b)

int transcode[LABELS_NUMBER];
int a,b;

  {
    if (!b) return;
    if (transcode[a]==0)
      transcode[a]=b;
    else
      {
	if (b!=transcode[a])
	  if (b>transcode[a])
	    mise_a_jour_transcode(transcode,b,transcode[a]);
	  else
	    {
	      mise_a_jour_transcode(transcode,transcode[a],b);
	      transcode[a]=b;
	    }
      }
  }
    
/****************************************************************************/

/* In the transcoding table, a label l1 may be associated with another label l2,
itself associated with a third label l3, etc... The following function simply
performs a recursive association of each label l1,l2,l3 with the smallest possible
label in the chain. */

void refresh(transcode,refresh_transcode,i,first)

int transcode[LABELS_NUMBER],refresh_transcode[LABELS_NUMBER];
int i;
int *first;  
    
  {
    if (transcode[i])
      {
	refresh(transcode,refresh_transcode,transcode[i],first);
	refresh_transcode[i]=*first;
      }
    else
      *first=i;
  }

/****************************************************************************/

void fconnected(In,line,col,FOREGROUND,NUMBER,not_writing_on,complement,connectivity)
float *In;
int line,col;
float FOREGROUND;
int *NUMBER;
char *not_writing_on,*complement,*connectivity;

  {
    int Line=line+2,Col=col+2;
    int *Output;
    register float *ptrin;
    register int *ptrout;
    register int dx,dy,i;
    int kernel[4];
    int label,first;
    int transcode[LABELS_NUMBER]; /* transcodage table obtained at first pass and constructed 
				     following the rule : if a point is connected with two different labels 
				     (more than two is impossible) one writes in the table the link
				     label_max->label_min if label_max is not yet linked, else, say we
				     already have label_max->label' one makes by recursivity 
				     label'->label_min or label_min->label' depending on 
				     max(label',label_min)
				     */
    int refresh_transcode[LABELS_NUMBER];  /* same as transcode except that each label is linked with
					      the smallest value in the chain to which it belongs */
    int normalize_transcode[LABELS_NUMBER];
    register int *ptr_trans,*ptr_ref,*ptr_norm;
    int norme;

    /* The frame used for processing is 2 pixels wider and higher than the input to avoid
       problems with borders */

    Output=(int*)malloc(Line*Col*sizeof(int));
    if (Output==NULL)
      mwerror(FATAL,1,"Not enough memory !\n");
    
    for (i=0,ptr_trans=transcode,ptr_ref=refresh_transcode,
	 ptr_norm=normalize_transcode;i<LABELS_NUMBER;i++,ptr_trans++,ptr_ref++,ptr_norm++)
      {*ptr_trans=0;*ptr_ref=0;*ptr_norm=0;}
    
    for (dy=0,ptrout=Output;dy<Col;dy++,ptrout++) *ptrout=0; 
    for (dx=2,ptrin=In;dx<Line;dx++)
      {
	*ptrout=0;ptrout++;
	for (dy=2;dy<Col;dy++,ptrout++,ptrin++)
	  if (complement)
	    if (*ptrin!=FOREGROUND) *ptrout=1;
	    else *ptrout=0;
	  else
	    if (*ptrin==FOREGROUND) *ptrout=1;
	    else *ptrout=0;
	*ptrout=0;ptrout++;
      }
    for (dy=0,ptrout=Output;dy<Col;dy++,ptrout++) *ptrout=0; 
   
    label=1;
    for (dx=2,ptrout=Output+(Col+1);dx<Line;dx++,ptrout+=2)
      for (dy=2;dy<Col;dy++,ptrout++)
	if (*ptrout)
	  {
	    if (label>LABELS_NUMBER)
	      mwerror(FATAL,1,"There are more than %d labels !\n",LABELS_NUMBER);
	    /* 'kernel' contains the four upper nearest neighbours of current point
	               * * *
		       * + x
		       x x x
	       if 4-connectivity is used then only 2 neighbors are taken into account */
	    kernel[0]=*(ptrout-1);
	    if (connectivity) kernel[1]=(-1);
	    else kernel[1]=*(ptrout-(Col+1));
	    kernel[2]=*(ptrout-Col);
	    if (connectivity) kernel[3]=(-1);
	    else kernel[3]=*(ptrout-(Col-1));
	    qsort(kernel,4,sizeof(int),u_compar_i);
	    /* kernel values are sorted decreasingly */
	    if (kernel[0]==0)
	      *ptrout=label++;   /* This point is given a new label */
	    else
	      {
		*ptrout=kernel[0];
		if (kernel[1])
		  if (kernel[1]!=kernel[0])
		    if (kernel[1]!=transcode[kernel[0]])
		      mise_a_jour_transcode(transcode,kernel[0],kernel[1]);
	      }
	  }
    first=0;
    for (i=label-1;i>0;i--)
      if ((transcode[i]) && (!(refresh_transcode[i])))
	refresh(transcode,refresh_transcode,i,&first);
    
    
    /* We want each connected component to be labelled between 1 and the number
       of components. This is done in the following step */
    norme=1;
    for (i=1,ptr_ref=refresh_transcode+1;i<label;i++,ptr_ref++){
      if (!(*ptr_ref))
	*ptr_ref=i;
      if (normalize_transcode[*ptr_ref]==0)
	{normalize_transcode[*ptr_ref]=norme++;}}
    *NUMBER=norme-1;

     /* Writing on image */
    if (!not_writing_on)
      {
	for (dx=2,ptrout=Output+(Col+1),ptrin=In;dx<Line;dx++,ptrout+=2)
	  for (dy=2;dy<Col;dy++,ptrout++,ptrin++)
	    {
	      if (*ptrout)
		*ptrin=(float)(normalize_transcode[refresh_transcode[*ptrout]]);
	      else
		*ptrin=0;
	    }
      }
    free(Output);
  }  




/****************************************************************************/

#define Square(u) (u)*(u)
#define Norm(u,v) sqrt(Square(u)+Square(v))
#define Min(u,v) (((u)<(v))?(u):(v))
#define Max(u,v) (((u)<(v))?(v):(u))
#define Sgn(u)   (((u)<(-1e-8))?(-1):1)
#define CPL 0.2

typedef struct a_jordan {  /* structure used to describe the T-junctions */
  double x,y; /* coordinates of a point of the Jordan curve */
  double vx,vy;  /* vector coordinates corresponding to the direction of the level line */
  char direction; /* direction of the associated line = either 0,2,4,6 (see below) */
  unsigned char gray1,gray2; /* values separated by current node (sorted counterclockwise) */
  struct a_jordan *junction;  /* jordan element it is linked to (NULL if no link) */
  struct a_jordan *previous,*next;
} jordan;

/* Numbering convention for direction

                      1 0 7  
                       \|/   
                     2-- --6 
                       /|\	
                      3 4 5
*/

	
/*  Called by function qsort for sorting unsigned char values by increasing values */
static int compar_uc(u,v)  
unsigned char *u,*v;
{
  if ((*u)>(*v)) return (1);
  if ((*u)==(*v)) return (0);
  return (-1);
}

int line_number,col_number;  /* image line and column numbers  */
int IONumber; /* Number of T-junctions asspciated with different level lines */

/* These three arrays are used for the computation of the optimal set of T-junctions
   (see ComputeOptimalSet) */
double **energy; 
int **correspond;
jordan **element;

char *anglep;
int energy_criterion;
int globx,globy;

static double CCost();
static void ComputeOptimalSet();

unsigned char *IImage,*OImage; /* Pointers to Input and Output images */
float *LImage; /* Pointer to Label image (image containing the only occlusions) */
float lvalue;  /* Current occlusion label */

typedef struct occlusionPoint{ /* structure used to describe the occlusion */
  int pos; /* coordinates of point */
  } occlusionPoint;
occlusionPoint *occlusion;
int numcolored; /* number of points already colorized */


typedef struct a_frontier { /* structure used to describe the vertices of the occlusion boundary */
  double x,y; /* coordinates of a point of the Jordan curve */
  struct a_frontier *previous,*next;
} frontier;

frontier *jb; /* First element of frontier */
int numfrontier=0;
char *pivots; /* A geodesic path can be made of different segments; 
		 at each segment apex (pivot) one needs the local conformation 
		 of the level line */
int *fpivots;

typedef struct Pxy_t {
    double x, y;
} Pxy_t;

typedef Pxy_t Ppoint_t;

typedef struct Ppoly_t {
    Ppoint_t *ps;
    int pn;
} Ppoly_t;

typedef Ppoly_t Ppolyline_t;

#define ISCCW 1
#define ISCW  2
#define ISON  3

#define DQ_FRONT 1
#define DQ_BACK  2

#ifndef TRUE
#define TRUE 1
#define FALSE 0
#endif

#define prerror(msg) mwerror(FATAL,1,"\nlibpath/%s:%d: %s\n", __FILE__, __LINE__, (msg))

#define POINTSIZE sizeof (Ppoint_t)

typedef struct pointnlink_t {
    Ppoint_t *pp;
    struct pointnlink_t *link;
} pointnlink_t;

#define POINTNLINKSIZE sizeof (pointnlink_t)
#define POINTNLINKPSIZE sizeof (pointnlink_t *)

typedef struct tedge_t {
    pointnlink_t *pnl0p;
    pointnlink_t *pnl1p;
    struct triangle_t *ltp;
    struct triangle_t *rtp;
} tedge_t;

typedef struct triangle_t {
    int mark;
    tedge_t e[3];
} triangle_t;

#define TRIANGLESIZE sizeof (triangle_t)

typedef struct deque_t {
    pointnlink_t **pnlps;
    int pnlpn, fpnlpi, lpnlpi, apex;
} deque_t;

pointnlink_t *pnls, **pnlps;
int pnll;

triangle_t *tris;
int trin, tril;

deque_t dq;

Ppoint_t *ops;

/****************************************************************************/

/* ccw test: counterclockwise, clockwise or collinear */
int ccw (p1p,p2p,p3p)
Ppoint_t *p1p,*p2p,*p3p;
  {
    double d;

    d = ((p1p->y - p2p->y) * (p3p->x - p2p->x)) -
            ((p3p->y - p2p->y) * (p1p->x - p2p->x));
    return (d > 0) ? ISCCW : ((d < 0) ? ISCW : ISON);
  }

/****************************************************************************/