XProfan

WindowTitle "Ritter's Fast Bounding-Ball Algorithm "+\
"(Errechnet eine ziemlich kleine Punktwolken-umfassende Kugel, hier in 2D,"+\
" im Schnitt 5% bis 20% größer als Minimalsphäre, aber sehr schnell!"
' <TransCoded to XProfan 11, (CL) CopyLeft 2014-07 by P. Specht@gmx.at>
' <Source: http:'geomalgorithms.com/a08-_containers.html#fastBall%28%29>
' Theory: https://en.wikipedia.org/wiki/Bounding_sphere
' =====================================================================
' Copyright 2001 softSurfer, 2012 Dan Sunday
' This code may be freely used and modified for any purpose providing that this
' copyright notice is included with it. SoftSurfer makes no warranty for this code,
' and cannot be held liable for any real or imagined damage resulting from its use.
' Users of this code must verify correctness for their application.
' <The Program is> based on the algorithm given by [Jack Ritter, 1990].
' ====================================================================================
' <ASSUMPTIONS ON AVAILABLE VECTOR MATH OVERRIDDEN BY SIMPLE XPROFAN MATH>
' Result: Ball (2D: Circle) B with Center and Radius <Remarks <..> by P.Specht>
'=====================================================================================
' TESTSYSTEM mit n& Punkten:
var n&=12:dec n&
'=====================================================================================
DECLARE BCenterX!,BCentery!,BRadius!' Übergabevariablen, global definiert
Declare Px![n&],Py![n&],i&
WindowStyle 24:Font 2:Randomize:Window 0,0-%maxx,%maxy-40
var xh&=width(%hwnd)/2:var yh&=height(%hwnd)/2

REPEAT

    Px![]=xh&*0.75+rnd(xh&/2):Py![]=yh&*1.5-rnd(xh&/2)
    CLS:usepen 0,5,255:whileloop 0,n&:i&=&Loop:line Px![i&],Py![i&] - Px![i&]+1,Py![i&]:endwhile
    FastBall Px![],Py![],n&
    Usepen 0,2,rgb(0,0,255):usebrush 0,rgb(0,255,0)
    Ellipse BCenterX!-BRadius!,BCenterY!+BRadius! - BCenterX!+BRadius!,BCenterY!-BRadius!
    waitinput 10000

until %key = 27

END

proc FastBall

    parameters Px![],Py![],n&
    ' (based on the algorithm given by [Jack Ritter, 1990])
    declare Cx!,Cy!,rad!,rad2!' Center of ball, radius and radius squared
    declare xmin!, xmax!, ymin!, ymax!' bounding box extremes
    declare Pxmin&,Pxmax&,Pymin&,Pymax&' index of  P[] at box extreme
    declare i&,dPxx!,dPxy!,dPyx!,dPyy!,dx2!,dy2!,dist!,dist2!,dPx!,dPy!
    ' Find a large diameter to start with
    xmin!=10^50:xmax!=-1*10^50:ymin!=10^50:ymax!=-1*10^50:Pxmin&=0:Pxmax&=0:Pymin&=0:Pymax&=0

    whileloop 0,n&:i&=&Loop

        if Px![i&]<xmin!: xmin!=Px![i&]:Pxmin&=i&:endif

            if Px![i&]>xmax!: xmax!=Px![i&]:Pxmax&=i&:endif

                if Py![i&]<ymin!: ymin!=Py![i&]:Pymin&=i&:endif

                    if Py![i&]>ymax!: ymax!=Py![i&]:Pymax&=i&:endif

                    endwhile

                    :::::::Usepen 0,1,rgb(0,0,255):usebrush 0,rgb(0,255,0)
                    :::::::rectangle Px![Pxmin&],Py![Pymin&] - Px![Pxmax&],Py![Pymax&]
                    cx!=(Px![Pxmin&]+Px![Pxmax&])/2
                    cy!=(Py![Pymin&]+Py![Pymax&])/2
                    rad2!=sqr(cx!-Px![Pxmin&])
                    case rad2!<sqr(cy!-Py![Pymin&]):rad2!=sqr(cy!-Py![Pymin&])
                    rad!=sqrt(rad2!)
                    :::::::usepen 0,9,0:line cx!,cy! - cx!+1,cy!
                    ' Now check that all points P[i] are inside the Ball
                    ' and if not, expand the ball just enough to include them
                    ' Vector dP; float dist, dist2;

                    whileloop 0,n&:i&=&Loop

                        dPx!=Px![i&]-Cx!
                        dPy!=Py![i&]-Cy!
                        dist2!=sqr(dPx!)+sqr(dPy!)'= norm2(dP);
                        case dist2!<=rad2!:continue' P[i] is inside the ball already
                        ' P[i] not in ball, so expand ball  to include it:
                        dist!=sqrt(dist2!)
                        rad!=(rad!+dist!)/2' enlarge radius just enough
                        rad2!=sqr(rad!)
                        Cx!=Cx!+((dist!-rad!)/dist!)*dPx!' shift Center toward P[i]
                        Cy!=Cy!+((dist!-rad!)/dist!)*dPy!

                    Endwhile

                    Bcenterx!=Cx!:Bcentery!=Cy!:BRadius!=rad!

                Endproc