Which structuring elements to detect branched points?

In a 2D skeleton obtained after thinning of a binary image, pixels having 3 or 4 (or more...) neighbors are branched-points. Hit-or-miss operator can be used to detect branched points (or junction points). Structuring elements, SE, can be 3x3 arrays representing the pixels configuration to be detected. Here different 3x3 structuring elements were explored on some simple test images with the hit-or-miss operator implemented in mahotas 1.0 .

• There are two SE to detect pixels with 4 neighbors, called here X0 and X1.
• There are two set of SE detecting pixels with at least 3 pixels, a Y-like configuration and a T-like configuration. Each of these configurations can be rotated by 45°, making 16 different SE.

18 SE were found , since don't care pixels were used a branched points can be detected by different SE, a pixel detected by X0 should be detected by T4 too. The 18 SE can be displayed as 3x3 image:

Black:background pixel (0), Blue: foreground pixel (1), Red:Don't care pixel (0 or 1)

Applied on a small skeleton-like test image (left), two branched-points  (junction) and four end-points are detected:

On a larger image, from the SIID images database:

Hit-or-miss transformation applied on the skeleton (up right), yields end-points (lower right). Since the same branched-points can be detected by different structuring elements, branched points appear of a variable color according to the detection occurence. Moreoever, branched-points can be connected yielding a branched domain of pixels.

ipython notebook

HMT-SE.ipynb

Detecting end points or branched points from the skeleton of a single chromosome

The end-points or the branched-points from the skeleton of a single chromosome can be extracted using mathematical morphology operators implemented in mahotas.

Two images of metaphasic human chromosomes hybridized with a Cy3-telomeric PNA probe and  counterstained with DAPI, were used.

 

False color DAPI image

False color telomeres image

The location of a chromosome was identified, some preprocessing step were performed to segment the chromosome (hi-pass filtering followed by a simple thresholding).

The end-points, the branching points are extracted from the skeleton with a hit-or-miss operator implemented in mahotas, then overlayed on the original chromosome image.

The contour of the segmented telomeres are also overlayed:

Depending on the amount of high-pass filtering the chromosome skeleton fit the chromosome axis or the chromatids:

Skeletonization after "low" high-pass filtering

Skeletonization after a stronger high-pass filtering

In both cases, the skeleton is displayed in red, the branched-points in green, the end-points in blue, the chromosome contour in magenta, and finally the telomeres in yellow.
To run the code (python 2.7), the path to the images must be modified and mahotas and pymorph must be installed.
Download code

Closed curve to set of ordered pixels: a python test

To transform a closed digital curve, remove one pixel, remember it and submit the opened curve to the following python code.
For example, the minimal closed curve:

gives:

closed curve
[[0 0 0 0 0]
 [0 0 1 0 0]
 [0 1 0 1 0]
 [0 1 0 1 0]
 [0 0 1 0 0]
 [0 0 0 0 0]]

index list
[array([1, 2]), 
array([2, 1]), 
array([3, 1]), 
array([4, 2]), 
array([3, 3]), 
array([2, 3]), 
array([1, 2])]
The fourth pixel:
[4 2]

Using functions, the code is a little bit cleaner, it relies only on numpy and scipy.ndimage for mathematical morphology (hit or miss), since only endpoints are needed here.If branched points are required, I will continue to use mahotas implementation of hit or miss operator because it handles the "don't care" pixels written with a "2" in the structuring element. Thus some of my scripts may not be run on windows, because mahotas has some installation problems due to freeimage on that plateform. I don't know if the scipy implementation of the hit or miss operator handle the "don't care" pixels.

# -*- coding: utf-8 -*-
"""
Created on Mon Jun 20 09:08:38 2011

@author: Jean-Patrick
"""
import time
import numpy as np
from scipy import ndimage as nd

