snaphu_util.c

phase unwrapping algorithm for SAR interferometry

 */
signed char IsFinite(double d){

  return(finite(d));
  /* return(isfinite(d)); */
  /* return(!(isnan(d) || isinf(d))); */
  /* return(TRUE) */
}


/* function: LRound()
 * ------------------
 * Rounds a floating point number to the nearest integer.  
 * The function takes a float and returns a long.
 */
long LRound(double a){

  return((long )rint(a));
}


/* function: Short2DRowColAbsMax()
 * -------------------------------
 * Returns the maximum of the absolute values of element in a 
 * two-dimensional short array.  The number of rows and columns 
 * should be passed in.
 */
long Short2DRowColAbsMax(short **arr, long nrow, long ncol){

  long row, col, maxval;

  maxval=0;
  for(row=0;row<nrow-1;row++){
    for(col=0;col<ncol;col++){
      if(labs(arr[row][col])>maxval){
	maxval=labs(arr[row][col]);
      }
    }
  }
  for(row=nrow-1;row<2*nrow-1;row++){
    for(col=0;col<ncol-1;col++){
      if(labs(arr[row][col])>maxval){
	maxval=labs(arr[row][col]);
      }
    }
  }
  return(maxval);
}


/* function: LinInterp1D()
 * -----------------------
 * Given an array of floats, interpolates at the specified noninteger
 * index.  Returns first or last array value if index is out of bounds.
 */
float LinInterp1D(float *arr, double index, long nelem){

  long intpart;
  double fracpart;

  intpart=(long )floor(index);
  fracpart=index-intpart;
  if(intpart<0){
    return(arr[0]);
  }else if(intpart>=nelem-1){
    return(arr[nelem-1]);
  }else{
    return(((1-fracpart)*arr[intpart]+fracpart*arr[intpart+1])/2.0);
  }
}


/* function: LinInterp2D()
 * -----------------------
 * Given a 2-D array of floats, interpolates at the specified noninteger
 * indices.  Returns first or last array values if index is out of bounds.
 */
float LinInterp2D(float **arr, double rowind, double colind , 
                  long nrow, long ncol){

  long rowintpart;
  double rowfracpart;

  rowintpart=(long )floor(rowind);
  rowfracpart=rowind-rowintpart;
  if(rowintpart<0){
    return(LinInterp1D(arr[0],colind,ncol));
  }else if(rowintpart>=nrow-1){
    return(LinInterp1D(arr[nrow-1],colind,ncol));
  }else{
    return(((1-rowfracpart)*LinInterp1D(arr[rowintpart],colind,ncol)
            +rowfracpart*LinInterp1D(arr[rowintpart+1],colind,ncol))/2.0);
  }
}


/* function: Despeckle()
 * ---------------------
 * Filters magnitude/power data with adaptive geometric filter to get rid of 
 * speckle.  Allocates 2D memory for ei.  Does not square before averaging.
 */
void Despeckle(float **mag, float ***ei, long nrow, long ncol){

  float **intensity;
  double ratio, ratiomax, wfull, wstick, w[NARMS+1];
  long row, col, i, j, k, Irow, Icol;
  short jmin[5]={2,2,0,1,2};
  short jmax[5]={2,3,4,3,2};
  enum{ C=0, T, B, R, L, TR, BL, TL, BR};

  /* get memory for output array */
  if(*ei==NULL){
    (*ei)=(float **)Get2DMem(nrow,ncol,sizeof(float *),sizeof(float));
  }

  /* pad magnitude and place into new array (don't touch original data) */
  intensity=MirrorPad(mag,nrow,ncol,ARMLEN,ARMLEN);
  if(intensity==mag){
    fprintf(sp0,"Despeckling box size too large for input array size\n"
	    "Abort\n");
    exit(ABNORMAL_EXIT);
  }


  /* do the filtering */
  for(row=0;row<nrow;row++){
    Irow=row+ARMLEN;
    for(col=0;col<ncol;col++){
      Icol=col+ARMLEN;
      
      /* filter only if input is nonzero so we preserve mask info in input */
      if(intensity[Irow][Icol]==0){

	(*ei)[row][col]=0;

      }else{

	for(k=0;k<NARMS+1;k++){
	  w[k]=0;
	}
	for(i=-1;i<=1;i++){
	  for(j=-1;j<=1;j++){
	    w[C]+=intensity[Irow+i][Icol+j];
	  }
	}
	for(i=-1;i<=1;i++){
	  for(j=2;j<ARMLEN+1;j++){
	    w[T]+=intensity[Irow-j][Icol+i];
	    w[B]+=intensity[Irow+j][Icol+i];
	    w[L]+=intensity[Irow+i][Icol-j];
	    w[R]+=intensity[Irow+i][Icol+j];
	  }
	}
	for(i=0;i<=4;i++){
	  for(j=jmin[i];j<=jmax[i];j++){
	    w[TR]+=intensity[Irow-i][Icol+j];
	    w[BR]+=intensity[Irow+i][Icol+j];
	    w[BL]+=intensity[Irow+i][Icol-j];
	    w[TL]+=intensity[Irow-i][Icol-j];
	  }
	}
	wfull=w[C]+w[T]+w[R]+w[B]+w[L];
	for(i=2;i<5;i++){
	  for(j=2;j<7-i;j++){
	    wfull+=intensity[Irow+i][Icol+j];
	    wfull+=intensity[Irow-i][Icol+j];
	    wfull+=intensity[Irow+i][Icol-j];
	    wfull+=intensity[Irow-i][Icol-j];
	  }
	} 
	ratiomax=1;
	for(k=1;k<=NARMS;k+=2){
	  wstick=w[0]+w[k]+w[k+1];
	  if((ratio=wstick/(wfull-wstick))<1){
	    ratio=1/ratio;
	  }
	  if(ratio>ratiomax){
	    ratiomax=ratio;
	    (*ei)[row][col]=wstick;
	  }
	}
      }
    }
  }   

  /* free memory */
  Free2DArray((void **)intensity,nrow+2*ARMLEN);
}



