utilities.hh

Go to the documentation of this file.

/*
 * Copyright (c) 1999-2003 Bert Kampes
 * Copyright (c) 1999-2003 Delft University of Technology, The Netherlands
 *
 * This file is part of Doris, the Delft o-o radar interferometric software.
 *
 * Doris program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * Doris is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 * Publications that contain results produced by the Doris software should
 * contain an acknowledgment. (For example: The interferometric processing
 * was performed using the freely available Doris software package developed
 * by the Delft Institute for Earth-Oriented Space Research (DEOS), Delft
 * University of Technology, or include a reference to: Bert Kampes and
 * Stefania Usai. \"Doris: The Delft Object-oriented Radar Interferometric
 * software.\" In: proceedings 2nd ITC ORS symposium, August 1999. (cdrom)).
 *
 */
/****************************************************************
 * $Source: /users/kampes/DEVELOP/DORIS/doris/src/RCS/utilities.hh,v $
 * $Revision: 3.12 $                                            *
 * $Date: 2005/04/11 13:47:45 $                                 *
 * $Author: kampes $                                            *
 *                                                              *
 * - Some utilities.                                            *
 ****************************************************************/


#ifndef UTILITIES_H                             // guard
#define UTILITIES_H

#include "constants.hh"                         // typedefs
#include "slcimage.hh"                          // my slc image class
#include "ioroutines.hh"                        // ERROR
#include "exceptions.hh"                        // exceptions class
#include <cmath>                                // floor, sqr


// === prototype requires below, so put first/.
// ______ Solve system of 3 equations ______
void solve33(
        matrix<real8> &Result,
        const matrix<real8> &rhs,
        const matrix<real8> &partials);


// ______ Solve set of 2 equations ______
matrix<real8> solve22(
        const matrix<real8> &rhs,
        const matrix<real8> &partials);




// ====== Inline ======
inline bool iseven(int16 w)             {return (w+1)%2;}
inline bool iseven(int32 w)             {return (w+1)%2;}
inline bool iseven(uint  w)             {return (w+1)%2;}
inline bool isodd (int16 w)             {return w%2;}
inline bool isodd (int32 w)             {return w%2;}
inline bool isodd (uint  w)             {return w%2;}

inline bool ispower2 (uint w)
  {if (w==2   || w==4    || w==8    || w==16   ||
       w==32  || w==64   || w==128  || w==256  ||
       w==512 || w==1024 || w==2048 || w==4096   )
     return true;
   return false;}

inline int32 Ncoeffs(int32 degree)
  {return int32(.5*(sqr(degree+1)+degree+1));}

inline int32 degree(int32 Ncoeffs)
  { return int32(.5*(-1 + int32(sqrt(real4(1+8*Ncoeffs)))))-1;}

inline real4 remainder(real4 number, real4 divisor)     
  {return number-floor(number/divisor)*divisor;}

inline real8 remainder(real8 number, real8 divisor)     
  {return number-floor(number/divisor)*divisor;}

inline real4 sinc(real4 x)      
  {return ((x==0) ? 1 : sin(PI*x)/(PI*x));}

//____RaffaeleNutricato START MODIFICATION SECTION 1
inline real8 sinc(real8 x)      
  {return ((x==0) ? 1 : sin(PI*x)/(PI*x));}
//____RaffaeleNutricato END MODIFICATION SECTION 1

inline real4 rect(real4 x)      
  {real4 ans = 0.;
   if (x<.5 && x>-.5) ans = 1;
   else if (x==.5 || x==-.5) ans = .5;
   return ans;}

inline real4 tri(real4 x)       
  {real4 ans = 0.;
   if (x<1. && x>-1.)
     {(x<0) ? ans=1+x : ans=1-x;}
   return ans;}

inline real8 onedecimal(real8 x)
  {return real8(int32(x*10+.5))/10.;}
inline real4 onedecimal(real4 x)
  {return real4(int32(x*10+.5))/10.;}


/****************************************************************
 * myrect                                                       *
 * rect window, lying vector                                    *
 * r(i) = 1 if abs(i)<=.5                                       *
 * used in filtering, see curlander, swabisch, geudtner         *
 *    Bert Kampes, 31-Mar-2000                                  *
 ****************************************************************/
inline matrix<real4> myrect(const matrix<real4> &X)
  {
  TRACE_FUNCTION("myrect (BK 31-Mar-2000)")
  #ifdef __DEBUG
    if (X.lines()!=1) 
      {
      PRINT_ERROR("myrect: only lying vectors.")
      throw(argument_error);
      }
  #endif
  matrix<real4> Res(1,X.pixels());      // init2zero
  real4 *pntX = X[0];
  real4 *pntR = Res[0];
  for (register int32 i=0; i<int32(X.pixels()); ++i)
    {
    if (abs(*pntX++)<=0.5)
      *pntR=1.;
    pntR++;
    }
  // BK 01-Nov-2000:  better: ??
  //for (real4 *pntX=&X[0][0]; pntX<=&X[X.lines()-1][X.pixels()-1]; ++pntX)
  //  *pntR++ = (abs(*pntX)<=0.5) ? 1 : 0;
  // rather something like: R(find(abs(X)<.5))=1;
  return Res;
  }


