Monday, March 31, 2014

wxOpenGL frame work for rendering numpy data



wxOpenGL framework

This frame work lets you render and dynamically modify numpy data (quad-, triangle-meshes) via provided wxOpenGL framework. 

Application functions (keyboard keys):
+ E - Toggle edge visibility

Application functions (viewport manipulation):
+ Middle mouse button & mouse motion - rotate opject
+ Mouse wheel - Zoom
+ Ctrl & Middle mouse button & mouse motion - rotate object
+ Shift & Middle mouse button & mouse motion - translate object

The framework can be downloaded at:
https://www.dropbox.com/sh/ip306uavm7jpekg/2adXHWODU1

The package contains the application (wx_opengl.py) and a render object class (render_object.py) which may be used to add triangle / quad data and meshed to the application viewport. 

Friday, March 28, 2014

Affine Image Registration (follow up, a feature-based approach)

This is a follow up on the previous rigid image registration approach. In the following example I am showing you how to implement an affine transformation registration based on image features with python. In our lab we quite often have to analyze pharmacokinetic parameters based on transient fluorescence imaging taken from tissue in an animal skin flap window. It is the nature of tissues to move during the imaging recording process especially when heat is applied to the tissues. To correct the movement between the individual images, a affine transformation approach is presented to correct the images. 

In a first step, the gradient images are derived from the original images and "landmark" features are detected via threshold criteria (Figure 1).


Figure 1

In a second step, affine transformation is applied to the matching image via optimization method until a defined distance criteria (nearest neighbor criteria with KDTree) between the landmarks reached a minimum (Figure 2).


Figure 2

The Python code:


       
'''
@author: Christian Rossmann, PhD
@license:  Public Domain
@blog: http://scientificcomputingco.blogspot.com/
'''

from scipy import ndimage
from scipy.spatial import cKDTree
from scipy.optimize import fmin

def register(im1,im2):
    
    # initialize the transformation matrix
    T = [[1.0,0.0],[0.0,1.0]]
    
    # extract features from first image  
    edges1 = ndimage.sobel(im1)
    features1 = dstack(where(edges1.T > edges1.mean()))[0]
    kd = cKDTree(features1)
    
    def errfunction(T,kd=kd,im2=im2):
        
        # transform the image
        img = ndimage.interpolation.affine_transform(im2,T.reshape(2,2))
        
        # extract features from transformed image 
        edges = ndimage.sobel(img)
        features = dstack(where(edges.T > edges.mean()))[0]
        
        dst,idx = kd.query([features])
        
        # measure the distances between features.
        return sum(dst**2)
        
    T = fmin(errfunction,T,xtol=1.e-6,ftol=1.e-6,maxfun=5000,maxiter=5000)
    
    return (T,ndimage.interpolation.affine_transform(im2,T.reshape(2,2)))

im1 = imread("sample_images/IMG1.png")[:,:,0]
im2 = imread("sample_images/IMG2.png")[:,:,0]

im1 = ndimage.filters.uniform_filter(im1,5)
im2 = ndimage.filters.uniform_filter(im2,5)

im1 -= im1.min(); im1 /= im1.max() 
im2 -= im2.min(); im2 /= im2.max()

(T,rimg) = register(im1,im2)

#%%
figure(1)
clf()
s = subplot(1,1,1)
title('Gradient Image with Detected Features')
imshow(edges,'bone')
edges = ndimage.sobel(im1)
features = dstack(where( (edges.T > 0.45) & (edges.T < 0.55) ))[0]
plot(features[:,0],features[:,1],'+',color='red')
xlim([0,im1.shape[0]])
ylim([0,im1.shape[1]])
s.axes.get_xaxis().set_visible(False)
s.axes.get_yaxis().set_visible(False)
#%%
figure(2)
clf()

s = subplot(2,2,1)
title('Original Image')
imshow(im1,'bone')
s.axes.get_xaxis().set_visible(False)
s.axes.get_yaxis().set_visible(False)

s = subplot(2,2,2)
title('Matching Image')
imshow(im2,'bone')
s.axes.get_xaxis().set_visible(False)
s.axes.get_yaxis().set_visible(False)

s = subplot(2,2,3)
title('Matching Image subtracted from Original Image')
imshow(im2-im1,'jet')
s.axes.get_xaxis().set_visible(False)
s.axes.get_yaxis().set_visible(False)

s = subplot(2,2,4)
title('Registred Image subtracted from Original Image')
imshow(rimg-im1,'jet')
s.axes.get_xaxis().set_visible(False)
s.axes.get_yaxis().set_visible(False)

show()