snaphu_util.c

phase unwrapping algorithm for SAR interferometry

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

  snaphu utility function 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: IsTrue()
 * ------------------
 * Returns TRUE if the string input is any of TRUE, True, true, 1, 
 * y, Y, yes, YES
 */
int IsTrue(char *str){

  if(!strcmp(str,"TRUE") || !strcmp(str,"true") || !strcmp(str,"True")
     || !strcmp(str,"1") || !strcmp(str,"y") || !strcmp(str,"Y")
     || !strcmp(str,"yes") || !strcmp(str,"YES") || !strcmp(str,"Yes")){
    return(TRUE);
  }else{
    return(FALSE);
  }
}


/* function: IsFalse()
 * ------------------
 * Returns FALSE if the string input is any of FALSE, False, false, 
 * 0, n, N, no, NO
 */
int IsFalse(char *str){

  if(!strcmp(str,"FALSE") || !strcmp(str,"false") || !strcmp(str,"False")
     || !strcmp(str,"0") || !strcmp(str,"n") || !strcmp(str,"N")
     || !strcmp(str,"no") || !strcmp(str,"NO") || !strcmp(str,"No")){
    return(TRUE);
  }else{
    return(FALSE);
  }
}


/* function: SetBoolenaSignedChar()
 * --------------------------------
 * Sets the value of a signed character based on the string argument passed.
 * Returns TRUE if the string was not a valid value, FALSE otherwise.
 */
signed char SetBooleanSignedChar(signed char *boolptr, char *str){

  if(IsTrue(str)){
    (*boolptr)=TRUE;
    return(FALSE);
  }else if(IsFalse(str)){
    (*boolptr)=FALSE;
    return(FALSE);
  }
  return(TRUE);
}


/* function: ModDiff()
 * -------------------
 * Computes floating point difference between two numbers.
 * f1 and f2 should be between [0,2pi).  The result is the 
 * modulo difference between (-pi,pi].  Assumes that
 * PI and TWOPI have been defined.  
 */  
double ModDiff(double f1, double f2){

  double f3;

  f3=f1-f2;
  if(f3>PI){
    f3-=TWOPI;
  }else if(f3<=-PI){
    f3+=TWOPI;
  }
  return(f3);
}


/* function: WrapPhase()
 * ---------------------
 * Makes sure the passed float array is properly wrapped into the [0,2pi)
 * interval.
 */
void WrapPhase(float **wrappedphase, long nrow, long ncol){

  long row, col;

  for(row=0;row<nrow;row++){
    for(col=0;col<ncol;col++){
      wrappedphase[row][col]-=TWOPI*floor(wrappedphase[row][col]/TWOPI);
    }
  }
}


/* function: CalcWrappedRangeDiffs()
 * ---------------------------------
 * Computes an array of wrapped phase differences in range (across rows).
 * Input wrapped phase array should be in radians.  Output is in cycles.
 */
void CalcWrappedRangeDiffs(float **dpsi, float **avgdpsi, float **wrappedphase,
			   long kperpdpsi, long kpardpsi,
			   long nrow, long ncol){
  long row, col;
  float **paddpsi;

  for(row=0;row<nrow;row++){
    for(col=0;col<ncol-1;col++){
      dpsi[row][col]=(wrappedphase[row][col+1]-wrappedphase[row][col])/TWOPI;
      if(dpsi[row][col]>=0.5){
	dpsi[row][col]-=1.0;
      }else if(dpsi[row][col]<-0.5){
	dpsi[row][col]+=1.0;
      }
    }
  }
  paddpsi=MirrorPad(dpsi,nrow,ncol-1,(kperpdpsi-1)/2,(kpardpsi-1)/2);
  if(paddpsi==dpsi){
    fprintf(sp0,"Wrapped-gradient averaging box too large "
	    "for input array size\nAbort\n");
    exit(ABNORMAL_EXIT);
  }
  BoxCarAvg(avgdpsi,paddpsi,nrow,ncol-1,kperpdpsi,kpardpsi);
  Free2DArray((void **)paddpsi,nrow+kperpdpsi-1);

}


/* function: CalcWrappedAzDiffs()
 * ---------------------------------
 * Computes an array of wrapped phase differences in range (across rows).
 * Input wrapped phase array should be in radians.  Output is in cycles.
 */
void CalcWrappedAzDiffs(float **dpsi, float **avgdpsi, float **wrappedphase,
			long kperpdpsi, long kpardpsi, long nrow, long ncol){
  long row, col;
  float **paddpsi;

  for(row=0;row<nrow-1;row++){
    for(col=0;col<ncol;col++){
      dpsi[row][col]=(wrappedphase[row][col]-wrappedphase[row+1][col])/TWOPI;
      if(dpsi[row][col]>=0.5){
	dpsi[row][col]-=1.0;
      }else if(dpsi[row][col]<-0.5){
	dpsi[row][col]+=1.0;
      }
    }
  }
  paddpsi=MirrorPad(dpsi,nrow-1,ncol,(kpardpsi-1)/2,(kperpdpsi-1)/2);
  if(paddpsi==dpsi){
    fprintf(sp0,"Wrapped-gradient averaging box too large "
	    "for input array size\nAbort\n");
    exit(ABNORMAL_EXIT);
  }
  BoxCarAvg(avgdpsi,paddpsi,nrow-1,ncol,kpardpsi,kperpdpsi);
  Free2DArray((void **)paddpsi,nrow-1+kpardpsi-1);

}


/* function: CycleResidue()
 * ------------------------
 * Computes the cycle array of a phase 2D phase array. Input arrays 
 * should be type float ** and signed char ** with memory pre-allocated.  
 * Numbers of rows and columns in phase array should be passed.  
 * Residue array will then have size nrow-1 x ncol-1.  Residues will
 * always be -1, 0, or 1 if wrapped phase is passed in.
 */
void CycleResidue(float **phase, signed char **residue, 
		  int nrow, int ncol){

  int row, col;
  float **rowdiff, **coldiff;

  rowdiff=(float **)Get2DMem(nrow-1,ncol,sizeof(float *),sizeof(float));
  coldiff=(float **)Get2DMem(nrow,ncol-1,sizeof(float *),sizeof(float));

  for(row=0;row<nrow-1;row++){
    for(col=0;col<ncol;col++){
      rowdiff[row][col]=ModDiff(phase[row+1][col],phase[row][col]);
    }
  }
  for(row=0;row<nrow;row++){
    for(col=0;col<ncol-1;col++){
      coldiff[row][col]=ModDiff(phase[row][col+1],phase[row][col]);
    }
  }

  for(row=0;row<nrow-1;row++){
    for(col=0;col<ncol-1;col++){
      residue[row][col]=(signed char)LRound((coldiff[row][col]
					     +rowdiff[row][col+1]
					     -coldiff[row+1][col]
					     -rowdiff[row][col])/TWOPI);
    }
  }

  Free2DArray((void **)rowdiff,nrow-1);
  Free2DArray((void **)coldiff,nrow);

}


/* function: CalcFlow()
 * --------------------
 * Calculates flow based on unwrapped phase data in a 2D array.  
 * Allocates memory for row and column flow arrays.
 */
void CalcFlow(float **phase, short ***flowsptr, long nrow, long ncol){

  long row, col;

  /* get memory for flow arrays */
  if((*flowsptr)==NULL){
    (*flowsptr)=(short **)Get2DRowColMem(nrow,ncol,
					 sizeof(short *),sizeof(short));
  }

  /* get row flows (vertical phase differences) */
  for(row=0;row<nrow-1;row++){
    for(col=0;col<ncol;col++){
      (*flowsptr)[row][col]=(short)LRound((phase[row][col]-phase[row+1][col])
					 /TWOPI);
    }
  }
  
  /* get col flows (horizontal phase differences) */
  for(row=0;row<nrow;row++){
    for(col=0;col<ncol-1;col++){
      (*flowsptr)[nrow-1+row][col]=(short)LRound((phase[row][col+1]
						  -phase[row][col])
						 /TWOPI);
    }
  }
}


/* function: IntegratePhase()
 * --------------------------
 * This function takes row and column flow information and integrates
 * wrapped phase to create an unwrapped phase field.  The unwrapped
 * phase field will be the same size as the wrapped field.  The array
 * rowflow should have size N-1xM and colflow size NxM-1 where the 
 * phase fields are NxM.  Output is saved to a file.
 */
