Fitting rectangle with maximum area inside polygon in PostGIS

There are many ways to get the desired result, the main thing is to understand what steps will lead you to the right result!

The main thing is all - the Idea, Method, Process, Result and Emotions!

In general, the Idea is in your fantasy (representation) and is as follows::

1. We are looking for the most important point that will help us build a inscribed rectangle into a polygon, in my example it is a table tblf;
2. Next we draw rectangular lines from this point is tblh and tbli;
3. Then we find 2 points on the borders of the polygon;
4. Next, we find the center of the rectangle;
5. Then we create a triangle on 3 points, expand the second triangle and unite them...

I've come all the way through, which is implemented as Postgre/PostGIS SQL code below:

WITH
    tbla(geom) AS (SELECT ST_BuildArea('POLYGON((529004.994 159096.691,529005.519 159096.992,529011.013 159100.145,529021.878 159106.109,529022.1 159105.95,529033.4 159097.85,529066.85 159073.7,529059.1 159060.25,529005.5 159096.35,529004.994 159096.691))')),
    tblb AS (SELECT ST_SETSrid(geom,3857) geom FROM tbla),
    tblc AS (SELECT geom1, geom2, CASE WHEN geom1>geom2 THEN geom1 WHEN geom1<=geom2 THEN geom2 END FROM (SELECT ST_MakeLine((ST_PointN(ST_Boundary(ST_OrientedEnvelope(geom)),1)),
    (ST_PointN(ST_Boundary(ST_OrientedEnvelope(geom)),2))) AS geom1, ST_MakeLine((ST_PointN(ST_Boundary(ST_OrientedEnvelope(geom)),2)), (ST_PointN(ST_Boundary(ST_OrientedEnvelope(geom)),3))) AS geom2 FROM tblb) AS foo),
    tbld AS (SELECT (ST_Dump(ST_Intersection(a.geom1, b.geom2))).geom FROM tblc a JOIN tblc b ON ST_Intersects(a.geom1, b.geom2)),
    tble AS (SELECT ST_ShortestLine(a.geom, ST_ExteriorRing(b.geom)) geom FROM tbld a, tblb b),
    tblf AS (SELECT ST_Intersection(a.geom, b.geom) geom FROM tble a JOIN tblb b ON ST_Intersects(a.geom, b.geom)),
    tblg AS (SELECT ST_MakeLine(ST_SetSrid(ST_MakePoint(0,0),3857), ST_Centroid(a.geom)) geom FROM tblf a),
    tblh AS (SELECT ST_Rotate(a.geom, pi()/5+Radians(ST_Azimuth(ST_EndPoint(b.geom1), ST_StartPoint(b.geom1))), ST_Centroid(c.geom)) geom FROM tblg a, tblc b, tblf c),
    tbli AS (SELECT ST_Rotate(a.geom,-pi()/2, ST_Centroid(b.geom)) geom FROM tblh a, tblf b),
    tblj AS (SELECT ST_Intersection(ST_ExteriorRing(a.geom), b.geom) geom FROM tblb a JOIN tblh b ON ST_Intersects(a.geom, b.geom)),
    tblk AS (SELECT ST_Intersection(ST_ExteriorRing(a.geom), b.geom) geom FROM tblb a JOIN tbli b ON ST_Intersects(a.geom, b.geom)),
    tbll AS (SELECT ((ST_Dump(ST_Difference(a.geom, b.geom))).geom) geom FROM tblj a JOIN tble b ON NOT ST_Disjoint(a.geom, ST_Buffer(b.geom,1))),                                                                                                                        
    tblm AS (SELECT((ST_Dump(ST_Difference(a.geom, b.geom))).geom) geom FROM tblk a JOIN tble b ON NOT ST_Disjoint(a.geom, ST_Buffer(b.geom,1))),
    tbln AS (SELECT ST_MakeLine(a.geom, b.geom) geom FROM tbll a, tblm b),
    tblo AS (SELECT ST_Centroid(geom) geom FROM (SELECT (a.geom) geom FROM tbln a, tbln b WHERE ST_Length(a.geom)>ST_Length(b.geom)) foo),                                                                                 
    tblp AS (SELECT geom FROM tblf UNION SELECT geom FROM tblj UNION SELECT geom FROM tblk),
    tblq AS (SELECT ((ST_Dump(ST_DelaunayTriangles(ST_Collect(geom)))).geom) geom FROM tblp),
    tblr AS (SELECT ST_Rotate(a.geom, -pi(), b.geom) geom FROM tblq a, tblo b),
    tbls AS (SELECT ST_Union(geom) geom FROM (SELECT geom FROM tblq UNION SELECT geom FROM tblr) foo)
    SELECT (ST_Dump(ST_Intersection(a.geom, b.geom))).geom FROM tbls a JOIN tblb b ON ST_Intersects(a.geom, b.geom)