def opencurveToPixelsList(im):
    def scipyEndPoints(imcurve):
        se1=np.array([[0,0,0],[0,1,0],[0,1,0]])
        se2=np.array([[0,0,0],[0,1,0],[0,0,1]])
        EndPoints=np.zeros(imcurve.shape,dtype=np.uint)
        #perform hit & miss 
        for i in (0,3):
            hit1=nd.morphology.binary_hit_or_miss(imcurve,se1)
            se1=np.rot90(se1)
            EndPoints=EndPoints+hit1
            hit2=nd.morphology.binary_hit_or_miss(imcurve,se2)
            se2=np.rot90(se1)
            EndPoints=EndPoints+hit2
        return EndPoints
    #total number of pixels in the curve
    pixel=np.array([-1,-1])
    pixN=np.sum(im==1)
    pixelsList=[]
    #get end point by H&Miss op by ndimage
    #print scipyEndPoints(im)
    #get end points by hit&miss operator provided by mahotas
    #ep=bep.get_endpoints(im)
    ep=scipyEndPoints(im)
    lab,n=nd.label(ep)
    #where the first endpoint is
    first_indices=np.where(lab==1)
    pixel[0]=first_indices[0][0]
    pixel[1]=first_indices[1][0]
    pixelsList.append(np.copy(pixel))
    #print "first point",first_indices," vector",walkingPixel
    firstEndPoint=np.uint8(lab==1)
    #keep an image of the second end point
    lastEndPoint=np.uint8(lab==2)
    last_Indices=np.where(lastEndPoint==1)
    #print "last point",last_Indices
    ######################################################
    ## walk on curve with a 3x3 neighborhood
    ######################################################
    #Init
    current_curve=np.copy(im)
    current_point=np.copy(firstEndPoint)
    ###start to walk on curve
    #remove the second last point
    #current_curve=current_curve-lastEndPoint##too heavy?, try indices
    current_curve[last_Indices[0][0],last_Indices[1][0]]=0
    c_point_ind=np.where(current_point==1)
    #print c_point_ind
    li=c_point_ind[0][0]
    col=c_point_ind[1][0]
    for i in range (0,pixN-2):    
        #3x3 neighborhood arround the endpoint
        neighbor=current_curve[li-1:li+2,col-1:col+2]
        neighbor[1,1]=0
        #################
        #Can only handle a curve
        #such np.where(neihgbor==1) must gives only one pixel
        #############
        nextPointIndices=np.where(neighbor==1)##vectO'M=nextPointIndices
        ##remove the first point from the curve
        current_curve[li,col]=0
        pixN=pixN-1
        #print current_curve
        ##compute nextPoint indices in the original base vectOM=OO'+O'M
        li=(li-1)+nextPointIndices[0][0]
        col=(col-1)+nextPointIndices[1][0]
        pixel[0]=li
        pixel[1]=col
        pixelsList.append(np.copy(pixel))
    #don't forget the last pixel
    pixel[0]=last_Indices[0][0]
    pixel[1]=last_Indices[1][0]
    pixelsList.append(np.copy(pixel))
    return pixelsList
############################
def closedcurveToPixelsList(im):
    '''
    given a closed curve (no test),
    not touching the image border(no test) as:
        [0,0,0,0,0],
        [0,0,1,0,0],
        [0,1,0,1,0],
        [0,1,0,1,0],
        [0,0,1,0,0],
        [0,0,0,0,0]
    The function returns an array list of ordered 
    points along the curve:
        [1,2],[2,3],[3,3],[4,2],[3,1],[2,1],[1,2]
    '''
    firstpoint=np.copy(np.where(im==1))
    print np.where(im==1)[1][0]
    print "first point",firstpoint
    openedcurve=np.copy(im)
    openedcurve[firstpoint[0][0],firstpoint[1][0]]=0
    print "now it's open"
    print openedcurve
    openlist=opencurveToPixelsList(openedcurve)
    #add the first point of the closed curved at the 
    #begiining and the end of the list:
    pixel=np.array([-1,-1])
    pixel[0]=firstpoint[0][0]
    pixel[1]=firstpoint[1][0]
    openlist.append(np.copy(pixel))
    openlist.insert(0,np.copy(pixel))
    return openlist
#tests
im=np.array([[0,0,0,0,0,0,0],
            [0,1,0,0,1,0,0],
            [0,0,1,0,0,1,0],
            [0,0,0,1,1,0,0],
            [0,0,0,0,0,0,0]])
liste=opencurveToPixelsList(im)
closedcurve=np.array([
        [0,0,0,0,0],
        [0,0,1,0,0],
        [0,1,0,1,0],
        [0,1,0,1,0],
        [0,0,1,0,0],
        [0,0,0,0,0]])
