disocclusion.c

image processing including fourier,wavelet,segmentation etc.

    jbi=(frontier*)NULL;
  }

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

/* Copies the frontier into a polygonal line whose structure is compatible 
   with Mitchell's shortest path algorithm. The frontier structure was built
   in such a way that it contains only those points of bifurcation of 
   the polygonal line enclosing the occlusion */

void copy_frontier(x)
int x;
  {
    frontier *jbi;
    int n;

    if (!(polypoints = (Ppoint_t *) malloc (POINTSIZE * numfrontier)))
      prerror ("cannot malloc polypoints");
    for (n=0,jbi=jb;n<numfrontier;n++,jbi=jbi->next)
      {polypoints[n].x=jbi->x+x;polypoints[n].y=jbi->y;}
  }

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

/* Colorization of pixels on both sides of the current geodesic path */

void geodfill(p1,p2,gray1,gray2,change)
int p1,p2;
unsigned char gray1,gray2;
char change;
  {
    /* LImage == (-1) means that the point has already been modified */

    if (LImage[p1]==lvalue) {OImage[p1]=gray1;LImage[p1]=(-1);occlusion[numcolored++].pos=p1;}
    if (LImage[p2]==lvalue) {OImage[p2]=gray2;LImage[p2]=(-1);occlusion[numcolored++].pos=p2;}
    if (!change)
      {if (LImage[p2]==(-1)) OImage[p2]=gray2;}
    else
      if (LImage[p1]==(-1)) OImage[p1]=gray1;
  }
	
/***************************************************************************/

/* Colorization of peculiar pixels in the neighborhood of geodesic path nodes  */

void pivotfill(p,gray,force)
int p;
unsigned char gray;
char force;

  {
    if (LImage[p]==lvalue) 
      {OImage[p]=gray;LImage[p]=(-1);occlusion[numcolored++].pos=p;return;}
    if (force)
      if (LImage[p]==(-1)) OImage[p]=gray;
  }

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

/* Colorization of peculiar pixels in the neighborhood of geodesic path nodes */

void fpivotfill(piv,p1,p2,gray1,gray2)
int piv,p1,p2;
unsigned char gray1,gray2;

  {
    if (p1<4) fpivots[piv+p1]=(int)gray1;
    if (p2<4) fpivots[piv+p2]=(int)gray2;
  }

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

/* Draws a two pixels wide straight line between two points. The algorithm is 
   very similar to Bresenham algorithm except that we force the line to be 
   two pixels wide and that we need some more work at each vertex of 
   the polygonal line. We use the following convention :
       xf=Horizontal coordinate of first point
       yf=Vertical coordinate of first point
       xe=Horizontal coordinate of end point
       ye=Vertical coordinate of end point */