/****************************************************************
 * myhamming                                                    *
 * hamming window, lying vector                                 *
 * w = (a + (1.-a).*cos((2.*pi/fs).*fr)) .* myrect(fr./Br);     *
 * scale/shift filter by g(x)=f((x-xo)/s)                       *
 * alpha==1 yields a myrect window                              *
 *    Bert Kampes, 31-Mar-2000                                  *
 ****************************************************************/
inline matrix<real4> myhamming(
        const matrix<real4> &fr,
        real8 RBW,
        real8 RSR,
        real8 alpha)
  {
  TRACE_FUNCTION("myhamming (BK 31-Mar-200)")
  #ifdef __DEBUG
    if (fr.lines()!=1)
      {
      PRINT_ERROR("myhamming: only lying vectors.")
      throw(argument_error);
      }
    if (alpha<0.0 || alpha>1.0)
      {
      PRINT_ERROR("myhamming: !alpha e{0..1}.")
      throw(argument_error);
      }
    if (RBW>RSR)
      {
      PRINT_ERROR("myhamming: RBW>RSR.")
      throw(argument_error);
      }
  #endif
  matrix<real4> Res(1,fr.pixels());
  real4 *pntfr = fr[0];
  real4 *pntR  = Res[0];
  for (register int32 i=0; i<int32(fr.pixels()); ++i)
    {
    if (abs((*pntfr)/RBW)<0.5)          // rect window
      *pntR = alpha+(1.-alpha)*cos((2*PI/RSR)*(*pntfr));
    pntfr++;
    pntR++;
    }
  // w = (a + (1.-a).*cos((2.*pi/fs).*fr)) .* myrect(fr./Br);
  return Res;
  }


/****************************************************************
 *    interpbilinear                                            *
 * bilinear interpolation with NORMALIZED data on grid.         *
 * f(row,col) = a00 + a01*c + a10*r +a11*r*c                    *
 * P1(0,0); P2(0,1); P3(1,0); P4(1,1); P[0:1,0:1];              *
 *                                                              *
 *         cols->                                               *
 *      r P1    P2                                              *
 *      o    .P                                                 *
 *      w P3    P4                                              *
 *      s                                                       *
 *                                                              *
 * input:                                                       *
 *  - normalized coordinates of point (row,col).                *
 *  - functional values at 4 corners.                           *
 * output:                                                      *
 *  - interpolated value                                        *
 *                                                              *
 *    Bert Kampes, 19-Jan-2000                                  *
 #%// BK 25-Sep-2000
 ****************************************************************/
inline real8 interpbilinear(
        const real8 r,
        const real8 c,
        const real8 P1,
        const real8 P2,
        const real8 P3,
        const real8 P4)         
  {return P1 + (P2-P1)*c + (P3-P1)*r + (P4-P3-P2+P1)*r*c;}


/****************************************************************
 *    interp_trilinear                                          *
 * trilinear interpolation using relative coordinates.          *
 *                                                              *
 * ====== Points A,B,C and do trigrid lin. interpolation at "+" *
 * ====== "+", origin of local coordinate system                *
 *                                                              *
 *     A                                                        *
 *                                                              *
 *      +                                                       *
 *  B                                                           *
 *                                                              *
 *              C                                               *
 *                                                              *
 * ====== Then, cn A.z = f(x,y) = a00+a10*A.x+a01*A.y;          *
 * ======       cn B.z = f(x,y) = a00+a10*B.x+a01*B.y;          *
 * ======       cn C.z = f(x,y) = a00+a10*C.x+a01*C.y;          *
 * ====== solve for a00; the Z value at origin, i.e., it!       *
 *                                                              *
 * input:                                                       *
 *  - normalized coordinates of point (row,col).                *
 *  - cn A,B,C with relative coordinates and functional values  *
 * output:                                                      *
 *  - interpolated value                                        *
 *                                                              *
 * Bert Kampes, 07-Apr-2005                                     *
 ****************************************************************/
