snaphu_tile.c

phase unwrapping algorithm for SAR interferometry

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

  snaphu tile-mode source file
  Written by Curtis W. Chen
  Copyright 2002 Board of Trustees, Leland Stanford Jr. University
  Please see the supporting documentation for terms of use.
  No warranty.

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

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <signal.h>
#include <limits.h>
#include <float.h>
#include <string.h>
#include <ctype.h>
#include <unistd.h>
#include <fcntl.h>
#include <sys/stat.h>
#include <sys/types.h>
#include <sys/wait.h>
#include <sys/time.h>
#include <sys/resource.h>

#include "snaphu.h"



/* function: SetupTile()
 * ---------------------
 * Sets up tile parameters and output file names for the current tile.
 */
void SetupTile(long nlines, long linelen, paramT *params, 
	       tileparamT *tileparams, outfileT *outfiles, 
	       outfileT *tileoutfiles, long tilerow, long tilecol){

  long ni, nj;
  char tempstring[MAXTMPSTRLEN], path[MAXSTRLEN], basename[MAXSTRLEN];
  char *tiledir;


  /* set parameters for current tile */
  ni=ceil((nlines+(params->ntilerow-1)*params->rowovrlp)
	  /(double )params->ntilerow);
  nj=ceil((linelen+(params->ntilecol-1)*params->colovrlp)
	  /(double )params->ntilecol);
  tileparams->firstrow=tilerow*(ni-params->rowovrlp);
  tileparams->firstcol=tilecol*(nj-params->colovrlp);
  if(tilerow==params->ntilerow-1){
    tileparams->nrow=nlines-(params->ntilerow-1)*(ni-params->rowovrlp);
  }else{
    tileparams->nrow=ni;
  }
  if(tilecol==params->ntilecol-1){
    tileparams->ncol=linelen-(params->ntilecol-1)*(nj-params->colovrlp);
  }else{
    tileparams->ncol=nj;
  }

  /* set output files */
  tiledir=params->tiledir;
  ParseFilename(outfiles->outfile,path,basename);
  sprintf(tempstring,"%s/%s%s_%ld_%ld.%ld",
	  tiledir,TMPTILEROOT,basename,tilerow,tilecol,tileparams->ncol);
  StrNCopy(tileoutfiles->outfile,tempstring,MAXSTRLEN);
  if(strlen(outfiles->initfile)){
    ParseFilename(outfiles->initfile,path,basename);
    sprintf(tempstring,"%s/%s%s_%ld_%ld.%ld",
	    tiledir,TMPTILEROOT,basename,tilerow,tilecol,tileparams->ncol);
    StrNCopy(tileoutfiles->initfile,tempstring,MAXSTRLEN);
  }else{
    StrNCopy(tileoutfiles->initfile,"",MAXSTRLEN);
  }
  if(strlen(outfiles->flowfile)){
    ParseFilename(outfiles->flowfile,path,basename);
    sprintf(tempstring,"%s/%s%s_%ld_%ld.%ld",
	    tiledir,TMPTILEROOT,basename,tilerow,tilecol,tileparams->ncol);
    StrNCopy(tileoutfiles->flowfile,tempstring,MAXSTRLEN);
  }else{
    StrNCopy(tileoutfiles->flowfile,"",MAXSTRLEN);
  }
  if(strlen(outfiles->eifile)){
    ParseFilename(outfiles->eifile,path,basename);
    sprintf(tempstring,"%s/%s%s_%ld_%ld.%ld",
	    tiledir,TMPTILEROOT,basename,tilerow,tilecol,tileparams->ncol);
    StrNCopy(tileoutfiles->eifile,tempstring,MAXSTRLEN);
  }else{
    StrNCopy(tileoutfiles->eifile,"",MAXSTRLEN);
  }
  if(strlen(outfiles->rowcostfile)){
    ParseFilename(outfiles->rowcostfile,path,basename);
    sprintf(tempstring,"%s/%s%s_%ld_%ld.%ld",
	    tiledir,TMPTILEROOT,basename,tilerow,tilecol,tileparams->ncol);
    StrNCopy(tileoutfiles->rowcostfile,tempstring,MAXSTRLEN);
  }else{
    StrNCopy(tileoutfiles->rowcostfile,"",MAXSTRLEN);
  }
  if(strlen(outfiles->colcostfile)){
    ParseFilename(outfiles->colcostfile,path,basename);
    sprintf(tempstring,"%s/%s%s_%ld_%ld.%ld",
	    tiledir,TMPTILEROOT,basename,tilerow,tilecol,tileparams->ncol);
    StrNCopy(tileoutfiles->colcostfile,tempstring,MAXSTRLEN);
  }else{
    StrNCopy(tileoutfiles->colcostfile,"",MAXSTRLEN);
  }
  if(strlen(outfiles->mstrowcostfile)){
    ParseFilename(outfiles->mstrowcostfile,path,basename);
    sprintf(tempstring,"%s/%s%s_%ld_%ld.%ld",
	    tiledir,TMPTILEROOT,basename,tilerow,tilecol,tileparams->ncol);
    StrNCopy(tileoutfiles->mstrowcostfile,tempstring,MAXSTRLEN);
  }else{
    StrNCopy(tileoutfiles->mstrowcostfile,"",MAXSTRLEN);
  }
  if(strlen(outfiles->mstcolcostfile)){
    ParseFilename(outfiles->mstcolcostfile,path,basename);
    sprintf(tempstring,"%s/%s%s_%ld_%ld.%ld",
	    tiledir,TMPTILEROOT,basename,tilerow,tilecol,tileparams->ncol);
    StrNCopy(tileoutfiles->mstcolcostfile,tempstring,MAXSTRLEN);
  }else{
    StrNCopy(tileoutfiles->mstcolcostfile,"",MAXSTRLEN);
  }
  if(strlen(outfiles->mstcostsfile)){
    ParseFilename(outfiles->mstcostsfile,path,basename);
    sprintf(tempstring,"%s/%s%s_%ld_%ld.%ld",
	    tiledir,TMPTILEROOT,basename,tilerow,tilecol,tileparams->ncol);
    StrNCopy(tileoutfiles->mstcostsfile,tempstring,MAXSTRLEN);
  }else{
    StrNCopy(tileoutfiles->mstcostsfile,"",MAXSTRLEN);
  }
  if(strlen(outfiles->corrdumpfile)){
    ParseFilename(outfiles->corrdumpfile,path,basename);
    sprintf(tempstring,"%s/%s%s_%ld_%ld.%ld",
	    tiledir,TMPTILEROOT,basename,tilerow,tilecol,tileparams->ncol);
    StrNCopy(tileoutfiles->corrdumpfile,tempstring,MAXSTRLEN);
  }else{
    StrNCopy(tileoutfiles->corrdumpfile,"",MAXSTRLEN);
  }
  if(strlen(outfiles->rawcorrdumpfile)){
    ParseFilename(outfiles->rawcorrdumpfile,path,basename);
    sprintf(tempstring,"%s/%s%s_%ld_%ld.%ld",
	    tiledir,TMPTILEROOT,basename,tilerow,tilecol,tileparams->ncol);
    StrNCopy(tileoutfiles->rawcorrdumpfile,tempstring,MAXSTRLEN);
  }else{
    StrNCopy(tileoutfiles->rawcorrdumpfile,"",MAXSTRLEN);
  }
  if(strlen(outfiles->costoutfile)){
    ParseFilename(outfiles->costoutfile,path,basename);
    sprintf(tempstring,"%s/%s%s_%ld_%ld.%ld",
	    tiledir,TMPTILEROOT,basename,tilerow,tilecol,tileparams->ncol);
    StrNCopy(tileoutfiles->costoutfile,tempstring,MAXSTRLEN);
  }else{
    sprintf(tempstring,"%s/%s%s%ld_%ld.%ld",
	    tiledir,TMPTILEROOT,TMPTILECOSTSUFFIX,tilerow,tilecol,
	    tileparams->ncol);
    StrNCopy(tileoutfiles->costoutfile,tempstring,MAXSTRLEN);
  }
  tileoutfiles->outfileformat=TMPTILEOUTFORMAT;

}


