1294 from images

galsim/chromatic.py

@@ -49,6 +49,7 @@
 from . import utilities
 from . import integ
 from . import dcr
+from .image import Image ##change
 
 
 class ChromaticObject:
@@ -1226,6 +1227,90 @@ def __init__(self, original, waves, oversample_fac=1.0, use_exact_sed=True,
 
         # Don't interpolate an interpolation.  Go back to the original.
         self.deinterpolated = original.deinterpolated
+
+    @classmethod
+    def from_images(self, images, waves, _force_stepk = None, _force_maxk = None, **kwargs):
+        """
+        Alternative initiazliation to InterpolatedChromaticObject from input images at specific wavelenghts. Any parameters accepted by
+        the InterpolatedImage class can be passed in kwargs. Note that stepk and maxk are parameters that can depend on the image size,
+        and therefore on the wavelength. If not given, as a list for every image, or a single number for all images, it will be caluclated
+        using the input image pixel scale and dimensions. This means it will be identical for all images. This will cause small differences
+        from the normal use of this class using chromatic objects whose stepk and maxk are wavelength-dependent. To avoid sanity checks
+        from method initialization, such as wavelength sorting, pixel scale and image dimensions compatibility, simply call the function
+        InterpolatedChromaticObject._from_images directly.
+
+        Parameters:
+            images:          list of Galsim Image objects to use as interpolated images. All images must have the same pixel scale and image
+                             dimensions.
+            wave:            list of wavelength values in nanometers.
+            _force_stepk:    list of step_k values to pass to InterpolatedImages for each image. Can also be single value
+                             to pass for all images alike. If not given stepk is calculated from image pixel scale and dimensions. [Default: None]
+            _force_maxk:     list of max_k values to pass to InterpolatedImages for each image. Can also be single value
+                             to pass for all images alike. If not given maxk is calculated from image pixel scale. [Default: None]
+
+        Returns:
+            An InterpolatedChromaticObject intialized from input images.
+        """ 
+        # check that all images have the same pixel scale and image dimensions. 
+        pix_scale = images[0].scale
+        img_dims = images[0].array.shape
+        for img in images[1:]:
+            if img.scale != pix_scale:
+                raise GalSimNotImplementedError("Cannot interpolate images with different pixel scales.")                         
+            if img.array.shape != img_dims:
+                raise GalSimNotImplementedError("Cannot interpolate images with different image dimensions.")
+        # sort wavelengths and apply sorting to image list
+        sorted_indices = np.argsort(waves)
+        sorted_waves = np.array(waves)[sorted_indices]
+        sorted_images = [images[i] for i in sorted_indices]
+        return self._from_images(sorted_images, sorted_waves, _force_stepk = _force_stepk, _force_maxk = _force_maxk, **kwargs)
+
+    @classmethod
+    def _from_images(cls, images, waves, _force_stepk = None, _force_maxk = None, **kwargs):
+        obj = cls.__new__(cls)  # Does not call __init__ 
+        # use_exact_sed set to false as input images won't have sed available. Ovsersample factor not relevant here, set to 1.0
+        obj.oversample = 1.0
+        obj.use_exact_sed = False
+        obj.separable = False
+        obj.interpolated = True
+        # image properties
+        obj.waves = np.array(waves)
+        pix_scale = images[0].scale 
+        img_dims = images[0].array.shape
+        n_img = len(images)
+
+        stepk = np.zeros(n_img)
+        maxk =  np.zeros(n_img)
+        mean_stepk, mean_maxk = 0.0, 0.0
+        # in the case _force_stepk and/or _force_maxk are passed as lists or single value, set up variables for dummy interpolated object
+        if _force_stepk is not None:
+            mean_stepk = np.mean(_force_stepk)
+        if _force_maxk is not None:
+            mean_maxk = np.mean(_force_maxk)
+
+        # set deinterpolated to a dummy interpolated image. Galsim uses this object to get gsparams and other properties
+        # so setting it to None will cause issues. Only obj.ims and obj.objs will affect future calculations if use_exact_sed = False.
+        # In here we set it to be the average profile in order to calculate the default stepk and maxk to use for all images
+        avg_image = np.zeros(img_dims, dtype=float)
+        for img in images:
+            avg_image += img.array
+        avg_image /= n_img
+        avg_img = Image(avg_image, scale=pix_scale)
+        obj.deinterpolated = InterpolatedImage(avg_img, _force_stepk = mean_stepk, _force_maxk= mean_maxk,  **kwargs)
+        stepk[:] = obj.deinterpolated._stepk
+        maxk[:] = obj.deinterpolated._maxk
+
+        # if stepk and maxk are provided as arguments, overwrite default from average profile
+        if _force_stepk is not None:
+            stepk[:] = _force_stepk
+        if _force_maxk is not None:
+            maxk[:] = _force_maxk
+
+        # set ims and objs manually using the input images
+        obj.ims = images
+        obj.objs = np.array([InterpolatedImage(images[i], _force_stepk=stepk[i], _force_maxk=maxk[i], **kwargs) for i in range(n_img) ])
+        return obj
+
 
     @property
     def sed(self):

tests/test_chromatic.py

@@ -1951,6 +1951,131 @@ def __repr__(self):
     # Also make sure that it ditched the interpolation.
     assert not hasattr(trans_interp_psf, 'waves')
 
+    # test alternate initialization method, "from_images()", of InterpolatedChromaticObject 
+    # that uses images at discrete wavelengths to initialize object.
+
+    # check sorting is done correctly for unsorted wavelengths/images
+    PSF_exact = ChromaticGaussian(sigma_0)
+    PSF = PSF_exact.interpolate(tricky_waves, oversample_fac=oversample_fac)
+    tricky_images = list(PSF.ims[::-1])
+    int_psf = galsim.InterpolatedChromaticObject.from_images(tricky_images, tricky_waves,_force_stepk =PSF.stepk_vals,
+                                                              _force_maxk = PSF.maxk_vals)
+    np.testing.assert_allclose(int_psf.waves, waves, atol=0,
+      err_msg='InterpolatedChromaticObject from_images initialization fails to sort wavelengths')
+
+    for i in range(len(int_psf.ims)):
+        np.testing.assert_allclose(int_psf.ims[i].array, PSF.ims[i].array, atol=1e-17,
+      err_msg='InterpolatedChromaticObject from_images initialization fails to sort images correctly')
+
+    # check error messages from inconsistent pixel scales and image dimensions
+    incorrect_ims = PSF.ims.copy()
+    incorrect_ims[0].scale += 0.01
+    assert_raises(galsim.GalSimNotImplementedError, galsim.InterpolatedChromaticObject.from_images, incorrect_ims, PSF.waves)
+    incorrect_ims[0].scale -= 0.01
+    smaller_image = incorrect_ims[0].array[1:]
+    incorrect_img = galsim.Image(smaller_image, scale=incorrect_ims[0].scale)
+    incorrect_ims[0] = incorrect_img
+    assert_raises(galsim.GalSimNotImplementedError, galsim.InterpolatedChromaticObject.from_images, incorrect_ims, PSF.waves)
+
+    # check input images are correctly intialized from underscored function directly
+    PSF = PSF_exact.interpolate(waves, oversample_fac=oversample_fac)
+    int_psf = galsim.InterpolatedChromaticObject._from_images(PSF.ims, list(PSF.waves),_force_stepk =PSF.stepk_vals,
+                                                              _force_maxk = PSF.maxk_vals)
+    assert isinstance(int_psf.waves, np.ndarray)
+    for i in range(len(int_psf.ims)):
+        np.testing.assert_allclose(int_psf.ims[i].array, PSF.ims[i].array, atol=1e-17,
+      err_msg='InterpolatedChromaticObject from_images initialization fails to initialize correct images')
+    test_obj = galsim.Convolve(int_psf, star)
+    true_obj = galsim.Convolve(PSF, star)
+    true_img = true_obj.drawImage(bandpass, scale=scale, method = 'auto')
+    im_interp = true_img.copy()
+    test_img = test_obj.drawImage(bandpass, image=im_interp, scale=scale)
+    # check drawing objects is identical when using same stepk and maxk as images
+    np.testing.assert_allclose(test_img.array, true_img.array, atol=1e-17,
+            err_msg='InterpolatedChromaticObject from_images initialization fails to reproduce default init. images')
+
+    # without specifying the same stepk and maxk for each image, stepk and maxk are caclulated
+    # based on the input image pixel scale and dimensions and have no wavelength dependance. 
+    int_psf = galsim.InterpolatedChromaticObject.from_images(PSF.ims, PSF.waves)
+    test_obj = galsim.Convolve(int_psf, star)
+    test_img = test_obj.drawImage(bandpass, image=im_interp, scale=scale, method = 'auto')
+    # check images are within 0.01% of the total flux, the flux and maximum within 0.1%
+    np.testing.assert_allclose(test_img.array, true_img.array, atol=1e-4*np.sum(true_img.array),
+        err_msg='InterpolatedChromaticObject images differ when using from_images initialization with default stepk, maxk')
+    np.testing.assert_allclose(test_img.array.sum(), true_img.array.sum(), rtol=1e-3,
+        err_msg='InterpolatedChromaticObject images flux differ when using from_images initialization with default stepk, maxk')
+    np.testing.assert_allclose(test_img.array.max(), true_img.array.max(), rtol=1e-3,
+        err_msg='InterpolatedChromaticObject images maximum differ when using from_images initialization with default stepk, maxk')
+    # check moments
+    truth_mom = galsim.hsm.FindAdaptiveMom(true_img)
+    test_mom = galsim.hsm.FindAdaptiveMom(test_img)
+
+    np.testing.assert_allclose(test_mom.moments_amp,
+                               truth_mom.moments_amp,
+                               rtol=1e-3, atol=0,
+                err_msg='InterpolatedChromaticObject from_images init. differs in moments amplitude ')
+    np.testing.assert_allclose(test_mom.moments_centroid.x,
+                               truth_mom.moments_centroid.x,
+                               rtol=0., atol=1e-2,
+                err_msg='InterpolatedChromaticObject from_images init. differs in moments centroid.x ')
+    np.testing.assert_allclose(test_mom.moments_centroid.y,
+                               truth_mom.moments_centroid.y,
+                               rtol=0., atol=1e-2,
+                err_msg='InterpolatedChromaticObject from_images init. differs in moments centroid.y ')
+    np.testing.assert_allclose(test_mom.moments_sigma,
+                               truth_mom.moments_sigma,
+                               rtol=1e-3, atol=0,
+                err_msg='InterpolatedChromaticObject from_images init. differs in moments sigma ')
+    np.testing.assert_allclose(test_mom.observed_shape.g1,
+                               truth_mom.observed_shape.g1,
+                               rtol=0, atol=1e-4,
+                err_msg='InterpolatedChromaticObject from_images init. differs in moments g1 ')
+    np.testing.assert_allclose(test_mom.observed_shape.g2,
+                               truth_mom.observed_shape.g2,
+                               rtol=0, atol=1e-4,
+                err_msg='InterpolatedChromaticObject from_images init. differs in moments g2 ')
+
+    # check that flux macthes when convolved with non-trivial light profile
+    gal = galsim.Exponential(half_light_radius = 2.*scale)
+    gal = gal.shear(g2 = 0.3)
+    gal = disk_SED*gal
+    obj_exact = galsim.Convolve(PSF, gal)
+    obj_interp = galsim.Convolve(int_psf, gal)
+    true_img = obj_exact.drawImage(bandpass_g, scale=scale)
+    im_interp = true_img.copy()
+    test_img = obj_interp.drawImage(bandpass_g, image=im_interp, scale=scale)
+    np.testing.assert_allclose(
+        test_img.array, true_img.array, atol=1.e-4,
+        err_msg='Interpolated ChromaticObject from_images initialization differs from default when drawing galaxy')
+    # check moments
+    truth_mom = galsim.hsm.FindAdaptiveMom(true_img)
+    test_mom = galsim.hsm.FindAdaptiveMom(test_img)
+
+    np.testing.assert_allclose(test_mom.moments_amp,
+                               truth_mom.moments_amp,
+                               rtol=1e-3, atol=0,
+                err_msg='InterpolatedChromaticObject from_images init. differs in moments amplitude ')
+    np.testing.assert_allclose(test_mom.moments_centroid.x,
+                               truth_mom.moments_centroid.x,
+                               rtol=0., atol=1e-2,
+                err_msg='InterpolatedChromaticObject from_images init. differs in moments centroid.x ')
+    np.testing.assert_allclose(test_mom.moments_centroid.y,
+                               truth_mom.moments_centroid.y,
+                               rtol=0., atol=1e-2,
+                err_msg='InterpolatedChromaticObject from_images init. differs in moments centroid.y ')
+    np.testing.assert_allclose(test_mom.moments_sigma,
+                               truth_mom.moments_sigma,
+                               rtol=1e-3, atol=0,
+                err_msg='InterpolatedChromaticObject from_images init. differs in moments sigma ')
+    np.testing.assert_allclose(test_mom.observed_shape.g1,
+                               truth_mom.observed_shape.g1,
+                               rtol=0, atol=1e-4,
+                err_msg='InterpolatedChromaticObject from_images init. differs in moments g1 ')
+    np.testing.assert_allclose(test_mom.observed_shape.g2,
+                               truth_mom.observed_shape.g2,
+                               rtol=0, atol=1e-4,
+                err_msg='InterpolatedChromaticObject from_images init. differs in moments g2 ')
+
 
 @timer
 def test_ChromaticOpticalPSF():