void drawgeod(xf,yf,xe,ye,gray1,gray2,piv)
int xf,yf,xe,ye;
unsigned char gray1,gray2;
int piv; /* Current pivot number */
  {
    int dx=xe-xf,dy=ye-yf,vect=0,d2x,d2y;
    int i,j,epsilon,oldi,oldj;
    char change=0; /* 1 in case orientation is changed */
    unsigned char g;

    if (dy<0){
      change=1;
      g=gray1;gray1=gray2;gray2=g;
      dx=-dx;
      epsilon=xe;xe=xf;xf=epsilon;
      dy=-dy;
      epsilon=ye;ye=yf;yf=epsilon;}
    d2x=2*dx;d2y=2*dy;
    
    /* vect=dx*(i-yf)-dy*(j-xf) is the vectorial product (careful to the orientation!!)*/

    if (abs(dx)>=dy) /* dx is the horizontal gap, dy the vertical gap */
      {
	/*epsilon=dy-abs(dx);*/
	epsilon=-abs(dx);
	if (dx>0)
	  {
	    for (j=xf,i=yf;j<=xe;j++){
	      if (j==xe)
		if (i!=oldi) pivots[8*(piv+1-change)]=1+change;
		else pivots[8*(piv+1-change)+7]=1+change;
	      if (2*(j/2)==j)
		if (2*(i/2)!=i)
		  {
		    if (j==xf+1) fpivotfill(4*(piv+change),1,2,gray1,gray2);
		    if (j==xe-1) fpivotfill(4*(piv+1-change),0,3,gray1,gray2);
		    geodfill((i-1)/2*col_number+j/2,(i+1)/2*col_number+j/2,gray1,gray2,change);
		  }
		else
		  if (vect<=0)
		    {
		      if (j==xf+1) fpivotfill(4*(piv+change),2,4,gray1,gray2);
		      if (j==xe-1) fpivotfill(4*(piv+1-change),0,3,gray1,gray2);
		      geodfill(i/2*col_number+j/2,(i+2)/2*col_number+j/2,gray1,gray2,change);
		    }
		  else 
		    {
		      if (j==xf+1) fpivotfill(4*(piv+change),1,2,gray1,gray2);
		      if (j==xe-1) fpivotfill(4*(piv+1-change),4,0,gray1,gray2);
		      geodfill((i-2)/2*col_number+j/2,i/2*col_number+j/2,gray1,gray2,change);
		    }
	      /*else
		if (j==xe && i==ye)
		geodfill((i-1)/2*col_number+(j+1)/2,(i+1)/2*col_number+(j-1)/2,gray1,gray2,change);*/
	      oldi=i;
	      epsilon+=d2y;vect-=dy;
	      if (epsilon>=0) {i++;epsilon-=d2x;vect+=dx;}
	      if (j==xf)
		if (i!=oldi) pivots[8*(piv+change)+4]=2-change;
		else pivots[8*(piv+change)+3]=2-change;}
	  }
	else
	  {
	    for (j=xf,i=yf;j>=xe;j--){
	      if (j==xe)
		if (i!=oldi) pivots[8*(piv+1-change)+2]=1+change;
		else pivots[8*(piv+1-change)+3]=1+change;
	      if (2*(j/2)==j)
		if (2*(i/2)!=i)
		  {
		    if (j==xf-1) fpivotfill(4*(piv+change),3,0,gray1,gray2);
		    if (j==xe+1) fpivotfill(4*(piv+1-change),2,1,gray1,gray2);
		    geodfill((i+1)/2*col_number+j/2,(i-1)/2*col_number+j/2,gray1,gray2,change);
		  }
		else
		  if (vect>=0) 
		    {
		      if (j==xf-1) fpivotfill(4*(piv+change),4,3,gray1,gray2);
		      if (j==xe+1) fpivotfill(4*(piv+1-change),2,1,gray1,gray2);
		      geodfill((i+2)/2*col_number+j/2,i/2*col_number+j/2,gray1,gray2,change);
		    }
		  else
		    {
		      if (j==xf-1) fpivotfill(4*(piv+change),3,0,gray1,gray2);
		      if (j==xe+1) fpivotfill(4*(piv+1-change),1,4,gray1,gray2);
		      geodfill(i/2*col_number+j/2,(i-2)/2*col_number+j/2,gray1,gray2,change);
		    }
	      /*else
		if (j==xe && i==ye)*/ /*if (2*(i/2)!=i && (j!=xf))*/
	      /*geodfill((i+1)/2*col_number+(j+1)/2,(i-1)/2*col_number+(j-1)/2,gray1,gray2,change);*/
	      oldi=i;
	      epsilon+=d2y;vect+=dy;
	      if (epsilon>=0) {i++;epsilon+=d2x;vect+=dx;}
	      if (j==xf)
		if (i!=oldi) pivots[8*(piv+change)+6]=2-change;
		else pivots[8*(piv+change)+7]=2-change;}
	  }
      }
    else 
      {
	/*epsilon=abs(dx)-dy;*/
	epsilon=-dy;
	if (dx>0)
	  {
	    for (j=xf,i=yf;i<=ye;i++){
	      if (i==ye)
		if (j!=oldj) pivots[8*(piv+1-change)]=1+change;
		else pivots[8*(piv+1-change)+1]=1+change;
	      if (2*(i/2)==i)
		if (2*(j/2)!=j)
		  {
		    if (i==yf+1) fpivotfill(4*(piv+change),2,3,gray1,gray2);
		    if (i==ye-1) fpivotfill(4*(piv+1-change),1,0,gray1,gray2);
		    geodfill(i/2*col_number+(j+1)/2,i/2*col_number+(j-1)/2,gray1,gray2,change);
		  }
		else
		  if (vect<=0) 
		    {
		      if (i==yf+1) fpivotfill(4*(piv+change),2,3,gray1,gray2);
		      if (i==ye-1) fpivotfill(4*(piv+1-change),0,4,gray1,gray2);
		      geodfill(i/2*col_number+j/2,i/2*col_number+(j-2)/2,gray1,gray2,change);
		    }
		  else 
		    {
		      if (i==yf+1) fpivotfill(4*(piv+change),4,2,gray1,gray2);
		      if (i==ye-1) fpivotfill(4*(piv+1-change),1,0,gray1,gray2);
		      geodfill(i/2*col_number+(j+2)/2,i/2*col_number+j/2,gray1,gray2,change);
		    }
	      /*else
		if (j==xe && i==ye)*/ /*if (2*(i/2)!=i && (i!=yf))*/
	      /*geodfill((i-1)/2*col_number+(j+1)/2,(i+1)/2*col_number+(j-1)/2,gray1,gray2,change);*/
	      oldj=j;
	      vect+=dx;epsilon+=d2x;
	      if (epsilon>=0) {j++;epsilon-=d2y;vect-=dy;}
	      if (i==yf)
		if (j!=oldj) pivots[8*(piv+change)+4]=2-change;
		else pivots[8*(piv+change)+5]=2-change;}
	  }
	else
	  {
	    for (j=xf,i=yf;i<=ye;i++){
	      if (i==ye)
		if (j!=oldj) pivots[8*(piv+1-change)+2]=1+change;
		else pivots[8*(piv+1-change)+1]=1+change;
	      if (2*(i/2)==i)
		if (2*(j/2)!=j)
		  {
		    if (i==yf+1) fpivotfill(4*(piv+change),2,3,gray1,gray2);
		    if (i==ye-1) fpivotfill(4*(piv+1-change),1,0,gray1,gray2);
		    geodfill(i/2*col_number+(j+1)/2,i/2*col_number+(j-1)/2,gray1,gray2,change);
		  }
		else
		  if (vect>=0)
		    {
		      if (i==yf+1) fpivotfill(4*(piv+change),2,3,gray1,gray2);
		      if (i==ye-1) fpivotfill(4*(piv+1-change),4,1,gray1,gray2);
		      geodfill(i/2*col_number+(j+2)/2,i/2*col_number+j/2,gray1,gray2,change);
		    }
		  else 
		    {
		      if (i==yf+1) fpivotfill(4*(piv+change),3,4,gray1,gray2);
		      if (i==ye-1) fpivotfill(4*(piv+1-change),1,0,gray1,gray2);
		      geodfill(i/2*col_number+j/2,i/2*col_number+(j-2)/2,gray1,gray2,change);
		    }
	      /*else
		if (j==xe && i==ye)*/ /*if (2*(i/2)!=i && (i!=yf))*/
	      /*geodfill((i+1)/2*col_number+(j+1)/2,(i-1)/2*col_number+(j-1)/2,gray1,gray2,change);*/
	      oldj=j;
	      vect+=dx;epsilon-=d2x;
	      if (epsilon>=0) {j--;epsilon-=d2y;vect+=dy;}
	      if (i==yf)
		if (j!=oldj) pivots[8*(piv+change)+6]=2-change;
		else pivots[8*(piv+change)+5]=2-change;}
	  }
      }
  }


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

