Logo Search packages:      
Sourcecode: saga version File versions  Download package

Morphometry.cpp


///////////////////////////////////////////////////////////
//                                                       //
//                         SAGA                          //
//                                                       //
//      System for Automated Geoscientific Analyses      //
//                                                       //
//                    Module Library:                    //
//                    ta_morphometry                     //
//                                                       //
//-------------------------------------------------------//
//                                                       //
//                    Morphometry.cpp                    //
//                                                       //
//                 Copyright (C) 2003 by                 //
//                      Olaf Conrad                      //
//                                                       //
//-------------------------------------------------------//
//                                                       //
// This file is part of 'SAGA - System for Automated     //
// Geoscientific Analyses'. SAGA 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; version 2 of the License.   //
//                                                       //
// SAGA 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.                                                  //
//                                                       //
//-------------------------------------------------------//
//                                                       //
//    e-mail:     oconrad@saga-gis.org                   //
//                                                       //
//    contact:    Olaf Conrad                            //
//                Institute of Geography                 //
//                University of Goettingen               //
//                Goldschmidtstr. 5                      //
//                37077 Goettingen                       //
//                Germany                                //
//                                                       //
///////////////////////////////////////////////////////////

//---------------------------------------------------------


///////////////////////////////////////////////////////////
//                                                                                   //
//                                                                                   //
//                                                                                   //
///////////////////////////////////////////////////////////

//---------------------------------------------------------
#include "Morphometry.h"


///////////////////////////////////////////////////////////
//                                                                                   //
//                                                                                   //
//                                                                                   //
///////////////////////////////////////////////////////////

//---------------------------------------------------------
CMorphometry::CMorphometry(void)
{
      Set_Name          (_TL("Local Morphometry"));

      Set_Author        (SG_T("(c) 2001 by O.Conrad"));

      Set_Description   (_TW(
            "Calculates local morphometric terrain attributes (i.e. slope, aspect and curvatures).\n\n"
            "References:\n\n"

            "Maximum Slope\n"
            "- Travis, M.R. / Elsner, G.H. / Iverson, W.D. / Johnson, C.G. (1975):\n"
            "    'VIEWIT: computation of seen areas, slope, and aspect for land-use planning',\n"
            "    USDA F.S. Gen. Tech. Rep. PSW-11/1975, 70p. Berkeley, California, U.S.A.\n\n"

            "Maximum Triangle Slope\n"
            "- Tarboton, D.G. (1997):\n"
            "    'A new method for the determination of flow directions and upslope areas in grid digital elevation models',\n"
            "    Water Ressources Research, Vol.33, No.2, p.309-319\n\n"

            "Least Squares or Best Fitted Plane\n"
            "- Horn, B. K. (1981):\n"
            "    'Hill shading and the relectance map',\n"
            "    Proceedings of the IEEE, v. 69, no. 1, p. 14-47.\n\n"

            "- Beasley, D.B. / Huggins, L.F. (1982):\n"
            "    'ANSWERS: Userís manual',\n"
            "    U.S. EPA-905/9-82-001, Chicago, IL. 54pp.\n\n"

            "- Costa-Cabral, M., and Burges, S.J., (1994):\n"
            "    'Digital Elevation Model Networks (DEMON): a model of flow over hillslopes for computation of contributing and dispersal areas',\n"
            "    Water Resources Research, v. 30, no. 6, p. 1681-1692.\n\n"

            "Fit 2.Degree Polynom\n"
            "- Bauer, J. / Rohdenburg, H. / Bork, H.-R. (1985):\n"
            "    'Ein Digitales Reliefmodell als Vorraussetzung fuer ein deterministisches Modell der Wasser- und Stoff-Fluesse',\n"
            "    Landschaftsgenese und Landschaftsoekologie, H.10, Parameteraufbereitung fuer deterministische Gebiets-Wassermodelle,\n"
            "    Grundlagenarbeiten zu Analyse von Agrar-Oekosystemen, (Eds.: Bork, H.-R. / Rohdenburg, H.), p.1-15\n\n"

            "- Heerdegen, R.G. / Beran, M.A. (1982):\n"
            "    'Quantifying source areas through land surface curvature',\n"
            "    Journal of Hydrology, Vol.57\n\n"

            "- Zevenbergen, L.W. / Thorne, C.R. (1987):\n"
            "    'Quantitative analysis of land surface topography',\n"
            "    Earth Surface Processes and Landforms, 12: 47-56.\n\n"

            "Fit 3.Degree Polynom\n"
            "- R.M. Haralick (1983):\n"
            "    'Ridge and valley detection on digital images',\n"
            "    Computer Vision, Graphics and Image Processing, Vol.22, No.1, p.28-38\n\n"
      ));


      //-----------------------------------------------------
      // Input...

      Parameters.Add_Grid(
            NULL  , "ELEVATION"     , _TL("Elevation"),
            _TL(""),
            PARAMETER_INPUT
      );


      //-----------------------------------------------------
      // Output...

      Parameters.Add_Grid(
            NULL  , "SLOPE"         , _TL("Slope"),
            _TL(""),
            PARAMETER_OUTPUT
      );

      Parameters.Add_Grid(
            NULL  , "ASPECT"        , _TL("Aspect"),
            _TL(""),
            PARAMETER_OUTPUT
      );

      Parameters.Add_Grid(
            NULL  , "CURV"          , _TL("Curvature"),
            _TL(""),
            PARAMETER_OUTPUT_OPTIONAL
      );

      Parameters.Add_Grid(
            NULL  , "HCURV"         , _TL("Plan Curvature"),
            _TL(""),
            PARAMETER_OUTPUT_OPTIONAL
      );

      Parameters.Add_Grid(
            NULL  , "VCURV"         , _TL("Profile Curvature"),
            _TL(""),
            PARAMETER_OUTPUT_OPTIONAL
      );


      //-----------------------------------------------------
      // Options...

      Parameters.Add_Choice(
            NULL  , "METHOD"        , _TL("Method"),
            _TL(""),
            CSG_String::Format(SG_T("%s|%s|%s|%s|%s|%s|%s|"),
                  _TL("Maximum Slope (Travis et al. 1975)"),
                  _TL("Maximum Triangle Slope (Tarboton 1997)"),
                  _TL("Least Squares Fitted Plane (Horn 1981, Costa-Cabral & Burgess 1996)"),
                  _TL("Fit 2.Degree Polynom (Bauer, Rohdenburg, Bork 1985)"),
                  _TL("Fit 2.Degree Polynom (Heerdegen & Beran 1982)"),
                  _TL("Fit 2.Degree Polynom (Zevenbergen & Thorne 1987)"),
                  _TL("Fit 3.Degree Polynom (Haralick 1983)")
            ), 5
      );
}

//---------------------------------------------------------
CMorphometry::~CMorphometry(void)
{}


///////////////////////////////////////////////////////////
//                                                                                   //
//                                                                                   //
//                                                                                   //
///////////////////////////////////////////////////////////

