diff --git a/ApplicationCode/ModelVisualization/Intersections/CMakeLists_files.cmake b/ApplicationCode/ModelVisualization/Intersections/CMakeLists_files.cmake index 5d2e68a6d6..1b0684cc7c 100644 --- a/ApplicationCode/ModelVisualization/Intersections/CMakeLists_files.cmake +++ b/ApplicationCode/ModelVisualization/Intersections/CMakeLists_files.cmake @@ -7,6 +7,7 @@ ${CMAKE_CURRENT_LIST_DIR}/RivHexGridIntersectionTools.h ${CMAKE_CURRENT_LIST_DIR}/RivIntersectionBoxGeometryGenerator.h ${CMAKE_CURRENT_LIST_DIR}/RivIntersectionBoxPartMgr.h ${CMAKE_CURRENT_LIST_DIR}/RivIntersectionBoxSourceInfo.h +${CMAKE_CURRENT_LIST_DIR}/RivSectionFlattner.h ) set (SOURCE_GROUP_SOURCE_FILES @@ -17,6 +18,7 @@ ${CMAKE_CURRENT_LIST_DIR}/RivHexGridIntersectionTools.cpp ${CMAKE_CURRENT_LIST_DIR}/RivIntersectionBoxGeometryGenerator.cpp ${CMAKE_CURRENT_LIST_DIR}/RivIntersectionBoxPartMgr.cpp ${CMAKE_CURRENT_LIST_DIR}/RivIntersectionBoxSourceInfo.cpp +${CMAKE_CURRENT_LIST_DIR}/RivSectionFlattner.cpp ) list(APPEND CODE_HEADER_FILES diff --git a/ApplicationCode/ModelVisualization/Intersections/RivSectionFlattner.cpp b/ApplicationCode/ModelVisualization/Intersections/RivSectionFlattner.cpp new file mode 100644 index 0000000000..ab4aa91419 --- /dev/null +++ b/ApplicationCode/ModelVisualization/Intersections/RivSectionFlattner.cpp @@ -0,0 +1,150 @@ +///////////////////////////////////////////////////////////////////////////////// +// +// 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 "RivSectionFlattner.h" +#include "cvfGeometryTools.h" + + +//-------------------------------------------------------------------------------------------------- +/// Returns the next index higher than idxToStartOfLineSegment that makes the line +// polyline[idxToStartOfLineSegment] .. polyline[nextIdx] not parallel to extrDir +/// +/// Returns size_t(-1) if no point is found +//-------------------------------------------------------------------------------------------------- +size_t RivSectionFlattner::indexToNextValidPoint(const std::vector& polyLine, const cvf::Vec3d extrDir, size_t idxToStartOfLineSegment) +{ + size_t lineCount = polyLine.size(); + if ( !(idxToStartOfLineSegment + 1 < lineCount) ) return -1; + + + cvf::Vec3d p1 = polyLine[idxToStartOfLineSegment]; + + for ( size_t lIdx = idxToStartOfLineSegment+1; lIdx < lineCount; ++lIdx ) + { + cvf::Vec3d p2 = polyLine[lIdx]; + cvf::Vec3d p1p2 = p2 - p1; + + if ( (p1p2 - (p1p2 * extrDir)*extrDir).length() > 0.1 ) + { + return lIdx; + } + } + + return -1; +} + +//-------------------------------------------------------------------------------------------------- +/// Returns one CS pr point, valid for the next segment +//-------------------------------------------------------------------------------------------------- +std::vector RivSectionFlattner::calculateFlatteningCSsForPolyline(const std::vector & polyLine, const cvf::Vec3d& extrusionDir, const cvf::Vec3d& startOffset, cvf::Vec3d* endOffset) +{ + CVF_ASSERT(endOffset); + size_t pointCount = polyLine.size(); + CVF_ASSERT(pointCount > 1); + + std::vector segmentTransforms; + segmentTransforms.reserve(pointCount); + + // Find initial transform, used if all is vertical + + cvf::Mat4d invSectionCS; + { + cvf::Vec3d p1 = polyLine[0]; + cvf::Vec3d p2 = polyLine[1]; + + cvf::Mat4d sectionLocalCS = calculateSectionLocalFlatteningCS(p1, p2, extrusionDir); + cvf::Mat4d invSectionCS = sectionLocalCS.getInverted(); + invSectionCS.setTranslation(invSectionCS.translation() + startOffset); + } + + cvf::Vec3d previousFlattenedSectionEndPoint = startOffset; + + size_t lIdx = 0; + while ( lIdx < pointCount ) + { + size_t idxToNextP = indexToNextValidPoint(polyLine, extrusionDir, lIdx); + + // If the rest is nearly parallel to extrusionDir, use the current inverse matrix for the rest of the points + + if ( idxToNextP == size_t(-1) ) + { + size_t inc = 0; + while ( (lIdx + inc) < pointCount ) + { + segmentTransforms.push_back(invSectionCS); + ++inc; + } + break; + } + + cvf::Vec3d p1 = polyLine[lIdx]; + cvf::Vec3d p2 = polyLine[idxToNextP]; + + cvf::Mat4d sectionLocalCS = calculateSectionLocalFlatteningCS(p1, p2, extrusionDir); + invSectionCS = sectionLocalCS.getInverted(); + cvf::Vec3d flattenedSectionEndPoint = p2.getTransformedPoint(invSectionCS); + + invSectionCS.setTranslation(invSectionCS.translation() + previousFlattenedSectionEndPoint); + + previousFlattenedSectionEndPoint += flattenedSectionEndPoint; + + // Assign the matrix to the points in between + + size_t inc = 0; + while ( (lIdx + inc) < idxToNextP ) + { + segmentTransforms.push_back(invSectionCS); + inc++; + } + + lIdx = idxToNextP; + } + + *endOffset = previousFlattenedSectionEndPoint; + return segmentTransforms; +} + +//-------------------------------------------------------------------------------------------------- +/// Origo in P1 +/// Ez in upwards extrusionDir +/// Ey normal to the section plane +/// Ex in plane along p1-p2 +//-------------------------------------------------------------------------------------------------- +cvf::Mat4d RivSectionFlattner::calculateSectionLocalFlatteningCS(const cvf::Vec3d& p1, const cvf::Vec3d& p2, const cvf::Vec3d& extrusionDir) +{ + using namespace cvf; + + Vec3d Ez = extrusionDir.z() > 0.0 ? extrusionDir: -extrusionDir; + + Vec3d sectionLineDir = p2 - p1; + + if ( cvf::GeometryTools::getAngle(sectionLineDir, extrusionDir) < 0.01 ) + { + sectionLineDir = Ez.perpendicularVector(); + } + + Vec3d Ey = Ez ^ sectionLineDir; + Ey.normalize(); + Vec3d Ex = Ey ^ Ez; + Ex.normalize(); + + return Mat4d(Ex[0], Ey[0], Ez[0], p1[0], + Ex[1], Ey[1], Ez[1], p1[1], + Ex[2], Ey[2], Ez[2], p1[2], + 0.0, 0.0, 0.0, 1.0); +} diff --git a/ApplicationCode/ModelVisualization/Intersections/RivSectionFlattner.h b/ApplicationCode/ModelVisualization/Intersections/RivSectionFlattner.h index c0973123b7..6d3a9288ae 100644 --- a/ApplicationCode/ModelVisualization/Intersections/RivSectionFlattner.h +++ b/ApplicationCode/ModelVisualization/Intersections/RivSectionFlattner.h @@ -1,146 +1,46 @@ +///////////////////////////////////////////////////////////////////////////////// +// +// 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 "cvfMatrix4.h" + +#include + class RivSectionFlattner { public: - //-------------------------------------------------------------------------------------------------- - /// Returns the next index higher than idxToStartOfLineSegment that makes the line - // polyline[idxToStartOfLineSegment] .. polyline[nextIdx] not parallel to extrDir - /// - /// Returns size_t(-1) if no point is found - //-------------------------------------------------------------------------------------------------- static size_t indexToNextValidPoint(const std::vector& polyLine, const cvf::Vec3d extrDir, - size_t idxToStartOfLineSegment) - { - size_t lineCount = polyLine.size(); - if ( !(idxToStartOfLineSegment + 1 < lineCount) ) return -1; + size_t idxToStartOfLineSegment); - - cvf::Vec3d p1 = polyLine[idxToStartOfLineSegment]; - - for ( size_t lIdx = idxToStartOfLineSegment+1; lIdx < lineCount; ++lIdx ) - { - cvf::Vec3d p2 = polyLine[lIdx]; - cvf::Vec3d p1p2 = p2 - p1; - - if ( (p1p2 - (p1p2 * extrDir)*extrDir).length() > 0.1 ) - { - return lIdx; - } - } - - return -1; - } - - //-------------------------------------------------------------------------------------------------- - /// Returns one CS pr point, valid for the next segment - //-------------------------------------------------------------------------------------------------- static std::vector calculateFlatteningCSsForPolyline(const std::vector & polyLine, const cvf::Vec3d& extrusionDir, const cvf::Vec3d& startOffset, - cvf::Vec3d* endOffset) - { - CVF_ASSERT(endOffset); - size_t pointCount = polyLine.size(); - CVF_ASSERT(pointCount > 1); - - std::vector segmentTransforms; - segmentTransforms.reserve(pointCount); - - // Find initial transform, used if all is vertical - - cvf::Mat4d invSectionCS; - { - cvf::Vec3d p1 = polyLine[0]; - cvf::Vec3d p2 = polyLine[1]; - - cvf::Mat4d sectionLocalCS = calculateSectionLocalFlatteningCS(p1, p2, extrusionDir); - cvf::Mat4d invSectionCS = sectionLocalCS.getInverted(); - invSectionCS.setTranslation(invSectionCS.translation() + startOffset); - } - - cvf::Vec3d previousFlattenedSectionEndPoint = startOffset; - - size_t lIdx = 0; - while ( lIdx < pointCount ) - { - size_t idxToNextP = indexToNextValidPoint(polyLine, extrusionDir, lIdx); - - // If the rest is nearly parallel to extrusionDir, use the current inverse matrix for the rest of the points - - if ( idxToNextP == size_t(-1) ) - { - size_t inc = 0; - while ( (lIdx + inc) < pointCount ) - { - segmentTransforms.push_back(invSectionCS); - ++inc; - } - break; - } - - cvf::Vec3d p1 = polyLine[lIdx]; - cvf::Vec3d p2 = polyLine[idxToNextP]; - - cvf::Mat4d sectionLocalCS = calculateSectionLocalFlatteningCS(p1, p2, extrusionDir); - invSectionCS = sectionLocalCS.getInverted(); - cvf::Vec3d flattenedSectionEndPoint = p2.getTransformedPoint(invSectionCS); - - invSectionCS.setTranslation(invSectionCS.translation() + previousFlattenedSectionEndPoint ); - - previousFlattenedSectionEndPoint += flattenedSectionEndPoint; - - // Assign the matrix to the points in between - - size_t inc = 0; - while ( (lIdx + inc) < idxToNextP ) - { - segmentTransforms.push_back(invSectionCS); - inc++; - } - - lIdx = idxToNextP; - } - - *endOffset = previousFlattenedSectionEndPoint; - return segmentTransforms; - } + cvf::Vec3d* endOffset); private: - //-------------------------------------------------------------------------------------------------- - /// Origo in P1 - /// Ez in upwards extrusionDir - /// Ey normal to the section plane - /// Ex in plane along p1-p2 - //-------------------------------------------------------------------------------------------------- static cvf::Mat4d calculateSectionLocalFlatteningCS(const cvf::Vec3d& p1, const cvf::Vec3d& p2, - const cvf::Vec3d& extrusionDir) - { - using namespace cvf; - - Vec3d Ez = extrusionDir.z() > 0.0 ? extrusionDir: -extrusionDir; - - Vec3d sectionLineDir = p2 - p1; - - if ( cvf::GeometryTools::getAngle(sectionLineDir, extrusionDir) < 0.01 ) - { - sectionLineDir = Ez.perpendicularVector(); - } - - Vec3d Ey = Ez ^ sectionLineDir; - Ey.normalize(); - Vec3d Ex = Ey ^ Ez; - Ex.normalize(); - - return Mat4d(Ex[0], Ey[0], Ez[0], p1[0], - Ex[1], Ey[1], Ez[1], p1[1], - Ex[2], Ey[2], Ez[2], p1[2], - 0.0, 0.0, 0.0, 1.0); - } + const cvf::Vec3d& extrusionDir); };