snaphu_tile.c

phase unwrapping algorithm for SAR interferometry

    scndryarcs[tilenum][i].to
      =&scndrynodes[tempnode->level][tempnode->incost];
    to=scndryarcs[tilenum][i].to;

    /* update secondary arc pointers in nodesupp strcutres */
    fromsupp=&nodesupp[from->row][from->col];
    j=0;
    while(fromsupp->neighbornodes[j]!=to){
      j++;
    }
    fromsupp->outarcs[j]=&scndryarcs[tilenum][i];
    tosupp=&nodesupp[to->row][to->col];
    j=0;
    while(tosupp->neighbornodes[j]!=from){
      j++;
    }
    tosupp->outarcs[j]=&scndryarcs[tilenum][i];
  }

  /* set outputs */
  nscndrynodes[tilenum]=nnewnodes;
  nscndryarcs[tilenum]=nnewarcs;
  totarclens[tilenum]=totarclen;

  /* free memory */
  Free2DArray((void **)primarynodes,nnrow);
  Free2DArray((void **)flows,2*nrow-1);
  Free2DArray((void **)rightedgeflows,nrow);
  Free2DArray((void **)leftedgeflows,nrow);
  Free2DArray((void **)upperedgeflows,1);
  Free2DArray((void **)loweredgeflows,1);
  Free2DArray((void **)rightedgecosts,nrow);
  Free2DArray((void **)leftedgecosts,nrow);
  Free2DArray((void **)upperedgecosts,1);
  Free2DArray((void **)loweredgecosts,1);
}


/* function: FindNumPathsOut()
 * ---------------------------
 * Check all outgoing arcs to see how many paths out there are. 
 */
long FindNumPathsOut(nodeT *from, paramT *params, long tilerow, long tilecol, 
		     long nnrow, long nncol, short **regions, 
		     short **nextregions, short **lastregions,
		     short *regionsabove, short *regionsbelow, long prevncol){

  long npathsout, ntilerow, ntilecol, fromrow, fromcol;

  /* initialize */
  ntilerow=params->ntilerow;
  ntilecol=params->ntilecol;
  fromrow=from->row;
  fromcol=from->col;
  npathsout=0;

  /* rightward arc */
  if(fromcol!=nncol-1){
    if(fromrow==0 || fromrow==nnrow-1 
       || regions[fromrow-1][fromcol]!=regions[fromrow][fromcol]){
      npathsout++;
    }
  }else{
    if(fromrow==0 || fromrow==nnrow-1 || 
       (tilecol!=ntilecol-1
        && nextregions[fromrow-1][0]!=nextregions[fromrow][0])){
      npathsout++;
    }
  }

  /* downward arc */
  if(fromrow!=nnrow-1){
    if(fromcol==0 || fromcol==nncol-1
       || regions[fromrow][fromcol]!=regions[fromrow][fromcol-1]){
      npathsout++;
    }
  }else{
    if(fromcol==0 || fromcol==nncol-1 ||
       (tilerow!=ntilerow-1 
        && regionsbelow[fromcol]!=regionsbelow[fromcol-1])){
      npathsout++;
    }
  }

  /* leftward arc */
  if(fromcol!=0){
    if(fromrow==0 || fromrow==nnrow-1 
       || regions[fromrow][fromcol-1]!=regions[fromrow-1][fromcol-1]){
      npathsout++;
    }
  }else{
    if(fromrow==0 || fromrow==nnrow-1 || 
       (tilecol!=0
        && (lastregions[fromrow][prevncol-1]
            !=lastregions[fromrow-1][prevncol-1]))){
      npathsout++;
    }
  }

  /* upward arc */
  if(fromrow!=0){
    if(fromcol==0 || fromcol==nncol-1
       || regions[fromrow-1][fromcol-1]!=regions[fromrow-1][fromcol]){
      npathsout++;
    }
  }else{
    if(fromcol==0 || fromcol==nncol-1 ||
       (tilerow!=0
        && regionsabove[fromcol-1]!=regionsabove[fromcol])){
      npathsout++;
    }
  }

  /* return number of paths out of node */
  return(npathsout);

}


/* function: RegionTraceCheckNeighbors()
 * -------------------------------------
 */
