;+
;
; Project   : Densitycube, reconstruction, 3D, and s3drs software
;                   
; Name      : simple3drv
;               
; Purpose   : used to manipulate 3D densitycubes and do 3D reconstructions
;               
; Explanation: sloppy little test routine to do a simple density sum
; comparison of view projections to get a zeroth-order
; result for the possible underlying reconstruction image.
; Optionally can return the 'timage' for processing

; Use       : standalone or via s3drs GUI
;    
; Inputs    : (many)
;               
; Outputs   : returns the solution image
;
; Keywords  : (many)
;               
; Common    : none
;               
; Restrictions: none
;               
; Side effects: 
;               
; Category    : 3D, reconstruction, density
;               
; Prev. Hist. : None.
;
; Written     : Sandy Antunes, NRC, March-April 2009
;               
;-            
; sloppy little test routine to do a simple density sum
; comparison of view projections to get a zeroth-order
; result for the possible underlying reconstruction image.
; Optionally can return the 'timage' for processing
;
; Called with no arguments, it does a simple 32x32x32 2-view demo routine
; via data generated by RTWL
;
; If given an 'rsun' in pixels (aka voxels), it will put a solid
; sun of that size and apply it as an LOS constraint --- NOT YET TESTED!!!
;
; If given the '/usemax' flag, it will enforce visual maxima, i.e. it
; will do a better reconstruction by constraining the amount in to
; the least amount from all contributing LOS.
;
; Can optionally take a Pixon 'pxwrapper' instead of a set of rotvects
;
; If given /backproject, does more accurate Thomson backprojection
; instead of summation.  Or, give it a mass image instead.
;
; /norm normalizes final solution to match data levels, recommended,
; but requires the 'pxwrapper' input structure.
;
; /carve carves out (removes) center sun, again needs 'pxwrapper' to
; work, recommended
;
; Added 'pxfile=file' so you can skip giving any inputs if you have a
; Pixon-format save file ready.
;
; Optional '/display' and '/contour' will plot results.
;

FUNCTION simple3drv,datum,rotvects,sloppy=sloppy,large=large,norm=norm,$
                    N=N,debug=debug,rsun=rsun,usemax=usemax,carve=carve,$
                    pxwrapper=pxwrapper,backproject=backproject,pxfile=pxfile,$
                    display=display,contour=contour,save=save,stem=stem

; ### CHECK ARGUMENTS

  sloppy=KEYWORD_SET(sloppy); choice of algorithm
  large=KEYWORD_SET(large); experimental large array size method
  display=KEYWORD_SET(display);
  contour=KEYWORD_SET(contour);
  debug=KEYWORD_SET(debug)
  usemax=KEYWORD_SET(usemax)
  backproject=KEYWORD_SET(backproject)
  norm=KEYWORD_SET(norm)
  carve=KEYWORD_SET(carve)
  save=KEYWORD_SET(save)
  if (n_elements(rsun) eq 0) then rsun=-1
  if (n_elements(stem) eq 0) then stem=''
  hasmax=KEYWORD_SET(hasmax); UNUSED SO FAR

  if (n_elements(pxfile) ne 0) then begin
      restore,/v,pxfile
      datum=diff*mset
  end

  if (n_elements(datum) eq 0) then begin

      if (n_elements(N) eq 0) then N=32
      fakedensitydata,datum,rotvects,N=N

  end else begin

      if (n_elements(rotvects) eq 0 and n_elements(pxwrapper) eq 0) then begin
          print,'Must include geometry (3-vect rotvects) of views, aborting.'
          return,fltarr(0,0,0)
      end

  end

