Added function for importing surfaces from OpenWorks XYZ files.

Added helper function for fuzzy comparing vectors.
Adjusted (old) code to avoid warnings.

ApplicationCode/Commands/SurfaceCommands/RicImportSurfacesFeature.cpp

@@ -54,7 +54,7 @@ void RicImportSurfacesFeature::onActionTriggered( bool isChecked )
     QStringList     fileNames  = RiuFileDialogTools::getOpenFileNames( Riu3DMainWindowTools::mainWindowWidget(),
                                                                   "Import Surfaces",
                                                                   defaultDir,
-                                                                  "Surface files (*.ptl *.ts);;All Files (*.*)" );
+                                                                  "Surface files (*.ptl *.ts *.dat);;All Files (*.*)" );
 
     if ( fileNames.isEmpty() ) return;
 

ApplicationCode/FileInterface/RifSurfaceImporter.cpp

@@ -152,7 +152,11 @@ void RifSurfaceImporter::readGocadFile( const QString& filename, RigGocadData* g
 
                 qstringLine.remove( "PROPERTIES" );
 
+#if QT_VERSION >= QT_VERSION_CHECK( 5, 15, 0 )
+                QStringList words = qstringLine.split( " ", Qt::SkipEmptyParts );
+#else
                 QStringList words = qstringLine.split( " ", QString::SkipEmptyParts );
+#endif
 
                 for ( auto w : words )
                 {
@@ -278,8 +282,8 @@ std::pair<std::vector<cvf::Vec3d>, std::vector<unsigned>> RifSurfaceImporter::re
 
     // Create full size grid matrix
 
-    size_t iCount = maxI - minI + 1;
-    size_t jCount = maxJ - minJ + 1;
+    size_t iCount = static_cast<size_t>( maxI ) - static_cast<size_t>( minI ) + 1;
+    size_t jCount = static_cast<size_t>( maxJ ) - static_cast<size_t>( minJ ) + 1;
 
     if ( iCount < 2 || jCount < 2 )
     {
@@ -294,7 +298,7 @@ std::pair<std::vector<cvf::Vec3d>, std::vector<unsigned>> RifSurfaceImporter::re
     {
         const auto& pointData = surfaceDataPoints[pIdx];
 
-        indexToPointData[pointData.i - minI][pointData.j - minJ] = pIdx;
+        indexToPointData[static_cast<size_t>( pointData.i ) - minI][static_cast<size_t>( pointData.j ) - minJ] = pIdx;
 
         vertices.push_back( pointData.point );
         // Todo: Move result values for each point into the
@@ -310,28 +314,329 @@ std::pair<std::vector<cvf::Vec3d>, std::vector<unsigned>> RifSurfaceImporter::re
     {
         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 ) && q4 != ( (unsigned)-1 ) && q3 != ( (unsigned)-1 ) )
-                {
-                    triangleIndices.push_back( q1 );
-                    triangleIndices.push_back( q4 );
-                    triangleIndices.push_back( q3 );
-                }
-            }
+            generateTriangleIndices( indexToPointData, iIdx, jIdx, triangleIndices );
         }
     }
 
     return std::make_pair( vertices, triangleIndices );
 }
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+std::pair<std::vector<cvf::Vec3d>, std::vector<unsigned>> RifSurfaceImporter::readOpenWorksXyzFile( const QString& filename )
+{
+    // Note: the assumption is that the grid will always be in the x-y plane.
+    std::ifstream stream( filename.toLatin1().data() );
+
+    struct vec2dCompare
+    {
+        bool operator()( const cvf::Vec2d& lhs, const cvf::Vec2d& rhs ) const { return lhs.length() < rhs.length(); }
+    };
+
+    std::vector<cvf::Vec3d> surfacePoints;
+
+    // Mapping normalized vectors to all vectors having their direction
+    std::map<cvf::Vec2d, double, vec2dCompare> axesVectorCandidates;
+    // Mapping axes vectors to their number of occurrence
+    std::map<cvf::Vec2d, unsigned, vec2dCompare> axesVectorCandidatesNum;
+    // Mapping axes vectors to their number of neighbored points
+    std::map<cvf::Vec2d, unsigned, vec2dCompare> numOfNeighbouredPoints;
+    // Vector of points where the vector between two following surface points changes.
+    std::vector<cvf::Vec2d> vectorChanges;
+    // Vector of potential end/start points.
+    std::vector<size_t> potentialEndPointIndices;
+
+    cvf::Vec2d lastVector;
+
+    double epsilon                = 0.001;
+    bool   ignoreNextVectorChange = false;
+
+    // Converts a 3d vector to a 2d vector
+    auto to2d = []( const cvf::Vec3d vector ) -> cvf::Vec2d { return cvf::Vec2d( vector.x(), vector.y() ); };
+
+    // Checks if the given vector is a possible new candidate for an axis vector and adds it to the given list of
+    // axesVectorCandidates. Also increases the number of occurrences of vector candidates.
+    auto maybeInsertAxisVectorCandidate =
+        [epsilon]( const cvf::Vec2d                              vector,
+                   std::map<cvf::Vec2d, double, vec2dCompare>&   axesVectorCandidates,
+                   std::map<cvf::Vec2d, unsigned, vec2dCompare>& axesVectorCandidatesNum,
+                   std::map<cvf::Vec2d, unsigned, vec2dCompare>& numOfNeighbouredPoints ) -> bool {
+        double     length           = vector.length();
+        cvf::Vec2d normalizedVector = vector.getNormalized();
+        for ( std::map<cvf::Vec2d, double, vec2dCompare>::iterator iter = axesVectorCandidates.begin();
+              iter != axesVectorCandidates.end();
+              ++iter )
+        {
+            if ( vectorFuzzyCompare( iter->first, normalizedVector, 0.1 ) )
+            {
+                normalizedVector = iter->first;
+                break;
+            }
+        }
+        if ( axesVectorCandidates.count( normalizedVector ) == 0 )
+        {
+            axesVectorCandidates.insert( std::pair<cvf::Vec2d, double>( normalizedVector, length ) );
+            axesVectorCandidatesNum.insert( std::pair<cvf::Vec2d, unsigned>( normalizedVector, 0 ) );
+            axesVectorCandidatesNum[normalizedVector] += 1;
+            return true;
+        }
+        else if ( length < axesVectorCandidates[normalizedVector] )
+        {
+            axesVectorCandidates[normalizedVector] = length;
+        }
+        if ( numOfNeighbouredPoints.count( normalizedVector ) == 0 )
+        {
+            numOfNeighbouredPoints.insert( std::pair<cvf::Vec2d, unsigned>( normalizedVector, 0 ) );
+        }
+        numOfNeighbouredPoints[normalizedVector] += 1;
+        axesVectorCandidatesNum[normalizedVector] += 1;
+        return false;
+    };
+
+    while ( stream.good() )
+    {
+        std::string line;
+        std::getline( stream, line );
+        std::istringstream lineStream( line );
+
+        double x( std::numeric_limits<double>::infinity() );
+        double y( std::numeric_limits<double>::infinity() );
+        double z( std::numeric_limits<double>::infinity() );
+
+        std::vector<double> values;
+
+        // First check if we can read all numbers
+        bool ok = true;
+        lineStream >> x;
+        ok &= lineStream.good();
+        lineStream >> y;
+        ok &= lineStream.good();
+        lineStream >> z;
+        ok &= lineStream.good() || lineStream.eof();
+        if ( ok ) // If we can, assume this line is a surface point
+        {
+            if ( x != std::numeric_limits<double>::infinity() && y != std::numeric_limits<double>::infinity() &&
+                 z != std::numeric_limits<double>::infinity() )
+            {
+                // Check for extra data
+                while ( lineStream.good() )
+                {
+                    double d;
+                    lineStream >> d;
+                    if ( lineStream.good() ) values.push_back( d );
+                }
+
+                // Add point
+                surfacePoints.push_back( {x, y, z} );
+
+                if ( surfacePoints.size() > 1 )
+                {
+                    cvf::Vec3d pointToPointVector = surfacePoints.back() - surfacePoints[surfacePoints.size() - 2];
+
+                    // Assuming that any change bigger than 45 degrees is a new column in the grid.
+                    if ( lastVector.dot( to2d( pointToPointVector ) ) < 0.5 && !ignoreNextVectorChange )
+                    {
+                        // Compare the current vector to all points where the vector has changed
+                        // in order to identify the secondary axis vector.
+                        for ( size_t i = 0; i < vectorChanges.size(); i++ )
+                        {
+                            cvf::Vec2d vector = to2d( surfacePoints.back() ) - vectorChanges[i];
+                            maybeInsertAxisVectorCandidate( vector,
+                                                            axesVectorCandidates,
+                                                            axesVectorCandidatesNum,
+                                                            numOfNeighbouredPoints );
+                        }
+                        vectorChanges.push_back( to2d( surfacePoints.back() ) );
+                        if ( !std::count( potentialEndPointIndices.begin(),
+                                          potentialEndPointIndices.end(),
+                                          surfacePoints.size() - 1 ) )
+                            potentialEndPointIndices.push_back( surfacePoints.size() - 1 );
+
+                        if ( !std::count( potentialEndPointIndices.begin(),
+                                          potentialEndPointIndices.end(),
+                                          surfacePoints.size() - 2 ) )
+                            potentialEndPointIndices.push_back( surfacePoints.size() - 2 );
+
+                        if ( vectorChanges.size() > 2 ) ignoreNextVectorChange = true;
+                    }
+                    else
+                    {
+                        if ( lastVector.dot( to2d( pointToPointVector ) ) < 0.5 && ignoreNextVectorChange )
+                        {
+                            ignoreNextVectorChange = false;
+                        }
+                        maybeInsertAxisVectorCandidate( to2d( pointToPointVector ),
+                                                        axesVectorCandidates,
+                                                        axesVectorCandidatesNum,
+                                                        numOfNeighbouredPoints );
+                    }
+
+                    lastVector = to2d( pointToPointVector ).getNormalized();
+                }
+                else
+                {
+                    vectorChanges.push_back( to2d( surfacePoints.back() ) );
+                }
+            }
+        }
+        else // Probably a comment line, skip
+        {
+        }
+    }
+
+    // Determine axes vectors
+    std::vector<std::pair<cvf::Vec2d, unsigned>> pairs;
+    for ( auto itr = axesVectorCandidatesNum.begin(); itr != axesVectorCandidatesNum.end(); ++itr )
+    {
+        pairs.push_back( *itr );
+    }
+
+    sort( pairs.begin(), pairs.end(), [=]( std::pair<cvf::Vec2d, unsigned>& a, std::pair<cvf::Vec2d, unsigned>& b ) {
+        return a.second > b.second;
+    } );
+
+    size_t primaryIndex = 0;
+    if ( numOfNeighbouredPoints[pairs[0].first] > numOfNeighbouredPoints[pairs[1].first] )
+    {
+        primaryIndex = 0;
+    }
+    else
+    {
+        primaryIndex = 1;
+    }
+
+    cvf::Vec2d primaryAxisVector   = pairs[primaryIndex].first * axesVectorCandidates[pairs[primaryIndex].first];
+    cvf::Vec2d secondaryAxisVector = pairs[1 - primaryIndex].first * axesVectorCandidates[pairs[1 - primaryIndex].first];
+
+    // Find starting and end point
+    cvf::Vec2d startPoint = to2d( surfacePoints.front() );
+    cvf::Vec2d endPoint   = to2d( surfacePoints.back() );
+
+    size_t row       = 0;
+    size_t column    = 0;
+    size_t maxColumn = ( potentialEndPointIndices.size() + 1 ) / 2;
+    size_t maxRow =
+        ( ( endPoint - secondaryAxisVector * ( maxColumn - 1 ) ) - startPoint ).length() / primaryAxisVector.length() + 1;
+
+    for ( size_t k = 0; k < potentialEndPointIndices.size(); k++ )
+    {
+        if ( column >= 2 )
+            column = ( k - 2 ) / 2;
+        else
+            column = 0;
+
+        // Bring all points on the same line, so that we are having a one-dimensional comparison
+        cvf::Vec2d currentProjected = to2d( surfacePoints[k] ) - secondaryAxisVector * column;
+        cvf::Vec2d endProjected     = endPoint - secondaryAxisVector * ( maxColumn - 1 );
+
+        // Check if pointing in same direction
+        if ( ( currentProjected - startPoint ).dot( endPoint - startPoint ) > 0.0 )
+        {
+            if ( ( currentProjected - startPoint ).length() > ( endProjected - startPoint ).length() )
+            {
+                // Found a new end point
+                endPoint = currentProjected + secondaryAxisVector * ( maxColumn - 1 );
+                maxRow   = ( currentProjected - startPoint ).length() / primaryAxisVector.length() + 1;
+                continue;
+            }
+        }
+        else
+        {
+            // Found new start point
+            maxRow += ( currentProjected - startPoint ).length() / primaryAxisVector.length();
+            startPoint = currentProjected;
+        }
+    }
+
+    // Fill list of index to point data in order to generate the triangles later.
+    cvf::Vec2d currentPoint = startPoint;
+
+    size_t                             index = 0;
+    std::vector<std::vector<unsigned>> indexToPointData;
+    indexToPointData.resize( maxColumn, std::vector<unsigned>( maxRow, -1 ) );
+    for ( column = 0; column < maxColumn; column++ )
+    {
+        currentPoint = startPoint + secondaryAxisVector * column;
+        for ( row = 0; row < maxRow; row++ )
+        {
+            if ( vectorFuzzyCompare( to2d( surfacePoints[index] ), currentPoint, epsilon ) )
+            {
+                indexToPointData[column][row] = static_cast<unsigned>( index );
+                currentPoint                  = to2d( surfacePoints[index] ) += primaryAxisVector;
+                index                         = std::min( index + 1, surfacePoints.size() - 1 );
+            }
+            else
+            {
+                currentPoint += primaryAxisVector;
+            }
+        }
+    }
+
+    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 )
+        {
+            generateTriangleIndices( indexToPointData, iIdx, jIdx, triangleIndices );
+        }
+    }
+
+    return std::make_pair( surfacePoints, triangleIndices );
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+bool RifSurfaceImporter::generateTriangleIndices( const std::vector<std::vector<unsigned>>& indexToPointData,
+                                                  const size_t&                             i,
+                                                  const size_t&                             j,
+                                                  std::vector<unsigned>&                    triangleIndices )
+{
+    unsigned q1 = indexToPointData[i + 0][j + 0];
+    unsigned q2 = indexToPointData[i + 0][j + 1];
+    unsigned q3 = indexToPointData[i + 1][j + 0];
+    unsigned q4 = indexToPointData[i + 1][j + 1];
+
+    if ( q1 != ( (unsigned)-1 ) && q2 != ( (unsigned)-1 ) && q4 != ( (unsigned)-1 ) && q1 != ( (unsigned)-1 ) &&
+         q4 != ( (unsigned)-1 ) && q3 != ( (unsigned)-1 ) )
+    {
+        triangleIndices.push_back( q1 );
+        triangleIndices.push_back( q2 );
+        triangleIndices.push_back( q4 );
+        triangleIndices.push_back( q1 );
+        triangleIndices.push_back( q4 );
+        triangleIndices.push_back( q3 );
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+bool RifSurfaceImporter::vectorFuzzyCompare( const cvf::Vec2d& vector1, const cvf::Vec2d& vector2, double epsilon )
+{
+    auto AlmostEqualRelativeAndAbs = [=]( double A, double B, double maxRelDiff ) -> bool {
+        // Check if the numbers are really close -- needed
+        // when comparing numbers near zero.
+        double diff = fabs( A - B );
+
+        A             = fabs( A );
+        B             = fabs( B );
+        float largest = ( B > A ) ? B : A;
+
+        if ( diff <= largest * maxRelDiff ) return true;
+        return false;
+    };
+    return ( AlmostEqualRelativeAndAbs( vector1.x(), vector2.x(), epsilon ) &&
+             AlmostEqualRelativeAndAbs( vector1.y(), vector2.y(), epsilon ) );
+}

ApplicationCode/FileInterface/RifSurfaceImporter.h

@@ -22,6 +22,7 @@
 
 #include <utility>
 #include <vector>
+#include <limits>
 
 #include <QString>
 
@@ -31,6 +32,15 @@ class RifSurfaceImporter
 {
 public:
     static void readGocadFile( const QString& filename, RigGocadData* gocadData );
-
     static std::pair<std::vector<cvf::Vec3d>, std::vector<unsigned>> readPetrelFile( const QString& filename );
+    static std::pair<std::vector<cvf::Vec3d>, std::vector<unsigned>> readOpenWorksXyzFile( const QString& filename );
+
+private:
+    static bool       generateTriangleIndices( const std::vector<std::vector<unsigned>>& indexToPointData,
+                                               const size_t&                             i,
+                                               const size_t&                             j,
+                                               std::vector<unsigned>&                    triangleIndices );
+    static bool       vectorFuzzyCompare( const cvf::Vec2d& vector1,
+                                          const cvf::Vec2d& vector2,
+                                          double            epsilon = std::numeric_limits<double>::epsilon() );
 };

ApplicationCode/ProjectDataModel/Surfaces/RimFileSurface.cpp

@@ -185,6 +185,10 @@ bool RimFileSurface::loadDataFromFile()
 
         surface = m_gocadData->gocadGeometry();
     }
+    else if ( filePath.endsWith( "dat", Qt::CaseInsensitive ) )
+    {
+        surface = RifSurfaceImporter::readOpenWorksXyzFile( filePath );
+    }
 
     m_vertices       = surface.first;
     m_tringleIndices = surface.second;

ApplicationCode/UnitTests/RifSurfaceImporter-Test.cpp

@@ -204,3 +204,59 @@ TEST( RifSurfaceImporter, ReadClippedPetrelData )
         EXPECT_TRUE( indices[i] != ( (unsigned)-1 ) );
     }
 }
