;+
;
; 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
;               
;-            
PRO triofaceoff
  ; a series of 45 degree rotations, we then compare within
  ; each rotation algorithm how much rotation perturbs the
  ; actual distribution.  Assumes symmetry.

  rot_triotest,timage,/slide,[45.0,45.0,45.0,45.0]
  restore,/v,'rot_triotest.sav'

  dptr=ptrarr(15)
  i=0

  dptr(i++)=ptr_new(dataset1[*,*,7]-dataset1[*,*,14])
  dptr(i++)=ptr_new(dataset1[*,*,4]-dataset1[*,*,10])
  dptr(i++)=ptr_new(dataset1[*,*,4]-dataset1[*,*,11])
  dptr(i++)=ptr_new(dataset1[*,*,1]-dataset1[*,*,8])
  dptr(i++)=ptr_new(dataset1[*,*,1]-dataset1[*,*,13])

  dptr(i++)=ptr_new(dataset2[*,*,7]-dataset2[*,*,14])
  dptr(i++)=ptr_new(dataset2[*,*,4]-dataset2[*,*,10])
  dptr(i++)=ptr_new(dataset2[*,*,4]-dataset2[*,*,11])
  dptr(i++)=ptr_new(dataset2[*,*,1]-dataset2[*,*,8])
  dptr(i++)=ptr_new(dataset2[*,*,1]-dataset2[*,*,13])

  dptr(i++)=ptr_new(datadiffs[*,*,7])
  dptr(i++)=ptr_new(datadiffs[*,*,4])
  dptr(i++)=ptr_new(datadiffs[*,*,4])
  dptr(i++)=ptr_new(datadiffs[*,*,1])
  dptr(i++)=ptr_new(datadiffs[*,*,1])

  tv_multi,fixed=128,/log,dptr,line=5

END



PRO minirot8interp

  angle=90.0

  xy=fltarr(8,8)
  xy(1,1)=3
  xy(2,2)=2
  xy(3,3)=1
;  tv_multi,fixed=64,/raw,xy,title='xy orig'

  xt1=xy*0.0
  xt1(1,6)=3
  xt1(2,5)=2
  xt1(3,4)=1
  xt2=xy*0.0
  xt2(1,2)=3
  xt2(2,1)=2
  xt2(3,0)=1

  xyr1=rot(xy,angle,/interp)
  xyr2=rot(xy,angle,/pivot,/interp)
  xyr3=rot(xy,angle,1.0,4,4,/interp)
  xyr4=rot(xy,angle,1.0,4,4,/pivot,/interp)
  xyr5=rot(xy,angle,1.0,3.5,3.5,/interp)
  xyr6=rot(xy,angle,1.0,3.5,3.5,/pivot,/interp)
  xyr7=rot(xy,angle,1.0,2,2,/interp)
  xyr8=rot(xy,angle,1.0,2,2,/pivot,/interp)
  xyr9=rot(xy,angle,1.0,1.5,1.5,/interp)
  xyra=rot(xy,angle,1.0,1.5,1.5,/pivot,/interp)

  xyset=[[[xy]],[[xt1]],[[xt2]],[[xyr1]],[[xyr2]],$
         [[xyr3]],[[xyr4]],[[xyr5]],[[xyr6]],[[xyr7]],[[xyr8]],$
         [[xyr9]],[[xyra]]]
  lab=['orig','90 theory','90off theory',$
       '90','90 p','90 ctr=4','90 ctr=4 p',$
       '90 ctr=3.5','90 ctr=3.5 p',$
       '90 off=2','90 off=2 p',$
      '90 off=1.5','90 off=1.5 p']

  tv_multi,fixed=64,line=2,/raw,xyset,/rev,lab=lab,/border
  tv_multi,/ps,fixed=64,title='rot 8x8, interpolating',$
    line=2,/raw,xyset,/rev,lab=lab,/border,file='rotinterp.ps'

END

PRO minirot8

  angle=90.0

  xy=fltarr(8,8)
  xy(1,1)=3
  xy(2,2)=2
  xy(3,3)=1
;  tv_multi,fixed=64,/raw,xy,title='xy orig'

  xt1=xy*0.0
  xt1(1,6)=3
  xt1(2,5)=2
  xt1(3,4)=1
  xt2=xy*0.0
  xt2(1,2)=3
  xt2(2,1)=2
  xt2(3,0)=1

  xyr1=rot(xy,angle)
  xyr2=rot(xy,angle,/pivot)
  xyr3=rot(xy,angle,1.0,4,4)
  xyr4=rot(xy,angle,1.0,4,4,/pivot)
  xyr5=rot(xy,angle,1.0,3.5,3.5)
  xyr6=rot(xy,angle,1.0,3.5,3.5,/pivot)
  xyr7=rot(xy,angle,1.0,2,2)
  xyr8=rot(xy,angle,1.0,2,2,/pivot)
  xyr9=rot(xy,angle,1.0,1.5,1.5)
  xyra=rot(xy,angle,1.0,1.5,1.5,/pivot)

  xyset=[[[xy]],[[xt1]],[[xt2]],[[xyr1]],[[xyr2]],$
         [[xyr3]],[[xyr4]],[[xyr5]],[[xyr6]],[[xyr7]],[[xyr8]],$
         [[xyr9]],[[xyra]]]
  lab=['orig','90 theory','90off theory',$
       '90','90 p','90 ctr=4','90 ctr=4 p',$
       '90 ctr=3.5','90 ctr=3.5 p',$
       '90 off=2','90 off=2 p',$
      '90 off=1.5','90 off=1.5 p']

  tv_multi,fixed=64,line=2,/raw,xyset,/rev,lab=lab,/border
  tv_multi,/ps,fixed=64,title='Rot 8x8, no interpolation',$
    line=2,/raw,xyset,/rev,lab=lab,/border,file='rotnointerp.ps'

