@ThomasG77 Yeah, they are of course tools doing this perfectly. We prefer however to implement ourselves as so we are looking for some efficient algorithms rather than packages/software. We do it now as we explained to Kirk Kuykendall using dict in Python. It does job but there is efficiency issue as we will apply it on millions of lines.

We are not sure about that either. In our previous attempts for another issue we noticed however Shapely was 10 times slower than our mixed Python-Fortran implementation.

Indeed, we have already implemented a similar idea in Python using defaultdict. This is not efficient as you need to update the dictionary entirely after joining every two segments.

Please note that we don't know which point belongs to which line. Line segments are not listed in any order. They are randomly located in an array. If you know which line is aligned (next) to a line, joining them is not an issue. But notice that we don't have any topology information of lines. No information between points and lines.

@DevdattaTengshe Right, even after splitting the will have points in common. The gaps above (part left) are only for demonstration. We don't have a record of belonging points to a specific line. Indeed, the inputs are segments (or lines for question 2) and points. We drew them here together for demonstration.

@PolyGeo None. We will code in Python. The points above are intersection points with many other lines that are removed later. So our interest is to split each line to many segments based on the intersecting points (part left). Our another interest is reversing the act, that is having a series of aligned (parallel) segments to join them as one line if possible (part right).

Thanks, finally we could understand your point. We added a demonstration of the problem too.

@whuber We added a complete solution that we are looking for as update 2.

@whuber & UffeKousgaard: The correct solution for your data is a rectangle (white) in the middle plus two triangles 3 and 4. These are largest ICHs (with no overlapping) and cover the area fully.

@UffeKousgaard As in the question, after finding the largest ICH it must be removed and the procedure needs to be repeated. We guess whuber's comment refers to this.

Note that this is not a complete solution but just an idea. As mentioned in the question, once the first largest convex found, it must be removed from data somehow and the procedure needs to be repeated as long as required. Even if you implement the complete solution (which is not an easy task, BTW) in addition to some incorrect results (sometimes) the method is not optimum.

The above illustration is the actual result of our implementation of Uffe's idea. Apparently, as we also mentioned there are some problems :(