; NAME:
;   interpix__init.por
;
; PURPOSE:
;   defines functions that are part of an interpix object
;
; CALLING SEQUENCE:
;   N/A - This entire file is compiled when an interpix
;   object is created.
;
; INPUTS:
;   N/A
;	
; KEYWORD PARAMETERS:
;   N/A
;
; OUTPUTS:
;   N/A
;
; EXAMPLE:
;   N/A
;
; MODIFICATION HISTORY:
;   Created by:
;       Ernie Morse, IDTL, October 7, 2005
;
Pro interpix::doPatchesInPair, im1, im2

    diffimg = im1 - im2
    sumimg = im1 + im2

    nrow = self.psize.w
    ncol = self.psize.h

    for i = 0, n_elements(*self.plocs) - 1 do begin
        c1 = (*self.plocs)[i].x
        c2 = c1 + ncol - 1
        r1 = (*self.plocs)[i].y
        r2 = r1 + ncol - 1
        *self.pdiff = diffimg[c1:c2, r1:r2]
        *self.psum = sumimg[c1:c2, r1:r2]
        if (self->goodPatch()) then begin
            self.p = self.p + 1
            self->correlatePatch
        endif
    endfor
End

Function interpix::goodPatch

    ; test the current patch
    ; if the current patch is OK, return 1
    ; missing thresholds means don't care.
    mthr = self.meanthresh
    varth = self.varthresh
    okmean = 0
    okvar = 0
    pdmean = mean(*self.pdiff)
    pdmin = min(*self.pdiff)
    pdmax = max(*self.pdiff)

    if (not mthr or (mthr and (-mthr lt pdmean and pdmean lt mthr))) then begin
       okmean = 1
    endif

    if (not varth or (varth and $
            (-varth lt pdmin and pdmin lt pdmax and pdmax lt varth))) then begin
       okvar = 1
    endif

    return, okmean and okvar
End

Pro interpix::correlatePatch

    ; Given a patch, add the center, vert, horiz and sig into the list.
    pdiff = *self.pdiff
    psum = *self.psum

    ; mean patch, squared, to offset light fluctuation
    pdmeansq = mean(pdiff)^2

    ; mean square noise pixel
    (*self.centers)[self.p] = mean(pdiff ^ 2) - pdmeansq

    for i = 0, n_elements(*self.checklist) - 1 do begin
        off1 = (*self.checklist)[i].off1
        off2 = (*self.checklist)[i].off2
        end1 = self.psize.w - (1 + off1)
        end2 = self.psize.h - (1 + off2)

        ; first order: off1=1, off2=0
        ; second order = 2,0, 1,0, 1,1
        ; for i,j = 1,1 that means row+i, col+j and row-j, col+i!!

        (*self.hors)[self.p, i] = mean(pdiff[0:end1, 0:end2] $
                * pdiff[off1:*, off2:*]) - pdmeansq

        (*self.vers)[self.p, i] = mean(pdiff[off2:*, 0:end1] $
                * pdiff[0:end2, off1:*]) - pdmeansq

    endfor

    bgo = self.order + 2
    bgend = self.psize.w - (bgo + 1)

    (*self.hbgs)[self.p] = mean(pdiff[0:bgend, *] * pdiff[bgo:*, *]) - pdmeansq
    (*self.vbgs)[self.p] = mean(pdiff[*, 0:bgend] * pdiff[*, bgo:*]) - pdmeansq

    (*self.means)[self.p] = mean(psum)

    ; next stuff for diagnostics only
    (*self.diffmeans)[self.p] = mean(pdiff)
    (*self.mins)[self.p] = min(pdiff)
    (*self.maxs)[self.p] = max(pdiff)

End


