A collection of resources and functions for performing image analysis and making figures in python.
Switch branches/tags
Nothing to show
Clone or download
Fetching latest commit…
Cannot retrieve the latest commit at this time.
Permalink
Failed to load latest commit information.
img
.gitignore
README.org
imageanalysis.py

README.org

Python functions for figures and analysis

This document is a work in progress and contains useful python functions I have discovered or written to aid in image analysis of microscopy data and making figures for presentations. Most of my work is done in python or R and I use Emacs Org-mode for literate programming. My goal for writing scientific papers is to work in org-mode and have all my figures generated programmatically when the document compiles to PDF. This makes the work more reproducible and easier to work with. I will be adding to this document as I work and learn. Most of it will not be org-mode specific, except where indicated.

This document is tangled into the accompanying python source doc in the repo and can be downloaded and used freely.

As a module

I am working on integrating these into a PyPi module. I made a class to contain images and the accompanying metadata which is here. I will keep updating this as I go/need.

Suggestions and contributions to either project are welcome!

Links and resources

Here is a collection of links and resources I have come across.

Making figures

# Emacs-plot venv
import os
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
if not 'img' in os.listdir("."):
    os.mkdir('img')

Basics

The first step is to read images. For this, we use a python wrapper around the OME Bioformats command line tool bfconvert. I wrote a python wrapper around bfconvert that will convert images or directories and store them in a dedicated folder for easier use and organization. This script and the documentation is available on github.

Next, images are read in using skimage like so

# using scikit-image to read images

import skimage.io
from skimage import img_as_float

def read_img(path, asfloat=False):
    img = skimage.io.imread(path)
    if not asfloat:
        return img
    img = img_as_float(img)
    return img

Parse metadata

Images are converted to OME-Tiff format because this is an open source data format that preserves all metadata and can be opened by virtually any image reader. The following class is from my imagetools repo. It uses the tifffile library to read the image metadata from the first page of the tiff file.

# main class 

import re
import tifffile
import skimage.io


class ImageInfo:
    """
    Class to hold tiff image and metadata.
    Arguments:
    ==========
    path: path to a valid tiff
    ==========
    Attempts to parse OME metadata. If it cannot, then it tries to
    parse imagej format. If this fails, None is returned for these
    attributes.

    """
    def __init__(self, path):
        self.__path = path
        self.__meta_page = self.parse_tif()
        self.__units_and_len = self.parse_ome_metadata()
        self.__pixel_size = self.__units_and_len[0]
        self.__unit = self.__units_and_len[1]
        self.__image = self.read_image()
        self.__shape = self.image.shape


    def parse_tif(self):
        with tifffile.TiffFile(self.__path) as tif:
            meta_page = tif[0]
            return meta_page


    def parse_ome_metadata(self):
        """accepts the path to an ome-tif file, or imageJ tif file.
        Attempts to validates image dimensions, returns pixel size
        and
        """
        meta_page = self.__meta_page
        omeXY = re.compile(r'PhysicalSize[XY]\s*\=\s*\"(\d+\.\d+)\"', re.I)
        omeXY_units = re.compile(r'PhysicalSize[XY]Unit\s*\=\s*\"(\D+)\"\s', re.I)
        blob = meta_page.image_description.decode('utf-8')
        findomeUnits = omeXY_units.findall(blob)
        XY = omeXY.findall(blob)
        if len(XY) and len(findomeUnits) == 2:
            try:
                assert XY[0] == XY[1]
                unit_as_float = float(XY[0].strip(' "'))
                return unit_as_float, findomeUnits[0]
            except AssertionError:
                print(f'OME parsing X resolution {XY[0]}!= {XY[1]}, returning None')
                return None, None
        else:
            print('Ome data not found, attmepting imagej parse')
            return self.parse_imagej_meta()


    def parse_imagej_meta(self):
        """
        Only called if Ome parsing fails
        method to parse imagej formats
        """
        print('WARNING ImagJ parsing is not as accurate!!')
        meta_page = self.__meta_page
        blob = self.__meta_page.image_description.decode('utf-8')
        try:
            XY = meta_page.x_resolution, meta_page.y_resolution
            find_ImageJ_units = re.findall(r'unit=(.+)',blob)
            if find_ImageJ_units[0] == 'micron':
                find_ImageJ_units = 'µm'
                if len(find_ImageJ_units) == 0:
                    print('ImageJ units were not found. Returning None')
                    find_ImageJ_units = None
                assert XY[0] == XY[1]
                return XY[0][0]/XY[0][1], find_ImageJ_units
        except AssertionError:
            print(f'ImageJ parsing X resolution {XY[0]}!={XY[1]}, returning None')
            return None, None
        except AttributeError:
            print('Could not parse, returning None,None')
            return None, None
        except Exception as e:
            print(f'unknown exception {e}')
            return None, None


    def read_image(self):
        """
        read image using skimage
        """
        return skimage.io.imread(self.__path)

    @property
    def image(self):
        return self.__image


    @property
    def pixel_size(self):
        """
        return the pixel size
        """
        return self.__pixel_size

    @property
    def pixel_unit(self):
        """
        return the pixel unit
        """
        return self.__unit

    @property
    def get_path(self):
        """
        return the original image path
        """
        return self.__path

    @property
    def get_meta_blob(self):
        """
        return the metadata blob taken from
        tifffile.TiffFile page 0 .image_description attribute
        """
        blob = self.__meta_page.image_description.decode('utf-8')

