/* CMD: Clay Morph
 *
 * Set up your control points in a layer by themselves and call ClayInit.
 * Then put your object to be morphed in the background layer.
 * Then move your control points and call ClayMorph.
 * That's all there is to it.
 *
 * Note that ClayMorph only moves points.  It doesn't make new polygons.
 * So make sure your object is subdivided enough.
 *
 */
        call addlib("LWModelerARexx.port",0)
        signal on error
        signal on syntax
        call addlib("rexxmathlib.library",0,-30,0)
        call addlib("rexxsupport.library",0,-30,0)
        call main
        exit

        syntax:
        error:
        t=Notify(1,'!Rexx Script Error','@'ErrorText(rc),'Line 'SIGL)
        if (mxx_add) then call remlib(mxx)
        exit

MAIN:

pi = 3.1415926385897932384626

call random(,,time('s'))

syscode = "Clay Morph"
claynam = "ENV:clay.tmp"

/* Store current and background layers and first empty layer. */
fg = curlayer()
bg = curblayer()
emp = emptylayers()
if (words(emp) = 0) then do
    ok = notify(2,'!No empty layer.','Morph will be placed','on current layer.','Original will be','!nuked.')
    if ok = 0 then return
    tlt = 0
end
if (words(emp) > 0) then do
    tl = word(emp,1)
    tlt = 1
end


/* Get original clay positions */

if (~open(state, claynam, 'R')) then break
clays = readln(state)
call meter_begin(clays*2, syscode, "Reading Clay")
ccx = 0
ccy = 0
ccz = 0
do i=1 to clays
    xyz = readln(state)
    parse var xyz ox.i oy.i oz.i
    ccx = ccx + ox.i
    ccy = ccy + oy.i
    ccz = ccz + oz.i
    call meter_step()
end i

do i=1 to clays
    x = ox.i-ccx/clays
    y = oy.i-ccy/clays
    z = oz.i-ccz/clays
    call GETROT
    orx.i = rox
    ory.i = roy
    orz.i = roz
    call meter_step()
end i

call meter_end()
call close(state)

n = xfrm_begin()
if (n ~= clays) then call notify(1,'Clay has been damaged.','Try again.')
if (n ~= clays) then return
call meter_begin(n, syscode, "Computing changes")
do i=1 to n
    xyz = xfrm_getpos(i)
    parse var xyz nx ny nz
    x = nx-ccx/clays
    y = ny-ccy/clays
    z = nz-ccz/clays
    call GETROT
    drx.i = rox-orx.i
    dry.i = roy-ory.i
    drz.i = roz-orz.i
    call meter_step()
end i
call meter_end()
call xfrm_end()


call setlayer(bg)
n = xfrm_begin()
if (n = 0) then return
call meter_begin(n, syscode, "Morphing Points")
do i=1 to n
    xyz = xfrm_getpos(i)
    parse var xyz x y z
    x = x-ccx/clays
    y = y-ccy/clays
    z = z-ccz/clays
    call GETROT
    x = x+ccx/clays
    y = y+ccy/clays
    z = z+ccz/clays
    cx = rox
    cy = roy
    cz = roz
    call Clay
    call SETROT
    call xfrm_setpos i,x y z
    call meter_step()
end i
call meter_end()
call xfrm_end()

CLAY:
dtot = 0
dmx = 0
dmy = 0
dmz = 0
do c=1 to clays
 if dtot ~= -1 then do
  dist = sqrt((x-ox.c)**2 + (y-oy.c)**2 + (z-oz.c)**2)
  if dist = 0 then do
   dtot = -1
   dmx = drx.c
   dmy = dry.c
   dmz = drz.c
  end
  if dist ~= 0 then do
   dtot = dtot + 1/dist
   dmx = dmx + drx.c/dist
   dmy = dmy + dry.c/dist
   dmz = dmz + drz.c/dist
  end
 end
end
/* normalize */
if dtot ~= -1 then do
 dmx = dmx/dtot
 dmy = dmy/dtot
 dmz = dmz/dtot
end
cx = cx+dmx
cy = cy+dmy
cz = cz+dmz
return

GETROT:
/* 3-D Radial Rotation */
/* takes x,y,z, returns with rox,roy,roz */
rox = 0
roy = 0
roz = sqrt(x*x+y*y+z*z)
/* First we adjust our heading */
yf=0
if z=0 then do
    if x>0 then roy=roy+pi/2
    if x<0 then roy=roy-pi/2
    yf = 1
end
if x=0 then do
    if z>0 then roy=roy+0
    if z<0 then roy=roy+pi
    yf = 1
end
if yf=0 then do
    actn=abs(x/z)
    call arctan
    if z>0 then do           /* quadrant 1 */
        if x>0 then roy=roy+th
    end
    if z>0 then do           /* quadrant 2 */
        if x<0 then roy=roy-th
    end
    if x<0 then do           /* quadrant 3 */
        if z<0 then roy=roy+pi+th
    end
    if z<0 then do           /* quadrant 4 */
        if x>0 then roy=roy+pi-th
    end
end
/* Now we pitch up */
nd=sqrt(z*z+x*x)
yf=0
if nd=0 then do
    if y>0 then rox=rox+pi/2
    if y<0 then rox=rox+-pi/2
    yf = 1
end
if y=0 then yf = 1
if yf=0 then do
    actn=abs(y/nd)
    call arctan
    if y>0 then rox=rox+th
    if y<0 then rox=rox-th
end
return

ARCTAN:
actu = actn
if actu < -1 then th = pi/2+atan(1/abs(actn))
if (actu < 1.00001)&(actu > -1.00001) then th = atan(actn)
if actu > 1 then th = pi/2-atan(1/actn)
return

SETROT:
/* takes cx,cy,cz & returns with x,y,z */
/* pitch */
x = 0
y = sin(cx)
z = cos(cx)
/* heading */
x = sin(cy)*z
y = y
z = cos(cy)*z
/* radius */
x = x * cz + ccx/clays
y = y * cz + ccy/clays
z = z * cz + ccz/clays
return