GammaLib

• GFft which performs a FFT on a n-dimensional array of type GNdarray
• GFftWavetable which is a helper class for GFft that contains the trigonometric coefficients for a factorisation

The class so far operates only on 1-dimensional arrays since the GSL does not provide support for more-dimensional arrays. Implementation of such support should however not be too complicated.

The GFft class now also supports 2-dimensional arrays.

Code has been merged into devel.

What remains is the implementation of the operators, and the usage of the GFft class for map smoothing in GSkyMap.

Here a test sequence to illustrate that smoothing of a 2-dimensional image works (this also illustrates how the smoothing kernel needs to be aligned):

Input image	Kernel	Smoothed imaged

And here the code that generated the smoothed image. Note that the code verifies the normalisation of the smoothed image. The sum over all pixels in the kernel needs to be unity.

        # Allocate and set 2-d array
        array = gammalib.GNdarray(10, 10)
        for i in range(3,7):
            array[i,i] = 1.0
        ref = array.sum()

        # Allocate and set 2-d kernel. The kernel needs to be normalized
        # to unity and the centre of the kernel needs to be at pixel [0,0],
        # and it needs to be wraped around to negative indices
        kernel = gammalib.GNdarray(10, 10)
        kernel[0,0] = 0.4
        # 0.4
        kernel[0,1] = 0.1
        kernel[1,0] = 0.1
        kernel[0,9] = 0.1
        kernel[9,0] = 0.1
        # 0.2
        kernel[1,1] = 0.05
        kernel[9,1] = 0.05
        kernel[9,9] = 0.05
        kernel[1,9] = 0.05

        # Smooth 2-d array using FFT
        fft_array  = gammalib.GFft(array)
        fft_kernel = gammalib.GFft(kernel)
        fft_smooth = fft_array * fft_kernel

        # Backtransform
        smooth = fft_smooth.backward()

        # Test sum
        sum = smooth.sum()
        self.test_value(sum, ref)

        # Store in sky map
        map = gammalib.GSkyMap('CAR','CEL',0.0,0.0,-1.0,1.0,10,10)
        for iy in range(10):
            for ix in range(10):
                map[ix+iy*10] = smooth[ix,iy]
        map.save('test_fft.fits', True)

The GSkyMap::smooth() method was added to accomplish the job.

So far the method supports smoothing using a uniform disk kernel and smoothing using a Gaussian kernel.