Re: arbitrary rotation of 3-d arrays

• Subject: Re: arbitrary rotation of 3-d arrays
• From: steinhh(at)ulrik.uio.no (Stein Vidar Hagfors Haugan)
• Date: 11 Jun 1999 09:31:59 GMT
• In-reply-to: "D. Mattes"'s message of Thu, 10 Jun 1999 10:33:45 -0700
• Newsgroups: comp.lang.idl-pvwave
• Organization: University of Oslo, Norway
• Xref: news.doit.wisc.edu comp.lang.idl-pvwave:15146

```

> has anybody out there in idl-land written or seen code to apply arbitrary
> rotations to 3-d arrays???

Although others have posted solutions using the t3d style
rotations, you might want to look at the procedures below.
I don't really like the idea of using a (global) system
variable designed for 3d *graphics* as "the" temporary
variable to accumulate all kinds of 3d manipulations...

I'm sorry for the complete lack of documentation ... this
was something I did to experiment myself towards an understanding
of such rotations..  ROT3DMATRIX returns an array that to
be applied like this

ROTATED_XYZ = ROT3DMATRIX([alphax,alphay,alphaz]) ## XYZ_ARRAY

(Note that [alphax,alphay,alphaz] is wrt. a fixed coordinate
system, the axes don't move after each partial rotation...
this may be different from the way t3d applies its input..)

Also, I wrote a test application (ROTATE_SCREW) that uses
direct graphics to manipulate 3D objects on screen with the
help of a trackball object.

Try clicking either left *and* middle buttons to rotate the
box, and the (two!) corkscrews inside the box. Double-clicking
the middle button changes the "sense" of how the "loose"
corkscrew is rotated - wrt the *box* coordinate system or wrt
the *screen* (sort of) coordinate system.... You may get the
difference if you e.g. turn the box 180 degrees around and
try manipulating the loose corkscrew...

Ok, here goes,

Stein Vidar

---------------------------------------------

;;
;; Return rotation matrix such that
;; rot3dmatrix([alphax,alphay,alphaz]) ## [X,Y,Z] gives your [X,Y,Z]
;; the y axis, and finally alphaz radians about the z axis. Note that
;; the axes are kept fixed. This allows an inverse operation to to
;; return alphax, alphay, alphaz from a given rotation matrix.

FUNCTION rot3dmatrix_angles,m

m = double(m)

;  From mathematica:
;
;  mm = [[Cy Cz, Cz Sx Sy - Cx Sz, Cx Cz Sy + Sx Sz   ], \$
;        [Cy Sz, Cx Cz + Sx Sy Sz, -(Cz Sx) + Cx Sy Sz],\$
;        [-Sy  , Sx Cy           ,     Cx Cy          ]]

cosyzero = (m(1,2) EQ 0 AND m(2,2) EQ 0)

IF NOT cosyzero THEN BEGIN

cys = 1

xfind = atan(cys*m(1,2),cys*m(2,2)) ;; May be wrong quadrants!
zfind = atan(cys*m(0,1),cys*m(0,0))

sinx = sin(xfind) & sxgood = (abs(sinx) GT 0.05)
cosx = cos(xfind) & cxgood = (abs(cosx) GT 0.05)

wt = double(sxgood+cxgood)
cosy = ((sxgood ? m(1,2)/sinx : 0)+(cxgood ? m(2,2)/cosx : 0))/wt

yfind = atan(-m(0,2),cosy)

return,[xfind,yfind,zfind]
END

;; Cos[y] == 0 => yfind = +/- Pi/2

yfind = atan(-m(0,2),0)

;; From mathematica we find that

;; MatrixForm[TrigFactor[r[x,+/-Pi/2,z]]]
;;
;; = [[ 0,  +/-Sin[x-z], +/-Cos[x+z]  ],\$
;;    [ 0,     Cos[x-z],    -Sin[x-z] ],\$
;;    [ -/+1,         0,           0  ]]

;; We thus arbitrarily set z = 0 and get:

zfind = 0
xfind = atan(-m(2,1),m(1,1))

return,[xfind,yfind,zfind]
END

FUNCTION rot3dmatrix,alpha,inverse=inverse

IF keyword_set(inverse) THEN return,rot3dmatrix_angles(alpha)

alpha = double(alpha)
ca = cos(alpha)
sa = sin(alpha)

mx = [[     1,     0,      0],\$
[     0, ca(0), -sa(0)],\$
[     0, sa(0),  ca(0)]]

my = [[ ca(1),     0,  sa(1)],\$
[     0,     1,      0],\$
[-sa(1),     0,  ca(1)]]

mz = [[ ca(2),-sa(2),      0],\$
[ sa(2), ca(2),      0],\$
[     0,     0,      1]]

return,mz ## (my ## mx)
END

;------------------------------------------------------------

PRO plotcube,xr,yr,zr

xi = [0,1,1,0,0]
yi = [0,0,1,1,0]
zi = [0,0,0,0,0]

plots,xr(xi),yr(yi),zr(zi),/t3d,/data
plots,xr(xi),yr(yi),zr(zi+1),/t3d,/data
FOR i=0,4 DO plots,xr(xi([i,i])),yr(yi([i,i])),zr([0,1]),/t3d,/data

END

PRO rotate_screw,rotation

;; Create widget draw window

xs = (ys=512) ;; Size of draw window

id = widget_base()
dummy = widget_draw(id,xsize=xs,ysize=ys,/button_ev,/motion)
widget_control,id,/realize
widget_control,dummy,get_value=win
wset,win

xrange = [(xmin=-10), (xmax=10)]
yrange = [(ymin=-10), (ymax=10)]
zrange = [(zmin=-10), (zmax=10)]

!x.s = [-xmin,1.0]/(xmax-xmin)
!y.s = [-ymin,1.0]/(ymax-ymin)
!z.s = [-zmin,1.0]/(zmax-zmin)

theta = findgen(120)/199.0*2*!PI
x = cos(20*theta)
y = sin(20*theta)
z = 5*theta

xyz = [[x],[y],[z]]
t3d,/reset,translate=-[.5,.5,.5] & xyzform = !P.t

;; X/Y/Z "axis" vectors for easy drawing

xa = 8*[[0,1],[0,0],[0,0]]
ya = 8*[[0,0],[0,1],[0,0]]
za = 8*[[0,0],[0,0],[0,1]]

h = [.5,.5,.5]
persp = 5
scale = .6

;; Build the initial viewing matrix

t3d,/reset,trans=-h
t3d,scale=.6*[1,1,1]
IF n_elements(rotation) EQ 3 THEN BEGIN
print,"Rotation"
END

t3d,trans=h

track = obj_new('trackball',[xs/2.,ys/2.0],0.25*xs)
track2 = obj_new('trackball',[xs/2.,ys/2.0],0.25*xs)

inspace = 1

REPEAT BEGIN
t = !P.t
t3d,trans=-h
angles = rot3dmatrix(!p.t(0:2,0:2),/inverse)
t3d,perspect=persp
t3d,trans=h

erase

plotcube,xrange,yrange,zrange

plots,transpose(xa),/t3d,/data
plots,transpose(ya),/t3d,/data,color=140
plots,transpose(za),/t3d,/data,color=100

xyouts,xa(1,0),xa(1,1),z=xa(1,2),"X",/t3d,/data
xyouts,ya(1,0),ya(1,1),z=ya(1,2),"Y",/t3d,/data
xyouts,za(1,0),za(1,1),z=za(1,2),"Z",/t3d,/data

plots,x,y,z,/t3d,/data
plots,transpose(xyz),/t3d,/data

xyouts,.1,.15,string(angles(0))+"!c"+string(angles(1))+\$
"!c"+string(angles(2))+"!c"+"insp:"+string(inspace),/normal

!P.t = t
empty
ev = widget_event(id)
IF ev.press EQ 2 AND ev.clicks EQ 2 THEN inspace = (inspace + 1) MOD 3

xformq = track->update(ev,transform=xform,mouse=1b)
IF xformq THEN BEGIN
t3d,translate=-h
!P.t = xform ## !p.t
t3d,translate=h
END

xformq = track2->update(ev,transform=xform,mouse=2b)
IF xformq THEN BEGIN
IF inspace EQ 1 THEN BEGIN

;; Now - rotate/shift the screw into the orientation it has on the
;; screen, then apply trackball rotation, then rotate/shift back
;; into it's own space.

xform = invert(!p.t) ## xform ## !p.t
END ELSE IF inspace EQ 2 THEN BEGIN
;; Rotate the screw "in its own data space".
;; Apply inverse transform on xyz,
;;
xform = xyzform ## xform ## invert(xyzform)
END

;; Apply the resulting transform on the screw. The xform really
;; includes a shift (since we placed it centered on the screen, not on
;; the coordinate axes), but this is not taken into account since
;; we're using only (0:2,0:2) of the resulting transform

xyz = xform(0:2,0:2) ## xyz

;; We need to keep track of the transform that has been applied
;; in order to do transformations in
xyzform = xform ## xyzform
END
END UNTIL ev.press EQ 4

widget_control,id

rotation = angles
END

```