disocclusion.c

image processing including fourier,wavelet,segmentation etc.

	    ptrpiv-=npi1;
	    npi2=0;
	    while (*ptrpiv!=2) {npi2++;ptrpiv++;}
	    ptrfpiv=fpivots+pi*4;
	    for (npi=0;npi<4;npi++)
	      if (*(ptrfpiv+npi)==(-1))
		{
		  switch (npi)
		    {
		    case 0:npipos=pos;break;
		    case 1:npipos=pos+1;break;
		    case 2:npipos=pos+col_number+1;break;
		    default:npipos=pos+col_number;break;
		    }
		  npi3=(npi+1)%4;
		  while (*(ptrfpiv+npi3)==(-1)) npi3=(npi3+1)%4;
		  npi4=(npi+3)%4;
		  while (*(ptrfpiv+npi4)==(-1)) npi4=(npi4+3)%4;
		  if (npi*2==npi1 || npi*2==npi2) prerror("Direction fpivots = probleme !!");
		  if (npi3!=npi4)
		    if (*(ptrfpiv+npi3)!=*(ptrfpiv+npi4))
		      {
			npi5=(2*npi+1)%8;
			/* The previous setting allows to switch from "fpivots" neighborhood to 
			   "pivots" neighborhood */
			while ((npi5!=npi1)&&(npi5)!=npi2) npi5=(npi5+1)%8;
			if (npi5==npi1) {pivotfill(npipos,gray2,(char)1);*(ptrfpiv+npi)=gray2;}
			else {pivotfill(npipos,gray1,(char)0);*(ptrfpiv+npi)=gray1;}
		      }
		}
	  }
	do jci=jci->next;
	while (jci->junction==NULL);
      }
    if (pivots) {free((void*)pivots);pivots=NULL;}
    if (fpivots) {free((void*)fpivots);fpivots=NULL;}
    FreeTPath();
    if (polypoints) 
      {free((void*)polypoints);polypoints=NULL;}
    if (endTjunctions) 
      {free((void*)endTjunctions);endTjunctions=NULL;}
    if (geodPoints){
      for (n=0;n<IONumber/2;n++)
	{free((void*)(geodPoints[n].ps));geodPoints[n].ps=NULL;}
      free((void*)geodPoints);
      geodPoints=NULL;}
  }

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


/* Once the geodesic paths have been drawn, it remains to let the values propagate
   within occlusion */

void propagation(occlusionArea)
int occlusionArea;

  {
    int k,l,pos;
    int oldnum=numcolored;
    unsigned char gray;


    /* First step : 4-connectivity geodesic dilation */
    for (k=0;k<oldnum;k++){
      gray=OImage[occlusion[k].pos];
      pos=occlusion[k].pos-col_number;
      for (l=0;l<4;l++){
	if (LImage[pos]==lvalue)
	  {OImage[pos]=gray;LImage[pos]=(-1);occlusion[numcolored++].pos=pos;}
	switch (l){
	case 0:case 2: pos+=col_number-1;break;
	case 1 : pos+=2;break;}}}
    
    /* Second step : 8-connectivity dilation */
    for (k=0;k<oldnum;k++){
      gray=OImage[occlusion[k].pos];
      pos=occlusion[k].pos-col_number-1;
      for (l=0;l<4;l++){
	if (LImage[pos]==lvalue)
	  {OImage[pos]=gray;LImage[pos]=(-1);occlusion[numcolored++].pos=pos;}
	switch (l){
	case 0:case 2: pos+=2;break;
	case 1 : pos+=2*col_number-2;break;}}}

    for (k=oldnum;k<occlusionArea;k++){
      gray=OImage[occlusion[k].pos];
      pos=occlusion[k].pos-col_number-1;
      for (l=0;l<8;l++){
	if (LImage[pos]==lvalue)
	  {OImage[pos]=gray;LImage[pos]=(-1);occlusion[numcolored++].pos=pos;}
	switch (l){
	case 0:case 1:case 5:case 6 : pos++;break;
	case 2:case 4 : pos+=col_number-2;break;
	case 3 : pos+=2;break;}}}
  }




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

/* This procedure is used to "clean" the jordan curve. Zero-length line are removed, 
   as well as lines that have the same endpoints */

jordan *update_jordan_curve(jc,ExternVal)
jordan *jc;
unsigned char *ExternVal;


  {
    jordan *jlook1,*jlook2,*jnext,*jprev;
    double ox1,oy1,ox2,oy2;
    char go_on;

    /* First step : find if possible a line with stricly positive length and start from this line */
    jlook1=jc;
    while ((jlook1->x==jlook1->junction->x) && (jlook1->y==jlook1->junction->y)){
      jlook1=jlook1->next;
      if (jlook1==jc){ /* All the junctions have zero-length ! */
	if (jc->vx*jc->junction->vy-jc->vy*jc->junction->vx>0) *ExternVal=jc->gray1;
	else *ExternVal=jc->gray2;
	IONumber=0;
	return ((jordan*)NULL);}}
    jc=jlook1;
    /* At that point the line jc -- jc->junction has a stricly positive length */

    /* Second step : remove all lines (and the in-between lines) whose length is zero */
    jlook1=jc->next;
    while (jlook1!=jc)
      if ((jlook1->x==jlook1->junction->x) && (jlook1->y==jlook1->junction->y)){
	jlook2=jlook1->junction->next;
	jlook1->previous->next=jlook2;
	jlook2->previous=jlook1->previous;
	do{
	  jnext=jlook1->next;
	  IONumber-=1;
	  free((void*)jlook1);
	  jlook1=jnext;}
	while (jlook1!=jlook2);}
      else jlook1=jlook1->next;

    /*  Third step : Merge lines which have the same endpoints 
	(but possibly with different directions at endpoints !!) */
    
    /* jlook1 equals jc */
    go_on=1;
    do{
      if (IONumber==2) return (jc);
      if (jlook1->next==jlook1->junction){
	jlook1=jlook1->junction;
	if (jlook1==jc) break;}
      ox1=jlook1->x;oy1=jlook1->y;ox2=jlook1->junction->x;oy2=jlook1->junction->y;
      jlook2=jlook1->next;
      if ((jlook2->x==ox1) && (jlook2->y==oy1) && 
	  (jlook2->junction->x==ox2) && (jlook2->junction->y==oy2)){
	jlook1->gray2=jlook2->gray2;
	jlook2->junction->next->gray1=jlook1->gray2;
	if (jlook2->next==jlook2->junction){
	  jlook1->next=jlook1->junction;
	  jlook1->junction->previous=jlook1;}
	else {
	  jnext=jlook2->next;
	  jprev=jlook2->junction->previous;
	  jlook1->next=jnext;
	  jnext->previous=jlook1;
	  jlook2->junction->next->previous=jprev;
	  jprev->next=jlook2->junction->next;}
	free((void*)(jlook2->junction));
	free((void*)(jlook2));
	IONumber-=2;}
      else{
	jlook1=jlook1->next;
	if (jlook1==jc) go_on=0;}}
    while (go_on);
    return (jc);
  }

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

/* In the same time that we construct the chain "jordan" of level lines, we construct the 
   jordan curve made of the vertices of the polygonal lines enclosing the occlusion. */

