;+
; 
; Name: thm_part_slice2d_2di.pro
; 
; Purpose:  Helper function for thm_part_slice2d.pro
;           Produces slice by interpolating the volume
;           in three dimensions then extracting a slice.
;
;-
pro thm_part_slice2d_3di, datapoints, velocities, resolution, $
                          orange=orange, displacement=displacement, $
                          part_slice=part_slice, $
                          xgrid=xgrid, ygrid=ygrid, $
                          fail=fail

    compile_opt idl2, hidden


  ;Error checks
  normal = [0,0,1.]
  xvec = [1.,0,0]

  ; Create cube grid
  velmm = [ [minmax(velocities[*,0])], [minmax(velocities[*,1])], [minmax(velocities[*,1])] ]
  xgrid = interpol(velmm[*,0], resolution)
  ygrid = interpol(velmm[*,1], resolution)
  zgrid = interpol(velmm[*,2], resolution)

  ; Get slice's center point
  center = ( (normal * displacement)/(velmm[1,*]-velmm[0,*]) + 0.5) * resolution
  if in_set(center le 0 or center gt resolution, 1) then begin
    fail = 'Slice displacement is outside the data range.'
    return
  endif  

  ; Must be copied to new variables for qhull
  xvels = velocities[*,0]
  yvels = velocities[*,1]
  zvels = temporary(velocities[*,2])
  
  qhull, xvels, yvels, zvels, th, /DELAUNAY
  
  ; Remove tetrahedra whose total velocity (centroid) is less than
  ; minimum velocity from distribution (prevents interpolation over
  ; lower energy limits)
  index = where( 1./16 * total(  xvels[th[0:3,*]] ,1 )^2 + $
                 1./16 * total(  yvels[th[0:3,*]] ,1 )^2 + $
                 1./16 * total(  zvels[th[0:3,*]] ,1 )^2  $
                  gt min(xvels^2+yvels^2+zvels^2), $
                  count, ncomplement=ncomp)
  if count gt 0 then begin
    if ncomp gt 0 then begin
      th=th[*,index]
    endif
  endif else begin
    fail = 'Range error in calculating triangulation. Cannot interpolate data.'
    return
  endelse 

  ; Interpolate data to regular 3D grid
  vol = qgrid3(xvels, yvels, zvels, datapoints, th, dimension=replicate(resolution,3))
  
  ; Remove erroneous data points (also helps prevent interpolation over gaps)
  derp = where(abs(vol) lt orange[0] or abs(vol) gt orange[1], nderp)
  if nderp gt 0 then begin
    vol[derp] = 0
  endif 

  ; Extract slice from regular grid
  part_slice=extract_slice(vol, resolution, resolution, $
                           center[0], center[1], center[2], $
                           normal, xvec, /sample)

;  ; For testing only: shows 3-D distribution of non-zero data (using valid scaling)
;  ww = where(vol gt 0, nww)
;  tmparr = bytscl(alog10(vol[ww]))
;  x_idxs = where(   vol[ww] eq min( (vol[ww])[where(vol[ww] ge orange[0])] )   )
;  x = float( tmparr[x_idxs[0]] )
;  m = -255./(x-255)
;  c = (255.*x)/(x-255)
;  tmparr = m*tmparr + c
;  ltz = where(tmparr lt 0, nltz)
;  if nltz gt 0 then tmparr[ltz] = 0
;  ;tmparr[where(tmparr gt 255)] = 255
;  tmparr = byte(round(tmparr))
;  x_idx = ww mod n_elements(xgrid) 
;  y_idx = (ww / n_elements(xgrid)) mod n_elements(ygrid)
;  z_idx = ww / (n_elements(xgrid)*n_elements(ygrid))
;  rg2 = [-max(xgrid),max(xgrid)]
;  iplot, xgrid[x_idx], ygrid[y_idx], zgrid[z_idx], linestyle=6, sym_index=1, $
;         xrange = rg2, yrange=rg2, zrange=rg2, $
;         rgb_table=13, vert_colors=tmparr


;  ; For testing only: View volume in IDL medical imaging software
;  mm = minmax(vol, min_value=0)
;  vol= temporary(vol)/mm[0]
;  vol = alog10(temporary(vol))
;  mm = minmax(vol, min_value=0)
;  vol = BYTSCL(temporary(vol), min=mm[0], max=mm[1])
;  hData = PTR_NEW(vol) ; slicer3 needs a pointer to the volume data 
;  slicer3, hData

end