snaphu.c

phase unwrapping algorithm for SAR interferometry

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

  snaphu main 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"



/* global (external) variable definitions */

/* flags used for signal handling */
char dumpresults_global;
char requestedstop_global;

/* ouput stream pointers */
/* sp0=error messages, sp1=status output, sp2=verbose, sp3=verbose counter */
FILE *sp0, *sp1, *sp2, *sp3;

/* node pointer for marking arc not on tree in apex array */
/* this should be treated as a constant */
nodeT NONTREEARC[1];

/* pointers to functions which calculate arc costs */
void (*CalcCost)();
long (*EvalCost)();

/* pointers to functions for tailoring network solver to specific topologies */
nodeT *(*NeighborNode)();
void (*GetArc)();


/***************************/
/* main program for snaphu */
/***************************/

int main(int argc, char **argv){

  /* variable declarations */
  infileT infiles[1];
  outfileT outfiles[1];
  paramT params[1];
  time_t tstart;
  double cputimestart;
  long linelen, nlines;


  /* get current wall clock and CPU time */
  StartTimers(&tstart,&cputimestart);

  /* set output stream pointers (may be reset after inputs parsed) */
  SetStreamPointers();

  /* print greeting */
  fprintf(sp1,"\n%s v%s\n",PROGRAMNAME,VERSION);

  /* set default parameters */
  SetDefaults(infiles,outfiles,params);
  ReadConfigFile(DEF_SYSCONFFILE,infiles,outfiles,&linelen,params);

  /* parse the command line inputs */
  ProcessArgs(argc,argv,infiles,outfiles,&linelen,params);

  /* set verbose output if specified */
  SetVerboseOut(params);

  /* set names of dump files if necessary */
  SetDumpAll(outfiles,params);

  /* get number of lines in file */
  nlines=GetNLines(infiles,linelen);

  /* check validity of parameters */
  CheckParams(infiles,outfiles,linelen,nlines,params);

  /* log the runtime parameters */
  WriteConfigLogFile(argc,argv,infiles,outfiles,linelen,params);

  /* unwrap, forming tiles and reassembling if necessary */
  Unwrap(infiles,outfiles,params,linelen,nlines);
    
  /* finish up */
  fprintf(sp1,"Program %s done\n",PROGRAMNAME);
  DisplayElapsedTime(tstart,cputimestart);
  exit(NORMAL_EXIT);

} /* end of main() */


/* function: Unwrap()
 * ------------------
 * Sets parameters for each tile and calls UnwrapTile() to do the
 * unwrapping.
 */
