Merge pull request #108 from starcksophie/deconvolve_binding

Deconvolve binding

pysap/extensions/sparse2d.py

@@ -29,14 +29,105 @@
 
 
 class Filter():
+    """ Define the structure that will be used to store the filter result.
+    """
     def __init__(self, **kwargs):
+        """ Define the filter.
+
+        Parameters
+        ----------
+
+        type_of_filtering: int
+        coef_detection_method: int
+        type_of_multiresolution_transform: float
+        type_of_filters: float
+        type_of_non_orthog_filters: float
+        sigma_noise: float
+        type_of_noise: int
+        number_of_scales: int
+        iter_max: double
+        epsilon: float
+        verbose: Boolean
+        tab_n_sigma: ndarray
+        suppress_isolated_pixels: Boolean
+
+        """
         self.data = None
         self.flt = pysparse.MRFilters(**kwargs)
 
     def filter(self, data):
+        """ Execute the filter operation.
+
+        Parameters
+        ----------
+        data: ndarray
+            the input data.
+        """
         self.data = pysap.Image(data=self.flt.filter(data))
 
     def show(self):  # pragma: no cover
+        """ Show the filtered data.
+        """
         if self.data is None:
             raise AttributeError("The data must be filtered first !")
         self.data.show()
+
+
+class Deconvolve():
+    """ Define the structure that will be used to
+        store the deconvolution result.
+    """
+    def __init__(self, **kwargs):
+        """ Define the deconvolution.
+
+        Parameters
+        ----------
+
+        type_of_deconvolution: int
+        type_of_multiresolution_transform: int
+        type_of_filters: int
+        number_of_undecimated_scales: int
+        sigma_noise: float
+        type_of_noise: int
+        number_of_scales: int
+        nsigma: float
+        number_of_iterations: int
+        epsilon: float
+        psf_max_shift: bool
+        verbose: bool
+        optimization: bool
+        fwhm_param: float
+        convergence_param: float
+        regul_param: float
+        first_guess: string
+        icf_filename: string
+        rms_map: string
+        kill_last_scale: bool
+        positive_constraint: bool
+        keep_positiv_sup: bool
+        sup_isol: bool
+        pas_codeur: float
+        sigma_gauss: float
+        mean_gauss: float
+        """
+        self.data = None
+        self.deconv = pysparse.MRDeconvolve(**kwargs)
+
+    def deconvolve(self, img, psf):
+        """ Execute the filter operation.
+
+        Parameters
+        ----------
+        img: ndarray
+            the input image.
+        psf: ndarray
+            the input psf
+        """
+        self.data = pysap.Image(data=self.deconv.deconvolve(img, psf))
+
+    def show(self):  # pragma: no cover
+        """ Show the deconvolved data.
+        """
+        if self.data is None:
+            raise AttributeError("The data must be deconvolved first !")
+        self.data.show()

pysap/extensions/tools.py

