Mesh rasterisation #2367

rinkk · 2019-02-18T10:07:34Z

Takes a 2D mesh and constructs an asc-raster by sampling elevation values at user-specified equidistant intervals.

endJunction · 2019-02-18T10:12:30Z

Applications/Utils/FileConverter/MeshToRaster.cpp

@@ -0,0 +1,205 @@
+/**
+ * \file moveMeshNodes.cpp


Wrong file name. Just use \file

endJunction · 2019-02-18T10:12:43Z

Applications/Utils/FileConverter/MeshToRaster.cpp

+ * 2012/03/07 KR Initial implementation
+ *
+ * \copyright
+ * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)


endJunction · 2019-02-18T10:13:27Z

Applications/Utils/FileConverter/MeshToRaster.cpp

+
+#include <memory>
+#include <string>
+#include <iostream>


Don't include iostream, use logog.

endJunction · 2019-02-18T10:14:59Z

Applications/Utils/FileConverter/MeshToRaster.cpp

+#include "MeshLib/Elements/Quad.h"
+#include "MeshLib/MeshSearch/MeshElementGrid.h"
+
+/// Returns the index of the element p is located in when projected onto a mesh.


Probably just me, but what is the "element p"?

endJunction · 2019-02-18T10:18:06Z

Applications/Utils/FileConverter/MeshToRaster.cpp

+                {
+                    MeshLib::Node const node(x+x_off[i], y+y_off[i], 0);
+                    std::size_t const elem_idx = getProjectedElement(elems, node);
+                    if (elem_idx != std::numeric_limits<std::size_t>::max())


~~This if condition seems to be always false, because it is in the else branch of the same test in line 169~~
Sorry, I'm wrong here. You tricked me by using the same variable name. Choose a different name, please.
Or, if you extract the "test for all of the pixels' corners" in own function, you can keep it as is; this would also be more readable.

endJunction · 2019-02-18T10:22:32Z

Applications/Utils/FileConverter/MeshToRaster.cpp

+    std::size_t const n_rows = static_cast<std::size_t>(std::ceil((max[1] - min[1]) / cellsize));
+
+    // raster header
+    std::ofstream out(output_arg.getValue(), std::ios::out);


ofstream is already having the default argument for the open mode set to std::ios::out, no need to mention it here.

endJunction · 2019-02-18T10:22:50Z

Applications/Utils/FileConverter/MeshToRaster.cpp

+                                               : mesh->getMinEdgeLength();
+    INFO ("Cellsize set to %f", cellsize);
+
+    std::size_t const nNodes (mesh->getNumberOfNodes());


endJunction · 2019-02-18T10:26:53Z

GeoLib/AnalyticalGeometry.h

+ * Tests if p is right of the line defined by a and b in 2D.
+ * (Note: This functionality is similar to getOrientation(). It is probably
+ * not quite as exact in numerically challenging situations but requires only
+ * about half the runtime.)


I'd recommend to prefer the correctness over runtime in such cases.

This is probably my computer graphics POV but nothing bad happens if this gives a result that is wrong by a measure of epsilon outside of the elevation also being wrong by an epsilon (which is most likely less of an error than any rounding artifacts, measurement errors, etc that happened before already). Should I still remove it and use the slower method?

@TomFischer your decision.

I would take the safe implementation.

endJunction · 2019-02-18T10:28:31Z

Applications/Utils/FileConverter/MeshToRaster.cpp

+    if (mesh->getDimension() != 2)
+    {
+        ERR("The programme requires a mesh containing two-dimensional elements "
+            "(i.e. triangals or quadrilaterals.");


endJunction · 2019-02-18T10:30:04Z

Applications/Utils/FileConverter/MeshToRaster.cpp

+    cmd.parse(argc, argv);
+
+    INFO("Rasterising mesh...");
+    if (!input_arg.isSet())


Isn't tclap already testing this? The third last argument to the constructor of the argument is the requirement, isn't it?

endJunction · 2019-02-18T10:31:48Z

Applications/Utils/FileConverter/MeshToRaster.cpp

+std::size_t getProjectedElement(std::vector<const MeshLib::Element*> const& elems,
+                                MeshLib::Node const& node)
+{
+   // std::vector<MeshLib::Element*> const& elems = mesh.getElements();


drop unused code...

codecov · 2019-02-18T15:32:16Z

Codecov Report

Merging #2367 into master will increase coverage by 0.25%.
The diff coverage is 0%.

@@            Coverage Diff             @@
##           master    #2367      +/-   ##
==========================================
+ Coverage   32.37%   32.63%   +0.25%     
==========================================
  Files         527      527              
  Lines       19935    19771     -164     
  Branches     9413     9260     -153     
==========================================
- Hits         6454     6452       -2     
+ Misses      10179    10032     -147     
+ Partials     3302     3287      -15

Impacted Files	Coverage Δ
MeshLib/MeshSearch/MeshElementGrid.h	`0% <0%> (ø)`	⬆️
GeoLib/Polyline.cpp	`35.86% <0%> (-0.61%)`	⬇️
Applications/ApplicationsLib/ProjectData.cpp	`0% <0%> (ø)`	⬆️

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 6382411...00ba9c1. Read the comment docs.

endJunction · 2019-02-18T15:59:02Z

Is it possible to have a small example to be included in the ctests?

TomFischer · 2019-02-19T07:29:42Z

Applications/Utils/FileConverter/MeshToRaster.cpp

+double getElevation(MeshLib::Element const& elem, MeshLib::Node const& node)
+{
+    MathLib::Vector3 const n = MeshLib::FaceRule::getSurfaceNormal(&elem).getNormalizedVector();
+    MeshLib::Node const orig = *elem.getNode(0);


Why making a copy here?
MeshLib::Node const& orig = *elem.getNode(0);

TomFischer · 2019-02-19T08:42:20Z

Applications/Utils/FileConverter/MeshToRaster.cpp

+    INFO ("Cellsize set to %f", cellsize);
+
+    std::vector<MeshLib::Node*> const& nodes_vec(mesh->getNodes());
+    GeoLib::AABB bounding_box(nodes_vec.begin(), nodes_vec.end());


TomFischer · 2019-02-19T08:59:28Z

Applications/Utils/FileConverter/MeshToRaster.cpp

+#include "MeshLib/MeshSearch/MeshElementGrid.h"
+
+/// Returns the index of the element the given node is located in when projected onto a mesh.
+std::size_t getProjectedElement(std::vector<const MeshLib::Element*> const& elems,


The function name doesn't express exactly what the comment says and what is returned.

TomFischer · 2019-02-19T09:02:41Z

Applications/Utils/FileConverter/MeshToRaster.cpp

+            // centre of the pixel is located within a mesh element
+            if (elem_idx != std::numeric_limits<std::size_t>::max())
+                out << getElevation(*(mesh->getElement(elem_idx)), node) << " ";
+            // centre of the pixel isn't located within a mesh element


From the above comment it is already clear that the centre of the pixel isn't located within a mesh element. Maybe, it is better to write what you'll do in the following (checking the corner points of the element).

There are multiple comments below this line, including one that mentions the corner points.

TomFischer · 2019-02-19T09:11:26Z

You could use GeoLib::Raster and Applications/FileIO/AsciiRasterInterface for the implementation to reuse existing code.

rinkk · 2019-02-19T10:28:15Z

The raster interface is hard to use in this case. For one thing, it'd require writing more code than I did. Most annoyingly, it'd require storing all the data first before I can start creating a raster object (currently I'm not storing anything!), then creating a header object and then writing all of it to a file.

rinkk · 2019-02-19T12:44:26Z

I'm using the slower (but presumably more correct) method now for determining if a point is in a triangle.

For the record: I think it's extremely annoying to do this. The runtime is noticably longer now but the result for the test data is exactly the same. Not a single byte was changed. We're basically hobbling our legs here, nobody was ever going to notice a difference in the results...

TomFischer · 2019-02-20T05:21:00Z

I would suggest to put the asc test file to git large filesystem.

rinkk · 2019-02-20T08:52:32Z

I did the test precisely as both Dima and Lars said I should do it.

endJunction · 2019-02-20T08:59:52Z

@rinkk For discussion an refactored version of the tool https://github.com/endJunction/ogs/tree/master2raster.

TomFischer · 2019-02-20T10:15:54Z

I don't insist to move the raster test data to git large file system.

copy happens element-wise, insert can do insertions in batches.

rinkk · 2019-02-20T13:34:06Z

👍

TomFischer

👍

ogsbot · 2020-06-19T12:58:09Z

OpenGeoSys development has been moved to GitLab.

See this pull request on GitLab.

endJunction reviewed Feb 18, 2019

View reviewed changes

TomFischer reviewed Feb 19, 2019

View reviewed changes

rinkk force-pushed the master2raster branch from b5b707e to 9213277 Compare February 19, 2019 12:18

endJunction mentioned this pull request Feb 19, 2019

Use exact predicates for orient2d #2371

Merged

endJunction force-pushed the master2raster branch from 69c89a5 to 4acd353 Compare February 20, 2019 12:26

rinkk added 3 commits February 20, 2019 13:28

implemented mesh rasterisation

8b26904

do not test for trailing whitespaces in asc-files

490f8dd

added test

126c245

endJunction added 2 commits February 20, 2019 13:28

[App/U] M2R; Use local lambda for is_right_of().

2f8dc47

[MeL] Use vector::insert which is faster than copy

b59dc9f

copy happens element-wise, insert can do insertions in batches.

endJunction force-pushed the master2raster branch from 4acd353 to 5863429 Compare February 20, 2019 12:28

[App/U] Refactor MeshToRaster utility.

00ba9c1

endJunction force-pushed the master2raster branch from 5863429 to 00ba9c1 Compare February 20, 2019 14:46

TomFischer approved these changes Feb 21, 2019

View reviewed changes

TomFischer merged commit 5685fd6 into ufz:master Feb 21, 2019

rinkk mentioned this pull request Feb 28, 2019

Mesh on mesh mapping fix #2390

Merged

rinkk deleted the master2raster branch April 12, 2019 15:30

Mesh rasterisation #2367

Mesh rasterisation #2367

Conversation

rinkk commented Feb 18, 2019

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

endJunction Feb 18, 2019 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

codecov bot commented Feb 18, 2019 • edited Loading

Codecov Report

endJunction commented Feb 18, 2019

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

TomFischer commented Feb 19, 2019

rinkk commented Feb 19, 2019

rinkk commented Feb 19, 2019

TomFischer commented Feb 20, 2019

rinkk commented Feb 20, 2019

endJunction commented Feb 20, 2019

TomFischer commented Feb 20, 2019

rinkk commented Feb 20, 2019

TomFischer left a comment

Choose a reason for hiding this comment

ogsbot commented Jun 19, 2020

endJunction Feb 18, 2019 •

edited

Loading

codecov bot commented Feb 18, 2019 •

edited

Loading