#5529 Surface : Add reader for GOCAD file format

wip

ApplicationCode/FileInterface/CMakeLists_files.cmake

@@ -50,6 +50,7 @@ ${CMAKE_CURRENT_LIST_DIR}/RifDerivedEnsembleReader.h
 ${CMAKE_CURRENT_LIST_DIR}/RifActiveCellsReader.h
 ${CMAKE_CURRENT_LIST_DIR}/RifCsvDataTableFormatter.h
 ${CMAKE_CURRENT_LIST_DIR}/RifEclipseInputPropertyLoader.h
+${CMAKE_CURRENT_LIST_DIR}/RifSurfaceReader.h
 
 # HDF5 file reader is directly included in ResInsight main CmakeList.txt
 #${CMAKE_CURRENT_LIST_DIR}/RifHdf5Reader.h
@@ -104,6 +105,8 @@ ${CMAKE_CURRENT_LIST_DIR}/RifActiveCellsReader.cpp
 ${CMAKE_CURRENT_LIST_DIR}/RifCsvDataTableFormatter.cpp
 ${CMAKE_CURRENT_LIST_DIR}/RifReaderEnsembleStatisticsRft.cpp
 ${CMAKE_CURRENT_LIST_DIR}/RifEclipseInputPropertyLoader.cpp
+${CMAKE_CURRENT_LIST_DIR}/RifSurfaceReader.cpp
+
 
 # HDF5 file reader is directly included in ResInsight main CmakeList.txt
 #${CMAKE_CURRENT_LIST_DIR}/RifHdf5Reader.cpp

ApplicationCode/FileInterface/RifSurfaceReader.cpp

