'3D skull from 2D Images - By David J. Baty 4/2019
'requires image file "skull.png"
SCREEN _NEWIMAGE(1600, 1000, 32): _FULLSCREEN: _DONTBLEND: _MOUSESHOW "CROSSHAIR"

'** change THESE values to use your own custom image *****************
tnf = 72 '                LINE INPUT "TOTAL NUMBER OF IMAGE FRAMES: ", tnf$: tnf = VAL(tnf$)
fpr = 10 '                LINE INPUT "TOTAL NUMBER OF FRAMES PER ROW: ", fpr$: fpr = VAL(fpr$)
frame_width = 799 '       LINE INPUT "PIXEL WIDTH OF EACH FRAME: ", frame_width$: frame_width = VAL(frame_width$)
frame_height = 660 '      LINE INPUT "PIXEL HEIGHT OF EACH FRAME: ", frame_height$: frame_height = VAL(frame_height$)
tnd = 180 '               LINE INPUT "TOTAL NUMBER OF DEGREES TO SCAN PER FRAME(1-180): ", tnd$: tnd = VAL(tnd$)
imagefile$ = "Skull.png" 'LINE INPUT "Enter Image File Name: ", imagefile$

CLS: PRINT "Loading Image: " + imagefile$
skull& = _LOADIMAGE(imagefile$) 'THIS SKULL image is a multi-image .jpg with 72 frames that are 799 x 660 pixels each
totrec = (tnf * tnd) + tnd
'*********************************************************************

DIM depth(totrec, 1) 'HOLDS CALCULATED 3D DEPTH VALUES(C2) (SORTED BACK TO FRONT)
DIM head#(totrec, 3) 'HOLDS LATITUDE, LONGITUDE, RADIUS
DIM puthead(totrec, 2) 'HOLDS 2D CALCULATED SCREEN COORDINATES

pointincriment# = (3.14159265359 / tnd) * (tnd / 180)
xcenter = 800: ycenter = 500 'MODELS X/Y CENTER on the PHYSICAL screen
xrot = 0: yrot = 320: zrot = 0 'Startup Rotation Angles
posize = .75 ' Percent of Size - SCALE THE MODEL
zobs = 200000 ' Z-OBSERVATION VALUE - use a low value(200 for example) for z-observation exaggeration!
cf# = 6.28318530718 / 360 'Pi*2 divided by 360 degrees
xpointsize = 9: ypointsize = 4 'this defines the SIZE of the filled boxs rendered at each point/ LARGER VALUES 'FILL-IN' THE MODEL
xpointer = 0: ypointer = 0: newframe& = _NEWIMAGE(frame_width, frame_height)
totrec = 0