void IntegratePhase(float **psi, float **phi, short **flows,
		    long nrow, long ncol){

  long row, col;
  short **rowflow, **colflow;
  rowflow=flows;
  colflow=&(flows[nrow-1]);
 
  /* set first element as seed */
  phi[0][0]=psi[0][0];

  /* integrate over first row */
  for(col=1;col<ncol;col++){
    phi[0][col]=phi[0][col-1]+(ModDiff(psi[0][col],psi[0][col-1])
      +colflow[0][col-1]*TWOPI);
  }

  /* integrate over columns */
  for(row=1;row<nrow;row++){
    for(col=0;col<ncol;col++){
      phi[row][col]=phi[row-1][col]+(ModDiff(psi[row][col],psi[row-1][col])
	-rowflow[row-1][col]*TWOPI);
    }
  }

}/* end IntegratePhase() */


/* function: ExtractFlow()
 * -----------------------
 * Given an unwrapped phase array, parse the data and find the flows.
 * Assumes only integer numbers of cycles have been added to get the 
 * unwrapped phase from the wrapped pase.  Gets memory and writes 
 * wrapped phase to passed pointer.  Assumes flows fit into short ints.
 */
float **ExtractFlow(float **unwrappedphase, short ***flowsptr, 
		    long nrow, long ncol){    

  long row, col;
  float **wrappedphase;
  
  /* get memory for wrapped phase array */
  wrappedphase=(float **)Get2DMem(nrow,ncol,sizeof(float *),sizeof(float));

  /* calculate wrapped phase */
  for(row=0;row<nrow;row++){
    for(col=0;col<ncol;col++){
      /* fmod() gives wrong results here (maybe because of float argument?) */
      /* wrappedphase[row][col]=fmod(unwrappedphase[row][col],TWOPI); */
      wrappedphase[row][col]=unwrappedphase[row][col]
	-TWOPI*floor(unwrappedphase[row][col]/TWOPI);
    }
  }

  /* calculate flows */
  CalcFlow(unwrappedphase,flowsptr,nrow,ncol);

  /* return pointer to wrapped phase array */
  return(wrappedphase);

}


/* function: FlipPhaseArraySign()
 * ------------------------------
 * Flips the sign of all values in a passed array if the flip flag is set.
 * Otherwise, does nothing.
 */
void FlipPhaseArraySign(float **arr, paramT *params, long nrow, long ncol){

  long row, col;

  if(params->flipphasesign){
    for(row=0;row<nrow;row++){
      for(col=0;col<ncol;col++){
	arr[row][col]*=-1;
      }
    }
  }
}


/* function: FlipFlowArraySign()
 * ------------------------------
 * Flips the sign of all values in a row-by-column array if the flip
 * flip flag is set.  Otherwise, does nothing.
 */
void FlipFlowArraySign(short **arr, paramT *params, long nrow, long ncol){

  long row, col, maxcol;

  if(params->flipphasesign){
    for(row=0;row<2*nrow-1;row++){
      if(row<nrow-1){
	maxcol=ncol;
      }else{
	maxcol=ncol-1;
      }
      for(col=0;col<maxcol;col++){
	arr[row][col]=-arr[row][col];
      }
    }
  }
}


/* function: Get2DMem()
 * --------------------
 * Allocates memory for 2D array.
 * Dynamically allocates memory and returns pointer of
 * type void ** for an array of size nrow x ncol.  
 * First index is row number: array[row][col]
 * size is size of array element, psize is size of pointer
 * to array element (eg sizeof(float *)).
 */
void **Get2DMem(int nrow, int ncol, int psize, size_t size){

  int row;
  void **array;

  if((array=malloc(nrow*psize))==NULL){
    fprintf(sp0,"Out of memory\n"); 
    exit(ABNORMAL_EXIT);
  }
  for(row=0; row<nrow; row++){
    if((array[row]=malloc(ncol*size))==NULL){
      fprintf(sp0,"Out of memory\n");
      exit(ABNORMAL_EXIT);
    }
  }
  return(array);
}


/* function: Get2DRowColMem()
 * --------------------------
 * Allocates memory for 2D array.  The array will have 2*nrow-1 rows.
 * The first nrow-1 rows will have ncol columns, and the rest will
 * have ncol-1 columns.
 */