@@ -0,0 +1,283 @@
+/////////////////////////////////////////////////////////////////////////////////
+//
+//  Copyright (C) 2020 Equinor 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 "RifSurfaceReader.h"
+
+#include "cvfVector3.h"
+
+#include <fstream>
+#include <limits>
+#include <map>
+#include <sstream>
+#include <vector>
+
+//--------------------------------------------------------------------------------------------------
+///
+/// Import vertices and triangle IDs from the first TFACE section in the file
+/// Returns vertices with z-value as depth, z is increasing downwards
+///
+///
+///
+///
+//--------------------------------------------------------------------------------------------------
+std::pair<std::vector<cvf::Vec3d>, std::vector<unsigned>> RifSurfaceReader::readGocadFile( const QString& filename )
+{
+    enum class GocadZPositive
+    {
+        Elevation,
+        Depth,
+        Unknown
+    };
+
+    std::vector<cvf::Vec3d> vertices;
+    std::map<int, unsigned> vertexIdToIndex;
+    std::vector<unsigned>   trianglesByIds;
+
+    {
+        std::ifstream stream( filename.toLatin1().data() );
+
+        bool isInTfaceSection = false;
+        while ( stream.good() )
+        {
+            std::string line;
+            std::getline( stream, line );
+
+            std::transform( line.begin(), line.end(), line.begin(), ::toupper );
+
+            std::istringstream lineStream( line );
+
+            GocadZPositive zDir = GocadZPositive::Unknown;
+
+            std::string firstToken;
+            lineStream >> firstToken;
+
+            if ( lineStream.good() ) // If we can, assume this line is a surface point
+            {
+                if ( isInTfaceSection )
+                {
+                    if ( firstToken.compare( "VRTX" ) == 0 )
+                    {
+                        int         vertexId = -1;
+                        double      x{std::numeric_limits<double>::infinity()};
+                        double      y{std::numeric_limits<double>::infinity()};
+                        double      z{std::numeric_limits<double>::infinity()};
+                        std::string endVertex;
+
+                        lineStream >> vertexId >> x >> y >> z >> endVertex;
+
+                        if ( vertexId > -1 )
+                        {
+                            if ( zDir == GocadZPositive::Elevation ) z = -z;
+
+                            vertices.emplace_back( cvf::Vec3d( x, y, z ) );
+                            vertexIdToIndex[vertexId] = static_cast<unsigned>( vertices.size() - 1 );
+                        }
+                    }
+                    else if ( firstToken.compare( "TRGL" ) == 0 )
+                    {
+                        int id1{-1};
+                        int id2{-1};
+                        int id3{-1};
+
+                        lineStream >> id1 >> id2 >> id3;
+
+                        if ( id1 >= 0 && id2 >= 0 && id3 >= 0 )
+                        {
+                            trianglesByIds.emplace_back( static_cast<unsigned int>( id1 ) );
+                            trianglesByIds.emplace_back( static_cast<unsigned int>( id2 ) );
+                            trianglesByIds.emplace_back( static_cast<unsigned int>( id3 ) );
+                        }
+                    }
+                    else if ( firstToken.compare( "END" ) == 0 )
+                    {
+                        isInTfaceSection = false;
+                    }
+                }
+                else if ( firstToken.compare( "TFACE" ) == 0 )
+                {
+                    isInTfaceSection = true;
+                }
+                else if ( firstToken.compare( "ZPOSITIVE" ) == 0 )
+                {
+                    std::string secondToken;
+                    lineStream >> secondToken;
+
+                    if ( secondToken == "DEPTH" )
+                    {
+                        zDir = GocadZPositive::Depth;
+                    }
+                    else if ( secondToken == "ELEVATION" )
+                    {
+                        zDir = GocadZPositive::Elevation;
+                    }
+                }
+            }
+            else // Probably a comment line, skip
+            {
+            }
+        }
+    }
+
+    std::vector<unsigned> triangleIndices;
+    if ( !trianglesByIds.empty() )
+    {
+        triangleIndices.reserve( trianglesByIds.size() );
+
+        for ( auto triangle : trianglesByIds )
+        {
+            auto vertexIndex = vertexIdToIndex[triangle];
+            triangleIndices.push_back( vertexIndex );
+        }
+    }
+
+    return std::make_pair( vertices, triangleIndices );
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+std::pair<std::vector<cvf::Vec3d>, std::vector<unsigned>> RifSurfaceReader::readPetrelFile( const QString& filename )
+{
+    std::ifstream stream( filename.toLatin1().data() );
+
+    struct SurfacePointData
+    {
+        int                 i;
+        int                 j;
+        cvf::Vec3d          point;
+        std::vector<double> values;
+    };
+
+    std::vector<SurfacePointData> surfaceDataPoints;
+    int                           minI = std::numeric_limits<int>::max();
+    int                           minJ = std::numeric_limits<int>::max();
+    int                           maxI = std::numeric_limits<int>::min();
+    int                           maxJ = std::numeric_limits<int>::min();
+
+    while ( stream.good() )
+    {
+        std::string line;
+        std::getline( stream, line );
+        std::istringstream lineStream( line );
+
+        double              x( HUGE_VAL ), y( HUGE_VAL ), z( HUGE_VAL );
+        int                 i( -1 ), j( -1 );
+        std::vector<double> values;
+
+        // First check if we can read a number
+
+        lineStream >> x;
+        if ( lineStream.good() ) // If we can, assume this line is a surface point
+        {
+            lineStream >> y >> z >> i >> j;
+
+            if ( x != HUGE_VAL && y != HUGE_VAL && z != HUGE_VAL && i != -1 && j != -1 )
+            {
+                // Check for extra data
+                while ( lineStream.good() )
+                {
+                    double d;
+                    lineStream >> d;
+                    if ( lineStream.good() ) values.push_back( d );
+                }
+
+                // Add point
+
+                surfaceDataPoints.push_back( {i, j, {x, y, z}, values} );
+
+                minI = std::min( minI, i );
+                minJ = std::min( minJ, j );
+                maxI = std::max( maxI, i );
+                maxJ = std::max( maxJ, j );
+            }
+        }
+        else // Probably a comment line, skip
+        {
+        }
+    }
+
+    // clang-format off
+    if ( surfaceDataPoints.empty() 
+        || minI == std::numeric_limits<int>::max() 
+        || minJ == std::numeric_limits<int>::max() 
+        || maxI == std::numeric_limits<int>::min() 
+        || maxJ == std::numeric_limits<int>::min() )
+    {
+        return {};
+    }
+    // clang-format on
+
+    // Create full size grid matrix
+
+    size_t iCount = maxI - minI + 1;
+    size_t jCount = maxJ - minJ + 1;
+
+    if ( iCount < 2 || jCount < 2 )
+    {
+        return {};
+    }
+
+    std::vector<std::vector<unsigned>> indexToPointData;
+    indexToPointData.resize( iCount, std::vector<unsigned>( jCount, -1 ) );
+    std::vector<cvf::Vec3d> vertices;
+
+    for ( unsigned pIdx = 0; pIdx < surfaceDataPoints.size(); ++pIdx )
+    {
+        const auto& pointData = surfaceDataPoints[pIdx];
+
+        indexToPointData[pointData.i - minI][pointData.j - minJ] = pIdx;
+
+        vertices.push_back( pointData.point );
+        // Todo: Move result values for each point into the
+    }
+
+    std::vector<unsigned> triangleIndices;
+    if ( indexToPointData.size() < 2 )
+    {
+        return {};
+    }
+
+    for ( size_t iIdx = 0; iIdx < indexToPointData.size() - 1; ++iIdx )
+    {
+        for ( size_t jIdx = 0; jIdx < indexToPointData[iIdx].size() - 1; ++jIdx )
+        {
+            {
+                unsigned q1 = indexToPointData[iIdx + 0][jIdx + 0];
+                unsigned q2 = indexToPointData[iIdx + 0][jIdx + 1];
+                unsigned q3 = indexToPointData[iIdx + 1][jIdx + 0];
+                unsigned q4 = indexToPointData[iIdx + 1][jIdx + 1];
+
+                if ( q1 != ( (unsigned)-1 ) && q2 != ( (unsigned)-1 ) && q4 != ( (unsigned)-1 ) )
+                {
+                    triangleIndices.push_back( q1 );
+                    triangleIndices.push_back( q2 );
+                    triangleIndices.push_back( q4 );
+                }
+
+                if ( q1 != ( (unsigned)-1 ) && q2 != ( (unsigned)-1 ) && q3 != ( (unsigned)-1 ) )
+                {
+                    triangleIndices.push_back( q1 );
+                    triangleIndices.push_back( q4 );
+                    triangleIndices.push_back( q3 );
+                }
+            }
+        }
+    }
+
+    return std::make_pair( vertices, triangleIndices );
+}

ApplicationCode/FileInterface/RifSurfaceReader.h

@@ -0,0 +1,34 @@
+/////////////////////////////////////////////////////////////////////////////////
+//
+//  Copyright (C) 2020 Equinor 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 "cvfVector3.h"
+
+#include <utility>
+#include <vector>
+
+#include <QString>
+
+class RifSurfaceReader
+{
+public:
+    static std::pair<std::vector<cvf::Vec3d>, std::vector<unsigned>> readGocadFile( const QString& filename );
+
+    static std::pair<std::vector<cvf::Vec3d>, std::vector<unsigned>> readPetrelFile( const QString& filename );
+};