disocclusion.c

image processing including fourier,wavelet,segmentation etc.

    Ppoint_t eps1,eps2;
    int trii, ftrii, ltrii;
    int ei;
    pointnlink_t epnls[2], *lpnlp, *rpnlp, *pnlp;
    triangle_t *trip;
    int splitindex;
    double gdistance;
    double gdx1,gdy1,gdx2,gdy2;
    double v1x=jc1->vx,v1y=jc1->vy,v2x=jc2->vx,v2y=jc2->vy;
    double vx,vy,vxref,vyref; 

    eps1.x=jc1->x+globx;eps1.y=jc1->y;
    eps2.x=jc2->x+globx;eps2.y=jc2->y;
    
    /* find first and last triangles */
    for (trii = 0; trii < tril; trii++)
      if (pointintri (trii, &eps1))
	break;
    if (trii == tril) {
      prerror ("source point not in any triangle");
      return -1;}
    ftrii = trii;
    for (trii = 0; trii < tril; trii++)
      if (pointintri (trii, &eps2))
            break;
    if (trii == tril) {
      prerror ("destination point not in any triangle");
      return -1;}
    ltrii = trii;
    
    /* mark the strip of triangles from eps1 to eps2 */
    if (!marktripath (ftrii, ltrii)) {
      prerror ("cannot find triangle path");
      return -1;}

    /* if endpoints in same triangle, compute the euclidean distance */
    if (ftrii == ltrii) 
      {
	for (trii = 0; trii < tril; trii++) 
	  tris[trii].mark=0; /* unmark all triangles */
	vx=eps2.x-eps1.x;vy=eps2.y-eps1.y;
	*TVangle=fabs(atan2(v1x*vy-v1y*vx,v1x*vx+v1y*vy))+fabs(atan2(v2y*vx-v2x*vy,-v2x*vx-v2y*vy));
	return (Norm(vx,vy));
      }
    dq.fpnlpi = dq.pnlpn / 2, dq.lpnlpi = dq.fpnlpi - 1;
    /* build funnel and shortest path linked list (in add2dq) */
    epnls[0].pp = &eps1, epnls[0].link = NULL;
    epnls[1].pp = &eps2, epnls[1].link = NULL;
    add2dq (DQ_FRONT, &epnls[0]);
    dq.apex = dq.fpnlpi;
    trii = ftrii;
    while (trii != -1) {
      trip = &tris[trii];
      trip->mark = 2;
      
      /* find the left and right points of the exiting edge */
      for (ei = 0; ei < 3; ei++)
	if (trip->e[ei].rtp && trip->e[ei].rtp->mark == 1)
	  break;
      if (ei == 3) { /* in last triangle */
	if (ccw (&eps2, dq.pnlps[dq.fpnlpi]->pp,dq.pnlps[dq.lpnlpi]->pp) == ISCCW)
	  {lpnlp = dq.pnlps[dq.lpnlpi]; rpnlp = &epnls[1];}
	else
	  {lpnlp = &epnls[1]; rpnlp = dq.pnlps[dq.lpnlpi];}}
      else {
	pnlp = trip->e[(ei + 1) % 3].pnl1p;
	if (ccw (trip->e[ei].pnl0p->pp, pnlp->pp,trip->e[ei].pnl1p->pp) == ISCCW)
	  {lpnlp = trip->e[ei].pnl1p; rpnlp = trip->e[ei].pnl0p;}
	else
	  {lpnlp = trip->e[ei].pnl0p; rpnlp = trip->e[ei].pnl1p;}}

      /* update deque */
      if (trii == ftrii) {
	add2dq (DQ_BACK, lpnlp);
	add2dq (DQ_FRONT, rpnlp);}
      else {
	if (dq.pnlps[dq.fpnlpi] != rpnlp && dq.pnlps[dq.lpnlpi] != rpnlp) {
	  /* add right point to deque */
	  splitindex = finddqsplit (rpnlp);
	  splitdq (DQ_BACK, splitindex);
	  add2dq (DQ_FRONT, rpnlp);
	  /* if the split is behind the apex, then reset apex */
	  if (splitindex > dq.apex)
	    dq.apex = splitindex;}
	else {
	  /* add left point to deque */
	  splitindex = finddqsplit (lpnlp);
	  splitdq (DQ_FRONT, splitindex);
	  add2dq (DQ_BACK, lpnlp);
	  /* if the split is in front of the apex, then reset apex */
	  if (splitindex < dq.apex)
	    dq.apex = splitindex;}}
      
      trii = -1;
      for (ei = 0; ei < 3; ei++)
	if (trip->e[ei].rtp && trip->e[ei].rtp->mark == 1) {
	  trii = trip->e[ei].rtp - tris;
	  break;}
    }
     
    pnlp=&epnls[1];
    vxref=v2x;vyref=v2y;(*TVangle)=0.0;
    gdx1=pnlp->pp->x;gdy1=pnlp->pp->y;
    mwdebug("Path : %f %f\n",gdx1,gdy1);
    gdistance=0.0;
    while (pnlp->link)
      {
	pnlp=pnlp->link;
	gdx2=pnlp->pp->x;gdy2=pnlp->pp->y;
	vx=gdx2-gdx1;vy=gdy2-gdy1;
	gdistance+=Norm(vx,vy);
	(*TVangle)+=fabs(atan2(vxref*vy-vyref*vx,vxref*vx+vyref*vy));
	vxref=vx;vyref=vy;
	if (!(pnlp->link)){
	    vx=-v1x;vy=-v1y;
	    (*TVangle)+=fabs(atan2(vxref*vy-vyref*vx,vxref*vx+vyref*vy));
	    if ((eps1.x!=gdx2)||(eps1.y!=gdy2))
	      mwerror(FATAL,1,"Erreur dans GeodesicPath\n");}
	gdx1=gdx2;gdy1=gdy2;
	mwdebug("Path : %f %f\n",gdx1,gdy1);
      }
    for (trii = 0; trii < tril; trii++) 
      tris[trii].mark=0; /* unmark all triangles */

    return gdistance;
  }

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

int Pshortestpath (eps1,eps2,output)
Ppoint_t eps1,eps2;
Ppolyline_t *output;
  {
    int pi;
    int trii, ftrii, ltrii;
    int ei;
    pointnlink_t epnls[2], *lpnlp, *rpnlp, *pnlp;
    triangle_t *trip;
    int splitindex;
    
    /*-------------------------------------------------------
      The shortest path between T-junctions 'eps1' and 'eps2'
       is computed
    ---------------------------------------------------------*/
    
    /* find first and last triangles */
    for (trii = 0; trii < tril; trii++)
      if (pointintri (trii, &eps1))
	break;
    if (trii == tril) {
      prerror ("source point not in any triangle");
      return -1;}
    ftrii = trii;
    for (trii = 0; trii < tril; trii++)
      if (pointintri (trii, &eps2))
            break;
    if (trii == tril) {
      prerror ("destination point not in any triangle");
      return -1;}
    ltrii = trii;
    
    /* mark the strip of triangles from eps1 to eps2 */
    if (!marktripath (ftrii, ltrii)) {
      prerror ("cannot find triangle path");
      return -1;}
    
    /* if endpoints in same triangle, use a single line */
    if (ftrii == ltrii) 
      {
	growops (2);
	output->pn = 2;
	ops[0] = eps1, ops[1] = eps2;
	output->ps = ops;
	for (trii = 0; trii < tril; trii++) 
	  tris[trii].mark=0; /* unmark all triangles */
	return 0;
      }
    
    dq.fpnlpi = dq.pnlpn / 2, dq.lpnlpi = dq.fpnlpi - 1;
    /* build funnel and shortest path linked list (in add2dq) */
    epnls[0].pp = &eps1, epnls[0].link = NULL;
    epnls[1].pp = &eps2, epnls[1].link = NULL;
    add2dq (DQ_FRONT, &epnls[0]);
    dq.apex = dq.fpnlpi;
    trii = ftrii;
    while (trii != -1) {
      trip = &tris[trii];
      trip->mark = 2;
      
      /* find the left and right points of the exiting edge */
      for (ei = 0; ei < 3; ei++)
	if (trip->e[ei].rtp && trip->e[ei].rtp->mark == 1)
	  break;
      if (ei == 3) { /* in last triangle */
	if (ccw (&eps2, dq.pnlps[dq.fpnlpi]->pp,
		 dq.pnlps[dq.lpnlpi]->pp) == ISCCW)
	  {lpnlp = dq.pnlps[dq.lpnlpi]; rpnlp = &epnls[1];}
	else
	  {lpnlp = &epnls[1]; rpnlp = dq.pnlps[dq.lpnlpi];}}
      else {
	pnlp = trip->e[(ei + 1) % 3].pnl1p;
	if (ccw (trip->e[ei].pnl0p->pp, pnlp->pp,
		 trip->e[ei].pnl1p->pp) == ISCCW)
	  {lpnlp = trip->e[ei].pnl1p; rpnlp = trip->e[ei].pnl0p;}
	else
	  {lpnlp = trip->e[ei].pnl0p; rpnlp = trip->e[ei].pnl1p;}}
      
      /* update deque */
      if (trii == ftrii) {
	add2dq (DQ_BACK, lpnlp);
	add2dq (DQ_FRONT, rpnlp);}
      else {
	if (dq.pnlps[dq.fpnlpi] != rpnlp && dq.pnlps[dq.lpnlpi] != rpnlp) {
	  /* add right point to deque */
	  splitindex = finddqsplit (rpnlp);
	  splitdq (DQ_BACK, splitindex);
	  add2dq (DQ_FRONT, rpnlp);
	  /* if the split is behind the apex, then reset apex */
	  if (splitindex > dq.apex)
	    dq.apex = splitindex;}
	else {
	  /* add left point to deque */
	  splitindex = finddqsplit (lpnlp);
	  splitdq (DQ_FRONT, splitindex);
	  add2dq (DQ_BACK, lpnlp);
	  /* if the split is in front of the apex, then reset apex */
	  if (splitindex < dq.apex)
	    dq.apex = splitindex;}}
      
      trii = -1;
      for (ei = 0; ei < 3; ei++)
	if (trip->e[ei].rtp && trip->e[ei].rtp->mark == 1) {
	  trii = trip->e[ei].rtp - tris;
	  break;}
    }
    
    mwdebug("polypath");
    for (pnlp = &epnls[1]; pnlp; pnlp = pnlp->link)
      mwdebug("%d %d  %f %f\n",pnlp,pnlp->pp,pnlp->pp->x, pnlp->pp->y);
    mwdebug("\n\n");
    
    for (pi = 0, pnlp = &epnls[1]; pnlp; pnlp = pnlp->link)  pi++;
    growops(pi);
    output->pn = pi;
    for (pi = pi - 1, pnlp = &epnls[1]; pnlp; pi--, pnlp = pnlp->link)
      ops[pi] = *pnlp->pp;
    output->ps = ops;
    
    for (trii = 0; trii < tril; trii++) 
      tris[trii].mark=0; /* unmark all triangles */
    
    return 0;
}
/****************************************************************************/

void FreeTPath()
  {
    if (tris) {free((void*)tris);tris=NULL;}
    if (pnls) {free((void*)pnls);pnls=NULL;}
    if (pnlps) {free((void*)pnlps);pnlps=NULL;}
    if (dq.pnlps) {free((void*)(dq.pnlps));dq.pnlps=NULL;}
  }

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

Ppoint_t *polypoints; 
    /* Copy of the frontier as a polygonal line whose structure 
       is compatible with Mitchell's shortest path algorithm */
Ppoint_t *endTjunctions; /* All the geodesics endpoints */

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

/* Adds a new element to the list of vertices of occlusion boundary */

void add_to_frontier(x,y)
double x,y;

  {
    frontier *inter;
  
    inter=jb;
    jb=(frontier*)malloc((size_t)sizeof(frontier));
    if (!jb) mwerror(FATAL,1,"Not enough memory (1)!\n");
    numfrontier++;
    jb->x=x;jb->y=y;
    jb->next=inter;
    if (inter) inter->previous=jb;
    jb->previous=(frontier*)NULL;
  }
/****************************************************************************/

/* Adds a new element to the list of T-junctions if current point is a T-junction
   and in any case, adds a new element to the list of occlusion frontier points
   Note that level line direction is at this level of algorithm either equal to 
   0, 2, 4 or 6 (see convention above) */