; ### PREPARE VARIABLES

  N = (size(datum))[1]
  nviews = (size(datum))[3]

  if (n_elements(pxwrapper) ne 0) then begin
      ; extract rotvects from pxwrapper
      phi=pxwrapper.sr.phi
      rho=pxwrapper.sr.rho
      rotvects=fltarr(3,nviews)
      stem=stem+pxwrapper.sc.gname
      for m=0,nviews-1 do rotvects[*,m] = [0,0.0-phi[m],rho[m]]
  end

  if (n_elements(rotvects) eq nviews) then begin
      ; special munging here, if we just give it one angle (presumed
      ; position along ecliptic) instead of a proper 3-angle for each
      ; view.  We alter the 'rotvects' to then have our proper form.
      ; useful only for testing.
      rho=rotvects
      rotvects=fltarr(3,nviews)
      for m=0,nviews-1 do rotvects[*,m] = [0,0,rho[m]]
  end


  ; NOW ALLOCATED HELPERS

  ; create each solution cube as simple average projection of data
  ; and then rotate to co-align and sum
  if (debug) then print,'allocating 2 large arrays of size ',N

  if (large) then begin

      sumimage = ptrarr(N, /allocate_heap)
      for m=0,N-1 do *sumimage[m]=fltarr(N,N)
      stemp = ptrarr(N, /allocate_heap)
      for m=0,N-1 do *stemp[m]=fltarr(N,N)

  end else begin

      sumimage=fltarr(N,N,N)
      stemp=fltarr(N,N,N)

  end

  if (debug) then print,'allocation done.'

; ### NOW LOOP THROUGH VIEWS AND CREATE SOLUTION

  for i=0,nviews-1 do begin
      if (debug) then print,'starting view ',i+1

      if (rsun ne -1) then begin
          ; we force the 'large' style for convenience later
          cube_los,N,rsun,los,loscount,large=large
      end else begin
          ; no LOS to worry about
          cube_los,N,rsun,los,loscount,large=large,/null
          loscount=float(N-1); just a scalar
      end

      if (backproject) then begin

          ; to fit with our summation scheme, we remove any
          ; rotational bits but keep the other spatial info
          pxtemp=pxwrapper
          pxtemp.sr.gamma[i]=0.0
          pxtemp.sr.rho[i]=0.0
          pxtemp.sr.theta[i]=0.0
          pxtemp.sr.phi[i]=0.0
          stemp=pr_backproject(pxtemp,datum[*,*,i],i)

      end else begin
          ; simple summation
          dtemp=datum(*,*,i)/loscount ; appropriate average

          if (debug) then print,'projecting data into temp image'

          if (large) then begin
              for m=0,N-1 do *stemp[m]=dtemp * (*los[m])
          end else begin
              for j=0,N-1 do stemp[j,*,*]=dtemp * los[j,*,*]
          end

      end

;;      if (debug) then print,'reordering image'
;;      if (large) then begin
;;          ;convert large array ijk to pixon schema
;;          stemp=reverse(temporary(stemp))
;          stemp=rearrange(stemp,[1,2,-3]) ; convert ijk to pixon schema
;;          if (rsun ne -1) then los=reverse(temporary(los))
;;      end else begin
          ;convert ijk to pixon schema
;;          stemp=rearrange(temporary(stemp),[1,2,-3])
;;          if (rsun ne -1) then los=rearrange(temporary(los),[1,2,-3])
;;      end

      if (debug) then print,'rotating image'

      ; modified to handle regular or large arrhandle
      cuberot3d,stemp,rotvects[*,i] 

      if (usemax) then $
        cuberot3d,los,rotvects[*,i]

      if (debug) then $
        threeview,stemp,rescale=128/N,title='frame '+int2str(i),$
        /pixon,/log

