Feature/get turn info add retrieve policy

[barendgehrels]

This PR (actually my first) adds a retrieve policy, and passes it to get_turn_info.

This policy retrieves a point known as "k" (pk / qk), the first non-duplicate point coming consecutively after the points "i" and "j" (two points from two segments p and q, part of the intersection of two segments, so these points are pi,pj and qi,qj).

In the past points pk/qk were just passed. But now we need another point (pl/ql) and it is not convenient to add that. Therefore the policy is introduced.

An additional advantage is that point "k" is now retrieved lazily, on demand. Often that point is not necessary. It is actually only necessary if there are collinearies. And often that is not the case. So in all those cases, performance should increase just a little bit.

This PR adds the policy and retrieves point "k" using the method .get(). In another PR I will add an argument to get, to be able to retrieve also the subsequente point, as described above. But for now this PR is large enough.

It is also related to side calculation, the place where that point is used, to check whether it is located left/right or collinear to a segment.

To my surprise, also non-typical usage of side calculation were introduced in the implementations for linear intersections, and the typical order is changed there. Therefore I copied the originals there - they cannot use the policy because they don't use the "k" point. Actually there should be looked at, because it is looking suspicious, but we can do that afterwards.

I'm curious to your findings.

[awulkiew]

Interesting. Regarding future plans, what for will be additional points pl and ql needed?

The PR looks good to me. I'm only wondering why you've made certain naming decisions:

• The name of the policy. Would it be possible to rename it to know what is retrieved exactly? What is it for on a higher level of abstraction? What is stored and why? AFAIU it's for optional, additional data of an intersection segment: next point and flags for the current segment (is first or last segment). Or are the flags for first/last point? Is that correct? OptionalPointDataPolicy, AdditionalSidePointPolicy, OptionalSegmentInfoPolicyetc.?
• The name of function get() may be confusing because it's not clear what is it returning exactly. Furthermore we use get() for coordinates. Why not get_point() or even explicitly get_point_k() (and later get_point_q())?

[barendgehrels]

Right, it is ready again, all tests pass (now in debug mode - to verify assertions)

[awulkiew]

It is probably me who wrote it. I don't remember exactly why turns are checked this way, I'd have to analyse the whole file and maybe debug. However I may point you in a right direction right now. This is probably for turns where linear geometries are involved, correct? With linear geometries the turns have slightly different meaning than in areal/areal case. For instance the equality doesn't mean that the whole segment was the same but that the intersection point was equal to one of the points of the geometry. Turns hold additional data indicating whether prior to the intersection point the segments were collinear or are crossing at the intersection point, etc. AFAIR this had something to do with the fact that linear geometries are not ordered like areal and that from various geometrically equal sets of segments I needed the same set of turns, no matter the direction of the segments. E.g. for areal/areal one relative direction could result in i/i operations but with reverse segments the operations would be x/x. And this is true for all other operations. But linear turns would be the same.

So AFAIR parts of this code are there to transform ordered areal turns into unordered linear turns. Sometimes additional info had to be calculated. I didn't check this for sure but this might have been such case.

[awulkiew]

Furthermore, this is the code calculating turns for endpoints, i.e. the first or the last point of a linestring. There is no such thing in areal geometries. However in this case we have to know how the other geometry behaves at the endpoint, when one geometry ends or begins. And since the geometry either begins or ends at intersection point maybe this is why "non-segment" points from other geometries are considered.

[awulkiew]

So unless we know what's happening exactly we should probably be careful about changing this code. But if all of the tests pass with this change it should be ok I think.

[awulkiew]

Out of curiosity, why did you decide to change the names of variables to less meaningful? Above I see a comment referencing specific points so this condition looks unrelated now.

[awulkiew]

Looks good, I'm only worried about the replacement of side_calculator_for_endpoint with manually calculated sides. AFAIU the whole point of it was to have this code used several times encapsulated as a tool. Why have you decided to remove it and replace with code looking less clear?

[barendgehrels]

IMO it is not less clear. The code was using named sides, but substituting them with completely other sides, even from different sources. I was really surprised by that and consider it as wrong.

In this new version you can a bit better see what is happening there, which is still not 100% clear, but IMO more clear than just call encapsulated magic in the wrong way.

I had to do it in really small steps to avoid breaking anything, it took me quite some time but all is equivalent as it was before, and all tests pass of course.

[barendgehrels]

@awulkiew About your question, related to the new method of using things as range_q.at(0):

Is this new way of expressing points more readable to you (e.g. instead of i(), j(), k())?

I was quite used to i,j and k - but I noticed that @vissarion kind of objected to this way of working, and I think I agree with him. So having range-like methods as at(0) (I also considered just subscripting) make it more clear which point is meant, then just alphabetical order i,j,k,l (l, so L, would have been new, but would have been really confusing with 1).

So to answer your question directly, yes, I think that at(0) etc is more readable to me than .i()

[barendgehrels]

About your question

Out of curiosity, why did you decide to change the names of variables to less meaningful?

Because the names were wrong. It used names like "pk" - while, as you mentioned, there is indeed no "pk" in end-points. It is another point, and which one depends on the situation. It is called 5 or 6 times in different ways, I could not come up with a common name, also because I don't understand all the tiny details of this code.

[barendgehrels]

Finally, about your comment:

So AFAIR parts of this code are there to transform ordered areal turns into unordered linear turns.

I did not realize they were unordered, and I don't think they are. Is that really the case? But for sure, I understand that on the end-points there is different behaviour, and different decisions should be made. What I don't understand (at least not understood) was that side calculator was used, a points were inserted for the k-points, and all looked like there were general calculations happening, but behind the screens there were sides calculated totally different from the normal cases (and called like the normal cases). That was my incentive to change it, besides, of course, the need to pass a range instead of points (which I could not do here) and to access the ranges differently, and pass them differently, as was coming from review remarks...

Also, this code was called from both linear/linear and linear/areal.

Are there bulk tests for these cases, inside or outside (at MySQL side) our codebase? I'm curious if it just always works...

Of course I'm willing to give more meaningful names if I know them, and to understand the code better, and refactor some more things if necessary - but for now I thought this was already quite a step and, as a step itself, consistent.

[awulkiew]

We have the tests here:
https://github.com/boostorg/geometry/blob/develop/test/algorithms/overlay/get_turns_linear_linear.cpp
https://github.com/boostorg/geometry/blob/develop/test/algorithms/overlay/get_turns_linear_areal.cpp

So debugging them, stepping inside and seeing in which scenarios which sides are calculated here and which points are taken into account should give us the answer. I'll check it out too. Btw, if you're using Visual Studio you could use the extension I wrote to easily see the geometries, I don't know if you know about it (https://marketplace.visualstudio.com/items?itemName=AdamWulkiewicz.GraphicalDebugging).

[awulkiew]

Ok, so this is the code handling a case when the intersection point is at the first or last point of the linestring (geometry P) and the second point of a segment of the other linestring or ring (geometry Q). So as usual we want to check which segment has to be followed when intersection or union is calculated. Because there is no 3rd point of the first geometry (and I probably wanted to use existing utility for this) we have to generate some point and it's sufficient to take the first point of Q.

