#!/usr/bin/env python
# coding: utf-8

# # 2. Numpy with images

# All images can be seen as matrices where each element represents the value of one pixel. Multiple images (or matrices) can be stacked together into a "multi-dimensional matrix" or tensor. The different dimensions and the pixel values can have very different meanings depending on the type of image considered (e.g. time-lapse, volume, etc.), but the structure remains the same.
# 
# Python does not allow by default to gracefully handle such multi-dimensional data. In particular it is not designed to handle matrix operations. [Numpy](https://numpy.org/) was developed to fill this blank and offers a clean way to handle creation and operations on such data called here **```arrays```**. Numpy is underlying a large number of packages and has become absolutely essential to Python scientific programming. In particular it underlies the functions of [scikit-image](https://scikit-image.org/), the main package used in this course for image processing. Scikit-image is then itself an important basis of other "higher-level" software like [CellProfiler](https://cellprofiler.org/). So Numpy really sits at the basis of the Python scientific world, and it is thus essential to have an understanding of it.
# 
# Instead of introducing Numpy in an abstract way, we are going here to present it through the lens of image processing in order to focus on the most useful features in the context of this course. For the moment, let's just import Numpy:

# In[1]:


import numpy as np


# We also import scikit-image and matplotlib. You will learn more on these packages in later chapters. At the moment we just want to be able to load and display a an image (i.e. a matrix) to understand Numpy.

# In[2]:


import skimage.io
import skimage
import matplotlib.pyplot as plt
plt.gray();


# ## 2.1 Exploring an image

# We start by importing an image using scikit-image (see the next notebook [03-Image_import](03-Image_import.ipynb) for details). We have some data already available in the Data folder:

# In[3]:


image = skimage.io.imread('../Data/neuron.tif')


# ### 2.1.1 Image size

# The first thing we can do with the image is simply look at the output:

# In[4]:


image


# We see that Numpy tells us we have an array and we don't have a simple list of pixels, but a sort of *list of lists* representing the fact that we are dealing with a two-dimensional object. Each list represents one row of pixels. Jupyter smartly only shows us the first/last rows/columns. As mentioned before, the array is the fundamental data container implemented by Numpy, and **ALL** the images we are going to process in this course will be Numpy arrays.
# 
# We have seen before that all Python variables are also objects with methods and attributes. Numpy arrays are no exception. For example we can get information such as the ```shape``` of the array, i.e. number of rows, columns, planes etc.

# In[5]:


image.shape


# This means that we have an image of 1024 rows and 1280 columns. We can now look at it using matplotlib (see next chapter for details) using the ```plt.imshow()``` function:

# In[6]:


plt.imshow(image);


# ### 2.1.2 Image type

# In[7]:


image


# In the output above we see that we have one additional piece of information: the array has ```dtype = uint8``` , which means that the image is of type *unsigned integer 8 bit*. We can also get the type of an array by using:

# In[8]:


image.dtype


# Standard formats we are going to see are 8bit (uint8), 16bit (uint16) and non-integers (usually float64). The image type sets what values pixels can take. For example 8bit means values from $0$ to $2^8 -1= 256-1 = 255$. Just like for example in Fiji, one can change the type of the image. If we know we are going to do operations requiring non-integers we can turn the pixels into floats trough the ```.astype()``` method.

# In[9]:


image_float = image.astype(np.float64)


# In[10]:


image_float.dtype


# The importance of the image type goes slightly against Python's philosophy of dynamics typing (no need to specify a type when creating a variable), but is a necessity when handling scientific images. We are going to see now what sort of operations we can do with arrays, and the importance of *types* is going to be more obvious.

# ## 2.2 Operations on arrays

# ### 2.2.1 Arithmetics on arrays

# Numpy is written in a smart way such that it is able to handle operations between arrays of different sizes. In the simplest case, one can combine a scalar and an array, for example through an addition:

# In[11]:


image+10


# Here Numpy automatically added the scalar 10 to **each** element of the array. Beyond the scalar case, operations between arrays of different sizes are also possible through a mechanism called broadcasting. This is an advanced (and sometimes confusing) features that we won't use in this course but about which you can read for example [here](https://jakevdp.github.io/PythonDataScienceHandbook/02.05-computation-on-arrays-broadcasting.html).
# 
# The only case we are going to consider here are operations between arrays of same size. For example we can multiply the image by itself. We use first the float version of the image:

# In[12]:


image_sq = image_float*image_float


# In[13]:


image_sq


# In[14]:


image_float


# Looking at the first row we see $32^2 = 1024$ and $31^2=961$ etc. which means that the multiplication operation has happened **pixel-wise**. Note that this is **NOT** a classical matrix multiplication. We can also see that the output has the same size as the original arrays:

# In[15]:


image_sq.shape


# In[16]:


image_float.shape


