Introduction

There are tons of software packages for working with the volumetric image data collected from confocal microscopes. For example Fiji/ImageJ is becoming a clear standard for image processing in biology, and there are also many commercial packages e.g. Metamorph. In spite of (or maybe because of) this richness of resources, a few years ago I ran into a image-processing and visualization problem where I didn't find an out of the box solution that perfectly fit my needs, so I ended up rolling my own. As a little bit of background, we were studying a collection of tiny muscles that attach to the fruit fly wing hinge, and found that if we expressed a fluorescent calcium sensor in the muscles we could visualize their activity straight through the intact cuticle. We wanted to use wide-field microscopy for speed and cost considerations; however, to make sense of the signals coming out of the thorax, we needed a better map of the anatomy in Drosophila than currently existed. I had the good fortune to work with a fantastic histologist, https://www.linkedin.com/in/anne-sustar-437401149> Anne Sustar , on this project and she developed a whole-mount protocol to stain fixed muscles with the fluorescently-labeled actin-binding stain Phalloidin. After clearing the tissue, she imaged the tissue on a confocal, also collecting a series of transmitted light images in register with the confocal slices. This was important because it allowed us to reconstruct the muscle anatomy with respect to easily-identifiable landmarks on the external cuticle. I ended up using these data to build a quantitative model of the muscle morphology using a method I'll describe in another post; however the procedure I used to segment the muscles from one another allowed me to create a useful visualization of the anatomy in python which I describe here.

First, I used the https://imagej.net/Extended_Depth_of_Field> extended depth of field plugin for ImageJ to project the stack of transmitted light images into a sharp image of the cuticular landmarks.

Next, I needed to segmentation each muscle from the surrounding tissue. I briefly spent some time with a few of the ImageJ segmentation tools, but I found that most of these required that I use a lasso tool to select the ROI on each slice. In my case, this was problematic for a couple of reasons. First, drawing the borders around the muscles this in this manner was not terribly mistake-tolerant. Second, sometimes the borders of a muscle in a particular slice was rather tortuous, this was especially true when many individual fibers of a single muscle were present on a single slice. Finally, I didn't like the fact that I was required to make all-or-none decisions about the location of tissue margins: often times this was not entirely clear from the micrographs, and thinking about where to put the border slowed me down. Since I had taken a digital imaging class for photography a while back, I knew that common way that photographers select one part of an image from the rest in a program like photoshop was to use the airbrush tool to paint in the selection on a mask layer. Unfortunately, Photoshop is an expensive bit of software, and I knew that I would want to come up with a system to automate parts of the procedure so I turned to Gimp. It's free, open-source and has a scripting and plug-in architecture.

Tissue segmentation in GIMP

The basic approach I took to segment the tissue was to import the series of confocal images into separate layers, and then create an interleaving set of layers where I would airbrush-in the selection for the underlying layer. There were 73 images in the confocal stack, so performing this operation using only the menu and GUI tools was bound to be very laborious. Without much experience with the gimp API I was able to cook up a few commands for command line that made this process much easier. The console is a bit clunky to use, but I imagine it would not be hard to convert these to scripts that could be quickly executed via a menu item or keyboard shortcut. I followed this basic procedure:

  • Import an 8bit tiff stack of a single channel into gimp: choose "Select All" and "Open pages as Layers" in the dialog box.

  • Convert to RGB (image>mode>RGB)

  • Rename the background image as "Page 0"

  • Save the imported stack as an .xcf (native gimp format) to use as a template.

  • open a python console within gimp: filters>python-fu>Console

  • get a reference to the image object:

>>> img = gimp.image_list()[0]

  • add a blank layer in-between each existing layer note that the naming convention for the interviening layers is defined by the c-string: "lr%s"%(lnum)

g = gimp.pdb w = g.gimp_image_width(img) h = g.gimp_image_height(img) for lnum,layer in enumerate(reversed(img.layers)): g.gimp_image_set_active_layer(img,layer) new_layer = g.gimp_layer_new(img, w, h, RGBA_IMAGE, "lr%s"%(lnum), 1, NORMAL_MODE) img.add_layer(new_layer)

  • Change the opacity of the segmentation layers to 0.2 so you can see what you are segmenting.

for layer in img.layers: if not('Page' in layer.name): layer.opacity = 0.2

  • now start the segmentation process first by hidng all layers