;      print,rotvects[*,i]
;      threeview,stemp

      if (debug) then print,'doing summation reconstruction'

      if (i eq 0) then begin
          ; first solution is of course unbounded and sets initial state

          if (large) then begin
              for m=0,N-1 do *sumimage[m] = *stemp[m]
          end else begin
              sumimage = stemp
          end

      end else if (sloppy) then begin
          ; sloppy 'summation of all possible solutions' method
          ; (alternately, prepares 'usemax' during first iteration
          ;  or prepares 'sum of most restrictive solutions' first iteration

          if (large) then begin
              for m=0,N-1 do *sumimage[m] = *sumimage[m] + *stemp[m]
          end else begin
              sumimage = temporary(sumimage) + stemp
          end

      end else if (usemax) then begin
          ; make a cube just giving the max value possible for
          ; each voxel, where max = minimum full density based
          ; on multiple views

          if (large) then begin
              for m=0,N-1 do *sumimage[m] = *sumimage[m] > *stemp[m]
          end else begin
;              c=where (los gt 0)
;              if (c[0] ne -1) then sumimage = min(sumimage[c],stemp[c])
              if (debug) then $
                threeview,title='sum'+int2str(i),/pixon,sumimage,/log,fixed=128
              if (debug) then $
                threeview,title='stemp'+int2str(i),/pixon,stemp,/log,fixed=128
              sumimage = stemp < sumimage
              if (debug) then $
                threeview,title='min'+int2str(i),/pixon,sumimage,/log,fixed=128
          end

      end else if (hasmax) then begin
          ; we have a max cube from earlier, so we can now solve


      end else begin
          ; mildly more accurate 'sum of most restrictive solutions'
          ; or summation only where there is intersection

          if (debug) then print,'(adding to solution)'

          if (large) then begin
              for m=0,N-1 do $
                *sumimage[m] = imageintersection(*sumimage[m],*stemp[m])
          end else begin
              sumimage = imageintersection(temporary(sumimage),stemp)
          end

      end

  end

  ; at this point, we have either a 'max' image or a 'sumimage' solution.
  if (usemax) then begin
    ; having made our max, we now recursively call to make a solution

  end

  if (debug) then print,'Done solution!  Now normalizing.  Sumimage info:'

; ### NORMALIZE

  if (debug) then help,sumimage

  if (large) then begin
      for m=0,N-1 do *sumimage[m]=*sumimage[m]/float(nviews)
  ; and normalize
      localmax=fltarr(N)
      for m=0,N-1 do localmax[m] = max(*sumimage[m])
      summax=max(localmax)
;      print,'max is ',summax
      for m=0,N-1 do *sumimage[m] = *sumimage[m]/summax
  end else begin
      sumimage=temporary(sumimage)/float(nviews)
  ; and normalize
      smax=max(sumimage)
      sumimage = temporary(sumimage)/smax
  end

  if (norm) then begin
      ; normalize to data levels, recommended
      if (n_elements(pxwrapper) eq 0) then begin
          print,'No pxwrapper, cannot normalize to data levels, sorry'
      end else begin
          normfactor=max(datum)/max(pr_render(pxwrapper,sumimage))
          if (large) then begin
              for m=0,N-1 do *sumimage[m]=*sumimage[m]/normfactor
          end else begin
              sumimage = sumimage * normfactor
          end
      end
  end

  if (carve) then begin
      ; removes sun-shaped region from solution, recommended
      if (n_elements(pxwrapper) eq 0) then begin
          print,'No pxwrapper, cannot carve out sun region, sorry'
      end else begin
          m=p_getmask('COR2',pxwrapper.sr.ncube,rpixels=rpixels)
          mask3D=cube_sphere(pxwrapper.sr.ncube,rpixels[0],r_outer=rpixels[1])
          if (large) then begin
              for m=0,N-1 do *sumimage[m]=*sumimage[m]/mask3D[*,*,m]
          end else begin
              sumimage = sumimage * mask3D
          end
      end
  end

; ### VISUALIZE

  if (display or contour) then begin
      threeview,sumimage,/pixon,oom=2,trio,/noplot
      zview=trio[*,*,2]
      if (contour) then begin
          image_cont,zview,/window_scale
          write_gif,stem+'.gif',tvrd()
      end
      if (display) then begin
          tv_multi,datum,rescale=128/N,title='original data',/log
          threeview,sumimage,rescale=128/N,title='simple reconstruction',$
            /pixon,/log
      end
  end else begin
      ;; force write/save of contours
      ;threeview,sumimage,/pixon,oom=2,trio,/noplot
      ;zview=trio[*,*,2]
      ;image_cont,zview
      ;write_gif,stem+'.gif',tvrd()
  end

  if (save) then begin
      threeview,sumimage,/pixon,oom=2,trio,/noplot
      zview=trio[*,*,2]
      save,file='simple_'+stem+'.sav',datum,sumimage,zview,trio
  end

;  if (debug) then help,sumimage

  if (large) then ptr_free,stemp
  return,sumimage

END