snaphu_cost.c

phase unwrapping algorithm for SAR interferometry

  double costscale;
  double nshortcycle, nshortcyclesq;
  float **dpsi, **avgdpsi;
  smoothcostT **costs, **rowcost, **colcost;

  /* get memory and set cost array pointers */
  costs=(smoothcostT **)Get2DRowColMem(nrow,ncol,sizeof(smoothcostT *),
				       sizeof(smoothcostT));
  rowcost=(smoothcostT **)costs;
  colcost=(smoothcostT **)&costs[nrow-1];

  /* set up */
  rho0=(params->rhosconst1)/(params->ncorrlooks)+(params->rhosconst2);
  defocorrthresh=params->defothreshfactor*rho0;
  rhopow=2*(params->cstd1)+(params->cstd2)*log(params->ncorrlooks)
    +(params->cstd3)*(params->ncorrlooks);
  sigsqrhoconst=2.0/12.0;
  sigsqcorr=params->sigsqcorr;
  sigsqshortmin=params->sigsqshortmin;
  kperpdpsi=params->kperpdpsi;
  kpardpsi=params->kpardpsi;
  costscale=params->costscale; 
  nshortcycle=params->nshortcycle;
  nshortcyclesq=nshortcycle*nshortcycle;

  /* get memory for wrapped difference arrays */
  dpsi=(float **)Get2DMem(nrow,ncol,sizeof(float *),sizeof(float));
  avgdpsi=(float **)Get2DMem(nrow,ncol,sizeof(float *),sizeof(float));

  /* build array of mean wrapped phase differences in range */
  /* simple average of phase differences is biased, but mean phase */
  /*   differences usually near zero, so don't bother with complex average */
  fprintf(sp2,"Building range cost arrays\n");
  CalcWrappedRangeDiffs(dpsi,avgdpsi,wrappedphase,kperpdpsi,kpardpsi,
			nrow,ncol);

  /* build colcost array (range slopes) */
  for(col=0;col<ncol-1;col++){
    for(row=0;row<nrow;row++){

      /* see if we have a masked pixel */
      if(mag[row][col]==0 || mag[row][col+1]==0){
	  
	/* masked pixel */
	colcost[row][col].offset=0;
	colcost[row][col].sigsq=LARGESHORT;

      }else{

	/* deformation-mode costs */
	
	/* calculate variance due to decorrelation */
	/* need symmetry for range if deformation */
	rho=(corr[row][col]+corr[row][col+1])/2.0;
	if(rho<defocorrthresh){
	  rho=0;
	}
	sigsqrho=(sigsqrhoconst*pow(1-rho,rhopow)+sigsqcorr)*nshortcyclesq;

	/* set cost paramaters in terms of flow, represented as shorts */
	if(rho>0){
	  colcost[row][col].offset=nshortcycle*
	    (dpsi[row][col]-avgdpsi[row][col]);
	}else{
	  colcost[row][col].offset=nshortcycle*
	    (dpsi[row][col]-0.5*avgdpsi[row][col]);       
	}
	colcost[row][col].sigsq=sigsqrho/(costscale*colweight[row][col]);
	if(colcost[row][col].sigsq<sigsqshortmin){
	  colcost[row][col].sigsq=sigsqshortmin;
	}
      }

      /* shift PDF to account for flattening by coarse unwrapped estimate */
      if(unwrappedest!=NULL){
	colcost[row][col].offset+=(nshortcycle/TWOPI*
				   (unwrappedest[row][col+1]
				    -unwrappedest[row][col]));
      }
    }
  }  /* end of range gradient cost calculation */

  /* build array of mean wrapped phase differences in azimuth */
  /* biased, but not much, so don't bother with complex averaging */
  fprintf(sp2,"Building azimuth cost arrays\n");
  CalcWrappedAzDiffs(dpsi,avgdpsi,wrappedphase,kperpdpsi,kpardpsi,
		     nrow,ncol);

  /* build rowcost array */
  for(col=0;col<ncol;col++){
    for(row=0;row<nrow-1;row++){

      /* see if we have a masked pixel */
      if(mag[row][col]==0 || mag[row+1][col]==0){
	  
	/* masked pixel */
	rowcost[row][col].offset=0;
	rowcost[row][col].sigsq=LARGESHORT;

      }else{

	/* deformation-mode costs */

	/* variance due to decorrelation */
	/* get correlation and clip small values because of estimator bias */
	rho=(corr[row][col]+corr[row+1][col])/2.0;
	if(rho<defocorrthresh){
	  rho=0;
	}
        sigsqrho=(sigsqrhoconst*pow(1-rho,rhopow)+sigsqcorr)*nshortcyclesq;

	/* set cost paramaters in terms of flow, represented as shorts */
	if(rho>0){
	  rowcost[row][col].offset=nshortcycle*
	    (dpsi[row][col]-avgdpsi[row][col]);
	}else{
	  rowcost[row][col].offset=nshortcycle*
	    (dpsi[row][col]-0.5*avgdpsi[row][col]);
	}
	rowcost[row][col].sigsq=sigsqrho/(costscale*rowweight[row][col]);
	if(rowcost[row][col].sigsq<sigsqshortmin){
	  rowcost[row][col].sigsq=sigsqshortmin;
	}
      }

      /* shift PDF to account for flattening by coarse unwrapped estimate */
      if(unwrappedest!=NULL){
	rowcost[row][col].offset+=(nshortcycle/TWOPI*
				   (unwrappedest[row+1][col]
				    -unwrappedest[row][col]));
      }
    }
  } /* end of azimuth cost calculation */

  /* free temporary arrays */
  Free2DArray((void **)corr,nrow);
  Free2DArray((void **)dpsi,nrow);
  Free2DArray((void **)avgdpsi,nrow);

  /* return pointer to costs arrays */
  return((void **)costs);

}


