如何避免用循环加速Python中FITS图像环形统计计算？

Hey there! I totally get your frustration—looping over 4 million pixels in Python is going to be slow no matter how fast your single function is. The good news is we can replace that entire nested loop with vectorized convolution operations, which are implemented in optimized C code and will speed things up by orders of magnitude.

Let's Break Down the Problem First #

Your snr function calculates two key values for each pixel's ring region:

1. Nmask: The number of pixels in the ring (inside rad2, outside rad1)
2. sum(array[mask]): The total intensity of those ring pixels

Then it computes 0.6052697 * abs(Nmask * flux - sum_ring). Instead of calculating this for each pixel one by one, we can compute these values for the entire array at once using convolution.

Step 1: Create a Reusable Ring Mask Template #

First, we make a fixed-size mask that represents our ring shape. This mask will be used as a convolution kernel to compute sums and counts across the entire image.

import numpy as np
from scipy.signal import convolve2d

def create_ring_kernel(rad1, rad2):
    # Create a square kernel with side length 2*rad2 + 1 (centered at rad2)
    kernel_size = 2 * rad2 + 1
    y, x = np.ogrid[:kernel_size, :kernel_size]
    center = rad2
    # Calculate squared distance from center
    dist_sq = (x - center)**2 + (y - center)**2
    # Mask pixels between rad1 and rad2 (inclusive)
    ring_mask = (dist_sq >= rad1**2) & (dist_sq <= rad2**2)
    return ring_mask.astype(float)  # Convert to float for convolution

Step 2: Compute Ring Sums and Pixel Counts with Convolution #

Convolution lets us slide our ring kernel over every pixel in the image, calculating the sum of overlapping pixels (for sum_ring) and the number of valid ring pixels (for Nmask) in one go.

# Load your FITS data as before
from astropy.io import fits  # Assuming you're using astropy for FITS handling

frame1 = fits.open(in_frame, mode='readonly')
data1 = frame1[ext].data
ny, nx = data1.shape
r1 = 5
r2 = 7

# Create the ring kernel
ring_kernel = create_ring_kernel(r1, r2)

# Compute sum of ring pixels for every position
sum_ring = convolve2d(data1, ring_kernel, mode='same', boundary='fill', fillvalue=0)

# Compute number of valid ring pixels for every position (handles edge cases)
nmask_array = convolve2d(np.ones_like(data1), ring_kernel, mode='same', boundary='fill', fillvalue=0)

Step 3: Calculate the Noise Array Vectorially #

Now we can compute the entire noise1 array in a single line, no loops needed:

# Apply your noise formula to the entire array
noise1 = 0.6052697 * np.abs(nmask_array * data1 - sum_ring)

Why This Is So Much Faster #

• No Python Loop Overhead: The convolution operations run in optimized C code, avoiding the slow Python-level loop that was processing 4 million pixels one by one.
• Reusable Kernel: We only create the ring mask once, instead of recreating it 4 million times in your original function.
• Vectorized Calculation: The final noise calculation uses NumPy's vectorized operations, which are far faster than per-pixel function calls.

Notes on Edge Behavior #

This code matches the behavior of your original function:

• For pixels near the image edge, the convolution automatically ignores (fills with 0) pixels outside the image bounds, which is equivalent to how your original mask only includes valid pixels within the image.
• The mode='same' parameter ensures the output noise1 array has the same dimensions as your input data1.

Performance Comparison #

On a typical machine, this approach should run in seconds instead of the minutes your original loop would take—we're talking about a 100-1000x speedup depending on your hardware.

内容的提问来源于stack exchange，提问作者M Neeser