diff --git a/ApplicationCode/Application/Tools/CMakeLists_files.cmake b/ApplicationCode/Application/Tools/CMakeLists_files.cmake index eea4fd4aeb..f5f755e976 100644 --- a/ApplicationCode/Application/Tools/CMakeLists_files.cmake +++ b/ApplicationCode/Application/Tools/CMakeLists_files.cmake @@ -28,6 +28,7 @@ ${CMAKE_CURRENT_LIST_DIR}/RiaTimeHistoryCurveMerger.h ${CMAKE_CURRENT_LIST_DIR}/RiaCurveDataTools.h ${CMAKE_CURRENT_LIST_DIR}/RiaTimeHistoryCurveResampler.h ${CMAKE_CURRENT_LIST_DIR}/RiaStatisticsTools.h +${CMAKE_CURRENT_LIST_DIR}/RiaPolyArcLineSampler.h ) set (SOURCE_GROUP_SOURCE_FILES @@ -59,6 +60,7 @@ ${CMAKE_CURRENT_LIST_DIR}/RiaTimeHistoryCurveMerger.cpp ${CMAKE_CURRENT_LIST_DIR}/RiaCurveDataTools.cpp ${CMAKE_CURRENT_LIST_DIR}/RiaTimeHistoryCurveResampler.cpp ${CMAKE_CURRENT_LIST_DIR}/RiaStatisticsTools.cpp +${CMAKE_CURRENT_LIST_DIR}/RiaPolyArcLineSampler.cpp ) list(APPEND CODE_HEADER_FILES diff --git a/ApplicationCode/Application/Tools/RiaPolyArcLineSampler.cpp b/ApplicationCode/Application/Tools/RiaPolyArcLineSampler.cpp new file mode 100644 index 0000000000..a04af42980 --- /dev/null +++ b/ApplicationCode/Application/Tools/RiaPolyArcLineSampler.cpp @@ -0,0 +1,191 @@ +///////////////////////////////////////////////////////////////////////////////// +// +// Copyright (C) 2018- Statoil ASA +// +// ResInsight 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 3 of the License, or +// (at your option) any later version. +// +// ResInsight 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 at +// for more details. +// +///////////////////////////////////////////////////////////////////////////////// + +#include "RiaPolyArcLineSampler.h" +#include "cvfGeometryTools.h" +#include "cvfMatrix4.h" + + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +RiaPolyArcLineSampler::RiaPolyArcLineSampler(const std::vector& lineArcEndPoints) + : m_lineArcEndPoints(lineArcEndPoints) + , m_samplingsInterval(0.15) + , m_isResamplingLines(true) + , m_totalMD(0.0) + , m_points(nullptr) + , m_meshDs(nullptr) +{ + m_lineArcEndPoints = lineArcEndPoints; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- + +void RiaPolyArcLineSampler::sampledPointsAndMDs(double sampleInterval, + bool isResamplingLines, + std::vector* points, + std::vector* mds) +{ + CVF_ASSERT(sampleInterval > 0.0); + + m_samplingsInterval = sampleInterval; + m_isResamplingLines = isResamplingLines; + + double startMD = 0.0; + points->clear(); + mds->clear(); + + if (m_lineArcEndPoints.size() < 2) return ; + + m_points = points; + m_meshDs = mds; + m_totalMD = startMD; + + cvf::Vec3d p1 = m_lineArcEndPoints[0]; + cvf::Vec3d p2 = m_lineArcEndPoints[1]; + + m_points->push_back(p1); + m_meshDs->push_back(m_totalMD); + + cvf::Vec3d t2; + + sampleLine(p1, p2, &t2); + + for (size_t pIdx = 1; pIdx < m_lineArcEndPoints.size() - 1 ; ++pIdx) + { + sampleSegment(t2, m_lineArcEndPoints[pIdx], m_lineArcEndPoints[pIdx + 1] , &t2); + } + + + return ; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RiaPolyArcLineSampler::sampleSegment(cvf::Vec3d t1, cvf::Vec3d p1, cvf::Vec3d p2, cvf::Vec3d* endTangent) +{ + cvf::Vec3d p1p2 = p2 - p1; + if (cvf::GeometryTools::getAngle(t1, p1p2) < 1e-9) + { + sampleLine(p1, p2, endTangent); + } + else // resample arc + { + sampleArc(t1, p1, p2, endTangent); + } +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RiaPolyArcLineSampler::sampleLine(cvf::Vec3d p1, cvf::Vec3d p2, cvf::Vec3d* endTangent ) +{ + cvf::Vec3d p1p2 = p2 - p1; + + double p1p2Length = p1p2.length(); + if ( p1p2Length > m_samplingsInterval && m_isResamplingLines ) + { + cvf::Vec3d tp1p2 = p1p2 / p1p2Length; + double mdInc = m_samplingsInterval; + while ( mdInc < p1p2Length ) + { + cvf::Vec3d ps = p1 + mdInc * tp1p2; + m_points->push_back(ps); + m_meshDs->push_back(m_totalMD + mdInc); + mdInc += m_samplingsInterval; + } + } + m_totalMD += p1p2Length; + m_points->push_back(p2); + m_meshDs->push_back(m_totalMD); + + (*endTangent) = p1p2.getNormalized(); +} + +std::pair calculateArcCSAndRadius(cvf::Vec3d t1, cvf::Vec3d p1, cvf::Vec3d p2); + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RiaPolyArcLineSampler::sampleArc(cvf::Vec3d t1, cvf::Vec3d p1, cvf::Vec3d p2, cvf::Vec3d* endTangent) +{ + // Find arc CS + auto CS_rad = calculateArcCSAndRadius(t1, p1, p2); + + // Find sampleLength angle + + double angleInc = m_samplingsInterval/ CS_rad.second; + + cvf::Vec3d C = CS_rad.first.translation(); + cvf::Vec3d N(CS_rad.first.col(2)); + cvf::Vec3d tr2 = (C - p2).getNormalized(); + cvf::Vec3d t2 = tr2 ^ N; + + // Sample arc by + // Rotate vector an increment, and transform to arc CS + + double arcAngle = cvf::GeometryTools::getAngle(N, p1-C, p2-C); + + for ( double angle = angleInc; angle < arcAngle; angle += angleInc ) + { + cvf::Vec3d C_to_incP = cvf::Vec3d::X_AXIS; + C_to_incP *= CS_rad.second; + C_to_incP.transformVector(cvf::Mat3d::fromRotation(cvf::Vec3d::Z_AXIS, angle)); + + C_to_incP.transformPoint(CS_rad.first); + + m_points->push_back(C_to_incP); + m_meshDs->push_back(m_totalMD + angle * CS_rad.second); + + } + m_totalMD += arcAngle*CS_rad.second; + m_points->push_back(p2); + m_meshDs->push_back(m_totalMD); + + (*endTangent) = t2; +} + +//-------------------------------------------------------------------------------------------------- +/// + p1 +/// t1 // +/// | + C +/// \ +/// + p2 +//-------------------------------------------------------------------------------------------------- +std::pair calculateArcCSAndRadius(cvf::Vec3d t1, cvf::Vec3d p1, cvf::Vec3d p2) +{ + t1.normalize(); + cvf::Vec3d p1p2 = p2 - p1; + cvf::Vec3d t12 = p1p2.getNormalized(); + cvf::Vec3d N = (t1 ^ t12).getNormalized(); + cvf::Vec3d tr1 = (N ^ t1).getNormalized(); + double radius = 0.5*p1p2.length()/(tr1.dot(t12)); + cvf::Vec3d C = p1 + radius * tr1; + + cvf::Vec3d nTr1 = -tr1; + cvf::Mat4d CS = cvf::Mat4d::fromCoordSystemAxes(&nTr1, &t1, &N); + CS.setTranslation(C); + + return std::make_pair(CS, radius); +} + + diff --git a/ApplicationCode/Application/Tools/RiaPolyArcLineSampler.h b/ApplicationCode/Application/Tools/RiaPolyArcLineSampler.h new file mode 100644 index 0000000000..89c9c3c370 --- /dev/null +++ b/ApplicationCode/Application/Tools/RiaPolyArcLineSampler.h @@ -0,0 +1,50 @@ +///////////////////////////////////////////////////////////////////////////////// +// +// Copyright (C) 2018- Statoil ASA +// +// ResInsight 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 3 of the License, or +// (at your option) any later version. +// +// ResInsight 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 at +// for more details. +// +///////////////////////////////////////////////////////////////////////////////// +#pragma once + +#include "cvfBase.h" +#include "cvfVector3.h" + +#include + +class RiaPolyArcLineSampler +{ +public: + RiaPolyArcLineSampler(const std::vector& lineArcEndPoints); + + + void sampledPointsAndMDs(double sampleInterval, + bool isResamplingLines, + std::vector* points, + std::vector* mds); + +private: + void sampleLine(cvf::Vec3d p1, cvf::Vec3d p2, cvf::Vec3d* endTangent); + void sampleArc(cvf::Vec3d t1, cvf::Vec3d p1, cvf::Vec3d p2, cvf::Vec3d* endTangent); + void sampleSegment(cvf::Vec3d t1, cvf::Vec3d p1, cvf::Vec3d p2, cvf::Vec3d* endTangent); + + std::vector* m_points; // Internal temporary pointers to collections beeing filled. + std::vector* m_meshDs; + + double m_samplingsInterval; + bool m_isResamplingLines; + double m_totalMD; + + std::vector m_lineArcEndPoints; +}; + diff --git a/ApplicationCode/UnitTests/CMakeLists_files.cmake b/ApplicationCode/UnitTests/CMakeLists_files.cmake index c0e628e837..bbc4df74bf 100644 --- a/ApplicationCode/UnitTests/CMakeLists_files.cmake +++ b/ApplicationCode/UnitTests/CMakeLists_files.cmake @@ -45,6 +45,7 @@ ${CMAKE_CURRENT_LIST_DIR}/RigWellLogExtractor-Test.cpp ${CMAKE_CURRENT_LIST_DIR}/RifEclipseSummaryAddress-Test.cpp ${CMAKE_CURRENT_LIST_DIR}/RiaTimeHistoryCurveTools-Test.cpp ${CMAKE_CURRENT_LIST_DIR}/SolveSpaceSolver-Test.cpp +${CMAKE_CURRENT_LIST_DIR}/RiaPolyArcLineSampler-Test.cpp ) list(APPEND CODE_HEADER_FILES diff --git a/ApplicationCode/UnitTests/RiaPolyArcLineSampler-Test.cpp b/ApplicationCode/UnitTests/RiaPolyArcLineSampler-Test.cpp new file mode 100644 index 0000000000..e6d46dc544 --- /dev/null +++ b/ApplicationCode/UnitTests/RiaPolyArcLineSampler-Test.cpp @@ -0,0 +1,33 @@ +#include "gtest/gtest.h" + +#include "RiaPolyArcLineSampler.h" +#include + +//-------------------------------------------------------------------------------------------------- +TEST(RiaPolyArcLineSampler, Basic) +{ + std::vector points { {0,0,0}, {0,0,-10}, {0,10,-20},{0,20,-20}, {0,30, -30} }; + + RiaPolyArcLineSampler sampler(points); + + std::vector sampledPoints; + std::vector mds; + + sampler.sampledPointsAndMDs(2, true, &sampledPoints, &mds); + #if 0 + for (size_t pIdx = 0; pIdx < sampledPoints.size(); ++pIdx) + { + std::cout << sampledPoints[pIdx].x() << " " + << sampledPoints[pIdx].y() << " " + << sampledPoints[pIdx].z() << " md: " << mds[pIdx] << std::endl; + } + #endif + EXPECT_EQ(27, (int)sampledPoints.size()); + EXPECT_NEAR(51.4159, mds.back(), 1e-4); + + EXPECT_NEAR(points[2].y(), sampledPoints[12].y(), 2); + EXPECT_NEAR(points[2].z(), sampledPoints[12].z(), 2); + + EXPECT_NEAR(points[4].y(), sampledPoints[25].y(), 2); + EXPECT_NEAR(points[4].z(), sampledPoints[25].z(), 2); +}