snaphu_tile.c

phase unwrapping algorithm for SAR interferometry

  /* free memory */
  free(nextnodelist);

}


/* function: RenumberRegion()
 * --------------------------
 * 
 */
void RenumberRegion(nodeT **nodes, nodeT *source, long newnum, 
		    long nrow, long ncol){

  long nextnodelistlen, nextnodelistnext, arcnum, arcrow, arccol, regionnum;
  nodeT *from, *to, **nextnodelist;

  
  /* initialize */
  nextnodelistlen=INITARRSIZE;
  nextnodelist=(nodeT **)MAlloc(nextnodelistlen*sizeof(nodeT **));
  nextnodelist[0]=source;
  nextnodelistnext=1;
  regionnum=source->incost;


  /* find all nodes in current region and switch their regions */
  while(nextnodelistnext){
    from=nextnodelist[--nextnodelistnext];
    from->incost=newnum;
    arcnum=0;
    while((to=RegionsNeighborNode(from,&arcnum,nodes,
				  &arcrow,&arccol,nrow,ncol))!=NULL){
      if(to->incost==regionnum){
	if(nextnodelistnext>=nextnodelistlen){
	  nextnodelistlen+=INITARRSIZE;
	  nextnodelist=(nodeT **)ReAlloc(nextnodelist,
					 nextnodelistlen*sizeof(nodeT *));
	}
	nextnodelist[nextnodelistnext++]=to;
      }
    }
  }

  /* free memory */
  free(nextnodelist);

}


/* function: AssembleTiles()
 * -------------------------
 */
void AssembleTiles(outfileT *outfiles, paramT *params, 
		   long nlines, long linelen){

  long tilerow, tilecol, ntilerow, ntilecol, ntiles, rowovrlp, colovrlp;
  long i, j, k, ni, nj, dummylong, costtypesize;
  long nrow, ncol, prevnrow, prevncol, nextnrow, nextncol;
  long n, ncycle, nflowdone, nflow, candidatelistsize, candidatebagsize;
  long nnodes, maxnflowcycles, arclen, narcs, sourcetilenum, flowmax;
  long *totarclens;
  long ***scndrycosts;
  double avgarclen;
  float **unwphase, **nextunwphase, **lastunwphase, **tempunwphase;
  float *unwphaseabove, *unwphasebelow;
  void **costs, **nextcosts, **lastcosts, **tempcosts;
  void *costsabove, *costsbelow;
  short **scndryflows, **bulkoffsets, **regions, **nextregions, **lastregions;
  short **tempregions, *regionsbelow, *regionsabove;
  short *nscndrynodes, *nscndryarcs;
  incrcostT **incrcosts;
  totalcostT totalcost, oldtotalcost;
  nodeT *source;
  nodeT **scndrynodes, ***scndryapexes;
  signed char **iscandidate;
  signed char notfirstloop;
  candidateT *candidatebag, *candidatelist;
  nodesuppT **nodesupp;
  scndryarcT **scndryarcs;
  bucketT *bkts;
  char filename[MAXSTRLEN];


  /* set up */
  fprintf(sp1,"Assembling tiles\n"); 
  ntilerow=params->ntilerow;
  ntilecol=params->ntilecol;
  ntiles=ntilerow*ntilecol;
  rowovrlp=params->rowovrlp;
  colovrlp=params->colovrlp;
  ni=ceil((nlines+(ntilerow-1)*rowovrlp)/(double )ntilerow);
  nj=ceil((linelen+(ntilecol-1)*colovrlp)/(double )ntilecol);
  nrow=0;
  ncol=0;
  flowmax=params->scndryarcflowmax;
  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 */
  regions=(short **)Get2DMem(ni,nj,sizeof(short *),sizeof(short));
  nextregions=(short **)Get2DMem(ni,nj,sizeof(short *),sizeof(short));
  lastregions=(short **)Get2DMem(ni,nj,sizeof(short *),sizeof(short));
  regionsbelow=(short *)MAlloc(nj*sizeof(short));
  regionsabove=(short *)MAlloc(nj*sizeof(short));
  unwphase=(float **)Get2DMem(ni,nj,sizeof(float *),sizeof(float));
  nextunwphase=(float **)Get2DMem(ni,nj,sizeof(float *),sizeof(float));
  lastunwphase=(float **)Get2DMem(ni,nj,sizeof(float *),sizeof(float));
  unwphaseabove=(float *)MAlloc(nj*sizeof(float));
  unwphasebelow=(float *)MAlloc(nj*sizeof(float));
  scndrynodes=(nodeT **)MAlloc(ntiles*sizeof(nodeT *));
  nodesupp=(nodesuppT **)MAlloc(ntiles*sizeof(nodesuppT *));
  scndryarcs=(scndryarcT **)MAlloc(ntiles*sizeof(scndryarcT *));
  scndrycosts=(long ***)MAlloc(ntiles*sizeof(long **));
  nscndrynodes=(short *)MAlloc(ntiles*sizeof(short));
  nscndryarcs=(short *)MAlloc(ntiles*sizeof(short));
  totarclens=(long *)MAlloc(ntiles*sizeof(long));
  bulkoffsets=(short **)Get2DMem(ntilerow,ntilecol,sizeof(short *),
				 sizeof(short));
  costs=(void **)Get2DRowColMem(ni+2,nj+2,sizeof(void *),costtypesize);
  nextcosts=(void **)Get2DRowColMem(ni+2,nj+2,sizeof(void *),costtypesize);
  lastcosts=(void **)Get2DRowColMem(ni+2,nj+2,sizeof(void *),costtypesize);
  costsabove=(void *)MAlloc(nj*costtypesize);
  costsbelow=(void *)MAlloc(nj*costtypesize);


  /* trace regions and parse secondary nodes and arcs for each tile */
  bulkoffsets[0][0]=0;
  for(tilerow=0;tilerow<ntilerow;tilerow++){
    for(tilecol=0;tilecol<ntilecol;tilecol++){

      /* read region, unwrapped phase, and flow data */
      if(tilecol==0){
	ReadNextRegion(tilerow,0,nlines,linelen,outfiles,params,
		       &nextregions,&nextunwphase,&nextcosts,
		       &nextnrow,&nextncol);
	prevnrow=nrow;
	nrow=nextnrow;
      }
      prevncol=ncol;
      ncol=nextncol;
      tempregions=lastregions;
      lastregions=regions;
      regions=nextregions;
      nextregions=tempregions;
      tempunwphase=lastunwphase;
      lastunwphase=unwphase;
      unwphase=nextunwphase;
      nextunwphase=tempunwphase;
      tempcosts=lastcosts;
      lastcosts=costs;
      costs=nextcosts;
      nextcosts=tempcosts;
      if(tilecol!=ntilecol-1){
	ReadNextRegion(tilerow,tilecol+1,nlines,linelen,outfiles,params,
		       &nextregions,&nextunwphase,&nextcosts,
		       &nextnrow,&nextncol);
      }
      ReadEdgesAboveAndBelow(tilerow,tilecol,nlines,linelen,params,
			     outfiles,regionsabove,regionsbelow,
			     unwphaseabove,unwphasebelow,
			     costsabove,costsbelow);

      /* trace region edges to form nodes and arcs */
      TraceRegions(regions,nextregions,lastregions,regionsabove,regionsbelow,
		   unwphase,nextunwphase,lastunwphase,unwphaseabove,
		   unwphasebelow,costs,nextcosts,lastcosts,costsabove,
		   costsbelow,prevnrow,prevncol,tilerow,tilecol,
		   nrow,ncol,scndrynodes,nodesupp,scndryarcs,
		   scndrycosts,nscndrynodes,nscndryarcs,totarclens,
		   bulkoffsets,params);

    }
  }

  /* scale costs based on average number of primary arcs per secondary arc */
  arclen=0;
  narcs=0;
  for(i=0;i<ntiles;i++){
    arclen+=totarclens[i];
    narcs+=nscndryarcs[i];
  }
  avgarclen=arclen/narcs;

  /* may need to adjust scaling so fewer costs clipped */
  for(i=0;i<ntiles;i++){
    for(j=0;j<nscndryarcs[i];j++){
      if(scndrycosts[i][j][2*flowmax+1]!=ZEROCOSTARC){
	for(k=1;k<=2*flowmax;k++){
	  scndrycosts[i][j][k]=(long )ceil(scndrycosts[i][j][k]/avgarclen);
	}
	scndrycosts[i][j][2*flowmax+1]=LRound(scndrycosts[i][j][2*flowmax+1]
					      /avgarclen);
	if(scndrycosts[i][j][2*flowmax+1]<1){
	  scndrycosts[i][j][2*flowmax+1]=1;
	}
      }
    }
  }

  /* free some intermediate memory */
  Free2DArray((void **)regions,ni);
  Free2DArray((void **)nextregions,ni);
  Free2DArray((void **)lastregions,ni);
  Free2DArray((void **)unwphase,ni);
  Free2DArray((void **)nextunwphase,ni);
  Free2DArray((void **)lastunwphase,ni);
  Free2DArray((void **)costs,ni);
  Free2DArray((void **)nextcosts,ni);
  Free2DArray((void **)lastcosts,ni);
  free(costsabove);
  free(costsbelow);
  free(unwphaseabove);
  free(unwphasebelow);
  free(regionsabove);
  free(regionsbelow);


  /* get memory for nongrid arrays of secondary network problem */
  scndryflows=(short **)MAlloc(ntiles*sizeof(short *));
  iscandidate=(signed char **)MAlloc(ntiles*sizeof(signed char*));
  scndryapexes=(nodeT ***)MAlloc(ntiles*sizeof(nodeT **));
  incrcosts=(incrcostT **)MAlloc(ntiles*sizeof(incrcostT *));
  nnodes=0;
  for(i=0;i<ntiles;i++){
    scndryflows[i]=(short *)CAlloc(nscndryarcs[i],sizeof(short));
    iscandidate[i]=(signed char *)MAlloc(nscndryarcs[i]*sizeof(signed char));
    scndryapexes[i]=(nodeT **)MAlloc(nscndryarcs[i]*sizeof(nodeT *));
    incrcosts[i]=(incrcostT *)MAlloc(nscndryarcs[i]*sizeof(incrcostT));
    nnodes+=nscndrynodes[i];
  }

  /* set up network for secondary solver */
  InitNetwork(scndryflows,&dummylong,&ncycle,&nflowdone,&dummylong,&nflow,
	      &candidatebagsize,&candidatebag,&candidatelistsize,
	      &candidatelist,NULL,NULL,&bkts,&dummylong,NULL,NULL,NULL,
	      NULL,NULL,NULL,NULL,ntiles,0,&notfirstloop,&totalcost,params);


  /* set pointers to functions for nongrid secondary network */
  CalcCost=CalcCostNonGrid;
  EvalCost=EvalCostNonGrid;
  NeighborNode=NeighborNodeNonGrid;
  GetArc=GetArcNonGrid;


  /* solve the secondary network problem */
  /* main loop: loop over flow increments and sources */
  fprintf(sp1,"Running optimizer for secondary network\n");
  maxnflowcycles=LRound(nnodes*params->maxcyclefraction);
  while(TRUE){ 
 
    fprintf(sp1,"Flow increment: %ld  (Total improvements: %ld)\n",
            nflow,ncycle);

    /* set up the incremental (residual) cost arrays */
    SetupIncrFlowCosts((void **)scndrycosts,incrcosts,scndryflows,nflow,ntiles,
		       ntiles,nscndryarcs,params); 

    /* set the tree root (equivalent to source of shortest path problem) */
    sourcetilenum=(long )ntilecol*floor(ntilerow/2.0)+floor(ntilecol/2.0);
    source=&scndrynodes[sourcetilenum][0];

    /* run the solver, and increment nflowdone if no cycles are found */
    n=TreeSolve(scndrynodes,nodesupp,NULL,source,&candidatelist,&candidatebag,
                &candidatelistsize,&candidatebagsize,bkts,scndryflows,
		(void **)scndrycosts,incrcosts,scndryapexes,iscandidate,0,
		nflow,NULL,NULL,NULL,ntiles,nscndrynodes,ntiles,nscndryarcs,
		ntiles,0,NULL,params);
    
    /* evaluate and save the total cost (skip if first loop through nflow) */
    if(notfirstloop){
      oldtotalcost=totalcost;
      totalcost=EvaluateTotalCost((void **)scndrycosts,scndryflows,ntiles,0,
				  nscndryarcs,params);
      if(totalcost>oldtotalcost || (n>0 && totalcost==oldtotalcost)){
        fprintf(sp0,"Unexpected increase in total cost.  Breaking loop\n");
        break;
      }
    }

    /* consider this flow increment done if not too many neg cycles found */
    ncycle+=n;
    if(n<=maxnflowcycles){
      nflowdone++;
    }else{
      nflowdone=1;
    }

    /* break if we're done with all flow increments or problem is convex */
    if(nflowdone>=params->maxflow){
      break;
    }

    /* update flow increment */
    nflow++;
    if(nflow>params->maxflow){
      nflow=1;
      notfirstloop=TRUE;
    }

  } /* end loop until no more neg cycles */

  /* free some memory */
  for(i=0;i<ntiles;i++){
    for(j=0;j<nscndryarcs[i];j++){
      free(scndrycosts[i][j]);
    }
  }
  Free2DArray((void **)scndrycosts,ntiles);
  Free2DArray((void **)scndryapexes,ntiles);
  Free2DArray((void **)iscandidate,ntiles);
  Free2DArray((void **)incrcosts,ntiles);
  free(candidatebag);
  free(candidatelist);  
  free(bkts->bucketbase);

  /* integrate phase from secondary network problem */
  IntegrateSecondaryFlows(linelen,nlines,scndrynodes,nodesupp,scndryarcs,
			  nscndryarcs,scndryflows,bulkoffsets,outfiles,params);

  /* free remaining memory */
  for(i=0;i<ntiles;i++){
    for(j=0;j<nscndrynodes[i];j++){
      free(nodesupp[i][j].neighbornodes);
      free(nodesupp[i][j].outarcs);
    }
  }
  Free2DArray((void **)nodesupp,ntiles);
  Free2DArray((void **)scndrynodes,ntiles);
  Free2DArray((void **)scndryarcs,ntiles);
  Free2DArray((void **)scndryflows,ntiles);
  free(nscndrynodes);
  free(nscndryarcs);
  Free2DArray((void **)bulkoffsets,ntilerow);

  /* remove temporary tile log files and tile directory */
  if(params->rmtmptile){
    for(tilerow=0;tilerow<ntilerow;tilerow++){
      for(tilecol=0;tilecol<ntilecol;tilecol++){
	sprintf(filename,"%s/%s%ld_%ld",
		params->tiledir,LOGFILEROOT,tilerow,tilecol);
	unlink(filename);
      }
    }
    rmdir(params->tiledir);
  }

}