So far I have solved this problem in this way, because other functions such as ST_MinimumBoundingCircle() behave approximately in my opinion, too, i.e. do not give 2(3) common points at all...

Maybe I or somebody will finalize it or write my ?...

EDIT 2

WITH
    tbla(geom) AS (SELECT ST_BuildArea('POLYGON((529004.994 159096.691,529005.519 159096.992,529011.013 159100.145,529021.878 159106.109,529022.1 159105.95,529033.4 159097.85,529066.85 159073.7,529059.1 159060.25,529005.5 159096.35,529004.994 159096.691))')),
    tblb AS (SELECT ST_SETSrid(geom,3857) geom FROM tbla),
    tblc AS (SELECT geom1, geom2, CASE WHEN geom1>geom2 THEN geom1 WHEN geom1<=geom2 THEN geom2 END FROM (SELECT ST_MakeLine((ST_PointN(ST_Boundary(ST_OrientedEnvelope(geom)),1)),
    (ST_PointN(ST_Boundary(ST_OrientedEnvelope(geom)),2))) AS geom1, ST_MakeLine((ST_PointN(ST_Boundary(ST_OrientedEnvelope(geom)),2)), (ST_PointN(ST_Boundary(ST_OrientedEnvelope(geom)),3))) AS geom2 FROM tblb) AS foo),
    tbld AS (SELECT (ST_Dump(ST_Intersection(a.geom1, b.geom2))).geom FROM tblc a JOIN tblc b ON ST_Intersects(a.geom1, b.geom2)),
    tble AS (SELECT ST_ShortestLine(a.geom, ST_ExteriorRing(b.geom)) geom FROM tbld a, tblb b),
    tblf AS (SELECT ST_Intersection(a.geom, b.geom) geom FROM tble a JOIN tblb b ON ST_Intersects(a.geom, b.geom)),
    tblg AS (SELECT ST_Rotate(a.geom1, -pi(), ST_Centroid(b.geom)) geom FROM tblc a, tblb b),
    tblh AS (SELECT ST_ShortestLine(a.geom, b.geom) geom FROM tblf a, tblg b),
    tbli AS (SELECT ST_Intersection(ST_ExteriorRing(a.geom), b.geom) geom FROM tblb a JOIN tblh b ON ST_Intersects(a.geom, b.geom)),
    tblj AS (SELECT ST_LongestLine(a.geom, b.geom) geom FROM tblf a, tblg b),
    tblq AS (SELECT ST_Intersection(ST_ExteriorRing(a.geom), b.geom) geom FROM tblb a JOIN tblj b ON ST_Intersects(a.geom, b.geom)),
    tbll AS (SELECT ST_Rotate(a.geom2, -pi(), ST_Centroid(b.geom)) geom FROM tblc a, tblb b),
    tblm AS (SELECT ST_ShortestLine(a.geom, b.geom) geom FROM tblf a, tbll b),
    tbln AS (SELECT ST_Intersection(ST_ExteriorRing(a.geom), b.geom) geom FROM tblb a JOIN tblm b ON ST_Intersects(a.geom, b.geom)),
    tblo AS (SELECT geom FROM tblf UNION SELECT geom FROM tbli UNION SELECT geom FROM tblq UNION SELECT geom FROM tbln)
    (SELECT ST_Union(geom) geom FROM (SELECT ((ST_Dump(ST_DelaunayTriangles(ST_Collect(geom)))).geom) geom FROM tblo) foo)

EDIT 3

Now the script has started to acquire the behavior of a new experimental function that will work for figures like yours, but it will take a little thought and effort to make it universal?...