void RegionTraceCheckNeighbors(nodeT *from, nodeT **nextnodeptr, 
			       nodeT **primarynodes, short **regions, 
			       short **nextregions, short **lastregions, 
			       short *regionsabove, short *regionsbelow,  
			       long tilerow, long tilecol, long nnrow, 
			       long nncol, nodeT **scndrynodes, 
			       nodesuppT **nodesupp, scndryarcT **scndryarcs, 
			       long *nnewnodesptr, long *nnewarcsptr, 
			       long flowmax, long nrow, long ncol, 
			       long prevnrow, long prevncol, paramT *params, 
			       void **costs, void **rightedgecosts, 
			       void **loweredgecosts, void **leftedgecosts, 
			       void **upperedgecosts, short **flows, 
			       short **rightedgeflows, short **loweredgeflows,
			       short **leftedgeflows, short **upperedgeflows,
			       long ***scndrycosts, 
			       nodeT ***updatednontilenodesptr, 
			       long *nupdatednontilenodesptr, 
			       long *updatednontilenodesizeptr,
			       short **inontilenodeoutarcptr, 
			       long *totarclenptr){

  long fromrow, fromcol;
  nodeT *to, *nextnode;


  /* initialize */
  fromrow=from->row;
  fromcol=from->col;
  nextnode=(*nextnodeptr);


  /* check rightward arc */
  if(fromcol!=nncol-1){
    to=&primarynodes[fromrow][fromcol+1];
    if(fromrow==0 || fromrow==nnrow-1 
       || regions[fromrow-1][fromcol]!=regions[fromrow][fromcol]){
      if(to!=from->pred){
	to->pred=from;
	if(to->group==NOTINBUCKET){
	  to->group=INBUCKET;
	  to->next=nextnode;
	  nextnode=to;
	}else if(to->group==ONTREE && (fromrow!=0 || tilerow==0)){
	  TraceSecondaryArc(to,scndrynodes,nodesupp,scndryarcs,scndrycosts,
			    nnewnodesptr,nnewarcsptr,tilerow,tilecol,flowmax,
			    nrow,ncol,prevnrow,prevncol,params,costs,
			    rightedgecosts,loweredgecosts,leftedgecosts,
			    upperedgecosts,flows,rightedgeflows,
			    loweredgeflows,leftedgeflows,upperedgeflows,
			    updatednontilenodesptr,nupdatednontilenodesptr,
			    updatednontilenodesizeptr,inontilenodeoutarcptr,
			    totarclenptr);
	}
      }
    }
  }


  /* check downward arc */
  if(fromrow!=nnrow-1){
    to=&primarynodes[fromrow+1][fromcol];
    if(fromcol==0 || fromcol==nncol-1
       || regions[fromrow][fromcol]!=regions[fromrow][fromcol-1]){
      if(to!=from->pred){
	to->pred=from;
	if(to->group==NOTINBUCKET){
	  to->group=INBUCKET;
	  to->next=nextnode;
	  nextnode=to;
	}else if(to->group==ONTREE && (fromcol!=0 || tilecol==0)){
	  TraceSecondaryArc(to,scndrynodes,nodesupp,scndryarcs,scndrycosts,
			    nnewnodesptr,nnewarcsptr,tilerow,tilecol,flowmax,
			    nrow,ncol,prevnrow,prevncol,params,costs,
			    rightedgecosts,loweredgecosts,leftedgecosts,
			    upperedgecosts,flows,rightedgeflows,
			    loweredgeflows,leftedgeflows,upperedgeflows,
			    updatednontilenodesptr,nupdatednontilenodesptr,
			    updatednontilenodesizeptr,inontilenodeoutarcptr,
			    totarclenptr);
	}
      }
    }
  }

      
  /* check leftward arc */
  if(fromcol!=0){
    to=&primarynodes[fromrow][fromcol-1];
    if(fromrow==0 || fromrow==nnrow-1 
       || regions[fromrow][fromcol-1]!=regions[fromrow-1][fromcol-1]){
      if(to!=from->pred){
	to->pred=from;
	if(to->group==NOTINBUCKET){
	  to->group=INBUCKET;
	  to->next=nextnode;
	  nextnode=to;
	}else if(to->group==ONTREE && (fromrow!=0 || tilerow==0)){
	  TraceSecondaryArc(to,scndrynodes,nodesupp,scndryarcs,scndrycosts,
			    nnewnodesptr,nnewarcsptr,tilerow,tilecol,flowmax,
			    nrow,ncol,prevnrow,prevncol,params,costs,
			    rightedgecosts,loweredgecosts,leftedgecosts,
			    upperedgecosts,flows,rightedgeflows,
			    loweredgeflows,leftedgeflows,upperedgeflows,
			    updatednontilenodesptr,nupdatednontilenodesptr,
			    updatednontilenodesizeptr,inontilenodeoutarcptr,
			    totarclenptr);
	}
      }
    }
  }


  /* check upward arc */
  if(fromrow!=0){
    to=&primarynodes[fromrow-1][fromcol];
    if(fromcol==0 || fromcol==nncol-1
       || regions[fromrow-1][fromcol-1]!=regions[fromrow-1][fromcol]){
      if(to!=from->pred){
	to->pred=from;
	if(to->group==NOTINBUCKET){
	  to->group=INBUCKET;
	  to->next=nextnode;
	  nextnode=to;
	}else if(to->group==ONTREE && (fromcol!=0 || tilecol==0)){
	  TraceSecondaryArc(to,scndrynodes,nodesupp,scndryarcs,scndrycosts,
			    nnewnodesptr,nnewarcsptr,tilerow,tilecol,flowmax,
			    nrow,ncol,prevnrow,prevncol,params,costs,
			    rightedgecosts,loweredgecosts,leftedgecosts,
			    upperedgecosts,flows,rightedgeflows,
			    loweredgeflows,leftedgeflows,upperedgeflows,
			    updatednontilenodesptr,nupdatednontilenodesptr,
			    updatednontilenodesizeptr,inontilenodeoutarcptr,
			    totarclenptr);
	}
      }
    }
  }


  /* set return values */
  *nextnodeptr=nextnode;

}