//---------------------------------------------------------
bool CMorphometry::On_Execute(void)
{
      int               x, y, Method;
      CSG_Colors  Colors;

      //-----------------------------------------------------
      Method            = Parameters("METHOD"   )->asInt();

      pDTM        = Parameters("ELEVATION")->asGrid();
      pSlope            = Parameters("SLOPE"    )->asGrid();
      pAspect           = Parameters("ASPECT"   )->asGrid();
      pCurvature  = Parameters("CURV"     )->asGrid();
      pCurv_Horz  = Parameters("HCURV"    )->asGrid();
      pCurv_Vert  = Parameters("VCURV"    )->asGrid();
      pCurv_Tang  = NULL;//Parameters("TCURV"    )->asGrid();

      //-----------------------------------------------------
      pSlope->Set_ZFactor     (M_RAD_TO_DEG);
      pSlope->Set_Unit  (_TL("Degree"));
      DataObject_Set_Colors(pSlope  , 100, SG_COLORS_YELLOW_RED);

      pAspect->Set_ZFactor(M_RAD_TO_DEG);
      pAspect->Set_Unit (_TL("Degree"));
      Colors.Set_Count(3);
      Colors.Set_Color(0,     SG_GET_RGB(255, 255, 190));
      Colors.Set_Color(1, SG_GET_RGB( 64,   0,   0));
      Colors.Set_Color(2, SG_GET_RGB(255, 255, 190));
      Colors.Set_Count(100);
      DataObject_Set_Colors(pAspect , Colors);

      DataObject_Set_Colors(pCurvature, 100, SG_COLORS_RED_GREY_BLUE, true);
      DataObject_Set_Colors(pCurv_Vert, 100, SG_COLORS_RED_GREY_BLUE, true);
      DataObject_Set_Colors(pCurv_Horz, 100, SG_COLORS_RED_GREY_BLUE, true);

      //-----------------------------------------------------
      _DX_2       = Get_Cellsize() * Get_Cellsize();
      _4DX_2            = 4.0 * _DX_2;
      _6DX        = 6.0 * Get_Cellsize();
      _2DX        = 2.0 * Get_Cellsize();

      //-----------------------------------------------------
      for(y=0; y<Get_NY() && Set_Progress(y); y++)
      {
            for(x=0; x<Get_NX(); x++)
            {
                  switch( Method )
                  {
                  case 0:
                        Do_MaximumSlope(  x, y );
                        break;

                  case 1:
                        Do_Tarboton(            x, y );
                        break;

                  case 2:
                        Do_LeastSquare(         x, y );
                        break;

                  case 3:
                        Do_FD_BRM(              x, y );
                        break;

                  case 4:
                        Do_FD_Heerdegen(  x, y );
                        break;

                  case 5:
                        Do_FD_Zevenbergen(      x, y );
                        break;

                  case 6:
                        Do_FD_Haralick(         x, y );
                        break;
                  }
            }
      }

      return( true );
}


///////////////////////////////////////////////////////////
//                                                                                   //
//                                                                                   //
//                                                                                   //
///////////////////////////////////////////////////////////

//---------------------------------------------------------
inline void CMorphometry::Set_Parameters(int x, int y, double Slope, double Aspect, double Curv, double vCurv, double hCurv, double tCurv)
{
      if( pSlope )
      {
            pSlope->Set_Value(            x, y, Slope );
      }

      if( pAspect )
      {
            pAspect->Set_Value(           x, y, Aspect );
      }

      if( pCurvature )
      {
            pCurvature->Set_Value(  x, y,  Curv );
      }

      if( pCurv_Vert )
      {
            pCurv_Vert->Set_Value(  x, y, vCurv );
      }

      if( pCurv_Horz )
      {
            pCurv_Horz->Set_Value(  x, y, hCurv );
      }

      if( pCurv_Tang )
      {
            pCurv_Tang->Set_Value(  x, y, tCurv );
      }
}

