08-30-2025, 11:59 AM
(08-30-2025, 09:33 AM)Torkel M. Jodalen Wrote: Anyone who cares to post the final iteration of the "clockwise sort of coordinates" source code, which was recently discussed in the PB forums? I waited too long to store a local copy.
Thanks.
Here's my final including the Moments of Inertia:
'
Code:
'Calculate Area, Centroid and Moments of Inertia
' of a polygon with clockwise ordered vertices
#COMPILE EXE
#DIM ALL
'
'MACRO RadToDeg(r) = (r * 57.2957795130823209##)
MACRO Pi = 3.14159265358979324##
MACRO TwoPi = 6.283185307179586##
TYPE CartPt
s AS LONG 'sequence number
x AS LONG ' x coord
y AS LONG ' y coord
a AS EXT 'angle in degrees or radians about centre
d AS EXT 'distance from centre
END TYPE
GLOBAL TotPoints AS LONG
GLOBAL Minx AS EXT
GLOBAL Maxx AS EXT
GLOBAL Miny AS EXT
GLOBAL MaxY AS EXT
GLOBAL gsCoords AS STRING
FUNCTION PBMAIN() AS LONG
LOCAL hGrf,hFont AS DWORD
LOCAL x AS LONG
LOCAL CentreX, CentreY,Ixx,Iyy,Ixy AS EXT
LOCAL totarea AS EXT
LOCAL pts() AS CartPt
FONT NEW "ARIAL",12,1 TO hFont
'Set parameters
totpoints = 6
minx = 0
miny = 0
maxx = 20
maxy = 20
totarea = (maxx-minx) * (maxy-miny)
'Initialise a random set of points
DIM pts(1 TO TotPoints)
RANDOMIZE TIMER
FOR x = 1 TO TotPoints
pts(x).s = x
pts(x).x = RND() * (Maxx - minx) + minx
pts(x).Y = RND() * (Maxy - miny) + miny
NEXT
'Sort points clockwise to from an ordered polygon
BuildPolygon pts()
'Calculate the Centroid (centre of mass) of the sorted polygon
PolygonCentroidInt Pts(),CentreX,CentreY
PolygonInertia(pts(), CentreX ,CentreY , Ixx, Iyy, Ixy )
'Display it
GRAPHIC WINDOW NEW "Points - Hit any key to cycle and exit", 50,50,600,600 TO hGrf
GRAPHIC SET FONT hFont
' GRAPHIC WINDOW STABILIZE
GRAPHIC SCALE (minx,miny) - (maxx,maxy)
PlotPoints(Pts())
GRAPHIC WAITKEY$
GRAPHIC SET POS (minX,MinY)
DrawPoly(Pts())
GRAPHIC WAITKEY$
GRAPHIC PAINT BORDER(CentreX,CentreY) , %RGB_GREENYELLOW, %RGB_BLACK
GRAPHIC WAITKEY$
Radiate(Pts(),CentreX,CentreY)
GRAPHIC COLOR %RGB_BLACK,-2
DrawCentre(CentreX,CentreY)
GRAPHIC WAITKEY$
GRAPHIC SET POS (minX,MinY)
GRAPHIC COLOR %RGB_BLACK,-2
GRAPHIC PRINT "Range: " & CHR$("(",STR$(minx),STR$(minY),") - (",STR$(maxX),",",STR$(maxY),")")
GRAPHIC PRINT CHR$("Total Area: ",STR$(totarea))
GRAPHIC PRINT CHR$("Polygon Area:",STR$(ShoelaceArea(pts())))
GRAPHIC PRINT CHR$("Percent coverage:",FORMAT$(ShoelaceArea(pts())/totarea*100,"0.0"))
GRAPHIC PRINT "Boundary Lattice (Integer) Points: " & STR$(BoundaryPoints(pts()))
GRAPHIC PRINT "Interior Lattice (Integer) Points: " & STR$((ShoelaceArea(pts()) -BoundaryPoints(pts())/2)+1)
GRAPHIC PRINT CHR$("Centroid: ",FORMAT$(Centrex,"0.0")," , ",FORMAT$(centreY," 0.0"))
GRAPHIC PRINT CHR$("Moments of Inertia: Ix ",FORMAT$(Ixx,"0.0")," ,Iy ",FORMAT$(Iyy," 0.0")," ,Ixy ",FORMAT$(Ixy," 0.0"))
GRAPHIC WAITKEY$
gsCoords &= "Centerpoint: " & FORMAT$(Centrex,"0.0") & "," & FORMAT$(Centrey,"0.0") & $LF
gsCoords &= CHR$("Moments of Inertia: Ix ",FORMAT$(Ixx,"0.0")," ,Iy ",FORMAT$(Iyy," 0.0")," ,Ixy ",FORMAT$(Ixy," 0.0"))
? gsCoords,,"Polygon"
GRAPHIC WAITKEY$
GRAPHIC WINDOW END
END FUNCTION
'==== Drawing Functions =====
FUNCTION Radiate(pts() AS cartPt,cx AS EXT, cy AS EXT )AS LONG
LOCAL x AS LONG
FOR x = LBOUND(pts()) TO UBOUND(pts())
GRAPHIC LINE (cx,cy) - (pts(x).x,pts(x).y),%RGB_WHITE
NEXT
END FUNCTION
FUNCTION PlotPoints(pts() AS CartPt) AS LONG
LOCAL x AS LONG
LOCAL xsize,ysize AS EXT
xsize = (maxx-minx)/200
ysize = (maxy-miny)/200
FOR x = LBOUND(pts()) TO UBOUND(pts())
GRAPHIC ELLIPSE (pts(x).x -xsize , pts(x).y-ysize) - (pts(x).x +xsize , pts(x).y +ysize),%RGB_RED,%RGB_RED
NEXT
END FUNCTION
FUNCTION DrawPoly(pts() AS CartPt) AS LONG
LOCAL x AS LONG
GRAPHIC SET POS (pts(LBOUND(pts())).x,pts(LBOUND(pts())).y)
FOR x = LBOUND(pts()) + 1 TO UBOUND(pts())
GRAPHIC LINE (pts(x-1).x,pts(x-1).y)- (pts(x).x,pts(x).y)
NEXT
GRAPHIC LINE (pts(UBOUND(pts())).x,pts(UBOUND(pts())).y) - (pts(LBOUND(pts())).x,pts(LBOUND(pts())).y)
END FUNCTION
FUNCTION DrawCentre(CentreX AS EXT,CentreY AS EXT) AS LONG
LOCAL tX,tY AS EXT
GRAPHIC TEXT SIZE "X" TO tX, tY
GRAPHIC SET POS (CentreX- tx/2,CentreY-ty/2)
GRAPHIC PRINT "X"
END FUNCTION
'==== Polygon calculation Functions ===============
FUNCTION BuildPolygon(pts() AS CartPt) AS LONG
'Create a convex polygon from an array of unsorted points (vertices)
LOCAL cx,cy,OffsetX,OffsetY AS EXT
LOCAL SUmx,SumY,x,Totpoints AS LONG
Totpoints = UBOUND(pts)
FOR x = 1 TO TotPoints
SumX += pts(x).x
SumY += pts(x).y
NEXT
CX = SumX/TotPoints
CY = SumY/TotPoints
'Calculate polar coordinate from arithmentic Mean Centre point
FOR x = 1 TO Totpoints
OffsetX = pts(x).x - CX
Offsety = pts(x).y - CY
pts(x).a = ATN(Offsety / Offsetx) ' if you want to use degrees: RadToDeg(ATN(Offsety / Offsetx))
IF Offsetx < 0 THEN 'Quadrant 2 or 3
pts(x).a += Pi '180 if RadToDeg
ELSEIF Offsety < 0 THEN ' Quadrant 4
pts(x).a += TwoPi '360 if RadToDeg
END IF
pts(x).d = SQR(Offsety ^ 2 + offsetx ^ 2)
NEXT
ARRAY SORT pts(), CALL SortClockwise
FOR x = 1 TO totpoints
gsCoords &=FORMAT$(pts(x).x,"0.0") & "," & FORMAT$(pts(x).y,"0.0") & $LF
NEXT
END FUNCTION
FUNCTION SortClockwise(p1 AS Cartpt,p2 AS Cartpt) AS LONG
'Used by CartPt Array Sort
IF p1.a < p2.a THEN FUNCTION = -1 : EXIT FUNCTION
IF p1.a > p2.a THEN FUNCTION = 1 : EXIT FUNCTION
'same angle - sort "inside" one first - arbitrary decision
IF p1.d < p2.d THEN FUNCTION = 1 : EXIT FUNCTION
IF p1.d > p2.d THEN FUNCTION = -1 : EXIT FUNCTION
END FUNCTION
FUNCTION PolygonCentroidInt(pts() AS CartPt, BYREF cx AS EXT, BYREF cy AS EXT) AS LONG
LOCAL i,j, n AS LONG
LOCAL A2 AS QUAD ' Twice the area (integer)
LOCAL sumX AS QUAD, sumY AS QUAD
LOCAL factor AS QUAD
FOR i = 1 TO UBOUND(pts)
j = IIF(i = UBOUND(pts), 1, i + 1)
factor =(pts(i).x * pts(j).y) - (pts(j).x * pts(i).y)
A2 = A2 + factor
sumX += (pts(i).x + pts(j).x) * factor
sumY += (pts(i).y + pts(j).y) * factor
NEXT
IF A2 = 0 THEN EXIT FUNCTION ' degenerate polygon
' Centroid = (sumX / (3*A2), sumY / (3*A2))
' Note: A2 = 2*Area, so denominator = 6*Area
cx = sumX / (3.0 * A2)
cy = sumY / (3.0 * A2)
FUNCTION = -1 ' success
END FUNCTION
'Pick's vTheorem
FUNCTION GCD(BYVAL a AS LONG, BYVAL b AS LONG) AS LONG
'Greatest Common Divisor of end points
'To understand the GCD's role, consider a line segment connecting two lattice points.
'The slope of this line segment can be represented as a fraction, say a/b.
'The GCD of a and b (written as gcd(a, b)) helps determine how many lattice points lie on that line segment.
'Specifically, the number of lattice points on the segment (excluding the endpoints) is gcd(a, b) - 1
'In essence, the GCD helps determine the spacing between lattice points along a line segment with a given slope,
'which is essential for accurately counting boundary points in Pick's theorem
LOCAL t AS LONG
a = ABS(a)
b = ABS(b)
DO WHILE b <> 0
t = b
b = a MOD b
a = t
LOOP
FUNCTION = a
END FUNCTION
FUNCTION ShoelaceArea(pts() AS CartPt) AS EXT
'The "shoelace formula," also known AS Gauss's area formula or the surveyor's formula,
'Is a mathematical algorithm used to calculate the area of a POLYGON given the Cartesian
'coordinates of its vertices.
'It's called the shoelace formula due to the cross-multiplication pattern of coordinates,
'resembling the lacing of shoes
LOCAL i,j AS LONG
LOCAL Area AS EXT
FOR i = 1 TO UBOUND(pts)
j = i + 1
IF j > UBOUND(pts) THEN j = 1
area = area + (pts(i).x * pts(j).y - pts(j).x * pts(i).y)
NEXT
FUNCTION = ABS(area) / 2
END FUNCTION
FUNCTION BoundaryPoints(pts() AS CartPt) AS LONG
'integer points that lie on the perimeter
LOCAL b,i,j AS LONG
FOR i = 1 TO UBOUND(pts)
j = i + 1
IF j > UBOUND(pts) THEN j = 1
B = B + GCD(pts(j).x - pts(i).x, pts(j).y - pts(i).y)
NEXT
FUNCTION = B
END FUNCTION
FUNCTION PolygonInertia (pts() AS cartpt, _
BYVAL XC AS EXT, BYVAL YC AS EXT, _
IXc AS EXT, IYc AS EXT, IXYc AS EXT) AS DOUBLE
LOCAL i AS LONG, j, n AS LONG
LOCAL a, Atotal AS EXT
LOCAL Ix, Iy, Ixy AS EXT
N = UBOUND(pts())
FOR i = 1 TO N
j = IIF(i = N, 1, i + 1)
a = pts(i).x * pts(j).y - pts(j).x * pts(i).y ' cross product term
Atotal = Atotal + a
Ix = Ix + a * (pts(i).y^2 + pts(i).y * pts(j).y + pts(j).y^2)
Iy = Iy + a * (pts(i).x^2 + pts(i).x * pts(j).x + pts(j).x^2)
Ixy = Ixy + a * (pts(i).x * pts(j).y + 2 * pts(i).x * pts(i).y + 2 * pts(j).x * pts(j).y +pts(j).x * pts(i).y)
NEXT
' Raw (about origin)
Ix = Ix / 12
Iy = Iy / 12
Ixy = Ixy / 24
' Apply parallel axis theorem to shift to centroid
IxC = Ix - Atotal / 2 * YC^2
IyC = Iy - Atotal / 2 * XC^2
IxyC = Ixy - Atotal / 2 * XC * YC
FUNCTION = Atotal / 2 ' Return signed area (can use ABS if needed)
END FUNCTION
'