Search code examples
pythonnumpyindexingmaskscikit-image

Mask an image (np.ndarray) section which lies between two curves

Borrowing this example from astropy:

import numpy as np
from matplotlib import pyplot as plt
from astropy.io import fits
from astropy.wcs import WCS
from astropy.utils.data import download_file

fits_file = 'http://data.astropy.org/tutorials/FITS-images/HorseHead.fits'
image_file = download_file(fits_file, cache=True)
hdu = fits.open(image_file)[0]
wcs = WCS(hdu.header)

fig = plt.figure()
ax = fig.add_subplot(111)
plt.imshow(hdu.data, origin='lower', cmap='cubehelix')
plt.xlabel('X')
plt.ylabel('Y')

x_array = np.arange(0, 1000)
line_1 = 1 * x_array + 20 * np.sin(0.05*x_array)
line_2 = x_array - 100 + 20 * np.sin(0.05*x_array)

plt.plot(x_array, line_1, color='red')
plt.plot(x_array, line_2, color='red')

ax.set_xlim(0, hdu.shape[1])
ax.set_ylim(0, hdu.shape[0])

plt.show()

I would like to calculate the median pixel value (in y direction for example) which lies between the two curves:

enter image description here

I believe the clever thing would be to create a mask for the region of interest.

Is there a way to generate this mask without looping across the image pixels?

Edit 1: Modified question improve understanding

Edit 2: I changed the example to better represent the question title (curves instead of straight lines)

Solution

  • A fairly straight forward way of generating such a mask would be to use the numpy.mgrid thingy that essentially gives you x and y-coordinate arrays (in that 2D case) that can then be used to compute the mask with the equation of the lines .

    Edit : Provided you can express your mask using an equation (like f(x,y)<0 where f is any function you want) or a combination of those equations, you can do anything you want. Here is an example with your new mask along with some extra art pieces: 

    import numpy as np
import matplotlib.pyplot as plt

y,x=np.mgrid[0:1000,0:1000]

#a wiggly ramp
plt.subplot(221)
mask0=(y < x + 20 * np.sin(0.05*x)) & (y > x - 100 + 20 * np.sin(0.05*x))
plt.imshow(mask0,origin='lower',cmap='gray')

#a ramp
plt.subplot(222)
mask1=(y < x) & (x < y+100)
plt.imshow(mask1,origin='lower',cmap='gray')

#a disk
plt.subplot(223)
mask2=(200**2>(x-500)**2+(y-500)**2)
plt.imshow(mask2,origin='lower',cmap='gray')

#a ying-yang attempt
plt.subplot(224)
mask3= (mask2 & (0 < np.sin(3.14*x/250)*100 + 500 - y) & (30**2 < (x-620)**2+(y-500)**2) )| (30**2 > (x-380)**2+(y-500)**2)
plt.imshow(mask3,origin='lower',cmap='gray')

plt.show()

    Output :

    Output