Add calculateSigma method

(GalSim-developers#308)

galsim/base.py

@@ -403,6 +403,56 @@ def calculateHLR(self, size=None, scale=None):
         return hlr
 
 
+    def calculateSigma(self, size=None, scale=None):
+        """Returns the sqrt of the radial second moment of the object, sigma = sqrt(T/2), where
+        T = Irr = Ixx+Iyy.
+
+        If the profile has a sigma attribute (only true for a Gaussian), it will just return that,
+        but in the general case, we draw the profile and estimate T directly.
+
+        The function optionally takes size and scale values to use for the image drawing.
+        The default scale is the nyquist scale, which generally produces results accurate
+        to about 1 decimal place.  Using a smaller scale will be more accurate at the expense
+        of speed.  The default size is None, which means drawImage will choose a size designed
+        to contain around 99.5% of the flux.  Using a larger size will again be more accurate
+        at the expense of speed.
+
+        Note: The results from this calculation should be taken as approximate at best.
+              They should usually be acceptable for things like testing that a galaxy has a
+              reasonable resolution, but they should not be trusted for very fine grain
+              discriminations.  For a more accurate estimate, see galsim.hlm.FindAdaptiveMom.
+
+        @param size         If given, the stamp size to use for the drawn image. [default: None,
+                            which will let drawImage choose the size automatically]
+        @param scale        If given, the pixel scale to use for the drawn image. [default:
+                            self.nyquistScale()]
+
+        @returns an estimate of sigma = sqrt(T/2)
+        """
+        if hasattr(self, 'sigma'): 
+            return self.sigma
+
+        if scale is None:
+            scale = self.nyquistScale()
+
+        # Draw the image.  Note: need a method that integrates over pixels to get flux right.
+        im = self.drawImage(nx=size, ny=size, scale=scale)
+
+        # Use radii at centers of pixels as approximation to the radial integral
+        x,y = np.meshgrid(range(im.array.shape[0]), range(im.array.shape[1]))
+        x = x - (im.trueCenter().x - im.bounds.xmin)
+        y = y - (im.trueCenter().y - im.bounds.ymin)
+        rsq = x*x + y*y
+
+        # Accumulate Irr
+        Irr = np.sum(rsq * im.array) / self.flux
+
+        # This has all been done in pixels.  So normalize according to the pixel scale.
+        sigma = np.sqrt(Irr/2.) * im.scale
+
+        return sigma
+
+
     @property
     def flux(self): return self.getFlux()
     @property

tests/test_calc.py

@@ -101,5 +101,64 @@ def test_hlr():
     print 'time for %s = %.2f'%(funcname(),t2-t1)
 
 
+def test_sigma():
+    """Test the calculateSigma method.
+    """
+    import time
+    t1 = time.time()
+
+    # Compare the calculation for a simple Gaussian.
+    g1 = galsim.Gaussian(sigma=5, flux=1.7)
+
+    print 'g1 native sigma = ',g1.sigma
+    print 'g1.calculateSigma = ',g1.calculateSigma()
+    # These should be exactly equal.
+    np.testing.assert_equal(g1.sigma, g1.calculateSigma(),
+                            err_msg="Gaussian.calculateSigma() returned wrong value.")
+
+    # Check for a convolution of two Gaussians.  Should be equivalent, but now will need to 
+    # do the calculation.
+    g2 = galsim.Convolve(galsim.Gaussian(sigma=3, flux=1.3), galsim.Gaussian(sigma=4, flux=23))
+    test_sigma = g2.calculateSigma()
+    print 'g2.calculateSigma = ',test_sigma
+    print 'ratio - 1 = ',test_sigma/g1.sigma-1
+    np.testing.assert_almost_equal(test_sigma/g1.sigma, 1.0, decimal=1,
+                                   err_msg="Gaussian.calculateSigma() is not accurate.")
+
+    # The default scale and size is only accurate to around 1 dp.# Using scale = 0.1 is accurate
+    # to 3 dp.
+    test_sigma = g2.calculateSigma(scale=0.1)
+    print 'g2.calculateSigma(scale=0.1) = ',test_sigma
+    print 'ratio - 1 = ',test_sigma/g1.sigma-1
+    np.testing.assert_almost_equal(test_sigma/g1.sigma, 1.0, decimal=3,
+                                   err_msg="Gaussian.calculateSigma(scale=0.1) is not accurate.")
+
+    # Next, use an Exponential profile
+    e1 = galsim.Exponential(scale_radius=5, flux=1.7)
+
+    # The true "sigma" for this is analytic, but not an attribute.
+    e1_sigma = np.sqrt(3.0) * e1.scale_radius
+    print 'true e1 sigma = sqrt(3) * e1.scale_radius = ',e1_sigma
+
+    # Test with the default scale and size.
+    test_sigma = e1.calculateSigma()
+    print 'e1.calculateSigma = ',test_sigma
+    print 'ratio - 1 = ',test_sigma/e1_sigma-1
+    np.testing.assert_almost_equal(test_sigma/e1_sigma, 1.0, decimal=1,
+                                   err_msg="Exponential.calculateSigma() is not accurate.")
+
+    # The default scale and size is only accurate to around 1 dp.  This time we have to both
+    # decrease the scale and also increase the size to get 3 dp of precision.
+    test_sigma = e1.calculateSigma(scale=0.1, size=2000)
+    print 'e1.calculateSigma(scale=0.1) = ',test_sigma
+    print 'ratio - 1 = ',test_sigma/e1_sigma-1
+    np.testing.assert_almost_equal(test_sigma/e1_sigma, 1.0, decimal=3,
+                                   err_msg="Exponential.calculateSigma(scale=0.1) is not accurate.")
+
+    t2 = time.time()
+    print 'time for %s = %.2f'%(funcname(),t2-t1)
+
+
 if __name__ == "__main__":
     test_hlr()
+    test_sigma()