inline real8 interp_trilinear(
        const cn &A,
        const cn &B,
        const cn &C)
  {
  // --- check if points are on line or are all same... ---
  // --- maybe simply add eps to A.x, A.y and subtract it from B.x, B.y?
  // --- since how else can we guarantee solve33 won't crash?
  //if (A.x==B.x && A.y==B.y)
  //  if (A.x==C.x && A.y==C.y) 
  //    return (A.z+B.z+C.z)/3.0;
  //  else // 2 same points?
  //const real8 EPS = 0.00001;// assume points are pixels
  // --- Vandermonde matrix ---
//  matrix<real8> DESIGNMAT(3,3); 
//  DESIGNMAT(0,0)=1.0;  DESIGNMAT(0,1)=A.x + EPS;  DESIGNMAT(0,2) = A.y + EPS;
//  DESIGNMAT(1,0)=1.0;  DESIGNMAT(1,1)=B.x - EPS;  DESIGNMAT(1,2) = B.y - EPS;
//  DESIGNMAT(2,0)=1.0;  DESIGNMAT(2,1)=C.x;        DESIGNMAT(2,2) = C.y;
//  // --- Right hand side ---
//  matrix<real8> RHS(3,1); 
//  RHS(0,0)=A.z;
//  RHS(1,0)=B.z;
//  RHS(2,0)=C.z;
//  // --- Solve system using library for all coefficients ---
//  matrix<real8> RESULT(3,1);// returned
//  solve33(RESULT,RHS,DESIGNMAT);// return result = a00,a10,a01
//  const real8 a00 = RESULT(0,0);
//  return a00;
//
// since first the last element is solved, turn it around and do it locally:
//* ====== Then, cn A.z = f(x,y) = [ya,xa,1][a01]
//* ======       cn B.z = f(x,y) = [yb,xb,1][a10]
//* ======       cn C.z = f(x,y) = [yc,xc,1][a00] <--- solve this on by LU
//* ====== solve for a00; the Z value at origin, i.e., it!
// ______ LU decomposition ______
const real8 one = real8(1.0);
const real8 EPS = 0.00001;// assume points are pixels
const real8 L10 = (B.y-EPS)/(A.y+EPS);// DESIGN(1,0)/DESIGN(0,0);
const real8 L20 = (C.y)/(A.y+EPS);// DESIGN(2,0)/DESIGN(0,0);
const real8 U11 = (B.x-EPS)-L10*(A.x+EPS);// DESIGN(1,1)-L10*DESIGN(0,1);
const real8 L21 = ((C.x)-(A.x+EPS))*L20/U11;// (DESIGN(2,1)-(DESIGN(0,1)*L20))/U11;
const real8 U12 = one-L10;// DESIGN(1,2)-L10*DESIGN(0,2);
const real8 U22 = one-L20-L21*U12;// DESIGN(2,2)-L20*DESIGN(0,2)-L21*U12;

// ______Solution: forward substitution______
const real8 b0  = A.z;// rhs(0,0);
const real8 b1  = B.z-b0*L10;// rhs(1,0)-b0*L10;
const real8 b2  = C.z-b0*L20-b1*L21;// rhs(2,0)-b0*L20-b1*L21;

// ______Solution: backwards substitution______
const real8 a00 = b2/U22;//RESULT(2,0) = b2/U22;
//RESULT(1,0) = (b1-U12*RESULT(2,0))/U11;
//RESULT(0,0) = (b0-A(0,1)*RESULT(1,0)-A(0,2)*RESULT(2,0))/A(0,0);
return a00;
  }


/****************************************************************
 *    ProjectPointOnLine                                        *
 * Computes Projection of a 3D point M                          *
 * on to a vector in 3D. The vector is directed                 *
 * from  P to S                                                 *
 * this can be useed in baseline computation                    *
 #%// KKM 09-March-2005                                         *
 ****************************************************************/
inline cn ProjectPointOnLine(
        cn M,
        cn P,
        cn S)
  {
  cn U            = P.min(S);// direction cosines of PS
  const real8 U_2 = U.in(U);
  U               = U.scale(1.0/sqrt(U_2));
  // --- Umat is the matrix representation of U to compute U*U' ---
  matrix<real4> Umat(3,1);
  Umat(0,0)       = U.x;
  Umat(1,0      ) = U.y;
  Umat(2,0)       = U.z;
  const matrix<real4> UU = matxmatT(Umat,Umat);// [3x3]
  // --- X is returned, base of perpendicular from M to PS
  matrix<real4> UmatT(1,3);
  UmatT(0,0)       = M.x-S.x;
  UmatT(0,1)       = M.y-S.y;
  UmatT(0,2)       = M.z-S.z;
  UmatT            = UmatT*UU;// [1x3]
  // --- X is returned, base of perpendicular from M to PS
  //cn X;
  //X.x              = UmatT(0,0)+S.x;
  //X.y              = UmatT(0,1)+S.y;
  //X.z              = UmatT(0,2)+S.z;
  //return X;
  return cn(UmatT(0,0)+S.x, UmatT(0,1)+S.y, UmatT(0,2)+S.z);
  }



