Histogram Equalization Ignoring zero in Python I wrote a small python code for histogram equalization (ignoring zero values) for an image but it is taking too long to run. Code is written with reference to wikipedia article on histogram equalization. 

def hist_eq_updated(img):
    img_flatten = img.flatten()
    img_flat_nnzero = img_flatten[np.where(img_flatten!=0)]
    image_histogram, bins = np.histogram(img_flat_nnzero,256) 
    #print(image_histogram,bins)
    cdf = image_histogram.cumsum() # cumulative distribution function
    #print(cdf)

    eq_img_flatten = np.zeros_like(img_flatten)
    cdf_min = min(cdf)
    cdf_max = max(cdf)
    for i in range(len(img_flatten)):
         val = img_flatten[i]
         ind = 0
         flag=False
         for j in range(len(bins)):
             if val==0:
                 flag=True
             if (val&gt;=bins[j]) and (val&lt;=bins[j+1]) or flag==True:
                 break;
             ind = ind + 1
         if flag==True:
             eq_img_flatten[i] = 0
         else:
             eq_img_flatten[i] = (((cdf[ind] - cdf_min)/(cdf_max-cdf_min))*254) + 1
    return eq_img_flatten.reshape(img.shape)


I am new to image processing in python and not familiar with any other "pythonic" way to do same. Thanks in advance !!!  

EDIT - Sorry, I forgot to mention that image size is ~4000 * ~4000, image type is float32.
Actually, this is not a signal processing quesiton, and would have been better at stackoverflow. Anyway, you can use a look-up-table for this problem. The output value of a pixel does only depend on the pixel value, not its position in the image. So, it makes no sense to go through each pixel. Instead, just go through each possible pixel value. Then, later on transform your original pixels to the new pixels using this look-up table

import skimage
def hist_eq_updated(img):
    img_flatten = img.flatten()
    img_flat_nnzero = img_flatten[np.where(img_flatten!=0)]
    image_histogram, bins = np.histogram(img_flat_nnzero,256) 
    #print(image_histogram,bins)
    cdf = image_histogram.cumsum() # cumulative distribution function
    #print(cdf)

    eq_img_flatten = np.zeros_like(img_flatten)
    cdf_min = min(cdf)
    cdf_max = max(cdf)
    LUT = np.zeros(256, dtype=img.dtype)
    for val in range(256):
         ind = 0
         flag=False
         for j in range(len(bins)):
             if val==0:
                 flag=True
             if (val&gt;=bins[j]) and (val&lt;=bins[j+1]) or flag==True:
                 break;
             ind = ind + 1
         if flag==True:
             LUT[val] = 0
         else:
             LUT[val] = (((cdf[ind] - cdf_min)/(cdf_max-cdf_min))*254) + 1
    return LUT[img]

Img = skimage.io.imread("http://upload.wikimedia.org/wikipedia/commons/thumb/4/4a/Portrait_elisha_gray.jpg/220px-Portrait_elisha_gray.jpg")

plt.figure()
plt.subplot(121)
plt.imshow(Img, interpolation='none', cmap='gray')
plt.subplot(122)
plt.imshow(hist_eq_updated(Img), interpolation='none', cmap='gray')