//---------------------------------------------------------
inline void CMorphometry::Set_Parameters_Derive(int x, int y, double D, double E, double F, double G, double H)
{
      double      k1, k2, Curv, vCurv, hCurv, tCurv;

      k1          = F * G * H;
      k2          = G*G + H*H;

      if( k2 )
      {
            Curv  = -2.0 * (E + D);
            vCurv = -2.0 * (D * G*G + E * H*H + k1) /  k2;
            hCurv = -2.0 * (D * H*H + E * G*G - k1) /  k2;
            tCurv = -2.0 * (D * H*H + E * G*G - k1) / (k2 * sqrt(k2 + 1.0));

//          Curv  = -2.0 * (E + D + F*F);
//          vCurv = -2.0 * (D * G*G + E * H*H + k1) / (k2 *  pow(k2 + 1.0, 3.0 / 2.0));
//          hCurv = -2.0 * (D * H*H + E * G*G - k1) / (      pow(k2      , 3.0 / 2.0));
//          tCurv = -2.0 * (D * H*H + E * G*G - k1) / (k2 * sqrt(k2 + 1.0)           );
      }
      else
      {
            Curv  = vCurv     = hCurv     = tCurv     = 0.0;
      }

      if( G != 0.0 )
      {
            Set_Parameters(x, y, atan(sqrt(k2)), M_PI_180 + atan2(H, G) , Curv, vCurv, hCurv, tCurv);
      }
      else if( H > 0.0 )
      {
            Set_Parameters(x, y, atan(sqrt(k2)), M_PI_270                     , Curv, vCurv, hCurv, tCurv);
      }
      else if( H < 0.0 )
      {
            Set_Parameters(x, y, atan(sqrt(k2)), M_PI_090                     , Curv, vCurv, hCurv, tCurv);
      }
      else
      {
            Set_Parameters(x, y, atan(sqrt(k2)), 0.0                          , Curv, vCurv, hCurv, tCurv);

            if( pAspect )
            {
                  pAspect->Set_NoData(x, y);
            }
      }
}

//---------------------------------------------------------
inline void CMorphometry::Set_Parameters_NoData(int x, int y, bool bCompletely)
{
      if( bCompletely )
      {
            if( pSlope )
            {
                  pSlope->Set_NoData(           x, y );
            }

            if( pAspect )
            {
                  pAspect->Set_NoData(    x, y );
            }

            if( pCurvature )
            {
                  pCurvature->Set_NoData( x, y );
            }

            if( pCurv_Vert )
            {
                  pCurv_Vert->Set_NoData( x, y );
            }

            if( pCurv_Horz )
            {
                  pCurv_Horz->Set_NoData( x, y );
            }

            if( pCurv_Tang )
            {
                  pCurv_Tang->Set_NoData( x, y );
            }
      }
      else
      {
            Set_Parameters(x, y, 0.0, 0.0, 0.0, 0.0, 0.0);

            if( pAspect )
            {
                  pAspect->Set_NoData(    x, y );
            }
      }
}


///////////////////////////////////////////////////////////
//                                                                                   //
//                                                                                   //
//                                                                                   //
///////////////////////////////////////////////////////////

//---------------------------------------------------------
// Indexing of the Submatrix:
//
//  +-------+    +-------+
//  | 7 0 1 |    | 2 5 8 |
//  | 6 * 2 | => | 1 4 7 |
//  | 5 4 3 |    | 0 3 6 |
//  +-------+    +-------+
//
//---------------------------------------------------------
inline bool CMorphometry::Get_SubMatrix3x3(int x, int y, double SubMatrix[9])
{
      static int  iSub[]      = { 5, 8, 7, 6, 3, 0, 1, 2 };

      int         i, ix, iy;
      double      z;

      if( pDTM->is_NoData(x, y) )
      {
            Set_Parameters_NoData(x, y, true);
      }
      else
      {
            z                       = pDTM->asDouble(x, y);
            SubMatrix[4]      = 0.0;

            for(i=0; i<8; i++)
            {
                  ix    = Get_xTo(i, x);
                  iy    = Get_yTo(i, y);

                  if( pDTM->is_InGrid(ix, iy) )
                  {
                        SubMatrix[iSub[i]]      = pDTM->asDouble(ix, iy) - z;
                  }
                  else
                  {
                        ix    = Get_xTo(i + 4, x);
                        iy    = Get_yTo(i + 4, y);

                        if( pDTM->is_InGrid(ix, iy) )
                        {
                              SubMatrix[iSub[i]]      = z - pDTM->asDouble(ix, iy);
                        }
                        else
                        {
                              SubMatrix[iSub[i]]      = 0.0;
                        }
                  }
            }

            return( true );
      }

      return( false );
}