I recently updated this function to also parse ImageJ encoded tifs. When you save a tif from imageJ, it encodes only the essential info in the file. By parsing the tif as explained in the documentation like so:

# example of metadata returned form an imageJ tif

import tifffile  
neun_path_example = '/Volumes/EXTENSION/RESTREPOLAB/images/neuronavigation/macklin_zeiss/2017-08-01/figures/MAX_2017-08-01_H001-017_img006.tif'
with tifffile.TiffFile(neun_path_example) as tif:
    images = tif.asarray()
    for page in tif:
        for tag in page.tags.values():
            t = tag.name, tag.value
            print(t)

You see a different set of results. So in ome-tif files, all the metadata is a xml blob in the image_description tag, while in a ImageJ encoded tif, you have to extract it from a series of top level tags. It turns out that these tags exist in the ome-tif ones too, but I think the xml is better to stick with if it is available. The imageJ parsing is less robust and informative, I am not exactly sure what all the tags mean, but I spot check it and it works for now. I think spacing in image_description refers to z-step size but I am not sure. Anyways this seems to work for parsing two tif encoding formats.

add scalebar

use the matplotlib-scalebar class. Here are some common arguments I like.

# example of a scalebar

import matplotlib.pyplot as plt
from matplotlib_scalebar.scalebar import ScaleBar

scalebar = ScaleBar(pixelLength, units, location = 'lower right', 
                   fixed_value = 25, color = 'black', frameon = False)

Here is a function to use it in a figure.

# function for plotting an image with a scalebar

import matplotlib.pyplot as plt
from matplotlib_scalebar.scalebar import ScaleBar

def scale_plot(img, imageSize, scale, units, scalebar_length, color):
    plt.figure(figsize=imageSize)
    plt.imshow(img)
    plt.axis('off')
    scalebar = ScaleBar(scale, units, location = 'lower right', 
                        fixed_value = scalebar_length, color = color, frameon = False)
    plt.gca().add_artist(scalebar)

panel figures

## example data
from matplotlib_scalebar.scalebar import ScaleBar
import numpy as np # for image example

np.random.seed(0)
example_image = np.random.rand(250,250,3)

fig = plt.figure(figsize=(10,10))
fig.suptitle('Experiment title here', fontsize=15)
one = fig.add_subplot(131)
plt.imshow(example_image[:,:,0])
one.set_title('Time one', fontsize=15)
plt.axis('off')
two = fig.add_subplot(132)
two.set_title('Time two', fontsize=15)
plt.imshow(example_image[:,:,1], cmap='gray')
plt.axis('off')
three = fig.add_subplot(133)
plt.imshow(example_image[:,:,2])
three.set_title('Time three', fontsize=15)
scalebar = ScaleBar(1,'um',location='lower right', fixed_value=25, color='black',frameon=True)
plt.gca().add_artist(scalebar)
plt.axis('off')
#plt.tight_layout()
plt.subplots_adjust(wspace=0.01)
plt.subplots_adjust(top=1.35)
plt.savefig('img/three_panel.png',bbox_inches='tight')

Gives the following figure

img/three_panel.png

cooler

import skimage.measure
# make line profiles

np.random.seed(0)
example_image = np.random.rand(250,250,3)