jordan *add_to_jordan(nouv,x,y,direction,gray1,gray2,dofrontier)
jordan *nouv;
double x,y;
char direction;
unsigned char gray1,gray2;
char dofrontier;

  {
    jordan *inter;
  
    if (gray1!=gray2){
      inter=nouv;
      nouv=(jordan*)malloc((size_t)sizeof(jordan));
      if (!nouv) mwerror(FATAL,1,"Not enough memory (2)!\n");
      nouv->x=x;nouv->y=y;nouv->direction=direction;
      nouv->gray1=gray1;nouv->gray2=gray2;
      nouv->junction=(jordan*)NULL;
      switch (direction){
      case 0:nouv->vx=1;nouv->vy=0;break;
      case 2:nouv->vx=0;nouv->vy=1;break;
      case 4:nouv->vx=(-1);nouv->vy=0;break;
      case 6:nouv->vx=0;nouv->vy=(-1);break;  
      default:mwerror(FATAL,1,"Impossible dans add_to_jordan\n");break;}
      /* Remember convention : x denotes vertical axis and y horizontal axis */
      nouv->next=inter;
      if (inter) inter->previous=nouv;
      nouv->previous=(jordan*)NULL;}
    if (dofrontier) add_to_frontier(x,y);
    /* dofrontier is useful to avoid redundancy : recall that a point can be associated 
       with several level lines. It is also used for taking only those points which  are vertices
       of the polygonal line describing the occlusion boundary */
    return nouv;
  }

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

/* Inserts a new element in the chain defining the structuring element*/

jordan *insert_jordan(jc,gray1,gray2)
jordan *jc;
unsigned char gray1,gray2;

  {
    jordan *new;  /* new is inserted between jc and jc->next */
  
    new=(jordan*)malloc((size_t)sizeof(jordan));
    if (!new) mwerror(FATAL,1,"Not enough memory (3)!\n");
    new->x=jc->x;new->y=jc->y;
    new->vx=jc->vx;new->vy=jc->vy;
    new->gray1=gray1;new->gray2=gray2;
    new->next=jc->next;
    if (jc->next) jc->next->previous=new;
    new->direction=jc->direction;
    new->junction=(jordan*)NULL;
    jc->next=new;
    new->previous=jc;
    return new;
  }
/****************************************************************************/

/* Deletes an element in the chain defining the structuring element*/

jordan *delete_jordan(jc)
jordan *jc;

  {
    if (jc->previous) jc->previous->next=jc->next;
    if (jc->next) jc->next->previous=jc->previous;
    free((void*)jc);
    return (jordan*)NULL;
  }

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

/*  Frees the chain defining the T-junctions */
/*  Be careful that this function is valid only if jordan curve is "broken" */

void free_jordan(jc)
jordan *jc;
  
  {
    if (jc->next) free_jordan(jc->next);
    free((void*)jc);
    jc=(jordan*)NULL;
  }

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

/*  Frees the frontier chain  */
/*  Be careful that this function is valid only if frontier curve is "broken" */

void free_frontier(jbi)
frontier *jbi;
  
  {
    if (jbi->next) free_frontier(jbi->next);
    free((void*)jbi);