FOR frame = 1 TO tnf
    _FREEIMAGE newframe&: newframe& = _NEWIMAGE(frame_width, frame_height)
    _PUTIMAGE , skull&, newframe&, (xpointer * frame_width, ypointer * frame_height)-((xpointer * frame_width) + frame_width - 1, (ypointer * frame_height) + frame_height - 1)
    xpointer = xpointer + 1: IF xpointer >= fpr THEN xpointer = 0: ypointer = ypointer + 1
    _DEST 0: _PUTIMAGE (0, 0), newframe&: LOCATE 1, 1: PRINT "Analyzing Frame"; frame; "of"; tnf
    _SOURCE newframe&



    FOR degrees# = 0 TO tnd * .017444444 STEP pointincriment#
        FOR size = 300 TO 1 STEP -.01 '300 = search radius extending from center of image frame
            lookx = INT(COS(degrees# - 1.570796326795) * size + (frame_width * .5))
            looky = INT(SIN(degrees# - 1.570796326795) * size + (frame_height * .5))
            WHILE POINT(lookx, looky) <> _RGB32(0, 0, 0)
                totrec = totrec + 1
                head#(totrec, 1) = (360 / tnf) * frame '         Latitude
                head#(totrec, 2) = degrees# / .0174532916666667 'Longitude
                head#(totrec, 3) = size '                        Radius of THIS point (distance out)
                CIRCLE (lookx, looky), 1, _RGB32(0, 255, 0)
                GOTO skip1
            WEND
        NEXT size
        skip1:
    NEXT degrees#
NEXT frame

_DEST 0
colorgradiant = totrec / 255: colorscale = 255 / colorgradiant 'values of 0 to colorgradiant determine the depth shading front to back
rbright = .5: gbright = .2: bbright = 0 ' CHANGES THE RGB DEPTH SHADING / BRIGHTNESS OF THE SKULL BROKEN DOWN INTO SEPARATE RGB VALUES

'========================= MAIN LOOP ============================
MainLoop:
mousecheck:
vermove = 0: hormove = 0: bigmove = 0: mb1 = 0: mb2 = 0: mb3 = 0
IF NOT _MOUSEINPUT THEN GOTO controls
vermove = _MOUSEMOVEMENTY: hormove = _MOUSEMOVEMENTX
mb1 = _MOUSEBUTTON(1): mb2 = _MOUSEBUTTON(2): mb3 = _MOUSEBUTTON(3)
WHILE _MOUSEINPUT: WEND
IF ABS(hormove) > ABS(vermove) THEN bigmove = hormove ELSE bigmove = vermove

controls:
mspeed = .1 + ABS(bigmove) * .1 'DEFAULT, mspeed smoothes out the model rotation in FREE and X,Y,Z modes
IF _KEYHIT = 27 THEN SYSTEM '[ESC] TO END THE PROGRAM
IF mb3 = -1 THEN posize = posize + (vermove * .01) 'PRESS/HOLD mouse wheel (button 3) to scale skull size up / down

freespin:
IF mb1 = -1 AND mb2 = -1 THEN zrot = zrot + bigmove * mspeed: GOTO checkit 'Z axis rotation
IF mb1 = -1 THEN xrot = xrot + (vermove * mspeed): yrot = yrot - (hormove * mspeed) 'X and Y axis rotation
checkit:
IF xrot >= 360 THEN xrot = xrot - 360 ELSE IF xrot < 0 THEN xrot = 360 + xrot
IF yrot >= 360 THEN yrot = yrot - 360 ELSE IF yrot < 0 THEN yrot = 360 + yrot
IF zrot >= 360 THEN zrot = zrot - 360 ELSE IF zrot < 0 THEN zrot = 360 + zrot

'Calculate 3d points (longitude/latitude) and translate into 2D X/Y screen coordinates
cz = COS(cf# * zrot): sz = SIN(cf# * zrot)
cx = COS(cf# * xrot): sx = SIN(cf# * xrot)
cy = COS(cf# * yrot): sy = SIN(cf# * yrot)
FOR tri = 1 TO totrec
    I = cf# * head#(tri, 1): j = cf# * head#(tri, 2): r = head#(tri, 3) * posize: a = r * SIN(I) * SIN(j): b = r * COS(j): c = r * COS(I) * SIN(j)
    a1 = a * cy - c * sy: c1 = a * sy + c * cy: b2 = b * cx - c1 * sx: c2 = b * sx + c1 * cx: A3 = a1 * cz - b2 * sz: b3 = a1 * sz + b2 * cz
    dr = c2 / (zobs - c2) + 1 'Distance ratio for computing perspective
    puthead(tri, 1) = INT(A3 * dr + xcenter): puthead(tri, 2) = INT(b3 * -dr + ycenter) 'contains the physical 2D screen x,y coordinates for 3D points
    depth(tri, 0) = c2: depth(tri, 1) = tri ' sort array index contains the point number(tri,1) and the 3d depth(tri,0) of that point
NEXT tri
start = 1: finish = totrec: CALL QuickSort(start, finish, depth()) 'SORT ALL POINTS FROM BACK TO FRONT

'render the sorted skull points
CLS
LOCATE 1, 1: PRINT "X,Y rotation = Left mouse"
LOCATE 2, 1: PRINT "Z rotation   = Left + Right mouse"
LOCATE 3, 1: PRINT "Scale / Size = Press/hold mouse wheel and move mouse up/down (currently"; STR$(posize); ")"
LOCATE 4, 1: PRINT "[esc]        = Exit"
FOR show = 1 TO totrec
    LINE (puthead(depth(show, 1), 1) - xpointsize, puthead(depth(show, 1), 2) - ypointsize)-(puthead(depth(show, 1), 1) + xpointsize, puthead(depth(show, 1), 2) + ypointsize), _RGB32(INT(show / colorgradiant) * rbright, INT(show / colorgradiant) * gbright, INT(show / colorgradiant) * bbright), BF
NEXT show
_DISPLAY
GOTO MainLoop

'*** SORT POINT DEPTHS ********* BACK TO FRONT ***** this section adapted from QB64 examples
SUB QuickSort (start AS INTEGER, finish AS INTEGER, depth() AS SINGLE)
    DIM Hi AS INTEGER, Lo AS INTEGER, Middle AS SINGLE
    Hi = finish: Lo = start
    Middle = depth((Lo + Hi) / 2, 0) 'find middle of array
    DO
        DO WHILE depth(Lo, 0) < Middle: Lo = Lo + 1: LOOP
        DO WHILE depth(Hi, 0) > Middle: Hi = Hi - 1: LOOP
        IF Lo <= Hi THEN
            SWAP depth(Lo, 0), depth(Hi, 0): SWAP depth(Lo, 1), depth(Hi, 1)
            Lo = Lo + 1: Hi = Hi - 1
        END IF
    LOOP UNTIL Lo > Hi
    IF Hi > start THEN CALL QuickSort(start, Hi, depth())
    IF Lo < finish THEN CALL QuickSort(Lo, finish, depth())
END SUB