/* Here we draw a geodesic path between each optimal couple of T-junctions
   Each geodesic path is actually a two-pixels wide polygonal line where the 
   two gray values associated with the level line are represented.
   In order to draw this polygonal line, we simply draw straight lines between 
   the vertices using Bresenham algorithm. However, since this line is two pixels
   wide and since we want afterwards to propagate the values, we have to be careful
   at each vertex. */

/* The geodesic paths are successively drawn with respect to the following rule : 
   gray values are set on both sides of the first path, which can be any path whose 
   starting point differs from the starting point of the path immediately before. 
   Then, taking the starting points one after the other when the occlusion boundary 
   is walked counterclockwise, we draw the corresponding paths by forcing the value 
   on the "right" of the path and setting it on the left only if it has not been set yet. */

void geodesic(jc,x)
jordan *jc;
int x;

  {
    jordan *jci;
    Ppolyline_t *geodPoints;
    int n,pi,npivots,npi,npi1,npi2,npi3,npi4,npi5,pos,npipos;
    int gx1,gy1,gx2,gy2;
    char *ptrpiv;
    int *ptrfpiv;
    unsigned char gray1,gray2;
 
    /* Be careful : at the end of this function, jci->junction->junction=NULL !! */

    /* The two following pointers are needed for the computation of geodesic paths.
     These paths have actually already been computed in the dynamic programing stage but
     it would have been much too expensive to try to save them at that time (due to the
     necessity of performing successive updates). */
    if (!(geodPoints = (Ppolyline_t *) malloc ((IONumber/2) * sizeof(Ppolyline_t))))
      prerror ("cannot malloc geodPoints");
    if (!(endTjunctions = (Ppoint_t *) malloc (POINTSIZE * IONumber)))
      prerror ("cannot malloc endTjunctions");
    
    /* This is to ensure that the first examined level line is the first line 
       starting at the associated T-junction (ie the first line reached when 
       curve is described counterclockwise). */
    while ((jc->x==jc->previous->x) && (jc->y==jc->previous->y)) jc=jc->next;
    
    /* To avoid a double examination of the same level line */
     jci=jc;n=0;
    do{
      if (jci->junction){
	endTjunctions[n].x=jci->x+x;endTjunctions[n++].y=jci->y;
	endTjunctions[n].x=jci->junction->x+x;endTjunctions[n++].y=jci->junction->y;
	jci->junction->junction=NULL;}
      jci=jci->next;}
    while (jci!=jc);
    
    pivots=NULL;npivots=0;
    for (n=0;n<IONumber;n+=2)
      {
	Pshortestpath(endTjunctions[n],endTjunctions[n+1], &geodPoints[n/2]);
	gray1=jci->gray1;gray2=jci->gray2;
	if (geodPoints[n/2].pn>npivots)
	  {
	    npivots=geodPoints[n/2].pn;
	    if (!pivots)
	      {
		/* Each point of the line has 8 neighbors in the semi-grid */
		/* This array gives the local conformation of the line */
		if (!(pivots = (char *) malloc (npivots*8*sizeof(char)))) 
		  prerror ("cannot malloc pivots");
		/* Each point of the line has 4 integer neighbors */
		/* This array says where to put which gray value */
		if (!(fpivots = (int *) malloc (npivots*4*sizeof(int)))) 
		  prerror ("cannot malloc fpivots");
	      }
	    else
	      {
		if (!(pivots = (char *) realloc((void*)pivots,npivots*8*sizeof(char))))
		  prerror ("cannot realloc pivots");
		if (!(fpivots = (int *) realloc((void*)fpivots,npivots*4*sizeof(int))))
		  prerror ("cannot realloc fpivots");
	      }
	  }
	else npivots=geodPoints[n/2].pn;
	for (npi=0,ptrpiv=pivots;npi<npivots*8;npi++,ptrpiv++) *ptrpiv=0;
	for (npi=0,ptrfpiv=fpivots;npi<npivots*4;npi++,ptrfpiv++) *ptrfpiv=(-1);
	
	/* Pivots neighborhood      0  1  2 
	                            7  x  3 
		                    6  5  4  */
	/* "pivots[i]" takes value 1 if the line is entering at point i, 2 it the
	   line is outgoing at i. */

	/* fpivots neighborhood     0  1
	                            3  2  */
	switch (jci->direction)
	  {
	  case 0:pivots[1]=1;fpivots[1]=gray1;fpivots[0]=gray2;break;
	  case 2:pivots[7]=1;fpivots[0]=gray1;fpivots[3]=gray2;break;
	  case 4:pivots[5]=1;fpivots[3]=gray1;fpivots[2]=gray2;break;
	  case 6:pivots[3]=1;fpivots[2]=gray1;fpivots[1]=gray2;break;
	  default: prerror("Probleme de direction");break;
	  }
	switch (jci->junction->direction)
	  {
	  case 0:pivots[(npivots-1)*8+1]=2;fpivots[(npivots-1)*4]=gray1;
	    fpivots[(npivots-1)*4+1]=gray2;break;
	  case 2:pivots[(npivots-1)*8+7]=2;fpivots[(npivots-1)*4+3]=gray1;
	    fpivots[(npivots-1)*4]=gray2;break;
	  case 4:pivots[(npivots-1)*8+5]=2;fpivots[(npivots-1)*4+2]=gray1;
	    fpivots[(npivots-1)*4+3]=gray2;break;
	  case 6:pivots[(npivots-1)*8+3]=2;fpivots[(npivots-1)*4+1]=gray1;
	    fpivots[(npivots-1)*4+2]=gray2;break;
	  default:prerror("Probleme de direction");break;
	  }
	gx1=(int)(rint(2*geodPoints[n/2].ps[0].y));
	gy1=(int)(rint(2*geodPoints[n/2].ps[0].x));
	for (pi = 1; pi < geodPoints[n/2].pn; pi++)
	  {
	    gx2=(int)(rint(2*geodPoints[n/2].ps[pi].y));
	    gy2=(int)(rint(2*geodPoints[n/2].ps[pi].x));
	    drawgeod(gx1,gy1,gx2,gy2,gray1,gray2,pi-1);
	    gx1=gx2;gy1=gy2;
	  }

	/* Now we examine each node of the polygonal line and check if it is necessary
	   to colorize furthermore. This is required by the propagation function if we
	   want to avoid conflicts between gray values */
	for (pi = 0; pi < geodPoints[n/2].pn; pi++)
	  {
	    pos=(int)(rint(geodPoints[n/2].ps[pi].x-0.5))*col_number+(int)(rint(geodPoints[n/2].ps[pi].y-0.5));
	    ptrpiv=pivots+8*pi;
	    npi1=0;
	    while (*ptrpiv!=1) {npi1++;ptrpiv++;}