/* function: GrowRegions()
 * -----------------------
 * Grows contiguous regions demarcated by arcs whose residual costs are
 * less than some threshold.  Numbers the regions sequentially from 0.
 */
void GrowRegions(void **costs, short **flows, long nrow, long ncol, 
		 incrcostT **incrcosts, outfileT *outfiles, paramT *params){

  long i, row, col, maxcol;
  long arcrow, arccol, arcnum, fromdist, arcdist;
  long regioncounter, *regionsizes, regionsizeslen, *thisregionsize;
  long closestregiondist, closestregion, lastfromdist;
  long costthresh, minsize, maxcost;
  short **regions;
  nodeT **nodes;
  nodeT *source, *from, *to, *ground;
  char regionfile[MAXSTRLEN];
  bucketT bkts[1];


  /* error checking */
  fprintf(sp1,"Growing reliable regions\n");
  minsize=params->minregionsize;
  costthresh=params->tilecostthresh;
  if(minsize>nrow*ncol){
    fprintf(sp0,"Minimum region size cannot exceed tile size\nAbort\n");
    exit(ABNORMAL_EXIT);
  }

  /* loop over all arcs */
  for(arcrow=0;arcrow<2*nrow-1;arcrow++){
    if(arcrow<nrow-1){
      maxcol=ncol;
    }else{
      maxcol=ncol-1;
    }
    for(arccol=0;arccol<maxcol;arccol++){

      /* compute incremental costs of unit flows in either direction */
      ReCalcCost(costs,incrcosts,flows[arcrow][arccol],
		 arcrow,arccol,1,nrow,params);

      /* store lesser of incremental costs in first field */
      if(incrcosts[arcrow][arccol].negcost<incrcosts[arcrow][arccol].poscost){
	incrcosts[arcrow][arccol].poscost=incrcosts[arcrow][arccol].negcost;
      }

      /* subtract costthresh and take negative of costs, then clip to zero */
      incrcosts[arcrow][arccol].poscost
	=-(incrcosts[arcrow][arccol].poscost-costthresh);
      if(incrcosts[arcrow][arccol].poscost<0){
	incrcosts[arcrow][arccol].poscost=0;
      }	
    }
  }

  /* thicken the costs arrays; results stored in negcost field */
  maxcost=ThickenCosts(incrcosts,nrow,ncol);
 
  /* initialize nodes and buckets for region growing */
  ground=NULL;
  nodes=(nodeT **)Get2DMem(nrow,ncol,sizeof(nodeT *),sizeof(nodeT));
  InitNodeNums(nrow,ncol,nodes,ground);
  InitNodes(nrow,ncol,nodes,ground);
  bkts->size=maxcost+2;
  bkts->minind=0;
  bkts->maxind=bkts->size-1;
  bkts->curr=0;
  bkts->wrapped=FALSE;
  bkts->bucketbase=(nodeT **)MAlloc(bkts->size*sizeof(nodeT *));
  bkts->bucket=bkts->bucketbase;
  for(i=0;i<bkts->size;i++){
    bkts->bucket[i]=NULL;
  }

  /* initialize region variables */
  regioncounter=-1;
  regionsizeslen=INITARRSIZE;
  regionsizes=(long *)MAlloc(regionsizeslen*sizeof(long));
  for(row=0;row<nrow;row++){
    for(col=0;col<ncol;col++){
      nodes[row][col].incost=-1;
    }
  }

  /* loop over all nodes (pixels) to make sure each is in a group */
  for(row=0;row<nrow;row++){
    for(col=0;col<ncol;col++){

      /* see if node is not in a group */
      if(nodes[row][col].incost<0){

	/* clear the buckets */
	ClearBuckets(bkts);

	/* make node source and put it in the first bucket */
	source=&nodes[row][col];
	source->next=NULL;
	source->prev=NULL;
	source->group=INBUCKET;
	source->outcost=0;
	bkts->bucket[0]=source;
	bkts->curr=0;
	lastfromdist=0;

	/* increment the region counter */
	if(++regioncounter>=regionsizeslen){
	  regionsizeslen+=INITARRSIZE;
	  regionsizes=(long *)ReAlloc(regionsizes,
				       regionsizeslen*sizeof(long));
	}
	thisregionsize=&regionsizes[regioncounter];

	/* set up */
	(*thisregionsize)=0;
	closestregiondist=VERYFAR;

	/* loop to grow region */
	while(TRUE){

	  /* set from node to closest node in circular bucket structure */
	  from=ClosestNode(bkts);
	  
	  /* break if we can't grow any more and the region is big enough */
	  if(from==NULL){
	    if(*thisregionsize>=minsize){

	      /* no more nonregion nodes, and current region is big enough */
	      break;

	    }else{

	      /* no more nonregion nodes, but current region still too small */
	      /* merge with another region */
	      MergeRegions(nodes,source,regionsizes,closestregion,nrow,ncol);
	      regioncounter--;
	      break;

	    }
	  }else{
	    fromdist=from->outcost;
	    if(fromdist>lastfromdist){
	      if(regionsizes[regioncounter]>=minsize){

		/* region grown to all nodes within mincost, is big enough */
		break;

	      }
	      if(fromdist>closestregiondist){

		/* another region closer than new node, so merge regions */
		MergeRegions(nodes,source,regionsizes,closestregion,nrow,ncol);
		regioncounter--;
		break;
	      }
	    }
	  }

	  /* make from node a part of the current region */
	  from->incost=regioncounter;
	  (*thisregionsize)++;
	  lastfromdist=fromdist;

	  /* scan from's neighbors */
	  arcnum=0;
	  while((to=RegionsNeighborNode(from,&arcnum,nodes,
					&arcrow,&arccol,nrow,ncol))!=NULL){
	   
	    /* get cost of arc to the to node */
	    arcdist=incrcosts[arcrow][arccol].negcost;

	    /* see if to node is already in another region */
	    if(to->incost>=0){

	      /* keep track of which neighboring region is closest */
	      if(to->incost!=regioncounter && arcdist<closestregiondist){
		closestregiondist=arcdist;
		closestregion=to->incost;
	      }

	    }else{

	      /* to node is not in another region */
	      /* compare distance of new nodes to temp labels */
	      if(arcdist<(to->outcost)){

		/* if to node is already in a (circular) bucket, remove it */
		if(to->group==INBUCKET){
		  BucketRemove(to,to->outcost,bkts);
		}
                
		/* update to node */
		to->outcost=arcdist;
		to->pred=from;

		/* insert to node into appropriate (circular) bucket */
		BucketInsert(to,arcdist,bkts);
		if(arcdist<bkts->curr){
		  bkts->curr=arcdist;
		}
	      }
	    }
	  }
	}
      }
    }
  }
  fprintf(sp2,"Tile partitioned into %ld regions\n",regioncounter+1);

  /* write regions array */
  /* write as shorts if multiple tiles */
  if(params->ntilerow > 1 || params->ntilecol>1){
    regions=(short **)Get2DMem(nrow,ncol,sizeof(short *),sizeof(short));
    for(row=0;row<nrow;row++){
      for(col=0;col<ncol;col++){
	if(nodes[row][col].incost>LARGESHORT){
	  fprintf(sp0,
		  "Number of regions in tile exceeds max allowed\nAbort\n");
	  exit(ABNORMAL_EXIT);
	}
	regions[row][col]=nodes[row][col].incost;
      }
    }
    sprintf(regionfile,"%s%s",outfiles->outfile,REGIONSUFFIX);
    fprintf(sp2,"Writing region data to file %s\n",regionfile);
    Write2DArray((void **)regions,regionfile,nrow,ncol,sizeof(short));
  }

  /* free memory */
  Free2DArray((void **)nodes,nrow);
  Free2DArray((void **)regions,nrow);
  free(bkts->bucketbase);

}


/* function: GrowConnCompMask()
 * ----------------------------
 * Grows contiguous regions demarcated by arcs whose residual costs are
 * less than some threshold.  Numbers the regions sequentially from 1.
 * Writes out byte file of connected component mask, with 0 for any pixels
 * not assigned to a component.
 */
void GrowConnCompsMask(void **costs, short **flows, long nrow, long ncol, 
		       incrcostT **incrcosts, outfileT *outfiles, 
		       paramT *params){

  long i, row, col, maxcol;
  long arcrow, arccol, arcnum;
  long regioncounter, *regionsizes, regionsizeslen, *thisregionsize;
  long *sortedregionsizes;
  long costthresh, minsize, maxncomps, ntied, newnum;
  nodeT **nodes;