/* function: SetUpperEdge()
 * ------------------------
 */
void SetUpperEdge(long ncol, long tilerow, long tilecol, void **voidcosts, 
		  void *voidcostsabove, float **unwphase, 
		  float *unwphaseabove, void **voidupperedgecosts, 
		  short **upperedgeflows, paramT *params, short **bulkoffsets){

  long col, reloffset;
  double dphi, dpsi;
  costT **upperedgecosts, **costs, *costsabove;
  smoothcostT **upperedgesmoothcosts, **smoothcosts, *smoothcostsabove;
  long nshortcycle;


  /* typecast generic pointers to costT pointers */
  upperedgecosts=(costT **)voidupperedgecosts;
  costs=(costT **)voidcosts;
  costsabove=(costT *)voidcostsabove;
  upperedgesmoothcosts=(smoothcostT **)voidupperedgecosts;
  smoothcosts=(smoothcostT **)voidcosts;
  smoothcostsabove=(smoothcostT *)voidcostsabove;
  
  /* see if tile is in top row */
  if(tilerow!=0){

    /* set up */
    nshortcycle=params->nshortcycle;
    reloffset=bulkoffsets[tilerow-1][tilecol]-bulkoffsets[tilerow][tilecol];

    /* loop over all arcs on the boundary */
    for(col=0;col<ncol;col++){
      dphi=(unwphaseabove[col]-unwphase[0][col])/TWOPI;
      upperedgeflows[0][col]=(short )LRound(dphi)-reloffset;
      dpsi=dphi-floor(dphi);
      if(dpsi>0.5){
	dpsi-=1.0;
      }
      if(params->costmode==TOPO || params->costmode==DEFO){
	upperedgecosts[0][col].offset=nshortcycle*dpsi;
	upperedgecosts[0][col].sigsq=ceil((costs[0][col].sigsq
					   +costsabove[col].sigsq)/2.0);
	if(costs[0][col].dzmax>costsabove[col].dzmax){
	  upperedgecosts[0][col].dzmax=costs[0][col].dzmax;
	}else{
	  upperedgecosts[0][col].dzmax=costsabove[col].dzmax;
	}
	if(costs[0][col].laycost<costsabove[col].laycost){
	  upperedgecosts[0][col].laycost=costs[0][col].laycost;
	}else{
	  upperedgecosts[0][col].laycost=costsabove[col].laycost;
	}
      }else if(params->costmode==SMOOTH){
	upperedgesmoothcosts[0][col].offset=nshortcycle*dpsi;
	upperedgesmoothcosts[0][col].sigsq=
	  ceil((smoothcosts[0][col].sigsq+smoothcostsabove[col].sigsq)/2.0);
      }else{
	fprintf(sp0,"Illegal cost mode in SetUpperEdge().  This is a bug.\n");
	exit(ABNORMAL_EXIT);
      }
    }

  }else{
    if(params->costmode==TOPO || params->costmode==DEFO){
      for(col=0;col<ncol;col++){
	upperedgecosts[0][col].offset=LARGESHORT/2;
	upperedgecosts[0][col].sigsq=LARGESHORT;
	upperedgecosts[0][col].dzmax=LARGESHORT;
	upperedgecosts[0][col].laycost=0;
      }
    }else if(params->costmode==SMOOTH){
      for(col=0;col<ncol;col++){
	upperedgesmoothcosts[0][col].offset=0;
	upperedgesmoothcosts[0][col].sigsq=LARGESHORT;
      }
    }else{
      fprintf(sp0,"Illegal cost mode in SetUpperEdge().  This is a bug.\n");
      exit(ABNORMAL_EXIT);
    }
  }
}