/* function: MirrorPad()
 * ---------------------
 * Returns pointer to 2D array where passed array is in center and
 * edges are padded by mirror reflections.  If the pad dimensions are
 * too large for the array size, a pointer to the original array is
 * returned.
 */
float **MirrorPad(float **array1, long nrow, long ncol, long krow, long kcol){

  long row, col;
  float **array2;

  /* get memory */
  array2=(float **)Get2DMem(nrow+2*krow,ncol+2*kcol,
			    sizeof(float *),sizeof(float));
  
  /* center array1 in new array */
  for(row=0;row<nrow;row++){
    for(col=0;col<ncol;col++){
      array2[row+krow][col+kcol]=array1[row][col];
    }
  }

  /* make sure pad is less than image dimensions */
  /* if array is too small for pad dimensions, return original array */
  /* requires checking by calling function */
  if(krow>nrow || kcol>ncol){
    return(array1);
  }

  /* mirror reflect edges */
  for(row=0;row<krow;row++){
    for(col=0;col<kcol;col++){
      array2[row][col]=array2[2*krow-row][2*kcol-col];
      array2[row][ncol+kcol+col]
        =array2[2*krow-row][ncol+kcol-2-col];
      array2[nrow+krow+row][col]
        =array2[nrow+krow-2-row][2*kcol-col];
      array2[nrow+krow+row][ncol+kcol+col]
        =array2[nrow+krow-2-row][ncol+kcol-2-col];
    }
  }
  for(row=krow;row<nrow+krow;row++){
    for(col=0;col<kcol;col++){
      array2[row][col]=array2[row][2*kcol-col];
      array2[row][ncol+kcol+col]
        =array2[row][ncol+kcol-2-col];
    }
  }
  for(col=kcol;col<ncol+kcol;col++){
    for(row=0;row<krow;row++){
      array2[row][col]=array2[2*krow-row][col];
      array2[nrow+krow+row][col]
        =array2[nrow+krow-2-row][col];
    }
  }
  
  /* return a pointer to the padded array */
  return(array2);

}


/* function: BoxCarAvg()
 * ---------------------
 * Takes in 2-D array, convolves with boxcar filter of size specified.  
 * Uses a recursion technique (but the function does not actually call
 * itself recursively) to compute the result, so there may be roundoff 
 * errors.
 */
void BoxCarAvg(float **avgarr, float **padarr, long nrow, long ncol, 
	       long krow, long kcol){

  long i, row, col, n;
  double window;

  /* loop over all rows */
  for(row=0;row<nrow;row++){

    /* calculate first cell */
    window=0;
    for(i=row;i<row+krow;i++){
      for(col=0;col<kcol;col++){
	window+=padarr[i][col];
      }
    }
    avgarr[row][0]=(float )window;

    /* convolve window with row, using result of last cell */
    for(col=1;col<ncol;col++){
      for(i=row;i<row+krow;i++){
	window-=padarr[i][col-1];
	window+=padarr[i][col+kcol-1];
      }
      avgarr[row][col]=(float )window;
    }
  }

  /* normalize */
  n=krow*kcol;
  for(row=0;row<nrow;row++){
    for(col=0;col<ncol;col++){
      avgarr[row][col]/=n;
    }
  }
}