//---------------------------------------------------------
inline bool CMorphometry::Get_SubMatrix5x5(int x, int y, double SubMatrix[25])
{
      int         i, ix, iy, jx, jy;
      double      z;

      if( pDTM->is_NoData(x, y) )
      {
            Set_Parameters_NoData(x, y, true);
      }
      else
      {
            z     = pDTM->asDouble(x,y);

            for(i=0, iy=y-2; iy<=y+2; iy++)
            {
                  jy    = iy < 0 ? 0 : ( iy >= Get_NY() ? Get_NY() - 1 : iy );

                  for(ix=x-2; ix<=x+2; ix++, i++)
                  {
                        jx    = ix < 0 ? 0 : ( ix >= Get_NX() ? Get_NX() - 1 : ix );

                        SubMatrix[i]      = pDTM->is_InGrid(jx, jy) ? pDTM->asDouble(jx, jy) - z : 0.0;
                  }
            }

            return( true );
      }

      return( false );
}


///////////////////////////////////////////////////////////
//                                                                                   //
//                            The Methods                                      //
//                                                                                   //
///////////////////////////////////////////////////////////

//---------------------------------------------------------
// Maximum Slope (Travis et al., 1975, Peucker & Douglas, 1975))
//
// Travis, M.R., Elsner, G.H., Iverson, W.D., and Johnson, C.G. 1975:
//          VIEWIT: computation of seen areas, slope, and aspect for land-use planning.
//          USDA F.S. Gen. Tech. Rep. PSW-11/1975, 70p. Berkeley, California, U.S.A.
//
//---------------------------------------------------------
void CMorphometry::Do_MaximumSlope(int x, int y)
{
      int         i, ix, iy, j, Aspect;

      double      z, zm[8], Slope, Curv, hCurv, a, b;

      //-----------------------------------------------------
      if( pDTM->is_NoData(x, y) )
      {
            Set_Parameters_NoData(x, y, true);
      }
      else
      {
            //-------------------------------------------------
            z           = pDTM->asDouble(x, y);
        Slope     = Curv      = 0.0;

            for(Aspect=-1, i=0; i<8; i++)
            {
                  ix          = Get_xTo(i, x);
                  iy          = Get_yTo(i, y);

                  if( !pDTM->is_InGrid(ix, iy) )
                  {
                        zm[i] = 0.0;
                  }
                  else
                  {
                        zm[i] = atan((z - pDTM->asDouble(ix, iy)) / Get_Length(i));
                        Curv  += zm[i];

                        if( zm[i] > Slope )
                        {
                              Aspect      = i;
                              Slope = zm[i];
                        }
                  }
            }

            //-------------------------------------------------
            if( Aspect < 0.0 )
            {
                  Set_Parameters_NoData(x, y);
            }
            else
            {
                  //---------------------------------------------
                  // Let's now estimate the plan curvature...

                  for(i=Aspect+1, j=0, a=0.0; i<Aspect+8; i++, j++)
                  {
                        if( zm[i % 8] < 0.0 )
                        {
                              a     = j + zm[(i - 1) % 8] / (zm[(i - 1) % 8] - zm[i % 8]);
                              break;
                        }
                  }

                  if( a != 0.0 )
                  {
                        for(i=Aspect+7, j=0, b=0.0; i>Aspect; i--, j++)
                        {
                              if( zm[i % 8] < 0.0 )
                              {
                                    b     = j + zm[(i + 1) % 8] / (zm[(i + 1) % 8] - zm[i % 8]);
                                    break;
                              }
                        }

                        hCurv =  45.0 * (a + b) - 180.0;
                  }
                  else
                  {
                        hCurv =  180.0;
                  }

                  //---------------------------------------------
                  Set_Parameters(x, y,
                        Slope,
                        Aspect * M_PI_045,
                        Curv,
                        zm[Aspect] + zm[(Aspect + 4) % 8],
                        hCurv
                  );
            }
    }
}

