Volume Transformation

QUESTION: I have two volume data sets and I am trying to rotate and translate them to match each other. It would be very nice if I could use something like "XVOLUME_ROTATE, /T3D" or "XVOLUME_ROTATE, MATRIX=!P.T." (Of course this isn't possible). What is the best way to accomplish this in the absence of such a program?

ANSWER: There have been numerous discussions on volume rotations and transformations on the IDL newsgroup over the years. Here is a volume transformation program that has been recommended recently. It was written by Martin Downing in the UK. The transformation can be performed in "chunks", which improves memory efficiency.

The short program is printed here (without it's accompaning documentation). You can download the entire program by clicking the link.

```FUNCTION Transform_Volume, volume, Rotation=rotation, \$
Scale=scale, Translate=translate, Centre_Rotation=centre_rotation, \$
T3Dmat=t3dmat, Inplace=inplace, Buffer_Size=buffer_size

; Error handling.

Catch, theError
IF theError NE 0 THEN BEGIN
Catch, /Cancel
ok = Dialog_Message(!Error.State.Msg)
IF N_Elements(volume) EQ 0 THEN volume = Temporary(vol_t)
RETURN, -1
ENDIF

; Find the dimensions of the volume.

s = Size(volume)
sx=s[1] & sy=s[2] & sz=s[3]
st = sx*sy*sz

IF Keyword_Set(inplace) THEN vol_t = Temporary(volume) ELSE vol_t = volume

; Create a transform matrix, if one is not provided.

IF N_Elements(t3dmat) EQ 0 THEN begin

IF N_Elements(rotation) EQ 0 THEN rotation =[0,0,0]
IF N_Elements(centre_rotation) EQ 0  THEN centre_rotation=[(sx-1)/2.0,(sy-1)/2.0,(sz-1)/2.0]
IF N_Elements(translate) EQ 0 THEN translate =[0,0,0]
IF N_Elements(scale) EQ 0 THEN scale =[1,1,1]

T3D, /Reset, Translate = -centre_rotation
T3D, Rotate=rotation
T3D, Translate= centre_rotation + translate, Scale=scale
t3dmat = !P.T

ENDIF

; Check buffer size. The size 128 is optimim on my system, You may
; want to try other values.

IF N_Elements(buffer_size) EQ 0 THEN buffer_size = 128
IF buffer_size LE 0 THEN buffer_size = st

; Perform the transformations.

FOR j=0L,(st-1),buffer_size DO BEGIN

; Account for possible odd last chunk.

bufsize = buffer_size < (st-j)

; Generate volume coordinates by interpolating temporary array of volume indices.

i = j + Lindgen(bufsize)
coords = [ [(i MOD sx)],[((i / sx) MOD (sy))], [(i / (sx * sy))], [Replicate(1b, bufsize)]]
coords = Temporary(coords) # t3dmat
vol_t[j:j+bufsize-1] = Interpolate(volume, coords[*,0], coords[*,1], coords[*,2], Missing=0)
ENDFOR

; Return the transformed volume.

RETURN, vol_t
END
```

 Web Coyote's Guide to IDL Programming