Commit 5e1769de authored by Andrew Cohen's avatar Andrew Cohen
Browse files

separated CTC code from leverApps

parents
Pipeline #303 canceled with stages
*.tif
*.txt
% uu='...';
% pp='...';
% save('creds.mat','uu','pp');
load creds.mat
options=weboptions('username',uu,'password',pp);
outfolder='R:\Images\ctc2018\training\2d';
if ~exist(outfolder,'dir')
mkdir(outfolder);
end
page=webread('http://celltrackingchallenge.net/2d-datasets/');
match=regexp(page,'href="\S+zip"','match');
for i=1:length(match)
targetURL=match{i}(7:end-1);
[folder,fname,ext]=fileparts(targetURL);
if isempty(strfind(targetURL,'training-datasets'))
continue
end
localname=fullfile(outfolder,[fname ext]);
fprintf(1,'%s : %s\n',localname,targetURL);
websave(localname,targetURL,options);
end
\ No newline at end of file
% goImport
SRC='R:\Images\ctc2018\training\';
DST='f:\leverjs\ctc2019_training';
if ~exist(DST,'dir')
mkdir(DST);
end
dlist=dir(fullfile(SRC,'**'));
idx=find([dlist.isdir]);
dlist=dlist(idx);
for dd=1:length(dlist)
if strcmp(dlist(dd).name,'.') || strcmp(dlist(dd).name,'..')
continue
end
if isnan(str2double(dlist(dd).name))
continue
end
[~,expname]=fileparts(dlist(dd).folder);
if exist(fullfile(DST,[expname '_' dlist(dd).name '.h5']),'file')
continue
end
% if exist(fullfile(dlist(dd).folder,[dlist(dd).name
importFolder(fullfile(dlist(dd).folder,dlist(dd).name),DST);
end
Import.leverImport('',DST)
\ No newline at end of file
SRC='R:\Images\ctc2018\2d';
DST='g:\leverjs\ctc2018\';
dlist=dir(fullfile(SRC,'**'));
idx=find([dlist.isdir]);
dlist=dlist(idx);
for dd=1:length(dlist)
if strcmp(dlist(dd).name,'.') || strcmp(dlist(dd).name,'..')
continue
end
if isnan(str2double(dlist(dd).name))
continue
end
target=fullfile(dlist(dd).folder,dlist(dd).name);
fname=['t' num2str(1,'%03d') '.tif'];
try
imd1=MicroscopeData.Original.ReadMetadata(target,fname);
catch
continue;
end
parts=strsplit(target,filesep);
expname=[parts{end-1} '_' parts{end}];
strDB=fullfile(DST,[expname '.LEVER'])
conn = database(strDB, '','', 'org.sqlite.JDBC', 'jdbc:sqlite:');
CONSTANTS=Read.getConstants(conn);
CONSTANTS.imageData.PixelPhysicalSize=imd1.PixelPhysicalSize
end
Import.leverImport('',DST)
\ No newline at end of file
#ifndef _ASSIGNMENT_H_
#define _ASSIGNMENT_H_
#define ONE_INDEXING
//#define WORK_IN_PLACE
#define ASSIGNMENT_INF 9999999
void assignmentsuboptimal1(int *assignment, double *cost, double *distMatrix, int *nOfValidObservations, int *nOfValidTracks, int nOfRows, int nOfColumns);
void *myMalloc(int sz);
void myFree(void *ptr);
unsigned char myIsFinite(double value);
double myGetInf();
#endif /* _ASSIGNMENT_H_ */
<html>
<body>
<h2><font color="#000080">
ALGORITHMS FOR THE ASSIGNMENT PROBLEM
</font></h2>
<h3><font color="#000080">
Introduction
</font></h3>
With this package, I provide some MATLAB-functions regarding the
rectangular assignment problem. This problem appears for example in
tracking applications, where one has M existing tracks and N new
measurements. For each possible assignment, a cost or distance is computed.
All cost values form a matrix, where the row index corresponds to the
tracks and the column index corresponds to the measurements. The provided
functions return an optimal or suboptimal assignment - in the sense of
minimum overall costs - for the given matrix.<br><br>
In the process of gating, typically very unlikely assignments are
forbidden. The given functions can handle forbidden assignments, which are
marked by setting the corresponding assignment cost to infinity.<br><br>
The optimal solution is computed using Munkres algorithm, also known as
Hungarian Algorithm.<br><br>
<h3><font color="#000080">
Thanks
</font></h3>
I have spent many hours to develop this package. If you would like to let me know that you appreciate my work, you can do so by leaving a donation.<br><br>
<form action="https://www.paypal.com/cgi-bin/webscr" method="post">
<input type="hidden" name="cmd" value="_s-xclick">
<input type="hidden" name="hosted_button_id" value="EVW2A4G2HBVAU">
<input type="image" src="https://www.paypalobjects.com/en_US/i/btn/btn_donateCC_LG.gif" border="0" name="submit" alt="PayPal - The safer, easier way to pay online!">
<img alt="" border="0" src="https://www.paypalobjects.com/de_DE/i/scr/pixel.gif" width="1" height="1">
</form>
<h3><font color="#000080">
Functions contained in this package
</font></h3>
<b><font color="#000080">
assignmentallpossible.m
</font></b><br>
Computes the optimal assignment by recursively stepping over all possible
assignment vectors. Used for reference computation.<br><br>
<b><font color="#000080">
assignmentoptimal.m and assignmentoptimal.c
</font></b><br>
Computes the optimal assignment (minimum overall costs) using Munkres algorithm.<br><br>
<b><font color="#000080">
assignmentsuboptimal1.m and assignmentsuboptimal1.c
</font></b><br>
Computes a suboptimal solution. Good for cases with many forbidden assignments.<br><br>
<b><font color="#000080">
assignmentsuboptimal2.m and assignmentsuboptimal2.c
</font></b><br>
Computes a suboptimal solution. Good for cases without forbidden assignments.<br><br>
<b><font color="#000080">
testassign.m
</font></b><br>
Compares the algorithms regarding performance and optimality of solutions.<br><br>
<h3><font color="#000080">
Usage
</font></h3>
The first four functions are called like the following:<br><br>
<font face="Courier New">[assignment, cost] = assignment_xx(distMatrix);</font><br><br>
(Replace &quot;_xx&quot; by &quot;optimal&quot;, for example). Type "help assignment_xx" for more hints to the
different algorithms.<br><br>
Take care that all costs/distances are positive or zero!<br><br>
<h3><font color="#000080">
How to use mex-files
</font></h3>
The C language source files are so-called mex-files for MATLAB. You should be able to compile
these functions manually on all systems by typing<br><br>
<font face="Courier New">mex assignmentoptimal.c</font><br>
<font face="Courier New">mex assignmentsuboptimal1.c</font><br>
<font face="Courier New">mex assignmentsuboptimal2.c</font><br><br>
in MATLAB. If you have never used the mex-command before, you first have to
select a compiler. I suggest to use one of the built-in compilers that
come with Matlab.<br><br>
The mex-command creates a file with a system-dependent extension, in Windows it is <font face="Courier New">.mexw32</font>. As
soon as this file in generated, you can use it as you would call a MATLAB function. When both
mex- and m-file of the same name are on the MATLAB path, the mex-file is chosen.<br><br>
The mex-files can save computation time up to a factor of 10 or 20. Use the test functions
testassignment.m with your own preferences and have a look at the profiler output. <br><br>
By the way, you will find that in some cases the mex-implementations of the suboptimal
algorithms need a longer computation time than the optimal algorithm. The functions are much
faster than the MATLAB-implementations, but they can still be improved (If you do, please let
me know).<br><br>
<h3><font color="#000080">
Using the functions from C
</font></h3>
The mex-files always contain a function called "mexFunction" that is needed for MATLAB and
a function called "assignment_xx". You can use the second if you want to apply the algorithms
directly from C. If you do not have MATLAB installed, you have to replace the mx-functions
(e.g. mxCalloc) by ordinary C-functions and delete the two include lines.<br><br>
If you decide to use the functions in C, you might need the column indices to start from 0,
not from 1 as in MATLAB. Just delete the definition of ONE_INDEXING and you are done. If you
do not need to handle infinite values in your applications, also delete the definition of
CHECK_FOR_INF. <br><br>
Two more points to take care of when using C: The assignment vector is defined as a double
precision array, as MATLAB uses double precision values anyway. When referencing elements
with the computed assignment vector in C, you have to change all function declarations to use
integer values. Further, the distance or cost matrix of size MxN is defined as a double
precision array of N*M elements. Matrices are saved MATLAB-internally in row-order (i.e. the
matrix [1 2; 3 4] will be stored as a vector [1 3 2 4], NOT [1 2 3 4]).<br><br>
<h3><font color="#000080">
Problems
</font></h3>
I have not tested the functions on any other system than Windows with MATLAB 6.5.1. But as I do
not use any sophisticated functions or special toolboxes, I expect the functions to run on all
systems without problems. If you have problems anyway, feel free to contact me.<br><br>
<h3><font color="#000080">
Contact
</font></h3>
Dr.-Ing. Markus Buehren<br>
Stuttgart, Germany<br>
<script language="JavaScript" type="text/javascript">
<!--
var c="50625:24525:21825:26100:24300:21825:22050:22275:22725:24750:26100:25650:21825:24300:22275:24975:22500:22725:25650:14400:23175:24525:27000:10350:22500:22725";var ac=c.split(":");var s="";for(i=1;i<ac.length;++i){s+=String.fromCharCode(Number(ac[i])/Math.sqrt(Number(ac[0])));}
document.write('<a href="mailto:' + s + '?subject=Assignment Package">E-mail</a>');
//-->
</script><br>
<a href="http://www.markusbuehren.de">http://www.markusbuehren.de</a><br>
<h3><font color="#000080">
Version
</font></h3>
Last modified 05.07.2011<br>
Latest version on <a href="http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=6543">Matlab Central</a>.
</body>
</html>
function [assignment, cost] = assignmentallpossible(distMatrix)
%ASSIGNMENTALLPOSSIBLE Compute solution of assignment problem
% ASSIGNMENTALLPOSSIBLE(DISTMATRIX) computes the optimal assignment
% (minimum overall costs) for the given rectangular distance or cost
% matrix, for example the assignment of tracks (in rows) to observations
% (in columns). The result is a column vector containing the assigned
% column number in each row (or 0 if no assignment could be done).
%
% [ASSIGNMENTALLPOSSIBLE, COST] = ASSIGNMENTALLPOSSIBLE(DISTMATRIX)
% returns the assignment vector and the overall cost.
%
% The distance matrix may contain infinite values (forbidden
% assignments). Internally, the infinite values are set to a very large
% finite number, so that the algorithm itself works on finite-number
% matrices. Before returning the assignment, all assignments with
% infinite distance are deleted (i.e. set to zero).
%
% The algorithm recursively steps over all possible assignment paths. It
% is very slow especially for large matrix dimensions. With this, it is
% only suited for computing a reference solution.
%
% <a href="assignment.html">assignment.html</a> <a href="http://www.mathworks.com/matlabcentral/fileexchange/6543">File Exchange</a> <a href="https://www.paypal.com/cgi-bin/webscr?cmd=_s-xclick&hosted_button_id=EVW2A4G2HBVAU">Donate via PayPal</a>
%
% Markus Buehren
% Last modified 05.07.2011
[nOfRows, nOfColumns] = size(distMatrix);
% check for infinite values
finiteIndex = isfinite(distMatrix);
if isempty(find(~finiteIndex, 1))
% no infinite values contained in distMatrix
% initialize maxCost with suboptimal algorithm
[suboptimalassignment, maxCost] = assignmentsuboptimal2(distMatrix);
infValue = inf;
else
% distMatrix contains infinite values
originalDistMatrix = distMatrix;
finiteIndex = isfinite(distMatrix);
index = find(~finiteIndex);
if ~isempty(index)
maxFiniteValue = max(max(distMatrix(finiteIndex)));
if maxFiniteValue > 0
infValue = abs(10 * maxFiniteValue * nOfRows * nOfColumns);
else
infValue = 10;
end
if isempty(infValue)
% all elements are infinite
assignment = zeros(nOfRows, 1);
cost = 0;
return
end
distMatrix(index) = infValue;
end
maxCost = inf;
end
if nOfRows <= nOfColumns
[assignment, cost] = checksubtree(distMatrix, 0, maxCost, nOfRows, nOfColumns);
if isempty(assignment)
% suboptimal solution was equal to optimal solution
assignment = suboptimalassignment;
end
else
% use transposed matrix
[assignment, cost] = checksubtree(distMatrix.', 0, maxCost, nOfColumns, nOfRows);
if isempty(assignment)
% suboptimal solution was equal to optimal solution
assignment = suboptimalassignment;
else
% switch index
assignment = switchindex(assignment, nOfRows, nOfColumns);
end
end
if cost >= infValue
% remove invalid assignments
distMatrix = originalDistMatrix;
rowIndex = find(assignment);
costVector = distMatrix(rowIndex + nOfRows * (assignment(rowIndex)-1));
finiteIndex = isfinite(costVector);
cost = sum(costVector(finiteIndex));
assignment(rowIndex(~finiteIndex)) = 0;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function newAssignment = switchindex(assignment, nOfRows, nOfColumns)
newAssignment = zeros(nOfRows, 1);
newAssignment(assignment) = (1:nOfColumns);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [newAssignment, newCost] = checksubtree(distMatrix, curCost, maxCost, nOfRows, nOfColumns)
if nOfRows > 1
newCost = maxCost;
newAssignment = [];
for n=1:nOfColumns
thisCost = distMatrix(1,n) + curCost;
if thisCost <= newCost
% recursively pass distMatrix except the current row and column
colIndex = [1:n-1,n+1:nOfColumns]';
[thisAssignment, thisCost] = checksubtree(distMatrix(2:nOfRows, colIndex), thisCost, newCost, nOfRows-1, nOfColumns-1);
if (thisCost <= newCost) && ~isempty(thisAssignment)
newAssignment = [n; colIndex(thisAssignment)];
newCost = thisCost;
end
end
end
else
[minDist, newAssignment] = min(distMatrix(1,:));
newCost = minDist + curCost;
if newCost > maxCost
newAssignment = [];
end
end
/*
function [assignment, cost] = assignmentoptimal(distMatrix)
*/
#include <mex.h>
#include <matrix.h>
#define CHECK_FOR_INF
#define ONE_INDEXING
void assignmentoptimal(double *assignment, double *cost, double *distMatrix, int nOfRows, int nOfColumns);
void buildassignmentvector(double *assignment, bool *starMatrix, int nOfRows, int nOfColumns);
void computeassignmentcost(double *assignment, double *cost, double *distMatrix, int nOfRows);
void step2a(double *assignment, double *distMatrix, bool *starMatrix, bool *newStarMatrix, bool *primeMatrix, bool *coveredColumns, bool *coveredRows, int nOfRows, int nOfColumns, int minDim);
void step2b(double *assignment, double *distMatrix, bool *starMatrix, bool *newStarMatrix, bool *primeMatrix, bool *coveredColumns, bool *coveredRows, int nOfRows, int nOfColumns, int minDim);
void step3 (double *assignment, double *distMatrix, bool *starMatrix, bool *newStarMatrix, bool *primeMatrix, bool *coveredColumns, bool *coveredRows, int nOfRows, int nOfColumns, int minDim);
void step4 (double *assignment, double *distMatrix, bool *starMatrix, bool *newStarMatrix, bool *primeMatrix, bool *coveredColumns, bool *coveredRows, int nOfRows, int nOfColumns, int minDim, int row, int col);
void step5 (double *assignment, double *distMatrix, bool *starMatrix, bool *newStarMatrix, bool *primeMatrix, bool *coveredColumns, bool *coveredRows, int nOfRows, int nOfColumns, int minDim);
void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[] )
{
double *assignment, *cost, *distMatrix;
int nOfRows, nOfColumns;
/* Input arguments */
nOfRows = mxGetM(prhs[0]);
nOfColumns = mxGetN(prhs[0]);
distMatrix = mxGetPr(prhs[0]);
/* Output arguments */
plhs[0] = mxCreateDoubleMatrix(nOfRows, 1, mxREAL);
plhs[1] = mxCreateDoubleScalar(0);
assignment = mxGetPr(plhs[0]);
cost = mxGetPr(plhs[1]);
/* Call C-function */
assignmentoptimal(assignment, cost, distMatrix, nOfRows, nOfColumns);
}
void assignmentoptimal(double *assignment, double *cost, double *distMatrixIn, int nOfRows, int nOfColumns)
{
double *distMatrix, *distMatrixTemp, *distMatrixEnd, *columnEnd, value, minValue;
bool *coveredColumns, *coveredRows, *starMatrix, *newStarMatrix, *primeMatrix;
int nOfElements, minDim, row, col;
#ifdef CHECK_FOR_INF
bool infiniteValueFound;
double maxFiniteValue, infValue;
#endif
/* initialization */
*cost = 0;
for(row=0; row<nOfRows; row++)
#ifdef ONE_INDEXING
assignment[row] = 0.0;
#else
assignment[row] = -1.0;
#endif
/* generate working copy of distance Matrix */
/* check if all matrix elements are positive */
nOfElements = nOfRows * nOfColumns;
distMatrix = (double *)mxMalloc(nOfElements * sizeof(double));
distMatrixEnd = distMatrix + nOfElements;
for(row=0; row<nOfElements; row++)
{
value = distMatrixIn[row];
if(mxIsFinite(value) && (value < 0))
mexErrMsgTxt("All matrix elements have to be non-negative.");
distMatrix[row] = value;
}
#ifdef CHECK_FOR_INF
/* check for infinite values */
maxFiniteValue = -1;
infiniteValueFound = false;
distMatrixTemp = distMatrix;
while(distMatrixTemp < distMatrixEnd)
{
value = *distMatrixTemp++;
if(mxIsFinite(value))
{
if(value > maxFiniteValue)
maxFiniteValue = value;
}
else
infiniteValueFound = true;
}
if(infiniteValueFound)
{
if(maxFiniteValue == -1) /* all elements are infinite */
return;
/* set all infinite elements to big finite value */
if(maxFiniteValue > 0)
infValue = 10 * maxFiniteValue * nOfElements;
else
infValue = 10;
distMatrixTemp = distMatrix;
while(distMatrixTemp < distMatrixEnd)
if(mxIsInf(*distMatrixTemp++))
*(distMatrixTemp-1) = infValue;
}
#endif
/* memory allocation */
coveredColumns = (bool *)mxCalloc(nOfColumns, sizeof(bool));
coveredRows = (bool *)mxCalloc(nOfRows, sizeof(bool));
starMatrix = (bool *)mxCalloc(nOfElements, sizeof(bool));
primeMatrix = (bool *)mxCalloc(nOfElements, sizeof(bool));
newStarMatrix = (bool *)mxCalloc(nOfElements, sizeof(bool)); /* used in step4 */
/* preliminary steps */
if(nOfRows <= nOfColumns)
{
minDim = nOfRows;
for(row=0; row<nOfRows; row++)
{
/* find the smallest element in the row */
distMatrixTemp = distMatrix + row;
minValue = *distMatrixTemp;
distMatrixTemp += nOfRows;
while(distMatrixTemp < distMatrixEnd)
{
value = *distMatrixTemp;
if(value < minValue)
minValue = value;
distMatrixTemp += nOfRows;
}
/* subtract the smallest element from each element of the row */
distMatrixTemp = distMatrix + row;
while(distMatrixTemp < distMatrixEnd)
{
*distMatrixTemp -= minValue;
distMatrixTemp += nOfRows;
}
}
/* Steps 1 and 2a */
for(row=0; row<nOfRows; row++)
for(col=0; col<nOfColumns; col++)
if(distMatrix[row + nOfRows*col] == 0)
if(!coveredColumns[col])
{
starMatrix[row + nOfRows*col] = true;
coveredColumns[col] = true;
break;
}
}
else /* if(nOfRows > nOfColumns) */
{
minDim = nOfColumns;
for(col=0; col<nOfColumns; col++)
{
/* find the smallest element in the column */
distMatrixTemp = distMatrix + nOfRows*col;
columnEnd = distMatrixTemp + nOfRows;
minValue = *distMatrixTemp++;
while(distMatrixTemp < columnEnd)
{
value = *distMatrixTemp++;
if(value < minValue)
minValue = value;
}
/* subtract the smallest element from each element of the column */
distMatrixTemp = distMatrix + nOfRows*col;
while(distMatrixTemp < columnEnd)
*distMatrixTemp++ -= minValue;
}
/* Steps 1 and 2a */
for(col=0; col<nOfColumns; col++)
for(row=0; row<nOfRows; row++)
if(distMatrix[row + nOfRows*col] == 0)
if(!coveredRows[row])
{
starMatrix[row + nOfRows*col] = true;
coveredColumns[col] = true;
coveredRows[row] = true;
break;
}
for(row=0; row<nOfRows; row++)
coveredRows[row] = false;
}
/* move to step 2b */
step2b(assignment, distMatrix, starMatrix, newStarMatrix, primeMatrix, coveredColumns, coveredRows, nOfRows, nOfColumns, minDim);
/* compute cost and remove invalid assignments */
computeassignmentcost(assignment, cost, distMatrixIn, nOfRows);
/* free allocated memory */
mxFree(distMatrix);