OilfieldToolsNet/ResInsight
4
0
Fork 0
mirror of https://github.com/OPM/ResInsight.git synced 2025-02-25 18:55:39 -06:00
Code Issues Packages Projects Releases Wiki Activity

#2825 Fix curve stitching issues against 3dwell log curve background.

Browse Source
This commit is contained in:
Gaute Lindkvist
2018-05-02 10:21:27 +02:00
parent 1f44ea1ea3
commit 6e3d989b2f
9 changed files with 381 additions and 124 deletions
Download Patch File Download Diff File Expand all files Collapse all files

207
ApplicationCode/ModelVisualization/Riv3dWellLogCurveGeomertyGenerator.cpp
View File

@@ -25,6 +25,7 @@
#include "RigWellPath.h"
#include "RigWellPathGeometryTools.h"
#include "cafLine.h"
#include "cafDisplayCoordTransform.h"
#include "cvfPrimitiveSetIndexedUInt.h"
@@ -49,18 +50,19 @@ void Riv3dWellLogCurveGeometryGenerator::createCurveDrawables(const caf::Display
const cvf::BoundingBox& wellPathClipBoundingBox,
const Rim3dWellLogCurve* rim3dWellLogCurve,
double planeOffsetFromWellPathCenter,
double planeWidth)
double planeWidth,
const std::vector<cvf::Vec3f>& gridVertices)
{
// Make sure all drawables are cleared in case we return early to avoid a
// previous drawable being "stuck" when changing result type.
clearCurvePointsAndGeometry();
std::vector<double> resultValues;
std::vector<double> resultMds;
rim3dWellLogCurve->curveValuesAndMds(&resultValues, &resultMds);
m_planeWidth = planeWidth;
// Make sure all drawables are cleared in case we return early to avoid a
// previous drawable being "stuck" when changing result type.
clearCurvePointsAndGeometry();
if (!wellPathGeometry()) return;
if (wellPathGeometry()->m_wellPathPoints.empty()) return;
if (!wellPathClipBoundingBox.isValid()) return;
@@ -80,7 +82,7 @@ void Riv3dWellLogCurveGeometryGenerator::createCurveDrawables(const caf::Display
if (wellPathCollection->wellPathClip)
{
double clipZDistance = wellPathCollection->wellPathClipZDistance;
clipLocation = wellPathClipBoundingBox.max() + clipZDistance * cvf::Vec3d(0, 0, 1);
clipLocation = wellPathClipBoundingBox.max() + clipZDistance * cvf::Vec3d(0, 0, 1);
}
clipLocation = displayCoordTransform->transformToDisplayCoord(clipLocation);
@@ -90,14 +92,15 @@ void Riv3dWellLogCurveGeometryGenerator::createCurveDrawables(const caf::Display
wellPathPoint = displayCoordTransform->transformToDisplayCoord(wellPathPoint);
}
std::vector<cvf::Vec3d> wellPathCurveNormals = RigWellPathGeometryTools::calculateLineSegmentNormals(wellPathPoints, rim3dWellLogCurve->drawPlaneAngle());
std::vector<cvf::Vec3d> wellPathCurveNormals =
RigWellPathGeometryTools::calculateLineSegmentNormals(wellPathPoints, rim3dWellLogCurve->drawPlaneAngle());
std::vector<cvf::Vec3d> interpolatedWellPathPoints;
std::vector<cvf::Vec3d> interpolatedCurveNormals;
// Iterate from bottom of well path and up to be able to stop at given Z max clipping height
for (auto md = resultMds.rbegin(); md != resultMds.rend(); md++)
{
cvf::Vec3d point = wellPathGeometry()->interpolatedVectorAlongWellPath(wellPathPoints, *md);
cvf::Vec3d point = wellPathGeometry()->interpolatedVectorAlongWellPath(wellPathPoints, *md);
cvf::Vec3d normal = wellPathGeometry()->interpolatedVectorAlongWellPath(wellPathCurveNormals, *md);
if (point.z() > clipLocation.z()) break;
@@ -111,10 +114,8 @@ void Riv3dWellLogCurveGeometryGenerator::createCurveDrawables(const caf::Display
std::reverse(interpolatedCurveNormals.begin(), interpolatedCurveNormals.end());
// The result values for the part of the well which is not clipped off, matching interpolatedWellPathPoints size
m_curveValues = std::vector<double>(resultValues.end() - interpolatedWellPathPoints.size(),
resultValues.end());
m_curveMeasuredDepths = std::vector<double>(resultMds.end() - interpolatedWellPathPoints.size(),
resultMds.end());
m_curveValues = std::vector<double>(resultValues.end() - interpolatedWellPathPoints.size(), resultValues.end());
m_curveMeasuredDepths = std::vector<double>(resultMds.end() - interpolatedWellPathPoints.size(), resultMds.end());
double maxClampedResult = -HUGE_VAL;
double minClampedResult = HUGE_VAL;
@@ -145,27 +146,26 @@ void Riv3dWellLogCurveGeometryGenerator::createCurveDrawables(const caf::Display
if (RigCurveDataTools::isValidValue(m_curveValues[i], false))
{
scaledResult =
planeOffsetFromWellPathCenter + (m_curveValues[i] - minClampedResult) * plotRangeToResultRangeFactor;
scaledResult = planeOffsetFromWellPathCenter + (m_curveValues[i] - minClampedResult) * plotRangeToResultRangeFactor;
}
cvf::Vec3d curvePoint(interpolatedWellPathPoints[i] + scaledResult * interpolatedCurveNormals[i]);
m_curveVertices.push_back(cvf::Vec3f(curvePoint));
}
createNewVerticesAlongTriangleEdges(gridVertices);
projectVerticesOntoTriangles(gridVertices);
std::vector<cvf::uint> indices;
indices.reserve(interpolatedWellPathPoints.size());
for (size_t i = 0; i < m_curveValues.size() - 1; ++i)
indices.reserve(m_curveVertices.size() * 2);
for (size_t i = 0; i < m_curveVertices.size() - 1; ++i)
{
if (RigCurveDataTools::isValidValue(m_curveValues[i], false) &&
RigCurveDataTools::isValidValue(m_curveValues[i + 1], false))
{
indices.push_back(cvf::uint(i));
indices.push_back(cvf::uint(i + 1));
}
indices.push_back(cvf::uint(i));
indices.push_back(cvf::uint(i + 1));
}
cvf::ref<cvf::PrimitiveSetIndexedUInt> indexedUInt = new cvf::PrimitiveSetIndexedUInt(cvf::PrimitiveType::PT_LINES);
cvf::ref<cvf::UIntArray> indexArray = new cvf::UIntArray(indices);
cvf::ref<cvf::UIntArray> indexArray = new cvf::UIntArray(indices);
m_curveDrawable = new cvf::DrawableGeo();
@@ -173,7 +173,7 @@ void Riv3dWellLogCurveGeometryGenerator::createCurveDrawables(const caf::Display
m_curveDrawable->addPrimitiveSet(indexedUInt.p());
cvf::ref<cvf::Vec3fArray> vertexArray = new cvf::Vec3fArray(m_curveVertices);
m_curveDrawable->setVertexArray(vertexArray.p());
m_curveDrawable->setVertexArray(vertexArray.p());
}
//--------------------------------------------------------------------------------------------------
@@ -212,10 +212,10 @@ bool Riv3dWellLogCurveGeometryGenerator::findClosestPointOnCurve(const cvf::Vec3
double* valueAtClosestPoint) const
{
cvf::Vec3f globalIntersectionFloat(globalIntersection);
float closestDistance = m_planeWidth * 0.1;
*closestPoint = cvf::Vec3d::UNDEFINED;
*measuredDepthAtPoint = cvf::UNDEFINED_DOUBLE;
*valueAtClosestPoint = cvf::UNDEFINED_DOUBLE;
float closestDistance = m_planeWidth * 0.1;
*closestPoint = cvf::Vec3d::UNDEFINED;
*measuredDepthAtPoint = cvf::UNDEFINED_DOUBLE;
*valueAtClosestPoint = cvf::UNDEFINED_DOUBLE;
if (m_curveVertices.size() < 2u) false;
CVF_ASSERT(m_curveVertices.size() == m_curveValues.size());
for (size_t i = 1; i < m_curveVertices.size(); ++i)
@@ -224,27 +224,154 @@ bool Riv3dWellLogCurveGeometryGenerator::findClosestPointOnCurve(const cvf::Vec3
RigCurveDataTools::isValidValue(m_curveValues[i - 1], false);
if (validCurveSegment)
{
cvf::Vec3f a = m_curveVertices[i - 1];
cvf::Vec3f b = m_curveVertices[i];
cvf::Vec3f a = m_curveVertices[i - 1];
cvf::Vec3f b = m_curveVertices[i];
cvf::Vec3f ap = globalIntersectionFloat - a;
cvf::Vec3f ab = b - a;
// Projected point is clamped to one of the end points of the segment.
float distanceToProjectedPointAlongAB = ap * ab / (ab * ab);
float clampedDistance = cvf::Math::clamp(distanceToProjectedPointAlongAB, 0.0f, 1.0f);
cvf::Vec3f projectionOfGlobalIntersection = a + clampedDistance * ab;
float distance = (projectionOfGlobalIntersection - globalIntersectionFloat).length();
float distanceToProjectedPointAlongAB = ap * ab / (ab * ab);
float clampedDistance = cvf::Math::clamp(distanceToProjectedPointAlongAB, 0.0f, 1.0f);
cvf::Vec3f projectionOfGlobalIntersection = a + clampedDistance * ab;
float distance = (projectionOfGlobalIntersection - globalIntersectionFloat).length();
if (distance < closestDistance)
{
*closestPoint = cvf::Vec3d(projectionOfGlobalIntersection);
*closestPoint = cvf::Vec3d(projectionOfGlobalIntersection);
closestDistance = distance;
*measuredDepthAtPoint = m_curveMeasuredDepths[i - 1] * (1.0f - clampedDistance) + m_curveMeasuredDepths[i] * clampedDistance;
*measuredDepthAtPoint =
m_curveMeasuredDepths[i - 1] * (1.0f - clampedDistance) + m_curveMeasuredDepths[i] * clampedDistance;
*valueAtClosestPoint = m_curveValues[i - 1] * (1.0f - clampedDistance) + m_curveValues[i] * clampedDistance;
}
}
}
if (closestPoint->isUndefined())
return false;
if (closestPoint->isUndefined()) return false;
return true;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void Riv3dWellLogCurveGeometryGenerator::createNewVerticesAlongTriangleEdges(const std::vector<cvf::Vec3f>& gridVertices)
{
std::vector<cvf::Vec3f> expandedCurveVertices;
std::vector<double> expandedMeasuredDepths;
std::vector<double> expandedValues;
size_t estimatedNumberOfPoints = m_curveVertices.size() + gridVertices.size();
expandedCurveVertices.reserve(estimatedNumberOfPoints);
expandedMeasuredDepths.reserve(estimatedNumberOfPoints);
expandedValues.reserve(estimatedNumberOfPoints);
for (size_t i = 0; i < m_curveVertices.size() - 1; i += 2)
{
if (RigCurveDataTools::isValidValue(m_curveValues[i], false) &&
RigCurveDataTools::isValidValue(m_curveValues[i + 1], false))
{
expandedCurveVertices.push_back(m_curveVertices[i]);
expandedMeasuredDepths.push_back(m_curveMeasuredDepths[i]);
expandedValues.push_back(m_curveValues[i]);
// Find segments that intersects the triangle edge
caf::Line<float> curveLine(m_curveVertices[i], m_curveVertices[i + 1]);
for (size_t j = 0; j < gridVertices.size() - 1; ++j)
{
caf::Line<float> gridLine(gridVertices[j], gridVertices[j + 1]);
bool withinSegments = false;
caf::Line<float> connectingLine = curveLine.findLineBetweenNearestPoints(gridLine, &withinSegments);
if (withinSegments)
{
cvf::Vec3f closestGridPoint = connectingLine.end();
double measuredDepth;
double valueAtPoint;
cvf::Vec3d closestPoint(closestGridPoint);
cvf::Vec3d dummyArgument;
// Interpolate measured depth and value
bool worked = findClosestPointOnCurve(closestPoint, &dummyArgument, &measuredDepth, &valueAtPoint);
if (worked)
{
expandedCurveVertices.push_back(closestGridPoint);
expandedMeasuredDepths.push_back(measuredDepth);
expandedValues.push_back(valueAtPoint);
}
}
}
// Next original segment point
expandedCurveVertices.push_back(m_curveVertices[i + 1]);
expandedMeasuredDepths.push_back(m_curveMeasuredDepths[i + 1]);
expandedValues.push_back(m_curveValues[i + 1]);
}
}
m_curveVertices.swap(expandedCurveVertices);
m_curveMeasuredDepths.swap(expandedMeasuredDepths);
m_curveValues.swap(expandedValues);
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void Riv3dWellLogCurveGeometryGenerator::projectVerticesOntoTriangles(const std::vector<cvf::Vec3f>& gridVertices)
{
for (size_t i = 0; i < m_curveVertices.size(); ++i)
{
for (size_t j = 0; j < gridVertices.size() - 2; j += 1)
{
cvf::Vec3f triangleVertex1, triangleVertex2, triangleVertex3;
if (j % 2 == 0)
{
triangleVertex1 = gridVertices[j];
triangleVertex2 = gridVertices[j + 1];
triangleVertex3 = gridVertices[j + 2];
}
else
{
triangleVertex1 = gridVertices[j];
triangleVertex2 = gridVertices[j + 2];
triangleVertex3 = gridVertices[j + 1];
}
bool wasInsideTriangle = false;
cvf::Vec3f projectedPoint = projectPointOntoTriangle(
m_curveVertices[i], triangleVertex1, triangleVertex2, triangleVertex3, &wasInsideTriangle);
if (wasInsideTriangle)
{
m_curveVertices[i] = projectedPoint;
}
}
}
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
cvf::Vec3f Riv3dWellLogCurveGeometryGenerator::projectPointOntoTriangle(const cvf::Vec3f& point,
const cvf::Vec3f& triangleVertex1,
const cvf::Vec3f& triangleVertex2,
const cvf::Vec3f& triangleVertex3,
bool* wasInsideTriangle)
{
*wasInsideTriangle = false;
cvf::Vec3f e1 = triangleVertex2 - triangleVertex1;
cvf::Vec3f e2 = triangleVertex3 - triangleVertex1;
cvf::Vec3f n = (e1.getNormalized() ^ e2.getNormalized()).getNormalized();
// Project vertex onto triangle plane
cvf::Vec3f av = point - triangleVertex1;
cvf::Vec3f projectedPoint = point - (av * n) * n;
// Calculate barycentric coordinates
float areaABC = n * (e1 ^ e2);
float areaPBC = n * ((triangleVertex2 - projectedPoint) ^ (triangleVertex3 - projectedPoint));
float areaPCA = n * ((triangleVertex3 - projectedPoint) ^ (triangleVertex1 - projectedPoint));
float u = areaPBC / areaABC;
float v = areaPCA / areaABC;
float w = 1.0 - u - v;
if (u >= 0.0 && v >= 0.0 && w >= 0.0)
{
*wasInsideTriangle = true;
}
return projectedPoint;
}

12
ApplicationCode/ModelVisualization/Riv3dWellLogCurveGeomertyGenerator.h
View File

@@ -51,8 +51,9 @@ public:
const cvf::BoundingBox& wellPathClipBoundingBox,
const Rim3dWellLogCurve* rim3dWellLogCurve,
double planeOffsetFromWellPathCenter,
double planeWidth);
double planeWidth,
const std::vector<cvf::Vec3f>& gridVertices);
void clearCurvePointsAndGeometry();
const RigWellPath* wellPathGeometry() const;
@@ -65,6 +66,13 @@ public:
double* valueAtClosestPoint) const;
private:
void createNewVerticesAlongTriangleEdges(const std::vector<cvf::Vec3f>& gridVertices);
void projectVerticesOntoTriangles(const std::vector<cvf::Vec3f>& gridVertices);
static cvf::Vec3f projectPointOntoTriangle(const cvf::Vec3f& point,
const cvf::Vec3f& triangleVertex1,
const cvf::Vec3f& triangleVertex2,
const cvf::Vec3f& triangleVertex3,
bool* wasInsideTriangle);
caf::PdmPointer<RimWellPath> m_wellPath;
double m_planeWidth;

18
ApplicationCode/ModelVisualization/Riv3dWellLogPlanePartMgr.cpp
View File

@@ -72,8 +72,15 @@ void Riv3dWellLogPlanePartMgr::appendPlaneToModel(cvf::ModelBasicList*
for (Rim3dWellLogCurve* rim3dWellLogCurve : m_wellPath->rim3dWellLogCurveCollection()->vectorOf3dWellLogCurves())
{
appendGridToModel(model, displayCoordTransform, wellPathClipBoundingBox, rim3dWellLogCurve, planeWidth());
append3dWellLogCurveToModel(model, displayCoordTransform, wellPathClipBoundingBox, rim3dWellLogCurve);
if (rim3dWellLogCurve->isShowingCurve())
{
appendGridToModel(model, displayCoordTransform, wellPathClipBoundingBox, rim3dWellLogCurve, planeWidth());
append3dWellLogCurveToModel(model,
displayCoordTransform,
wellPathClipBoundingBox,
rim3dWellLogCurve,
m_3dWellLogGridGeometryGenerator->vertices());
}
}
}
@@ -83,10 +90,10 @@ void Riv3dWellLogPlanePartMgr::appendPlaneToModel(cvf::ModelBasicList*
void Riv3dWellLogPlanePartMgr::append3dWellLogCurveToModel(cvf::ModelBasicList* model,
const caf::DisplayCoordTransform* displayCoordTransform,
const cvf::BoundingBox& wellPathClipBoundingBox,
Rim3dWellLogCurve* rim3dWellLogCurve)
Rim3dWellLogCurve* rim3dWellLogCurve,
const std::vector<cvf::Vec3f>& gridVertices)
{
CVF_ASSERT(rim3dWellLogCurve);
if (!rim3dWellLogCurve->isShowingCurve()) return;
cvf::ref<Riv3dWellLogCurveGeometryGenerator> generator = rim3dWellLogCurve->geometryGenerator();
if (generator.isNull())
@@ -99,7 +106,8 @@ void Riv3dWellLogPlanePartMgr::append3dWellLogCurveToModel(cvf::ModelBasicList*
wellPathClipBoundingBox,
rim3dWellLogCurve,
wellPathCenterToPlotStartOffset(rim3dWellLogCurve),
planeWidth());
planeWidth(),
gridVertices);
cvf::ref<cvf::DrawableGeo> curveDrawable = generator->curveDrawable();

3
ApplicationCode/ModelVisualization/Riv3dWellLogPlanePartMgr.h
View File

@@ -59,7 +59,8 @@ private:
void append3dWellLogCurveToModel(cvf::ModelBasicList* model,
const caf::DisplayCoordTransform* displayCoordTransform,
const cvf::BoundingBox& wellPathClipBoundingBox,
Rim3dWellLogCurve* rim3dWellLogCurve);
Rim3dWellLogCurve* rim3dWellLogCurve,
const std::vector<cvf::Vec3f>& gridVertices);
void appendGridToModel(cvf::ModelBasicList* model,
const caf::DisplayCoordTransform* displayCoordTransform,