# Let's see now what happens when we square the original 8bit image:

# In[17]:


image*image


# We see that we don't get at all the expected result. Since we multiplied two 8bit images, Numpy assumes we want an 8bit output. And therefore the values are bound between 0-255. For example the second value of the first row is just the remainder of the modulo 256:

# In[18]:


961%256


# The same thing happens e.g. if we add an integer scalar to the matrix. We can also visualize the result:

# In[19]:


print(image+230)


# In[20]:


plt.imshow(image+230)


# Clearly something went wrong as we get values that are smaller than 230. The image also has a strange aspect, with dark regions suddenly brighter. Again any value **over-flowing** above 255 restarts at 0.
# 
# This problem can be alleviated in different ways. For example we can combine a integer array with a float scalar and Numpy will automatically give a result using the "most complex" type:

# In[21]:


image_plus_float = image+230.0


# In[22]:


print(image_plus_float)


# To be on the safe side we can also explicitely change the type when we know we might run into this kind of trouble. This can be done via the ```.astype()``` method: 

# In[23]:


image_float = image.astype(float)


# In[24]:


image_float.dtype


# Again, if we combine floats and integers the output is going to be a float: 

# In[25]:


image_float+230


# It is a very common error to use the wrong image type for a certain application. When you do series of operations, it's always a good idea to regularly look at the image itself to catch such problems that remain "silent" (no error message).

# ### 2.2.2 Logical operations

# A set of important operations when processing images are logical (or boolean) operations that allow e.g. to create object masks. Those have a very simple syntax in Numpy. For example, let's compare pixel intensities to some value ```threshold```:

# In[26]:


threshold = 100


# In[27]:


image > threshold


# We see that the result is again a pixel-wise comparison with ```threshold```, generating in the end a boolean or logical matrix. We can directly assign this logical matrix to a variable and verify its shape and type and plot it:

# In[28]:


image_threshold = image > threshold


# In[29]:


image_threshold.shape


# In[30]:


image_threshold.dtype


# In[31]:


image_threshold


# In[32]:


plt.imshow(image_threshold);


# Of course other logical operator can be used (<, >, ==, !=) and the resulting boolean matrices combined:

# In[33]:


threshold1 = 70
threshold2 = 100
image_threshold1 = image > threshold1
image_threshold2 = image < threshold2


# In[34]:


image_AND = image_threshold1 & image_threshold2
image_XOR = image_threshold1 ^ image_threshold2


# In[35]:


fig, ax = plt.subplots(nrows=1,ncols=4,figsize=(15,15))
ax[0].imshow(image_threshold1)
ax[1].imshow(image_threshold2)
ax[2].imshow(image_AND)
ax[3].imshow(image_XOR);


# ## 2.3 Numpy functions
# 
# To broadly summarize, one can say that Numpy offers four types of operations: 1. Creation of various types of arrays, 2. Pixel-wise modifications of arrays, 3. Operations changing array dimensions, 4. Combinations of arrays.

# ### 2.3.1 Array creation

# Often we are going to create new arrays that later transform them. Functions creating arrays usually take arguments specifying both the content of the array and its dimensions. 
# 
# Some of the most useful functions create 1D arrays of ordered values. For example to create a sequence of numbers separated by a given step size:

# In[36]:


np.arange(0,20,2)


# Or to create an array with a given number of equidistant values:

# In[37]:


np.linspace(0,20,5)


# In higher dimensions, the simplest example is the creation of arrays full of ones or zeros. In that case one only has to specify the dimensions. For example to create a 3x5 array of zeros:

# In[38]:


np.zeros((3,5))


# Same for an array filled with ones:

# In[39]:


np.ones((3,5))


# Until now we have only created one-dimensional lists of 2D arrays. However Numpy is designed to work with arrays of arbitrary dimensions. For example we can easily create a three-dimensional "ones-array" of dimension 5x8x4:

# In[40]:


array3D = np.ones((2,6,5))


# In[41]:


array3D


# In[42]:


array3D.shape


# And all operations that we have seen until now and the following ones apply to such high-dimensional arrays exactly in the same way as before:

# In[43]:


array3D*5


# We can also create more complex arrays. For example an array filled with numbers drawn from a normal distribution:

# In[44]:


np.random.standard_normal((3,5))


# As mentioned before, some array-creating functions take additional arguments. For example we can draw samples from a gaussian distribution whose mean and variance we can specify. 

# In[45]:


np.random.normal(loc=10, scale=2, size=(5,2))


# ### 2.3.2 Pixel-wise operations

# Numpy has a large trove of functions to do all common mathematical operations matrix-wise. For example you can take the cosine of all elements of your array:

# In[46]:


angles = np.random.random_sample(5)
angles


# In[47]:


np.cos(angles)


# Or to calculate exponential values:

# In[48]:


np.exp(angles)


