Here's an example image to make the problem very clear:

Crop of Hubble R136 image

I've used OpenCV's blob_log() as well as the rank.otsu threshold approaches. The problem with the blob detection is that the blob radii are quantized, e.g., 3, 7, 11, etc. The otsu local adaptive threshold works fairly well but needs severe hyperparameter tuning -- in other words, it works in some locations but not others.

I don't know how I can ask the following question: "find all circular blobs in the image. If the blob is large (meaning it is overexposed), expand the search to include diffraction spikes (at known fixed angles) plus the Airy disk."

My goal is to remove the stars and their optical effects from the image, then use an inpainting algorithm to fill in the missing pixels. I have several inpainting candidates to try.

There is noise in the image, and I'm fine with missing any stars that are indistinguishable from noise.

This is one channel -- the other image channels do add additional information, and the actual removal of stars/effects will be based on the logical or of all three analyzed channels.

I am aware that there are other methods for star removal that pop up in google search. The best automatic one I've used is starnet, e.g. in https://github.com/nekitmm/starnet . The issue I have with it is it doesn't separate out the star mask but goes ahead and does the inpainting -- unfortunately, very poorly in large regions. And it doesn't handle spikes well, per the author's commentary and my observations.

I wanted to avoid a machine learning approach. Instead, I wanted to solve the problem directly from first principles. It seems feasible to build a structural model of the star/spike/airy disk with some effort(e.g., as central disk at max intensity becomes larger, spikes grow larger and so does Airy disk.) But how do I match this model to the objects in the image?

Thanks in advance.

The desired output is a 1 bit mask image. Bit set to 1 if this is a star pixel, 0 otherwise. I don't have a problem converting a list of (y,x,r) values to a bit mask, and already have code to do that.

The blob detectors I have used (a) don't compute accurate blob sizes (b) generate blobs for some of the noisy areas.

Regarding finding the isolated blobs being a snap, give me an idea of what you consider is a "snap." My very first attempt looked for local maxima and attempted to grow a circular region. It was very slow and had problems with noise.

The region-growing approaches all failed miserably. A different 2d filtering approach turned out to be more promising. The basic idea is to:

  1. Capture all of the overly bright areas via thresholding, e.g., against around 248 in the bottom image. That gives the first mask = image > 248, e.g.
  2. Apply an emboss filter at 45 and 135 degree angles (filter kernels are [-1 -1 0][-1 0 1][0 1 1] plus 128 for the 45 degree angle.)
  3. Threshold values more than 45 away from grey (128). mask = (embossed < 128-45) | (embossed > 128 + 45), e.g.
  4. Or together all of the three thresholds.

That gives the following result:

After simple processing

That approach has two magic numbers, yes. But it works fairly well for such a simple algorithm.

I'm aware the blobs are a bit too big in places. Also, the blue and red channels of the original image don't work quite as well due to more noise in them. Time to experiment further.

  • $\begingroup$ What is your desired output? A list of "blobs" and characteristics? Say (x,y, r1,r2) Where r1 might be your blob radius and r2 the length of the diffraction spikes? Finding the isolate little blobs is a snap, it's your overlap and cloudy areas that will be tougher. $\endgroup$ – Cedron Dawg Aug 19 '20 at 0:25
  • $\begingroup$ If you can put a bounding box around a blob that is black on the edges, you simply calculate the weighted average of the pixel values in regards to your coordinates and the mean will be the center and you can get std dev values that will give you width parameters. Being "a snap" conceptually does not necessarily mean computationally efficient. Is this a commercial situation or personal research? $\endgroup$ – Cedron Dawg Aug 19 '20 at 14:38
  • $\begingroup$ Just personal research. Yes, I had already done the weighted mean to compute the center of a small round blob. I ran into issues when blobs overlapped. So I thought perhaps a structural matching approach would work better. E.g., here are, say, Airy-shaped blobs ranging from 1 pixel in size to 30. What is the algorithm that chooses the correct blob size and works even in overlapping blob cases? I just have the feeling I am not seeing a whole class of obvious algorithms. This matching problem has to have been solved, right? $\endgroup$ – Walt Donovan Aug 19 '20 at 19:27

Here is a first stab. It is essentially a doubly smoothed contour map.

enter image description here

  • $\begingroup$ Thanks, interesting idea. I'm going to try the algorithm in section 7.1 of people.kth.se/~tony/papers/Lindeberg-PhD-thesis-24-May-1991.pdf since that seems to match my intuitive understanding of what a blob is. $\endgroup$ – Walt Donovan Aug 20 '20 at 6:39
  • $\begingroup$ @WaltDonovan That was a big file. From first appearances, what I did is sort of a quick and dirty version. Coincidentally there seemed to be a large overlap with the Gaussian kernel which is my focus of study right now, finding the exact discrete version vs using a sampled continuous version. Note, the levels I choose were somewhat arbitrary. You could use a much finer gradation. At some level, two adjacent humps will separate into separate contour lines (always orthogonal to the gradient). $\endgroup$ – Cedron Dawg Aug 21 '20 at 0:00

Here's the code

import cv2 as cv
import numpy as np

img = cv.imread("r136-crop.png")
b, g, r = cv.split(img)
b = b.astype(np.uint16)     # so we can accumulate in 16 bits

# filter kernels
diag45 = np.array(( [-1, -1, 0], [-1, 0, 1], [0, 1, 1]))
diag135 = np.array(( [0, -1, -1], [1, 0, -1], [1, 1, 0]))

delta = 35      # delta from grey for an emblss signal
overbright = 240  # threshold for "white"

# take mean of three channels
# this also cuts down noise a bit
b = b + g + r + 2
chan = (b//3).astype(np.uint8)

cv.createTrackbar('delta', 'out', delta, 128, nothing)
cv.createTrackbar('overbright', 'out', overbright, 255, nothing)

d45out = np.zeros(chan.shape, np.uint8)
d135out = np.zeros(chan.shape, np.uint8)

cv.filter2D(chan, -1, diag45, d45out, delta=128, borderType=cv.BORDER_REPLICATE)
cv.filter2D(chan, -1, diag135, d135out, delta=128, borderType=cv.BORDER_REPLICATE)

delta = cv.getTrackbarPos('delta', 'out')
overbright = cv.getTrackbarPos('overbright', 'out')

# get mask of brightest pixels -- apply to unblurred image
mask = (chan >= overbright).astype(np.uint8)

# find extremal values in the embossed images
mask45 = ((d45out < 128 - delta) | (d45out > 128 + delta)).astype(np.uint8)
mask135 = ((d135out < 128 - delta) | (d135out > 128 + delta)).astype(np.uint8)

mask = ((mask | mask45 | mask135) * 255).astype(np.uint8)

# dilate the mask to cover the little spots we missed
element = cv.getStructuringElement(cv.MORPH_RECT, (3,3))
mask = cv.dilate(mask, element)

masked_img = img
# mask out the stars in the image
for ch in range(3):
    masked_img[:,:,ch] = masked_img[:,:,ch] & (255 - mask)

dst = cv.inpaint(masked_img, mask, 6, cv.INPAINT_NS )

cv.imshow('out', dst)
cv.imshow('masked', masked_img)
cv.imshow('mask', mask)

Original Image Original Image

Masked Image Masked Image

Infilled Image Infilled image

The result is somewhat acceptable, but removing the halo around the brightest stars is going to require a very different approach and a rethinking of what it is exactly I want to accomplish here.

It's good to see that a fairly simple algorithm is able to do fairly well what seemed to be a fairly simple task.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.