<jupyter_start><jupyter_text>2. Numpy with images All images are essentially matrices with a variable number of dimensions where each element represents the value of one pixel. The different dimensions and the pixel values can have very different meanings depending on the type of image considered, but the structure is the same.Python does not allow by default to gracefully handle multi-dimensional data. In particular it is not desgined to handle matrix operations. Numpy was developed to fill in this blank and offers a very similar framework as the one offered by Matlab. It is underlying a large number of packages and has become abolsutely essential to Python scientific programming. In particular it underlies the functions of scikit-image. The latter in turn forms the basis of other software like CellProfiler. It is thus essential to have a good understanding of Numpy to  proceed.Instead of introducing Numpy in an abstract way, we are going here to present it through the lense of image processing in order to focus on the most useful features in the context of this course.  2.1 Exploring an image Some test images are provided directly in skimage, so let us look at one (we'll deal with the details of image import later). First let us import the necessary packages.<jupyter_code>import numpy as np
import skimage
import matplotlib.pyplot as plt
image = skimage.data.coins()<jupyter_output><empty_output><jupyter_text>The first thing we can do is check the type of the image:<jupyter_code>type(image)<jupyter_output><empty_output><jupyter_text>We thus see that the imported image is indeed a Numpy array. We can also find out what is the type of each pixel:<jupyter_code>image.dtype<jupyter_output><empty_output><jupyter_text>We see here that the image is an 8bit image. Classical format we are going to see are 8bit (uint8), 16bit (uint16) and non-integers (usually float64). The type of the image pixels set what values they can take. For example 8bit means values from $0$ to $2^8 -1= 256-1 = 255$. Since typing is dynamic in Python, if we for example divide our uint8 image by 2 it ends up being a float.<jupyter_code>image
image2 = image/2
image2
image2.dtype<jupyter_output><empty_output><jupyter_text>Let us check what are the image dimensions using the shape method:<jupyter_code>image.shape<jupyter_output><empty_output><jupyter_text>We see that the image has 303 rows and 384 columns. Let us have a look at it:<jupyter_code>plt.imshow(image, cmap = 'gray')
plt.show()<jupyter_output><empty_output><jupyter_text>Let us also look at the actual matrix behind that image:<jupyter_code>image<jupyter_output><empty_output><jupyter_text>Jupyter only shows us a part of this large matrix, adding ... between the first and last few rows/columns. Each value here corresponds to a single pixel of the image. The type of data contained in the array, here 8bit integer, is also displayed.  2.2 Operations on arrays  2.2.2 Arithmetic operations Two types of simple arithmetic operations can be done on arrays. One can modify them using a scalar, for example by adding a constant:<jupyter_code>image
image+10
print(image+10)<jupyter_output>[[ 57 133 143 ...  24  13  22]
 [103 154 155 ...  22  17  17]
 [136 157 153 ...  12  23  13]
 ...
 [ 91  89  84 ...  16  14  17]
 [ 98  92  84 ...  15  17  18]
 [101  89  78 ...  14  20  17]]<jupyter_text>One has to be very careful with the type of the image. For example let's add a large number:<jupyter_code>print(image+230)<jupyter_output>[[ 21  97 107 ... 244 233 242]
 [ 67 118 119 ... 242 237 237]
 [100 121 117 ... 232 243 233]
 ...
 [ 55  53  48 ... 236 234 237]
 [ 62  56  48 ... 235 237 238]
 [ 65  53  42 ... 234 240 237]]<jupyter_text>Clearly something went wrong as we get values that are smaller than 230. As the array has 8bit type, with vlaues ranging from 0 to 255, any value above 255 limit gets reassigned at the start of the range.It is important to note that when doing an operation, Python adopts the most complex type of the variables involved. Hence if we add a float instead of an integer:<jupyter_code>print(image+230.0)<jupyter_output>[[277. 353. 363. ... 244. 233. 242.]
 [323. 374. 375. ... 242. 237. 237.]
 [356. 377. 373. ... 232. 243. 233.]
 ...
 [311. 309. 304. ... 236. 234. 237.]
 [318. 312. 304. ... 235. 237. 238.]
 [321. 309. 298. ... 234. 240. 237.]]<jupyter_text>the result is now correct and is of type float. If needed, arrays can be converted to approriate formats using the astype() method:<jupyter_code>image