/* function: GetIntensityAndCorrelation()
 * --------------------------------------
 * Reads amplitude and correlation info from files if specified.  If ampfile
 * not given, uses interferogram magnitude.  If correlation file not given,
 * generates correlatin info from interferogram and amplitude.
 */
void GetIntensityAndCorrelation(float **mag, float **wrappedphase, 
				float ***pwrptr, float ***corrptr, 
				infileT *infiles, long linelen, long nlines,
				long nrow, long ncol, outfileT *outfiles, 
				paramT *params, tileparamT *tileparams){

  long row, col, krowcorr, kcolcorr, iclipped;
  float **pwr, **corr;
  float **realcomp, **imagcomp;
  float **padreal, **padimag, **avgreal, **avgimag;
  float **pwr1, **pwr2, **padpwr1, **padpwr2, **avgpwr1, **avgpwr2; 
  double rho0, rhomin, biaseddefaultcorr;

  /* initialize */
  pwr=NULL;
  corr=NULL;
  pwr1=NULL;
  pwr2=NULL;
  
  /* read intensity, if specified */
  if(strlen(infiles->ampfile)){
    ReadIntensity(&pwr,&pwr1,&pwr2,infiles,linelen,nlines,params,tileparams);
  }else{
    if(params->costmode==TOPO){
      fprintf(sp1,"No brightness file specified.  ");
      fprintf(sp1,"Using interferogram magnitude as intensity\n");
    }
    pwr=(float **)Get2DMem(nrow,ncol,sizeof(float *),sizeof(float));
    for(row=0;row<nrow;row++){
      for(col=0;col<ncol;col++){
	pwr[row][col]=mag[row][col];
      }
    }
  }

  /* read corrfile, if specified */
  if(strlen(infiles->corrfile)){
    ReadCorrelation(&corr,infiles,linelen,nlines,tileparams); 
  }else if(pwr1!=NULL && pwr2!=NULL && params->havemagnitude){

    /* generate the correlation info from the interferogram and amplitude */
    fprintf(sp1,"Generating correlation from interferogram and intensity\n");

    /* get the correct number of looks, and make sure its odd */
    krowcorr=1+2*floor(params->ncorrlooksaz/(double )params->nlooksaz/2);
    kcolcorr=1+2*floor(params->ncorrlooksrange/(double )params->nlooksrange/2);

    /* calculate equivalent number of independent looks */
    params->ncorrlooks=(kcolcorr*(params->dr/params->rangeres))
      *(krowcorr*(params->da/params->azres))*params->nlooksother;
    fprintf(sp1,"   (%.1f equivalent independent looks)\n",
	    params->ncorrlooks);
    
    /* get real and imaginary parts of interferogram */
    realcomp=(float **)Get2DMem(nrow,ncol,sizeof(float *),sizeof(float));
    imagcomp=(float **)Get2DMem(nrow,ncol,sizeof(float *),sizeof(float));
    for(row=0;row<nrow;row++){
      for(col=0;col<ncol;col++){
	realcomp[row][col]=mag[row][col]*cos(wrappedphase[row][col]);
	imagcomp[row][col]=mag[row][col]*sin(wrappedphase[row][col]);
      }
    }

    /* do complex spatial averaging on the interferogram */
    padreal=MirrorPad(realcomp,nrow,ncol,(krowcorr-1)/2,(kcolcorr-1)/2);
    padimag=MirrorPad(imagcomp,nrow,ncol,(krowcorr-1)/2,(kcolcorr-1)/2);
    if(padreal==realcomp || padimag==imagcomp){
      fprintf(sp0,"Correlation averaging box too large for input array size\n"
	      "Abort\n");
      exit(ABNORMAL_EXIT);
    }
    avgreal=realcomp;
    BoxCarAvg(avgreal,padreal,nrow,ncol,krowcorr,kcolcorr);
    avgimag=imagcomp;
    BoxCarAvg(avgimag,padimag,nrow,ncol,krowcorr,kcolcorr);
    Free2DArray((void **)padreal,nrow);
    Free2DArray((void **)padimag,nrow);

    /* spatially average individual SAR power images */
    padpwr1=MirrorPad(pwr1,nrow,ncol,(krowcorr-1)/2,(kcolcorr-1)/2);
    avgpwr1=pwr1;
    BoxCarAvg(avgpwr1,padpwr1,nrow,ncol,krowcorr,kcolcorr);
    padpwr2=MirrorPad(pwr2,nrow,ncol,(krowcorr-1)/2,(kcolcorr-1)/2);
    avgpwr2=pwr2;
    BoxCarAvg(avgpwr2,padpwr2,nrow,ncol,krowcorr,kcolcorr);
    Free2DArray((void **)padpwr1,nrow);
    Free2DArray((void **)padpwr2,nrow);

    /* build correlation data */
    corr=(float **)Get2DMem(nrow,ncol,sizeof(float *),sizeof(float));
    for(row=0;row<nrow;row++){
      for(col=0;col<ncol;col++){
	if(avgpwr1[row][col]<=0 || avgpwr2[row][col]<=0){
	  corr[row][col]=0.0;
	}else{
	  corr[row][col]=sqrt((avgreal[row][col]*avgreal[row][col]
			       +avgimag[row][col]*avgimag[row][col])
			      /(avgpwr1[row][col]*avgpwr2[row][col]));
	}
      }
    }

    /* free temporary arrays */
    Free2DArray((void **)avgreal,nrow);
    Free2DArray((void **)avgimag,nrow);
    Free2DArray((void **)avgpwr1,nrow);
    Free2DArray((void **)avgpwr2,nrow);
    pwr1=NULL;
    pwr2=NULL;

  }else{

    /* no file specified: set corr to default value */
    /* find biased default correlation using */
    /* inverse of unbias method used by BuildCostArrays() */
    corr=(float **)Get2DMem(nrow,ncol,sizeof(float *),sizeof(float));
    fprintf(sp1,"No correlation file specified.  Assuming correlation = %g\n",
	   params->defaultcorr);
    rho0=(params->rhosconst1)/(params->ncorrlooks)+(params->rhosconst2);
    rhomin=params->rhominfactor*rho0;
    if(params->defaultcorr>rhomin){
      biaseddefaultcorr=params->defaultcorr;
    }else{
      biaseddefaultcorr=0.0;
    }
    for(row=0;row<nrow;row++){
      for(col=0;col<ncol;col++){
	corr[row][col]=biaseddefaultcorr;
      }
    }
  }

  /* dump correlation data if necessary */
  if(strlen(outfiles->rawcorrdumpfile)){
    Write2DArray((void **)corr,outfiles->rawcorrdumpfile,
		 nrow,ncol,sizeof(float)); 
  }

  /* check correlation data validity */
  iclipped=0;
  for(row=0;row<nrow;row++){
    for(col=0;col<ncol;col++){
      if(!IsFinite(corr[row][col])){
	fprintf(sp0,"NaN or infinity found in correlation data\nAbort\n");
	exit(ABNORMAL_EXIT);
      }else if(corr[row][col]>1.0){
	if(corr[row][col]>1.001){
	  iclipped++;               /* don't warn for minor numerical errors */
	}
	corr[row][col]=1.0;
      }else if(corr[row][col]<0.0){
	if(corr[row][col]<-0.001){
	  iclipped++;               /* don't warn for minor numerical errors */
	}
	corr[row][col]=0.0;
      }
    }
  }
  if(iclipped){
    fprintf(sp0,"WARNING: %ld illegal correlation values clipped to [0,1]\n",
	    iclipped);
  }
  
  /* dump correlation data if necessary */
  if(strlen(outfiles->corrdumpfile)){
    Write2DArray((void **)corr,outfiles->corrdumpfile,
		 nrow,ncol,sizeof(float)); 
  }

  /* free memory and set output pointers */
  if(pwr1!=NULL){
    Free2DArray((void **)pwr1,nrow);
  }
  if(pwr2!=NULL){
    Free2DArray((void **)pwr2,nrow);
  }
  if(params->costmode==DEFO && pwr!=NULL){
    Free2DArray((void **)pwr,nrow);
    pwr=NULL;
  }
  *pwrptr=pwr;
  *corrptr=corr;

}


