Keck-DataReductionPipelines / MosfireDRP

http://keck-datareductionpipelines.github.io/MosfireDRP
10 stars 13 forks source link

High error values in the final rectified spectrum: bad pixel interpolation of the sigma spectrum is not handled correctly? #42

Open themiyan opened 8 years ago

themiyan commented 8 years ago

I think there is something wrong in the way how bad pixel interpolation is handled during rectification. In the standard deviation spectrum for first order we expect the error of a bad pixel to be ~2X higher than the surrounding pixels. Even with rectification this should be approximately correct for pixels around the bad pixel.

In the attached image, the white dots in the top left frame corresponds to bad pixels. The value between a white pixel and the neighbouring pixel can vary even more than X10 times! This should not be correct?

It seems this occurs from the interpolation step, for which I've shown the code below. Using r_interpol() function works fine for the eps and itime frames, but I think not so much for the ivar (perhaps because it use 1d interpolation). When I change code to interpolate the variance or the standard deviation instead ( as shown in the other two frames in the screenshot) these high peak values disappear. Now the values are what we expect it to be (~2 times higher) and the bad pixels outside the common dithering space has a value of 0 (These are the dark black pixels in the top right and bottom left frames).

This may be related to issue https://github.com/Keck-DataReductionPipelines/MosfireDRP/issues/35 as well.

screenshot 2016-02-04 22 46 46

for shift in shifts:
        output = r_interpol(ll, eps, fidl, tops, top, shift_pix=shift/0.18,
            pad=[mnshift, mxshift], fill_value = np.nan)
        epss.append(sign * output)

        ivar = 1/vv
        bad = np.where(np.isfinite(ivar) ==0)
        ivar[bad] = 0.0
        #TN: change code here to interpolate var or sd and not ivar?
        output = r_interpol(ll,ivar, fidl, tops, top, shift_pix=shift/0.18,
            pad=[mnshift, mxshift], fill_value=np.nan) 
        ivss.append(output)

        output = r_interpol(ll, it, fidl, tops, top, shift_pix=shift/0.18,
            pad=[mnshift, mxshift], fill_value=np.nan) 
        itss.append(output)

        sign *= -1

    it_img = np.nansum(np.array(itss), axis=0)
    eps_img = np.nanmean(epss, axis=0)

    # Remove any NaNs or infs from the variance array
    ivar_img = []
    for ivar in ivss:
        bad = np.where(np.isfinite(ivar) == 0)
        ivar[bad] = 0.0

        ivar_img.append(ivar)

    IV = np.array(ivar_img)
    bad = np.isclose(IV,0)
    IV[bad] = np.inf
    var_img = np.nanmean(1/np.array(IV), axis=0)
    sd_img = np.sqrt(var_img)
csteidel commented 8 years ago

This is actually a rather subtle problem, and one that we should think about carefully. By “bad pixel” I assume you mean a pixel that is in the bad pixel mask. do you mean a pixel that is in the bad pixel mask? Assuming that is the case, it is a very subtle problem what to do with pixels that (after rectification) that are partially affected by a “bad” pixel. Strictly speaking, one should probably use a “drizzle” type algorithm to make sure that the new pixel grid properly accounts (photometrically) for pixels that contain portions of a masked pixel in the original grid. The conservative approach would be to mask out all such pixels, but with only 2 spatial positions we cannot afford to do that. Anyway, what to do with the sigma (or variance or inverse variance) is a good question- I don’t think interpolation is the right thing to do, in any case- perhaps one should think of masked pixels as “NaNs” (rather than pixels having zero inverse variance) which should be completely ignored in creating the new (remapped) pixel grid. Since the units of the frames are e-/second at the point of rectifying (is that correct?), a final pixel should probably be ignored (assigned NaN?) if more than half of its area was inherited from a masked pixel, and assigned the value of the dominant “parent” pixel if less than half of its area was once a bad pixel. I must remind myself of the units being used at this point to be sure- it makes a difference whether the units are total electrons or electrons/sec…

On Feb 4, 2016, at 4:13 AM, Themiya Nanayakkara notifications@github.com wrote:

I think there is something wrong in the way how bad pixel interpolation is handled during rectification. In the standard deviation spectrum for first order we expect the error of a bad pixel to be ~2X higher than the surrounding pixels. Even with rectification this should be approximately correct for pixels around the bad pixel.

In the attached image, the white dots in the top left frame corresponds to bad pixels. The value between a white pixel and the neighbouring pixel can vary even more than X10 times! This should not be correct?

It seems this occurs from the interpolation step, for which I've shown the code below. Using r_interpol() function works fine for the eps and itime frames, but I think not so much for the ivar (perhaps because it use 1d interpolation). When I change code to interpolate the variance or the standard deviation instead ( as shown in the other two frames in the screenshot) these high peak values disappear. Now the values are what we expect it to be (~2 times higher) and the bad pixels outside the common dithering space has a value of 0 (These are the dark black pixels in the top right and bottom left frames).

This may be related to issue #35 https://github.com/Keck-DataReductionPipelines/MosfireDRP/issues/35 as well.

https://cloud.githubusercontent.com/assets/5274688/12814178/3bae2d36-cb91-11e5-8a2f-935a9749d68c.png for shift in shifts: output = r_interpol(ll, eps, fidl, tops, top, shift_pix=shift/0.18, pad=[mnshift, mxshift], fill_value = np.nan) epss.append(sign * output)

    ivar = 1/vv
    bad = np.where(np.isfinite(ivar) ==0)
    ivar[bad] = 0.0
    #TN: change code here to interpolate var or sd and not ivar?
    output = r_interpol(ll,ivar, fidl, tops, top, shift_pix=shift/0.18,
        pad=[mnshift, mxshift], fill_value=np.nan) 
    ivss.append(output)

    output = r_interpol(ll, it, fidl, tops, top, shift_pix=shift/0.18,
        pad=[mnshift, mxshift], fill_value=np.nan) 
    itss.append(output)

    sign *= -1

it_img = np.nansum(np.array(itss), axis=0)
eps_img = np.nanmean(epss, axis=0)

# Remove any NaNs or infs from the variance array
ivar_img = []
for ivar in ivss:
    bad = np.where(np.isfinite(ivar) == 0)
    ivar[bad] = 0.0

    ivar_img.append(ivar)

IV = np.array(ivar_img)
bad = np.isclose(IV,0)
IV[bad] = np.inf
var_img = np.nanmean(1/np.array(IV), axis=0)
sd_img = np.sqrt(var_img)

— Reply to this email directly or view it on GitHub https://github.com/Keck-DataReductionPipelines/MosfireDRP/issues/42.

themiyan commented 8 years ago

This is actually a rather subtle problem, and one that we should think about carefully. By “bad pixel” I assume you mean a pixel that is in the bad pixel mask. do you mean a pixel that is in the bad pixel mask?

yes, thats right.

I don’t think interpolation is the right thing to do, in any case- perhaps one should think of masked pixels as “NaNs” (rather than pixels having zero inverse variance) which should be completely ignored in creating the new (remapped) pixel grid. Since the units of the frames are e-/second at the point of rectifying (is that correct?), a final pixel should probably be ignored (assigned NaN?) if more than half of its area was inherited from a masked pixel, and assigned the value of the dominant “parent” pixel if less than half of its area was once a bad pixel.

Yes, this is a more reasonable thing to do, as long as the exposure frames, itime and sd is treated similarly. However at the moment any NAN or INF error is assigned a value of 0 by the code.

The units of the exposure frames at this point is in e-/s but the standard deviation frames are in total e-.