snaphu_tile.c

phase unwrapping algorithm for SAR interferometry

  nodeT *source, *from, *to, *ground;
  unsigned char **components;
  bucketT bkts[1];


  /* error checking */
  fprintf(sp1,"Growing connected component mask\n");
  minsize=params->minconncompfrac*nrow*ncol;
  maxncomps=params->maxncomps;
  costthresh=params->conncompthresh;
  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 */
  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=1;
  bkts->minind=0;
  bkts->maxind=0;
  bkts->wrapped=FALSE;
  bkts->bucketbase=(nodeT **)MAlloc(sizeof(nodeT *));
  bkts->bucket=bkts->bucketbase;
  bkts->bucket[0]=NULL;
  
  /* initialize region variables */
  regioncounter=0;
  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;

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

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

	/* 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(regionsizes[regioncounter]>=minsize){

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

	    }else{

	      /* no more nonregion nodes, but current region still too small */
	      /* zero out the region */
	      RenumberRegion(nodes,source,0,nrow,ncol);
	      regioncounter--;
	      break;

	    }
	  }

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

	  /* scan from's neighbors */
	  arcnum=0;
	  while((to=RegionsNeighborNode(from,&arcnum,nodes,
					&arcrow,&arccol,nrow,ncol))!=NULL){
	   
	    /* see if to can be reached */
	    if(to->incost<0 && incrcosts[arcrow][arccol].negcost==0 
	       && to->group!=INBUCKET){

	      /* update to node */
	      to->pred=from;
	      BucketInsert(to,0,bkts);

	    }
	  }
	}
      }
    }
  }
  fprintf(sp2,"%ld connected components formed\n",regioncounter);

  /* make sure we don't have too many components */
  if(regioncounter>maxncomps){

    /* copy regionsizes array and sort to find new minimum region size */
    fprintf(sp2,"Keeping only %ld connected components\n",maxncomps);
    sortedregionsizes=(long *)MAlloc(regioncounter*sizeof(long));
    for(i=0;i<regioncounter;i++){
      sortedregionsizes[i]=regionsizes[i+1];
    }
    qsort((void *)sortedregionsizes,regioncounter,sizeof(long),LongCompare);
    minsize=sortedregionsizes[regioncounter-maxncomps];

    /* see how many regions of size minsize still need zeroing */
    ntied=0;
    i=regioncounter-maxncomps-1;
    while(i>=0 && sortedregionsizes[i]==minsize){
      ntied++;
      i--;
    }

    /* zero out regions that are too small */
    newnum=-1;
    for(row=0;row<nrow;row++){
      for(col=0;col<ncol;col++){
	i=nodes[row][col].incost;
	if(i>0){
	  if(regionsizes[i]<minsize 
	     || (regionsizes[i]==minsize && (ntied--)>0)){

	    /* region too small, so zero it out */
	    RenumberRegion(nodes,&(nodes[row][col]),0,nrow,ncol);

	  }else{

	    /* keep region, assign it new region number */
	    /* temporarily assign negative of new number to avoid collisions */
	    RenumberRegion(nodes,&(nodes[row][col]),newnum--,nrow,ncol);

	  }
	}
      }
    }

    /* flip temporary negative region numbers so they are positive */
    for(row=0;row<nrow;row++){
      for(col=0;col<ncol;col++){
	nodes[row][col].incost=-nodes[row][col].incost;
      }
    }
  }

  /* write components array */
  components=(unsigned char **)Get2DMem(nrow,ncol,sizeof(unsigned char *),
					sizeof(unsigned char));
  for(row=0;row<nrow;row++){
    for(col=0;col<ncol;col++){
      if(nodes[row][col].incost>255){
	fprintf(sp0,"Number of connected components too large for byte data\n"
		"Abort\n");
	exit(ABNORMAL_EXIT);
      }
      components[row][col]=(unsigned char )(nodes[row][col].incost);
    }
  }
  fprintf(sp1,"Writing connected components to file %s\n",
	  outfiles->conncompfile);
  Write2DArray((void **)components,outfiles->conncompfile,nrow,ncol,
	       sizeof(unsigned char));

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

}


/* function: ThickenCosts()
 * ------------------------
 */
long ThickenCosts(incrcostT **incrcosts, long nrow, long ncol){

  long row, col, templong, maxcost;
  double n;
  

  /* initialize variable storing maximum cost */
  maxcost=-LARGELONG;

  /* loop over row arcs and convolve */
  for(row=0;row<nrow-1;row++){
    for(col=0;col<ncol;col++){
      templong=2*incrcosts[row][col].poscost;
      n=2.0;
      if(col!=0){
	templong+=incrcosts[row][col-1].poscost;
	n+=1.0;
      }
      if(col!=ncol-1){
	templong+=incrcosts[row][col+1].poscost;
	n+=1.0;
      }
      templong=LRound(templong/n);
      if(templong>LARGESHORT){
        fprintf(sp0,"WARNING: COSTS CLIPPED IN ThickenCosts()\n");
	incrcosts[row][col].negcost=LARGESHORT;
      }else{
	incrcosts[row][col].negcost=templong;
      }
      if(incrcosts[row][col].negcost>maxcost){
	maxcost=incrcosts[row][col].negcost;
      }
    }
  }

  /* loop over column arcs and convolve */
  for(row=nrow-1;row<2*nrow-1;row++){
    for(col=0;col<ncol-1;col++){
      templong=2*incrcosts[row][col].poscost;
      n=2.0;
      if(row!=nrow-1){
	templong+=incrcosts[row-1][col].poscost;
	n+=1.0;
      }
      if(row!=2*nrow-2){
	templong+=incrcosts[row+1][col].poscost;
	n+=1.0;
      }
      templong=LRound(templong/n);
      if(templong>LARGESHORT){
        fprintf(sp0,"WARNING: COSTS CLIPPED IN ThickenCosts()\n");
	incrcosts[row][col].negcost=LARGESHORT;
      }else{      
	incrcosts[row][col].negcost=templong;
      }
      if(incrcosts[row][col].negcost>maxcost){
	maxcost=incrcosts[row][col].negcost;
      }
    }
  }

  /* return maximum cost */
  return(maxcost);

}


/* function: RegionsNeighborNode()
 * -------------------------------
 * Return the neighboring node of the given node corresponding to the
 * given arc number.
 */
nodeT *RegionsNeighborNode(nodeT *node1, long *arcnumptr, nodeT **nodes, 
			   long *arcrowptr, long *arccolptr, 
			   long nrow, long ncol){
  
  long row, col;

  row=node1->row;
  col=node1->col;

  while(TRUE){
    switch((*arcnumptr)++){
    case 0:
      if(col!=ncol-1){
	*arcrowptr=nrow-1+row;
	*arccolptr=col;
	return(&nodes[row][col+1]);
      }
      break;
    case 1:
      if(row!=nrow-1){
      *arcrowptr=row;
      *arccolptr=col;
	return(&nodes[row+1][col]);
      }
      break;
    case 2:
      if(col!=0){
	*arcrowptr=nrow-1+row;
	*arccolptr=col-1;
	return(&nodes[row][col-1]);
      }
      break;
    case 3:
      if(row!=0){
	*arcrowptr=row-1;
	*arccolptr=col;
	return(&nodes[row-1][col]);
      }
      break;
    default:
      return(NULL);
    }
  }
}


/* function: ClearBuckets()
 * ------------------------
 * Removes any nodes in the bucket data structure passed, and resets
 * their distances to VERYFAR.  Assumes bukets indexed from 0.
 */
void ClearBuckets(bucketT *bkts){

  nodeT *currentnode, *nextnode;
  long i;

  /* loop over all buckets */
  for(i=0;i<bkts->size;i++){

    /* clear the bucket */
    nextnode=bkts->bucketbase[i];
    while(nextnode!=NULL){
      currentnode=nextnode;
      nextnode=currentnode->next;
      currentnode->group=NOTINBUCKET;
      currentnode->outcost=VERYFAR;
      currentnode->pred=NULL;
    }
    bkts->bucketbase[i]=NULL;
  }

  /* reset bucket parameters */
  bkts->minind=0;
  bkts->maxind=bkts->size-1;
  bkts->wrapped=FALSE;
}


/* function: MergeRegions()
 * ------------------------
 * 
 */
void MergeRegions(nodeT **nodes, nodeT *source, long *regionsizes, 
		  long closestregion, 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=closestregion;
    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;
      }
    }
  }

  /* update size of region to which we are merging */
  regionsizes[closestregion]+=regionsizes[regionnum];