END

PRO minirot

  xy=fltarr(4,4)
  xy(0,0)=2
  xy(1,1)=1
;  tv_multi,fixed=64,/raw,xy,title='xy orig'

  xt1=xy*0.0
  xt1(0,3)=2
  xt1(1,2)=1
  xt2=xy*0.0
  xt2(0,1)=2
  xt2(1,0)=1

  xyr1=rot(xy,90.0)
  xyr2=rot(xy,90.0,/pivot)
  xyr3=rot(xy,90.0,1.0,2,2)
  xyr4=rot(xy,90.0,1.0,2,2,/pivot)
  xyr5=rot(xy,90.0,1.0,1,1)
  xyr6=rot(xy,90.0,1.0,1,1,/pivot)

  xyset=[[[xy]],[[xt1]],[[xt2]],$
         [[xyr1]],[[xyr2]],[[xyr3]],[[xyr4]],[[xyr5]],[[xyr6]]]
  lab=['orig','90 theory','90off theory',$
       '90base','90base p','90ctr','90ctr p','90off','90off p']

  tv_multi,fixed=64,line=2,/raw,xyset,/rev,lab=lab,/border

END


pro loadme, timage
on_error,2
restore,/v,'/home/antunes/RUNS/data/dialme/fluxyB.img'
timage=fluxyB
render_rot_gui,timage

threeview,timage,/log,title='original, not rotated',/pixon

end

pro fourcube,timage

  timage=fltarr(4,4,4)
  timage[0,3,*]=10.0^5
  timage[1,2,*]=10.0^5

end

PRO offcenter_rottest,timage

  ; as before, but rotating where the sun center aka coord center
  ; is elsewhere

;  rot_triotest,timage,/slide,$
;    [0.0,45.0,90.0,45.0,0.0],[0.0,0.0,0.0,45.0,45.0];

  rot_triotest,timage,/slide,$
    [90.0,90.0,90.0,90.0],xyzcenter=[3,3,2]


END