/****************************************************************
 * matrix<real4> I = indgen(real4(10))                          *
 * matrix<int32> I = indgen(int32(10)), etc.                    *
 * create lying matrix [0,1,2,..N-1]                            *
 * Bert Kampes, 08-Apr-2005                                     *
 ****************************************************************/
template <class Type>
inline 
matrix<Type> indgen(
        const Type N) // use Type of N to create appropriate matrix
  {
  matrix<Type> e(1,N); 
  for(int32 i=0; i<N; ++i) e(0,i)=Type(i);
  return e;
  }
/****************************************************************
 * matrix<real4> I = linspace(x0,xn,N)                          *
 * create lying matrix [x0:xn] equally space N samples          *
 * Bert Kampes, 08-Apr-2005                                     *
 ****************************************************************/
template <class Type>
inline 
matrix<Type> linspace(
        const Type x0,
        const Type xN,
        const int32 N) 
  {
  return x0+indgen(Type(N))*((xN-x0)/Type(N-1));// use Type(N) to get ok.
  }
/****************************************************************
 * matrix<real4> E = ones(10,10)                                *
 * create matrix filled with ones.                              *
 * Bert Kampes, 08-Apr-2005                                     *
 ****************************************************************/
template <class Type>
inline
matrix<Type>  ones(
        const Type Ny,
        const Type Nx)// use Type of N to create appropriate matrix 
  {
  matrix<Type> E(Ny,Nx); 
  E.setdata(Type(1));
  return E;
  }
/****************************************************************
 * matrix<real4> I = eye(10)                                    *
 * create identity matrix with ones on main diagonal.           *
 * Bert Kampes, 08-Apr-2005                                     *
 ****************************************************************/
template <class Type>
inline
matrix<Type>  eye(
        const Type N)// use Type of N to create appropriate matrix
  {
  matrix<Type> I(N,N); 
  for(int32 i=0; i<N; ++i) I(i,i)=Type(1);
  return I;
  }




// ====== Prototypes ============================================
// ______ Call getorb for precise orbits ______
void getorb(
        const slcimage &image,
        const input_pr_orbits &inputorb,
        int16 FILEID);


// ______ Converts output to desired format ______
void convertgetorbout(
        int16 FILEID,
        const char* file="scratchorbit");

// ______ Return next power of 2 ______
uint nextpow2 (
        real8 w);


// ______ Evaluate 2d polynomial on grid ______
matrix<real4> polyval(
        const matrix<real4> &x,
        const matrix<real4> &y,
        const matrix<real8> &coeff,
        int32 degree = -1);


// ______ Evaluate 2d polynomial ______
real8 polyval(
        real8 x,
        real8 y,
        const matrix<real8> &coeff);


// ______ Evaluate 2d polynomial ______
real8 polyval(
        real8 x,
        real8 y,
        const matrix<real8> &coeff,
        int32 degree);


// ______ Evaluate 1d polynomial ______
real8 polyval1d(
        real8 x,
        const matrix<real8> &coeff);
//        int32 degree);


// ______ Normalize a matrix data -> [-2,2] ______
void normalize(
        matrix<real4> &data,
        real8 min,
        real8 max
        );


// ______ Normalize a matrix data -> [-2,2] ______
void normalize(
        matrix<real8> &data,
        real8 min,
        real8 max
        );


// ______ Normalize data from [min,max] to [-2,2] ______
inline real4 normalize(real4 data, real8 min, real8 max)
  {return real4((data-(.5*(min+max)))/(.25*(max-min)));}


// ______ Normalize data from [min,max] to [-2,2] ______
inline real8 normalize(real8 data, real8 min, real8 max)
  {return (data-(.5*(min+max)))/(.25*(max-min));}


// ______ Use this to compute Bperp for a few points, slow! ______
void BBparBperp(
        real8 &B, real8 &Bpar, real8 &Berp,     // returned
        const cn Master, const cn Point, const cn Slave);

// ______ Sign not computed, slow! ______
void BBhBv(
        real8 &B, real8 &Bh, real8 &Bv,         // returned
        const cn Master, const cn Slave);

// ______ Temporal baseline in days returned ______
// ______ (positive for master earlier than slave) ______
int32 Btemp(
        const char *utc_master, const char *utc_slave);

// ______ Baseline parametrizations returned ______
void BalphaBhBvBparBperpTheta(
        real8 &B, real8 &alpha, real8 &Bh, real8 &Bv,
        real8 &Bpar, real8 &Bperp, real8 &theta,
        const cn M, const cn P, const cn S);

// ______ Shift azimuth spectrum from fDC to zero, and vv. ______
void shiftazispectrum(
        matrix<complr4> &data,  // slcdata sin space domain
        const slcimage  &slave, // polynomial fDC
        const real4      shift);// abs(shift==firstpix in slave system)
                                // if shift>0 then shift from zero to fDC

#endif // UTILITIES_H

---