start_coords = [50,50]
stop_coords = [150,150]

# make line profiles


start_y_line = skimage.measure.profile_line(example_image[:,:,0], start_coords, stop_coords)
middle_y_line = skimage.measure.profile_line(example_image[:,:,1], start_coords, stop_coords)
last_y_line = skimage.measure.profile_line(example_image[:,:,2], start_coords, stop_coords)
linescan_dist = (np.linalg.norm(np.array(start_coords) - np.array(stop_coords)))
line_axis = np.linspace(0,linescan_dist+1,len(start_y_line))

# column 1
fig = plt.figure(figsize=(10,8))
fig.suptitle('1040nm exposure', fontsize=15)
one = fig.add_subplot(231)
plt.imshow(example_image[:,:,0], cmap='gray')
plt.plot([start_coords[0],stop_coords[0]], [start_coords[1],stop_coords[1]],'r-', linewidth=4)
one.set_title('Start exposure', fontsize=15)
plt.axis('off')
onescan = fig.add_subplot(234)
plt.plot(line_axis, start_y_line,'-', color='black')
onescan.spines['right'].set_visible(False)
onescan.spines['top'].set_visible(False)
plt.ylabel('Fluorescence intensity (AU)', fontsize=15)
plt.xlabel(r'Distance ($\mu{}m$)', fontsize=15)


#column2
two = fig.add_subplot(232)
two.set_title('22 s', fontsize=15)
plt.imshow(example_image[:,:,1], cmap='gray')
plt.plot([start_coords[0],stop_coords[0]], [start_coords[1],stop_coords[1]],'r-', linewidth=4)
plt.axis('off')
middlescan = fig.add_subplot(235)
plt.plot(line_axis, middle_y_line,'-', color='black')
plt.axis('off')


#column3
three = fig.add_subplot(233)
plt.imshow(example_image[:,:,2], cmap='gray')
plt.plot([start_coords[0],stop_coords[0]], [start_coords[1],stop_coords[1]],'r-', linewidth=4)
three.set_title('45 s', fontsize=15)
plt.axis('off')
scalebar = ScaleBar(1,'um',location='lower right', fixed_value=25,color = 'black', frameon=True)
plt.gca().add_artist(scalebar)
lastscan = fig.add_subplot(236)
plt.plot(line_axis, last_y_line, '-', color='black')
plt.axis('off')
plt.subplots_adjust(wspace=0.01)
plt.subplots_adjust(top=.9)
fig.savefig('img/three_panel_with_scans.png', bbox_inches='tight')

That code generates the figure below

img/three_panel_with_scans.png

In a row

np.random.seed(0)
example_image = np.random.rand(250,250,3)

roi_stim_coords_start = [50,50]
roi_stim_coords_end = [150,150]

# make line profiles


pre_exposure_y_line = skimage.measure.profile_line(example_image[:,:,0], start_coords, stop_coords)
post_exposure_y_line = skimage.measure.profile_line(example_image[:,:,1], start_coords, stop_coords)
linescan_dist = (np.linalg.norm(np.array(roi_stim_coords_start) - np.array(roi_stim_coords_end)))
stim_line_axis = np.linspace(0,linescan_dist+1,len(pre_exposure_y_line))

# TP 1
fig = plt.figure(figsize=(25,10))
one = fig.add_subplot(141)
plt.imshow(example_image[:,:,0], cmap='gray')
one.set_title('Pre-exposure',fontsize=20)
plt.plot([roi_stim_coords_start[0],roi_stim_coords_end[0]], 
         [roi_stim_coords_start[1],roi_stim_coords_end[1]], 'r', linewidth=4)
plt.axis('off')

# TP 2
two = fig.add_subplot(142)
two.set_title('During exposure',fontsize=20)
plt.imshow(example_image[:,:,1], cmap='gray')
plt.plot([roi_stim_coords_start[0],roi_stim_coords_end[0]], 
        [roi_stim_coords_start[1],roi_stim_coords_end[1]], 'r', linewidth=4)
plt.axis('off')

# TP 3 
three = fig.add_subplot(143)
plt.imshow(example_image[:,:,2], cmap='gray')
three.set_title('post-exposure',fontsize=20)
plt.axis('off')
scalebar = ScaleBar(1,'um',location='lower left', fixed_value=25,color = 'black', frameon=True)
plt.gca().add_artist(scalebar)
plt.plot([roi_stim_coords_start[0],roi_stim_coords_end[0]], 
         [roi_stim_coords_start[1],roi_stim_coords_end[1]], 'r', linewidth=4)

