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