@@ -104,20 +104,30 @@ def mr_deconv(
     """ Wrap the Sparse2d 'mr_deconv'.
     """
     # Generate the command
-    cmd = [
-        "mr_deconv",
-        "-d", type_of_deconvolution,
-        "-t", type_of_multiresolution_transform,
-        "-T", type_of_filters,
-        "-m", type_of_noise,
-        "-n", number_of_scales,
-        "-s", nsigma,
-        "-i", number_of_iterations,
-        "-e", epsilon,
-        "-G", regul_param]
+    cmd = ["mr_deconv"]
+
+    if (type_of_deconvolution != 3):
+        cmd += ["-d", type_of_deconvolution]
+    if (type_of_multiresolution_transform != 2):
+        cmd += ["-t", type_of_multiresolution_transform]
+    if (type_of_filters != 1):
+        cmd += ["-T", type_of_filters]
+    if (type_of_noise != 1):
+        cmd += ["-m", type_of_noise]
+    if (number_of_scales != 4):
+        cmd += ["-n", number_of_scales]
+    if (nsigma != 3):
+        cmd += ["-s", nsigma]
+    if (number_of_iterations != 500):
+        cmd += ["-i", number_of_iterations]
+    if (epsilon != 0.001):
+        cmd += ["-e", epsilon]
+    if (regul_param != 0):
+        cmd += ["-G", regul_param]
+
     for key, value in [
             ("-P", suppress_positive_constraint),
-            ("-S", do_not_auto_shift_max_psf),
+            ("-S", no_auto_shift_max_psf),
             ("-p", detect_only_positive_structure),
             ("-k", suppress_isolated_pixels),
             ("-K", suppress_last_scale),
@@ -136,6 +146,7 @@ def mr_deconv(
             ("-O", optimization)]:
         if value is not None:
             cmd += [key, value]
+
     cmd += [in_image, in_psf, out_image]
 
     # Execute the command

pysap/test/test_binding.py

@@ -36,6 +36,10 @@ def setUp(self):
             get_sample_data(dataset_name="mri-slice-nifti"),
             get_sample_data(dataset_name="astro-ngc2997")
             ]
+        self.deconv_images = [
+            get_sample_data(dataset_name="astro-galaxy"),
+            get_sample_data(dataset_name="astro-psf")
+        ]
         print("[info] Image loaded for test: {0}.".format(
             [i.data.shape for i in self.images]))
         transforms_struct = pysap.wavelist(["isap-2d", "isap-3d"])
@@ -151,13 +155,13 @@ def test_accessors(self):
             band_array[:, :] = 10
             self.assertTrue(numpy.allclose(transform[0, 0], band_array))
 
-    def test_init_filter(self):
+    def test_filter_init(self):
         flt = sp.Filter()
         data = numpy.copy(self.images[0])
         flt.filter(data)
         assert(flt.data is not None)
 
-    def test_default_filter(self):
+    def test_filter_default(self):
         # filter with binding
         flt = sp.Filter()
         data = numpy.copy(self.images[1])
@@ -249,7 +253,7 @@ def test_filter_options_f3_n5(self):
             diff = flt.data - image
             self.assertFalse(diff.all())
 
-    def test_noise_value_error(self):
+    def test_filter_noise_valueerror(self):
         data = numpy.copy(self.images[1])
         with assert_raises(ValueError):
             flt = sp.Filter(epsilon_poisson=5, type_of_noise=2)
@@ -264,6 +268,156 @@ def test_noise_value_error(self):
             flt = sp.Filter(rms_map=in_image, type_of_noise=6)
             flt.filter(data)
 
+    def test_deconv_init(self):
+        deconv = sp.Deconvolve()
+        data = self.deconv_images[0].data
+        psf = self.deconv_images[1].data
+        deconv.deconvolve(data, psf)
+        assert(deconv.data is not None)
+
+    def test_deconv_default(self):
+        deconv = sp.Deconvolve()
+        data = self.deconv_images[0].data
+        psf = self.deconv_images[1].data
+        deconv.deconvolve(data, psf)
+        image = 0
+        # deconvolve with wrapper
+        with pysap.TempDir(isap=True) as tmpdir:
+            in_image = os.path.join(tmpdir, "in.fits")
+            in_psf = os.path.join(tmpdir, "in_psf.fits")
+            out_file = os.path.join(tmpdir, "out.fits")
+            pysap.io.save(data, in_image)
+            pysap.io.save(psf, in_psf)
+
+            pysap.extensions.mr_deconv(in_image, in_psf, out_file)
+            image = numpy.copy(pysap.io.load(out_file))
+            diff = deconv.data.data - image
+            self.assertFalse(diff.all())
+
+    def test_deconv_poisson(self):
+        deconv = sp.Deconvolve(type_of_multiresolution_transform=14,
+                               type_of_noise=2)
+        data = self.deconv_images[0].data
+        psf = self.deconv_images[1].data
+        deconv.deconvolve(data, psf)
+        image = 0
+        # deconvolve with wrapper
+        with pysap.TempDir(isap=True) as tmpdir:
+            in_image = os.path.join(tmpdir, "in.fits")
+            in_psf = os.path.join(tmpdir, "in_psf.fits")
+            out_file = os.path.join(tmpdir, "out.fits")
+            pysap.io.save(data, in_image)
+            pysap.io.save(psf, in_psf)
+
+            pysap.extensions.mr_deconv(in_image, in_psf, out_file,
+                                       type_of_multiresolution_transform=14,
+                                       type_of_noise=2)
+            image = numpy.copy(pysap.io.load(out_file))
+            diff = deconv.data.data - image
+            self.assertFalse(diff.all())
+
+    def test_deconv_d1_t20(self):
+        deconv = sp.Deconvolve(type_of_deconvolution=1,
+                               type_of_multiresolution_transform=20)
+        data = self.deconv_images[0].data
+        psf = self.deconv_images[1].data
+        deconv.deconvolve(data, psf)
+        image = 0
+        # deconvolve with wrapper
+        with pysap.TempDir(isap=True) as tmpdir:
+            in_image = os.path.join(tmpdir, "in.fits")
+            in_psf = os.path.join(tmpdir, "in_psf.fits")
+            out_file = os.path.join(tmpdir, "out.fits")
+            pysap.io.save(data, in_image)
+            pysap.io.save(psf, in_psf)
+
+            pysap.extensions.mr_deconv(in_image, in_psf, out_file,
+                                       type_of_deconvolution=1,
+                                       type_of_multiresolution_transform=20)
+            image = numpy.copy(pysap.io.load(out_file))
+            diff = deconv.data.data - image
+            self.assertFalse(diff.all())
+
+    def test_deconv_u3_d5_n2_p_g5(self):
+        deconv = sp.Deconvolve(number_of_undecimated_scales=3,
+                               type_of_deconvolution=5,
+                               number_of_scales=2,
+                               positive_constraint=False,
+                               sigma_noise=5.)
+        data = self.deconv_images[0].data
+        psf = self.deconv_images[1].data
+        deconv.deconvolve(data, psf)
+        image = 0
+        # deconvolve with wrapper
+        with pysap.TempDir(isap=True) as tmpdir:
+            in_image = os.path.join(tmpdir, "in.fits")
+            in_psf = os.path.join(tmpdir, "in_psf.fits")
+            out_file = os.path.join(tmpdir, "out.fits")
+            pysap.io.save(data, in_image)
+            pysap.io.save(psf, in_psf)
+
+            pysap.extensions.mr_deconv(in_image, in_psf, out_file,
+                                       number_of_undecimated_scales=3,
+                                       type_of_deconvolution=5,
+                                       number_of_scales=2,
+                                       suppress_positive_constraint=False,
+                                       sigma=5.)
+            image = numpy.copy(pysap.io.load(out_file))
+            diff = deconv.data.data - image
+            self.assertFalse(diff.all())
+
+    def test_deconv_d2_m5_S_K(self):
+        deconv = sp.Deconvolve(type_of_noise=5,
+                               psf_max_shift=False,
+                               kill_last_scale=True,
+                               type_of_deconvolution=2)
+        data = self.deconv_images[0].data
+        psf = self.deconv_images[1].data
+        deconv.deconvolve(data, psf)
+        image = 0
+        # deconvolve with wrapper
+        with pysap.TempDir(isap=True) as tmpdir:
+            in_image = os.path.join(tmpdir, "in.fits")
+            in_psf = os.path.join(tmpdir, "in_psf.fits")
+            out_file = os.path.join(tmpdir, "out.fits")
+            pysap.io.save(data, in_image)
+            pysap.io.save(psf, in_psf)
+
+            pysap.extensions.mr_deconv(in_image, in_psf, out_file,
+                                       type_of_noise=5,
+                                       no_auto_shift_max_psf=False,
+                                       suppress_last_scale=True,
+                                       type_of_deconvolution=2)
+            image = numpy.copy(pysap.io.load(out_file))
+            diff = deconv.data.data - image
+            self.assertFalse(diff.all())
+
+    def test_deconv_p_G5_P_f2(self):
+        deconv = sp.Deconvolve(regul_param=5,
+                               positive_constraint=False,
+                               keep_positiv_sup=True,
+                               fwhm_param=2)
+        data = self.deconv_images[0].data
+        psf = self.deconv_images[1].data
+        deconv.deconvolve(data, psf)
+        image = 0
+        # deconvolve with wrapper
+        with pysap.TempDir(isap=True) as tmpdir:
+            in_image = os.path.join(tmpdir, "in.fits")
+            in_psf = os.path.join(tmpdir, "in_psf.fits")
+            out_file = os.path.join(tmpdir, "out.fits")
+            pysap.io.save(data, in_image)
+            pysap.io.save(psf, in_psf)
+
+            pysap.extensions.mr_deconv(in_image, in_psf, out_file,
+                                       regul_param=5,
+                                       suppress_positive_constraint=False,
+                                       detect_only_positive_structure=True,
+                                       icf_fwhm=2)
+            image = numpy.copy(pysap.io.load(out_file))
+            diff = deconv.data.data - image
+            self.assertFalse(diff.all())
+
 
 if __name__ == "__main__":
     unittest.main()