/* function: RemoveMean()
 * -------------------------
 * Divides intensity by average over sliding window.
 */
void RemoveMean(float **ei, long nrow, long ncol, 
		       long krowei, long kcolei){

  float **avgei, **padei;
  long row, col;

  /* make sure krowei, kcolei are odd */
  if(!(krowei % 2)){
    krowei++;
  }
  if(!(kcolei % 2)){
    kcolei++;
  }

  /* get memory */
  avgei=(float **)Get2DMem(nrow,ncol,sizeof(float *),sizeof(float));

  /* pad ei in new array */
  padei=MirrorPad(ei,nrow,ncol,(krowei-1)/2,(kcolei-1)/2);
  if(padei==ei){
    fprintf(sp0,"Intensity-normalization averaging box too large "
	    "for input array size\nAbort\n");
    exit(ABNORMAL_EXIT);
  }

  /* calculate average ei by using sliding window */
  BoxCarAvg(avgei,padei,nrow,ncol,krowei,kcolei);

  /* divide ei by avgei */
  for(row=0;row<nrow;row++){
    for(col=0;col<ncol;col++){
      if(avgei){
	ei[row][col]/=(avgei[row][col]);
      }
    }
  }

  /* free memory */
  Free2DArray((void **)padei,nrow+krowei-1);
  Free2DArray((void **)avgei,nrow);

}


