Locating the solar disk
Uses
-
-
Step
returns
We must restrict the search for sunspots to the solar disk itself, to avoid picking up noise from the background. The intensity analysis shows that an intensity of 6000 is a good threshold to separate the solar disk from the dark current of the CCD chip.
solar_disk_threshold = 750
As a precaution, let's check that this threshold corresponds to a minimum in the intensity histogram.
counts = np.bincount(raw_image.flatten())
max_count = counts[solar_disk_threshold-100:solar_disk_threshold+100].max()
if max_count > 100:
print(f'Warning: many pixels are near the solar disk intensity threshold ({max_count})', file=sys.stderr)
We take all pixels above the threshold to be inside the solar disk. We can then compute its center as the centroid of the surface defined by these pixels, and the radius from the radius of gyration of the disk.
temp_solar_disk_mask = (raw_image >= solar_disk_threshold).astype(np.uint8) temp_solar_disk_pixel_count = np.sum(temp_solar_disk_mask) x_indices = np.arange(raw_image.shape[1])[np.newaxis, :] y_indices = np.arange(raw_image.shape[0])[:, np.newaxis] x_center = np.sum(temp_solar_disk_mask * x_indices) / temp_solar_disk_pixel_count y_center = np.sum(temp_solar_disk_mask * y_indices) / temp_solar_disk_pixel_count radial_distance_sq = temp_solar_disk_mask * ((x_indices-x_center)**2 + (y_indices-y_center)**2) radius = np.sqrt(2.*np.sum(radial_distance_sq) / temp_solar_disk_pixel_count) (x_center,y_center, radius)
Some functions used later require integer parameters for the solar disk.
x_center = int(round(x_center)) y_center = int(round(y_center)) radius = int(round(radius))
A sanity check: the radius must be smaller than half the image size.
if radius >= raw_image.shape[0] // 2:
print(f'Warning: solar disk radius cannot be {radius}', file=sys.stderr)
Next, we make a new perfectly circular mask from the center and radius, plus a new solar disk image in which all background pixles are set to zero.
import skimage.draw solar_disk_mask = np.zeros(raw_image.shape, dtype=np.uint8) yy, xx = skimage.draw.disk((y_center, x_center), radius, shape=solar_disk_mask.shape) solar_disk_mask[yy, xx] = 1 solar_disk = solar_disk_mask*raw_image