disocclusion.c

image processing including fourier,wavelet,segmentation etc.

/* is pbp between pap and pcp */
int between (pap,pbp,pcp)
Ppoint_t *pap,*pbp,*pcp;
  {
    Ppoint_t p1, p2;

    p1.x = pbp->x - pap->x, p1.y = pbp->y - pap->y;
    p2.x = pcp->x - pap->x, p2.y = pcp->y - pap->y;
    if (ccw (pap, pbp, pcp) != ISON)
        return FALSE;
    return (p2.x * p1.x + p2.y * p1.y >= 0) &&
            (p2.x * p2.x + p2.y * p2.y <= p1.x * p1.x + p1.y * p1.y);
  }

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

/* line to line intersection */
int intersects (pap,pbp,pcp,pdp)
Ppoint_t *pap,*pbp,*pcp,*pdp;
  {
    int ccw1, ccw2, ccw3, ccw4;

    if (ccw (pap, pbp, pcp) == ISON || ccw (pap, pbp, pdp) == ISON ||
            ccw (pcp, pdp, pap) == ISON || ccw (pcp, pdp, pbp) == ISON) {
        if (between (pap, pbp, pcp) || between (pap, pbp, pdp) ||
                between (pcp, pdp, pap) || between (pcp, pdp, pbp))
            return TRUE;
    } else {
        ccw1 = (ccw (pap, pbp, pcp) == ISCCW) ? 1 : 0;
        ccw2 = (ccw (pap, pbp, pdp) == ISCCW) ? 1 : 0;
        ccw3 = (ccw (pcp, pdp, pap) == ISCCW) ? 1 : 0;
        ccw4 = (ccw (pcp, pdp, pbp) == ISCCW) ? 1 : 0;
        return (ccw1 ^ ccw2) && (ccw3 ^ ccw4);
    }
    return FALSE;
  }

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

/* check if (i, i + 2) is a diagonal */
int isdiagonal (pnli,pnlip2,pnln)
int pnli,pnlip2,pnln;
  {
    int pnlip1, pnlim1, pnlj, pnljp1, res;

    /* neighborhood test */
    pnlip1 = (pnli + 1) % pnln;
    pnlim1 = (pnli + pnln - 1) % pnln;
    /* If P[pnli] is a convex vertex [ pnli+1 left of (pnli-1,pnli) ]. */
    if (ccw (pnlps[pnlim1]->pp, pnlps[pnli]->pp, pnlps[pnlip1]->pp) == ISCCW)
        res = (ccw (pnlps[pnli]->pp, pnlps[pnlip2]->pp,
                pnlps[pnlim1]->pp) == ISCCW) &&
                (ccw (pnlps[pnlip2]->pp, pnlps[pnli]->pp,
                pnlps[pnlip1]->pp) == ISCCW);
    /* Assume (pnli - 1, pnli, pnli + 1) not collinear. */
    else
        res = (ccw (pnlps[pnli]->pp, pnlps[pnlip2]->pp,
                pnlps[pnlip1]->pp) == ISCW);
    if (!res)
        return FALSE;

    /* check against all other edges */
    for (pnlj = 0; pnlj < pnln; pnlj++) {
        pnljp1 = (pnlj + 1) % pnln;
        if (!((pnlj == pnli) || (pnljp1 == pnli) ||
                (pnlj == pnlip2) || (pnljp1 == pnlip2)))
            if (intersects (pnlps[pnli]->pp, pnlps[pnlip2]->pp,
                    pnlps[pnlj]->pp, pnlps[pnljp1]->pp))
                return FALSE;
    }
    return TRUE;
  }

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

void growtris (newtrin)
int newtrin;
  {
    if (newtrin <= trin) return;
    if (!tris) {
        if (!(tris = (triangle_t *) malloc (TRIANGLESIZE * newtrin))) {
            prerror ("cannot malloc tris");
            abort ();
        }
    } else {
        if (!(tris = (triangle_t *) realloc ((void *) tris,
                TRIANGLESIZE * newtrin))) {
            prerror ("cannot realloc tris");
            abort ();
        }
    }
    trin = newtrin;
  }

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

void loadtriangle (pnlap,pnlbp,pnlcp)
pointnlink_t *pnlap,*pnlbp,*pnlcp;
  {
    triangle_t *trip;
    int ei;

    /* make space */
    if (tril >= trin) growtris (trin + 20);
    trip = &tris[tril++];
    trip->mark = 0;
    trip->e[0].pnl0p = pnlap, trip->e[0].pnl1p = pnlbp, trip->e[0].rtp = NULL;
    trip->e[1].pnl0p = pnlbp, trip->e[1].pnl1p = pnlcp, trip->e[1].rtp = NULL;
    trip->e[2].pnl0p = pnlcp, trip->e[2].pnl1p = pnlap, trip->e[2].rtp = NULL;
    for (ei = 0; ei < 3; ei++)
        trip->e[ei].ltp = trip;
}

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