jordan *allocation_jordan_curve(instart,jc,previous_direction,next_direction,i,j)
unsigned char *instart;
jordan *jc;
char previous_direction,next_direction;
int i,j;

  {
    unsigned char *ptrin;
    unsigned char gray1,gray2;

    /* The eight neighbors of a pixel are numbered using the following convention
                                1 0 7  
                                 \|/   
                               2-- --6 
                                 /|\	
                                3 4 5

       Each pixel has four vertices whose coordinates are half-integer. Level lines pass
       through this type of points
                            s1              s2
                     (i-0.5,j-0.5)  (i-0.5,j+0.5)
                                *----*
                                |    |
                                |    |
                                *----*
                     (i+0.5,j-0.5)  (i+0.5,j+0.5)
                       	    s3             s4
    */
    
    switch (previous_direction)
      {
      case 0:case 7:
	if (next_direction==6) {mwerror(FATAL,1,"Impossible 1\n");break;}
	if (previous_direction==0)
	  if (next_direction==7) {mwerror(FATAL,1,"Impossible 2\n");break;}
	/* jc=add_to_jordan(jc,(double)i-0.5,(double)j+0.5,6);*/ /* predecessor already made */
	if (next_direction==5) {add_to_frontier((double)i+0.4999,(double)j+0.5001);break;}
	ptrin=instart+(i*col_number+j+1);gray2=*ptrin;ptrin+=col_number;gray1=*ptrin;
	jc=add_to_jordan(jc,(double)i+0.5,(double)j+0.5,(char)6,gray1,gray2,(char)0);
	if (next_direction==4) break;
	gray2=gray1;ptrin--;gray1=*ptrin;
	jc=add_to_jordan(jc,(double)i+0.5,(double)j+0.5,(char)4,gray1,gray2,(char)1);
	if (next_direction==3) {add_to_frontier((double)i+0.5001,(double)j-0.4999);break;}
	gray2=gray1;ptrin--;gray1=*ptrin;
	jc=add_to_jordan(jc,(double)i+0.5,(double)j-0.5,(char)4,gray1,gray2,(char)0);
	if (next_direction==2) break;
	gray2=gray1;ptrin-=col_number;gray1=*ptrin;
	jc=add_to_jordan(jc,(double)i+0.5,(double)j-0.5,(char)2,gray1,gray2,(char)1);
	if (next_direction==1) {add_to_frontier((double)i-0.4999,(double)j-0.5001);break;}
	gray2=gray1;ptrin-=col_number;gray1=*ptrin;
	jc=add_to_jordan(jc,(double)i-0.5,(double)j-0.5,(char)2,gray1,gray2,(char)0);
	if (next_direction==0) break;
	gray2=gray1;ptrin++;gray1=*ptrin;
	jc=add_to_jordan(jc,(double)i-0.5,(double)j-0.5,(char)0,gray1,gray2,(char)1);
	add_to_frontier((double)i-0.5001,(double)j+0.4999); 
	/* This slight translation is used to obtain a simple curve */
	break; /* next_direction and previous_direction equal 7 */
      case 5:case 6:
	if (next_direction==4) {mwerror(FATAL,1,"Impossible 3\n");break;}
	if (previous_direction==6)
	  if (next_direction==5) {mwerror(FATAL,1,"Impossible 4\n");break;}
	/*jc=add_to_jordan(jc,(double)i+0.5,(double)j+0.5,4);*/ /* predecessor already made */
	if (next_direction==3) {add_to_frontier((double)i+0.5001,(double)j-0.4999);break;}
	ptrin=instart+((i+1)*col_number+j);gray2=*ptrin;ptrin--;gray1=*ptrin;
	jc=add_to_jordan(jc,(double)i+0.5,(double)j-0.5,(char)4,gray1,gray2,(char)0);
	if (next_direction==2) break;
	gray2=gray1;ptrin-=col_number;gray1=*ptrin;
	jc=add_to_jordan(jc,(double)i+0.5,(double)j-0.5,(char)2,gray1,gray2,(char)1);
	if (next_direction==1) {add_to_frontier((double)i-0.4999,(double)j-0.5001);break;}
	gray2=gray1;ptrin-=col_number;gray1=*ptrin;
	jc=add_to_jordan(jc,(double)i-0.5,(double)j-0.5,(char)2,gray1,gray2,(char)0);
	if (next_direction==0) break;
	gray2=gray1;ptrin++;gray1=*ptrin;
	jc=add_to_jordan(jc,(double)i-0.5,(double)j-0.5,(char)0,gray1,gray2,(char)1);
	if (next_direction==7) {add_to_frontier((double)i-0.5001,(double)j+0.4999);break;}
	gray2=gray1;ptrin++;gray1=*ptrin;
	jc=add_to_jordan(jc,(double)i-0.5,(double)j+0.5,(char)0,gray1,gray2,(char)0);
	if (next_direction==6) break;
	gray2=gray1;ptrin+=col_number;gray1=*ptrin;
	jc=add_to_jordan(jc,(double)i-0.5,(double)j+0.5,(char)6,gray1,gray2,(char)1);
	add_to_frontier((double)i+0.4999,(double)j+0.5001);
	break;  /* next_direction and previous_direction equal 5 */
      case 3:case 4:
	if (next_direction==2) {mwerror(FATAL,1,"Impossible 5\n");break;}
	if (previous_direction==4)
	  if (next_direction==3) {mwerror(FATAL,1,"Impossible 5\n");break;}
	/*jc=add_to_jordan(jc,(double)i+0.5,(double)j-0.5,2);*/ /* predecessor already made */
	if (next_direction==1) {add_to_frontier((double)i-0.4999,(double)j-0.5001);break;}
	ptrin=instart+(i*col_number+j-1);gray2=*ptrin;ptrin-=col_number;gray1=*ptrin;
	jc=add_to_jordan(jc,(double)i-0.5,(double)j-0.5,(char)2,gray1,gray2,(char)0);
	if (next_direction==0) break;
	gray2=*ptrin;ptrin++;gray1=*ptrin;
	jc=add_to_jordan(jc,(double)i-0.5,(double)j-0.5,(char)0,gray1,gray2,(char)1);
	if (next_direction==7) {add_to_frontier((double)i-0.5001,(double)j+0.4999);break;}
	gray2=*ptrin;ptrin++;gray1=*ptrin;
	jc=add_to_jordan(jc,(double)i-0.5,(double)j+0.5,(char)0,gray1,gray2,(char)0);
	if (next_direction==6) break;
	gray2=*ptrin;ptrin+=col_number;gray1=*ptrin;
	jc=add_to_jordan(jc,(double)i-0.5,(double)j+0.5,(char)6,gray1,gray2,(char)1);
	if (next_direction==5) {add_to_frontier((double)i+0.4999,(double)j+0.5001);break;}
	gray2=*ptrin;ptrin+=col_number;gray1=*ptrin;
	jc=add_to_jordan(jc,(double)i+0.5,(double)j+0.5,(char)6,gray1,gray2,(char)0);
	if (next_direction==4) break;
	gray2=*ptrin;ptrin--;gray1=*ptrin;
	jc=add_to_jordan(jc,(double)i+0.5,(double)j+0.5,(char)4,gray1,gray2,(char)1);
	add_to_frontier((double)i+0.5001,(double)j-0.4999);
	break;  /* next_direction and previous_direction equal 3 */
      case 1:case 2:
	if (next_direction==0) {mwerror(FATAL,1,"Impossible 7\n");break;}
 	if (previous_direction==2)
	  if (next_direction==1) {mwerror(FATAL,1,"Impossible 8\n");break;}
	/*jc=add_to_jordan(jc,(double)i-0.5,(double)j-0.5,0);*/ /* predecessor already made */
	if (next_direction==7) {add_to_frontier((double)i-0.5001,(double)j+0.4999);break;}
	ptrin=instart+((i-1)*col_number+j);gray2=*ptrin;ptrin++;gray1=*ptrin;
	jc=add_to_jordan(jc,(double)i-0.5,(double)j+0.5,(char)0,gray1,gray2,(char)0);
	if (next_direction==6) break;
	gray2=*ptrin;ptrin+=col_number;gray1=*ptrin;
	jc=add_to_jordan(jc,(double)i-0.5,(double)j+0.5,(char)6,gray1,gray2,(char)1);
	if (next_direction==5) {add_to_frontier((double)i+0.4999,(double)j+0.5001);break;}
	gray2=*ptrin;ptrin+=col_number;gray1=*ptrin;
	jc=add_to_jordan(jc,(double)i+0.5,(double)j+0.5,(char)6,gray1,gray2,(char)0);
	if (next_direction==4) break;
	gray2=*ptrin;ptrin--;gray1=*ptrin;
	jc=add_to_jordan(jc,(double)i+0.5,(double)j+0.5,(char)4,gray1,gray2,(char)1);
	if (next_direction==3) {add_to_frontier((double)i+0.5001,(double)j-0.4999);break;}
	gray2=*ptrin;ptrin--;gray1=*ptrin;
	jc=add_to_jordan(jc,(double)i+0.5,(double)j-0.5,(char)4,gray1,gray2,(char)0);
	if (next_direction==2) break;
	gray2=*ptrin;ptrin-=col_number;gray1=*ptrin;
	jc=add_to_jordan(jc,(double)i+0.5,(double)j-0.5,(char)2,gray1,gray2,(char)1);
	add_to_frontier((double)i-0.4999,(double)j-0.5001);
	break;  /* next_direction and previous_direction equal 1 */
      default:
	mwerror(FATAL,1,"Impossible global\n");
	break;
      }
    return jc;
  }
