Fix of Circle.intersection and other suggestions...

[richardalligier]

Hello,

I have discovered a bug several months ago on Circle.intersection. I fixed it in: https://github.com/richardalligier/geoml/commit/b3c59e341e23a9bdf9916a1db79c8e11c43e356b/src/circle.ml

I am not "git fluent" and the commit containing this fix also include others things so I cannot easily make a pull request so please feel free to make a copy/pasta... sorry...

Just a few additional suggestions: _As the doc is automatically generated, I do not feel like it should be included in the git _I feel like that modeling lines as (y=ax+b or x=c) is prone to numerical issues, so i have changed as (ax+by+c=0) (see https://github.com/richardalligier/geoml/commit/b3c59e341e23a9bdf9916a1db79c8e11c43e356b/src/line.ml ). Likewise, the Line.intersection uses the package lacaml to solve the linear system. I think it will do a better job than me at dealing with numerical issues. However it introduces a dependence... If you are interested in all the things I suggested, you can merge the commit https://github.com/richardalligier/geoml/commit/b3c59e341e23a9bdf9916a1db79c8e11c43e356b

Best, Richard

[ghilesZ]

Hi richard, thanks for your time and suggestions, i will take a look at it next week. Just to know, what was the bug precisely?

[richardalligier]

Do not remember well, but by looking at it right now, I would say that the code of Circle.intersection return the right result only if the circles have same radius.

If you look at the formula you can convince yourself that it is wrong: the barycenter of c1 (center of c) and c2 (center of c') is computed and then the intersections are supposed to be inside the perpendicular at c1_c2 passing through this barycenter. However with the positive weights, this barycenter always lies inside the segment [c1;c2]. As a consequence Circle.intersection cannot return the right answer for some situations where the projections of the correct intersections on the line (c1;c2) are outside [c1;c2]. Such case happens where the radius of one is much larger than the other one as depicted in https://mathalino.com/sites/default/files/images/006-internally-tangent-circles.jpg

You can test it by yourself by checking that the points returned by Circle.intersection are inside both circles. We just have to compute the distance between the alleged intersections and the centers of the circles. Then we can try different radius and see if the results remain correct.

[ghilesZ]

Ok thanks again i'll try that asap (i dont have access to a computer until monday).

As the doc is automatically generated, I do not feel like it should be included in the git

Totally agree.

Line.intersection uses the package lacaml to solve the linear system

I'm not a big fan of dependencies when they can be avoided, especially if it is just for just one function.

I feel like that modeling lines as (y=ax+b or x=c) is prone to numerical issues, so i have changed as (ax+by+c=0)

We though about this one, and we hesitated between the two for a while. The main reasons we avoided the latter are:

• when a = b = 0 the equation defines something different than a line (either the whole plane when c = 0, either the empty set when c <> 0)
• a given line could have several representation (ex 1x+1y+1=0, 2x+2y+2=0 etc) which would make structural equality not equivalent to semantic equality which is a bit of a shame.

To solve these issues we could have used a sanitize/normalize function that only allow for normalized value (only 1x + 1y + 1 = 0 would be allowed for instance and not 2x + 2y + 2 = 0) but we felt that it was more elegant to use a "smaller" type (caridinality wise) with only right values than using a bigger type and forbid some of its values. That way, we kinda use the type system and pattern matching exhaustivity to make sure that the two problem mentioned never occur. This however has other drawbacks: code is arguably larger, some operations might be slightly slower ... but at the end this is a matter of taste.

[richardalligier]

We though about this one, and we hesitated between the two for a while. The main reasons we avoided the latter are:

• when a = b = 0 the equation defines something different than a line (either the whole plane when c = 0, either the empty set when c <> 0)

This problem remains with your choice. In your code, if a = b = 0, then the line will be X(infinity) or X(-infinity) or X(nan) depending on the sign of c. It will be good to add an exception to be raised in such case.

• a given line could have several representation (ex 1x+1y+1=0, 2x+2y+2=0 etc) which would make structural equality not equivalent to semantic equality which is a bit of a shame.

Yes this is a valid concern if you want to use structural equality. As you said, one could sanitize/normalize.

That way, we kinda use the type system and pattern matching exhaustivity to make sure that the two problem mentioned never occur.

As the line is a private type, you can control how a line is created and ensure that only "normalized" lines are used. That way you also solve the two problem. As a drawback, you need to be sure that every function in line.ml produces "normalized" lines.

I just thought of another model that could solve these issues: type t = {theta:float;c:float}. The coefficient a and b are encoded with theta: a = cos(theta) and b=sin(theta). That way, we have structural equality <=> semantic equality, we are sure that we are not in any degenerate case (a=b=0) and that every produced lines are "normalized". A possible drawback is the computation speed as we will have to call cos and sin.

[ghilesZ]

I've implemented a fix for Circle.intersection which reuses some of your code @richardalligier. It's a shame that you don't appear as a contributor though.

Regarding the discussion about the Line.t type, this seems to be an othogonal subject hence, i'm closing this issue. However, feel free to open other ones if you feel the Line.t type must be changed. Otherwise i'm planning on simply adding a check in Line.make to make sure the case a = b = 0 is rejected.