E.g. in this case: "LINESTRING(1 0,1 1)", "LINESTRING(0 0,1 0,2 0)" (https://github.com/boostorg/geometry/blob/develop/test/algorithms/overlay/get_turns_linear_linear.cpp#L40) these are the initial segments (P - 1 red segment, Q - 2 green segments):

and these are the new 2 blue segments for P (let's call it P') with the added first point from Q and 2 segments of Q (green):

So now we can check the sides "the regular way" using the existing tools, i.e. side_calculator to check sides of P' and Q and IntersectionInfo to check spikes of Q.

With that said in L/L case (2 linear geometries) 2 sides checks would be enough: Pi or Pj WRT Q1 and Q2 (so side_pk_q2 and side_pk_p (P') because this is the same as side_pk_q in the original code) followed by spike check to see if the side is opposite. There are additional conditions in operations_of_equal because this function is used also explicitly in the L/A (linear and areal) case (e.g. here https://github.com/boostorg/geometry/blob/develop/include/boost/geometry/algorithms/detail/overlay/get_turn_info_la.hpp#L780), so if the segment of linear P (or second segment of P') is on the left of Q, it's outside Q which means union on P, and if it's right of Q this means intersection on P. Because the areal geometries are ordered clockwise. And by left of Q I mean P' k (which is Pj for first and Pi for last segment of a linestring) turning more left than Qk WRT Q1 (in operations_of_equal named _p, so of P'). So yes:

• the naming is messed up,
• for L/L this additional side calculation could be removed which implies different implementation for L/A.

[awulkiew]

I edited the comment to make it more clear what I'm talking about ;)

I don't suggest that this code should be refactored/cleaned as a part of this PR. I just wanted to help understanding it.

Anyway, I agree with merging it together with possible change of names of variables.

So this would probably be the key:

• first1
  • side_pk_q2 -> side_pj_q2
  • side_pk_p -> side_pj_q
  • side_qk_p -> side_qk_q
• last1
  • side_pk_q2 -> side_pi_q2
  • side_pk_p -> side_pi_q
  • side_qk_p -> side_qk_q

[awulkiew]

One more clarification, by "unordered" I meant that linear geometries are not CW nor CCW and that there is no interior or exterior on left or right side of linear geometry. That regardless of the relative direction of this geometry the methods and operations are the same. They allow to simply follow the linear geometry and generate intersection or union or stop at x. Previously when the regular/areal turns were used for linear geometries different turns were generated based on the side of the geometries, like i/i vs x/x etc. Maybe "unordered" was not a good name, my point was that the operations and methods may mean something different than you're used to. E.g. operation x means "the end of linestring" and method e may be generated for the first endpoint of linestring because it means "intersection point equal to points of both geometries", not "the segments before arriving at IP are equal" as in A/A case, etc.

[barendgehrels]

Thanks a lot for your clarifications! I'll continue later.

[barendgehrels]

Thanks again for your comments and clear pictures.

Sure, I know your debugging tools, certainly. But I'm always using QtCreator now - and not Visual Studio. I remember that you did have a tool for that too, but (at that time) it did slow things down, sometimes giving issues, and besides that I usually have the need to add more and more debug information, so therefore I rely on SVG creation and do not visualize geometries themselves in the IDE.

It now becomes clear to me why a point of Q was substituted for P. But I still think it was mainly a trick to be able to use the side calculator. That need is gone now, because (by applying ranges) it became clear that it was messing names up (I don't mean implementation, but just the names, as you confirmed).
If that substitution is not necessary anymore, it is trivial to calculate just the point of p (pj in this case) w.r.t. q1 or q2.

For which you provide the key, as far as I can see, and I could use those names (in the calling functions, so i/o x. But note that they are called as, for example the last one:
int const side_qk_x = sides.apply(range2.at(0), range1.at(1), range2.at(2));
If this is side qk_q - why is the second point then coming from range 1? I tried to make it consistent with your naming (so use range2.at(1)) but... that gives me errors. So it is not side_qk_q but the side of qk w.r.t. a segment made from qi to pj.

Or maybe I did not understand you completely...

[awulkiew]

I'm playing with the code in develop and replacing the second parameter of side_calc at line 432 (first1) from ri1 to rj2 and at line 485 (last1) from rj1 to rj2 doesn't change the results for me. At least all of the set operations tests pass. So in case of the first segment (first1) the P' is changed from ri2, ri1, rj1 to ri2, rj2, rj1, or in the new notation from rng2.at(0), rng1.at(0), rng1.at(1) to rng2.at(0), rng2.at(1), rng1.at(1)

Anyway, I also tried to remove side_calculator at all and at line 428 put this:

auto side = inters.get_side_strategy();

int const side_pj_q2 = side.apply(rj2, rk2, rj1);
int const side_pj_q1 = side.apply(ri2, rj2, rj1);
int const side_qk_q1 = side.apply(ri2, rj2, rk2);

std::pair<operation_type, operation_type>
    operations = operations_of_equal(side_pj_q2, side_pj_q1, side_qk_q1);

and at line 482:

auto side = inters.get_side_strategy();

int const side_pi_q2 = side.apply(rj2, rk2, ri1);
int const side_pi_q1 = side.apply(ri2, rj2, ri1);
int const side_qk_q1 = side.apply(ri2, rj2, rk2);

std::pair<operation_type, operation_type>
    operations = operations_of_equal(side_pi_q2, side_pi_q1, side_qk_q1);

of course in operations_of_equal take (int side_px_q2, int side_px_q1, int side_qk_q1), where px means that it's either pi or pj so the endpoint of the segment P which is not the intersection point. After this change get_turns_linear_linear test passes as well.

[barendgehrels]

Right, I think I tried the same thing (but with errors...)
I propose to merge it as it is now, and I don't object, of course, if you change the variable names and/or the substituted points afterwards.

[awulkiew]

Overlay tests also pass for me. In get_turns_la I adapted operations_of_equal to new arguments, e.g. like that:

std::pair<operation_type, operation_type>
    operations = get_info_e::operations_of_equal(side_calc.pk_wrt_q2(), side_calc.pk_wrt_p1(), side_calc.qk_wrt_p1());

I also changed the names of points passed to side_calculator when it was possible. Note that at line 793 the intersection point is in the middle of Q1 and is equal to pi which means that a new segment qi->pi has to be considered so in this case sides are calculated WRT this new segment and not Q1. And later the same for last, the intersection point may be in the middle. But this is because side_calculator is used, if it were not used and sides calculated manually we could probably instead simply use Q1 because it has the same direction as the segment from qi to the intersection point at the Q1.

[awulkiew]

Do I understand correctly that you do not get the errors with the code as it is now, correct? Only after making changes? And you suggest I play with it after the merging? If that's the case I agree.

[barendgehrels]

Sure, that is what I suggested. So I did not get the errors as it is now.
Thanks for merging! And thanks again for reviewing

[awulkiew]

Follow-up: #540
It worked for me out of the box. Maybe you had problems because of some typo?
Anyway, let's move the discussion there.

[barendgehrels]

All right - I was tired, probably not a typo but tried one of them, not the pair together. So obviously that didn't work... Thanks, will go to that new PR