WITH
    tbla AS (SELECT (ST_Dump(geom)).geom geom FROM polygon),
    tblb AS (SELECT geom1, geom2, CASE WHEN geom1>geom2 THEN geom1 WHEN geom1<=geom2 THEN geom2 END FROM (SELECT ST_MakeLine((ST_PointN(ST_Boundary(ST_OrientedEnvelope(geom)),1)),
    (ST_PointN(ST_Boundary(ST_OrientedEnvelope(geom)),2))) AS geom1, ST_MakeLine((ST_PointN(ST_Boundary(ST_OrientedEnvelope(geom)),2)), (ST_PointN(ST_Boundary(ST_OrientedEnvelope(geom)),3))) AS geom2 FROM tbla) AS foo),
    tblc AS (SELECT (ST_Dump(ST_Intersection(a.geom1, b.geom2))).geom FROM tblb a JOIN tblb b ON ST_Intersects(a.geom1, b.geom2)),
    tbld AS (SELECT ST_ShortestLine(a.geom, ST_ExteriorRing(b.geom)) geom FROM tblc a, tbla b),
    tble AS (SELECT ST_Intersection(a.geom, b.geom) geom FROM tbld a JOIN tbla b ON ST_Intersects(a.geom, b.geom)),
    tblf AS (SELECT ST_Rotate(a.geom1, -pi(), ST_Centroid(b.geom)) geom FROM tblb a, tbla b),
    tblg AS (SELECT ST_ShortestLine(a.geom, b.geom) geom FROM tble a, tblf b),
    tblh AS (SELECT (ST_Dump(ST_Intersection(ST_ExteriorRing(a.geom), b.geom))).geom geom FROM tbla a JOIN tblg b ON ST_Intersects(a.geom, b.geom)),
    tbli AS (SELECT ST_Rotate(a.geom2, -pi(), ST_Centroid(b.geom)) geom FROM tblb a, tbla b),
    tblj AS (SELECT DISTINCT ST_ShortestLine(a.geom, b.geom) geom FROM tbli a JOIN LATERAL (SELECT (geom) AS geom FROM tblh) AS b ON true),
    tblk AS (SELECT ST_Intersection(ST_ExteriorRing(a.geom), b.geom) geom FROM tbla a JOIN tblj b ON ST_Intersects(a.geom, b.geom)),                                                                                                               
    tbll AS (SELECT (ST_Dump(geom)).geom geom FROM tblh UNION SELECT (ST_Dump(geom)).geom geom FROM tblk)
    (SELECT ST_Union(geom) geom FROM (SELECT ((ST_Dump(ST_DelaunayTriangles(ST_Collect(geom)))).geom) geom FROM tbll) foo)

EDIT 4

Here, as an option, may look a new custom function that tries to fit a rectangle as much as possible in the wrong trapezoidal rectangle, circle, triangle, correct polygon:

CREATE OR REPLACE FUNCTION ST_MaximumAreaInscribedRectangleInPolygon(
geom GEOMETRY
)
RETURNS GEOMETRY AS
$BODY$
WITH
tbl_rigth AS (WITH 
    tbla AS (SELECT (ST_Dump(geom)).geom geom),
    tblb AS (SELECT geom1, geom2, CASE WHEN geom1>geom2 THEN geom1 WHEN geom1<=geom2 THEN geom2 END FROM (SELECT ST_MakeLine((ST_PointN(ST_Boundary(ST_OrientedEnvelope(geom)),1)),
    (ST_PointN(ST_Boundary(ST_OrientedEnvelope(geom)),2))) AS geom1, ST_MakeLine((ST_PointN(ST_Boundary(ST_OrientedEnvelope(geom)),2)), (ST_PointN(ST_Boundary(ST_OrientedEnvelope(geom)),3))) AS geom2 FROM tbla) AS foo),
    tblc AS (SELECT ST_OffsetCurve(geom2, -1) geom FROM tblb UNION SELECT ST_OffsetCurve(ST_Rotate(a.geom2, -pi(), ST_Centroid(b.geom)), -1) geom FROM tblb a, tbla b),
    tbld AS (SELECT ST_OffsetCurve(geom1, -1) geom FROM tblb UNION SELECT ST_OffsetCurve(ST_Rotate(a.geom1, -pi(), ST_Centroid(b.geom)), -1) geom FROM tblb a, tbla b),
    tble AS (SELECT (ST_DumpPoints(ST_OffsetCurve(geom2, -1))).geom geom FROM tblb),
    tblf AS (SELECT ST_ShortestLine(a.geom, ST_ExteriorRing(b.geom)) geom FROM tble a, tbla b),
    tblg AS (SELECT (a.geom) geom FROM tblf a, tblf b WHERE ST_Length(a.geom)>ST_Length(b.geom)),
    tblh AS (SELECT ((ST_Dump(ST_Intersection(a.geom, b.geom))).geom) geom FROM tblg a JOIN tbla b ON ST_Intersects(a.geom, b.geom)),
    tbli AS (SELECT ST_ShortestLine(a.geom, b.geom) geom FROM tblh a JOIN LATERAL (SELECT (geom) AS geom FROM tbld) AS b ON true),
    tblj AS (SELECT ((ST_Dump(ST_Intersection(ST_ExteriorRing(a.geom), b.geom))).geom) geom FROM tbla a JOIN tbli b ON ST_Intersects(a.geom, b.geom)),
    tblk AS (SELECT ((ST_Dump(ST_Difference(a.geom, b.geom))).geom) geom FROM tblj a JOIN tblh b ON ST_Disjoint(a.geom, ST_Buffer(b.geom,0.1))),
    tbll AS (SELECT ST_ShortestLine((ST_Dump(a.geom)).geom, b.geom) geom FROM tblk a, tblc b),
    tblm AS (SELECT ((ST_Dump(ST_Intersection(ST_ExteriorRing(a.geom), b.geom))).geom) geom FROM tbla a JOIN tbll b ON ST_Intersects(a.geom, b.geom)),
    tbln AS (SELECT ((ST_Dump(ST_Difference(a.geom, b.geom))).geom) geom FROM tblm a JOIN tblk b ON ST_Disjoint(a.geom, ST_Buffer(b.geom,0.1))),
    tblo AS (SELECT ST_ShortestLine(a.geom, b.geom) geom FROM tbln a JOIN LATERAL (SELECT (geom) AS geom FROM tbld) AS b ON true),
    tblp AS (SELECT ST_ShortestLine((ST_Dump(a.geom)).geom, b.geom) geom FROM tblh a, tblc b),
    tblq AS (SELECT ((ST_Dump(ST_Intersection(a.geom, b.geom))).geom) geom FROM tblo a JOIN tblp b ON ST_Intersects(a.geom, b.geom)),                             
    tblr AS (SELECT (geom) geom FROM tblh UNION SELECT (geom) geom FROM tblk 
            UNION SELECT (geom) geom FROM tbln UNION SELECT (geom) geom FROM tblq)
    SELECT ST_Union(geom) geom FROM (SELECT ((ST_Dump(ST_DelaunayTriangles(ST_Collect(geom)))).geom) geom FROM tblr) foo),