/* triangulate polygon */
void triangulate (pnln)
int pnln; 
  {
    int pnli, pnlip1, pnlip2;

    if (pnln > 3) {
        for (pnli = 0; pnli < pnln; pnli++) {
            pnlip1 = (pnli + 1) % pnln;
            pnlip2 = (pnli + 2) % pnln;
	    if (isdiagonal (pnli, pnlip2,  pnln)) {
		loadtriangle (pnlps[pnli], pnlps[pnlip1], pnlps[pnlip2]);
		for (pnli = pnlip1; pnli < pnln - 1; pnli++)
                    pnlps[pnli] = pnlps[pnli + 1];
                triangulate (pnln - 1);
                return;
            }
        }
	prerror ("the polygon frontier is not a simple curve");
        abort ();
    } else
        loadtriangle (pnlps[0], pnlps[1], pnlps[2]);
  }

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

/* connect a pair of triangles at their common edge (if any) */
void connecttris (tri1,tri2)
int tri1,tri2;
  {
    triangle_t *tri1p, *tri2p;
    int ei, ej;

    for (ei = 0 ; ei < 3; ei++) {
        for (ej = 0; ej < 3; ej++) {
            tri1p = &tris[tri1], tri2p = &tris[tri2];
            if ((tri1p->e[ei].pnl0p->pp == tri2p->e[ej].pnl0p->pp &&
                    tri1p->e[ei].pnl1p->pp == tri2p->e[ej].pnl1p->pp) ||
                    (tri1p->e[ei].pnl0p->pp == tri2p->e[ej].pnl1p->pp &&
                    tri1p->e[ei].pnl1p->pp == tri2p->e[ej].pnl0p->pp))
                tri1p->e[ei].rtp = tri2p, tri2p->e[ej].rtp = tri1p;
        }
    }
  }

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

/* find and mark path from trii, to trij */
int marktripath (trii,trij)
int trii,trij;
  {
    int ei;

    if (tris[trii].mark)
        return FALSE;
    tris[trii].mark = 1;
    if (trii == trij)
        return TRUE;
    for (ei = 0; ei < 3; ei++)
        if (tris[trii].e[ei].rtp &&
                marktripath (tris[trii].e[ei].rtp - tris, trij))
            return TRUE;
    tris[trii].mark = 0;
    return FALSE;
}

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

/* add a new point to the deque, either front or back */
void add2dq (side,pnlp)
int side;
pointnlink_t *pnlp;
  {
    if (side == DQ_FRONT) {
        if (dq.lpnlpi - dq.fpnlpi >= 0)
            pnlp->link = dq.pnlps[dq.fpnlpi]; /* shortest path links */
        dq.fpnlpi--;
        dq.pnlps[dq.fpnlpi] = pnlp;
    } else {
        if (dq.lpnlpi - dq.fpnlpi >= 0)
            pnlp->link = dq.pnlps[dq.lpnlpi]; /* shortest path links */
        dq.lpnlpi++;
        dq.pnlps[dq.lpnlpi] = pnlp;
    }
  }

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

void splitdq (side,index)
int side,index;
 {
    if (side == DQ_FRONT)
        dq.lpnlpi = index;
    else
        dq.fpnlpi = index;
 }

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

int finddqsplit (pnlp)
pointnlink_t *pnlp;
  {
    int index;

    for (index = dq.fpnlpi; index < dq.apex; index++)
        if (ccw (dq.pnlps[index + 1]->pp, dq.pnlps[index]->pp, pnlp->pp) == ISCCW)
            return index;
    for (index = dq.lpnlpi; index > dq.apex; index--)
        if (ccw (dq.pnlps[index - 1]->pp, dq.pnlps[index]->pp, pnlp->pp) == ISCW)
            return index;
    return dq.apex;
  }

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

int pointintri (trii,pp)
int trii;
Ppoint_t *pp;
  {
    int ei, sum;

    for (ei = 0, sum = 0; ei < 3; ei++)
      if (ccw (tris[trii].e[ei].pnl0p->pp,
	       tris[trii].e[ei].pnl1p->pp, pp) != ISCW)
	sum++;
    return (sum == 3 || sum == 0);
  } 

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

void growpnls (newpnln)
int newpnln;
  {
    if (!pnls) {
      if (!(pnls = (pointnlink_t *) malloc (POINTNLINKSIZE * newpnln))) {
	prerror ("cannot malloc pnls");
	abort ();}
      if (!(pnlps = (pointnlink_t **) malloc (POINTNLINKPSIZE * newpnln))) {
	prerror ("cannot malloc pnlps");
	abort ();}}
  }

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

