Commit c1e35b29 authored by gonciarz's avatar gonciarz
Browse files

Initial commit

#include "mex.h"
#include <math.h>
// finds out if a point is inside a polygon
// to be compiled using mex
void mexFunction( int nlhs, mxArray *plhs[],
int nrhs, const mxArray *prhs[] )
// general remarks: this is a standard mex-function setup. Arguments are passed as
// arrays of pointers to the input and output, respectively (*plhs[] means pointer
// left hand side, thus output, *prhs[] is for input). We use the functions mxGetPr
// and mxGetScalar to get local copies of the data (thzat is at least how I
// understand it)
// pointers to input:
mxArray *p_poly_x, *p_poly_y;
// input scalars:
int n_pol, n_test;
mxArray *p_test_x, *p_test_y;
// pointers to output
double *p_in;
// array input data, all 1D arrays:
double *poly_x, *poly_y, *test_x, *test_y;
// copy input pointers
p_poly_x = prhs[0];
p_poly_y = prhs[1];
p_test_x = prhs[2];
p_test_y = prhs[3];
// get matrices out of mxArrays
poly_x = mxGetPr(p_poly_x);
poly_y = mxGetPr(p_poly_y);
test_x = mxGetPr(p_test_x);
test_y = mxGetPr(p_test_y);
// get scalars
n_pol = (int)mxGetScalar(prhs[4]);
n_test = (int)mxGetScalar(prhs[5]);
// allocate memory and assign output pointer
plhs[0] = mxCreateDoubleMatrix(n_test, 1, mxREAL); //mxReal is our data-type
// get a pointer to the data space in our newly allocated memory
p_in = mxGetPr(plhs[0]);
// check if the points are inside or not
// sorry for the cryptic loop, see
// for background info on the idea of the checking method
int i, j, k, c;
for (k = 0; k < n_test; k++) {
c = 0;
for (i = 0, j = n_pol-1; i < n_pol; j = i++) {
if ((((poly_y[i] <= test_y[k]) && (test_y[k] < poly_y[j])) ||
((poly_y[j] <= test_y[k]) && (test_y[k] < poly_y[i]))) &&
(test_x[k] < (poly_x[j] - poly_x[i]) * (test_y[k] - poly_y[i]) / (poly_y[j] - poly_y[i]) + poly_x[i]))
c = !c;
p_in[k] = (double)c;
\ No newline at end of file
File added
File added
File added
File added
1. Set up the tool
- All code should be in the main folder. Either set a path to the folder or work directly from within the folder.
- There are a few functions coded in c that might have to be recompiled to be used in MATLAB. There are compiled versions for G4 Macs and Intel Macs.
- Some MATLAB toolboxes are used (splines, ...)
2. Play with the tool
Two functions/scripts are most important for the beginning:
- job_specs.m sets all parameters and stores them in the workspace
- estimateShapes.m is the main function. Use it with an image (e.g. "toughSampleImage.tif") with double values in [0,255] as the first argument. Leave the argument S empty (i.e. use "[]"). Use the parameters struct "params" as the "params" argument.
Consider setting the plottingPeriod parameters to low number, then one can see the outlines evolving.
3.) Use the tool properly
- The tool uses two data-files: PSF.mat and Iscaling.mat. They contain a 1D PSF (radial profile) and the intensity calibration information. They are created with the output of estimatePSF.m and intensityCalibrationCurve.m. For each data set they need to be recomputed to account for the specific imaging conditions. The two files given match the image "toughSampleImage.tif".
- The many parameters need to be set properly. The first bunch of parameters is important for the segmentation, the second for selection of objects, the third for the refinement. Parameters are usually explained in the functions where they are used.
Appendix 1: List of cleaned functions
\ No newline at end of file
% ADDPOISSONNOISE: add poisson distributed noise to a microscope image
% INPUT: I: a gray-scale image
% FWC: the full-well capacity, namely the number of
% photons detected corresponding to the max signal
% OUTPUT: I_noisy: the image with added noise
% NOTE: Poisson noise comes from processes with countable events. One
% photon on the detector is one event, one gray-sclae value is not an
% event. Hence one needs to estimate how many photons correspond to a
% gray-scale value to get a good idea of the poisson noise level.
% Jo Helmuth, Mosaic, ETHZ, 11.02.09
function I_noisy = addPoissonNoise(I,FWC)
% re-scale to electrons
electrons = FWC/255 * I;
% add poisson noise
I_noisy = poissrnd(electrons);
% re-scale to intensity
I_noisy = (255/FWC) * I_noisy;
% saturate
Ind_sat = find(I_noisy>=255);
I_noisy(Ind_sat) = 255;
% convert
I_noisy = uint8(round(I_noisy));
\ No newline at end of file
% CALCPEAKSHAPE: measure the radial intensity profile about the intensity
% centroid of a bright object
% INPUT: I: image matrix (double)
% c_init: position of local intensity maximum
% OUTPUT: Ir: radial intensity profile
% r: sampling radii
% peakIntensity: peak intensity of the bright object
% c: intensity centroid of bright object
% For PSF estimation it does not matter, for intensity scaling "I" should be
% background free. I should not be scaled for intensity scaling, for PSF
% estimation it doesn't matter
% The centroid estimation is the improved version of the centroid
% estimation in Ivo's 2005 SPT tracking paper
% A clarification regarding the coordinate system:
% - pixel centers are on integer positions
% - the uppermost and leftmost pixel is at (1,1)
% - the first index in an image matrix is y, the second is x.
% - the (+y)-direction is downward
% Jo Helmuth, Mosaic, ETHZ, 11.02.09
function [Ir,r,peakIntensity,c] = calcPeakShape(I,c_init,params)
% This function is used for two things: PSF estimation and measurement of
% radial intenity profiles for the estimation of the intensity scaling
% parameter. The parameter names are motivated by the first application.
% set parameters
dr = params.drPSFmeasurement; % spacing for the sampling in PIXELS
r_max = params.radiusPSFmeasurement; % maximum sampling radius in PIXELS
pointSize = params.pointSize; % approx. point size in PIXELS
% The following parameters are arbitrary hard-coded choices
phi = linspace(0,18/10*pi,10); % angles for the sampling
r = 0:dr:r_max; % sampling radii
% get gauss kernel for image smoothing, used for centroid estimation
g = fspecial('gauss',[3 3],0.5);
% STEP 0: image preprocessing
% smooth the image
% keep the original for later measurement of radial intensity profile
I_orig = I;
I = imfilter(I,g);
% STEP 1: refine centroid estimate
itCnt = 0;
c = c_init;
% this "20" is a empirically found value that gives good results
while (itCnt < 20)
% integer position (nearest pixel center)
c_int = round(c);
% offset from pixel center
c_off = c - c_int;
itCnt = itCnt + 1;
% ensure that we do not try to sample intensities outside the image
% domain
max_radius = pointSize;
max_radius = min(max_radius,size(I,2) - c_int(2));
max_radius = min(max_radius,c_int(2)-1);
max_radius = min(max_radius,size(I,1) - c_int(1));
max_radius = min(max_radius,c_int(1)-1);
radius = max_radius;
% x and y distance matrices
a = -radius:1:radius;
% Distances between pixel centers
Apc = repmat(a,length(a),1);
Bpc = Apc';
% Distances between current centroid and pixel centers
A = Apc - c_off(2);
B = Bpc - c_off(1);
% compute weights for pixels based on their overlap with a circular
% mask
W = pointSize - sqrt(A.*A + B.*B) + 0.5;
W(W<0) = 0;
W(W>1) = 1;
% take the part of the image around the centroid that is relevant
miniIm = I((c_int(1)-radius):(c_int(1)+radius),(c_int(2)-radius):(c_int(2)+radius));
weightedMiniIm = miniIm.*W;
% total intensity of that image
sumWeightedMiniIm = sum(weightedMiniIm(:));
% compute the centroid offset
dc(1) = sum(sum(B.*weightedMiniIm))/sumWeightedMiniIm;
dc(2) = sum(sum(A.*weightedMiniIm))/sumWeightedMiniIm;
% move to new centroid
c = c + dc;
% STEP 2: measure the intensity around the centroid
% put sample points on concentric circles
R = repmat(r,length(phi),1);
PHI = repmat(phi',1,length(r));
X = c(2) + R.*cos(PHI);
Y = c(1) + R.*sin(PHI);
% sample
I_sampled = zeros(size(X,1),size(X,2));
I_sampled(:) = interp2(I_orig,X(:),Y(:),'spline');
% average on the circles
Ir = mean(I_sampled,1);
% STEP 3: finalize the output
% the peak Intensity at this bright point
peakIntensity = max(Ir);
% do some test whether the profile makes any sense:
% 1.) Is the peak kind of close to zero?
i_max = find(Ir==peakIntensity,1,'first');
i_n = length(r);
ok = i_max < 0.5*i_n;
% 2.) Is the peak non-negative?
ok = ok && peakIntensity==max(abs(Ir));
% 3.) No strong minimum between i_0 and i_max?
ok = ok && isempty(find(Ir(1:i_max)<(0.5*peakIntensity)));
warning('The intensity profile is non-normal. Returning empty output.');
Ir = [];
r = [];
peakIntensity = [];
c = [];
% scale and brush the profile
% this assumes the image has a background that is kind of homogenous. It is
% a bit sensitive to funny things in the background and nearby objects.
%minIr = min(Ir);
% this assumes zero background and is more robust, since too small sampling
% radii do not matter. However, if there is still a background, especially
% PSF estimation will behave weird
minIr = 0;
% scale maximum to 1, and set first minimum and everything beyond to 0
maxIr = max(Ir);
ind = find(Ir == minIr,1,'first');
Ir(ind:end) = minIr;
Ir = (Ir-minIr)/(maxIr-minIr);
% CLUSTEROBJECTS: form groups of neighboring objects based on simple
% minimum distance threshold
% INPUT: Object: Array of object structs
% W: Segmentation
% Im: Gray-scale image on which the segmentation was based
% OUTPUT: ObjectClusters: Array of cluster structs
% This function groups neighboring objects in clusters. This allows to
% later improve the outline estimates for multiple objects simultaneously,
% which is important when objects are so close that they influence
% identical pixels in the image.
% The clusters are formed by adding objects to clusters if they are close
% to any cluster and merging clusters into single clusters when at least
% one object of one cluster is close to one object of the other cluster.
% The distance threshold is an arbritrary choice but reflects the typical
% width of a PSF. The PSF is, however, not explicitely taken into account,
% which should be done at some point. On the other hand, the threshold is
% set very conservative (meaning high).
% NOTE 1: This function runs slow.
% NOTE 2: Although it looks awful and complicated, it work pretty well.
% NOTE 3: Any other algorithm for forming groups would be fine.
% Jo Helmuth, Mosaic, ETHZ, 12.02.09
function ObjectClusters = clusterObjects(Objects,W,Im)
% The distance threshold
% TODO: calculate this based on the width of the actual PSF
distThresh = 13;
% label of background
% TODO: make this a function parameter at some point
backgroundLabel = W(1,1);
% image size
s = size(W);
% find highest label and store labels in a vector for convenience
maxLabel = -1;
for i = 1:length(Objects)
labels(i) = Objects(i).label;
maxLabel = max(labels(i),maxLabel);
% STEP 1: set a cluster label for each object
% allocate memory for the cluster labels
clusterLabels = zeros(length(Objects),1);
% loop through all objects and find potential neighbors and check cluster
% condition
clusterCnt = 0;
for i = 1:length(Objects)
% pick object and define ROI
snake = Objects(i).snake;
% this the upper left corner of the ROI in real image coordinates
ul = Objects(i).ul - 1;
ul(1) = max(ul(1),1);
ul(2) = max(ul(2),1);
% this the lower right
lr = [ul(1)+max(snake(:,1))+2 ul(2)+max(snake(:,2))+2];
lr(1) = min(lr(1),s(1));
lr(2) = min(lr(2),s(2));
% now see what is inside the ROI
ROI = W(ul(1):lr(1),ul(2):lr(2));
objectPixels = find(ROI==labels(i));
ROI(objectPixels) = 0;
ROI(find(ROI==backgroundLabel)) = 0;
% these are the indices of anything that is not zero at this point
otherObjInd = find(ROI~=0);
otherLabels = [];
% check what labels are present
for j = 1:length(otherObjInd)
otherLabels(end+1) = ROI(otherObjInd(j));
% now check distances to all present objects and build clusters
for j = 1:length(otherLabels)
[obPRow,obPCol] = ind2sub(size(ROI),objectPixels);
[otherObRow,otherObCol] = find(ROI==otherLabels(j));
minDist = computeDistance(obPRow,obPCol,otherObRow,otherObCol);
% build clusters
% associate the object to the current cluster
% "i" is the index of the first object, find the index of the
% second involved object
indexOtherObj = find(labels==otherLabels(j));
% check if objects are already part of a cluster
if(clusterLabels(i)==0 && clusterLabels(indexOtherObj)==0)
% none of the two is in a cluster, make a new one
clusterLabels(i) = labels(i);
clusterLabels(indexOtherObj) = labels(i);
elseif(clusterLabels(i)~=0 && clusterLabels(indexOtherObj)==0)
% the first object is in a cluster, assign the second one
% the same cluster label
clusterLabels(indexOtherObj) = clusterLabels(i);
elseif(clusterLabels(i)~=0 && clusterLabels(indexOtherObj)~=0)
% both are in a cluster, assign all objects in the cluster
% of the second object the cluster label of the first object
c2 = clusterLabels(indexOtherObj);
c1 = clusterLabels(i);
clusterLabels(find(clusterLabels==c2)) = c1;
% now the clusterLabels vector has an entry for each object that is in a
% cluster with more than one object
% STEP 2: set cluster label for isolated objects to object label
% fill it to have a reasonable entry for all objects
for i = 1:length(Objects)
clusterLabels(i) = labels(i);
% STEP 3: construct a cluster struct for each cluster of object
% loop over all clusters
for i = min(clusterLabels):max(clusterLabels)
inClusterObj = find(clusterLabels==i);
nObj = length(inClusterObj);
% check dimension of image that contains all objects
ulMin = [10000000 10000000];
lrMax = [0 0];
for j = 1:nObj
% check dimension of image that contains all objects
ul = Objects(inClusterObj(j)).ul;
lr = Objects(inClusterObj(j)).lr;
ulMin(1) = min(ulMin(1),ul(1));
ulMin(2) = min(ulMin(2),ul(2));
lrMax(1) = max(lrMax(1),lr(1));
lrMax(2) = max(lrMax(2),lr(2));
Obs = cell(nObj,1);
for j = 1:nObj
object = Objects(inClusterObj(j));
ul = object.ul;
snake = object.snake;
% move to new coordinate system of the bigger image
snake(:,1) = snake(:,1) + (ul(1)-ulMin(1));
snake(:,2) = snake(:,2) + (ul(2)-ulMin(2));
object.snake = snake;
Obs{j} = object;
% get the bigger image that contains all objects in the cluster
miniIm = Im(ulMin(1):lrMax(1),ulMin(2):lrMax(2));
miniW = W(ulMin(1):lrMax(1),ulMin(2):lrMax(2));
% make struct with the cluster data
newCluster = struct('objects',{Obs},'miniIm',miniIm,'miniW', ...
% append to the array of clusters
clusterCnt = clusterCnt + 1;
ObjectClusters(clusterCnt) = newCluster;
function minDist = computeDistance(obPRow,obPCol,otherObRow,otherObCol)
minDist = 1000;
for i = 1:length(obPRow)
for j = 1:length(otherObRow)
d_sq = (obPRow(i)-otherObRow(j))^2 + (obPCol(i)-otherObCol(j))^2;
minDist = min(minDist,d_sq);
#include "mex.h"
#include <math.h>
// computes distances of points to a polygon
// to be compiled using mex
/* The function does not require anything from the polygon, except that it
must be closed, thus the first point is repeated as the last point.
The function computes the distance to the polygon points, takes the
shortest as candidate and than checks if any point-line distance is
shorter. If it is, it is further checked if the intersection of the line
and a line normal to it through the point lies between the two supporting
polygon points of the line */
/* this function is usually NOT called directly, rather dist2poly.m is used,
which also checks the sign of the distance to account for inside resp.
outside points*/
void mexFunction( int nlhs, mxArray *plhs[],
int nrhs, const mxArray *prhs[] )
// pointers to input
mxArray *p_poly_x, *p_poly_y, *p_normal_x, *p_normal_y,*p_q_x, *p_q_y;
// input scalars
int n_p, n_q;
// pointers to output
double *d;
// input data
double *p_x, *p_y, *q_x, *q_y, *n_x, *n_y;
// get all the input data etc.
// copy input pointers
p_poly_x = prhs[0];
p_poly_y = prhs[1];
p_normal_x = prhs[2];
p_normal_y = prhs[3];
p_q_x = prhs[4];
p_q_y = prhs[5];
// get matrices out of mxArrays
p_x = mxGetPr(p_poly_x);
p_y = mxGetPr(p_poly_y);
n_x = mxGetPr(p_normal_x);
n_y = mxGetPr(p_normal_y);
q_x = mxGetPr(p_q_x);
q_y = mxGetPr(p_q_y);
// get scalars
n_p = (int)mxGetScalar(prhs[6]);
n_q = (int)mxGetScalar(prhs[7]);
// organize stuff for output
// allocate memory and assign output pointer
plhs[0] = mxCreateDoubleMatrix(n_q, 1, mxREAL); //mxReal is our data-type
// get a pointer to the data space in our newly allocated memory
d = mxGetPr(plhs[0]);
// space for internal variables
double *p2p_x, *p2p_y, *absp2p, *q2p_x, *q2p_y;
p2p_x = (double*)malloc(sizeof(double)*n_p-1);
p2p_y = (double*)malloc(sizeof(double)*n_p-1);
absp2p = (double*)malloc(sizeof(double)*n_p-1);
q2p_x = (double*)malloc(sizeof(double)*n_p-1);
q2p_y = (double*)malloc(sizeof(double)*n_p-1);
int ip, iq, i;
double d_min;
double d_temp, proj, absqp;
// do the actual work
// compute point to point distances within polygon