/* function: StrNCopy()
 * --------------------
 * Just like strncpy(), but terminates string by putting a null 
 * character at the end (may overwrite last character).
 */
char *StrNCopy(char *dest, const char *src, size_t n){
 
  char *s;

  s=strncpy(dest,src,n-1);
  dest[n-1]='\0';
  return(s);
}


/* function: FlattenWrappedPhase()
 * -------------------------------
 * Given a wrapped phase array and an unwrapped phase array, subtracts
 * the unwrapped data elementwise from the wrapped data and stores
 * the result, rewrapped to [0,2pi), in the wrapped array.
 */
void FlattenWrappedPhase(float **wrappedphase, float **unwrappedest, 
			 long nrow, long ncol){
 
  long row, col;

  /* loop to subtract, rewrap, store in wrapped array. */
  for(row=0;row<nrow;row++){
    for(col=0;col<ncol;col++){
      wrappedphase[row][col]-=unwrappedest[row][col];
      wrappedphase[row][col]=fmod(wrappedphase[row][col],TWOPI);
      if(wrappedphase[row][col]<0){
	wrappedphase[row][col]+=TWOPI;
      }
    }
  }
}


/* function: Add2DFloatArrays()
 * ----------------------------
 * Addes the values of two 2-D arrays elementwise.
 */
void Add2DFloatArrays(float **arr1, float **arr2, long nrow, long ncol){

  long row, col;
 
  /* loop over all rows and columns, add and store result in first array */
  for(row=0;row<nrow;row++){
    for(col=0;col<ncol;col++){
      arr1[row][col]+=arr2[row][col];
    }
  }

}


/* function: StringToDouble()
 * --------------------------
 * Uses strtod to convert a string to a double, but also does error checking.
 * If any part of the string is not converted, the function does not make
 * the assignment and returns TRUE.  Otherwise, returns FALSE.
 */
int StringToDouble(char *str, double *d){

  double tempdouble;
  char *endp;
  
  endp=str;
  tempdouble=strtod(str,&endp);
  if(strlen(endp) || tempdouble>=HUGE_VAL || tempdouble<=-HUGE_VAL){
    return(TRUE);
  }else{
    *d=tempdouble;
    return(FALSE);
  }
}


/* function: StringToLong()
 * ------------------------
 * Uses strtol to convert a string to a base-10 long, but also does error 
 * checking.  If any part of the string is not converted, the function does 
 * not make the assignment and returns TRUE.  Otherwise, returns FALSE.
 */
int StringToLong(char *str, long *l){

  long templong;
  char *endp;
  
  endp=str;
  templong=strtol(str,&endp,10);
  if(strlen(endp) || templong==LONG_MAX || templong==LONG_MIN){
    return(TRUE);
  }else{
    *l=templong;
    return(FALSE);
  }
}


/* function: CatchSignals()
 * ------------------------
 * Traps common signals that by default cause the program to abort.  
 * Sets (pointer to function) Handler as the signal handler for all.
 * Note that SIGKILL usually cannot be caught.  No return value.
 */
void CatchSignals(void (*SigHandler)(int)){

  signal(SIGHUP,SigHandler);
  signal(SIGINT,SigHandler);
  signal(SIGQUIT,SigHandler);
  signal(SIGILL,SigHandler);
  signal(SIGABRT,SigHandler);
  signal(SIGFPE,SigHandler);
  signal(SIGSEGV,SigHandler);
  signal(SIGPIPE,SigHandler);
  signal(SIGALRM,SigHandler);
  signal(SIGTERM,SigHandler);
  signal(SIGBUS,SigHandler);
}


/* function: SetDump()
 * -------------------
 * Set the global variable dumpresults_global to TRUE if SIGINT or SIGHUP
 * signals recieved.  Also sets requestedstop_global if SIGINT signal 
 * received.  This function should only be called via signal() when 
 * a signal is caught.
 */ 
void SetDump(int signum){

  if(signum==SIGINT){

    /* exit if we receive another interrupt */
    signal(SIGINT,exit);

    /* print nice message and set global variables so program knows to exit */
    fprintf(sp0,"\n\nSIGINT signal caught.  Please wait for graceful exit\n");
    fprintf(sp0,"(One more interrupt signal halts job)\n");
    dumpresults_global=TRUE;
    requestedstop_global=TRUE;

  }else if(signum==SIGHUP){

    /* make sure the hangup signal doesn't revert to default behavior */
    signal(SIGHUP,SetDump);

    /* print a nice message, and set the dump variable */
    fprintf(sp0,"\n\nSIGHUP signal caught.  Dumping results\n");
    dumpresults_global=TRUE;

  }else{
    fprintf(sp0,"WARNING: Invalid signal (%d) passed to signal handler\n",
	    signum);
  }
}