tbl_lefth AS (WITH 
    tbla AS (SELECT (ST_Dump(geom)).geom geom),
    tblb AS (SELECT geom1, geom2, CASE WHEN geom1>geom2 THEN geom1 WHEN geom1<=geom2 THEN geom2 END FROM (SELECT ST_MakeLine((ST_PointN(ST_Boundary(ST_OrientedEnvelope(geom)),1)),
    (ST_PointN(ST_Boundary(ST_OrientedEnvelope(geom)),2))) AS geom1, ST_MakeLine((ST_PointN(ST_Boundary(ST_OrientedEnvelope(geom)),2)), (ST_PointN(ST_Boundary(ST_OrientedEnvelope(geom)),3))) AS geom2 FROM tbla) AS foo),
    tblc AS (SELECT ST_OffsetCurve(geom2, -1) geom FROM tblb UNION SELECT ST_OffsetCurve(ST_Rotate(a.geom2, -pi(), ST_Centroid(b.geom)), -1) geom FROM tblb a, tbla b),
    tbld AS (SELECT ST_OffsetCurve(geom1, -1) geom FROM tblb UNION SELECT ST_OffsetCurve(ST_Rotate(a.geom1, -pi(), ST_Centroid(b.geom)), -1) geom FROM tblb a, tbla b),
    tble AS (SELECT (ST_DumpPoints(ST_OffsetCurve(geom2, -1))).geom geom FROM tblb),
    tblf AS (SELECT ST_ShortestLine(a.geom, ST_ExteriorRing(b.geom)) geom FROM tble a, tbla b),
    tblg AS (SELECT (a.geom) geom FROM tblf a, tblf b WHERE ST_Length(a.geom)<ST_Length(b.geom)),
    tblh AS (SELECT ((ST_Dump(ST_Intersection(a.geom, b.geom))).geom) geom FROM tblg a JOIN tbla b ON ST_Intersects(a.geom, b.geom)),
    tbli AS (SELECT ST_ShortestLine(a.geom, b.geom) geom FROM tblh a JOIN LATERAL (SELECT (geom) AS geom FROM tbld) AS b ON true),
    tblj AS (SELECT ((ST_Dump(ST_Intersection(ST_ExteriorRing(a.geom), b.geom))).geom) geom FROM tbla a JOIN tbli b ON ST_Intersects(a.geom, b.geom)),
    tblk AS (SELECT ((ST_Dump(ST_Difference(a.geom, b.geom))).geom) geom FROM tblj a JOIN tblh b ON ST_Disjoint(a.geom, ST_Buffer(b.geom,0.1))),
    tbll AS (SELECT ST_ShortestLine((ST_Dump(a.geom)).geom, b.geom) geom FROM tblk a, tblc b),
    tblm AS (SELECT ((ST_Dump(ST_Intersection(ST_ExteriorRing(a.geom), b.geom))).geom) geom FROM tbla a JOIN tbll b ON ST_Intersects(a.geom, b.geom)),
    tbln AS (SELECT ((ST_Dump(ST_Difference(a.geom, b.geom))).geom) geom FROM tblm a JOIN tblk b ON ST_Disjoint(a.geom, ST_Buffer(b.geom,0.1))),
    tblo AS (SELECT ST_ShortestLine(a.geom, b.geom) geom FROM tbln a JOIN LATERAL (SELECT (geom) AS geom FROM tbld) AS b ON true),
    tblp AS (SELECT ST_ShortestLine((ST_Dump(a.geom)).geom, b.geom) geom FROM tblh a, tblc b),
    tblq AS (SELECT ((ST_Dump(ST_Intersection(a.geom, b.geom))).geom) geom FROM tblo a JOIN tblp b ON ST_Intersects(a.geom, b.geom)),                             
    tblr AS (SELECT (geom) geom FROM tblh UNION SELECT (geom) geom FROM tblk 
            UNION SELECT (geom) geom FROM tbln UNION SELECT (geom) geom FROM tblq)
    SELECT ST_Union(geom) geom FROM (SELECT ((ST_Dump(ST_DelaunayTriangles(ST_Collect(geom)))).geom) geom FROM tblr) foo)
    SELECT geom FROM tbl_rigth UNION SELECT geom FROM tbl_lefth
    $BODY$
LANGUAGE SQL

Run

SELECT ST_MaximumAreaInscribedRectangleInPolygon(geom) geom FROM <name_table>

Use the function for your solutions if necessary, but don't forget the author himself ?...

Do not forget to adjust the outer boundary lines with the ST_OffsetCurve() function, based on the shape of the original polygon...

Of course, remember the unexpected result but sometimes a slight turn of a figure can bring the expected result in order?...- it's automation ?...

Original solutions,

Now you know that my style is **IMPRE**?

Good luck in learning ...

Translated with www.DeepL.com/Translator (free version)

Well, it looks like a really complicated question. I don't think that there is an optimum algorithm, but maybe you can find one that fits your particular needs. For exemple, if your polygons can be concave, it will be way more difficult (the blue rectangle would be the answer?):

So you should look at the specificities (only convexe polygon?) of your problem, and the degree of precision that you can accept.

If you have only convexe polygon, I think (I may be wrong) that you can test the sides one after another, and try to build the biggest rectangle possible using this side as a side of your rectangle (for exemple by projecting perpendiculars until you cross another side). Not sure if it always works, but it should give you an answer at least. But you still need to handle the case when you have an alternance between angles >90° and <90°. You also probably can try to make rectangular triangle with all the <90° angles before, maybe using the smallest side of this angle, like that:

And from there use the new points (the intersection between the dotted lines and the sides) as if they were points of your polygons, and do a I said before.

For concave polygons, I think it is way more difficult. Maybe you will need to cut the lines of your polygons to make intermediates points (in the first exemple above, the blue rectangle is between points that does not exists, so you first need to cut the sides of your polygons to add points).