for layer in img.layers: layer.visible = False

  • starting from page 0 make the layers visible until you have reached the deepest slice of tissue you want to segment. Make that page and the segmentation layer (lr..) above it visible.
  • Make sure you have selected the segmentation layer as the active layer and paint over tissue to create a mask. It is best if you pick a single color to use for all the segmentation (I like red). The utility of gimp is that you can now use the all of the drawing tools and brushes to do the segmentation. I found that plugging in a drawing tablet really helped speed things up.

  • Once you have finished segmenting a layer, move up by making the next Page and segmentation layer visible.

  • Saving the segmented stack requires a little gymnastics. First I delete all the 'Page' layers in two steps:

    • First make all layers visible:
    • Then hide layers not containing the string 'Page'

for layer in img.layers: layer.visible = True

for layer in img.layers: if not('Page' in layer.name): layer.visible = False

  • Merge the visible layers (Image>Merge Visible Layers)
  • Delete the merged layer
  • Now it is possible to save the remaining segmentation layers as separate tiff files using this gimp plugin: http://registry.gimp.org/node/28268
Defining masks in gimp.
Example of the 'fuzzy border problem'

Loading the images into Python

Thanks to support from the National Science Foundation, we were able to publish these segmented annatomical stacks on dryad as a companion to our paper in Current Biology. So I'll skip ahead and pull the data out of the oven by downloading from this repository.

In [1]:
from skimage.external import tifffile
import matplotlib.pyplot as plt
import numpy as np
import warnings
import figurefirst as fifi
warnings.filterwarnings('ignore')
%matplotlib inline
In [2]:
from subprocess import call
call(['wget','http://datadryad.org/bitstream/'+\
      'handle/10255/dryad.135236/muscle_segmentation.zip?sequence=1'])
call(['unzip','muscle_segmentation.zip?sequence=1'])
Out[2]:
0

Next I load the segmentation data: stacks of masks that segment each muscle from the rest of the tissue.

In [3]:
base_str = "muscle_segmentation/65G06_muscle_8bit_%(mu)s_layers/"+\
           "65G06_muscle_8bit_%(mu)s_layers%(num)s.tif"

segmented_data = {}
for muscle in ['b1','b2','b3','i1','i2',
               'iii1','iii24','iii3',
               'hg1','hg2','hg3','hg3','hg4']:
    
    f_names = [base_str%{'mu':muscle,'num':i} for i in range(74)]
    segmented_data[muscle] = tifffile.imread(f_names,key = 0)
    #print(muscle + ' loaded')

I then load the phalloidin-stained muscle images. This will be multipled with the mask to form the max intensity image for each muscle.

In [4]:
ph_staining = tifffile.imread('muscle_segmentation/' + \
                              'phalloidin_stained_hemithorax.tif',key = 0)
ph_staining = ph_staining[::-1] #reverse the stack order

Finally want to use the bright field image of the cuticle as a backdrop.

In [5]:
cuticle = tifffile.imread('muscle_segmentation/65G06_brightfield_cuticle.tif',key = 0)
In [6]:
plt.subplot(1,3,1)
plt.imshow(cuticle,cmap = plt.cm.gray)
plt.subplot(1,3,2)
# phalloidin slice #50
plt.imshow(ph_staining[50],cmap = plt.cm.gray)
plt.subplot(1,3,3)
plt.imshow(segmented_data['i1'][50,:,:,3],cmap = plt.cm.gray)
lns = [ax.axis('off') for ax in plt.gcf().axes]

Mixing-in the segmentation masks

My goal is to combine the segment masks with the confocal data to create a maximum intensity image for each muscle that only includes the segmented pixels. Since I created my mask using the airbrush tool in gimp, each slice of my mask has fuzzy margins that tell where the muscle begins and end. This is either a good thing or a bad thing depending on your perspective: on one hand my goal in segmenting the data is to define boundries would seem appropriate/ On the other hand, in many cases those bountries are not entierly clear from the images, so my segmentation shoud reflect that. It always possible to convert a greyscale image to binary image by setting a threshold, for example: binary_img = gray. So, I will retain the uncertanty in the segment boundries until it is nessesary to throw away that information. For the purposes of the visualization, I have found it is usefull to soften up the edges even further by applying a gaussian filter to the mask data.

In [7]:
from scipy.ndimage.filters import gaussian_filter
In [8]:
plt.subplot(1,2,1)
plt.imshow(segmented_data['i1'][50,:,:,3])
plt.subplot(1,2,2)
#the edge blur is more apparent with a non grayscale colormap.
plt.imshow(gaussian_filter(segmented_data['i1'][50,:,:,3],sigma = 15.))
Out[8]:

Now I want to create a stack with just the phalloidin staining of i1. To do this I elementwize multiply the phalloidin stack with the blurred stack of masks, but I first need to apply the gaussian filter to the masks. The gaussian_filter function I applied from scipy.ndimage.filters works in n-dimensions and can use an independent sigma for each dimension. I will primarily blurr in x and y but also apply a little across z.

In [9]:
blurred_mask = gaussian_filter(segmented_data['i1'][...,3],
                               sigma = (1.,15.0,15.0))
blurred_mask = blurred_mask.astype(float)
plt.imshow(blurred_mask[50])
#looks a little different
Out[9]:

I create the max intensity projection by applying np.max accross the product of the two stacks.

In [10]:
max_projection = np.max(blurred_mask*ph_staining,axis = 0)
max_projection /= np.max(max_projection)
plt.imshow(max_projection,cmap = plt.cm.gray)
Out[10]:
In [11]:
#Here is a closeup of the origin of the i1 filtered with and 
#without the gamma parameters.
plt.subplot(1,2,1)
plt.imshow(max_projection[650:850,100:300],cmap = plt.cm.gray)
plt.subplot(1,2,2)
sharp_max_projection = np.max(blurred_mask**5 *ph_staining,axis = 0)
sharp_max_projection /= np.max(sharp_max_projection)
plt.imshow(sharp_max_projection[650:850,100:300],cmap = plt.cm.gray)
Out[11]:
In [12]:
def max_segmented(muscle = 'i1',sig_z = 1.0,sig_xy = 5.0):
    blurred_mask = gaussian_filter(segmented_data[muscle][...,3],
                               sigma = (sig_z,sig_xy,sig_xy))
    blurred_mask = blurred_mask.astype(float)**8
    max_projection = np.max(blurred_mask*ph_staining,axis = 0)
    max_projection /= np.max(max_projection)
    alpha = np.max(blurred_mask,axis = 0)
    alpha -= np.min(alpha)
    alpha /= np.max(alpha)
    return max_projection,alpha

#maximg,alpha = max_segmented('i1')
#plt.imshow(maximg,cmap = plt.cm.gray)

Colorizing the histology

The matplotlib colomaps are a nice way to colorize the grayscale images but I need to rescale the values from 0-1 first. Once I've created a color image I want to use the values from segmentation mask to set the transparency. In this way if I plot the muscle over the cuticle image only the segmented pixels get included in the overlay.

I will keep the vertical and horizontal scale, but transform from pixels to um using a pixel pitch of 0.76 um.

In [15]:
layout = fifi.FigureLayout('thorax_img_layout.svg',
                           make_mplfigures=True)
ax = layout.axes['img']

PPCH = 0.76

#plt.figure(figsize = (10,10))
xextent = PPCH* np.shape(cuticle)[0]
yextent = PPCH* np.shape(cuticle)[1]
xtnt = [0,xextent,0,yextent]

ax.imshow(cuticle,
           cmap = plt.cm.gray,
           clim = (0,500),extent = xtnt)


def plot_muscle(muscle,
                sig_z = 0.2,
                sig_xy = 1.0,
                cmap = plt.cm.Blues_r,
                ax = plt.gca()):
    maximg,alpha = max_segmented(muscle,
                                 sig_z = sig_z,
                                 sig_xy = sig_xy)
    muscle_color = cmap(maximg)
    muscle_color[...,3] = 0.8*alpha
    ax.imshow(muscle_color,extent = xtnt)

plot_muscle('hg1',cmap = plt.cm.Oranges_r)

plot_muscle('i1',cmap = plt.cm.Purples_r)
plot_muscle('i2',cmap = plt.cm.Purples_r)

plot_muscle('iii1',cmap = plt.cm.Reds_r)
plot_muscle('iii24',cmap = plt.cm.Reds_r)
plot_muscle('iii3',cmap = plt.cm.Reds_r)

plot_muscle('b1',cmap = plt.cm.Blues_r)
plot_muscle('b3',cmap = plt.cm.Blues_r)
plot_muscle('b2',cmap = plt.cm.Blues_r)

plot_muscle('hg2',cmap = plt.cm.Oranges_r)
plot_muscle('hg3',cmap = plt.cm.Oranges_r)
plot_muscle('hg4',cmap = plt.cm.Oranges_r)

ax.set_xlabel('uM')

fifi.mpl_functions.set_spines(layout)
In [16]:
from IPython.display import display,SVG
In [17]:
layout.save('thorax_img.svg')
display(SVG('thorax_img.svg'))
image/svg+xml b2 i1 iii1 iii3 hg2 hg3 hg4 b1 b3 iii4 i2 hg1
In [ ]:
 

Comments

comments powered by Disqus