;+
;
; 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
;               
;-            
; given data and view angles, produces a 'probability cube'
; for the likely underlying distribution, using this method:
;
; Back project each data into a cube, then
; for each voxel, return the minimum value from each back projection
;
; This basically indicates the _brightest_ any single voxel can be,
; without 
;
; The default is a simple summing, which therefore returns a magnitude
; but not actual electron density values.
; The /backproject flag uses Paul's Thomson scattering back projection
; and does return electron density values.
; Assumes the back project function can return such a 'worst case
; solution'; if the back project function returns an averaged
; n_e density, this method is still valid but does not produce a
; limiting case, only a statistically likely average representation.
;

FUNCTION simpledata2cube,pxwrapper,datum,index,image=image

  ; simple takes an image and (in the Pixon axes convention) sums
  ; it along the cube to make a pseudo-backprojection, rotating
  ; the resulting image cube so the answer returned can be
  ; compared to projections of other angles.

  N=pxwrapper.sr.Ncube

  theta=pxwrapper.sr.theta[index]
  gamma=pxwrapper.sr.rho[index]
  fast=0;pxwrapper.sc.fast

  if (n_elements(image) eq 0) then image=fltarr(N,N,N) else image=image*0.0
  for i=0,N-1 do for j=0,N-1 do image[*,j,i]=datum[i,j,index]
  if (fast) then begin
      cuberot3d,image,[0.0,0.0-theta,0.0]
      cuberot3d,image,[0.0,0.0,gamma]
  end else begin
      image=turn_3d(image,0.0,0.0-theta,0.0-gamma,/resize)
  end

  return,image

END


FUNCTION prob_cube,pxwrapper,datum,backproject=backproject,debug=debug,$
                   nframes=nframes

  backproject=setyesno(backproject)
  debug=setyesno(debug)

  N=pxwrapper.sr.Ncube

  if (n_elements(nframes) eq 0) then nframes=pxwrapper.sr.Nd

  tempcube=fltarr(N,N,N); allocate here to avoid pr_backproject re-allocs
  imgpair=fltarr(N,N,N,2); stores both working versions and the solution

  ; assemble possible solutions and min on the fly
  for iframe=0,nframes-1 do begin

      tempcube=tempcube*0.0
      if (backproject) then imgpair[*,*,*,1]=$
        pr_backproject(pxwrapper,datum,iframe,image=tempcube) else $
        imgpair[*,*,*,1]=simpledata2cube(pxwrapper,datum,iframe,image=tempcube)

      if (debug) then threeview,fixed=128,imgpair[*,*,*,1],$
        title='prob_cube '+strtrim(string(iframe),2),/autocap,oom=2,/pixon,/raw

      if (iframe eq 0) then begin
          imgpair[*,*,*,0]=imgpair[*,*,*,1]; store first frame guess as answer
      end else begin
          imgpair[*,*,*,0]=min(imgpair,dimension=4)
      end

  end

  return,imgpair[*,*,*,0]

END