/* function: KillChildrenExit()
 * ----------------------------
 * Signal handler that sends a KILL signal to all processes in the group
 * so that children exit when parent exits.
 */
void KillChildrenExit(int signum){

  fprintf(sp0,"Parent received signal %d\nKilling children and exiting\n",
	  signum);
  fflush(NULL);
  signal(SIGTERM,SIG_IGN);
  kill(0,SIGTERM);
  exit(ABNORMAL_EXIT);

}


/* function: SignalExit()
 * ----------------------
 * Signal hanlder that prints message about the signal received, then exits.
 */
void SignalExit(int signum){
  
  signal(SIGTERM,SIG_IGN);
  fprintf(sp0,"Exiting with status %d on signal %d\n",ABNORMAL_EXIT,signum);
  fflush(NULL);
  exit(ABNORMAL_EXIT);

}


/* function: StartTimers()
 * -----------------------
 * Starts the wall clock and CPU timers for use in conjunction with
 * DisplayElapsedTime().
 */
void StartTimers(time_t *tstart, double *cputimestart){
  
  struct rusage usagebuf;

  *tstart=time(NULL);
  *cputimestart=-1.0;
  if(!getrusage(RUSAGE_SELF,&usagebuf)){
    *cputimestart=(double )(usagebuf.ru_utime.tv_sec
			    +(usagebuf.ru_utime.tv_usec/(double )1000000)
			    +usagebuf.ru_stime.tv_sec
			    +(usagebuf.ru_stime.tv_usec/(double )1000000));
    if(!getrusage(RUSAGE_CHILDREN,&usagebuf)){
      *cputimestart+=(double )(usagebuf.ru_utime.tv_sec
			       +(usagebuf.ru_utime.tv_usec/(double )1000000)
			       +usagebuf.ru_stime.tv_sec
			       +(usagebuf.ru_stime.tv_usec/(double )1000000));
    }
  }
}


/* function: DisplayElapsedTime()
 * ------------------------------
 * Displays the elapsed wall clock and CPU times for the process and its
 * children.  Times should be initialized at the start of the program with
 * StartTimers().  The code is written to show the total processor time
 * for the parent process and all of its children, but whether or not 
 * this is actually done depends on the implementation of the system time 
 * functions.
 */
void DisplayElapsedTime(time_t tstart, double cputimestart){

  double cputime, walltime, seconds;
  long hours, minutes;
  time_t tstop;
  struct rusage usagebuf;

  cputime=-1.0;
  if(!getrusage(RUSAGE_CHILDREN,&usagebuf)){
    cputime=(double )(usagebuf.ru_utime.tv_sec
		       +(usagebuf.ru_utime.tv_usec/(double )1000000)
		       +usagebuf.ru_stime.tv_sec
		       +(usagebuf.ru_stime.tv_usec/(double )1000000));
    if(!getrusage(RUSAGE_SELF,&usagebuf)){
      cputime+=(double )(usagebuf.ru_utime.tv_sec
			 +(usagebuf.ru_utime.tv_usec/(double )1000000)
			 +usagebuf.ru_stime.tv_sec
			 +(usagebuf.ru_stime.tv_usec/(double )1000000));
    }
  }
  tstop=time(NULL);
  if(cputime>0 && cputimestart>=0){
    cputime-=cputimestart;
    hours=(long )floor(cputime/3600);
    minutes=(long )floor((cputime-3600*hours)/60);
    seconds=cputime-3600*hours-60*minutes;
    fprintf(sp1,"Elapsed processor time:   %ld:%02ld:%05.2f\n",
	    hours,minutes,seconds);
  }
  if(tstart>0 && tstop>0){
    walltime=tstop-tstart;
    hours=(long )floor(walltime/3600);
    minutes=(long )floor((walltime-3600*hours)/60);
    seconds=walltime-3600*hours-60*minutes;
    fprintf(sp1,"Elapsed wall clock time:  %ld:%02ld:%02ld\n",
	    hours,minutes,(long )seconds);
  }
}


/* function: LongCompare()
 * -----------------------
 * Compares two long integers.  For use with qsort().
 */
int LongCompare(const void *c1, const void *c2){
  return((*((long *)c1))-(*((long *)c2)));
}