|
BrowserUk has asked for the wisdom of the Perl Monks concerning the following question:
I have a (large) dataset that consists of and X,Y coordinate pair and a real number (0.0 -> 3.0) representing an attribute, call it a height, and I want plot contour lines connecting points of similar heights.
Any thoughts or pointers on how to start?
The data looks like this, but there is not enough (of the 2.4 millions points) here to do anything useful with: -69.282032302755084 40.000000000000014 0
-69.123493781255831 39.908467741935496 -1.4443382565142906e-006
-68.748013135538145 40.911009397420145 0
-68.964955259756593 39.816935483870985 -2.990721348858345e-006
-68.370049396495517 40.298149116372635 -6.3944096502804299e-006
-68.202015544462682 41.814890597403163 0
-67.658222977894795 40.432681542784742 -1.1025156665695599e-005
-68.80641673825734 39.725403225806467 -4.7264009945253546e-006
-68.647878216758102 39.633870967741949 -6.7019561401411285e-006
-68.489339695258849 39.542338709677438 -8.9173860088063258e-006
-67.970337613779577 39.892083135646438 -1.3463345566431798e-005
-67.644134662322315 42.711486110711832 0
-66.168460737542475 40.816251240724036 -7.7940171180040021e-006
-66.403424513471577 39.677023978360658 -3.2568451287380103e-005
-67.417482219894339 39.751269670319942 -2.2905275056155327e-005
-68.330801173759596 39.45080645161292 -1.132337127325882e-005
-68.172262652260343 39.359274193548401 -1.4263293915609116e-005
-68.013724130761105 39.267741935483883 -1.7815964293169953e-005
-67.074467692415936 43.600639717542869 0
-65.320216546227712 42.631856022495803 -2.5121171527202102e-005
-64.452638549546535 41.120575875191086 4.5678927263972832e-005
-65.2313333965173 39.144242118219182 -1.7497211352324817e-005
-66.154489394499308 39.00980515543462 -6.3610401599566315e-005
-66.84644386289331 39.27582098868416 -4.1213274977535648e-005
-67.855185609261852 39.176209677419372 -2.1907142377233984e-005
-67.696647087762614 39.084677419354854 -2.6613602927452167e-005
-67.538108566263361 38.993145161290336 -3.1880809503868505e-005
-66.493113891610946 44.482196494745864 0
-63.455014677037695 42.328563579380834 -6.7791580746250923e-005
-63.743298666018838 43.86849215765006 -0.00012339986352570304
-63.010618308126652 40.142191039199616 0.00028666487947822099
-63.74890107758462 38.983319224147046 0.00022091344000053716
-62.272335538668685 41.301062854252187 9.8592328446381225e-005
...
I also have another set of data that connects the points as triangles. The numbers are indexes into the dataset above, and the indexes are "given in counter-clockwise direction". Whether that means that the first node of each triangle is the right-most I'm not sure, but otherwise I'm not quite sure why that is specified, beyond completeness.
My thoughts so far are:
- First, find the lowest and highest 'heights'; and the x/y bounds.
- The establish the period of the contours: either divide the range into an arbitrary subdivisions; or some equally arbitrary numerical subdivisions 0.0-0.999999, 0.1-0.199999, 0.2->0,29999 etc.
- Then ... I'm kinda blank at the moment.
Note: I'm not actually wanting to plot the mesh and contours at this point; though that's a possibility later. What I need to determine is the angle, relative to the origin, of each of the (short) lines between nodes that make up each contour. That's easy once I have the set of points and their order than make up each contour.
Some contours will be partial -- starting and/or ending at the edge of the data space.
I'm certainly not expecting code; but just clues on how to begin. Thanks.
With the rise and rise of 'Social' network sites: 'Computers are making people easier to use everyday'
Examine what is said, not who speaks -- Silence betokens consent -- Love the truth but pardon error.
In the absence of evidence, opinion is indistinguishable from prejudice.
|
Re: Contour mapping?
by pryrt (Abbot) on Apr 27, 2016 at 19:55 UTC
|
Normally, I just leave my 3d plotting to gnuplot. That said, it's probably not what you wanted. :-) So I googled a bit, and found a couple of articles that may be of help -- the first, as a rather technical paper, and the second that gives the outline for a class lesson the subject
It's at least a place to start for ideas, even if it doesn't give everything you want right away. But if you do get a solution, make sure you at least post an outline, if not all the code; I'd be interested in the end result.
| [reply] |
Re: Contour mapping?
by pryrt (Abbot) on Apr 27, 2016 at 20:07 UTC
|
Alternately, assuming a smoothish/continuous surface
- take the gradient (dz/dx,dz/dy)
- find (x,y) for local minima and maxima (slopes=0)
- work your way up or down your particular valley or hill, and things that are approximately the same height (within your bands, or whatever) will be connected as a contour band; work clockwise around the min/max
- look for intersections with the bands from other minima or maxima, and possibly join bands when necessary
Not sure how well this will work... and I cannot immediately see whether or not your triangle data set will help. Are the triangles all supposed to be in the same contour band, or are those just surface polygons which could be all over, height-wise?
edit: http://www.innovativegis.com/basis/mapanalysis/topic11/topic11.htm gives some interesting surfaces for GIS-based data, so some of the details there might help, since I started thinking about hills and valleys. (and the google:how to distinguish mountain or valley or saddle in 3d surface data that found that, which has some other interesting links)
| [reply] |
|
|
I cannot immediately see whether or not your triangle data set will help. Are the triangles all supposed to be in the same contour band, or are those just surface polygons which could be all over, height-wise?
No. They are just form the mesh via which the "heights" were calculated (FEM). The mesh is generated prior to the heights being calculated, so there's no immediate help there.
That said, if I start at a triangle on the periphery, check each of its node for the highest; then check each of the nodes of all the other triangle that contain that highest point; and so on, eventually it ought to get me to the highest point. But nah, .... I can find that by in a single pass of the nodes data.
So, start at the highest point, follow each of edges of each of the triangles that contain that point, recursively ... to what end? Where am I going?...
Thanks for all the links. Looks like I've a lot of reading to do.
And no, this isn't about plotting; its about calculating and collating all these damned angles.
With the rise and rise of 'Social' network sites: 'Computers are making people easier to use everyday'
Examine what is said, not who speaks -- Silence betokens consent -- Love the truth but pardon error.
In the absence of evidence, opinion is indistinguishable from prejudice.
| [reply] |
|
|
| [reply]
[d/l] |
|
|
And no, this isn't about plotting; its about calculating and collating all these damned angles.
Actually, since plotting isn't your goal, I'm not sure you'll need the angle-sort I suggested... or maybe even the local maxima search-connected. It really depends: do you need to keep two different "islands" separate?
6 g f e d
7 6 c b
7 7 6 7 7 7 a
7 9 7 6 7 7
7 7 6 7 8 7
6 7 7 7 7
7
If you have a grid as pictured, would you want to keep the two rings of height=7 separate, or would they all be the same "contour", in terms of the angles you need to calculate? And for calculating angles, assuming the alphabet ridge were all at the same height, would you need to know that a was adjacent to b, and b adjacent to c, but a not adjacent to c (etc)? Or would you just not care about ordering and adjacency, because you just have to compute so many annoying angles?
- If adjacency doesn't matter, and you just need bands of approximately equal height, whether or not they are part of the same ring, then you don't need the local min/max searches, and you don't need the angle sorting
- If adjacency doesn't matter, but you do need to treat the two rings of 7s separately, then you'll still have to do the min/max search, but not the angle sorting
- If adjacency does matter and the rings need to be separate, then the procedure I outlined before is probably a good starting thought.
| [reply]
[d/l] |
|
|
Re: Contour mapping?
by graff (Chancellor) on Apr 28, 2016 at 01:55 UTC
|
It's odd that you start by saying "an X,Y coordinate pair and a real number (0.0 -> 3.0)", but then most of the points in your sample have negative values in the third column. (Granted, all the negative values are very close to zero, but still...)
I realize the sample is too small to draw any conclusions, but it's worth noting that all the X values are unique, as are all the Y values. Also, if I sort on either X or Y, and look at intervals between the points on the given axis, there's very little repetition (but the little repetition I see appears near the middle of the range from shortest to longest interval). In other words, the distribution of sample points over the XY plane appears to be pretty well randomized, and is likely to be fairly dense.
It might be interesting to look at a sample of the "triangle" data, along with the corresponding rows of XYZ coordinates that it refers to. I'd be curious whether the triangles evince some relevant property, like (maybe): no XY coordinate is contained within any such triangle, or the XY plane is fully covered without overlap by the set of triangles, or the Z values satisfy some condition for each set of three points, or whatever. (There was presumably some motivation for defining the triangles, and one can hope that it would be relevant to the task at hand.)
Apart from that, if your primary goal is "to determine ... the angle, relative to the origin, of each of the ... lines between nodes that make up each contour", I'd be tempted to set "arbitrary subdivisions" of the X and Y ranges as well, to quantize and simplify (at least for a first pass) the process of interpolating endpoints for those lines whose slopes are to be reported.
Or maybe the triangle data will provide some way to simplify that process? E.g. get the XYZ centroids of the triangles (if they don't overlap), or maybe even interpolate the Z value at a uniformly quantized XY point that happens to fall within the triangle. I don't know...
Of course, just getting a visual display of (some portion of) the surface via gnuplot or whatever would be a good start, too. | [reply] |
|
|
I'd be curious whether the triangles evince some relevant property, like (maybe): no XY coordinate is contained within any such triangle, or the XY plane is fully covered without overlap by the set of triangles
The triangles are a non-overlapping mesh (Delaunay or complient Delaunay Trangulation); the data is from a Finite Element Modeling program.
And yes; I'm also pretty sure the triangles are the key. Having played a little with my stratification idea; it is clear that it is going no where -- producing a tangled web of string rather than nice flowing contours.
This is what I now think I need to do:
Starting with a simple triangular mesh of points, with the heights at the nodes:
/4\
/ | \
/ | \
/ | \
/ | \
/ | \
/ | \
/ | \
/ | \
/ | \
/ | \
3--------------------6.5--------------------5
/|\ | /|
+\
/ | \ | / |
+ \
/ | \ | / |
+ \
/ | \ | / |
+ \
/ | \ | / |
+ \
/ | \ | / |
+ \
/ | \ | / |
+ \
/ | \ | . / |
+ \
/ | \ | / |
+ \
/ | \ | / . |
+ \
/ | \|/ |
+ \
2--------------------5.5--------------------9--------------------7.5
+---------------------6
\ | /|\ |
+ /
\ | / | \ |
+ /
\ | / | \ |
+ /
\ | / | . \ |
+ /
\ | / | \ |
+ /
\ | / | \ |
+ /
\ | / | \ |
+ /
\ | / | \ |
+ /
\ | / | \ |
+ /
\ | / | \ |
+ /
\|/ | \|
+/
3--------------------6.5--------------------6
\ | /
\ | /
\ | /
\ | /
\ | /
\ | /
\ | /
\ | /
\ | /
\ | /
\4/
I pick my contour values -- integers 3 to 8 are convenient for this examples -- and process the edges of the triangles linearly (being careful not to process edges twice where they are part of two triangles), and interpolate the positions on those edges where any of my chosen contour values will cross.
Eg. On the edge between value 2 and 5.5, I find the positions are which 3, 4, 5 cross that edge; and store those coordinates in separate arrays for each value. I end up with something as below:
/4\
/ | \
/ | \
/ | \
/ | \
/ 5 \
/ | \
/ | \
/ | \
/ 6 \
/ | \
3-----4------5-----6-6.5-----6--------------5
/|\ | /|
+\
/ | 4 | / |
+ \
/ | \ 7 6 |
+ \
/ | 5 | / |
+ \
/ 4 \ | / 6
+ \
/ | 6 | 7 |
+ \
/ | \ 8 / |
+ \
/ | 7 | . / |
+ \
/ 5 \ | 8 |
+ \
/ | 8 | / . 7
+ \
/ | \|/ |
+ \
2-----3------4----5 -5.5--6-----7-----8-----9-------------8------7.5
+------7--------------6
\ | /|\ |
+ /
\ | 8 | \ |
+ /
\ 5 / | \ |
+ /
\ | 7 | . \ 7
+ /
\ | / 8 8 |
+ /
\ | 6 | \ |
+ /
\ 4 / | \ |
+ /
\ | 5 | \ |
+ /
\ | / 7 7 |
+ /
\ | 4 | \ |
+ /
\|/ | \|
+/
3-----4------5-----6-6.5--------------------6
\ | /
\ | /
\ 6 /
\ | /
\ | /
\ | 5
\ 5 /
\ | /
\ | /
\ | /
\4/
Then it just a case of ordering those arrays of points to produce the polylines that make up the contours something like this:
/4\
/ *| \
/ * | \
/ * | \
/ * | \
/ * 5 \
/ * * | * \
/ * * | * \
/ * * | * \
/ * * 6* * \
/ * * * | * * \
3-----4------5-----6-6.5-----6--------------5
/|\ * * * | * /|\
/ *| 4 * * | * / |
+ \
/ * | * \ * * *7 * * 6 |
+ \
/ * |* 5 * * | * / ** |
+ \
/ * 4 * \ * * | * / 6
+ \
/ * *| * 6 * | 7* |
+ * \
/ * * | * * \ * 8 / * |
+ * \
/ * * |* * 7 *| * . / * |
+ * \
/ * * 5 * * \ * | 8 * |
+ * \
/ * * * | * * 8 | / * . 7
+* * \
/ * * * | * * * \|/ * |
+ * * \
2-----3------4----5 -5.5--6-----7-----8-----9-------------8------7.5
+------7--------------6
\ * * * | * * * /|\ * |
+ * *
\ * * * | * * 8 | \ * |
+ * *
\ * * 5 * * / * | \ * |
+* *
\ * * |* * 7 *| . \ * 7
+ *
\ * * | * * / * 8********8 *|
+ *
\ * *| * *6 * | \ * |
+ *
\ * 4 * / * * | \ * |
+ *
\ * |* 5 * * | \ * |
+ *
\ * | * / * * 7****************7 |
+ *
\ *| 4 * * | \ |
+ *
\|/ * * * | \|*
+
3-----4------5-----6-6.5--------6.25--------6
\ * * * | *** /
\ * * *| *** /
\ * * 6 *** /
\ * * | /
\ * * | /
\ * *| * *5
\ * 5 * * /
\ * | /
\ * | /
\ *| /
\4/
I don't have a full handle on the ordering process; but I've a few ideas to try. (Actually, I don't need to order them to do the intersections with the boundaries; though it might be nice to draw the contours to see that they are correct before going on.)
Then all I have to do :) is work out which of those contours intersect the polylines that form the boundaries and work out where and at what angle they intersect.
And then I can start on trying to collate the results in a way that is meaningful...the driving force behind the whole problem.
With the rise and rise of 'Social' network sites: 'Computers are making people easier to use everyday'
Examine what is said, not who speaks -- Silence betokens consent -- Love the truth but pardon error.
In the absence of evidence, opinion is indistinguishable from prejudice.
| [reply]
[d/l]
[select] |
|
|
Wow.I'd say don't worry about ordering (until later). Take one triangle at a time, and look at its three edges in whatever sequence is given.
- If one edge crosses a given whole-number Z value, then exactly one other edge in that triangle will cross that same value, so you'll get exactly one iso-contour line segment per whole-number Z per triangle, and its position and angle in the overall XY coordinate space can be worked out.
- While a single vertex (XY point) can be shared among a variable number of triangles (2 to ...8? 9?), a given vertex pair (i.e. an edge) will appear in at most two triangles. (Edges appearing in only one triangle define the outer boundary of the XY plane being charted.)
- Putting those two things together, as you march through each triangle, watch for the case where you come to an edge you've seen before (of course, its end points will be presented in opposite order, relative to the first time you saw this edge). If you've already seen it, you already know if it crosses any whole-number Z value(s), and you have the coordinates of the end points for that iso-contour line in the adjacent (previously done) triangle.
- When you find the other edge(s) of the current triangle crossing the same whole-number Z value(s) -- these also may be edges you've seen before -- you have what you need to maintain a linked list of the nodes for the given iso-contour line, in addition to knowing the position and angle of each line segment.
For 2.4 million vertices (how many triangles?), it's a lot to keep track of, but it seems like just a matter of looping over triangles, and for each of those, looping over edges, keeping track of the edges you've done, and keeping a properly indexed list of iso-contour line segments as they come up, so that you know when you encounter each interior edge the second time, and can traverse the connections of linked contour lines.
Regarding that first point above, the logical possibilities are:
- None of the three edges contain a whole-number Z value
- One or more of the vertices is itself a whole-number Z value
- An edge contains just one whole-number Z value, so that value appears on one other edge (or may be the at the opposite vertex)
- An edge contains two or more whole-number Z values; either all those values will also appear on one other edge, or else they'll be divvied up among the other two edges (or one of them may be at the opposite vertex)
I hope you have a generous schedule (and are getting paid enough) for doing this. Sounds like a lot of work.
(updated wording in a couple bullet points)
| [reply] |
|
|
Re: Contour mapping?
by Anonymous Monk on Apr 27, 2016 at 21:43 UTC
|
This ought to take forever... Pick a target height, say 1.0000, then go through all triangles that have points on both sides of the target. Plot a point in the middle of the triangle. After all triangles have been examined and the points plotted, zoom way back, and there is your contour.
| [reply] |
Re: Contour mapping?
by Anonymous Monk on Apr 28, 2016 at 14:17 UTC
|
You have triangles. Each triangle determines a plane in three space. Find the intersection of that plane with a horizontal (constant height) plane, this will be a line. The direction of that line is the angle you seek.
| [reply] |
|
|
You have triangles. Each triangle determines a plane in three space. Find the intersection of that plane with a horizontal (constant height) plane, this will be a line. The direction of that line is the angle you seek.
Hm. If that is true; and I can't honestly dismiss it on first (nor second) reading, it is a brilliant observation and would save huge amounts of work. I would only need to consider those triangle that border the boundaries I'm interested in.
What's puzzling me is that I wasn't aware that I had described the angles I'm seeking in sufficient detail for anyone else to understand; but despite that, you've hit the nail on the head and saved me a huge amount of processing in the bargain.
Any pair of points of equal height, on any two of the three sides will define the angle of all pairs of equal height points across the triangle. so if I follow my boundary polylines, I only need consider those triangles on either side of them, and finding the angles is simple geometry.
I simply do not know how to thank you enough.
(And there are people who want to ban Anonymonk....)
With the rise and rise of 'Social' network sites: 'Computers are making people easier to use everyday'
Examine what is said, not who speaks -- Silence betokens consent -- Love the truth but pardon error.
In the absence of evidence, opinion is indistinguishable from prejudice.
| [reply] |
|
|
(P1 - P0) crossproduct (P2 - P0) crossproduct (0, 0, 1)
Where P0 P1 P2 are points of the triangle,
and (0, 0, 1) is the vertical vector
| [reply]
[d/l] |
|
|
|
|
| |
|
|
|
| [reply] |
|
|
| [reply] |
|
|
