pro p_insert_retro,pxwrapper,datum,oPixon,iele=iele,$
                  sigma=sigma,poisson=poisson,masks=masks,nonorm=nonorm,$
                  autoscale=autoscale,retro=retro,image=image
;
;+
; $Id: p_insert_retro.pro,v 1.1 2009/04/22 18:31:32 antunes Exp $
;
; Project   : developer area for STEREO SECCHI
;                   
; Name      : p_insert_retro
;             
; Purpose   : inserts datum into a pixon object, optionally also
;             inserting a provided noise (sigma)
;               
; Explanation: 
;              
;               
; Use       : IDL> p_insert_data,pxwrapper,datum,oPixon
;    
; Inputs    : 
;               
; Outputs   : 
;
; Keywords  : 
;
; Calls from LASCO : none 
;
; Common    : none 
;               
; Restrictions: none
;               
; Side effects: note this routine modifies three of its inputs as outputs
;               
; Category    : Simulation
;               
; Prev. Hist. : See Paul Reiser's code and notes.
;
; Written     : Sandy Antunes, NRL, Nov-Dec 2005
;               
;
; if given (iele), inserts the ptr into the object at that
; position, otherwise the default is to insert the dataset
; as a whole.
; When inserting a single item, pass only the single datum, not
; an array of datum.  When inserting a set, obviously you have
; to pass in an array of datum.
; Note the pixon object must already exist and be of the
; appropriate size-- ptr array of data has proper number of elements.
;
; Note that, in general, we multiply our input data by sr.normd
; before inserting.  However, if you set /nonorm, this will be skipped.
;
; /retro uses the earlier Reiser method, which does not allow for pointers
;

if (n_elements(iele) eq 0) then iele=-1
nonorm=KEYWORD_SET(nonorm)
retro=KEYWORD_SET(retro)
autoscale=KEYWORD_SET(autoscale)

if (iele eq -1) then single=0 else single=1

if (n_elements(image) ne 0) then force_precision,image,pxwrapper.sc.doubleprec
force_precision,datum,pxwrapper.sc.doubleprec

; fragile routine here, for backwards compatibility with earlier Reiser sims
if (retro) then begin
  Nd=n_elements(datum[0,0,*])
  if (n_elements(image) ne 0) then $
    oPixon->Set,image=image ; insert starter image
  if (n_elements(sigma) ne 0) then oPixon->Set,sigma=sigma
  oData = oPixon->get(/oData)
  for m=0,Nd-1 do oData[m]->Set, data=PTR_NEW(datum[*,*,m])
  return; early exit
end

; this is here just for compatability with earlier test cases,
; where both data and image were set together.
if (n_elements(image) ne 0) then $
  oPixon->Setimage,image=image  ; insert starter image

if (iele eq -1) then Nd=(size(datum))[3] else Nd=n_elements(pxwrapper.sr.gamma)

if (Nd ne n_elements(pxwrapper.sr.gamma)) then begin
    print,'Insufficient geometrical info for dataset, exiting p_insert_data'
    stop
end

pixontype=pxwrapper.object


; Note potential memory leak, as we do not free the old pointer before
; assigning this anew.  Since we typically insert data just once and
; at the beginning, this should rarely actually cause leaks.
; Even if done several time, leak is small relative to
; other data use and thus we do not code around it.

; This bit looks up the desired Pixon range pxwrapper.sp.autorange
; and then sets the normalization constant pxwrapper.sr.normd
; so that the data will be rescaled to fit that range.  Best
; used when ingesting a set of data at a time, if used when
; ingesting individual files, you may end up with idiosyncratic
; scaling, so we disable that.  You can always just run
;  p_autorange,datum,pxwrapper   manually to set the scaling constants.
; Note p_autorange only sets the normalization constant, it does not
; apply it.  We apply it here in p_insert_data (unless /nonorm is set).
;

if (autoscale) then begin
    if (single) then begin
        print,'Warning, will not set autoscale with single data, skipping autoscaling'
    end else begin
        p_autorange,datum,pxwrapper
    end
end

if (nonorm) then normd=1.0 else normd=pxwrapper.sr.normd

print,'autoscale is ',autoscale,' and norm is ',normd; debug

; just ensure they are pointers and move on
; Note that 'force_ptrs' can apply a normalization that ends up
;   applying to the returned data but does not alter the original
;   data, so in general 'force_ptrs' actually returns a copy of
;   the original data, for safety.
;

if (iele ne -1) then begin
    dptr=force_ptrs(datum,1,/single,norm=normd)
end else begin
    dptr=force_ptrs(datum,Nd,norm=normd)
end

;help,dptr; debug