# linescans
four = fig.add_subplot(144)
plt.plot(stim_line_axis, pre_exposure_y_line, '--', color='blue',linewidth=2, label='pre-exposure')
plt.plot(stim_line_axis, post_exposure_y_line, '-', color='black',linewidth=2, label='post-exposure')
four.spines['right'].set_visible(False)
four.spines['top'].set_visible(False)
four.legend(loc='lower center', fontsize=15)
plt.ylabel('Fluorescence intensity (a.u.)',fontsize=15)
plt.xlabel(r'Distance ($\mu{}m$)',fontsize=15)
plt.subplots_adjust(wspace=None)
fig.savefig('img/in_a_row.png', bbox_inches='tight')

How about in a row?

img/in_a_row.png

Subplots

The following format works well for splitting two channels and merging. Some of my microscopy images are only two channels. To plot these, I add an extra channel of zeros. To do this, do the following:

# create three channel image from 2 channel


import numpy as np

two_channel_image = np.random.rand(250,250,2)
print('Original image shape is {}'.format(two_channel_image.shape))

# make it three

now_three =np.dstack((two_channel_image[:,:,0], two_channel_image[:,:,1],
                   np.zeros_like(two_channel_image[:,:,0])))
print('Three channel image shape is {}'.format(now_three.shape))
two_channel_image = np.random.rand(250,250,2)
now_three =np.dstack((two_channel_image[:,:,0], two_channel_image[:,:,1], 
                      np.zeros_like(two_channel_image[:,:,0])))

# plot 2 channels of an image with scalebar
fig = plt.figure(figsize=(10,10))
one = fig.add_subplot(131)
plt.imshow(now_three[:,:,0], cmap="Greens_r") # note colormap
one.axis('off')
one.set_title('Channel 1',size=15)
two = fig.add_subplot(132)
plt.imshow(now_three[:,:,1] ,cmap="Reds_r") # note colormap
two.set_title('Channel 2',size=15)
two.axis('off')
scalebar = ScaleBar(1, units, location = 'lower right', 
                        fixed_value = 25, color = 'black', frameon = True)
three = fig.add_subplot(133)
plt.imshow(now_three)
plt.gca().add_artist(scalebar)
three.set_title('Merge', size=15)
three.axis('off')
plt.tight_layout()
fig.savefig('img/fake_channels.png', bbox_inches='tight')

Easy way to work with channels.

img/fake_channels.png

more to come…

Processing

Max projections

Shocked there was not already a method for this?

## max project

import numpy as np


def max_project(image, start_slice = 0, stop_slice = None):
    """ takes ONE CHANNEL nd array image
        optional args = start_slice, stop_slice
        range to max project into
        returns new projection"""
    if stop_slice is None:
        stop_slice = image.shape[0]
    print(stop_slice)
    max_proj = [image[i,:,:] for i in range(start_slice,stop_slice)]
    return np.maximum.reduce(max_proj)

I will expand this for multichannel images and also implement mean projection, sum projection etc. Example coming soon.

Analysis

Interactive line profiles are cool!

# draw a line profile interactively


import matplotlib
matplotlib.use('TKAgg') # I don't have matplotlib installed as a framework so I need this..
from skimage import data
from skimage.viewer import ImageViewer
from skimage.viewer.plugins.lineprofile import LineProfile

def make_profile(image):
    """ 
    Takes a 2D image, gives an PyQt image
    viewer that you can make a ROI on. 
    returns line profile values
    """
    viewer = ImageViewer(image)
    viewer += LineProfile()
    _, line = zip(*viewer.show())
    return line

Dray your profile line then close the image. This returns a list of length 1 containing a tuple. The tuple contains the image array and the line profile. Very annoying, but I used argument unpacking with zip(*args) to fix it. I used the tutorial and put them in a function.

How similar are your images? You can use structured similarity index (SSIM) to get an indication.

# compare image channels with SSIM

import skimage.measure

score, diff = skimage.measure.compare_ssim(image[0,:,:], image[1,:,:], full=True,
                                           gaussian_weights=True, sigma=1.5, 
                                           use_sample_covariance=False)

a score of 1 is most similar and a score of -1 is least similar. See the skimage documentation and the paper.