qgsderivativefilter.cpp

Go to the documentation of this file.

/***************************************************************************
                          qgsderivativefilter.cpp  -  description
                          -----------------------
    begin                : August 7th, 2009
    copyright            : (C) 2009 by Marco Hugentobler
    email                : marco dot hugentobler at karto dot baug dot ethz dot ch
 ***************************************************************************/
 
/***************************************************************************
 *                                                                         *
 *   This 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.                                   *
 *                                                                         *
 ***************************************************************************/
 
#include "qgsderivativefilter.h"
 
QgsDerivativeFilter::QgsDerivativeFilter( const QString &inputFile, const QString &outputFile, const QString &outputFormat )
  : QgsNineCellFilter( inputFile, outputFile, outputFormat )
{
 
}
 
float QgsDerivativeFilter::calcFirstDerX( float *x11, float *x21, float *x31, float *x12, float *x22, float *x32, float *x13, float *x23, float *x33 ) const
{
  //the basic formula would be simple, but we need to test for nodata values...
  //return (( (*x31 - *x11) + 2 * (*x32 - *x12) + (*x33 - *x13) ) / (8 * mCellSizeX));
 
  int weight = 0;
  double sum = 0;
 
  //first row
  if ( *x31 != mInputNodataValue && *x11 != mInputNodataValue ) //the normal case
  {
    sum += ( *x31 - *x11 );
    weight += 2;
  }
  else if ( *x31 == mInputNodataValue && *x11 != mInputNodataValue && *x21 != mInputNodataValue ) //probably 3x3 window is at the border
  {
    sum += ( *x21 - *x11 );
    weight += 1;
  }
  else if ( *x11 == mInputNodataValue && *x31 != mInputNodataValue && *x21 != mInputNodataValue ) //probably 3x3 window is at the border
  {
    sum += ( *x31 - *x21 );
    weight += 1;
  }
 
  //second row
  if ( *x32 != mInputNodataValue && *x12 != mInputNodataValue ) //the normal case
  {
    sum += 2 * ( *x32 - *x12 );
    weight += 4;
  }
  else if ( *x32 == mInputNodataValue && *x12 != mInputNodataValue && *x22 != mInputNodataValue )
  {
    sum += 2 * ( *x22 - *x12 );
    weight += 2;
  }
  else if ( *x12 == mInputNodataValue && *x32 != mInputNodataValue && *x22 != mInputNodataValue )
  {
    sum += 2 * ( *x32 - *x22 );
    weight += 2;
  }
 
  //third row
  if ( *x33 != mInputNodataValue && *x13 != mInputNodataValue ) //the normal case
  {
    sum += ( *x33 - *x13 );
    weight += 2;
  }
  else if ( *x33 == mInputNodataValue && *x13 != mInputNodataValue && *x23 != mInputNodataValue )
  {
    sum += ( *x23 - *x13 );
    weight += 1;
  }
  else if ( *x13 == mInputNodataValue && *x33 != mInputNodataValue && *x23 != mInputNodataValue )
  {
    sum += ( *x33 - *x23 );
    weight += 1;
  }
 
  if ( weight == 0 )
  {
    return mOutputNodataValue;
  }
 
  return sum / ( weight * mCellSizeX ) * mZFactor;
}
 
float QgsDerivativeFilter::calcFirstDerY( float *x11, float *x21, float *x31, float *x12, float *x22, float *x32, float *x13, float *x23, float *x33 ) const
{
  //the basic formula would be simple, but we need to test for nodata values...
  //return (((*x11 - *x13) + 2 * (*x21 - *x23) + (*x31 - *x33)) / ( 8 * mCellSizeY));
 
  double sum = 0;
  int weight = 0;
 
  //first row
  if ( *x11 != mInputNodataValue && *x13 != mInputNodataValue ) //normal case
  {
    sum += ( *x11 - *x13 );
    weight += 2;
  }
  else if ( *x11 == mInputNodataValue && *x13 != mInputNodataValue && *x12 != mInputNodataValue )
  {
    sum += ( *x12 - *x13 );
    weight += 1;
  }
  else if ( *x31 == mInputNodataValue && *x11 != mInputNodataValue && *x12 != mInputNodataValue )
  {
    sum += ( *x11 - *x12 );
    weight += 1;
  }
 
  //second row
  if ( *x21 != mInputNodataValue && *x23 != mInputNodataValue )
  {
    sum += 2 * ( *x21 - *x23 );
    weight += 4;
  }
  else if ( *x21 == mInputNodataValue && *x23 != mInputNodataValue && *x22 != mInputNodataValue )
  {
    sum += 2 * ( *x22 - *x23 );
    weight += 2;
  }
  else if ( *x23 == mInputNodataValue && *x21 != mInputNodataValue && *x22 != mInputNodataValue )
  {
    sum += 2 * ( *x21 - *x22 );
    weight += 2;
  }
 
  //third row
  if ( *x31 != mInputNodataValue && *x33 != mInputNodataValue )
  {
    sum += ( *x31 - *x33 );
    weight += 2;
  }
  else if ( *x31 == mInputNodataValue && *x33 != mInputNodataValue && *x32 != mInputNodataValue )
  {
    sum += ( *x32 - *x33 );
    weight += 1;
  }
  else if ( *x33 == mInputNodataValue && *x31 != mInputNodataValue && *x32 != mInputNodataValue )
  {
    sum += ( *x31 - *x32 );
    weight += 1;
  }
 
  if ( weight == 0 )
  {
    return mOutputNodataValue;
  }
 
  return sum / ( weight * mCellSizeY ) * mZFactor;
}