if (n_elements(sigma) ne 0) then begin

    if (iele eq -1) then begin

        ; easy case, put them all in
;sdebug        oPixon->Set,sigma=force_ptrs(sigma,Nd); fails
;        print,'*** Doing sigma'; debug
;        help,sigma

;; used to work but now does not
;;;        sigma=normd*sigma
;;;        print,'max sigma is ',max(sigma)
;;;        oPixon->Set,sigma=sigma ; works

;;;;;  worked as sigma, but discarded to use new wt method instead
;;;;;        oPixon->Set,sigma=force_ptrs(sigma,norm=normd)

        psigma=force_ptrs(sigma,norm=normd)
        nsig=n_elements(psigma)
        pmasks=ptrarr(nsig)
        if (n_elements(masks) eq 0) then begin
            for m=0,nsig-1 do begin
                N=pxwrapper.sr.ndata[m]
                pmasks[m]=ptr_new(make_array(/float,N,N,nsig,value=1.0))
            end
        end else begin
            pmasks=force_ptrs(masks)
        end
        
        wt=ptrarr(nsig)
        for m=0,nsig-1 do begin
            wt[m] = ptr_new((*(pmasks[m])) / ((*(psigma[m])) * normd))
        end
        oPixon->Set,wt=wt

    end else begin

        ; single element, so extract set, update, and replace
        ; note we always save sigma as pointers, even if it is not
        ; yet one already
        tempset=force_ptrs(oPixon->Get(/wt))
        if (n_elements(masks) eq 0) then begin
            N=pxwrapper.sr.ndata[iele]
            mask=ptr_new(make_array(/float,N,N,value=1.0))
        end else begin
            ; can ingest the full mask set then select, or handle a solo mask
            mask=force_ptrs(masks)
            if (n_elements(mask) ne 1) then mask=mask[iele]
        end
        
        wt = ptr_new ((*(mask)) / ((*(tempset[iele])) * normd))
        tempset[iele]=wt
        oPixon->Set,wt=tempset

    end

end else begin
  ; if no sigma given, set to default of 1
  oPixon->Set,wt=1.0/normd
end


if (pixontype eq 't') then begin

    data  = PTRARR(1,Nd)        ; TETRA
    if (iele ne -1) then begin
        ; first copy existing data, even if blank
        odata=oPixon->Get(/data)
        for m=0,Nd-1 do data[0,m]=odata[0,m]

        ; now insert the one datum
        data[0,iele]=dptr

        ns=n_elements((*dptr)[0,*])
        pxwrapper.sr.ndata[iele]=ns; update its size
    end else begin
        ; insert the whole set
;        for m=0,Nd-1 do data[0,m]=ptr_new(datum[*,*,m]) ; TETRA
        for m=0,Nd-1 do data[0,m]=dptr[m]
    end
    oPixon->Replace,data=data  ; TETRA

end else if (pixontype eq 'c') then begin

    oData = oPixon->get(/oData)
    if (iele ne -1) then begin
        ; insert one datum
        ns=n_elements((*dptr)[0,*])
        pxwrapper.sr.ndata[iele]=ns; update its size
;        ratid = float(pxwrapper.sr.Ncube)/float(pxwrapper.sr.Ndata[iele])
;        ratid = float(pxwrapper.sr.Ndata[iele])/float(pxwrapper.sr.Ncube)
;        print,'setting ratio as ',ratid,pxwrapper.sr.Ncube,pxwrapper.sr.Ndata[iele]
;        oData[iele]->Set, data=PTR_NEW(datum)
;        oData[iele]->Set, ratid=ptr_new(ratid)
        oData[iele]->Set, data=dptr
    end else begin
        ; insert the whole set
;        help,oData; debug
;        print,'*** Doing datum'
;        for m=0,Nd-1 do oData[m]->Set, data=PTR_NEW(datum[*,*,m]); works
;        for m=0,Nd-1 do oData[m]->Set, data=ptr_new(*dptr[m]); works
        for m=0,Nd-1 do oData[m]->Set, data=dptr[m]; works, and simplest
        for m=0,Nd-1 do print,'max data is ',max(*dptr[m])
    end

end else begin

    print,'Cannot determine pixon version, exiting p_insert_data'
    stop

end

if (n_elements(poisson) ne 0) then begin

    if (iele eq -1) then begin

        ; easy case, put them all in
        oPixon->Replace,poisson=force_ptrs(poisson,Nd)

    end else begin

        ; single element, so extract set, update, and replace
        ; note we always save poisson as pointers, even if it is not
        ; yet one already
        tempset=force_ptrs(oPixon->Get(/poisson))
        tempset[iele]=force_ptrs(poisson,1,/single)
        oPixon->Replace,poisson=tempset

    end

end


END