image.astype(float)
print(image.astype(float)+230)<jupyter_output>[[277. 353. 363. ... 244. 233. 242.]
 [323. 374. 375. ... 242. 237. 237.]
 [356. 377. 373. ... 232. 243. 233.]
 ...
 [311. 309. 304. ... 236. 234. 237.]
 [318. 312. 304. ... 235. 237. 238.]
 [321. 309. 298. ... 234. 240. 237.]]<jupyter_text>The second type of operation is element-wise by using two arrays of the same dimension. Let's change the image type to avoid the type of issue mentionded above.<jupyter_code>image = image.astype(float)
image.dtype<jupyter_output><empty_output><jupyter_text>Now we can for example multiply the image by itself:<jupyter_code>print(image*image)
image.shape
image3 = image*image
image3.shape<jupyter_output><empty_output><jupyter_text>2.2.3 Logical operations A set of important operations when processing images are logical (or boolean) operations. Those have a very simple syntax in Numpy. For example, let's compare pixel intensities to some value *a*:<jupyter_code>a = 100
image4 = image > a
image4.shape
image4.dtype
image4
image4.astype(int)
plt.imshow(image4)
plt.show()
aa = 100
print(image)
print(image>a)<jupyter_output>[[ 47 123 133 ...  14   3  12]
 [ 93 144 145 ...  12   7   7]
 [126 147 143 ...   2  13   3]
 ...
 [ 81  79  74 ...   6   4   7]
 [ 88  82  74 ...   5   7   8]
 [ 91  79  68 ...   4  10   7]]
[[False  True  True ... False False False]
 [False  True  True ... False False False]
 [ True  True  True ... False False False]
 ...
 [False False False ... False False False]
 [False False False ... False False False]
 [False False False ... False False False]]<jupyter_text>We see that the returned object is a logical matrix (containg True and False) of the same size as the image. That returned object can be directly assigned to a variable, leading to a very compact formulation:<jupyter_code>image_threshold = image > a

print (image_threshold)<jupyter_output>[[False  True  True ... False False False]
 [False  True  True ... False False False]
 [ True  True  True ... False False False]
 ...
 [False False False ... False False False]
 [False False False ... False False False]
 [False False False ... False False False]]<jupyter_text>Of course other logical operator can be used (, ==, !=) and the resulting boolean matrices combined:<jupyter_code>a = 70
b = 100
image_threshold = (image>a)*(image<b)

print (image_threshold)<jupyter_output>[[False False False ... False False False]
 [ True False False ... False False False]
 [False False False ... False False False]
 ...
 [ True  True  True ... False False False]
 [ True  True  True ... False False False]
 [ True  True False ... False False False]]<jupyter_text>2.2.4 Functional operations Numpy has a large trove of functions to act on arrays. Just like the arithmetic operations described above, those functions are designed to be applied to entire arrays, and can either generate a new array or a single value. For example the mean, maximum and minmal pixel values can be abtained like this:<jupyter_code>help(np.mean)