+
+TEST( RifSurfaceImporter, ReadTinyOpenWorksXyzFile )
+{
+    QDir baseFolder( TEST_DATA_DIR );
+
+    QString filename( "RifSurfaceImporter/tiny-test.dat" );
+    QString filePath = baseFolder.absoluteFilePath( filename );
+    EXPECT_TRUE( QFile::exists( filePath ) );
+
+    auto surface = RifSurfaceImporter::readOpenWorksXyzFile( filePath );
+
+    auto vertices = surface.first;
+    auto indices  = surface.second;
+
+    EXPECT_EQ( (size_t)13, vertices.size() );
+    EXPECT_EQ( (size_t)24, indices.size() );
+
+    if ( indices.size() > 0 )
+    {
+        EXPECT_EQ( (size_t)0, indices.front() );
+        EXPECT_EQ( (size_t)11, indices.back() );
+
+        for ( size_t i = 0; i < indices.size(); i++ )
+        {
+            EXPECT_TRUE( indices[i] != ( (unsigned)-1 ) );
+        }
+    }
+}
+
+TEST( RifSurfaceImporter, ReadLargeOpenWorksXyzFile )
+{
+    QDir baseFolder( TEST_DATA_DIR );
+
+    QString filename( "RifSurfaceImporter/norne_xyz.dat" );
+    QString filePath = baseFolder.absoluteFilePath( filename );
+    EXPECT_TRUE( QFile::exists( filePath ) );
+
+    auto surface = RifSurfaceImporter::readOpenWorksXyzFile( filePath );
+
+    auto vertices = surface.first;
+    auto indices  = surface.second;
+
+    EXPECT_EQ( (size_t)60805, vertices.size() );
+    EXPECT_EQ( (size_t)360396, indices.size() );
+
+    if ( indices.size() > 0 )
+    {
+        EXPECT_EQ( (size_t)0, indices.front() );
+        EXPECT_EQ( (size_t)60801, indices.back() );
+
+        for ( size_t i = 0; i < indices.size(); i++ )
+        {
+            EXPECT_TRUE( indices[i] != ( (unsigned)-1 ) );
+        }
+    }
+}

ApplicationCode/UnitTests/TestData/RifSurfaceImporter/tiny-test.dat

@@ -0,0 +1,13 @@
+0.0 0.0 0.0
+1.0 1.0 0.0
+2.0 2.0 0.0
+4.0 4.0 0.0
+2.0 0.0 0.0
+3.0 1.0 0.0
+4.0 2.0 0.0
+5.0 3.0 0.0
+5.0 1.0 0.0
+6.0 2.0 0.0
+6.0 0.0 0.0
+7.0 1.0 0.0
+8.0 2.0 0.0