/* function: ReadNextRegion()
 * --------------------------
 */
void ReadNextRegion(long tilerow, long tilecol, long nlines, long linelen,
		    outfileT *outfiles, paramT *params, 
		    short ***nextregionsptr, float ***nextunwphaseptr,
		    void ***nextcostsptr, 
		    long *nextnrowptr, long *nextncolptr){

  long nexttilelinelen, nexttilenlines, costtypesize;
  tileparamT nexttileparams[1];
  outfileT nexttileoutfiles[1];
  char nextfile[MAXSTRLEN], tempstring[MAXTMPSTRLEN];
  char path[MAXSTRLEN], basename[MAXSTRLEN];
  
  /* 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);
  }

  /* use SetupTile() to set filenames only; tile params overwritten below */
  SetupTile(nlines,linelen,params,nexttileparams,outfiles,nexttileoutfiles,
	    tilerow,tilecol);
  nexttilenlines=nexttileparams->nrow;
  nexttilelinelen=nexttileparams->ncol;

  /* set tile parameters, overwriting values set by SetupTile() above */
  SetTileReadParams(nexttileparams,nexttilenlines,nexttilelinelen,
		    tilerow,tilecol,nlines,linelen,params);

  /* read region data */
  ParseFilename(outfiles->outfile,path,basename);
  sprintf(tempstring,"%s/%s%s_%ld_%ld.%ld%s",
	  params->tiledir,TMPTILEROOT,basename,tilerow,tilecol,
	  nexttilelinelen,REGIONSUFFIX);
  StrNCopy(nextfile,tempstring,MAXSTRLEN);
  Read2DArray((void ***)nextregionsptr,nextfile,
	      nexttilelinelen,nexttilenlines,
	      nexttileparams,sizeof(short *),sizeof(short));

  /* read unwrapped phase data */
  if(TMPTILEOUTFORMAT==ALT_LINE_DATA){
    ReadAltLineFilePhase(nextunwphaseptr,nexttileoutfiles->outfile,
			 nexttilelinelen,nexttilenlines,nexttileparams);
  }else if(TMPTILEOUTFORMAT==FLOAT_DATA){
    Read2DArray((void ***)nextunwphaseptr,nexttileoutfiles->outfile,
		nexttilelinelen,nexttilenlines,nexttileparams,
		sizeof(float *),sizeof(float));
  }else{
    fprintf(sp0,"Cannot read format of unwrapped phase tile data\nAbort\n");
    exit(ABNORMAL_EXIT);
  }

  /* read cost data */
  if(params->p<0){
    Read2DRowColFile((void ***)nextcostsptr,nexttileoutfiles->costoutfile,