np.mean()
np.mean(image, axis =0).shape
image.shape
np.min(image)<jupyter_output><empty_output><jupyter_text>Alternatively one can take the square root of each array element by doing:<jupyter_code>np.sqrt(image)<jupyter_output><empty_output><jupyter_text>It is impossible to give an exhaustive presentation of Numpy functions. We will simply describe functions that we use a we go through course material.  2.3 Slicing and indexing Let us now use another image to explore how to select specific parts of an image.<jupyter_code>image = skimage.data.chelsea()
image.shape<jupyter_output><empty_output><jupyter_text>We see that the image has three dimensions, probably it's a stack of three images of size 300x400. Let us try to have a look at this image hoping that dimensions are handled gracefully:<jupyter_code>plt.imshow(image)
plt.show()<jupyter_output><empty_output><jupyter_text>So we have an image of a cat with dimensions 300x400. The image being in natural colors, the three dimensions probably indicate an RGB (red, green, blue) format, and the plotting function just knows what to do in that case.  2.3.1 Array slicing Let us now just look at one of the three planes composing the image. To do that, we are going the select a portion of the image array by slicing it. One can give:- a single index e.g. 0 for the first element- a range e.g. 0:10 for the first 10 elements- take all elements using a semi-column :What portion is selected has to be specified for each dimensions of an array. In our particular case, we want to select all rows, all columns and a single plane of the image:<jupyter_code>image.shape
image[:,:,1].shape
image[:,0,0]
plt.imshow(image[:,:,0],cmap='gray')
plt.title('First plane: Red')
plt.show()<jupyter_output><empty_output><jupyter_text>We see now the red layer of the image. We can do the same for the others by specifying planes 0, 1, and 2:<jupyter_code>fig,ax = plt.subplots(figsize=(10,10))
plt.subplot(1,3,1)
plt.imshow(image[:,:,0],cmap='gray')
plt.title('First plane: Red')
plt.subplot(1,3,2)
plt.imshow(image[:,:,1],cmap='gray')
plt.title('Second plane: Green')
plt.subplot(1,3,3)
plt.imshow(image[:,:,2],cmap='gray')
plt.title('Third plane: Blue')
plt.show()<jupyter_output><empty_output><jupyter_text>Logically intensities are high for the red channel and low for the blue channel as the image has red/brown patterns. We can confirm that by measuring the mean of each plane. To do that we use the same function as above but apply it to a singel sliced plane:<jupyter_code>image4 = image[:,:,0]
np.mean(image4)<jupyter_output><empty_output><jupyter_text>and for all planes using a comprehension list:<jupyter_code>np.mean(image[:,:,0])
a = np.zeros((3,4,2,8,4,5))
a.shape
[np.mean(image[:,:,i]) for i in range(3)]<jupyter_output><empty_output><jupyter_text>To look at some more details let us focus on a smaller portion of the image e.g. one of the cat's eyes. For that we are going to take a slice of the red image and store it in a new variable and display the selection. We consider pixel rows from 80 to 150 and columns from 130 to 210 of the first plane (0).<jupyter_code>image_red = image[80:150,130:210,0]
plt.imshow(image_red,cmap='gray')
plt.show()<jupyter_output><empty_output><jupyter_text>There are different ways to select parts of an array. For example one can select every n'th element by giving a step size. In the case of an image, this subsamples the data:<jupyter_code>image_subsample = image[80:150:3,130:210:3,0]
plt.imshow(image_subsample,cmap='gray')
plt.show()<jupyter_output><empty_output><jupyter_text>Or use negative indices or steps to start from the end of the array. In our case this flips dimensions:<jupyter_code>image_reverse = image[150:80:-1,210:130:-1,0]
plt.imshow(image_reverse,cmap='gray')
plt.show()<jupyter_output><empty_output><jupyter_text>2.3.2 Array indexing In addition to slicing an array, we can also select specific values out of it. There are [many](https://docs.scipy.org/doc/numpy-1.13.0/reference/arrays.indexing.html) different ways to achieve that, but we focus here on two main ones.First one might have a list of pixel positions and one wishes to get the values of those pixels. By passing two lists of the same size containing the rows and columns positions of those pixels, one can recover them:<jupyter_code>image_red[0:5,0:5]
image_red.shape
image_red[row_position,col_position]
row_position = [0,1,2,3]
col_position = [0,1,0,1]

print(image_red[0:5,0:5])

image_red[row_position,col_position]<jupyter_output>[[166 162 169 174 185]
 [183 192 185 183 173]
 [179 178 168 175 176]
 [187 184 187 189 185]
 [195 192 187 181 169]]<jupyter_text>Alternatively, one can pass a logical array of the same dimensions as the original array, and only the True pixels are selected. For example, let us create a logical array by picking values above a threshold:<jupyter_code>threshold_image = image_red>120<jupyter_output><empty_output><jupyter_text>Let's visualize it. Matplotlib handles logical arrays simply as a binary image:<jupyter_code>plt.imshow(threshold_image*image_red,cmap = 'gray')
plt.title('Thresholded logical image')
plt.show()
np.where(threshold_image)<jupyter_output><empty_output><jupyter_text>We can recover the value of all the "white" (True) pixels in the original image by using:<jupyter_code>threshold_image.shape
image_red[threshold_image]
selected_pixels.shape
selected_pixels = image_red[threshold_image]
print(selected_pixels)<jupyter_output>[166 162 169 ... 148 137 132]<jupyter_text>And now ask how many pixels are above threshold and what their average value is.<jupyter_code>len(selected_pixels)
np.mean(selected_pixels)<jupyter_output><empty_output><jupyter_text>We now know that there are 2585 pixels above the threshold and that their mean is 153.6  2.4 Creating arrays Often we are going to create new arrays that we are going to then transform. There are a few common ways of doing that.Creating a zero array of dimensions 3x5:<jupyter_code>np.zeros((3,5))<jupyter_output><empty_output><jupyter_text>Same for an array filled with ones:<jupyter_code>np.ones((3,5))<jupyter_output><empty_output><jupyter_text>Create an array filled with random uniform (0-1) values:<jupyter_code>np.random.rand((3,5))<jupyter_output><empty_output><jupyter_text>Clearly we did something wrong. The argument that was passed was not accepted by the functions. This is a good occasion to look at the help for a more complex function, in order to know how to pass arguments.<jupyter_code>help(np.random.rand)<jupyter_output>Help on built-in function rand:

rand(...) method of mtrand.RandomState instance
    rand(d0, d1, ..., dn)
    
    Random values in a given shape.
    
    Create an array of the given shape and populate it with
    random samples from a uniform distribution
    over ``[0, 1)``.
    
    Parameters
    ----------
    d0, d1, ..., dn : int, optional
        The dimensions of the returned array, should all be positive.
        If no argument is given a single Python float is returned.
    
    Returns
    -------
    out : ndarray, shape ``(d0, d1, ..., dn)``
        Random values.
    
    See Also
    --------
    random
    
    Notes
    -----
    This is a convenience function. If you want an interface that
    takes a shape-tuple as the first argument, refer to
    np.random.random_sample .
    
    Examples
    --------
    >>> np.random.rand(3,2)
    array([[ 0.14022471,  0.96360618],  #random
           [ 0.37601032,  0.25528411],  #random
           [ 0.49313049,  0.94909878]][...]<jupyter_text>So we understand now that we can just pass a list of numbers:<jupyter_code>np.zeros((3,5))
np.random.rand((3,5))<jupyter_output><empty_output><jupyter_text>2.5 Assembling arrays It might happen that we have separate images that need to be assembled into one larger matrix. Numpy provides several functions to achieve that. For example we can concatenate images, which means that we append an image to another. Let's grab each channel of the cat image:<jupyter_code>image_red = image[:,:,0]
image_green = image[:,:,1]
image_blue = image[:,:,2]<jupyter_output><empty_output><jupyter_text>And concateante them:<jupyter_code>concatenated_image = np.concatenate((image_red,image_green))
plt.imshow(concatenated_image,cmap = 'gray')
plt.show()<jupyter_output><empty_output><jupyter_text>We can specify along which axis they should be assembled:<jupyter_code>concatenated_image = np.concatenate((image_red,image_green),axis = 1)
plt.imshow(concatenated_image,cmap = 'gray')
plt.show()<jupyter_output><empty_output><jupyter_text>Alternatively, we might want to create a stack of image from separate images. Let's recreate the cat image by switching the red and blue channel. We assemble the separate color images and specific that they should be stacked along the third axis (axis =2).<jupyter_code>switch_image = np.stack((image_blue,image_green,image_red),axis = 2)<jupyter_output><empty_output><jupyter_text>We verify that the dimensions are the expected ones:<jupyter_code>switch_image.shape
plt.imshow(switch_image)
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09571652027433522, 'fraction_numerical': 0.0527012393213813, 'mean_word_length': 4.464497041420119, 'pattern_counts': {'":': 0, '<': 160, '<?xml version=': 0, '>': 167, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '42', 'dir': 'jupyter-structured-clean-dedup', 'id': '331041', 'license': 'MIT', 'n_tokens_mistral': 5832, 'n_tokens_neox': 5114, 'n_words': 2060, 'path': '02-Numpy_images.ipynb', 'repo_name': 'guiwitz/Python_image_processing', 'size': '16622'}