/* function: SetLowerEdge()
 * ------------------------
 */
void SetLowerEdge(long nrow, long ncol, long tilerow, long tilecol, 
		  void **voidcosts, void *voidcostsbelow, 
		  float **unwphase, float *unwphasebelow, 
		  void **voidloweredgecosts, short **loweredgeflows, 
		  paramT *params, short **bulkoffsets){

  long *flowhistogram;
  long col, iflow, reloffset, nmax;
  long flowlimhi, flowlimlo, maxflow, minflow, tempflow;
  double dphi, dpsi;
  costT **loweredgecosts, **costs, *costsbelow;
  smoothcostT **loweredgesmoothcosts, **smoothcosts, *smoothcostsbelow;
  long nshortcycle;

  /* typecast generic pointers to costT pointers */
  loweredgecosts=(costT **)voidloweredgecosts;
  costs=(costT **)voidcosts;
  costsbelow=(costT *)voidcostsbelow;
  loweredgesmoothcosts=(smoothcostT **)voidloweredgecosts;
  smoothcosts=(smoothcostT **)voidcosts;
  smoothcostsbelow=(smoothcostT *)voidcostsbelow;

  /* see if tile is in bottom row */
  if(tilerow!=params->ntilerow-1){
  
    /* set up */
    nshortcycle=params->nshortcycle;
    flowlimhi=LARGESHORT;
    flowlimlo=-LARGESHORT;
    flowhistogram=(long *)CAlloc(flowlimhi-flowlimlo+1,sizeof(long));
    minflow=flowlimhi;
    maxflow=flowlimlo;

    /* loop over all arcs on the boundary */
    for(col=0;col<ncol;col++){
      dphi=(unwphase[nrow-1][col]-unwphasebelow[col])/TWOPI;
      tempflow=(short )LRound(dphi);
      loweredgeflows[0][col]=tempflow;
      if(tempflow<minflow){
	if(tempflow<flowlimlo){
	  fprintf(sp0,"Overflow in tile offset\nAbort\n");
	  exit(ABNORMAL_EXIT);
	}
	minflow=tempflow;
      }
      if(tempflow>maxflow){
	if(tempflow>flowlimhi){
	  fprintf(sp0,"Overflow in tile offset\nAbort\n");
	  exit(ABNORMAL_EXIT);
	}
	maxflow=tempflow;
      }
      flowhistogram[tempflow-flowlimlo]++;
      dpsi=dphi-floor(dphi);
      if(dpsi>0.5){
	dpsi-=1.0;
      }
      if(params->costmode==TOPO || params->costmode==DEFO){
	loweredgecosts[0][col].offset=nshortcycle*dpsi;
	loweredgecosts[0][col].sigsq=ceil((costs[nrow-2][col].sigsq
					   +costsbelow[col].sigsq)/2.0);
	if(costs[nrow-2][col].dzmax>cos