Compute a deg=3 spline as overscan correction (closes #258)

megaradrp/processing/tests/test_trimover.py

@@ -1,9 +1,5 @@
 from megaradrp.processing.trimover import OverscanCorrector
-from numina.core import DataFrame
-import shutil
-from tempfile import mkdtemp
-import numpy as np
-import astropy.io.fits as fits
+import numpy
 
 def test_OverscanCorrector():
 
@@ -31,5 +27,36 @@ def test_OverscanCorrector():
     assert oc.orow2 == (slice(2106, 2156, None),slice(50, 4146, None))
 
 
+def test_amp_1():
+
+    detconf = {
+        'trim1': [[0, 2056], [50, 4146]],
+        'trim2': [[2156, 4212], [50, 4146]],
+        'bng': [1, 1],
+        'overscan1': [[0, 2056], [4149, 4196]],
+        'overscan2': [[2156, 4212], [0, 50]],
+        'prescan1': [[0, 2056], [0, 50]],
+        'prescan2': [[2156, 4212], [4145, 4196]],
+        'middle1': [[2056, 2106], [50, 4146]],
+        'middle2': [[2106, 2156], [50, 4146]]
+    }
+
+    oc = OverscanCorrector(detconf)
+
+    data = 2000.0 + numpy.zeros(shape=(4212, 4196))
+
+    # k,c,deg = spl1._eval_args
+
+    fit1, spl1 = oc.eval_spline_amp1(data)
+    data[oc.trim1] -= fit1[:, numpy.newaxis]
+
+    assert fit1.shape == (2056,)
+
+    fit2, spl2 = oc.eval_spline_amp2(data)
+    data[oc.trim2] -= fit2[:, numpy.newaxis]
+
+    assert fit2.shape == (2056,)
+
+
 if __name__ == "__main__":
     test_OverscanCorrector()

megaradrp/processing/trimover.py

@@ -14,9 +14,10 @@
 
 import numpy as np
 import numpy
-import numpy.polynomial.polynomial as nppol
+from scipy.interpolate import LSQUnivariateSpline
 from astropy.io import fits
-import numina.array.trace.extract as extract
+
+
 from numina.processing import Corrector
 
 
@@ -238,32 +239,68 @@ def run(self, img):
         oc1 = np.median(data[self.ocol1])
         _logger.debug('median col overscan1 is %f', oc1)
         # avg = (p1 + or1 + oc1) / 3.0
-        avg1 = oc1
-        _logger.debug('overscan1 correction is %f', avg1)
-        data[self.trim1] -= avg1
-
-        # p2 = data[self.pcol2].mean()
-        # _logger.debug('prescan2 is %f', p2)
-        # or2 = data[self.orow2].mean()
-        # _logger.debug('row overscan2 is %f', or2)
+
         oc2 = np.median(data[self.ocol2])
         _logger.debug('median col overscan2 is %f', oc2)
-        # avg = (p2 + or2 + oc2) / 3.0
-        avg2 = oc2
-        _logger.debug('overscan2 correction is %f', avg2)
-        data[self.trim2] -= avg2
+
+        _logger.debug('compute spline for overscan1')
+        fit1, spl1 = self.eval_spline_amp1(data)
+        data[self.trim1] -= fit1[:, numpy.newaxis]
+
+        _logger.debug('compute spline for overscan2')
+        fit2, spl2 = self.eval_spline_amp2(data)
+        data[self.trim2] -= fit2[:, numpy.newaxis]
+
         hdr = img['primary'].header
         hdr['NUM-OVPE'] = self.calibid
         hdr['history'] = 'Overscan correction with {}'.format(self.calibid)
         hdr['history'] = 'Overscan correction time {}'.format(datetime.datetime.utcnow().isoformat())
         # hdr['history'] = 'Mean of prescan1 is %f' % p1
         hdr['history'] = 'Median of col overscan1 is {}'.format(oc1)
-        hdr['history'] = 'Overscan1 correction is {}'.format(avg1)
+        hdr['history'] = 'Overscan1 correction is {}'.format('spline3')
         # hdr['history'] = 'prescan2 is %f' %  p2
         hdr['history'] = 'Median of col overscan2 is {}'.format(oc2)
-        hdr['history'] = 'Overscan2 correction is {}'.format(avg2)
+        hdr['history'] = 'Overscan2 correction is {}'.format('spline3')
+
+        for spl, label in zip([spl1, spl2], ['overscan1', 'overscan2']):
+            k, c, deg = spl._eval_args
+            hdr['history'] = '{} deg {}'.format(label, deg)
+            hdr['history'] = '{} knots {}'.format(label, list(k))
+            hdr['history'] = '{} coeffs {}'.format(label, list(c))
+
         return img
 
+    def fit_spline_amp1(self, data):
+
+        region = self.ocol1
+        u = numpy.arange(region[0].start, region[0].stop)
+        v = data[region].mean(axis=1)
+        knots1 = [1200]
+        spl1 = LSQUnivariateSpline(u, v, knots1, k=3)
+        return spl1
+
+    def eval_spline_amp1(self, data):
+        spl1 = self.fit_spline_amp1(data)
+        region = self.ocol1
+        u = numpy.arange(region[0].start, region[0].stop)
+        fit1 = spl1(u)
+        return fit1, spl1
+
+    def fit_spline_amp2(self, data):
+        region = self.ocol2
+        u = numpy.arange(region[0].start, region[0].stop)
+        v = data[region].mean(axis=1)
+        knots2 = [3100]
+        spl2 = LSQUnivariateSpline(u, v, knots2, k=3)
+        return spl2
+
+    def eval_spline_amp2(self, data):
+        spl2 = self.fit_spline_amp2(data)
+        region = self.ocol2
+        u = numpy.arange(region[0].start, region[0].stop)
+        fit2 = spl2(u)
+        return fit2, spl2
+
 
 class TrimImage(Corrector):
     """A Corrector Node that trims MEGARA images."""