arguments
science/weight list
number of iteration
exponent
number of input science images
medname (None at first iteration)
read an input file consisting of:
science1 var1
science2 var2
etc...
Read all sciences and weights
replace all 0 in var(n) with 1.0E6
sigma = sqrt ( var )
var_sum += var
Calculate a starting mean image
foreach pixel (looking toward the stack)
if var_sum >= 1.0E6
dummy = science where sig_stack != 1.0E3 (sqrt( 1.0E6 ))
dummy_wgt = 1.0/(( sig_stack != 1.0E3 )**2)
else
dummy = science
dummy_wgt = 1.0/( sig_stack**2)
if dimension of dummy > 0 and <= 2
mean_image = min ( dummy )
else if dimension of dummy > 2
dummy_s = sort ( dummy)
dummy_wgt_s = dummy_wgt sorted like dummy
then find the mid point (where weight is 50% of total weight ) in the sorted array
(could use a linear interpolation if mid point is not an exact image number)
mean_image = mid_point
med_image = mean_image
Apply skeptical algorithm
for each iteration
apply equation
wgt_stack = choose(equal(sig_stack,1.e3),((1.0/(sig_stack**2))*(1.0/(1.0+((absolute(sci_stack-mean
_image)/sig_stack)**hexpon))), 0.0))
# mean image for current iteration is calculated here. Pixels with no valid values
# throughout the entire stack have a weight sum of zero. This zero weight sum
# would cause a division by zero error in the calculation for the mean. An
# appropriate trap is therefore set here as well, and bad pixel values are set to -1.0e6
num = sum(sci_stack*wgt_stack)
denum = sum(wgt_stack)
denum2 = where(equal(denum,0.0),1.0,denum)
mean_image = num/denum2
mean_image = where(equal(denum,0.0),-1.0e6,mean_image)
final combined weighted-average image (Computing final weighted combined image)
apply equation using the mean_image from last iteration
wgt_stack = choose(equal(sig_stack,1.e3),((1.0/(sig_stack**2))*(1.0/(1.0+((absolute(sci_stack-mean_image)/
sig_stack)**hexpon))), 0.0))
num = sum(sci_stack*wgt_stack)
denum = sum(wgt_stack)
denum2 = where(equal(denum,0.0),1.0,denum)
comb_image = num/denum2
comb_image = where(equal(denum,0.0),-1.0e6,comb_image)
final combined output weighted-variance image (Computing final weighted combined sigma image )
recompute the weight
num = sum((wgt_stack*sig_stack)**2)
denum = sum(wgt_stack)**2
denum2 = where(equal(denum,0.0),1.0,denum)
comb_image_var = num/denum2
comb_image_var = where(equal(denum,0.0),1.0e6,comb_image_var)
Use lanczos3 2D interpolation to remove bad pixels and flag them
nonetheless in weight image interpolation area is 9x9. The bad pixel
interpolation is use for cosmetic purpose only, the weght stays at 1.0e06
|
(Subscribe to new public observations)