/***************************************************************************/

/* Here we scan image for computing level lines input and output, construct the
jordan curve, compute the optimal set of junctions and computes the final result */

void perform_disocclusion(ptrin,gray_value,ptr_label,value,ptrout)

register unsigned char *ptrin,*ptrout;
unsigned char gray_value; /* value of a pixel connected to the occlusion but out of it */
register float *ptr_label;
float value;

  {
    int i,j;
    int istart,jstart;
    char next_direction,previous_direction;
    char previous_direction_start,next_direction_start;
    jordan *jc,*jci;
    frontier *jbi;
    unsigned char *instart;
    float *label_start,*central_label;
    int level_number,new_level_number;
    unsigned char *level_sort,*ptr_sort,*level_array;
    unsigned char gray1,gray2;
    char neighbors; /* equal to 1 if the first point of the jordan curve has neighbors, 0 else */
    char encore;
    unsigned char ExternVal=0;
    int Line;
    int baire;
    double bx,by,norm_v,bvalue;
    Ppoly_t poly;
    int occlusionArea; /* Number of occluding points */
    int Number=0;

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

    label_start=ptr_label;
    instart=ptrin;
    lvalue=value;

    /* Coordinates of level lines refer to the first line of the set 
       (equal to globx) and to column 0 */
    occlusionArea=0;
    for (i=0;i<line_number-globx;i++){
      encore=0; /* if there is no pixel of occlusion on a line, it is unnecessary to 
		   study the next line */
      for (j=0;j<col_number;j++,ptr_label++)
	if (*ptr_label==value) {encore=1;occlusionArea++;}
      if (!encore) break;}

    /*
      We are going to compute the set of T-junctions, starting from the first pixel encountered 
      when the image is scanned from top to bottom and left to right. This point is used to 
      initialize the jordan curve. One looks for its predecessor and successor on the 
      occlusion boundary walked along clockwise, and computes all nodes associated to this pixel
      that could be input or output nodes for a level lines (there are at most 4 nodes for 
      each pixel, see figure below).
                            s1              s2
                     (i-0.5,j-0.5)  (i-0.5,j+0.5)
                                *----*
                                |    |
                                |    |