Automatically remove hot / dead pixels from image in python

Question

I use numpy and scipy to handle multiple images taken with a CCD camera. These images have several hot (and dead) pixels with very large (or small) values. They interfere with other image processing, so they must be removed. Unfortunately, although some pixels are stuck at 0 or 255 and always have the same value in all images, there are some pixels that temporarily get stuck in other values ​​for several minutes (data takes many hours).

I am wondering if there is a way to identify (and remove) hot pixels already implemented in python. If not, I wonder what would be an effective method for this. Hot / dead pixels are relatively easy to identify by comparing them with neighboring pixels. I could see how to write a loop that looks at each pixel, compares its value with its 8 nearest neighbors. Or, it seems better to use some kind of convolution to create a smoother image, and then subtract it from the image containing hot pixels, which makes their identification easier.

I tried this “blur method” in the code below and it works fine, but I doubt it is the fastest. In addition, it gets confused at the edge of the image (possibly because the gaussian_filter function takes a convolution, and the convolution becomes weird near the edge). So, is there a better way to do this?

Code example:

import numpy as np import matplotlib.pyplot as plt import scipy.ndimage plt.figure(figsize=(8,4)) ax1 = plt.subplot(121) ax2 = plt.subplot(122) #make a sample image x = np.linspace(-5,5,200) X,Y = np.meshgrid(x,x) Z = 255*np.cos(np.sqrt(x**2 + Y**2))**2 for i in range(0,11): #Add some hot pixels Z[np.random.randint(low=0,high=199),np.random.randint(low=0,high=199)]= np.random.randint(low=200,high=255) #and dead pixels Z[np.random.randint(low=0,high=199),np.random.randint(low=0,high=199)]= np.random.randint(low=0,high=10) #Then plot it ax1.set_title('Raw data with hot pixels') ax1.imshow(Z,interpolation='nearest',origin='lower') #Now we try to find the hot pixels blurred_Z = scipy.ndimage.gaussian_filter(Z, sigma=2) difference = Z - blurred_Z ax2.set_title('Difference with hot pixels identified') ax2.imshow(difference,interpolation='nearest',origin='lower') threshold = 15 hot_pixels = np.nonzero((difference>threshold) | (difference<-threshold)) #Don't include the hot pixels that we found near the edge: count = 0 for y,x in zip(hot_pixels[0],hot_pixels[1]): if (x != 0) and (x != 199) and (y != 0) and (y != 199): ax2.plot(x,y,'ro') count += 1 print 'Detected %i hot/dead pixels out of 20.'%count ax2.set_xlim(0,200); ax2.set_ylim(0,200) plt.show() 

And the conclusion:

Answer 1

Basically, I believe that the fastest way to deal with hot pixels is to simply use an average size filter. Then, poof, your hot pixels disappear, and you also kill all sorts of other high-frequency sensor noises from your camera.

If you really want to remove ONLY hot pixels, then substituting you can subtract the median filter from the original image, as it was in the question, and replace only these values ​​with the values ​​from the middle filtered image. This does not work very well around the edges, so if you can ignore the pixels along the edge, this will greatly simplify the work.