/* function: BuildDZRCritLookupTable()
 * -----------------------------------
 * Builds a 1-D lookup table of dzrcrit values indexed by incidence angle 
 * (in rad).
 */
float *BuildDZRCritLookupTable(double *nominc0ptr, double *dnomincptr, 
			       long *tablesizeptr, tileparamT *tileparams, 
			       paramT *params){

  long tablesize, k;
  double nominc, nominc0, nomincmax, dnominc;
  double a, re, slantrange;
  float *dzrcrittable;
  
  /* compute nominal spherical earth incidence angle for near and far range */
  a=params->orbitradius;
  re=params->earthradius;
  slantrange=params->nearrange+params->dr*tileparams->firstcol;
  nominc0=acos((a*a-slantrange*slantrange-re*re)/(2*slantrange*re));
  slantrange+=params->dr*tileparams->ncol;
  nomincmax=acos((a*a-slantrange*slantrange-re*re)/(2*slantrange*re));
  if(!IsFinite(nominc0) || !IsFinite(nomincmax)){
    fprintf(sp0,"Geometry error detected.  " 
	    "Check altitude, near range, and earth radius parameters\n"
	    "Abort\n");
    exit(ABNORMAL_EXIT);
  }

  /* build lookup table */
  dnominc=params->dnomincangle;
  tablesize=(long )floor((nomincmax-nominc0)/dnominc)+1;
  dzrcrittable=MAlloc(tablesize*sizeof(float));
  nominc=nominc0;
  for(k=0;k<tablesize;k++){
    dzrcrittable[k]=(float )SolveDZRCrit(sin(nominc),cos(nominc),params,
				 params->threshold);
    nominc+=dnominc;
    if(nominc>PI/2.0){
      nominc-=dnominc;
    }
  }
  
  /* set return variables */
  (*nominc0ptr)=nominc;
  (*dnomincptr)=dnominc;
  (*tablesizeptr)=tablesize;
  return(dzrcrittable);

}


/* function: SolveDZRCrit()
 * ------------------------