liste2=closedcurveToPixelsList(closedcurve)
print liste
print liste[0]
print "closed curve"
print closedcurve
print liste2
print liste2[3]

Converting an open curve to a list of ordered pixels: a python test code

In some conditions, a particle skeleton is an open curve with no branched point, according to a 4-connected neighborhood:

Test curve

Getting the pixels indices according to their order on the curve may be usefull (computing direction, bending,fitting...).

The following SkelToPixel_2 script, depends on two other scripts:BranchedEndPoints and StructElem. The scripts must be stored in the same folder.

The idea is simple, the script takes the curve at one end, walk on it searching a neighbor pixel up to the ...antepenultimate pixel (that works).Opencv yields the coordinates of pixels belonging to the contour of a particle, i.e. a closed curve. I don't know how to, but I suppose it could be possible to do the same thing with opencv for an open curve.

# -*- coding: utf-8 -*-
"""
Created on Tue Jun  7 10:35:17 2011
                      SkelToPixel_2
A code to convert an open curve stored in a numpy array
into a list of neihbor pixels

@author: jean-Patrick Pommier
http://dip4fish.blogspot.com
"""

import numpy as np
from scipy import ndimage as nd
import BranchedEndPoints as bep

im=np.array([[0,0,0,0,0,0,0],
            [0,1,0,0,0,0,0],
            [0,0,1,1,0,0,0],
            [0,0,0,0,1,1,0],
            [0,0,0,0,0,0,0]]) 
#total number of pixels in the curve
pixN=np.sum(im==1)
pixelsList=[]
#get end points by hit&miss operator
ep=bep.get_endpoints(im)
lab,n=nd.label(ep)
#where the first endpoint is
first_indices=np.where(lab==1)
firstEndPoint=np.uint8(lab==1)
#keep an image of the second end point
lastEndPoint=np.uint8(lab==2)
pixelsList.append(first_indices)
current_curve=np.copy(im)
current_point=firstEndPoint
print current_curve
#start to walk on the curve
while pixN<>2:
    current_curve=current_curve-current_point
    pixN=pixN-1
    #this is not very smart, using a 3x3 neighborhood
    #arround the current point would reduce the research
    current_point=bep.get_endpoints(current_curve)-lastEndPoint
    pixelsList.append(np.where(current_point==1))
#add the last pixel
pixelsList.append(np.where(lastEndPoint==1))
for i in range(0,len(pixelsList)):
    print "i=",i,pixelsList[i]

Run from spyder, the script output is:

[[0 0 0 0 0 0 0]
 [0 1 0 0 0 0 0]
 [0 0 1 1 0 0 0]
 [0 0 0 0 1 1 0]
 [0 0 0 0 0 0 0]]
i= 0 (array([1]), array([1]))
i= 1 (array([2]), array([2]))
i= 2 (array([2]), array([3]))
i= 3 (array([3]), array([4]))
i= 4 (array([3]), array([5]))

The script doesn't yield directly pixels indices.
Feel free to propose a nicer, faster code.

The result with another curve

And the output console:

[[0 0 0 0 0 0 0]
 [0 1 0 0 1 0 0]
 [0 0 1 0 0 1 0]
 [0 0 0 1 1 0 0]
 [0 0 0 0 0 0 0]]
i= 0 (array([1]), array([1]))
i= 1 (array([2]), array([2]))
i= 2 (array([3]), array([3]))
i= 3 (array([3]), array([4]))
i= 4 (array([2]), array([5]))
i= 5 (array([1]), array([4]))

Detecting end points in a skeleton

In order to detect overlapping chromosomes, or to analyse the shape of a chromosome it could be useful to analyse the shape of the particle skeleton.

In the previous post hit-or-miss morphological operator, with the proper structuring element, detects branched points of a skeleton. Those points are indicated to classify overlapping chromosomes, even if pruning the skeleton might be necessary.

Here, with other structuring elements, hit-or-miss operator detects the end points of particle constituted of overlapping and touching chromosomes.

Download python code

 Lower image: green points (branched points of the skeleton), blue points (end points of the skeleton)