//---------------------------------------------------------
// Maximum Triangle Slope
//
// Tarboton, D.G. (1997):
//          'A new method for the determination of flow directions and upslope areas in grid digital elevation models',
//          Water Ressources Research, Vol.33, No.2, p.309-319
//
//---------------------------------------------------------
void CMorphometry::Do_Tarboton(int x, int y)
{
      int         i, ix, iy, j;
      double      z, zm[8], iSlope, iAspect, Slope, Aspect, G, H;

      //-----------------------------------------------------
      if( pDTM->is_NoData(x, y) )
      {
            Set_Parameters_NoData(x, y, true);
      }
      else
      {
            z           = pDTM->asDouble(x, y);

            for(i=0; i<8; i++)
            {
                  ix          = Get_xTo(i, x);
                  iy          = Get_yTo(i, y);

                  if( pDTM->is_InGrid(ix, iy) )
                  {
                        zm[i] =  pDTM->asDouble(ix, iy);
                  }
                  else
                  {
                        ix          = Get_xTo(i + 4, x);
                        iy          = Get_yTo(i + 4, y);

                        if( pDTM->is_InGrid(ix, iy) )
                        {
                              zm[i] =  z - (pDTM->asDouble(ix, iy) - z);
                        }
                        else
                        {
                              zm[i] =  z;
                        }
                  }
            }

            //---------------------------------------------
        Slope     =  0.0;
            Aspect      = -1.0;

            for(i=0, j=1; i<8; i++, j=(j+1)%8)
            {
                  if( i % 2 ) // i => diagonal
                  {
                        G           = (z        - zm[j]) / Get_Cellsize();
                        H           = (zm[j]    - zm[i]) / Get_Cellsize();
                  }
                  else        // i => orthogonal
                  {
                        G           = (z        - zm[i]) / Get_Cellsize();
                        H           = (zm[i]    - zm[j]) / Get_Cellsize();
                  }

                  if( H < 0.0 )
                  {
                        iAspect     = 0.0;
                        iSlope      = G;
                  }
                  else if( H > G )
                  {
                        iAspect     = M_PI_045;
                        iSlope      = (z - zm[i % 2 ? i : j]) / (sqrt(2.0) * Get_Cellsize());
                  }
                  else
                  {
                        iAspect     = atan(H / G);
                        iSlope      = sqrt(G*G + H*H);
                  }

                  if( iSlope > Slope )
                  {
                        Aspect      = i * M_PI_045 + (i % 2 ? M_PI_045 - iAspect : iAspect);
                        Slope = iSlope;
                  }
            }

            //---------------------------------------------
            if( Aspect < 0.0 )
            {
                  Set_Parameters_NoData(x, y);
            }
            else
            {
                  Set_Parameters(x, y, atan(Slope), Aspect);
            }
      }
}

//---------------------------------------------------------
// Least Squares or Best Fit Plane (Horn 1981, Beasley & Huggins 1982, Costa-Cabral & Burgess 1994)
//
// Horn, B. K. (1981):
//      Hill shading and the relectance map.
//      Proceedings of the IEEE, v. 69, no. 1, p 14-47.
//
// Beasley, D.B. and Huggins, L.F. 1982:
//          ANSWERS: Userís manual.
//          U.S. EPA-905/9-82-001, Chicago, IL. 54pp.
//
// Costa-Cabral, M., and Burges, S.J., 1994:
//          Digital Elevation Model Networks (DEMON): a model of flow over hillslopes for computation of contributing and dispersal areas
//          Water Resources Research, v. 30, no. 6, p. 1681-1692.
//
//---------------------------------------------------------
void CMorphometry::Do_LeastSquare(int x, int y)
{
      double      zm[9], a, b;

      if( Get_SubMatrix3x3(x, y, zm) )
      {
            a = ((zm[2] + 2 * zm[5] + zm[8]) - (zm[0] + 2 * zm[3] + zm[6])) / (8 * Get_Cellsize());
            b = ((zm[6] + 2 * zm[7] + zm[8]) - (zm[0] + 2 * zm[1] + zm[2])) / (8 * Get_Cellsize());

            if( a != 0.0 )
            {
                  Set_Parameters(x, y, atan( sqrt(a*a + b*b) ), M_PI_180 + atan2(b, a));
            }
            else if( b > 0.0 )
            {
                  Set_Parameters(x, y, atan( sqrt(a*a + b*b) ), M_PI_270);
            }
            else if( b < 0.0 )
            {
                  Set_Parameters(x, y, atan( sqrt(a*a + b*b) ), M_PI_090);
            }
            else
            {
                  Set_Parameters_NoData(x, y);
            }
      }
}

//---------------------------------------------------------
// Quadratic Function Approximation (Bauer, Rohdenburg & Bork, 1985)
//
// Bauer, J. / Rohdenburg, H. / Bork, H.-R., (1985):
//          'Ein Digitales Reliefmodell als Vorraussetzung fuer ein deterministisches Modell der Wasser- und Stoff-Fluesse',
//          Landschaftsgenese und Landschaftsoekologie, H.10, Parameteraufbereitung fuer deterministische Gebiets-Wassermodelle,
//          Grundlagenarbeiten zu Analyse von Agrar-Oekosystemen, (Eds.: Bork, H.-R. / Rohdenburg, H.), p.1-15
//
//---------------------------------------------------------
void CMorphometry::Do_FD_BRM(int x, int y)
{
      double      zm[9], D, E, F, G, H;

      if( Get_SubMatrix3x3(x, y, zm) )
      {
            D     = ( (zm[0] + zm[2] + zm[3] + zm[5] + zm[6] + zm[8]) - 2 * (zm[1] + zm[4] + zm[7]) ) / _DX_2;
            E     = ( (zm[0] + zm[6] + zm[1] + zm[7] + zm[2] + zm[8]) - 2 * (zm[3] + zm[4] + zm[5]) ) / _DX_2;
        F   = (  zm[8] + zm[0] - zm[7] )                                                                                / _4DX_2;
            G     = ( (zm[2] - zm[0]) + (zm[5] - zm[3]) + (zm[8]-zm[6]) )                                               / _6DX;
            H     = ( (zm[6] - zm[0]) + (zm[7] - zm[1]) + (zm[8]-zm[2]) )                                               / _6DX;

            Set_Parameters_Derive(x, y, D, E, F, G, H);
      }
}

//---------------------------------------------------------
// Quadratic Function Approximation (Heerdegen & Beran, 1984)
//
// Heerdegen, R.G. / Beran, M.A. (1982):
//          'Quantifying source areas through land surface curvature',
//          Journal of Hydrology, Vol.57
//
//---------------------------------------------------------
void CMorphometry::Do_FD_Heerdegen(int x, int y)
{
      double      zm[9], D, E, F, G, H, a, b;

      //-----------------------------------------------------
      if( Get_SubMatrix3x3(x, y, zm) )
      {
            a     =   zm[0] + zm[2] + zm[3] + zm[5] + zm[6] + zm[8];
            b     =   zm[0] + zm[1] + zm[2] + zm[6] + zm[7] + zm[8];
            D     = (0.3 * a - 0.2 * b)                                             / _DX_2;
            E     = (0.3 * b - 0.2 * a)                                             / _DX_2;
            F     = ( zm[0] - zm[2]                 - zm[6] + zm[8])    / _4DX_2;
            G     = (-zm[0] + zm[2] - zm[3] + zm[5] - zm[6] + zm[8])    / _6DX;
        H   = (-zm[0] - zm[1] - zm[2] + zm[6] + zm[7] + zm[8])    / _6DX;

            Set_Parameters_Derive(x, y, D, E, F, G, H);
      }
}

//---------------------------------------------------------
// Quadratic Function Approximation (Zevenbergen und Thorne, 1986)
//
// Zevenbergen, L.W. and C.R. Thorne. 1987:
//          Quantitative analysis of land surface topography
//          Earth Surface Processes and Landforms, 12: 47-56.
//
//---------------------------------------------------------
void CMorphometry::Do_FD_Zevenbergen(int x, int y)
{
      double      zm[9], D, E, F, G, H;

      //-----------------------------------------------------
      if( Get_SubMatrix3x3(x, y, zm) )
      {
            D     = ((zm[3] + zm[5]) / 2.0 - zm[4])   / _DX_2;
            E     = ((zm[1] + zm[7]) / 2.0 - zm[4])   / _DX_2;
            F     =  (zm[0] - zm[2] - zm[6] + zm[8])  / _4DX_2;
            G     =  (zm[5] - zm[3])                              / _2DX;
        H   =  (zm[7] - zm[1])                              / _2DX;

            Set_Parameters_Derive(x, y, D, E, F, G, H);
      }
}

//---------------------------------------------------------
// Cubic Function Approximation (Haralick, 1991)
//
// R.M. Haralick (1983):
//          'Ridge and Valley Detection on digital images',
//          Computer Vision, Graphics and Image Processing, Vol.22, No.1, p.28-38
//
//---------------------------------------------------------
void CMorphometry::Do_FD_Haralick(int x, int y)
{
      //-----------------------------------------------------
      // Finite Differenzen Methode Matrizen...

      const int   Mtrx[][5][5]      = {
      {     { 31,- 5,-17,- 5, 31}, {-44,-62,-68,-62,-44}, {  0,  0,  0,  0,  0}, { 44, 62, 68, 62, 44}, {-31,  5, 17,  5,-31} },
      {     { 31,-44,  0, 44,-31}, {- 5,-62,  0, 62,  5}, {-17,-68,  0, 68, 17}, {- 5,-62,  0, 62,  5}, { 31,-44,  0, 44,-31} },
      {     {  2,  2,  2,  2,  2}, {- 1,- 1,- 1,- 1,- 1}, {- 2,- 2,- 2,- 2,- 2}, {- 1,- 1,- 1,- 1,- 1}, {  2,  2,  2,  2,  2} },
      {     {  4,  2,  0,- 2,- 4}, {  2,  1,  0,- 1,- 2}, {  0,  0,  0,  0,  0}, {- 2,- 1,  0,  1,  2}, {- 4,- 2,  0,  2,  4} },
      {     {  2,- 1,- 2,- 1,  2}, {  2,- 1,- 2,- 1,  2}, {  2,- 1,- 2,- 1,  2}, {  2,- 1,- 2,- 1,  2}, {  2,- 1,- 2,- 1,  2} },    };

      const int   QMtrx[]                 = { 4200, 4200, 700, 1000, 700 };

      //-----------------------------------------------------
      int         i, ix, iy, n;
      double      Sum, zm[25], k[5];

      //-----------------------------------------------------
      if( Get_SubMatrix5x5(x, y, zm) )
      {
            for(i=0; i<5; i++)
            {
                  for(n=0, Sum=0.0, iy=0; iy<5; iy++)
                  {
                        for(ix=0; ix<5; ix++, n++)
                        {
                              Sum   += zm[n] * Mtrx[i][ix][iy];
                        }
                  }

                  k[i]  = Sum / QMtrx[i];
            }

            Set_Parameters_Derive(x, y, k[4], k[2], k[3], k[1], k[0]);
    }
}

Generated by  Doxygen 1.6.0   Back to index