; given a set of rotations, consecutively rotates.  A good
; test is give it a series of rotations plus their inverses
; to see how well it reconstrues the original.
;
; The toggle /pr_backproject runs the calculation as per
; 'pr_backproject' on just the 1st set of angles, to verify
; that we have the right order on our rotations.
; /pr_render does the same as per our pr_render routine.
; So /pr_render will rotate, then rotate back (a check to make
; sure we restore our original.
; /pr_backproject will do just the backward rotation (cube-centric)
;

PRO rot_triotest,timage,rhoset,phiset,N=N,res=res,$
                 pr_render=pr_render,pr_backproject=pr_backproject,$
                 slideshow=slideshow,ps=ps,pixon=pixon,xyzcenter=xyzcenter

  if (n_elements(timage) eq 0) then $
    restore,/v,'/home/antunes/RUNS/data/dialme/fluxyB.img'

  pr_render=setyesno(pr_render)
  pr_backproject=setyesno(pr_backproject)
  slideshow=setyesno(slideshow)
  ps=setyesno(ps)
  pixon=setyesno(pixon)

  if (n_elements(rhoset) eq 0) then rhoset=[45.0]
  if (n_elements(phiset) eq 0) then phiset=rhoset*0.0

  if (n_elements(timage) eq 0) then loadme,timage

  if (n_elements(N) ne 0) then timage=congrid(timage,N,N,N)

  image1=timage
  image2=timage
  image3=timage

  if (pr_render or pr_backproject) then begin
      ; these tests work only one angle set
      rhoset=[rhoset[0]]
      phiset=[phiset[0]]
  end

  threeview,fixed=64,/raw,image1,dataset1,/noplot,pixon=pixon

  starttime=systime(/seconds)

  phitally=0
  rhotally=0

  for i=0,n_elements(rhoset)-1 do begin

      xang=0.0                  ; always 0
      yang=phiset[i]                  ; height angle above plane
      zang=rhoset[i]              ; CCW around z

      phitally=phitally+yang
      rhotally=rhotally+zang

      if (pr_render) then begin
          cuberot3d,image1,[0.0,0.0,0.0-zang],xyzcenter
          cuberot3d,image1,[xang,yang,0.0],xyzcenter
          threeview,fixed=64,/raw,image1,view1,/noplot,pixon=pixon
          dataset1=[[[view1]],[[dataset1]]]
          cuberot3d,image1,[xang,0.0-yang,0.0],xyzcenter
          cuberot3d,image1,[0.0,0.0,zang],xyzcenter
      end else if (pr_backproject) then begin
          cuberot3d,image1,[xang,0.0-yang,0.0],xyzcenter
          cuberot3d,image1,[0.0,0.0,zang],xyzcenter
      end else begin
          print,phitally,rhotally
          cuberot3d,image1,[0.0,0.0,0.0-zang],xyzcenter
          cuberot3d,image1,[xang,yang,0.0],xyzcenter
      end

      threeview,fixed=64,/raw,image1,view1,/noplot,pixon=pixon
      if (slideshow or pr_render) then dataset1=[[[view1]],[[dataset1]]] else $
        dataset1=view1

  end
  print,'cuberot3d ran in ',systime(/seconds)-starttime,' seconds'

; turn_3d runs slower than cuberot3d (6x slower) but yields
; better results.
; also, turn_3d with the /reverse mod runs only 4x slower

  threeview,fixed=64,/raw,image2,dataset2,/noplot,pixon=pixon

  starttime=systime(/seconds)

  for i=0,n_elements(rhoset)-1 do begin

      xang=0.0                  ; always 0
      yang=phiset[i]                  ; height angle above plane
      zang=rhoset[i]              ; CCW around z

; turn_3d expands by 50% to preserve info
; double pass, to preserve order
;;      image2=turn_3d(image2,0.0,0.0,0.0-zang)
;;      image2=image2[33:160,33:160,33:160]
;;      image2=turn_3d(image2,xang,yang,0.0)
;;      image2=image2[33:160,33:160,33:160]
; single pass, for speed
      print,xang,yang,zang,xyzcenter

      if (pr_render) then begin
          image2=turn_3d(image2,xang,yang,zang,/reverse,/resize,$
                         xyzcenter=xyzcenter,/pivot,/interp)
          ; include intermediate frame
          threeview,fixed=64,/raw,image2,view2,/noplot,pixon=pixon
          dataset2=[[[view2]],[[dataset2]]]
          image2=turn_3d(image2,xang,0.0-yang,0.0-zang,/resize,$
                         xyzcenter=xyzcenter,/pivot,/interp)
      end else if (pr_backproject) then begin
          image2=turn_3d(image2,xang,0.0-yang,0.0-zang,/resize,$
                         xyzcenter=xyzcenter,/pivot,/interp)
      end else begin
          image2=turn_3d(image2,xang,yang,zang,/reverse,/resize,$
                         xyzcenter=xyzcenter,/pivot,/interp)
      end

;      image2=image2[33:160,33:160,33:160]

      threeview,fixed=64,/raw,image2,view2,/noplot,pixon=pixon
      if (slideshow or pr_render) then dataset2=[[[view2]],[[dataset2]]] else $
        dataset2=view2

  end
  print,'turn_3d ran in ',systime(/seconds)-starttime,' seconds'


  datadiffs=dataset1-dataset2
  
  tv_multi,/border,/raw,fixed=64,dataset1,line=3,/log,$
    title='cuberot3d',res=res
  tv_multi,/border,/raw,fixed=64,dataset2,line=3,/log,$
    title='turn_3d',res=res
  tv_multi,/border,/raw,fixed=64,datadiffs,$
    line=3,/log,title='difference',res=res

  if (ps) then begin
      tv_multi,/border,/raw,fixed=64,dataset1,line=3,/log,$
        title='cuberot3d',res=res,/ps,file='s.ps'
      tv_multi,/border,/raw,fixed=64,dataset2,line=3,/log,$
        title='turn_3d',res=res,/ps,file='t.ps'
      tv_multi,/border,/raw,fixed=64,datadiffs,$
        line=3,/log,title='difference',res=res,/ps,file='r.ps'
  end

  save,dataset1,dataset2,datadiffs,file='rot_triotest.sav'



; We removed transform_volume because it was slow (8x cuberot3d)
; but had the same inaccuracies as cuberot3d.
;  starttime=systime(/seconds)
;  for i=0,n_elements(rhoset)-1 do begin

;      xang=0.0                  ; always 0
;      yang=phiset[i]                  ; height angle above plane
;      zang=0.0-rhoset[i]              ; CCW around z

;      image3=transform_volume(image3,rotation=[0.0,0.0,0.0-zang])
;      image3=transform_volume(image3,rotation=[yang,0.0,0.0])

;  end
;  print,'transform_volume ran in ',systime(/seconds)-starttime,' seconds'

;  threeview,fixed=64,/raw,image3,/log,title='transform_volume'


end

pro ignore

; very odd behavior, hard to parameterize, skip
image4=timage
interpolate3d,image4,[xang,yang,zang]
threeview,fixed=64,/raw,image4,/log,title='interpolate3d'

; did not work
;rotme=rot3dmatrix([xang,yang,zang]*!dtor)
;image5 = rotme ## timage
;threeview,fixed=64,/raw,image5,/log,title='rot3dmatrix'

end