Function interpix::conversionFactor

    ; return an estimate of conversion factor, given accumulated data.
    ; Assume all data at the same light level
    signal = mean((*self.means)[0:self.p])
    uncompensated = mean((*self.centers)[0:self.p]) - self.readnoise
    noise_squared = uncompensated

    for i = 0, n_elements(*self.checklist) - 1 do begin
        noise_squared = noise_squared + 2 * mean((*self.hors)[0:self.p, i])
        noise_squared = noise_squared + 2 * mean((*self.vers)[0:self.p, i])

        ; remove background flutter by going a few pixels away.
        noise_squared = noise_squared - 2 * mean((*self.hbgs)[0:self.p])
        noise_squared = noise_squared - 2 * mean((*self.vbgs)[0:self.p])
    endfor

    return, [signal / uncompensated, signal / noise_squared]

End

Function interpix::conversionFactorFromSlope

    ; return an estimate of conversion factor, given accumulated data.
    ; Assume we have data at varying intensities, and fit noise-squared
    ; to signal with a curve.
    signal = (*self.means)[0:self.p]
    noise_squared = (*self.centers)[0:self.p]
    uncompensated = noise_squared

    for i = 0, n_elements(*self.checklist) - 1 do begin
        noise_squared = noise_squared + 2 * (*self.hors)[0:self.p, i]
        noise_squared = noise_squared + 2 * (*self.vers)[0:self.p, i]

        ; remove background flutter by going a few pixels away.
        noise_squared = noise_squared - 2 * (*self.hbgs)[0:self.p]
        noise_squared = noise_squared - 2 * (*self.vbgs)[0:self.p]
   endfor
   uslope = ( (mean(signal * uncompensated) - mean(signal) * mean(uncompensated)) $
           / (mean(signal * signal) - mean(signal)^2) )
   slope = ( (mean(signal * noise_squared) - mean(signal) * mean(noise_squared)) $
           / (mean(signal * signal) - mean(signal)^2) )

    return, [1 / uslope, 1 / slope]

End

Pro interpix::cleanup
    ptr_free, $
            self.centers, $
            self.means, $
            self.hbgs, $
            self.vbgs, $
            self.diffmeans, $
            self.mins, $
            self.maxs, $
            self.hors, $
            self.vers, $
            self.pdiff, $
            self.psum, $
            self.checklist
End

Function interpix::init, psize, plocs, numpairs, readnoise = readnoise, $
        meanthresh = meanthresh, varthresh = varthresh, $
        order = order

    ; Set default values for input
    if (not keyword_set(readnoise)) then begin
        readnoise = 0.0
    endif

    if (not keyword_set(order)) then begin
        order = 1
     endif

    if (not keyword_set(meanthresh)) then begin
        meanthresh = 0
     endif

    if (not keyword_set(varthresh)) then begin
        varthresh = 0
     endif

    n = n_elements(plocs) * numpairs

    ; Set object properties
    self.psize.w = psize[0]
    self.psize.h = psize[1]
    self.plocs = ptr_new(plocs)
    self.readnoise = 2 * (readnoise ^ 2)
    self.meanthresh = meanthresh
    self.varthresh = varthresh
    self.order = order
    self.p = -1l

    self.checklist = ptr_new(replicate({off1:0, off2:0}, $
            total(indgen(order + 1))))

    index = 0
    for i = 0, order do begin
        for j = 0, i - 1 do begin
            (*self.checklist)[index] = {off1:order + 1 - i, off2:j}
            index = index + 1
        endfor
    endfor

    self.centers = ptr_new(dblarr(n))
    self.means = ptr_new(dblarr(n))
    self.hbgs = ptr_new(dblarr(n))
    self.vbgs = ptr_new(dblarr(n))
    self.diffmeans = ptr_new(dblarr(n))
    self.mins = ptr_new(dblarr(n))
    self.maxs = ptr_new(dblarr(n))

    self.hors = ptr_new(dblarr(n, n_elements(*self.checklist)))
    self.vers = ptr_new(dblarr(n, n_elements(*self.checklist)))

    self.pdiff = ptr_new(dblarr(psize[0], psize[1]))
    self.psum = ptr_new(dblarr(psize[0], psize[1]))

    return, 1

End