void growdq (newdqn)
int newdqn;
  {
    if (newdqn <= dq.pnlpn) return;
    if (!dq.pnlps) {
      dq.pnlps = (pointnlink_t **)malloc(POINTNLINKPSIZE * newdqn);
      if (!dq.pnlps){
	prerror ("cannot malloc dq.pnlps");
	abort ();}}
    else
      if (!(dq.pnlps = (pointnlink_t **)realloc ((void *) dq.pnlps,
						 POINTNLINKPSIZE * newdqn))) {
	prerror ("cannot realloc dq.pnlps");
	abort ();}
    dq.pnlpn = newdqn;
  }

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

void growops (newopn)
int newopn;
  {
    ops=NULL;
    if (!(ops = (Ppoint_t *) malloc (POINTSIZE * newopn))) {
      prerror ("cannot malloc ops");
      abort ();}
  }

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

void Triangulation(polyp)
Ppoly_t *polyp;

  {
    int pi, minpi;
    double minx;
    Ppoint_t p1, p2, p3;
    int trii, trij;
    int ei;
    int pnli;
   
    /* make space */
    growpnls (polyp->pn);
    pnll = 0;
    tril = 0;
	
    /* make sure polygon is counterclockwise and load pnls array */
    for (pi = 0, minx = 1e+200, minpi = -1; pi < polyp->pn; pi++) {
        if (minx > polyp->ps[pi].x)
            minx = polyp->ps[pi].x, minpi = pi;
    }
    p2 = polyp->ps[minpi];
    p1 = polyp->ps[((minpi == 0) ? polyp->pn - 1: minpi - 1)];
    p3 = polyp->ps[((minpi == polyp->pn - 1) ? 0 : minpi + 1)];
    if (((p1.x == p2.x && p2.x == p3.x) && (p3.y > p2.y)) ||
            ccw (&p1, &p2, &p3) != ISCCW) {
        for (pi = polyp->pn - 1; pi >= 0; pi--) {
            if (pi < polyp->pn - 1 && polyp->ps[pi].x == polyp->ps[pi + 1].x &&
                    polyp->ps[pi].y == polyp->ps[pi + 1].y)
                continue;
            pnls[pnll].pp = &polyp->ps[pi];
            pnls[pnll].link = &pnls[pnll % polyp->pn];
            pnlps[pnll] = &pnls[pnll];
            pnll++;
        }
    } else {
        for (pi = 0; pi < polyp->pn; pi++) {
            if (pi > 0 && polyp->ps[pi].x == polyp->ps[pi - 1].x &&
                    polyp->ps[pi].y == polyp->ps[pi - 1].y)
                continue;
            pnls[pnll].pp = &polyp->ps[pi];
            pnls[pnll].link = &pnls[pnll % polyp->pn];
            pnlps[pnll] = &pnls[pnll];
            pnll++;
        }
    }

    mwdebug("points\n%d\n", pnll);
    for (pnli = 0; pnli < pnll; pnli++)
        mwdebug("%f %f\n", pnls[pnli].pp->x, pnls[pnli].pp->y);

    /* generate list of triangles */
    trin=0;
    triangulate (pnll);
    mwdebug("triangles\n%d\n", tril);
    for (trii = 0; trii < tril; trii++)
        for (ei = 0; ei < 3; ei++)
            mwdebug("%f %f\n", tris[trii].e[ei].pnl0p->pp->x,
                    tris[trii].e[ei].pnl0p->pp->y);

    /* connect all pairs of triangles that share an edge */
    for (trii = 0; trii < tril; trii++)
        for (trij = trii + 1; trij < tril; trij++)
            connecttris (trii, trij);
    dq.pnlpn=0;
    dq.pnlps=NULL;
    growdq (polyp->pn * 2);

  }
/****************************************************************************/
/* return an angle between -pi and pi */

double AngleFix(a)
double a;
{
    while (a > M_PI)
        a -= 2 * M_PI;
    while (a < -M_PI)
        a += 2 * M_PI;
    return a;
}
/****************************************************************************/

/* compare the current angle a with reference angle b */

double AngleCheck(a,b)
double a,b;
{
    if ((a+M_PI>=b+M_PI_2)&&(a+M_PI<=b+3*M_PI_2)) return a;
    return AngleFix(a+M_PI);
}


/****************************************************************************/
/* eps1 and eps2 are the T-junctions. 
   meandir1 is the geodesic curve mean direction at eps1
   meandir2 is the geodesic curve mean direction at eps2 */

double GeodesicDistance(jc1,jc2,TVangle)
jordan *jc1,*jc2;
double *TVangle;
  {