If you want to deal with edges, you can use the code below. However, this is not the fastest:

 import numpy as np import matplotlib.pyplot as plt import scipy.ndimage plt.figure(figsize=(10,5)) ax1 = plt.subplot(121) ax2 = plt.subplot(122) #make some sample data x = np.linspace(-5,5,200) X,Y = np.meshgrid(x,x) Z = 100*np.cos(np.sqrt(x**2 + Y**2))**2 + 50 np.random.seed(1) for i in range(0,11): #Add some hot pixels Z[np.random.randint(low=0,high=199),np.random.randint(low=0,high=199)]= np.random.randint(low=200,high=255) #and dead pixels Z[np.random.randint(low=0,high=199),np.random.randint(low=0,high=199)]= np.random.randint(low=0,high=10) #And some hot pixels in the corners and edges Z[0,0] =255 Z[-1,-1] =255 Z[-1,0] =255 Z[0,-1] =255 Z[0,100] =255 Z[-1,100]=255 Z[100,0] =255 Z[100,-1]=255 #Then plot it ax1.set_title('Raw data with hot pixels') ax1.imshow(Z,interpolation='nearest',origin='lower') def find_outlier_pixels(data,tolerance=3,worry_about_edges=True): #This function finds the hot or dead pixels in a 2D dataset. #tolerance is the number of standard deviations used to cutoff the hot pixels #If you want to ignore the edges and greatly speed up the code, then set #worry_about_edges to False. # #The function returns a list of hot pixels and also an image with with hot pixels removed from scipy.ndimage import median_filter blurred = median_filter(Z, size=2) difference = data - blurred threshold = 10*np.std(difference) #find the hot pixels, but ignore the edges hot_pixels = np.nonzero((np.abs(difference[1:-1,1:-1])>threshold) ) hot_pixels = np.array(hot_pixels) + 1 #because we ignored the first row and first column fixed_image = np.copy(data) #This is the image with the hot pixels removed for y,x in zip(hot_pixels[0],hot_pixels[1]): fixed_image[y,x]=blurred[y,x] if worry_about_edges == True: height,width = np.shape(data) ###Now get the pixels on the edges (but not the corners)### #left and right sides for index in range(1,height-1): #left side: med = np.median(data[index-1:index+2,0:2]) diff = np.abs(data[index,0] - med) if diff>threshold: hot_pixels = np.hstack(( hot_pixels, [[index],[0]] )) fixed_image[index,0] = med #right side: med = np.median(data[index-1:index+2,-2:]) diff = np.abs(data[index,-1] - med) if diff>threshold: hot_pixels = np.hstack(( hot_pixels, [[index],[width-1]] )) fixed_image[index,-1] = med #Then the top and bottom for index in range(1,width-1): #bottom: med = np.median(data[0:2,index-1:index+2]) diff = np.abs(data[0,index] - med) if diff>threshold: hot_pixels = np.hstack(( hot_pixels, [[0],[index]] )) fixed_image[0,index] = med #top: med = np.median(data[-2:,index-1:index+2]) diff = np.abs(data[-1,index] - med) if diff>threshold: hot_pixels = np.hstack(( hot_pixels, [[height-1],[index]] )) fixed_image[-1,index] = med ###Then the corners### #bottom left med = np.median(data[0:2,0:2]) diff = np.abs(data[0,0] - med) if diff>threshold: hot_pixels = np.hstack(( hot_pixels, [[0],[0]] )) fixed_image[0,0] = med #bottom right med = np.median(data[0:2,-2:]) diff = np.abs(data[0,-1] - med) if diff>threshold: hot_pixels = np.hstack(( hot_pixels, [[0],[width-1]] )) fixed_image[0,-1] = med #top left med = np.median(data[-2:,0:2]) diff = np.abs(data[-1,0] - med) if diff>threshold: hot_pixels = np.hstack(( hot_pixels, [[height-1],[0]] )) fixed_image[-1,0] = med #top right med = np.median(data[-2:,-2:]) diff = np.abs(data[-1,-1] - med) if diff>threshold: hot_pixels = np.hstack(( hot_pixels, [[height-1],[width-1]] )) fixed_image[-1,-1] = med return hot_pixels,fixed_image hot_pixels,fixed_image = find_outlier_pixels(Z) for y,x in zip(hot_pixels[0],hot_pixels[1]): ax1.plot(x,y,'ro',mfc='none',mec='r',ms=10) ax1.set_xlim(0,200) ax1.set_ylim(0,200) ax2.set_title('Image with hot pixels removed') ax2.imshow(fixed_image,interpolation='nearest',origin='lower',clim=(0,255)) plt.show() 

Output: