From: adv Date: Fri, 30 Aug 2013 07:33:16 +0000 (+0000) Subject: Calculation of the point's altitude from bathymetry (feature 9). X-Git-Tag: BR_hydro_v_0_1~91 X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=a14d55983f3a7f8f6573520d7a2c686eed8f54dd;p=modules%2Fhydro.git Calculation of the point's altitude from bathymetry (feature 9). --- diff --git a/src/HYDROData/HYDROData_Bathymetry.cxx b/src/HYDROData/HYDROData_Bathymetry.cxx index 31bd8cc7..afbb8d04 100644 --- a/src/HYDROData/HYDROData_Bathymetry.cxx +++ b/src/HYDROData/HYDROData_Bathymetry.cxx @@ -1,5 +1,6 @@ #include "HYDROData_Bathymetry.h" +#include "HYDROData_Tool.h" #include #include @@ -9,8 +10,12 @@ #include #include #include +#include #include +#define INVALID_ALTITUDE_VALUE -9999.0 + + IMPLEMENT_STANDARD_HANDLE(HYDROData_Bathymetry, HYDROData_Object) IMPLEMENT_STANDARD_RTTIEXT(HYDROData_Bathymetry, HYDROData_Object) @@ -77,6 +82,59 @@ void HYDROData_Bathymetry::RemoveAltitudePoints() aLab.ForgetAllAttributes(); } +void interpolateAltitudeForPoints( const gp_XY& thePoint, + const HYDROData_Bathymetry::AltitudePoint& theFirstPoint, + const HYDROData_Bathymetry::AltitudePoint& theSecPoint, + HYDROData_Bathymetry::AltitudePoint& theResPoint, + const bool& theIsVertical ) +{ + double aCoordX = thePoint.X(); + double aCoordY = thePoint.Y(); + + if ( theIsVertical ) + { + aCoordX = theFirstPoint.X(); + + if ( !ValuesEquals( theFirstPoint.X(), theSecPoint.X() ) ) + { + // Recalculate X coordinate by equation of line from two points + aCoordX = ( ( ( thePoint.Y() - theFirstPoint.Y() ) * ( theSecPoint.X() - theFirstPoint.X() ) ) / + ( theSecPoint.Y() - theFirstPoint.Y() ) ) + theFirstPoint.X(); + } + } + else + { + aCoordY = theFirstPoint.Y(); + + if ( !ValuesEquals( theFirstPoint.Y(), theSecPoint.Y() ) ) + { + // Recalculate y by equation of line from two points + aCoordY = ( ( ( thePoint.X() - theFirstPoint.X() ) * ( theSecPoint.Y() - theFirstPoint.Y() ) ) / + ( theSecPoint.X() - theFirstPoint.X() ) ) + theFirstPoint.Y(); + } + } + + theResPoint.SetX( aCoordX ); + theResPoint.SetY( aCoordY ); + + // Calculate coefficient for interpolation + double aLength = Sqrt( Pow( theSecPoint.Y() - theFirstPoint.Y(), 2 ) + + Pow( theSecPoint.X() - theFirstPoint.X(), 2 ) ); + + double aInterCoeff = 0; + if ( aLength != 0 ) + aInterCoeff = ( theSecPoint.Z() - theFirstPoint.Z() ) / aLength; + + + double aNewLength = Sqrt( Pow( theResPoint.Y() - theFirstPoint.Y(), 2 ) + + Pow( theResPoint.X() - theFirstPoint.X(), 2 ) ); + + // Calculate interpolated value + double aResVal = theFirstPoint.Z() + aInterCoeff * aNewLength; + + theResPoint.SetZ( aResVal ); +} + double HYDROData_Bathymetry::GetAltitudeForPoint( const QPointF& thePoint ) const { gp_XY aGpPoint( thePoint.x(), thePoint.y() ); @@ -88,8 +146,125 @@ double HYDROData_Bathymetry::GetAltitudeForPoint( const gp_XY& thePoint ) const double aResAltitude = -9999.90; AltitudePoints anAltitudePoints = GetAltitudePoints(); + if ( anAltitudePoints.isEmpty() ) + return aResAltitude; + + QPolygonF aBoundingRect; + + // Boundary plane + // [ 0 (top-left) ] [ 1 (top-right) ] + // thePoint + // [ 2 (bot-left) ] [ 3 (bot-right) ] + AltitudePoint aBounds[ 4 ] = { AltitudePoint( -DBL_MAX, -DBL_MAX, INVALID_ALTITUDE_VALUE ), + AltitudePoint( DBL_MAX, -DBL_MAX, INVALID_ALTITUDE_VALUE ), + AltitudePoint( -DBL_MAX, DBL_MAX, INVALID_ALTITUDE_VALUE ), + AltitudePoint( DBL_MAX, DBL_MAX, INVALID_ALTITUDE_VALUE ) }; + + AltitudePoints::const_iterator aListItBeg = anAltitudePoints.constBegin(); + AltitudePoints::const_iterator aListItEnd = anAltitudePoints.constEnd(); + for ( ; aListItBeg != aListItEnd; ++aListItBeg ) + { + const AltitudePoint& aPoint = *aListItBeg; + + double aDeltaX = Abs( aPoint.X() ) - Abs( thePoint.X() ); + double aDeltaY = Abs( aPoint.Y() ) - Abs( thePoint.Y() ); + + if ( ValuesEquals( aDeltaX, 0.0 ) ) // Both left and right sides + { + if ( ValuesEquals( aDeltaY, 0.0 ) ) // Both top and bottom sides + { + aResAltitude = aPoint.Z(); + return aResAltitude; + } + else if ( aDeltaY < 0 ) // top side + { + // top border + if ( ValuesMoreEquals( aPoint.X(), aBounds[ 0 ].X() ) && ValuesMoreEquals( aPoint.Y(), aBounds[ 0 ].Y() ) ) + aBounds[ 0 ] = aPoint; + if ( ValuesLessEquals( aPoint.X(), aBounds[ 1 ].X() ) && ValuesMoreEquals( aPoint.Y(), aBounds[ 1 ].Y() ) ) + aBounds[ 1 ] = aPoint; + } + else + { + // bottom border + if ( ValuesMoreEquals( aPoint.X(), aBounds[ 2 ].X() ) && ValuesLessEquals( aPoint.Y(), aBounds[ 2 ].Y() ) ) + aBounds[ 2 ] = aPoint; + if ( ValuesLessEquals( aPoint.X(), aBounds[ 3 ].X() ) && ValuesLessEquals( aPoint.Y(), aBounds[ 3 ].Y() ) ) + aBounds[ 3 ] = aPoint; + } + } + else if ( aDeltaX < 0 ) // left side + { + if ( ValuesEquals( aDeltaY, 0.0 ) ) + { + // Left border + if ( ValuesMoreEquals( aPoint.X(), aBounds[ 0 ].X() ) && ValuesMoreEquals( aPoint.Y(), aBounds[ 0 ].Y() ) ) + aBounds[ 0 ] = aPoint; + if ( ValuesMoreEquals( aPoint.X(), aBounds[ 2 ].X() ) && ValuesLessEquals( aPoint.Y(), aBounds[ 2 ].Y() ) ) + aBounds[ 2 ] = aPoint; + } + else if ( aDeltaY < 0 ) + { + // top left corner + if ( ValuesMoreEquals( aPoint.X(), aBounds[ 0 ].X() ) && ValuesMoreEquals( aPoint.Y(), aBounds[ 0 ].Y() ) ) + aBounds[ 0 ] = aPoint; + } + else + { + // bottom left corner + if ( ValuesMoreEquals( aPoint.X(), aBounds[ 2 ].X() ) && ValuesLessEquals( aPoint.Y(), aBounds[ 2 ].Y() ) ) + aBounds[ 2 ] = aPoint; + } + } + else // right side + { + if ( ValuesEquals( aDeltaY, 0.0 ) ) + { + // Right border + if ( ValuesLessEquals( aPoint.X(), aBounds[ 1 ].X() ) && ValuesMoreEquals( aPoint.Y(), aBounds[ 1 ].Y() ) ) + aBounds[ 1 ] = aPoint; + if ( ValuesLessEquals( aPoint.X(), aBounds[ 3 ].X() ) && ValuesLessEquals( aPoint.Y(), aBounds[ 3 ].Y() ) ) + aBounds[ 3 ] = aPoint; + } + else if ( aDeltaY < 0 ) + { + // top right corner + if ( ValuesLessEquals( aPoint.X(), aBounds[ 1 ].X() ) && ValuesMoreEquals( aPoint.Y(), aBounds[ 1 ].Y() ) ) + aBounds[ 1 ] = aPoint; + } + else + { + // bottom right corner + if ( ValuesLessEquals( aPoint.X(), aBounds[ 3 ].X() ) && ValuesLessEquals( aPoint.Y(), aBounds[ 3 ].Y() ) ) + aBounds[ 3 ] = aPoint; + } + } + + // Update bounding rectangle of our global grid + aBoundingRect << QPointF( aPoint.X(), aPoint.Y() ); + } - // TODO : implement + // Check if requested point is inside of our bounding rectangle + if ( !aBoundingRect.boundingRect().contains( thePoint.X(), thePoint.Y() ) ) + return INVALID_ALTITUDE_VALUE; + + // Calculate result altitude for point + AltitudePoint aFirstPoint( aBounds[ 0 ] ), aSecPoint( aBounds[ 1 ] ); + + // At first we merge top and bottom borders + if ( aBounds[ 0 ].Y() != aBounds[ 2 ].Y() || aBounds[ 0 ].X() != aBounds[ 2 ].X() ) + interpolateAltitudeForPoints( thePoint, aBounds[ 0 ], aBounds[ 2 ], aFirstPoint, true ); + + if ( aBounds[ 1 ].Y() != aBounds[ 3 ].Y() || aBounds[ 1 ].X() != aBounds[ 3 ].X() ) + interpolateAltitudeForPoints( thePoint, aBounds[ 1 ], aBounds[ 3 ], aSecPoint, true ); + + AltitudePoint aResPoint( aFirstPoint ); + + // At last we merge left and right borders + if ( aFirstPoint.Y() != aSecPoint.Y() || aFirstPoint.X() != aSecPoint.X() ) + interpolateAltitudeForPoints( thePoint, aFirstPoint, aSecPoint, aResPoint, false ); + + aResAltitude = aResPoint.Z(); return aResAltitude; } @@ -136,7 +311,7 @@ bool HYDROData_Bathymetry::importFromXYZFile( QFile& theFile, while ( !theFile.atEnd() ) { - QString aLine = theFile.readLine(); + QString aLine = theFile.readLine().simplified(); if ( aLine.isEmpty() ) continue;