snaphu_tile.c

phase unwrapping algorithm for SAR interferometry

		      nexttilelinelen,nexttilenlines,nexttileparams,
		      costtypesize);
  }else{
    fprintf(sp0,"Tile reassembly not enabled in Lp mode\nAbort\n");
    exit(ABNORMAL_EXIT);
  }

  /* flip sign of wrapped phase if flip flag is set */
  FlipPhaseArraySign(*nextunwphaseptr,params,
		     nexttileparams->nrow,nexttileparams->ncol);

  /* set outputs */
  (*nextnrowptr)=nexttileparams->nrow;
  (*nextncolptr)=nexttileparams->ncol;

}

/* function: SetTileReadParams()
 * -----------------------------
 * Set parameters for reading the nonoverlapping piece of each tile.  
 * ni and nj are the numbers of rows and columns in this particular tile.
 * The meanings of these variables are different for the last row 
 * and column.
 */
void SetTileReadParams(tileparamT *tileparams, long nexttilenlines, 
		       long nexttilelinelen, long tilerow, long tilecol, 
		       long nlines, long linelen, paramT *params){

  long rowovrlp, colovrlp;

  /* set temporary variables */
  rowovrlp=params->rowovrlp;
  colovrlp=params->colovrlp;

  /* row parameters */
  if(tilerow==0){
    tileparams->firstrow=0;
  }else{
    tileparams->firstrow=ceil(rowovrlp/2.0);
  }
  if(tilerow!=params->ntilerow-1){
    tileparams->nrow=nexttilenlines-floor(rowovrlp/2.0)-tileparams->firstrow;
  }else{
    tileparams->nrow=nexttilenlines-tileparams->firstrow;
  }

  /* column parameters */
  if(tilecol==0){
    tileparams->firstcol=0;
  }else{
    tileparams->firstcol=ceil(colovrlp/2.0);
  }
  if(tilecol!=params->ntilecol-1){
    tileparams->ncol=nexttilelinelen-floor(colovrlp/2.0)-tileparams->firstcol;
  }else{
    tileparams->ncol=nexttilelinelen-tileparams->firstcol;
  }
}


/* function: ReadEdgesAboveAndBelow()
 * ----------------------------------
 */