void Unwrap(infileT *infiles, outfileT *outfiles, paramT *params, 
	    long linelen, long nlines){

  long nexttilerow, nexttilecol, ntilerow, ntilecol, nthreads, nchildren;
  long sleepinterval;
  tileparamT tileparams[1];
  outfileT tileoutfiles[1];
  pid_t pid;
  int childstatus;
  double tilecputimestart;
  time_t tiletstart;


  /* set up */
  ntilerow=params->ntilerow;
  ntilecol=params->ntilecol;
  nthreads=params->nthreads;
  dumpresults_global=FALSE;
  requestedstop_global=FALSE;


  /* do the unwrapping */
  if(ntilerow==1 && ntilecol==1){

    /* only single tile */

    /* do the unwrapping */
    tileparams->firstrow=params->piecefirstrow;
    tileparams->firstcol=params->piecefirstcol;
    tileparams->nrow=params->piecenrow;
    tileparams->ncol=params->piecencol;
    UnwrapTile(infiles,outfiles,params,tileparams,nlines,linelen);

  }else{

    /* don't unwrap if in assemble-only mode */
    if(!params->assembleonly){

      /* make a temporary directory into which tile files will be written */
      MakeTileDir(params,outfiles);

      /* different code for parallel or nonparallel operation */
      if(nthreads>1){

	/* parallel code */

	/* initialize */
	nexttilerow=0;
	nexttilecol=0;
	nchildren=0;
	sleepinterval=LRound(nlines*linelen/ntilerow/ntilecol*SECONDSPERPIXEL);

	/* trap signals so children get killed if parent dies */
	CatchSignals(KillChildrenExit);

	/* loop until we're done unwrapping */
	while(TRUE){

	  /* unwrap next tile if there are free processors and tiles left */
	  if(nchildren<nthreads && nexttilerow<ntilerow){
	    
	    /* fork to create new process */
	    fflush(NULL);
	    pid=fork();

	    /* see if parent or child (or error) */
	    if(pid<0){

	      /* parent kills children and exits if there was a fork error */
	      fprintf(sp0,"Error while forking\nAbort\n");
	      kill(0,SIGKILL);
	      exit(ABNORMAL_EXIT);

	    }else if(pid==0){

	      /* child executes this code after fork */

	      /* reset signal handlers so that children exit nicely */
	      CatchSignals(SignalExit);

	      /* start timers for this tile */
	      StartTimers(&tiletstart,&tilecputimestart);

	      /* set up tile parameters */
	      pid=getpid();
	      fprintf(sp1,"Unwrapping tile at row %ld, column %ld (pid %ld)\n",
		      nexttilerow,nexttilecol,(long )pid);
	      SetupTile(nlines,linelen,params,tileparams,outfiles,tileoutfiles,
			nexttilerow,nexttilecol);
	      
	      /* reset stream pointers for logging */
	      ChildResetStreamPointers(pid,nexttilerow,nexttilecol,params);

	      /* unwrap the tile */
	      UnwrapTile(infiles,tileoutfiles,params,tileparams,
			 nlines,linelen);

	      /* log elapsed time */
	      DisplayElapsedTime(tiletstart,tilecputimestart);

	      /* child exits when done unwrapping */
	      exit(NORMAL_EXIT);

	    }
	      
	    /* parent executes this code after fork */

	    /* increment tile counters */
	    nchildren++;
	    if(++nexttilecol==ntilecol){
	      nexttilecol=0;
	      nexttilerow++;
	    }

	    /* wait a little while for file i/o before beginning next tile */
	    sleep(sleepinterval);

	  }else{

	    /* wait for a child to finish (only parent gets here) */
	    pid=wait(&childstatus);

	    /* make sure child exited cleanly */
	    if(!(WIFEXITED(childstatus)) || (WEXITSTATUS(childstatus))!=0){
	      fprintf(sp0,"Unexpected or abnormal exit of child process %ld\n"
		      "Abort\n",(long )pid);
	      signal(SIGTERM,SIG_IGN);
	      kill(0,SIGTERM);
	      exit(ABNORMAL_EXIT);
	    }

	    /* we're done if there are no more active children */
	    if(--nchildren==0){
	      break;
	    }

	  } /* end if free processor and tiles remaining */
	} /* end while loop */

	/* return signal handlers to default behavior */
	CatchSignals(SIG_DFL);

      }else{

	/* nonparallel code */

	/* loop over all tiles */
	for(nexttilerow=0;nexttilerow<ntilerow;nexttilerow++){
	  for(nexttilecol=0;nexttilecol<ntilecol;nexttilecol++){
	    
	    /* set up tile parameters */
	    fprintf(sp1,"Unwrapping tile at row %ld, column %ld\n",
		    nexttilerow,nexttilecol);
	    SetupTile(nlines,linelen,params,tileparams,outfiles,tileoutfiles,
		      nexttilerow,nexttilecol);
	    
	    /* unwrap the tile */
	    UnwrapTile(infiles,tileoutfiles,params,tileparams,nlines,linelen);

	  }
	}

      } /* end if nthreads>1 */

    } /* end if !params->assembleonly */

    /* reassemble tiles */
    AssembleTiles(outfiles,params,nlines,linelen);

  } /* end if multiple tiles */

} /* end of Unwrap() */


/* function: UnwrapTile()
 * ----------------------
 * This is the main phase unwrapping function for a single tile.
 */
void UnwrapTile(infileT *infiles, outfileT *outfiles, paramT *params,