# And many many more. If you need any specific one, look for it in the docs or just Google it.

# ### 2.2.3 Operations changing dimensions

# Some functions are accessible in the form of *methods*, i.e. they are called using the *dot-notation*. For example to find the maximum in an array:

# In[49]:


angles.max()


# Alternatively there's also a maximum function:

# In[50]:


np.max(angles)


# The ```max``` function like many others (min, mean, median etc.) can also be applied to a given axis. Let's imagine we have an image of dimensions 10x10x4 (e.g. 10 pixels high, 10 pixels wide and 4 planes):

# In[51]:


volume = np.random.random((10,10,4))


# If we want to do a maximum projection along the third axis, we can specify (remember that we start counting at 0 in Python). We expect to get a single 10x10px image out of this:

# In[52]:


projection = np.max(volume, axis = 2)


# In[53]:


projection.shape


# We see that we have indeed a new array with one dimension less because of the projection.

# ### 2.3.4 Combination of arrays

# Finally arrays can be combined in multiple ways. For example if we want to assemble two images with the same size into a stack, we can use the ```np.stack``` function:

# In[54]:


image1 = np.ones((4,4))
image2 = np.zeros((4,4))

stack = np.stack([image1, image2],axis = 2)


# In[55]:


stack.shape


# ## 2.3 Slicing and indexing

# Just like broadcasting, the selection of parts of arrays by slicing or indexing can become very sophisticated. We present here only the very basics to avoid confusion. There are often multiple ways to do slicing/indexing and we favor here easier to understand but sometimes less efficient solutions.
# 
# Let's look at a painting (Paul Klee, "Übermut"):

# In[56]:


image = skimage.io.imread('../Data/Klee.jpg')


# In[57]:


image.shape


# We see that the image has three dimensions, probably it's a stack of three images of size 643x471. Let us try to have a look at this image hoping that dimensions are handled gracefully:

# In[58]:


plt.imshow(image);


# Whenever an image has three planes in its third dimension like here, matplotlib assumes that it is dealing with a natural image where planes encode the colors Red, Green and Blue (RGB) and displays it accordingly.

# ### 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. For each dimension we can:
# - give a single index e.g. ```0``` to recover the first element
# - give a range e.g. ```0:10``` to recover the first 10 elements
# - recover **all** elements using the sign ```:```
# 
# In our particular case, we want to select all rows ```:``` for the first axis, all columns ```:``` for the second axis, and a single plane ```1``` for the third axis of the image:

# In[59]:


image[:,:,1].shape


# In[60]:


plt.imshow(image[:,:,0],cmap='gray')
plt.title('First plane: Red');


# We see now the red layer of the image, and indeed the bright regions in this gray scale image correspond with red regions. We can do the same for the others by specifying planes 0, 1, and 2:

# In[61]:


fig, ax = plt.subplots(nrows=1, ncols=4, figsize=(10,10))
ax[0].imshow(image[:,:,0],cmap='gray')
ax[0].set_title('First plane: Red')
ax[1].imshow(image[:,:,1],cmap='gray')
ax[1].set_title('Second plane: Green')
ax[2].imshow(image[:,:,2],cmap='gray')
ax[2].set_title('Third plane: Blue');
ax[3].imshow(image)
ax[3].set_title('Complete image');


# Logically intensities are high for the red channel (red, brown regions) and low for the blue channel (few blue shapes). 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 single sliced plane:

# In[62]:


image0 = image[:,:,0]


# In[63]:


np.mean(image0)


# For all planes, we can either use a comprehension list to go through all planes:

# In[64]:


[np.mean(image[:,:,i]) for i in range(3)]


# Or more elegantly, calculate the mean over the axis 0 and 1:

# In[65]:


np.mean(image,axis = (0,1))


# To look at some more details let us focus on a smaller portion of the image e.g. the head of the character. 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 170 to 300 and columns from 200 to 330 of the first plane (0).

# In[66]:


image_red = image[170:300,200:330,0]
plt.imshow(image_red,cmap='gray');


# 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:

# In[67]:


image_subsample = image[170:300:3,200:330:3,0]
plt.imshow(image_subsample,cmap='gray');


# ### 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:

# In[68]:


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]


# 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 below a threshold:

# In[69]:


threshold_image = image_red<120


# Let's visualize it. Matplotlib handles logical arrays simply as a binary image:

# In[70]:


plt.imshow(threshold_image)


# We can recover the value of all the "white" (True) pixels in the original image by **indexing one array with the other**:

# In[71]:


selected_pixels = image_red[threshold_image]
print(selected_pixels)


# And now ask how many pixels are above threshold and what their average value is.

# In[72]:


len(selected_pixels)


# In[73]:


np.mean(selected_pixels)


# We now know that there are 5054 pixels above the threshold and that their mean is 64.5.