void **Get2DRowColMem(long nrow, long ncol, int psize, size_t size){

  long row;
  void **array;

  if((array=malloc((2*nrow-1)*psize))==NULL){
    fprintf(sp0,"Out of memory\n"); 
    exit(ABNORMAL_EXIT);
  }
  for(row=0; row<nrow-1; row++){
    if((array[row]=malloc(ncol*size))==NULL){
      fprintf(sp0,"Out of memory\n");
      exit(ABNORMAL_EXIT);
    }
  }
  for(row=nrow-1; row<2*nrow-1; row++){
    if((array[row]=malloc((ncol-1)*size))==NULL){
      fprintf(sp0,"Out of memory\n");
      exit(ABNORMAL_EXIT);
    }
  }
  return(array);
}


/* function: Get2DRowColZeroMem()
 * ------------------------------
 * Allocates memory for 2D array.  The array will have 2*nrow-1 rows.
 * The first nrow-1 rows will have ncol columns, and the rest will
 * have ncol-1 columns.  Memory is initialized to zero.
 */
void **Get2DRowColZeroMem(long nrow, long ncol, int psize, size_t size){

  long row;
  void **array;

  if((array=malloc((2*nrow-1)*psize))==NULL){
    fprintf(sp0,"Out of memory\n"); 
    exit(ABNORMAL_EXIT);
  }
  for(row=0; row<nrow-1; row++){
    if((array[row]=calloc(ncol,size))==NULL){
      fprintf(sp0,"Out of memory\n");
      exit(ABNORMAL_EXIT);
    }
  }
  for(row=nrow-1; row<2*nrow-1; row++){
    if((array[row]=calloc((ncol-1),size))==NULL){
      fprintf(sp0,"Out of memory\n");
      exit(ABNORMAL_EXIT);
    }
  }
  return(array);
}


/* function: MAlloc()
 * --------------------
 * Has same functionality as malloc(), but exits if out of memory.
 */
void *MAlloc(size_t size){

  void *ptr;

  if((ptr=malloc(size))==NULL){
    fprintf(sp0,"Out of memory\n");
    exit(ABNORMAL_EXIT);
  }
  return(ptr);
}


/* function: CAlloc()
 * ------------------
 * Has same functionality as calloc(), but exits if out of memory.
 */
void *CAlloc(size_t nitems, size_t size){
  
  void *ptr;
  
  if((ptr=calloc(nitems,size))==NULL){
    fprintf(sp0,"Out of memory\n");
    exit(ABNORMAL_EXIT);
  }
  return(ptr);
}


/* function: ReAlloc()
 * -------------------
 * Has same functionality as realloc(), but exits if out of memory.
 */
void *ReAlloc(void *ptr, size_t size){
  
  void *ptr2;
  
  if((ptr2=realloc(ptr,size))==NULL){
    fprintf(sp0,"Out of memory\n");
    exit(ABNORMAL_EXIT);
  }
  return(ptr2);
}


/* function: Free2DArray()
 * -----------------------
 * This function frees the dynamically allocated memory for a 2D
 * array.  Pass in a pointer to a pointer cast to a void **.
 * The function assumes the array is of the form arr[rows][cols]
 * so that nrow is the number of elements in the pointer array.
 */
void Free2DArray(void **array, unsigned int nrow){

  int row;

  for(row=0; row<nrow; row++){
    free(array[row]);
  }
  free(array);
}


/* function: Set2DShortArray()
 * -------------------------
 * Sets all entries of a 2D array of shorts to the given value.  Assumes
 * that memory is already allocated.
 */
void Set2DShortArray(short **arr, long nrow, long ncol, long value){

  long row, col;

  for(row=0;row<nrow;row++){
    for(col=0;col<ncol;col++){
      arr[row][col]=value;
    }
  }
}


/* function: ValidDataArray()
 * --------------------------
 * Given a 2D floating point array, returns FALSE if any elements are NaN
 * or infinite, and TRUE otherwise (uses math library finite() function).
 */
signed char ValidDataArray(float **arr, long nrow, long ncol){

  long row, col;

  for(row=0;row<nrow;row++){
    for(col=0;col<ncol;col++){
      if(!IsFinite(arr[row][col])){
	return(FALSE);
      }
    }
  }
  return(TRUE);
}


/* function: IsFinite()
 * --------------------
 * This function takes a double and returns a nonzero value if 
 * the arguemnt is finite (not NaN and not infinite), and zero otherwise.
 * Different implementations are given here since not all machines have
 * these functions available.