<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Introduction to Correlation and Colocalisation with Python.
Reading images
Dominic Waithe 2015 (c)
Exercise
Step1:  The Pearson's test
Exercise
Step2:  Pearson's comparison of microscopy derived images
Step3:  Maybe remove so not to clash with Mark's.
Last challenge
Exercise
<ASSISTANT_TASK:>
Python Code:
#This line is very important: (It turns on the inline visuals!)
%pylab inline
a = [2,9,32,12,14,6,9,23,4,5,13,6,7,92,21,45];
b = [7,21,4,2,92,9,9,6,13,12,45,5,6,23,14,32];
#Please calculate the dot product of the vectors 'a' and 'b'.
#You may use any method you like. If get stuck. Check:
#http://docs.scipy.org/doc/numpy/reference/generated/numpy.dot.html
#If you rearrange the numbers in 'b', what sequence will give
#the highest dot-product magnitude?
Explanation: Introduction to Correlation and Colocalisation with Python.
Reading images
Dominic Waithe 2015 (c)
Exercise: See the similarities between the dot-product and correlation. Apply correlation to images to obtain a metric of colocalisation/similarity. Use colocalisation to assess the quality of registration.
We start with two lists of numbers (or two vectors or arrays as they are known). Please find the dot product of the two vectors. The dot product formula is a follows:<img src="dotProduct.png">
In python there is more than one way to find the dot product of two vectors. It can be performed using 'for loops' or through vectorised notation
End of explanation
#The cross-correlation algorithm is another name for the Pearson's test.
#Here it is written in code form and utilising the builtin functions:
c = [0,1,2]
d = [3,4,5]
rho = np.average((c-np.average(c))*(d-np.average(d)))/(np.std(c)*np.std(d))
print('rho',round(rho,3))
#equally you can write
rho = np.dot(c-np.average(c),d-np.average(d))/sqrt(((np.dot(c-np.average(c),c-np.average(c)))*np.dot(d-np.average(d),d-np.average(d))))
print('rho',round(rho,3))
#Why is the rho for c and d, 1.0?
#Edit the variables c and d and find the pearson's value for 'a' and 'b'.
#What happens when you correlate 'a' with 'a'?
#Here is an image from the Fiji practical
from tifffile import imread as imreadtiff
im = imreadtiff('neuron.tif')
print('image dimensions',im.shape, ' im dtype:',im.dtype)
subplot(2,2,1)
imshow(im[0,:,:],cmap='Blues_r')
subplot(2,2,2)
imshow(im[1,:,:],cmap='Greens_r')
subplot(2,2,3)
imshow(im[2,:,:],cmap='Greys_r')
subplot(2,2,4)
imshow(im[3,:,:],cmap='Reds_r')
Explanation: The Pearson's test
Exercise: See the similarities
The above example shows you how two number sequences can be compared with nothing more complicated than by using the dot product. This works as long as the sequences comprise of the same numbers but in a shuffled order.  To compare different sequences with the original we normalise by the magnitude of the vectors. To include this step. We use a more complicated equation:
<img src="eqn_full.gif">
https://en.wikipedia.org/wiki/Pearson_product-moment_correlation_coefficient
https://en.wikipedia.org/wiki/Cross-correlation
Hopefully you can see the top of this equation is very similar to the dot-product, except that it is centered on zero (subtraction of the mu, the mean) and the variance is normalised (division by standard deviation).
Because the equation is normalised, a perfectly correlated sequence yeilds a rho value of 1.0. A perfectly random comparison yields 0 and two anti-correlated sequences will yield a value of -1.0.
End of explanation
a = im[0,:,:].reshape(-1)
b = im[3,:,:].reshape(-1)
#Calculate the pearson's coefficent (rho) for the image channel 0, 3.
#You should hopefully obtain a value 0.829
#from tifffile import imread as imreadtiff
im = imreadtiff('composite.tif')
#The organisation of this file is not simple. It is also a 16-bit image.
print("shape of im: ",im.shape,"bit-depth: ",im.dtype)
#We can assess the image data like so.
CH0 = im[0,0,:,:]
CH1 = im[1,0,:,:]
#Single channels visualisation can handle 16-bit
subplot(2,2,1)
imshow(CH0,cmap='Reds_r')
subplot(2,2,2)
imshow(CH1,cmap='Greens_r')
subplot(2,2,3)
#RGB data have to range between 0 and 255 in each channel and be int (8-bit).
imRGB = np.zeros((CH0.shape[0],CH0.shape[1],3))
imRGB[:,:,0] = CH0/255.0
imRGB[:,:,1] = CH1/255.0
imshow(255-imRGB.astype(int8))
#What is the current Pearson's value for this image?
Explanation: Pearson's comparison of microscopy derived images
End of explanation
rho_max = 0
#This moves one of your images with respect to the other.
for c in range(1,40):
    for r in range(1,40):
        #We need to dynamically sample our image.
        temp = CH0[c:-40+c,r:-40+r].reshape(-1);
        #The -40 makes sure they are the same size.
        ref = CH1[:-40,:-40].reshape(-1);
        
        rho = np.dot(temp-np.average(temp),ref-np.average(ref))/sqrt(((np.dot(temp-np.average(temp),temp-np.average(temp)))*np.dot(ref-np.average(ref),ref-np.average(ref))))
        
        #You will need to work out the highest rho value is recorded.
        #You will then need to find the coordinates of this high rho.
        #You will then need to provide a visualisation with the image translated.
        
Explanation: Maybe remove so not to clash with Mark's.
Last challenge
Exercise: The above image is not registered. Can you devise a way of registering this image using the Pearson's test, as a measure for the similarity of the image in different positions. hint you will need to move one of the images relative to the other and measure the colocalisation in this position. The best localisation will have the highest rho value. Produce an image of your fully registered image.
End of explanation
<END_TASK>