;+
;
; Project   : s3drs reconstruction software
;                   
; Name      : (various test routines)
;               
; Purpose   : part of the s3drs test kit
;               
; Explanation: (archived here for completeness, not intended for users use)
;               
; Use       : (self-documented, see notes in code)
;    
; Written     : Sandy Antunes, NRC, 2005-2009
;               
;-            
;IDL users:  thanks for the lively discussion regarding this topic.  i
;never did find a canned procedure for me to use, but i implemented a 3D
;volume rotation class method using VERT_T3D and INTERPOLATE.  i submit it
;here for anybody who would like it, and for suggestions for optimization
;and improvement.  
;
;cheers,
;david mattes
;http://groups.google.com/group/comp.lang.idl-pvwave/browse_thread/thread/bce4883c35398a37/339e33701f8e88be?hl=en&lnk=st&q=#339e33701f8e88be

;pro ApplyRigidTransformation,image,U,t
pro interpolate3D,image,U,t

     ;U is a 3x3 unitary rotation matrix
     ;t is a 3x1 vector of x,y,z translations

     if (n_elements(U) eq 0) then U=make_array(3,3,/float,value=0.0);def=no rot
     if (n_elements(t) eq 0) then t=[0.0,0.0,0.0]; default=no translation


     Urot=U
     ; shorthand here, if we only get 3 elements, assume degrees and rotate
     if (n_elements(Urot) eq 3) then begin
         Urot=rot3dmatrix(U*!dtor)
     end

     xdim=n_elements(image[*,0,0])
     ydim=n_elements(image[0,*,0])
     zdim=n_elements(image[0,0,*])

     ;build 4x4 homogeneous transformation matrix
    localT=fltarr(4,4)
    localT(0:2,0:2)=Urot
    localT(3,0:2)=t
    localT(3,3)=1.

     ;build x,y,z interpolation points
    vert_size=LONG(xdim)*$
        LONG(ydim)*$
        LONG(zdim)
    verts=fltarr(3,vert_size)
    count=0L
    for i=0,zdim-1 do begin
        for j=0,ydim-1 do begin
            for k=0,xdim-1 do begin
                 ;notice index swapping here!!!!!!!!!!!!
                 ;for our image, the x coordinate is the most
                 ;quickly varying, so it must also be this way for the
                 ;interpolate locations
                verts(*,count)=[k,j,i]
                count=count+1
            endfor
        endfor
    endfor

     ;verts are the location points for which we want
     ;interpolates...really just the indices in the volume array.
     ;we transform these indices using the localT transformation
     ;matrix, and pass them to the interpolate function.
     ;notice: verts is 3 x (n*m*l) 2D matrix for a n x m x l array!!!
    verts=VERT_T3D(verts,MATRIX=localT,/NO_COPY,/NO_DIVIDE)

     ;right now, use missing=-1 for values outside input data range.
     ;cubic interpolation is not supported for 3-d case
    image=INTERPOLATE(TEMPORARY(image), $
                      verts(0,*),verts(1,*),verts(2,*),MISSING=-1)

     ;the interpolate function returns a 1-d array of interpolated
     ;points, which must be resized into the original array shape.
    image=REFORM(image, xdim, ydim, zdim,/OVERWRITE)

end