void ReadEdgesAboveAndBelow(long tilerow, long tilecol, long nlines, 
			    long linelen, paramT *params, outfileT *outfiles, 
			    short *regionsabove, short *regionsbelow,
			    float *unwphaseabove, float *unwphasebelow,
			    void *costsabove, void *costsbelow){

  long ni, nj, readtilelinelen, readtilenlines, costtypesize;
  long ntilerow, ntilecol, rowovrlp, colovrlp;
  tileparamT tileparams[1];
  outfileT outfilesabove[1], outfilesbelow[1];
  float **unwphaseaboveptr, **unwphasebelowptr;
  void **costsaboveptr, **costsbelowptr;
  short **regionsaboveptr, **regionsbelowptr;
  char tempstring[MAXTMPSTRLEN], readregionfile[MAXSTRLEN];
  char path[MAXSTRLEN], basename[MAXSTRLEN];

  /* set temporary variables */
  ntilerow=params->ntilerow;
  ntilecol=params->ntilecol;
  rowovrlp=params->rowovrlp;
  colovrlp=params->colovrlp;
  ni=ceil((nlines+(ntilerow-1)*rowovrlp)/(double )ntilerow);
  nj=ceil((linelen+(ntilecol-1)*colovrlp)/(double )ntilecol);

  /* size of the data type for holding cost data depends on cost mode */
  if(params->costmode==TOPO){
    costtypesize=sizeof(costT);
  }else if(params->costmode==DEFO){
    costtypesize=sizeof(costT);
  }else if(params->costmode==SMOOTH){
    costtypesize=sizeof(smoothcostT);
  }

  /* set names of files with SetupTile() */
  /* tile parameters set by SetupTile() will be overwritten below */
  if(tilerow!=0){
    SetupTile(nlines,linelen,params,tileparams,outfiles,outfilesabove,
	      tilerow-1,tilecol);
  }
  if(tilerow!=ntilerow-1){
    SetupTile(nlines,linelen,params,tileparams,outfiles,outfilesbelow,
	      tilerow+1,tilecol);
  }

  /* temporary pointers, so we can use Read2DArray() with 1D output array */
  unwphaseaboveptr=&unwphaseabove;
  unwphasebelowptr=&unwphasebelow;
  costsaboveptr=&costsabove;
  costsbelowptr=&costsbelow;
  regionsaboveptr=&regionsabove;
  regionsbelowptr=&regionsbelow;
  
  /* set some reading parameters */
  if(tilecol==0){
    tileparams->firstcol=0;
  }else{
    tileparams->firstcol=ceil(colovrlp/2.0);
  }
  if(tilecol!=params->ntilecol-1){
    readtilelinelen=nj;
    tileparams->ncol=readtilelinelen-floor(colovrlp/2.0)-tileparams->firstcol;
  }else{
    readtilelinelen=linelen-(ntilecol-1)*(nj-colovrlp);
    tileparams->ncol=readtilelinelen-tileparams->firstcol;
  }
  tileparams->nrow=1;

  /* read last line of tile above */
  readtilenlines=ni;
  if(tilerow!=0){
    tileparams->firstrow=readtilenlines-floor(rowovrlp/2.0)-1;

    /* read region data */
    ParseFilename(outfiles->outfile,path,basename);
    sprintf(tempstring,"%s/%s%s_%ld_%ld.%ld%s",
	    params->tiledir,TMPTILEROOT,basename,tilerow-1,tilecol,
	    readtilelinelen,REGIONSUFFIX);
    StrNCopy(readregionfile,tempstring,MAXSTRLEN);
    Read2DArray((void ***)&regionsaboveptr,readregionfile,
		readtilelinelen,readtilenlines,
		tileparams,sizeof(short *),sizeof(short));

    /* read unwrapped phase data */
    if(TMPTILEOUTFORMAT==ALT_LINE_DATA){
      ReadAltLineFilePhase(&unwphaseaboveptr,outfilesabove->outfile,
			   readtilelinelen,readtilenlines,tileparams);
    }else if(TMPTILEOUTFORMAT==FLOAT_DATA){
      Read2DArray((void ***)&unwphaseaboveptr,outfilesabove->outfile,
		  readtilelinelen,readtilenlines,tileparams,
		  sizeof(float *),sizeof(float));
    }

    /* flip sign of wrapped phase if flip flag is set */
    FlipPhaseArraySign(unwphaseaboveptr,params,
		       tileparams->nrow,tileparams->ncol);

    /* read costs data */
    tileparams->firstrow--;
    Read2DRowColFileRows((void ***)&costsaboveptr,outfilesabove->costoutfile,
			 readtilelinelen,readtilenlines,tileparams,
			 costtypesize);

    /* remove temporary tile cost file unless told to save it */
    if(params->rmtmptile && !strlen(outfiles->costoutfile)){
      unlink(outfilesabove->costoutfile);
    }
  }

  /* read first line of tile below */
  if(tilerow!=ntilerow-1){
    if(tilerow==params->ntilerow-2){
      readtilenlines=nlines-(ntilerow-1)*(ni-rowovrlp);
    }
    tileparams->firstrow=ceil(rowovrlp/2.0);

    /* read region data */
    ParseFilename(outfiles->outfile,path,basename);
    sprintf(tempstring,"%s/%s%s_%ld_%ld.%ld%s",
	    params->tiledir,TMPTILEROOT,basename,tilerow+1,tilecol,
	    readtilelinelen,REGIONSUFFIX);
    StrNCopy(readregionfile,tempstring,MAXSTRLEN);
    Read2DArray((void ***)&regionsbelowptr,readregionfile,
		readtilelinelen,readtilenlines,
		tileparams,sizeof(short *),sizeof(short));

    /* read unwrapped phase data */
    if(TMPTILEOUTFORMAT==ALT_LINE_DATA){
      ReadAltLineFilePhase(&unwphasebelowptr,outfilesbelow->outfile,
			   readtilelinelen,readtilenlines,tileparams);
    }else if(TMPTILEOUTFORMAT==FLOAT_DATA){
      Read2DArray((void ***)&unwphasebelowptr,outfilesbelow->outfile,
		  readtilelinelen,readtilenlines,tileparams,
		  sizeof(float *),sizeof(float));
    }

    /* flip the sign of the wrapped phase if flip flag is set */
    FlipPhaseArraySign(unwphasebelowptr,params,
		       tileparams->nrow,tileparams->ncol);

    /* read costs data */
    Read2DRowColFileRows((void ***)&costsbelowptr,outfilesbelow->costoutfile,
			 readtilelinelen,readtilenlines,tileparams,
			 costtypesize);

  }else{

    /* remove temporoary tile cost file for last row unless told to save it */
    if(params->rmtmptile && !strlen(outfiles->costoutfile)){
      SetupTile(nlines,linelen,params,tileparams,outfiles,outfilesbelow,
		tilerow,tilecol);
      unlink(outfilesbelow->costoutfile);
    }
  }
}


/* function: TraceRegions()
 * ------------------------
 * Trace edges of region data to form nodes and arcs of secondary
 * (ie, region-level) network problem.  Primary nodes and arcs are
 * those of the original, pixel-level network problem.  Flows along
 * edges are computed knowing the unwrapped phase values of edges
 * of adjacent tiles.  Costs along edges are approximated in that they
 * are calculated from combining adjacent cost parameters, not from 
 * using the exact method in BuildCostArrays().
 */
