;+
;
; Project   : Densitycube, reconstruction, 3D, and s3drs software
;                   
; Name      : cuberot3d
;               
; Purpose   : used to manipulate 3D densitycubes and do 3D reconstructions
;               
; Explanation: rotates a 3D array around an arbitrary center
;               
; Use       : standalone or via s3drs GUI
;    
; Inputs    : MRI2 (cube), phi (angles)
;               
; Outputs   : MRI2 (input cube is modified)
;
; Keywords  : trans, render, xyzcenter, reverse, doublepec
;               
; Common    : none
;               
; Restrictions: none
;               
; Side effects: none
;               
; Category    : 3d, density, rotation
;               
; Prev. Hist. : None.
;
; Written     : Sandy Antunes, NRC, March-April 2009
;               
;-            
; a fast but not ultra-precise rotation program
; from Usenet, by Marc Shellens, Sept 2001,
; transcribed and embedded by Sandy Antunes, Dec 2005
; Modified to allow different center than 'cube'
;
; If given /reverse, does rotation in reverse order
;
; Original source and discussion:
; http://groups.google.com/group/comp.lang.idl-pvwave/browse_thread/thread/6626c477c2eb4991/d26bf81e8cc5f77f?lnk=st&q=idl+fast+cube+array+rotation&rnum=5&hl=en#d26bf81e8cc5f77f
;
; Note it presumes cubes are floats (it casts each slice into an NxN float
; array), if you are using double cubes, this may cause a small loss
; in fidelity unless you include the /doubleprec flag.
;

PRO cuberot3d,MRI2,phi,trans=trans,render=render,xyzcenter,reverse=reverse,$
              doubleprec=doubleprec


;phi fltarr(3) the three angles to rotate (+/- is convention)
;trans intarr(3) the translation (in voxels)

  reverse=KEYWORD_SET(reverse)
  doubleprec=KEYWORD_SET(doubleprec)

  large=0

  if (datatype(MRI2) eq 'PTR') then begin
      ; experimental handling of a large array handling variant
      N=n_elements(MRI2)
      large=1
;      print,'using large handler'
  end else begin
      N=(size(MRI2))[1]
  end


  if (n_elements(xyzcenter) eq 0) then $
    xyzcenter=[(N-1)/2.0,(N-1)/2.0,(N-1)/2.0]

  if (n_elements(trans) eq 0) then trans=[0,0,0]

;print,'X...'
  if phi[0] ne 0.0 and reverse eq 0 then begin
      for x=0,N-1 do begin

          if (doubleprec) then mycut=dblarr(N,N) else mycut=fltarr(N,N)
          if (large) then begin
              for m=0,N-1 do mycut[*,m] = (*MRI2[m])[x,*]
              mycut = rot(/INTERP, reform(mycut,N,N),-phi[0],$
                          1.0,xyzcenter[1],xyzcenter[2],MISSING=0,/pivot)
              for m=0,N-1 do (*MRI2[m])[x,*] = mycut[*,m]
          end else begin
              MRI2[x,*,*]=rot(/INTERP,reform(MRI2[x,*,*],N,N),-phi[0],$
                              1.0,xyzcenter[1],xyzcenter[2],MISSING=0,/pivot)

          end

      endfor
  endif

;print,'Y...'
  if phi[1] ne 0.0 and reverse eq 0 then begin
      for y=0,N-1 do begin
          
          if (doubleprec) then mycut=dblarr(N,N) else mycut=fltarr(N,N)
          if (large) then begin
              for m=0,N-1 do mycut[*,m] = (*MRI2[m])[*,y]
              mycut = rot(/INTERP, reform(mycut,N,N),-phi[1],$
                          1.0,xyzcenter[0],xyzcenter[2],MISSING=0,/pivot)
              for m=0,N-1 do (*MRI2[m])[*,y] = mycut[*,m]
          end else begin
              MRI2[*,y,*]=rot(/INTERP,reform(MRI2[*,y,*],N,N),phi[1],$
                          1.0,xyzcenter[0],xyzcenter[2],MISSING=0,/pivot)
          end

      endfor
  endif

;print,'Z...'
  if phi[2] ne 0.0 then begin
      for z=0,N-1 do begin
          if (large) then begin
              (*MRI2[z])[*,*]=rot(/INTERP,(*MRI2[z])[*,*],-phi[2],$
                          1.0,xyzcenter[0],xyzcenter[1],MISSING=0,/pivot)
          end else begin
              MRI2[*,*,z]=rot(/INTERP,MRI2[*,*,z],-phi[2],$
                          1.0,xyzcenter[0],xyzcenter[1],MISSING=0,/pivot)
          end
      endfor
  endif


;print,'Y...'
  if phi[1] ne 0.0 and reverse eq 1 then begin
      for y=0,N-1 do begin
          
          if (doubleprec) then mycut=dblarr(N,N) else mycut=fltarr(N,N)
          if (large) then begin
              for m=0,N-1 do mycut[*,m] = (*MRI2[m])[*,y]
              mycut = rot(/INTERP, reform(mycut,N,N),-phi[1],$
                          1.0,xyzcenter[0],xyzcenter[2],MISSING=0,/pivot)
              for m=0,N-1 do (*MRI2[m])[*,y] = mycut[*,m]
          end else begin
              MRI2[*,y,*]=rot(/INTERP,reform(MRI2[*,y,*],N,N),phi[1],$
                          1.0,xyzcenter[0],xyzcenter[2],MISSING=0,/pivot)
          end

      endfor
  endif

;print,'X...'
  if phi[0] ne 0.0 and reverse eq 1 then begin
      for x=0,N-1 do begin

          if (doubleprec) then mycut=dblarr(N,N) else mycut=fltarr(N,N)
          if (large) then begin
              for m=0,N-1 do mycut[*,m] = (*MRI2[m])[x,*]
              mycut = rot(/INTERP, reform(mycut,N,N),-phi[0],$
                          1.0,xyzcenter[1],xyzcenter[2],MISSING=0,/pivot)
              for m=0,N-1 do (*MRI2[m])[x,*] = mycut[*,m]
          end else begin
              MRI2[x,*,*]=rot(/INTERP,reform(MRI2[x,*,*],N,N),-phi[0],$
                              1.0,xyzcenter[1],xyzcenter[2],MISSING=0,/pivot)

          end

      endfor
  endif

;print,'shift...'

  if (trans[0] ne 0 or trans[1] ne 0 or trans[2] ne 0) then begin
      print,'transing'
      if (large) then begin
          for m=0,N-1 do *MRI2[m]=shift(*MRI2[m],trans[0],trans[1])
          MRI2=shift(temporary(MRI2),trans[2])
      end else begin
          MRI2=shift(temporary(MRI2),trans[0],trans[1],trans[2]) 
      end
  end

END