Measuring weighted mean length from an electrophoresis gel image

Question

Background:

My question relates to extracting feature from an electrophoresis gel (see below). In this gel, DNA is loaded from the top and allowed to migrate under a voltage gradient. The gel has sieves so smaller molecules migrate further than longer molecules resulting in the separation of DNA by length. So higher up the molecule, the longer it is.

Question:

In this image there are 9 lanes each with separate source of DNA. I am interested in measuring the mean location (value on the y axis) of each lane. I am really new to image processing, but I do know MATLAB and I can get by with R with some difficulty. I would really appreciate it if someone can show me how to go about finding the mean of each lane.

Answer 1

Here's my try. It requires that the gels are nice (i.e. straight lanes and the gel should not be rotated), but should otherwise work fairly generically. Note that there are two image-size-dependent parameters that will need to be adjusted to make this work on images of different size.

%# first size-dependent parameter: should be about 1/4th-1/5th
%# of the lane width in pixels.
minFilterWidth = 10;

%# second size-dependent parameter for filtering the 
%# lane profiles
gaussWidth = 5;

%# read the image, normalize to 0...1
img = imread('http://img823.imageshack.us/img823/588/gele.png');
img = rgb2gray(img);
img = double(img)/255;

%# Otsu thresholding to (roughly) find lanes
thMsk = img < graythresh(img);

%# count the mask-pixels in each columns. Due to 
%# lane separation, there will be fewer pixels
%# between lanes
cts = sum(thMsk,1);

%# widen the local minima, so that we get a nice
%# separation between lanes
ctsEroded = imerode(cts,ones(1,minFilterWidth));

%# use imregionalmin to identify the separation 
%# between lanes. Invert to get a positive mask
laneMsk = ~repmat(imregionalmin(ctsEroded),size(img,1),1);

Image with lanes that will be used for analysis

%# for each lane, create an averaged profile
lblMsk = bwlabel(laneMsk);
nLanes = max(lblMsk(:));

profiles = zeros(size(img,1),nLanes);
midLane = zeros(1,nLanes);

for i = 1:nLanes
profiles(:,i) = mean(img.*(lblMsk==i),2);
midLane(:,i) = mean(find(lblMsk(1,:)==i));
end

%# Gauss-filter the profiles (each column is an
%# averaged intensity profile
G = exp(-(-gaussWidth*5:gaussWidth*5).^2/(2*gaussWidth^2));
G=G./sum(G);
profiles = imfilter(profiles,G','replicate'); %'

%# find the minima
[~,idx] = min(profiles,[],1);

%# plot
figure,imshow(img,[])
hold on, plot(midLane,idx,'.r')