void TraceRegions(short **regions, short **nextregions, short **lastregions, 
		  short *regionsabove, short *regionsbelow, float **unwphase, 
		  float **nextunwphase, float **lastunwphase, 
		  float *unwphaseabove, float *unwphasebelow, void **costs, 
		  void **nextcosts, void **lastcosts, void *costsabove, 
		  void *costsbelow, long prevnrow, long prevncol, long tilerow,
		  long tilecol, long nrow, long ncol, nodeT **scndrynodes,
		  nodesuppT **nodesupp, scndryarcT **scndryarcs, 
		  long ***scndrycosts, short *nscndrynodes, 
		  short *nscndryarcs, long *totarclens, short **bulkoffsets, 
		  paramT *params){

  long i, j, row, col, nnrow, nncol, tilenum, costtypesize;
  long nnewnodes, nnewarcs, npathsout, flowmax, totarclen;
  long nupdatednontilenodes, updatednontilenodesize, ntilecol;
  short **flows;
  short **rightedgeflows, **loweredgeflows, **leftedgeflows, **upperedgeflows;
  short *inontilenodeoutarc;
  void **rightedgecosts, **loweredgecosts, **leftedgecosts, **upperedgecosts;
  nodeT **primarynodes, **updatednontilenodes;
  nodeT *from, *to, *nextnode, *tempnode;
  nodesuppT *fromsupp, *tosupp;


  /* initialize */
  ntilecol=params->ntilecol;
  nnrow=nrow+1;
  nncol=ncol+1;
  primarynodes=(nodeT **)Get2DMem(nnrow,nncol,sizeof(nodeT *),sizeof(nodeT));
  for(row=0;row<nnrow;row++){
    for(col=0;col<nncol;col++){
      primarynodes[row][col].row=row;
      primarynodes[row][col].col=col;
      primarynodes[row][col].group=NOTINBUCKET;
      primarynodes[row][col].pred=NULL;
      primarynodes[row][col].next=NULL;
    }
  }
  nextnode=&primarynodes[0][0];
  tilenum=tilerow*ntilecol+tilecol;
  scndrynodes[tilenum]=NULL;
  nodesupp[tilenum]=NULL;
  scndryarcs[tilenum]=NULL;
  scndrycosts[tilenum]=NULL;
  nnewnodes=0;
  nnewarcs=0;
  totarclen=0;
  flowmax=params->scndryarcflowmax;
  updatednontilenodesize=INITARRSIZE;
  nupdatednontilenodes=0;

  /* size of the data type for holding cost data depends on cost mode */
  if(params->costmode==TOPO){
    costtypesize=sizeof(costT);
  }else if(params->costmode==DEFO){
    costtypesize=sizeof(costT);
  }else if(params->costmode==SMOOTH){
    costtypesize=sizeof(smoothcostT);
  }

  /* get memory */
  updatednontilenodes=(nodeT **)MAlloc(updatednontilenodesize*sizeof(nodeT *));
  inontilenodeoutarc=(short *)MAlloc(updatednontilenodesize*sizeof(short));
  flows=(short **)Get2DRowColMem(nrow+1,ncol+1,sizeof(short *),sizeof(short));
  rightedgeflows=(short **)Get2DMem(nrow,1,sizeof(short *),sizeof(short));
  leftedgeflows=(short **)Get2DMem(nrow,1,sizeof(short *),sizeof(short));
  upperedgeflows=(short **)Get2DMem(1,ncol,sizeof(short *),sizeof(short));
  loweredgeflows=(short **)Get2DMem(1,ncol,sizeof(short *),sizeof(short));
  rightedgecosts=(void **)Get2DMem(nrow,1,sizeof(void *),costtypesize);
  leftedgecosts=(void **)Get2DMem(nrow,1,sizeof(void *),costtypesize);
  upperedgecosts=(void **)Get2DMem(1,ncol,sizeof(void *),costtypesize);
  loweredgecosts=(void **)Get2DMem(1,ncol,sizeof(void *),costtypesize);

  /* parse flows for this tile */
  CalcFlow(unwphase,&flows,nrow,ncol);

  /* set up cost and flow arrays for boundaries */
  SetUpperEdge(ncol,tilerow,tilecol,costs,costsabove,unwphase,unwphaseabove,
	       upperedgecosts,upperedgeflows,params, bulkoffsets);
  SetLowerEdge(nrow,ncol,tilerow,tilecol,costs,costsbelow,unwphase,
	       unwphasebelow,loweredgecosts,loweredgeflows,
	       params,bulkoffsets);
  SetLeftEdge(nrow,prevncol,tilerow,tilecol,costs,lastcosts,unwphase,
	      lastunwphase,leftedgecosts,leftedgeflows,params, bulkoffsets);
  SetRightEdge(nrow,ncol,tilerow,tilecol,costs,nextcosts,unwphase, 
	       nextunwphase,rightedgecosts,rightedgeflows,
	       params,bulkoffsets);

  /* trace edges between regions */
  while(nextnode!=NULL){
 
    /* get next primary node from stack */
    from=nextnode;
    nextnode=nextnode->next;
    from->group=NOTINBUCKET;

    /* find number of paths out of from node */
    npathsout=FindNumPathsOut(from,params,tilerow,tilecol,nnrow,nncol,regions,
			      nextregions,lastregions,regionsabove,
			      regionsbelow,prevncol);

    /* secondary node exists if region edges fork */
    if(npathsout>2){

      /* mark primary node to indicate that secondary node exists for it */
      from->group=ONTREE;
      
      /* create secondary node if not already created in another tile */
      if((from->row!=0 || tilerow==0) && (from->col!=0 || tilecol==0)){

	/* create the secondary node */
	nnewnodes++;
	scndrynodes[tilenum]=(nodeT *)ReAlloc(scndrynodes[tilenum],
					      nnewnodes*sizeof(nodeT));
	nodesupp[tilenum]=(nodesuppT *)ReAlloc(nodesupp[tilenum],
					       nnewnodes*sizeof(nodesuppT));
	scndrynodes[tilenum][nnewnodes-1].row=tilenum;
	scndrynodes[tilenum][nnewnodes-1].col=nnewnodes-1;
	nodesupp[tilenum][nnewnodes-1].row=from->row;
	nodesupp[tilenum][nnewnodes-1].col=from->col;
	nodesupp[tilenum][nnewnodes-1].noutarcs=0;
	nodesupp[tilenum][nnewnodes-1].neighbornodes=NULL;
	nodesupp[tilenum][nnewnodes-1].outarcs=NULL;
      }      

      /* create the secondary arc to this node if it doesn't already exist */
      if(from->pred!=NULL
	 && ((from->row==from->pred->row && (from->row!=0 || tilerow==0))
	     || (from->col==from->pred->col && (from->col!=0 || tilecol==0)))){

	TraceSecondaryArc(from,scndrynodes,nodesupp,scndryarcs,scndrycosts,
			  &nnewnodes,&nnewarcs,tilerow,tilecol,flowmax,
			  nrow,ncol,prevnrow,prevncol,params,costs,
			  rightedgecosts,loweredgecosts,leftedgecosts,
			  upperedgecosts,flows,rightedgeflows,loweredgeflows, 
			  leftedgeflows,upperedgeflows,&updatednontilenodes,
			  &nupdatednontilenodes,&updatednontilenodesize,
			  &inontilenodeoutarc,&totarclen);
      }
    }

    /* scan neighboring primary nodes and place path candidates into stack */
    RegionTraceCheckNeighbors(from,&nextnode,primarynodes,regions,
			      nextregions,lastregions,regionsabove,
			      regionsbelow,tilerow,tilecol,nnrow,nncol,
			      scndrynodes,nodesupp,scndryarcs,&nnewnodes,
			      &nnewarcs,flowmax,nrow,ncol,prevnrow,prevncol,
			      params,costs,rightedgecosts,loweredgecosts,
			      leftedgecosts,upperedgecosts,flows,
			      rightedgeflows,loweredgeflows,leftedgeflows,
			      upperedgeflows,scndrycosts,&updatednontilenodes,
			      &nupdatednontilenodes,&updatednontilenodesize,
			      &inontilenodeoutarc,&totarclen);
  }


  /* reset temporary secondary node and arc pointers in data structures */
  /* secondary node row, col stored level, incost of primary node pointed to */

  /* update nodes in this tile */
  for(i=0;i<nnewnodes;i++){
    for(j=0;j<nodesupp[tilenum][i].noutarcs;j++){
      tempnode=nodesupp[tilenum][i].neighbornodes[j];
      nodesupp[tilenum][i].neighbornodes[j]
	=&scndrynodes[tempnode->level][tempnode->incost];
    }
  }

  /* update nodes not in this tile that were affected (that have new arcs) */
  for(i=0;i<nupdatednontilenodes;i++){
    row=updatednontilenodes[i]->row;
    col=updatednontilenodes[i]->col;
    j=inontilenodeoutarc[i];
    tempnode=nodesupp[row][col].neighbornodes[j];
    nodesupp[row][col].neighbornodes[j]
      =&scndrynodes[tempnode->level][tempnode->incost];
  }

  /* update secondary arcs */
  for(i=0;i<nnewarcs;i++){

    /* update node pointers in secondary arc structure */
    tempnode=scndryarcs[tilenum][i].from;
    scndryarcs[tilenum][i].from
      =&scndrynodes[tempnode->level][tempnode->incost];
    from=scndryarcs[tilenum][i].from;
    tempnode=scndryarcs[tilenum][i].to;