main_angle() function calculates wrong angles in some cases

[kannes]

Describe the bug
The new main_angle() function calculates wrong angles in some cases. In the following screenshot the top row has "linear" rectangles rotated manually at about 10° per step. The two rows below have a duplicated rectangle rotated at exact steps of 10°. The CRS is EPSG:32618.

The blue line points to the North, the red arrow goes from the centroid to a point in the direction of main_angle($geometry), the label shows main_angle($geometry), the line pattern is rotated 90-main_angle degrees (as rotation 0 is in the East).

The arrow line is drawn with this Geometry Generator:

with_variable(
  'centroid',
  centroid($geometry),
  make_line(
    @centroid,
    project(@centroid, 1, radians(main_angle($geometry)))
  )
)

You can see that at some rotations of the rectangles, the calculated main_angle is 90° off.

• For the 5th rectangle in the top row, the angle should be ~40° but is ~130°.
• For the 8th it should be ~70° but is ~160°.
• In the bottom row the angle for the 2nd rectangle should be 70° but is 160°.
• In the bottom row the angle for the last rectangle should be 90° but is 0°.

How to Reproduce
I attached the data and the style to reproduce this.
rectangles.wkt.txt
style.qml.gz

QGIS and OS versions
Master from today:

QGIS version	3.17.0-Master	QGIS code revision	24aa545
Compiled against Qt	5.15.2	Running against Qt	5.15.2
Compiled against GDAL/OGR	3.3.0dev-762a261f32-dirty	Running against GDAL/OGR	3.3.0dev-762a261f32-dirty
Compiled against GEOS	3.8.1-CAPI-1.13.3	Running against GEOS	3.8.1-CAPI-1.13.3
Compiled against SQLite	3.34.0	Running against SQLite	3.34.0
PostgreSQL Client Version	13.1	SpatiaLite Version	4.3.0a
QWT Version	6.1.6	QScintilla2 Version	2.11.6
Compiled against PROJ	6.3.2	Running against PROJ	Rel. 6.3.2, May 1st, 2020
OS Version	Arch Linux

[kannes]

This is caused by

QGIS/src/core/geometry/qgsinternalgeometryengine.cpp

Lines 1508 to 1510 in f327256

	// constrain angle to 0 - 180
	if ( angle > 180.0 )
	angle = std::fmod( angle, 180.0 );

not taking certain cases into account or maybe rather using faulty logic to fix them.

I am interested in trying to fix this myself, so unless someone feels super eager to do this bug, give me some time to find a calm hour. Should be easy(tm) and I can figure out how to add a test for this too. :)

[roya0045]

@kannes I fail to see how constraining the angle between 0 and 180 would cause this behaviour, the issue seems to be higher up in the code when generating the hull and iterating through vertices.

[kannes]

My hunch is that the angles are negative before the modulo and I have not thought about if that makes any sense to be honest. If it could be fixed higher up, that would be perfect of course.

[roya0045]

Any progress on the debug? I'm curious about this one.

[kannes]

Haven't had a moment to look into it yet, my build process was broken for quite a while. Sorry!
If someone else wants to check it out, please go ahead. :)

[roya0045]

I'll look into it as soon as I can fix this cursed msvc c++17 build.

[kannes]

I spoke too soon. Got my setup running and am looking into it now.

[kannes]

It is definitely not in that code snippet I posted, which does make sense... ;)

Polygon ((353534.44156514067435637 378966.25529708701651543, 353537.43783939041895792 378967.3458517276449129, 353537.78222506638849154 378966.39965985936578363, 353534.78595081664388999 378965.30910521873738617, 353534.44156514067435637 378966.25529708701651543))
angle before if: 340.000000
angle changed by if: 160.000000

-> should be 70

Polygon ((353542.63843526312848553 378974.92373468796722591, 353544.95808016933733597 378975.7369054548908025, 353545.27173175086500123 378974.84218527673510835, 353542.95208684465615079 378974.02901450981153175, 353542.63843526312848553 378974.92373468796722591))
angle before if: 340.681462
angle changed by if: 160.681462

-> should be ~70

I think that for "parallel" ("regular"?) geometries like rectangles the if ( currentArea < area ) condition will not work as intended. The geometry will continue to be rotated and the area will not shrink, leading to the correct angle value not being set properly.

With another if for a abs(currentArea - area) < small_epsilon with a break inside I can get it to exit early with the correct value but looking at the code and the results, I have no clue what it is supposed to do and thus what a proper fix without side effects might be. I think it was written by @nyalldawson, could you please add some inline comments or function "docstring" to the code explaining the concept and supposed function of the algorithm?

[roya0045]

I think that for "parallel" ("regular"?) geometries like rectangles the if ( currentArea < area ) condition will not work as intended. The geometry will continue to be rotated and the area will not shrink, leading to the correct angle value not being set properly.

The main issue is that a box will have the same minimal hull if rotated 90 degree as the width and height are inverted.

The only thing that should be done is if currentArea = area pick the one with the biggest height, don't ignore that parameter. That should solve it. No need to uncertainty terms or quitting early as you may only ever get the wide side first.

[kannes]

That might be a fix for the rectangles but are you sure there are no possible geometries that could have the same convex hull bounding box area for more than two subsequent rotations? It might just be a local minima I guess.

[kannes]

Current state with my testing data (the top row is two squares, then one vertex being moved to the top resp. right, first with a tiny step, then visibly), lots of printf-ing because I have no proper setup and probably remaining bugs. Should really be tested with better test cases I think.

main_angle.zip

The equality comparison is absolutely not ready for usage and there are remaining broken tests I just introduced this morning, yay...

I might continue next weekend but cannot promise it so if someone else wants to have a go, enjoy! :)

[roya0045]

The equality comparison is absolutely not ready for usage and there are remaining broken tests I just introduced this morning, yay...
Are you talking about your epsilon approach or the one I suggested?

[kannes]

The abs(a-b) < epsilon I used!

[roya0045]

try preferring height on stalemate and don't quick abort. Your approach will give varying result depending on the vertex order if you stop too soon.

[kannes]

stalewhat? :}

I would still like to understand what the algorithm is trying to do. I think I mostly got it but guessing seems like a waste of time and potentially introduce new bugs.

[roya0045]

The algorithm creates a convex hull of the geometry ( a minimal all encompassing geometry). From this you iterate over the edges by aligning the bounding box on each edge in order to find the best fit. This is do that way since aligning the box on one edge should produce the optimal result in theory.

In your case you hav 4 edged and the bounding box get the same area each time as the width and eight are the only two parameter changing.