Diffraction image visualization

This exercise is about single crystal diffraction images (protein crystallography) and highlights difficulties one can encounter when trying to display such image.

The image we will try to display, run2_1_00147.cbf, is part of the insulin test-set provided by the manufacturer (Dectris) to advertise the quality of its Pilatis 6M detectors.

The first difficulty one may encounter is the compressed form of this image. Fortunately, both PyMca and FabIO support the CIF Binary Format and the byte-offset compression scheme used by Dectris.



In [2]:

    
%pylab inline









    



Populating the interactive namespace from numpy and matplotlib



In [3]:

    
import fabio
img = fabio.open("run2_1_00148.cbf")
imshow(img.data)









    Out[3]:





<matplotlib.image.AxesImage at 0x7f3ee40d9e10>

When trying to display such image, it appears empty. One has to check the content to be sure it is not:



In [4]:

    
print(img.data.shape)
print(img.data.max())
print(img.data.min())
print(img.data.mean())
print(img.data.std())









    



(2527, 2463)
664762
-2
6.4045405841
473.992411656

Diffraction data exhibit a very large dynamical range, here from 0 to 7e5 but most datapoint are arround 0. One solution, at least, one part of the solution is to display the image in a logarithmic color-scale.

Logarithmic colors

Simply display the logarithm of the intensity (intensities <=0 will raise a warning and be transparent)



In [5]:

    
imshow(np.log(img.data))









    



-c:1: RuntimeWarning: divide by zero encountered in log
-c:1: RuntimeWarning: invalid value encountered in log






    Out[5]:





<matplotlib.image.AxesImage at 0x7f3ee402bb90>

Taking the log of the image helps a lot: one clearly sees now the 5x12 Pilatus modules with their gaps, the beamstop and the diffuse scattering attributed to water (ice-ring). But no diffraction spots are visible, actually a few of them, but not enough to assess the quality of diffraction signal.

Here the images displayed in the IPython Notebook is about 326x333 while the input image was 2463x2527, so one pixel on the screen is representing 10x10 pixels on the detector. Because the Pilatus is a pixel-detector with a point spread function of 1 pixel, only 1% of the information is displayed. Two solution exists::

* make peaks 10x larger
* bin the image by a factor 10

The latter solution has already been seen and is called "binning"

Image binning

To bin the image, it needs to be a full multiple of the binning factor, here 10. Moreover to avoid side effects with integers, we will convert the to float. We use here one of the two binning techics developped in the numpy exercise.



In [6]:

    
chunk = img.data[:2520,:2460].astype("float")
chunk.shape = 252, 10, 246, 10
binned = chunk.sum(axis=-1).sum(axis=1)
imshow(np.log(binned))









    



-c:4: RuntimeWarning: divide by zero encountered in log
-c:4: RuntimeWarning: invalid value encountered in log






    Out[6]:





<matplotlib.image.AxesImage at 0x7f3ee3f61990>

Note the gaps have almost disapeared due to the binning of the image ... but a dozen of diffraction spots are now clearly visible.

Increase the size of spots

This can be done via "grey-scale" morphology from scipy.ndimage. More details about the mathematics behind grayscale morphology can be found in: http://en.wikipedia.org/wiki/Mathematical_morphology We will use "grey_dilation" which actually will transform a peak of 1 pixel wide into a disc of 10 pixels.



In [7]:

    
from scipy import ndimage



In [8]:

    
imshow(np.log(ndimage.grey_dilation(img.data,10)))









    Out[8]:





<matplotlib.image.AxesImage at 0x7f3ee3e2e9d0>

In this picture, one can clearly see hundreeds of diffraction spots up to a very wide diffraction angle (exhibiting the very good resolution for the protein diffractionist).

The last step is now to resume the masked pixelsfrom the original image:



In [9]:

    
dil = ndimage.grey_dilation(img.data,10)
dil[img.data<0] = 0
imshow(np.log(dil))









    



-c:3: RuntimeWarning: divide by zero encountered in log






    Out[9]:





<matplotlib.image.AxesImage at 0x7f3ee3dc6710>

Alternative solution

Instead of using sum in the binning use the max to highlight peaks



In [10]:

    
chunk = img.data[:2520,:2460]
chunk.shape = 252, 10, 246, 10
binned = chunk.max(axis=-1).max(axis=1)
imshow(np.log(binned))









    Out[10]:





<matplotlib.image.AxesImage at 0x7f3ed7162fd0>

Conclusion

With such thumbnail image, of only 300x300, a protein crystalographer can actually assess the quality of the data-collection ongoing. Such technique is implemented in an EDNA plugin to generate the thumbnail displayed into ISPyB.



In [10]:



In [10]: