#3135 Create a class to do sampling of an poly-arc-line

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

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 <http://www.gnu.org/licenses/gpl.html> 
+//  for more details.
+//
+/////////////////////////////////////////////////////////////////////////////////
+
+#include "RiaPolyArcLineSampler.h"
+#include "cvfGeometryTools.h"
+#include "cvfMatrix4.h"
+
+
+//--------------------------------------------------------------------------------------------------
+/// 
+//--------------------------------------------------------------------------------------------------
+RiaPolyArcLineSampler::RiaPolyArcLineSampler(const std::vector<cvf::Vec3d>& 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<cvf::Vec3d>* points,
+                                                    std::vector<double>* 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<cvf::Mat4d, double> 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<cvf::Mat4d, double> 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);
+}
+
+

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 <http://www.gnu.org/licenses/gpl.html> 
+//  for more details.
+//
+/////////////////////////////////////////////////////////////////////////////////
+#pragma once
+
+#include "cvfBase.h"
+#include "cvfVector3.h"
+
+#include <vector>
+
+class RiaPolyArcLineSampler
+{
+public:
+    RiaPolyArcLineSampler(const std::vector<cvf::Vec3d>& lineArcEndPoints);
+
+
+    void sampledPointsAndMDs(double sampleInterval,
+                               bool isResamplingLines,
+                               std::vector<cvf::Vec3d>* points,
+                               std::vector<double>* 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<cvf::Vec3d>* m_points; // Internal temporary pointers to collections beeing filled.
+    std::vector<double>*     m_meshDs;
+
+    double                   m_samplingsInterval;
+    bool                     m_isResamplingLines;
+    double                   m_totalMD; 
+    
+    std::vector<cvf::Vec3d>  m_lineArcEndPoints;
+};
+

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

ApplicationCode/UnitTests/RiaPolyArcLineSampler-Test.cpp

@@ -0,0 +1,33 @@
+#include "gtest/gtest.h"
+
+#include "RiaPolyArcLineSampler.h"
+#include <iostream>
+
+//--------------------------------------------------------------------------------------------------
+TEST(RiaPolyArcLineSampler, Basic)
+{
+    std::vector<cvf::Vec3d> points {  {0,0,0}, {0,0,-10}, {0,10,-20},{0,20,-20}, {0,30, -30} }; 
+
+    RiaPolyArcLineSampler sampler(points);
+
+    std::vector<cvf::Vec3d> sampledPoints; 
+    std::vector<double> 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);
+}