diff --git a/matlab/GoPrintTracks.m b/matlab/GoPrintTracks.m
new file mode 100644
index 0000000000000000000000000000000000000000..2860f224e5ef392a17bee1cdc3326ae174ef4a50
--- /dev/null
+++ b/matlab/GoPrintTracks.m
@@ -0,0 +1,90 @@
+%inputroot = 'C:\Users\mwinter\Documents\Data\OHSU\BDNF';
+%outputroot = fullfile(inputroot, 'OutputTrackingResults');
+%inputroot = 'C:\Users\mwinter\Documents\Data\OHSU\ssNPY';
+%outputroot = fullfile(inputroot, 'OutputTrackingResults');
+inputroot = 'C:\Users\mwinter\Documents\Data\OHSU\BDNF';
+outputroot = fullfile(inputroot, 'OutputTrackingResults');
+
+if ( ~exist(outputroot, 'dir') )
+ mkdir(outputroot);
+end
+
+if ( ~exist('movielist', 'var') )
+ if ( exist(fullfile(inputroot, 'movielist.mat'), 'file') )
+ load(fullfile(inputroot, 'movielist.mat'));
+ else
+ movielist = MakeBDNFMovieList(inputroot);
+ save(fullfile(inputroot, 'movielist.mat'), 'movielist');
+ end
+end
+
+if ( ~exist('startidx', 'var') )
+ startidx = 1;
+end
+
+if ( ~exist('endidx', 'var') )
+ endidx = length(movielist);
+end
+
+disp(['Start time: ' datestr(clock())]);
+
+goodsegs = [];
+for i=startidx:endidx
+ for tryalpha = 0.01:0.01:2.0
+% if ( exist(fullfile(outputroot, [movielist(i).root movielist(i).cellname '_a' num2str(10*tryalpha,'%02d') '.mat']), 'file') )
+% system(['copy "' fullfile(outputroot, [movielist(i).root movielist(i).cellname '_a' num2str(10*tryalpha,'%02d') '.mat']) '" "' fullfile(outputroot, [movielist(i).root movielist(i).cellname '_a' mkalphastr(tryalpha) '.mat']) '"']);
+% delete(fullfile(outputroot, [movielist(i).root movielist(i).cellname '_a' num2str(10*tryalpha,'%02d') '.mat']));
+% end
+% if ( exist(fullfile(outputroot, [movielist(i).root movielist(i).cellname '_a' num2str(10*tryalpha,'%02d') '_seg.dat']), 'file') )
+% system(['copy "' fullfile(outputroot, [movielist(i).root movielist(i).cellname '_a' num2str(10*tryalpha,'%02d') '_seg.dat']) '" "' fullfile(outputroot, [movielist(i).root movielist(i).cellname '_a' mkalphastr(tryalpha) '_seg.dat']) '"']);
+% delete(fullfile(outputroot, [movielist(i).root movielist(i).cellname '_a' num2str(10*tryalpha,'%02d') '_seg.dat']));
+% end
+% if ( exist(fullfile(outputroot, [movielist(i).root movielist(i).cellname '_a' num2str(10*tryalpha,'%02d') '_mb.dat']), 'file') )
+% system(['copy "' fullfile(outputroot, [movielist(i).root movielist(i).cellname '_a' num2str(10*tryalpha,'%02d') '_mb.dat']) '" "' fullfile(outputroot, [movielist(i).root movielist(i).cellname '_a' mkalphastr(tryalpha) '_mb.dat']) '"']);
+% delete(fullfile(outputroot, [movielist(i).root movielist(i).cellname '_a' num2str(10*tryalpha,'%02d') '_mb.dat']));
+% end
+% if ( exist(fullfile(outputroot, [movielist(i).root movielist(i).cellname '_a' num2str(10*tryalpha,'%02d') '_tracks.txt']), 'file') )
+% system(['copy "' fullfile(outputroot, [movielist(i).root movielist(i).cellname '_a' num2str(10*tryalpha,'%02d') '_tracks.txt']) '" "' fullfile(outputroot, [movielist(i).root movielist(i).cellname '_a' mkalphastr(tryalpha) '_tracks.txt']) '"']);
+% delete(fullfile(outputroot, [movielist(i).root movielist(i).cellname '_a' num2str(10*tryalpha,'%02d') '_tracks.txt']));
+% end
+% if ( exist(fullfile(outputroot, [movielist(i).root movielist(i).cellname '_a' num2str(10*tryalpha,'%02d') '_dbginfo.txt']), 'file') )
+% system(['copy "' fullfile(outputroot, [movielist(i).root movielist(i).cellname '_a' num2str(10*tryalpha,'%02d') '_dbginfo.txt']) '" "' fullfile(outputroot, [movielist(i).root movielist(i).cellname '_a' mkalphastr(tryalpha) '_dbginfo.txt']) '"']);
+% delete(fullfile(outputroot, [movielist(i).root movielist(i).cellname '_a' num2str(10*tryalpha,'%02d') '_dbginfo.txt']));
+% end
+
+ if ( ~exist(fullfile(outputroot, [movielist(i).root movielist(i).cellname '_a' mkalphastr(tryalpha) '.mat']), 'file') )
+ continue;
+ end
+
+ goodsegs(end+1,:) = [i tryalpha];
+ end
+end
+
+if ( ~exist(fullfile(outputroot,'Prints'), 'dir') )
+ mkdir(fullfile(outputroot,'Prints'));
+end
+
+for i=1:size(goodsegs,1)
+ mvidx = goodsegs(i,1);
+ alpha = goodsegs(i,2);
+
+ if ( ~exist(fullfile(outputroot, [movielist(mvidx).root movielist(mvidx).cellname '_a' mkalphastr(alpha) '_tracks.txt']), 'file')...
+ || ~exist(fullfile(outputroot,[movielist(mvidx).root movielist(mvidx).cellname '_a' mkalphastr(alpha) '_dbginfo.txt']), 'file') )
+ continue;
+ end
+
+ outdatname = [movielist(mvidx).root movielist(mvidx).cellname '_a' mkalphastr(alpha) '_tracks.txt'];
+ dbgdatname = [movielist(mvidx).root movielist(mvidx).cellname '_a' mkalphastr(alpha) '_dbginfo.txt'];
+ load(fullfile(outputroot, [movielist(mvidx).root movielist(mvidx).cellname '_a' mkalphastr(alpha) '.mat']));
+ trackingStruct = readConnectData(trackingStruct, fullfile(outputroot, outdatname), fullfile(outputroot, dbgdatname));
+ save(fullfile(outputroot, [trackingStruct.fullname '_a' mkalphastr(alpha) '.mat']), 'trackingStruct');
+
+ try
+ plotConnect(trackingStruct,5);
+ print('-r300', '-dtiffn', [fullfile(outputroot,'Prints') filesep trackingStruct.fullname '_a' mkalphastr(alpha) '_tracks.tiff']);
+ catch ME
+ disp(['Problem plotting: ' fullfile(outputroot,'Prints') filesep trackingStruct.fullname '_a' mkalphastr(alpha) '_tracks.tiff']);
+ disp([' :' ME.message]);
+ end
+ close all;
+end
\ No newline at end of file
diff --git a/matlab/GoTrackFakes.m b/matlab/GoTrackFakes.m
new file mode 100644
index 0000000000000000000000000000000000000000..df8bbc6f6000b45b4ed415e972aeeca64e6be30d
--- /dev/null
+++ b/matlab/GoTrackFakes.m
@@ -0,0 +1,143 @@
+outdir = 'C:\Users\mwinter\Documents\MATLAB\faketest_singleaxon_three';
+mkdir(outdir);
+% numsdetlist = [10 25 50 75 100 125 150 200 250 300 350];
+numsdetlist = [10 20 40 60 80 100];
+falsedet = [0 10 20 40 60];
+numtrys = 20;
+
+trackertypes = {'_mta';'_sfa';'_Jaqaman'};
+% trackertypes = ['_sfa'];
+
+% % axons = [];
+% load(fullfile(outdir,'axons.mat'));
+% for idx=1:length(numsdetlist)
+% numdet = numsdetlist(idx);
+%
+% for idxfalse = 1:length(falsedet)
+% pctfalse = falsedet(idxfalse);
+% numfalse = round(numdet*(pctfalse/100));
+%
+% for i=1:numtrys
+% [ds, axons] = makePhantom(512, 512, 'numAxons',1, 'axons',axons, 'numDetections',numdet, 'speedVarProb',0, 'paramVarProb',0, 'falseDetections', numfalse);
+% % load(fullfile(outdir,['fakedata_' num2str(numdet,'%04d') 'det_' num2str(i,'%04d') 'try.mat']));
+%
+% for k=1:length(ds.rgDetect)
+% for j=1:size(ds.rgDetect{i},1)
+% ds.rgDetect{k}(j,:) = [ds.loc{k}(1,j) ds.loc{k}(2,j) ones(1,25)];
+% end
+% end
+%
+% fname_prefix = ['fakedata_' num2str(numdet,'%04d') 'det_' num2str(pctfalse,'%04d') 'false_' num2str(i,'%04d') 'try'];
+%
+% save(fullfile(outdir,[fname_prefix '.mat']), 'ds');
+%
+% for tkidx=1:size(trackertypes,1)
+% outdatname = [fname_prefix trackertypes(tkidx,:) '_tracks.txt'];
+% dbgdatname = [fname_prefix trackertypes(tkidx,:) '_dbginfo.txt'];
+% segdatname = [fname_prefix '_seg.dat'];
+% mbdatname = [fname_prefix '_mb.dat'];
+%
+% if ( ~exist(fullfile(outdir,outdatname), 'file') )
+% writeSegData(fullfile(outdir,segdatname), ds.rgDetect, ds.rgCC, [1 512; 1 512]);
+% writeMotionBlurData(fullfile(outdir,mbdatname), ds.rgMotionBlur);
+%
+% system(['KymoTracker_v2' trackertypes(tkidx) '.exe "' fullfile(outdir,segdatname) '" "' fullfile(outdir,mbdatname) '" "' fullfile(outdir,outdatname) '" "' fullfile(outdir,dbgdatname) '" &']);
+% end
+% end
+% clear ds;
+% end
+% end
+% end
+% save(fullfile(outdir,'axons.mat'), 'axons');
+
+%errstruct = [];
+for tkidx=1:length(trackertypes)
+% if ( ~exist('pctAdds','var') || ~all(size(pctAdds) == [length(numsdetlist) length(falsedet) numtrys]) )
+ pctAdds = nan*ones(length(numsdetlist),length(falsedet),numtrys);
+ pctMisses = nan*ones(length(numsdetlist),length(falsedet),numtrys);
+ pctErrs = nan*ones(length(numsdetlist),length(falsedet),numtrys);
+ avgPixelErrs = nan*ones(length(numsdetlist),length(falsedet),numtrys);
+% else
+% pctAdds = errstruct(tkidx,1).pctAdds;
+% pctMisses = errstruct(tkidx,1).pctMisses;
+% pctErrs = errstruct(tkidx,1).pctErrs;
+% avgPixelErrs = errstruct(tkidx,1).avgPixelErrs;
+% end
+
+ for idx=1:length(numsdetlist)
+ numdet = numsdetlist(idx);
+
+ for idxfalse = 1:length(falsedet)
+ pctfalse = falsedet(idxfalse);
+ numfalse = round(numdet*(pctfalse/100));
+ %numfalse = pctfalse;
+
+ for i=1:numtrys
+
+ fname_prefix = ['fakedata_' num2str(numdet,'%04d') 'det_' num2str(pctfalse,'%04d') 'false_' num2str(i,'%04d') 'try'];
+ %fname_prefix = ['fakedata_' num2str(numdet,'%04d') 'det_' num2str(i,'%04d') '_false' num2str(pctfalse,'%04d') '_try'];
+
+ disp(['Dataset: ' num2str(numdet) ' detections, ' num2str(pctfalse) ' false, ' num2str(i) 'th try'])
+
+ if ( isnan(pctAdds(idx,idxfalse,i)) )
+ if ( strcmpi(trackertypes{tkidx},'_Jaqaman') )
+ outdatname = [fname_prefix trackertypes{tkidx} '_tracks.mat'];
+ else
+ outdatname = [fname_prefix trackertypes{tkidx} '_tracks.txt'];
+ dbgdatname = [fname_prefix trackertypes{tkidx} '_dbginfo.txt'];
+ end
+
+ if ( ~exist(fullfile(outdir,outdatname), 'file') )
+ continue;
+ end
+
+ disp('Loading true phantom data...');
+ tic
+ load(fullfile(outdir,[fname_prefix '.mat']));
+ toc
+
+ if ( exist('dataset','var') )
+ ds = dataset;
+ clear dataset
+ save(fullfile(outdir,[fname_prefix '.mat']), 'ds');
+ end
+
+ intracking.gAssigned = [];
+ intracking.gConnect = [];
+ intracking.gCellData = [];
+ intracking.gCellDbgData = [];
+ intracking.trackList = [];
+
+ disp('Loading phantom tracked data...');
+ tic
+ try
+ if ( strcmpi(trackertypes{tkidx},'_Jaqaman') )
+ load(fullfile(outdir,outdatname));
+ intracking = trackingStruct;
+ else
+ intracking = readConnectData(intracking, fullfile(outdir,outdatname), fullfile(outdir,dbgdatname));
+ end
+ catch me
+ continue;
+ end
+ toc
+
+ disp('Calculating edge errors...');
+ tic
+ [pctAdds(idx,idxfalse,i), pctMisses(idx,idxfalse,i), pctErrs(idx,idxfalse,i), avgPixelErrs(idx,idxfalse,i)] = calcErrs(ds, intracking, numdet, numfalse);
+ toc
+ end
+
+ disp('---------------------');
+ end
+ end
+ end
+ errstruct(tkidx,1).name = trackertypes{tkidx};
+ errstruct(tkidx,1).pctAdds = pctAdds;
+ errstruct(tkidx,1).pctMisses = pctMisses;
+ errstruct(tkidx,1).pctErrs = pctErrs;
+ errstruct(tkidx,1).avgPixelErrs = avgPixelErrs;
+end
+
+save(['trackComparisonErrs_' num2str(max(numsdetlist),'%04d') 'detmax_' num2str(max(falsedet),'%04d') 'falsemax.mat'], 'errstruct', 'numsdetlist', 'falsedet', 'numtrys', 'trackertypes');
+%errSurfaces(100*pctErrs, falsedet, numsdetlist, 30)
diff --git a/matlab/buildMovieInfo.m b/matlab/buildMovieInfo.m
new file mode 100644
index 0000000000000000000000000000000000000000..f4a7dc454f1f6e601b0ec6280f6a25d04196d700
--- /dev/null
+++ b/matlab/buildMovieInfo.m
@@ -0,0 +1,8 @@
+function movieInfo = makeMovieInfo(ds)
+ for i=1:length(ds.rgDetect)
+ movieInfo(i).xCoord = [ds.rgDetect{i}(:,1) zeros(size(ds.rgDetect{i}(:,1)))];
+ movieInfo(i).yCoord = [ds.rgDetect{i}(:,2) zeros(size(ds.rgDetect{i}(:,2)))];
+ movieInfo(i).zCoord = [zeros(size(ds.rgDetect{i}(:,1))) zeros(size(ds.rgDetect{i}(:,1)))];
+ movieInfo(i).amp = [ds.rgDetect{i}(:,15) zeros(size(ds.rgDetect{i}(:,1)))];
+ end
+end
\ No newline at end of file
diff --git a/matlab/buildTrackingStruct.m b/matlab/buildTrackingStruct.m
new file mode 100644
index 0000000000000000000000000000000000000000..d9d94b0a6562ebe7ead8f41790779a0a5af23ab0
--- /dev/null
+++ b/matlab/buildTrackingStruct.m
@@ -0,0 +1,45 @@
+function trackingStruct = buildTrackingStruct(ds, tracksFinal, compWindowSize)
+ l2g = zeros(length(ds.rgDetect)+1,1);
+ g2l = zeros(size(ds.gConnectActual,1),2);
+ for j=1:length(ds.rgDetect)
+ l2g(j+1) = l2g(j) + size(ds.rgDetect{j},1);
+ g2l((l2g(j)+1):l2g(j+1),:) = [j*ones(size(ds.rgDetect{j},1),1) cumsum(ones(size(ds.rgDetect{j},1),1))];
+ end
+
+ tracksList = [];
+ numEdges = length([tracksFinal(:).tracksFeatIndxCG]);
+ gAssignedTracking = sparse([],[],[],size(ds.gConnectActual,1),size(ds.gConnectActual,2),numEdges);
+ gConnectTracking = sparse([],[],[],size(ds.gConnectActual,1),size(ds.gConnectActual,2),numEdges);
+ gCellDatatTracking = [];
+ for i=1:length(tracksFinal)
+ bNonMissedFrames = ~(tracksFinal(i).tracksFeatIndxCG == 0);
+ trackFrames = tracksFinal(i).seqOfEvents(1,1):tracksFinal(i).seqOfEvents(end,1);
+ tracksList{i} = tracksFinal(i).tracksFeatIndxCG + l2g(trackFrames)';
+ tracksList{i} = tracksList{i}(bNonMissedFrames);
+
+ r = tracksList{i}(1:end-1);
+ c = tracksList{i}(2:end);
+ linIdx = sub2ind(size(gAssignedTracking),r,c);
+
+ bAlreadyAssigned = any(gAssignedTracking(r,:),2);
+ if ( any(bAlreadyAssigned) )
+ error('We do not handle shared detections for the Jaqaman tracker yet');
+ end
+
+ gAssignedTracking(linIdx) = 1;
+
+ for j=1:length(r)
+ newCellIdx = length(gCellDatatTracking)+1;
+
+ gConnectTracking(r(j),c(j)) = newCellIdx;
+
+ gCellDatatTracking(newCellIdx).trackID = i;
+ gCellDatatTracking(newCellIdx).path = tracksList{i}(j:min(j+compWindowSize,length(tracksList{i})));
+ end
+ end
+
+ trackingStruct.gConnect = gConnectTracking;
+ trackingStruct.gCellData = gCellDatatTracking;
+ trackingStruct.gAssigned = gAssignedTracking;
+ trackingStruct.trackList = tracksList;
+end
\ No newline at end of file
diff --git a/matlab/calcErrs.m b/matlab/calcErrs.m
new file mode 100644
index 0000000000000000000000000000000000000000..0f4818ccf6ec4da3a4834776fdd14a052804c1b3
--- /dev/null
+++ b/matlab/calcErrs.m
@@ -0,0 +1,126 @@
+function [pctAdds, pctMisses, pctErrs, avgPixelErr, errs] = calcErrs(ds, intracking, numdet, numfalse)
+ l2g = zeros(length(ds.rgDetect)+1,1);
+ g2l = zeros(size(ds.gConnectActual,1),2);
+ for j=1:length(ds.rgDetect)
+ l2g(j+1) = l2g(j) + size(ds.rgDetect{j},1);
+ g2l((l2g(j)+1):l2g(j+1),:) = [j*ones(size(ds.rgDetect{j},1),1) cumsum(ones(size(ds.rgDetect{j},1),1))];
+ end
+
+ gAssignedActual = ds.gConnectActual > 0;
+ gAssignedTracking = intracking.gAssigned;
+ gAssignedExpandedShared = gAssignedTracking;
+
+ mintracklength = 6;
+
+ % Expand shared edges
+ % Note: Could also expand missed detections by interpolation
+ [r,c] = find(gAssignedTracking);
+ sharedidx = find(c > l2g(g2l(r,1)+2));
+ for i=1:length(sharedidx)
+ ridx = r(sharedidx(i));
+ cidx = c(sharedidx(i));
+ cellidx = intracking.gConnect(ridx,find(intracking.gAssigned(ridx,:)));
+ if ( length(cellidx) > 1 )
+ error(['Multiple edge assignments in gConnect row: ' num2str(i)]);
+ end
+% if ( isempty(cellidx) )
+% continue;
+% end
+ % Note: Could also expand missed detections by interpolation
+ % if path has no internal points (>ridx,<cidx) then interp
+
+ path = intracking.gCellData(cellidx).path;
+ path = path(path >= ridx & path <= cidx);
+ linidx = sub2ind(size(gAssignedExpandedShared), path(1:end-1),path(2:end));
+
+ gAssignedExpandedShared(linidx) = 1;
+ gAssignedExpandedShared(ridx,cidx) = 0;
+ end
+
+ edgecount = nnz(gAssignedActual);
+
+ misscount = 0;
+ errs = [];
+ swaperr = [];
+ swaperrcount = 1;
+ addcount = 0;
+ falseaddcount = 0;
+ for i=1:size(gAssignedActual,1)
+ t = g2l(i,1);
+ locidx = g2l(i,2);
+
+ % No error if assigned edges match perfectly
+ if ( all((gAssignedActual(i,:) - gAssignedTracking(i,:)) == 0) )
+ continue;
+ end
+
+ cellidx = intracking.gConnect(i,find(intracking.gAssigned(i,:)));
+ tracklen = mintracklength;
+ if ( ~isempty(cellidx) )
+ tracklen = length(intracking.trackList{intracking.gCellData(cellidx).trackID});
+ end
+
+ bIsFalse = (locidx > numdet);
+
+ % Ignore short track edge errors in false
+ if ( bIsFalse && tracklen < mintracklength )
+ continue;
+ end
+
+ comprow = any(gAssignedActual(i,:)) - any(gAssignedExpandedShared(i,:));
+
+ trueinidx = find(gAssignedActual(i,:));
+ trackedinidx = find(gAssignedExpandedShared(i,:));
+ if ( comprow > 0 && ~any(gAssignedActual(:,trueinidx)) )
+ misscount = misscount + 1;
+ elseif (comprow < 0 && ~any(any(gAssignedActual(:,trackedinidx))) )
+ addcount = addcount + 1;
+ if ( bIsFalse )
+ falseaddcount = falseaddcount + 1;
+ end
+ elseif ( any(gAssignedExpandedShared(i,:)) )
+ % Errors are associated with the row they occur on not column?
+ trackedinidx = find(gAssignedExpandedShared(i,:));
+ [trueoutidx,c] = find(gAssignedActual(:,trackedinidx));
+
+ trueoutidx = trueoutidx';
+
+ trueedges = [i*ones(size(trueinidx)) trueoutidx; trueinidx trackedinidx(c)];
+ trackedges = [i*ones(size(trackedinidx)); trackedinidx];
+
+ trackerr = 0;
+ disterr = Inf*ones(size(trueedges,2),1);
+ for j=1:size(trackedges,2)
+ for k=1:size(trueedges,2)
+ if ( any(trueedges(:,k) == trackedges(:,j)) )
+ truepts(:,1) = ds.loc{g2l(trueedges(1,k),1)}(:,g2l(trueedges(1,k),2));
+ truepts(:,2) = ds.loc{g2l(trueedges(2,k),1)}(:,g2l(trueedges(2,k),2));
+ trackpts(:,1) = ds.loc{g2l(trackedges(1,j),1)}(:,g2l(trackedges(1,j),2));
+ trackpts(:,2) = ds.loc{g2l(trackedges(2,j),1)}(:,g2l(trackedges(2,j),2));
+
+ disterr(k) = sum(sqrt((truepts(1,:)-trackpts(1,:)).^2 + (truepts(2,:)-trackpts(2,:)).^2));
+ end
+ end
+ trackerr = trackerr + min(disterr);
+ end
+
+ if ( ~isempty(trueinidx) )
+ errs{swaperrcount,1} = [trackerr i trueinidx trackedinidx];
+ else
+ errs{swaperrcount,1} = [trackerr i 0 trackedinidx];
+ end
+
+ swaperr(swaperrcount,1) = trackerr;
+ swaperrcount = swaperrcount + 1;
+ end
+ end
+
+ pctAdds = addcount / edgecount;
+ pctMisses = misscount / edgecount;
+ pctErrs = swaperrcount / edgecount;
+
+ avgPixelErr = 0;
+ if ( ~isempty(swaperr) )
+ avgPixelErr = mean(swaperr);
+ end
+end
\ No newline at end of file
diff --git a/matlab/compareErrSurf.m b/matlab/compareErrSurf.m
new file mode 100644
index 0000000000000000000000000000000000000000..4ede69e41c2de2a8260b5b3bbde42dbcc5c583be
--- /dev/null
+++ b/matlab/compareErrSurf.m
@@ -0,0 +1,46 @@
+function compareErrSurf(erra, errb, errc, gridx, gridy, interppts, interpmethod)
+
+ cmap = [summer(512);autumn(512);winter(512)];
+
+ [X,Y] = meshgrid(gridx, gridy);
+ [XI,YI] = meshgrid(linspace(gridx(1), gridx(end), interppts), linspace(gridy(1), gridy(end), interppts));
+
+ mederra = nanmedian(erra,3);
+ interpmederra = interp2(X,Y,mederra,XI,YI, interpmethod);
+ mederrb = nanmedian(errb,3);
+ interpmederrb = interp2(X,Y,mederrb,XI,YI, interpmethod);
+ mederrc = nanmedian(errc,3);
+ interpmederrc = interp2(X,Y,mederrc,XI,YI, interpmethod);
+
+ ca = 511*((mederra-min(mederra(:))) / (max(mederra(:))-min(mederra(:)))) + 1;
+ cb = 511*((mederrb-min(mederrb(:))) / (max(mederrb(:))-min(mederrb(:)))) + 513;
+ cc = 511*((mederrc-min(mederrc(:))) / (max(mederrc(:))-min(mederrc(:)))) + 1025;
+
+ colora = 511*((interpmederra-min(interpmederra(:))) / (max(interpmederra(:))-min(interpmederra(:)))) + 1;
+ colorb = 511*((interpmederrb-min(interpmederrb(:))) / (max(interpmederrb(:))-min(interpmederrb(:)))) + 513;
+ colorc = 511*((interpmederrc-min(interpmederrb(:))) / (max(interpmederrc(:))-min(interpmederrb(:)))) + 1025;
+
+
+
+ figure;colormap(cmap);
+ surf(X,Y,mederra,ca);
+ hold on;
+ surf(X,Y,mederrb,cb);
+ surf(X,Y,mederrc,cc);
+ xlabel('False Detections(% of Detections)');
+ ylabel('Detections per Frame');
+
+
+ figure;colormap(cmap);
+ h = surf(XI,YI,interpmederra,colora);
+ hold on;
+ h = surf(XI,YI,interpmederrb,colorb);
+ alpha(h, 0.85);
+ h = surf(XI,YI,interpmederrc,colorc);
+ alpha(h, 0.75);
+
+ xlabel('False Detections(% of Detections)');
+ ylabel('Detections per Frame');
+ %legend('Multitemporal Assignment','Single Frame Assignment');
+
+end
\ No newline at end of file
diff --git a/matlab/compareKymo.m b/matlab/compareKymo.m
new file mode 100644
index 0000000000000000000000000000000000000000..8d2d269e8a5d2e3086bc7443c9567ae7331e363e
--- /dev/null
+++ b/matlab/compareKymo.m
@@ -0,0 +1,269 @@
+function outColors = compareKymo(root, fname, trackingStruct, tracks, times)
+
+ rmkdir(root);
+
+ cmap = hsv(length(tracks)+1);
+ ctweak = [zeros(size(cmap,1),1) zeros(size(cmap,1),1) zeros(size(cmap,1),1)];
+ ctweak(4:6,3) = [0 0.3 0.4];
+ cmap = hsv2rgb(rgb2hsv(cmap)-ctweak);
+ trackColors = zeros(length(tracks),3);
+ for i=1:length(tracks)
+ trackColors(i,:) = cmap(i,:);
+ end
+ outColors = round(255*trackColors);
+
+ fullTrackList = BuildFullTrackList(trackingStruct);
+
+ kymosize = size(trackingStruct.imKymo);
+
+ tracelen = norm([trackingStruct.jPath(end) trackingStruct.iPath(end)] - [trackingStruct.jPath(1) trackingStruct.iPath(1)]);
+ trcppix = tracelen / kymosize(1);
+
+ umppix = 0.18;
+ spfrm = 0.650;
+
+ invfig = figure('Color','w');[rgbKymo rgbmap] = gray2ind(1-mat2gray(trackingStruct.imKymo), 255);
+ imagesc([0 spfrm*(kymosize(2)-1)], [0 umppix*(kymosize(1)-1)*trcppix],ind2rgb(rgbKymo, rgbmap));axis xy; truesize(invfig, round(2*(kymosize./[trcppix 1])))
+ frm = getframe(invfig);
+ [frmimg, map] = frame2im(frm);
+ if ( ~isempty(map) )
+ imwrite(frmimg, map, fullfile(root, [fname '_kymo.tiff']), 'tiff');
+ else
+ imwrite(frmimg, fullfile(root, [fname '_kymo.tiff']), 'tiff');
+ end
+
+ trackfig = figure('Color','w');[rgbKymo rgbmap] = gray2ind(trackingStruct.imKymo, 255);
+ imagesc([0 spfrm*(kymosize(2)-1)], [0 umppix*(kymosize(1)-1)*trcppix],ind2rgb(rgbKymo, rgbmap));axis xy;
+
+ vdata = [];
+ kvdata = [];
+
+ figure(trackfig);hold on;
+ for i=1:length(tracks)
+ trkt = trackingStruct.globNodeLUT(fullTrackList{tracks(i)},1)';
+ trkidx = trackingStruct.globNodeLUT(fullTrackList{tracks(i)},2)';
+
+ bValidId = trkt >= times(i,1) & trkt <= times(i,2);
+ trkt = trkt(bValidId);
+ trkidx = trkidx(bValidId);
+ trkpts = GetKymoPt(trackingStruct.rgDetect, trkidx', trkt', trackingStruct.LKymo, [trackingStruct.jPath';trackingStruct.iPath']);
+
+ plot(spfrm*(trkt-1), umppix*(trkpts-1)*trcppix, 'Color', trackColors(i,:), 'LineWidth',2, 'LineStyle','-');
+ text(spfrm*(trkt(round(end/2))-1), umppix*(trkpts(round(end/2))-1)*trcppix, ['\leftarrow ' num2str(i)], 'Color', trackColors(i,:), 'BackgroundColor','w');
+
+ lclSpeed = [];
+ vtan = [];
+ xldata = {i '';'Velocity' 'Kymograph Velocity';'' ''};
+ headeroffset = size(xldata,1);
+ arclen = 0;
+ for k=1:length(trkt)-1
+ lclVel = ((trackingStruct.rgDetect{trkt(k+1)}(trkidx(k+1),1:2) - trackingStruct.rgDetect{trkt(k)}(trkidx(k),1:2)) / (trkt(k+1)-trkt(k)) * (umppix / spfrm));
+ lclVel = [lclVel(2),lclVel(1)];
+
+ lclDelta = (trackingStruct.rgDetect{trkt(k+1)}(trkidx(k+1),1:2) - trackingStruct.rgDetect{trkt(k)}(trkidx(k),1:2));
+ lclDelta = [lclDelta(2),lclDelta(1)];
+
+ lclKymoVel = trkpts(k+1) - trkpts(k);
+ kymovec = [trackingStruct.jPath(trkpts(k+1)) trackingStruct.iPath(trkpts(k+1))] - [trackingStruct.jPath(trkpts(k)) trackingStruct.iPath(trkpts(k))];
+
+ lclSpeed(k,1) = sign(lclKymoVel)*norm(lclVel);
+
+ lclArcLen = dot(lclDelta, kymovec) / norm(kymovec);
+ vtan(k,1) = sign(lclKymoVel)*dot(lclVel, kymovec) / norm(kymovec);
+ if ( norm(kymovec) < 1e-7 )
+ vtan(k,1) = 0;
+ lclArcLen = 0;
+ end
+ xldata{k+headeroffset,1} = lclSpeed(k,1);
+ xldata{k+headeroffset,2} = vtan(k,1);
+
+ arclen = arclen + lclArcLen;
+ end
+ %xlswrite(fullfile(root, [fname '_stats.xlsx']), xldata, [char((2*(i-1)+1)+65) num2str(1) ':' char((2*(i-1)+2)+65) num2str(length(lclSpeed)+headeroffset)]);
+
+ avgSpeed = (norm((trackingStruct.rgDetect{trkt(end)}(trkidx(end),1:2) - trackingStruct.rgDetect{trkt(1)}(trkidx(1),1:2))) / (trkt(end)-trkt(1))) * (umppix / spfrm);
+ avgKymoSpeed = (arclen / (trkt(end)-trkt(1))) * (umppix / spfrm);
+
+ vdata = [vdata;lclSpeed];
+ kvdata = [kvdata;vtan];
+ end
+ hold off;
+ hleg = legend(arrayfun(@num2str,(1:length(tracks)), 'UniformOutput',0), 'Location','NorthEastOutside');
+ %legend(hleg, 'boxoff');
+ truesize(trackfig, round(3*(kymosize./[trcppix 1])));
+ %legpos = get(hleg, 'Position');
+ %set(hleg, 'Position',legpos - [0.04,0,0,0]);
+ frm = getframe(trackfig);
+ [frmimg, map] = frame2im(frm);
+ if ( ~isempty(map) )
+ imwrite(frmimg, map, fullfile(root, [fname '_tracks.tiff']), 'tiff');
+ else
+ imwrite(frmimg, fullfile(root, [fname '_tracks.tiff']), 'tiff');
+ end
+
+ xbins = -5.5:0.25:4.5;
+ figure();hist(vdata(abs(vdata)>0.2),xbins);
+ title(['ssNPY velocities']);
+ xlabel('Velocity [um/s] (+/-: anterograde/retrograde)');
+ h = findobj(gca,'Type','patch');
+ set(h, 'FaceColor','r');
+ set(h, 'FaceAlpha',0.5);
+ set(gca,'xminortick','on')
+ set(gca,'xtick',[-6:6])
+ print('-dtiff', fullfile(root, [fname '_velhistogram.tiff']));
+
+ xbins = -5.5:0.25:4.5;
+ figure();hist(kvdata(abs(kvdata)>0.2),xbins);
+ title(['ssNPY kymograph velocities']);
+ xlabel('Velocity [um/s] (+/-: anterograde/retrograde)');
+ h = findobj(gca,'Type','patch');
+ set(h, 'FaceColor','r');
+ set(h, 'FaceAlpha',0.5);
+ set(gca,'xminortick','on')
+ set(gca,'xtick',[-6:6])
+ print('-dtiff', fullfile(root, [fname '_kymovelhistogram.tiff']));
+
+ close all;
+
+ %movie
+ iminfo = imfinfo(fullfile(trackingStruct.moviepath, trackingStruct.moviename));
+ kymoinfo = imfinfo(fullfile(trackingStruct.moviepath, trackingStruct.moviename));
+ N = length(iminfo);
+
+ [tracemin, tracemax] = getTrackExtents(trackingStruct, N);
+
+ imwght = ones(iminfo(1).Height, iminfo(1).Width);
+
+ imheight = tracemax(2) - tracemin(2);
+ kymoheight = size(trackingStruct.imKymo,1);
+
+ bLeft = boolean(ones(length(tracks),1));
+ for i=1:length(tracks)
+ trkt = trackingStruct.globNodeLUT(fullTrackList{tracks(i)},1)';
+ trkidx = trackingStruct.globNodeLUT(fullTrackList{tracks(i)},2)';
+
+ bValidId = trkt >= times(i,1) & trkt <= times(i,2);
+ trkt = trkt(bValidId);
+ trkidx = trkidx(bValidId);
+
+ trkpts = GetKymoPt(trackingStruct.rgDetect, trkidx', trkt', trackingStruct.LKymo, [trackingStruct.jPath';trackingStruct.iPath']);
+ n2 = round(length(trkt)/2);
+ if ( trackingStruct.jPath(round(trkpts(n2))) < trackingStruct.rgDetect{trkt(n2)}(trkidx(n2),2) )
+ bLeft(i) = 0;
+ end
+ end
+
+ imfig = figure();
+ for t=1:N
+ %skymoFrame = drawKymo();
+ imFrame = drawTracking();
+
+ %fullFrame = kymoFrame;
+ %fullFrame.cdata = [kymoFrame.cdata imFrame.cdata];
+
+ %movfile = addframe(movfile,fullFrame);
+ [frmimg, map] = frame2im(imFrame);
+ if ( ~isempty(map) )
+ imwrite(frmimg, map, fullfile(root, [fname '_frm_' num2str(t,'%03d') '.tiff']), 'tiff');
+ else
+ imwrite(frmimg, fullfile(root, [fname '_frm_' num2str(t,'%03d') '.tiff']), 'tiff');
+ end
+ end
+
+ close all;
+
+ function frm = drawTracking()
+ figure(imfig);
+ im = imread(fullfile(trackingStruct.moviepath, trackingStruct.moviename), t);
+ im = mat2gray(im);
+ imagesc(im.*imwght);colormap(gray);hold on;
+ brighten(0.2);
+
+ plot(trackingStruct.jPath(1), trackingStruct.iPath(1), 'og');
+ plot(trackingStruct.jPath(end), trackingStruct.iPath(end), 'xg');
+
+ for i=1:length(tracks)
+ trkt = trackingStruct.globNodeLUT(fullTrackList{tracks(i)},1)';
+ trkidx = trackingStruct.globNodeLUT(fullTrackList{tracks(i)},2)';
+
+ bValidId = trkt >= times(i,1) & trkt <= times(i,2);
+ trkt = trkt(bValidId);
+ trkidx = trkidx(bValidId);
+
+ idx = find(trkt == t,1);
+ if ( isempty(idx) )
+ continue;
+ end
+
+ trkpts = GetKymoPt(trackingStruct.rgDetect, trkidx', trkt', trackingStruct.LKymo, [trackingStruct.jPath';trackingStruct.iPath']);
+
+ x=trackingStruct.rgDetect{t}(trkidx(idx),2);
+ y=trackingStruct.rgDetect{t}(trkidx(idx),1);
+
+ bwim = zeros(size(im));
+ r = trackingStruct.rgCC{t}{trkidx(idx)}(:,1);
+ c = trackingStruct.rgCC{t}{trkidx(idx)}(:,2);
+ linidx = sub2ind(size(bwim), r,c);
+ bwim(linidx) = 1;
+ bwbound = bwboundaries(bwim);
+ plot(bwbound{1}(:,2),bwbound{1}(:,1), 'Color',trackColors(i,:), 'LineStyle','-');
+
+ if ( bLeft(i) )
+ text(x+5,y,['\leftarrow ' num2str(i)], 'FontSize',10, 'HorizontalAlignment','Left', 'Color',trackColors(i,:));
+ else
+ text(x-5,y,[num2str(i) ' \rightarrow'], 'FontSize',10, 'HorizontalAlignment','Right', 'Color',trackColors(i,:));
+ end
+ end
+
+ %text(tracemin(1),tracemin(2), ['t=' num2str(t)], 'HorizontalAlign','Left', 'VerticalAlign','Top', 'FontSize',8, 'BackgroundColor',[0 0.6 0.5]);
+ text(2,2, ['t=' num2str(t)], 'HorizontalAlign','Left', 'VerticalAlign','Top', 'FontSize',8, 'BackgroundColor',[0 0.6 0.5]);
+
+ %xlim([tracemin(1) tracemax(1)]);
+ %ylim([tracemin(2) tracemax(2)]);
+
+ hold off;
+ drawnow;
+
+ %truesize(imfig, imsize);
+ imsize = size(im);
+ truesize(imfig, 2*imsize);
+ frm = getframe(gca);
+ end
+end
+
+function fullTrackList = BuildFullTrackList(trackingStruct)
+ fullTrackList = trackingStruct.trackList;
+
+ for i=1:length(trackingStruct.trackList)
+ trkt = trackingStruct.globNodeLUT(trackingStruct.trackList{i},1)';
+
+ deltat = diff(trkt);
+ sharedidx = find(deltat > 1);
+ tshift = 0;
+ for k =1 :length(sharedidx)
+ srcNode = trackingStruct.trackList{i}(sharedidx(k));
+ destNode = trackingStruct.trackList{i}(sharedidx(k)+1);
+
+ srcPath = trackingStruct.gCellData(trackingStruct.gConnect(srcNode,destNode)).path;
+
+ trkt = [trkt(1:sharedidx(k)+tshift) trackingStruct.globNodeLUT(srcPath(2:deltat(sharedidx(k))),1)' trkt(sharedidx(k)+tshift+1:end)];
+ fullTrackList{i} = [fullTrackList{i}(1:sharedidx(k)+tshift) srcPath(2:deltat(sharedidx(k))) fullTrackList{i}(sharedidx(k)+tshift+1:end)];
+
+ tshift = tshift + deltat(sharedidx(k))-1;
+ end
+ end
+end
+
+function [tracemin, tracemax] = getTrackExtents(trackingStruct, N)
+ tracemin = [Inf Inf];
+ tracemax = [0 0];
+ for t=1:N
+ xmin = min(trackingStruct.rgDetect{t}(:,2));
+ xmax = max(trackingStruct.rgDetect{t}(:,2));
+ ymin = min(trackingStruct.rgDetect{t}(:,1));
+ ymax = max(trackingStruct.rgDetect{t}(:,1));
+
+ tracemin = [min(tracemin(1), xmin) min(tracemin(2), ymin)];
+ tracemax = [max(tracemax(1), xmax) max(tracemax(2), ymax)];
+ end
+end
\ No newline at end of file
diff --git a/matlab/compareTracking.m b/matlab/compareTracking.m
new file mode 100644
index 0000000000000000000000000000000000000000..3fd17e6766f3f8cc4f7752ec776c0918f79dbcbd
--- /dev/null
+++ b/matlab/compareTracking.m
@@ -0,0 +1,162 @@
+function compareTracking(ds, trackingStruct, extents, selectTrack)
+ gAssignedActual = (ds.gConnectActual > 0);
+
+ if ( ~exist('selectTrack','var') )
+ selectTrack = [];
+ end
+
+ cmap = [cool(256); spring(256); summer(256); winter(256)];
+ trackColorsActual = cmap(round(1023*rand(length(ds.tracksActual),1)+1),:);
+ trackColorsTest = cmap(round(1023*rand(length(trackingStruct.trackList),1)+1),:);
+
+ fig = figure;
+ set(fig, 'WindowScrollWheelFcn', @scrollWheelCallback);
+ set(fig, 'KeyPressFcn', @keyReleaseCallback);
+
+ t = 1;
+ NFRAMES = length(ds.rgDetect);
+
+ testlims = extents;
+ actuallims = extents;
+
+ l2g = zeros(length(ds.rgDetect)+1,1);
+ g2l = zeros(size(ds.gConnectActual,1),2);
+ for i=1:length(ds.rgDetect)
+ l2g(i+1) = l2g(i) + size(ds.rgDetect{i},1);
+ g2l((l2g(i)+1):l2g(i+1),:) = [i*ones(size(ds.rgDetect{i},1),1) cumsum(ones(size(ds.rgDetect{i},1),1))];
+ end
+
+ %Error types
+ bCompRows = (any(gAssignedActual,2) - any(trackingStruct.gAssigned,2));
+
+ bMisses = (bCompRows > 0);
+ bAdds = (bCompRows < 0);
+ bGood = (bCompRows == 0);
+
+ bMistakes = bGood;
+ bMistakes(bGood) = any(gAssignedActual(bGood,:) - trackingStruct.gAssigned(bGood,:),2);
+
+ dispFrame(testlims, actuallims);
+
+
+ function scrollWheelCallback(src, evt)
+ if ( evt.VerticalScrollCount < 0 )
+ t = t - 1;
+ else
+ t = t + 1;
+ end
+
+ if ( t > NFRAMES )
+ t = NFRAMES;
+ elseif ( t < 1 )
+ t = 1;
+ end
+
+ h = subplot(1,2,1);
+ actuallims = [xlim(h);ylim(h)];
+ h = subplot(1,2,2);
+ testlims = [xlim(h);ylim(h)];
+
+ dispFrame(actuallims, testlims);
+ end
+
+ function dispFrame(leftlims, rightlims)
+ i = [];
+ dotradius = 4;
+ traillength = 3;
+
+ pdfmax = normpdf(0,0,dotradius/3)^2;
+
+ fakeim = zeros(extents(2,2),extents(1,2));
+ for j=1:size(ds.loc{t},2)
+ fakeim = drawGaussianPt(fakeim, ds.loc{t}(1,j), ds.loc{t}(2,j), 1, dotradius, dotradius);
+ end
+
+ h = subplot(1,2,1);
+ cla;
+ imshow(0.7*mat2gray(fakeim));colormap(gray);hold on;
+ %imshow(fakeim);colormap(gray);hold on;
+ for j=1:length(ds.axons)
+ axnpts = ppval(ds.axons(j).pp, 0:ds.axons(j).len+2);
+ plot(axnpts(1,:), axnpts(2,:), '-b');
+ end
+ plot(ds.loc{t}(1,:),ds.loc{t}(2,:), '.g');
+ for j=1:length(ds.tracksActual)
+ %for j = 7:7
+ trkt = g2l(ds.tracksActual{j},1);
+ trkidx = g2l(ds.tracksActual{j},2);
+ idx = find((trkt >= t-traillength) & (trkt <= t));
+ if ( ~isempty(idx) )
+ loc = [];
+ for k=1:length(idx)
+ loc(:,k) = ds.loc{trkt(idx(k))}(:,trkidx(idx(k)));
+ end
+ plot(loc(1,:), loc(2,:), '-', 'Color',trackColorsActual(j,:));
+ if ( trkt(idx(end)) == t )
+% plot(loc(1,end), loc(2,end), 'o', 'Color',trackColorsActual(j,:));
+ plot(loc(1,end), loc(2,end), '.r');
+% text(loc(1,end)+4, loc(2,end)+4, num2str(j), 'Color',trackColorsActual(j,:));
+ end
+
+% for k=1:length(idx)
+% if ( bMisses(ds.tracksActual{j}(idx(k))) || bMistakes(ds.tracksActual{j}(idx(k))) )
+% if ( k < length(idx) )
+% errloc = [ds.loc{trkt(idx(k))}(:,trkidx(idx(k))) ds.loc{trkt(idx(k+1))}(:,trkidx(idx(k+1)))];
+% plot(errloc(1,:), errloc(2,:), '-r', 'LineWidth',2);
+% else
+% errloc = ds.loc{trkt(idx(k))}(:,trkidx(idx(k)));
+% plot(errloc(1,:), errloc(2,:), 'xr', 'LineWidth',2);
+% end
+% end
+% end
+ end
+ end
+ hold off;
+ xlim(h,leftlims(1,:));
+ ylim(h,leftlims(2,:));
+
+ h = subplot(1,2,2);
+ cla;
+ imshow(mat2gray(fakeim));colormap(gray);hold on;
+ %imshow(fakeim);colormap(gray);hold on;
+% for j=1:length(ds.axons)
+% axnpts = ppval(ds.axons(j).pp, 0:ds.axons(j).len+2);
+% plot(axnpts(1,:), axnpts(2,:), '-b');
+% end
+% for j=1:length(trackingStruct.trackList)
+% trkt = g2l(trackingStruct.trackList{j},1);
+% trkidx = g2l(trackingStruct.trackList{j},2);
+% idx = find((trkt >= t-traillength) & (trkt <= t));
+% if ( ~isempty(idx) )
+%
+% loc = [];
+% for k=1:length(idx)
+% loc(:,k) = ds.loc{trkt(idx(k))}(:,trkidx(idx(k)));
+% end
+% plot(loc(1,:), loc(2,:), '-', 'Color',trackColorsTest(j,:));
+% if ( trkt(idx(end)) == t )
+% plot(loc(1,end), loc(2,end), 'o', 'Color',trackColorsTest(j,:));
+% text(loc(1,end)+4, loc(2,end)+4, num2str(j), 'Color',trackColorsTest(j,:));
+% end
+%
+% for k=1:length(idx)
+% if ( bAdds(trackingStruct.trackList{j}(idx(k))) || bMistakes(trackingStruct.trackList{j}(idx(k))) )
+% if ( k < length(idx) )
+% errloc = [ds.loc{trkt(idx(k))}(:,trkidx(idx(k))) ds.loc{trkt(idx(k+1))}(:,trkidx(idx(k+1)))];
+% plot(errloc(1,:), errloc(2,:), '-r', 'LineWidth',2);
+% else
+% errloc = ds.loc{trkt(idx(k))}(:,trkidx(idx(k)));
+% plot(errloc(1,:), errloc(2,:), 'xr', 'LineWidth',2);
+% end
+% end
+% end
+% end
+% end
+ hold off;
+ xlim(h,rightlims(1,:));
+ ylim(h,rightlims(2,:));
+
+ title(num2str(t));
+
+ end
+end
\ No newline at end of file
diff --git a/matlab/drawGaussianPt.m b/matlab/drawGaussianPt.m
new file mode 100644
index 0000000000000000000000000000000000000000..8cbf67aae69a8149d7fd8a7da3fa7fcd36b66617
--- /dev/null
+++ b/matlab/drawGaussianPt.m
@@ -0,0 +1,16 @@
+function im = drawGaussianPt(im, x, y, intensity, xRadius, yRadius)
+ grdx = round(x-xRadius-1):round(x+xRadius+1);
+ validx = (grdx >= 1 & grdx <= size(im,2));
+ grdx = grdx(validx);
+
+ grdy = round(y-yRadius-1):round(y+yRadius+1);
+ validy = (grdy >= 1 & grdy <= size(im,1));
+ grdy = grdy(validy);
+
+ pdfmax = normpdf(0,0,xRadius/3)*normpdf(0,0,yRadius/3);
+
+ [X,Y] = meshgrid(grdx, grdy);
+ alpha = normpdf(X,x,xRadius/3).*normpdf(Y,y,yRadius/3) / pdfmax;
+
+ im(grdy,grdx) = ((1-alpha).*im(grdy,grdx) + alpha*intensity);
+end
\ No newline at end of file
diff --git a/matlab/errSurfaces.m b/matlab/errSurfaces.m
new file mode 100644
index 0000000000000000000000000000000000000000..7fdd9c327254e17976af1c831c1df6c3de1f70ce
--- /dev/null
+++ b/matlab/errSurfaces.m
@@ -0,0 +1,16 @@
+function errSurfaces(err, gridx, gridy, interppts, cmap)
+ if ( ~exist('cmap','var') )
+ cmap = jet(1024);
+ end
+
+ [X,Y] = meshgrid(gridx, gridy);
+ [XI,YI] = meshgrid(linspace(gridx(1), gridx(end), interppts), linspace(gridy(1), gridy(end), interppts));
+ mederr = nanmedian(err,3);
+ interpmederr = interp2(X,Y,mederr,XI,YI, 'bicubic');
+
+ hc = round((length(cmap)-1)*(mederr - min(mederr(:)))/(max(mederr(:))-min(mederr(:)))) + 1;
+ hci = round((length(cmap)-1)*(interpmederr - min(mederr(:)))/(max(mederr(:))-min(mederr(:)))) + 1;
+
+ figure;surf(X,Y,mederr,hc);
+ figure;surf(XI,YI,interpmederr,hci);
+end
\ No newline at end of file
diff --git a/matlab/getTrackStats.m b/matlab/getTrackStats.m
new file mode 100644
index 0000000000000000000000000000000000000000..4c41f15d8a1627e7f2166c93dac658f1f1628dea
--- /dev/null
+++ b/matlab/getTrackStats.m
@@ -0,0 +1,8 @@
+function vs = getTrackStats(trackingStruct, trackIdx, names, stats)
+ dataidx = find(strcmp(names,trackingStruct.fullname));
+ if ( isempty(dataidx) )
+ return;
+ end
+
+ vs = stats{dataidx}.vs{trackIdx};
+end
\ No newline at end of file
diff --git a/matlab/makePhantom.m b/matlab/makePhantom.m
new file mode 100644
index 0000000000000000000000000000000000000000..9ed569619fe813703d115b916fb9386cbfd39d18
--- /dev/null
+++ b/matlab/makePhantom.m
@@ -0,0 +1,524 @@
+function [ds, axons, images] = makePhantom(xlen, ylen, varargin)
+
+ p = inputParser;
+ p.StructExpand = false;
+
+ p.addRequired('xlen');
+ p.addRequired('ylen');
+
+ p.addParamValue('numAxons',0);
+ p.addParamValue('axons',[]);
+ p.addParamValue('numDetections',100);
+ p.addParamValue('numFrames',120);
+
+ p.addParamValue('speedRange',[0 10]);
+ p.addParamValue('speedVarTime',[3 5]);
+ p.addParamValue('speedVarProb',0.1);
+
+ p.addParamValue('axonWidth',8);
+ p.addParamValue('paramVarTime',[7 10]);
+ p.addParamValue('paramVarProb',0.05);
+
+ p.addParamValue('falseDetections',0);
+
+ p.addParamValue('frameTime',0.650);
+ p.addParamValue('exposureTime',0.150);
+
+ p.addParamValue('maxCurve',20);
+
+ p.addParamValue('dotRadius',6);
+
+ p.parse(xlen, ylen, varargin{:});
+
+ args = p.Results;
+
+ tic;
+
+ expratio = args.exposureTime / args.frameTime;
+
+ %dot size?
+
+ extents = [1 xlen;1 ylen];
+
+ numframes = args.numFrames;
+ numdetections = args.numDetections;
+ numaxons = args.numAxons;
+
+ axonwidth = args.axonWidth;
+
+ falsedetpad = 3*args.falseDetections;
+ gConnectActual = sparse([], [], [], numframes*(numdetections+falsedetpad), numframes*(numdetections+falsedetpad), numframes*numdetections);
+
+ % Generate initial detections and associated parameters
+ detLoc{1} = [xlen;ylen]*ones(1,numdetections).*rand(2,numdetections);
+ if ( numaxons > 0 )
+ if ( ~isempty(args.axons) )
+ axons = args.axons;
+ else
+ axons = MakeAxons(extents, args.maxCurve);
+ end
+ detParam{1} = [ceil((numaxons-eps)*rand(1,numdetections) + eps); 2*axonwidth*rand(1,numdetections)-axonwidth; rand(1,numdetections); 2*round(rand(1,numdetections))-1];
+ detParam{1}(3,:) = detParam{1}(3,:).*[axons(detParam{1}(1,:)).len];
+
+ for j=1:numdetections
+ detLoc{1}(:,j) = GetAxonLoc(detParam{1}(:,j));
+ bValid = ((detLoc{1}(1,j) >= extents(1,1) & detLoc{1}(1,j) <= extents(1,2)) & (detLoc{1}(2,j) >= extents(2,1) & detLoc{1}(2,j) <= extents(2,2)));
+ checkiter = 0;
+ while ( ~bValid )
+ if ( checkiter > 200 )
+ error('Iteration limit exceeded for valid axon start point');
+ end
+ detParam{1}(3,j) = rand()*axons(detParam{1}(1,j)).len;
+ detLoc{1}(:,j) = GetAxonLoc(detParam{1}(:,j));
+ bValid = ((detLoc{1}(1,j) >= extents(1,1) & detLoc{1}(1,j) <= extents(1,2)) & (detLoc{1}(2,j) >= extents(2,1) & detLoc{1}(2,j) <= extents(2,2)));
+
+ checkiter = checkiter + 1;
+ end
+ end
+ else
+ detParam{1} = 2*pi*rand(1,numdetections);
+ end
+
+ detSpeed{1} = (args.speedRange(2)-args.speedRange(1))*rand(1,numdetections) + args.speedRange(1);
+ %detSpeed{1} = 5*ones(1,numdetections);
+
+ speedTarget = detSpeed{1};
+ if ( numaxons > 0 )
+ paramTarget = detParam{1}(2,:);
+ else
+ paramTarget = detParam{1};
+ end
+
+ speedTime = zeros(1,numdetections);
+ paramTime = zeros(1,numdetections);
+
+ [detFalse{1}, numfalsedet] = AddFalseDetections();
+
+ totaldetections(1,1) = 0;
+ totaldetections(2,1) = numdetections + numfalsedet;
+
+ g2l((totaldetections(1,1)+1):totaldetections(2,1),:) = [ones((numdetections+numfalsedet),1) cumsum(ones((numdetections+numfalsedet),1))];
+
+ bNew = boolean(zeros(1,numdetections));
+ for i=2:numframes
+ %Vary speed
+ speedTime(bNew) = 0;
+ speedTarget(bNew) = detSpeed{i-1}(bNew);
+ [speedTarget,speedTime,deltaSpeed] = UpdateVariation(speedTarget, speedTime, detSpeed{i-1}, args.speedVarProb, [-3 3], args.speedVarTime, 0);
+ speedTarget(speedTarget < args.speedRange(1)) = args.speedRange(1);
+ speedTarget(speedTarget > args.speedRange(2)) = args.speedRange(2);
+
+ %Vary lane or angle (parameters)
+ paramTime(bNew) = 0;
+ if ( numaxons > 0 )
+ paramTarget(bNew) = detParam{i-1}(2,bNew);
+ [paramTarget,paramTime,deltaParam] = UpdateVariation(paramTarget, paramTime, detParam{i-1}(2,:), args.paramVarProb, [-axonwidth axonwidth], args.paramVarTime, 1);
+ else
+ paramTarget(bNew) = detParam{i-1}(bNew);
+ [paramTarget,paramTime,deltaParam] = UpdateVariation(paramTarget, paramTime, detParam{i-1}, args.paramVarProb, [-(2*pi/8) (2*pi/8)], args.paramVarTime, 0);
+ end
+
+ [detLoc{i} detParam{i} detSpeed{i}] = MoveDetections(detLoc{i-1}, detParam{i-1}, detSpeed{i-1}+deltaSpeed, deltaParam, extents);
+ [detLoc{i} detParam{i} detSpeed{i} bNew] = AddDetections(detLoc{i}, detParam{i}, detSpeed{i}, extents);
+
+ [detFalse{i}, numfalsedet] = AddFalseDetections();
+
+ totaldetections(i+1,1) = totaldetections(i,1) + size(detLoc{i},2) + numfalsedet;
+ g2l((totaldetections(i,1)+1):totaldetections(i+1,1),:) = [i*ones((numdetections+numfalsedet),1) cumsum(ones((numdetections+numfalsedet),1))];
+ end
+
+ gConnectActual = gConnectActual(1:totaldetections(end),1:totaldetections(end));
+
+ for i=1:numframes
+ detLoc{i} = [detLoc{i} detFalse{i}];
+ end
+
+ dotradius = args.dotRadius;
+ tracksActual = {};
+ for i=1:numframes-1
+ for j=1:size(detLoc{i},2)
+ inidx = j + totaldetections(i);
+ outidx = find(gConnectActual(inidx,:),1);
+ if ( ~isempty(outidx) )
+ if ( gConnectActual(inidx,outidx) < 0 )
+ gConnectActual(inidx,outidx) = length(tracksActual)+1;
+ tracksActual{gConnectActual(inidx,outidx),1} = inidx;
+ end
+ trackidx = gConnectActual(inidx,outidx);
+ tracksActual{trackidx,1} = [tracksActual{trackidx} outidx];
+
+ nextidx = find(gConnectActual(outidx,:),1);
+ if ( ~isempty(nextidx) )
+ gConnectActual(outidx,nextidx) = trackidx;
+ end
+ end
+ end
+ end
+
+% detections = [];
+% for i=1:numframes
+% ddx = (ones(size(detLoc{i},2),1)*detLoc{i}(1,:) - ((ones(size(detLoc{i},2),1)*detLoc{i}(1,:))')).^2;
+% ddy = (ones(size(detLoc{i},2),1)*detLoc{i}(2,:) - ((ones(size(detLoc{i},2),1)*detLoc{i}(2,:))')).^2;
+% dd = sqrt(ddx+ddy);
+% dd(boolean(tril(ones(size(dd))))) = Inf;
+% [r,c] = find(dd < dotradius);
+% detections(i).actualmap = num2cell(setdiff([1:size(detLoc{i},2)]',union(r,c)));
+% detections(i).loc = detLoc{i}(:,[detections(i).actualmap{:}]);
+% while( ~isempty(r) )
+% mergepts = [r(1) c(1)];
+% r = r(2:end);
+% c = c(2:end);
+% bchkmerge = ismember([r c], mergepts);
+% if ( ~isempty(bchkmerge) )
+% bchkmerge = bchkmerge(:,1) | bchkmerge(:,2);
+% end
+% while ( any(bchkmerge) )
+% mergepts = unique([mergepts r(bchkmerge)' c(bchkmerge)']);
+% r = r(~bchkmerge);
+% c = c(~bchkmerge);
+% bchkmerge = ismember([r c], mergepts);
+% if ( ~isempty(bchkmerge) )
+% bchkmerge = bchkmerge(:,1) | bchkmerge(:,2);
+% end
+% end
+% centroid = mean(detLoc{i}(:,mergepts),2);
+% detections(i).actualmap{end+1} = mergepts;
+% detections(i).loc(:,end+1) = centroid;
+% end
+% end
+
+ toc;
+
+
+ cmap = hsv(1024);
+ trackColors = zeros(length(tracksActual),3);
+ for i=1:length(tracksActual)
+ trackColors(i,:) = cmap(round(rand()*1023)+1,:);
+ end
+
+ images = [];
+% for i=1:numframes
+% fakeim = zeros(ylen,xlen);
+% for j=1:size(detLoc{i},2)
+% grdx = round(detLoc{i}(1,j)-dotradius-1):round(detLoc{i}(1,j)+dotradius+1);
+% validx = (grdx >= 1 & grdx <= xlen);
+% grdx = grdx(validx);
+%
+% grdy = round(detLoc{i}(2,j)-dotradius-1):round(detLoc{i}(2,j)+dotradius+1);
+% validy = (grdy >= 1 & grdy <= ylen);
+% grdy = grdy(validy);
+%
+% [X,Y] = meshgrid(grdx, grdy);
+% %fakeim(grdy,grdx) = max(fakeim(grdy,grdx), normpdf(X,detLoc{i}(1,j),dotradius/3).*normpdf(Y,detLoc{i}(2,j),dotradius/3));
+% fakeim(grdy,grdx) = fakeim(grdy,grdx) + normpdf(X,detLoc{i}(1,j),dotradius/3).*normpdf(Y,detLoc{i}(2,j),dotradius/3);
+% end
+% cla;
+% images(:,:,i) = mat2gray(fakeim);
+% imshow(fakeim);colormap(gray);hold on;
+% for j=1:numaxons
+% axnpts = ppval(axons(j).pp, 0:axons(j).len+2);
+% plot(axnpts(1,:), axnpts(2,:), '-b');
+% end
+%
+% % for j=1:length(tracksActual)
+% % trkt = g2l(tracksActual{j},1);
+% % trkidx = g2l(tracksActual{j},2);
+% % idx = find(trkt == i,1);
+% % if ( ~isempty(idx) )
+% % loc = detLoc{i}(:,trkidx(idx));
+% % plot(loc(1), loc(2), 'o', 'Color',trackColors(j,:));
+% % text(loc(1)+4, loc(2)+4, num2str(j), 'Color',trackColors(j,:));
+% % end
+% % end
+% end
+
+ ds.axons = axons;
+ ds.loc = detLoc;
+ ds.param = detParam;
+ ds.speed = detSpeed;
+
+ ds.gConnectActual = gConnectActual;
+ ds.tracksActual = tracksActual;
+
+ ds.rgCC = {};
+ ds.rgDetect = {};
+ ds.rgMotionBlur = {};
+
+ tic;
+
+ for i=1:numframes
+ for j=1:size(detLoc{i},2)
+ ds.rgDetect{i}(j,:) = [round(detLoc{i}(1,j)) round(detLoc{i}(2,j)) ones(1,25)];
+
+ grdx = round(detLoc{i}(1,j)-dotradius-1):round(detLoc{i}(1,j)+dotradius+1);
+ validx = (grdx >= 1 & grdx <= xlen);
+ grdx = grdx(validx);
+
+ grdy = round(detLoc{i}(2,j)-dotradius-1):round(detLoc{i}(2,j)+dotradius+1);
+ validy = (grdy >= 1 & grdy <= ylen);
+ grdy = grdy(validy);
+
+ [X,Y] = meshgrid(grdx, grdy);
+ bFG = (((X-detLoc{i}(1,j)).^2 + (Y-detLoc{i}(2,j)).^2) <= dotradius.^2);
+ r = X(bFG);
+ c = Y(bFG);
+
+ ds.rgCC{i}{j} = [r,c];
+ if ( i < numframes )
+ ds.rgMotionBlur{i}{j} = [];
+ for nxtidx=1:size(detLoc{i+1},2)
+ if ( sqrt(sum((detLoc{i+1}(:,nxtidx)-detLoc{i}(:,j)).^2)) <= 20 )
+ ds.rgMotionBlur{i}{j}(1:2,end+1) = [totaldetections(i+1)+nxtidx;0];
+ end
+ end
+ end
+ end
+ end
+
+ toc;
+
+
+
+ % Helper functions
+ function [loc param speed bValid] = MoveDetections(inLoc, inParam, inSpeed, deltaParam, extents)
+ param = inParam;
+ speed = inSpeed;
+
+ if ( numaxons > 0 )
+ bRailParam = (deltaParam > abs(inSpeed));
+ deltaParam(bRailParam) = sign(deltaParam(bRailParam)).*inSpeed(bRailParam);
+
+ param(2,:) = inParam(2,:) + deltaParam;
+ arclen = sqrt(deltaParam.^2);
+ testlen = inSpeed;
+
+ inc = param(4,:) .* (inSpeed / 20);
+
+ loc = inLoc;
+ for k=1:size(param,2)
+ loc(:,k) = GetAxonLoc(param(:,k));
+ end
+
+ bIncArc = arclen < testlen;
+ oldLoc = loc;
+ while ( any(bIncArc) )
+ for k=1:length(bIncArc)
+ if ( ~bIncArc(k) )
+ continue;
+ end
+ param(3,k) = param(3,k) + inc(k);
+
+ loc(:,k) = GetAxonLoc(param(:,k));
+ newlen = sqrt((loc(1,k) - oldLoc(1,k)).^2 + (loc(2,k) - oldLoc(2,k)).^2);
+ arclen(k) = arclen(k) + newlen;
+ end
+
+ bIncArc = (arclen < testlen);
+ oldLoc = loc;
+ end
+ else
+ param = param + deltaParam;
+ loc = inLoc + (ones(2,1)*inSpeed).*[cos(param);sin(param)];
+ end
+
+ bValid = ((loc(1,:) >= extents(1,1) & loc(1,:) <= extents(1,2)) & (loc(2,:) >= extents(2,1) & loc(2,:) <= extents(2,2)));
+
+ oldidx = find(bValid);
+ newidx = 1:nnz(bValid);
+ linidx = sub2ind(size(gConnectActual), (oldidx + totaldetections(i-1)), (newidx + totaldetections(i)));
+ gConnectActual(linidx) = -1;
+
+ if ( numaxons > 0 )
+ param = param(:,bValid);
+ else
+ param = param(bValid);
+ end
+ speed = speed(bValid);
+ loc = loc(:,bValid);
+ end
+
+ function [loc param speed bNew] = AddDetections(inLoc, inParam, inSpeed, extents)
+ numnew = numdetections - size(inLoc,2);
+
+ loc = inLoc;
+ param = inParam;
+ speed = inSpeed;
+
+ edge = [2:(xlen-1) xlen*ones(1,ylen) (xlen-1):-1:2 ones(1,ylen); ones(1,xlen-2) 1:ylen ylen*ones(1,xlen-2) ylen:-1:1];
+
+ if ( numaxons > 0 )
+ newparams = [ceil((numaxons-eps)*rand(1,numnew) + eps); 2*axonwidth*rand(1,numnew)-axonwidth; zeros(1,numnew); 2*round(rand(1,numnew))-1];
+ for k=1:numnew
+ if ( newparams(4,k) < 0 )
+ newparams(3,k) = axons(newparams(1,k)).len;
+ end
+ end
+ param = [param newparams];
+ speed = [speed 10*rand(1,numnew)];
+ oldpidx = size(inLoc,2);
+ for k=1:numnew
+ param(3,oldpidx+k) = FindStart(param(:,oldpidx+k), extents);
+ loc = [loc GetAxonLoc(param(:,oldpidx+k))];
+ end
+ else
+ idxs = round((size(edge,2)-1)*rand(1,numnew)) + 1;
+ loc = [loc edge(:,idxs)];
+ param = [param 2*pi*rand(1,numnew)];
+ speed = [speed 10*rand(1,numnew)];
+ end
+
+ bNew = boolean([zeros(1,size(inLoc,2)) ones(1,numnew)]);
+ end
+
+ function [loc numfalse] = AddFalseDetections()
+ loc = [];
+ stdvfalse = args.falseDetections / 4;
+ numfalse = round(stdvfalse*randn(1,1) + args.falseDetections);
+
+ if ( args.falseDetections > 0 && numfalse )
+ if ( numaxons > 0 )
+ parms = [ceil((numaxons-eps)*rand(1,numfalse) + eps); 2*axonwidth*rand(1,numfalse)-axonwidth; rand(1,numfalse); 2*round(rand(1,numfalse))-1];
+ parms(3,:) = parms(3,:).*[axons(parms(1,:)).len];
+
+ for k=1:numfalse
+ loc(:,k) = GetAxonLoc(parms(:,k));
+ bOnScreen = ((loc(1,k) >= extents(1,1) & loc(1,k) <= extents(1,2)) & (loc(2,k) >= extents(2,1) & loc(2,k) <= extents(2,2)));
+ checkiter = 0;
+ while ( ~bOnScreen )
+ if ( checkiter > 200 )
+ error('Iteration limit exceeded for valid axon start point');
+ end
+ parms(3,k) = rand()*axons(parms(1,k)).len;
+ loc(:,k) = GetAxonLoc(parms(:,k));
+ bOnScreen = ((loc(1,k) >= extents(1,1) & loc(1,k) <= extents(1,2)) & (loc(2,k) >= extents(2,1) & loc(2,k) <= extents(2,2)));
+
+ checkiter = checkiter + 1;
+ end
+ end
+ else
+ loc = [extents(1,2);extents(2,2)]*ones(1,numfalse).*rand(2,numfalse);
+ end
+ end
+ end
+
+ function [target time delta] = UpdateVariation(inTarget, inTime, inCurrent, pVariation, varRange, timeRange, bStaticVariation)
+ delta = zeros(1,length(inCurrent));
+
+ target = inTarget;
+ time = max(inTime - 1, zeros(size(inTime)));
+
+ bUpdateVaried = ((inTarget ~= inCurrent) & (inTime > 0));
+ bStopVaried = ((inTarget ~= inCurrent) & (inTime == 0));
+
+ delta(bUpdateVaried) = (inTarget(bUpdateVaried) - inCurrent(bUpdateVaried)) ./ inTime(bUpdateVaried);
+ target(bStopVaried) = inCurrent(bStopVaried);
+
+ bNewVaried = (rand(1,length(inCurrent)) < pVariation);
+
+ time(bNewVaried) = floor((timeRange(2) - timeRange(1)+1-eps)*rand(1,nnz(bNewVaried))) + timeRange(1);
+ if ( bStaticVariation )
+ target(bNewVaried) = (varRange(2)-varRange(1))*rand(1,nnz(bNewVaried)) + varRange(1);
+ else
+ target(bNewVaried) = inCurrent(bNewVaried) + ((varRange(2)-varRange(1))*rand(1,nnz(bNewVaried)) + varRange(1));
+ end
+
+ delta(bNewVaried) = (target(bNewVaried) - inCurrent(bNewVaried)) ./ time(bNewVaried);
+ end
+
+ function arclen = FindStart(param, extents)
+ widthtol = 50;
+ arcrange = [-widthtol axons(param(1)).len+widthtol];
+ arclen = param(3);
+ tstloc = GetAxonLoc(param, arclen);
+ edgedist = [tstloc(1)-extents(1,1) extents(1,2)-tstloc(1) tstloc(2)-extents(2,1) extents(2,2)-tstloc(2)];
+ [mindist sideidx] = min(abs(edgedist));
+
+ if ( any(edgedist < 0) )
+ stepinc = param(4)*0.1;
+ else
+ stepinc = -param(4)*0.1;
+ end
+
+ oldstate = any(edgedist < 0);
+ while ( (mindist > 0.5) || any(edgedist < 0) )
+ arclen = arclen + stepinc;
+ tstloc = GetAxonLoc(param, arclen);
+ edgedist = [tstloc(1)-extents(1,1) extents(1,2)-tstloc(1) tstloc(2)-extents(2,1) extents(2,2)-tstloc(2)];
+ [mindist sideidx] = min(abs(edgedist(:)));
+
+ if ( (mindist > 0.5) && oldstate ~= any(edgedist < 0) )
+ stepinc = -(stepinc / 2);
+ end
+
+ oldstate = any(edgedist < 0);
+
+ if ( arclen < arcrange(1) || arclen > arcrange(2) )
+ break;
+ end
+ end
+
+ end
+
+ function [loc, tan, nrm] = GetAxonLoc(param, arclen)
+ if ( ~exist('arclen','var') )
+ arclen = param(3);
+ end
+ sloc = ppval(axons(param(1)).pp, arclen);
+ tan = ppval(axons(param(1)).dpp, arclen);
+ nrm = [tan(2); -tan(1)] / norm(tan);
+ loc = sloc+param(2)*nrm;
+ end
+
+ function outAxons = MakeAxons(extents, maxcurve)
+ outAxons = struct();
+ extrapts = floor((3-eps)*rand(1,numaxons));
+
+ edge = [2:(xlen-1) xlen*ones(1,ylen) (xlen-1):-1:2 ones(1,ylen); ones(1,xlen-2) 1:ylen ylen*ones(1,xlen-2) ylen:-1:1];
+ sides = [1 xlen xlen+ylen+1 2*xlen+ylen+1 length(edge)+1];
+
+ axsides = zeros(2,numaxons);
+ axsides(1,:) = floor((4-eps)*rand(1,numaxons));
+ axsides(2,:) = mod(axsides(1,:) + floor((3-eps)*rand(1,numaxons))+1, 4);
+ axsides = axsides + 1;
+
+ tstextents = extents;
+ tstextents(:,1) = tstextents(:,1) + axonwidth;
+ tstextents(:,2) = tstextents(:,2) - axonwidth;
+
+ for k=1:numaxons
+ edgeidx = [round((sides(axsides(1,k)+1)-sides(axsides(1,k)))*rand())+sides(axsides(1,k)) round((sides(axsides(2,k)+1)-sides(axsides(2,k)))*rand())+sides(axsides(2,k))];
+ pts = edge(:,edgeidx);
+
+ while ( norm(pts(:,2) - pts(:,1)) < 20.0 )
+ edgeidx = [round((sides(axsides(1,k)+1)-sides(axsides(1,k)))*rand())+sides(axsides(1,k)) round((sides(axsides(2,k)+1)-sides(axsides(2,k)))*rand())+sides(axsides(2,k))];
+ pts = edge(:,edgeidx);
+ end
+
+ if ( extrapts(k) > 0)
+ curveoffset = cumsum(ones(1,extrapts(k)))+0.2;
+ curverat = 0.6*rand(1,extrapts(k)) + curveoffset;
+ curverat = ((1/extrapts(k))*curverat + (1-(1/extrapts(k)))*curverat(1)) * (1/3);
+
+ curvyness = 2*maxcurve*rand(1,extrapts(k))-maxcurve;
+
+ curvetan = (pts(:,2)-pts(:,1));
+ curvenrm = [curvetan(2); -curvetan(1)] / norm(curvetan);
+ curvepts = curvetan*curverat + pts(:,1)*ones(1,extrapts(k)) + curvenrm*curvyness;
+
+ curvepts(1,:) = max(curvepts(1,:), tstextents(1,1)*ones(1,extrapts(k)));
+ curvepts(1,:) = min(curvepts(1,:), tstextents(1,2)*ones(1,extrapts(k)));
+ curvepts(2,:) = max(curvepts(2,:), tstextents(2,1)*ones(1,extrapts(k)));
+ curvepts(2,:) = min(curvepts(2,:), tstextents(2,2)*ones(1,extrapts(k)));
+
+ pts = [pts(:,1) curvepts pts(:,2)];
+ end
+
+ apprxlen = sqrt(sum(sum(((pts(:,2:end)-pts(:,1:end-1)).^2))));
+
+ outAxons(k).pp = pchip(0:(apprxlen/(length(pts)-1)):apprxlen, pts);
+ [breaks,coefs,l,kp,d] = unmkpp(outAxons(k).pp);
+ outAxons(k).dpp = mkpp(breaks,repmat(kp-1:-1:1,d*l,1).*coefs(:,1:kp-1),d);
+ outAxons(k).len = apprxlen;
+ end
+ end
+end
\ No newline at end of file
diff --git a/matlab/makePhantomKymo.m b/matlab/makePhantomKymo.m
new file mode 100644
index 0000000000000000000000000000000000000000..9848abbe30e732699945d80c03d3a316adb642e9
--- /dev/null
+++ b/matlab/makePhantomKymo.m
@@ -0,0 +1,147 @@
+function kymos = makePhantomKymo(ds, trackingStruct)
+ kymos = [];
+
+ cmap = [cool(512); winter(512)];
+ trackColorsTest = cmap(round(1023*rand(length(trackingStruct.trackList),1)+1),:);
+
+ l2g = zeros(length(ds.rgDetect)+1,1);
+ g2l = zeros(size(ds.gConnectActual,1),2);
+ for j=1:length(ds.rgDetect)
+ l2g(j+1) = l2g(j) + size(ds.rgDetect{j},1);
+ g2l((l2g(j)+1):l2g(j+1),:) = [j*ones(size(ds.rgDetect{j},1),1) cumsum(ones(size(ds.rgDetect{j},1),1))];
+ end
+
+ numFrames = length(ds.rgDetect);
+
+ lineThickness = 2.5;
+ ptThickness = 14;
+
+ for axIdx = 1:length(ds.axons)
+ axonLen = floor(ds.axons(axIdx).len);
+
+ actfig = figure();
+ dotfig = figure();
+
+ [allKymoPts, noiseKymoPts] = getAllKymoPts(ds, ds.axons(axIdx));
+
+ figure(actfig);hold on;
+ plot(noiseKymoPts(:,1), noiseKymoPts(:,2), '.', 'MarkerSize',ptThickness, 'Color',[0.25,0.25,0.25]);
+
+ figure(dotfig);hold on;
+ plot(noiseKymoPts(:,1), noiseKymoPts(:,2), '.', 'MarkerSize',ptThickness, 'Color',[0.25,0.25,0.25]);
+
+ plotActualTracks(actfig, ds, g2l, axIdx, '.k', lineThickness, ptThickness);
+ plotActualTracks(dotfig, ds, g2l, axIdx, '-k', lineThickness, ptThickness);
+
+ plotTrackingResults(actfig, ds, g2l, allKymoPts, trackingStruct, trackColorsTest, 5);
+ plotTrackingResults(dotfig, ds, g2l, allKymoPts, trackingStruct, trackColorsTest, 5);
+
+ figure(dotfig);hold on;
+ xlim([1 numFrames]);
+ ylim([0 axonLen]);
+
+ figure(actfig);hold on;
+ xlim([1 numFrames]);
+ ylim([0 axonLen]);
+ end
+end
+
+function [allKymoPts noiseKymoPts] = getAllKymoPts(ds, axon)
+ numFrames = length(ds.rgDetect);
+
+ axonLen = floor(axon.len);
+
+ axonPts = zeros(2,axonLen+1);
+ for i=1:(axonLen+1)
+ axonPts(:,i) = ppval(axon.pp, axon.len*((i-1)/axonLen));
+ end
+
+ allKymoPts = [];
+ noiseKymoPts = [];
+ for t=1:numFrames
+ numTotal = length(ds.loc{t}(1,:));
+ numActual = length(ds.param{t}(1,:));
+
+ % True Detections
+ allKymoPts{t} = ds.param{t}(3,:)';
+
+ % Axon location of false detections
+ noisePts = ds.loc{t}(:,(numActual+1):end)';
+ [dist, segparam, segidx] = distToPath(noisePts(:,1), noisePts(:,2), axonPts);
+
+ allArcLen = axon.len*((segidx-1)/axonLen) + segparam;
+
+ bOnAxon = (dist <= 8);
+ allArcLen(~bOnAxon) = nan;
+ noiseArcLen = allArcLen(bOnAxon);
+
+ allKymoPts{t} = [allKymoPts{t}; allArcLen];
+
+ noiseKymoPts = [noiseKymoPts; t*ones(length(noiseArcLen),1) noiseArcLen];
+ end
+end
+
+function plotActualTracks(fig, ds, g2l, axIdx, linespec, linewidth, ptwidth)
+ figure(fig);hold on;
+
+ numTracks = length(ds.tracksActual);
+
+ for i=1:numTracks
+ trkidx = g2l(ds.tracksActual{i},2)';
+ trkt = g2l(ds.tracksActual{i},1)';
+
+ if ( ds.param{trkt(1)}(1,trkidx(1)) ~= axIdx )
+ continue;
+ end
+
+ trkVal = zeros(1,length(trkt));
+ for j=1:length(trkt)
+ trkVal(j) = ds.param{trkt(j)}(3,trkidx(j));
+ end
+
+ plot(trkt, trkVal, linespec, 'LineWidth',linewidth, 'MarkerSize',ptwidth, 'LineSmoothing','on');
+ end
+ hold off;
+end
+
+function plotTrackingResults(fig, ds, g2l, allKymoPts, trackingStruct, trackColors, minlength)
+ figure(fig);hold on;
+
+ gAssignedActual = ds.gConnectActual > 0;
+ gAssignedTracking = trackingStruct.gAssigned;
+
+ numTracks = length(trackingStruct.trackList);
+
+ for i=1:numTracks
+ trkidx = g2l(trackingStruct.trackList{i},2)';
+ trkt = g2l(trackingStruct.trackList{i},1)';
+
+ if ( length(trkt) < minlength )
+ continue;
+ end
+
+ trkPts = zeros(size(trkt));
+ for j=1:length(trkt)
+ trkPts(j) = allKymoPts{trkt(j)}(trkidx(j));
+ end
+
+ plot(trkt, trkPts, '-', 'Color',trackColors(i,:), 'LineWidth',1);
+ end
+
+ bErr = full(any((gAssignedActual - gAssignedTracking),2));
+ comprows = Inf*ones(size(bErr));
+
+ comprows(bErr) = full(any(gAssignedActual(bErr,:),2) - any(gAssignedTracking(bErr,:),2));
+ edgeErrs = find(comprows == 0);
+
+ errt = g2l(edgeErrs,1)';
+ erridx = g2l(edgeErrs,2)';
+
+ errPts = zeros(size(erridx));
+ for j=1:length(errt)
+ errPts(j) = allKymoPts{errt(j)}(erridx(j));
+ end
+ plot(errt,errPts,'or', 'MarkerSize',6);
+
+ hold off;
+end
diff --git a/matlab/makePhantomMovie.m b/matlab/makePhantomMovie.m
new file mode 100644
index 0000000000000000000000000000000000000000..ce95a4dda908d56644339b16725c27f74961db56
--- /dev/null
+++ b/matlab/makePhantomMovie.m
@@ -0,0 +1,152 @@
+function makePhantomMovie(fname, ds, trackingStruct, extents, zoom, fps)
+ gAssignedActual = (ds.gConnectActual > 0);
+
+ l2g = zeros(length(ds.rgDetect)+1,1);
+ g2l = zeros(size(ds.gConnectActual,1),2);
+ for i=1:length(ds.rgDetect)
+ l2g(i+1) = l2g(i) + size(ds.rgDetect{i},1);
+ g2l((l2g(i)+1):l2g(i+1),:) = [i*ones(size(ds.rgDetect{i},1),1) cumsum(ones(size(ds.rgDetect{i},1),1))];
+ end
+
+ cmap = [cool(256); spring(256); summer(256); winter(256)];
+ trackColorsActual = cmap(round(1023*rand(length(ds.tracksActual),1)+1),:);
+ trackColorsTest = cmap(round(1023*rand(length(trackingStruct.trackList),1)+1),:);
+
+ zoomsize = [zoom(1,2)-zoom(1,1)+1 zoom(2,2)-zoom(2,1)+1];
+
+ axExtents = getAxonExtents(ds, extents, 20);
+
+ imsize = [axExtents(1,2)-axExtents(1,1) axExtents(2,2)-axExtents(2,1)];
+ if ( imsize(2) < 500 )
+ imsize = imsize * 2;
+ elseif ( imsize(2) < 600 )
+ imsize = round(imsize * 1.5);
+ imsize = imsize + mod(imsize(2),2);
+ end
+
+ zoomratio = imsize(2) / zoomsize(2);
+
+ imfig = figure();
+ zoomfig = figure();
+
+ numFrames = length(ds.rgDetect);
+
+ movfile = [];
+
+ try
+ movfile = avifile(fname, 'Compression','None', 'fps',fps);
+
+ for t=1:numFrames
+ drawImage(imfig, t, ds, extents, 1, 1);
+ overlayBox(imfig, zoom, 1.5);
+ xlim(axExtents(1,:));
+ ylim(axExtents(2,:));
+
+ drawImage(zoomfig, t, ds, extents, 0.7, 0);
+ overlayBox(zoomfig, [zoom(1,1)+0.25 zoom(1,2)-0.3; zoom(2,1)+0.375 zoom(2,2)-0.125], 3);
+ overlayTracks(zoomfig, t, ds, g2l, ds.tracksActual, trackColorsActual, 3);
+ xlim(zoom(1,:));
+ ylim(zoom(2,:));
+
+ figure(imfig);
+ truesize(imfig, [imsize(2) imsize(1)]);
+ imfrm = getframe(gca);
+
+ figure(zoomfig);
+ truesize(zoomfig, round(zoomratio*[zoomsize(2) zoomsize(1)]));
+ zoomfrm = getframe(gca);
+
+ if ( size(imfrm.cdata,2) ~= size(zoomfrm.cdata,2) )
+ newsz = min(size(imfrm.cdata,1), size(zoomfrm.cdata,1));
+ imfrm.cdata = imfrm.cdata(1:newsz,:,:);
+ zoomfrm.cdata = zoomfrm.cdata(1:newsz,:,:);
+ end
+
+ fullframe = imfrm;
+ fullframe.cdata = [fullframe.cdata zoomfrm.cdata];
+
+ movfile = addframe(movfile,fullframe);
+
+ clf(imfig);
+ clf(zoomfig);
+ end
+ catch excp
+ if ( ~isempty(movfile) )
+ movfile = close(movfile);
+ end
+ rethrow (excp);
+
+ end
+
+ movfile = close(movfile);
+end
+
+function drawImage(fig, t, ds, extents, maxalpha, bDrawCenterline)
+ dotradius = 4;
+
+ fakeim = zeros(extents(2,2),extents(1,2));
+ for j=1:size(ds.loc{t},2)
+ fakeim = drawGaussianPt(fakeim, ds.loc{t}(1,j), ds.loc{t}(2,j), 1, dotradius, dotradius);
+ end
+
+ figure(fig);
+ imshow(maxalpha*mat2gray(fakeim));colormap(gray);hold on;
+
+ if ( bDrawCenterline )
+ for j=1:length(ds.axons)
+ axnpts = ppval(ds.axons(j).pp, 0:ds.axons(j).len+2);
+ plot(axnpts(1,:), axnpts(2,:), '-b');
+ end
+ end
+end
+
+function overlayBox(fig, boxedges, thickness)
+ figure(fig);
+ rectangle('Position',[boxedges(:,1)' (boxedges(:,2)-boxedges(:,1))'], 'Curvature',[0 0], 'EdgeColor','r', 'LineWidth',thickness);
+% plot(boxedges(1,:), ones(1,2)*boxedges(2,1), '-r', 'LineWidth',thickness);
+% plot(boxedges(1,:), ones(1,2)*boxedges(2,2), '-r', 'LineWidth',thickness);
+% plot(ones(1,2)*boxedges(1,1), boxedges(2,:), '-r', 'LineWidth',thickness);
+% plot(ones(1,2)*boxedges(1,2), boxedges(2,:), '-r', 'LineWidth',thickness);
+end
+
+function overlayTracks(fig, t, ds, g2l, trackList, trackColors, traillength)
+ figure(fig);hold on;
+
+ numDets = length(ds.param{t});
+ noisePts = ds.loc{t}(:,numDets+1:end);
+ plot(noisePts(1,:), noisePts(2,:), 'xm', 'MarkerSize',12);
+
+ for j=1:length(trackList)
+ trkt = g2l(trackList{j},1);
+ trkidx = g2l(trackList{j},2);
+
+ idx = find((trkt >= t-traillength) & (trkt <= t));
+
+ if ( ~isempty(idx) )
+ loc = [];
+ for k=1:length(idx)
+ loc(:,k) = ds.loc{trkt(idx(k))}(:,trkidx(idx(k)));
+ end
+ plot(loc(1,:), loc(2,:), '-', 'Color',trackColors(j,:), 'LineWidth',1.5);
+ if ( trkt(idx(end)) == t )
+ plot(loc(1,end), loc(2,end), '.r', 'MarkerSize',12);
+ end
+ end
+ end
+end
+
+
+function extents = getAxonExtents(ds, inextents, backoff)
+
+ extents = [0 0; 0 0];
+
+ for axIdx=1:length(ds.axons)
+ axlens = [0:ds.axons(axIdx).len-1 ds.axons(axIdx).len];
+ axPts = ppval(ds.axons(axIdx).pp, axlens);
+
+ extents(1,1) = max(min(axPts(1,:))-backoff, inextents(1,1));
+ extents(1,2) = min(max(axPts(1,:))+backoff, inextents(1,2));
+ extents(2,1) = max(min(axPts(2,:))-backoff, inextents(2,1));
+ extents(2,2) = min(max(axPts(2,:))+backoff, inextents(2,2));
+ end
+end
\ No newline at end of file
diff --git a/matlab/runfakejaqaman.m b/matlab/runfakejaqaman.m
new file mode 100644
index 0000000000000000000000000000000000000000..acc1a903167a0ec9afb1c7cbc73dfdc81262e807
--- /dev/null
+++ b/matlab/runfakejaqaman.m
@@ -0,0 +1,58 @@
+outdir = 'C:\Users\mwinter\Documents\MATLAB\faketest_singleaxon_three';
+mkdir(outdir);
+numsdetlist = [10 20 40 60 80 100];
+falsedet = [0 10 20 40 60];
+numtrys = 20;
+
+trackertypes = {'_Jaqaman'};
+
+% axons = [];
+load(fullfile(outdir,'axons.mat'));
+for idx=1:length(numsdetlist)
+ numdet = numsdetlist(idx);
+
+ for idxfalse = 1:length(falsedet)
+ pctfalse = falsedet(idxfalse);
+ numfalse = round(numdet*(pctfalse/100));
+
+ for i=1:numtrys
+ fname_prefix = ['fakedata_' num2str(numdet,'%04d') 'det_' num2str(pctfalse,'%04d') 'false_' num2str(i,'%04d') 'try'];
+ segdatname = [fname_prefix '_seg.dat'];
+ mbdatname = [fname_prefix '_mb.dat'];
+
+ if ( ~exist(fullfile(outdir,[fname_prefix '.mat']),'file') )
+% [ds, axons] = makePhantom(512, 512, 'numAxons',1, 'axons',axons, 'numDetections',numdet, 'speedVarProb',0, 'paramVarProb',0, 'falseDetections', numfalse);
+ for k=1:length(ds.rgDetect)
+ for j=1:size(ds.rgDetect{k},1)
+ ds.rgDetect{k}(j,:) = [ds.loc{k}(1,j) ds.loc{k}(2,j) ones(1,25)];
+ end
+ end
+% save(fullfile(outdir,[fname_prefix '.mat']), 'ds');
+% writeSegData(fullfile(outdir,segdatname), ds.rgDetect, ds.rgCC, [1 512; 1 512]);
+% writeMotionBlurData(fullfile(outdir,mbdatname),
+% ds.rgMotionBlur);
+ else
+ load(fullfile(outdir,[fname_prefix '.mat']));
+
+ for k=1:length(ds.rgDetect)
+ for j=1:size(ds.rgDetect{k},1)
+ ds.rgDetect{k}(j,:) = [ds.loc{k}(1,j) ds.loc{k}(2,j) ones(1,25)];
+ end
+ end
+ end
+
+ for tkidx=1:length(trackertypes)
+ outdatname = [fname_prefix trackertypes{tkidx} '_tracks.mat'];
+
+ if ( ~exist(fullfile(outdir,outdatname), 'file') )
+ movieInfo = buildMovieInfo(ds);
+ scriptTrackGeneral
+ trackingStruct = buildTrackingStruct(ds, tracksFinal, 6);
+ save(fullfile(outdir,outdatname), 'trackingStruct', 'tracksFinal', 'movieInfo');
+ end
+ end
+ clear ds;
+ end
+ end
+end
+save(fullfile(outdir,'axons.mat'), 'axons');
\ No newline at end of file
diff --git a/matlab/runfaketrack.m b/matlab/runfaketrack.m
new file mode 100644
index 0000000000000000000000000000000000000000..fd483fef7dbf1db64a5f4de8ec7ea2af36f721bf
--- /dev/null
+++ b/matlab/runfaketrack.m
@@ -0,0 +1,100 @@
+outdir = 'faketest_singleaxon_largesets';
+mkdir(outdir);
+%numsdetlist = [10 25 50 75 100 125 150 200 250 300 350];
+%numsdetlist = [10 20 40 60 80 100];
+%falsedet = [0 10 20 40 60];
+numsdetlist = [10 20 40 60 80 100 120 140];
+falsedet = [80 100];
+numtrys = 20;
+
+trackertypes = {'_mta';'_sfa'};
+%trackertypes = ['_sfa'];
+
+numwindows = 10;
+windowprocessing = [];
+windidx = 1;
+prctime = [];
+giveuptime = 54000;
+
+% axons = [];
+pause on;
+load(fullfile(outdir,'axons.mat'));
+for idx=1:length(numsdetlist)
+ numdet = numsdetlist(idx);
+
+ for idxfalse = 1:length(falsedet)
+ pctfalse = falsedet(idxfalse);
+ numfalse = round(numdet*(pctfalse/100));
+
+ for i=1:numtrys
+ fname_prefix = ['fakedata_' num2str(numdet,'%04d') 'det_' num2str(pctfalse,'%04d') 'false_' num2str(i,'%04d') 'try'];
+ segdatname = [fname_prefix '_seg.dat'];
+ mbdatname = [fname_prefix '_mb.dat'];
+
+ if ( ~exist(fullfile(outdir,[fname_prefix '.mat']),'file') )
+ [ds, axons] = makePhantom(512, 512, 'numAxons',1, 'axons',axons, 'numDetections',numdet, 'speedVarProb',0, 'paramVarProb',0, 'falseDetections', numfalse);
+ for k=1:length(ds.rgDetect)
+ for j=1:size(ds.rgDetect{k},1)
+ ds.rgDetect{k}(j,:) = [ds.loc{k}(1,j) ds.loc{k}(2,j) ones(1,25)];
+ end
+ end
+ save(fullfile(outdir,[fname_prefix '.mat']), 'ds');
+ writeSegData(fullfile(outdir,segdatname), ds.rgDetect, ds.rgCC, [1 512; 1 512]);
+ writeMotionBlurData(fullfile(outdir,mbdatname), ds.rgMotionBlur);
+ else
+ load(fullfile(outdir,[fname_prefix '.mat']));
+
+ for k=1:length(ds.rgDetect)
+ for j=1:size(ds.rgDetect{k},1)
+ ds.rgDetect{k}(j,:) = [ds.loc{k}(1,j) ds.loc{k}(2,j) ones(1,25)];
+ end
+ end
+ end
+
+ for tkidx=1:length(trackertypes)
+ outdatname = [fname_prefix trackertypes{tkidx} '_tracks.txt'];
+ dbgdatname = [fname_prefix trackertypes{tkidx} '_dbginfo.txt'];
+ segdatname = [fname_prefix '_seg.dat'];
+ mbdatname = [fname_prefix '_mb.dat'];
+
+ if ( ~exist(fullfile(outdir,outdatname), 'file') )
+ while ( 1 )
+ if ( isempty(prctime) )
+ prctime = tic();
+ end
+ deltat = toc(prctime);
+
+ keepidx = boolean(ones(size(windowprocessing)));
+ for chkwidx=1:length(windowprocessing)
+ windowprocessing(chkwidx).time = windowprocessing(chkwidx).time + deltat;
+
+ if ( ~exist(windowprocessing(chkwidx).name,'file') && windowprocessing(chkwidx).time < giveuptime )
+ continue;
+ end
+
+ keepidx(chkwidx) = 0;
+ end
+
+ windowprocessing = windowprocessing(keepidx);
+ windidx = nnz(keepidx)+1;
+
+ prctime = tic();
+ if ( windidx > numwindows )
+ pause(30);
+ else
+ break;
+ end
+ end
+
+ system(['KymoTracker_v2' trackertypes{tkidx} '.exe "' fullfile(outdir,segdatname) '" "' fullfile(outdir,mbdatname) '" "' fullfile(outdir,outdatname) '" "' fullfile(outdir,dbgdatname) '" &']);
+
+ windowprocessing(windidx).name = fullfile(outdir,outdatname);
+ windowprocessing(windidx).time = 0;
+ windidx = windidx + 1;
+ end
+ end
+ clear ds;
+ end
+ end
+end
+save(fullfile(outdir,'axons.mat'), 'axons');
\ No newline at end of file
diff --git a/matlab/sfaTrack.m b/matlab/sfaTrack.m
new file mode 100644
index 0000000000000000000000000000000000000000..e14e6b40848235c28c9ee1a4214602b02ba0b4e1
--- /dev/null
+++ b/matlab/sfaTrack.m
@@ -0,0 +1,53 @@
+%
+% This function takes a cell array of detections and outputs a matrix of
+% tracks.
+%
+function tracks = sfaTrack(rgDetect)
+ numFrames = length(rgDetect);
+ tracks=zeros(size(rgDetect{1},1), numFrames-1);
+ for t=1:numFrames-1
+ N = size(rgDetect{t},1);
+ M = size(rgDetect{t+1},1);
+ C = zeros(N, M);
+ trackAssgn = [];
+ for i=1:N
+ if ( t > 1 )
+ lastAssgn = find(tracks(:,t-1) == i);
+ if ( isempty(lastAssgn) )
+ trackAssgn(i) = size(tracks,1)+1;
+ tracks(end+1,:) = zeros(1, numFrames-1);
+ else
+ trackAssgn(i) = lastAssgn;
+ end
+ else
+ trackAssgn = 1:N;
+ end
+
+ for j=1:M
+ C(i,j) = GetCost(rgDetect{t}(i,:), rgDetect{t+1}(j,:));
+ end
+ end
+
+ [assign, cost] = assignmentoptimal(C);
+
+ for i=1:length(assign)
+ if ( assign(i) > 0 )
+ tracks(trackAssgn(i),t) = assign(i);
+ end
+ end
+ end
+end
+
+function cost = GetCost(p1, p2)
+ costEpsilon = 1e-7;
+ maxDistanceDelta = 100;
+
+ cost = inf;
+
+ distErr = (p1(1)-p2(1))^2 + (p1(2)-p2(2))^2;
+ if ( distErr > maxDistanceDelta )
+ return;
+ end
+
+ cost = distErr;
+end
\ No newline at end of file
diff --git a/matlab/testpathdist.m b/matlab/testpathdist.m
new file mode 100644
index 0000000000000000000000000000000000000000..57dc648af8faa4bcc2c745e1d25087a5b12e7104
--- /dev/null
+++ b/matlab/testpathdist.m
@@ -0,0 +1,24 @@
+path = [-3, 0, 2, 6, 9;0, 0, 0, 0, 0];
+seglen = 5;
+for i=1:length(path)-1
+ path(2,i+1) = sqrt(seglen^2 - (path(1,i+1)-path(1,i))^2) + path(2,i);
+end
+
+gridx = -30:0.05:30;
+gridy = -10:0.05:45;
+[X,Y] = meshgrid(gridx, gridy);
+dtst = zeros(size(X));
+dist = zeros(size(X));
+
+% for i=1:size(X,1)
+% for j=1:size(X,2)
+% [dist(i,j) dtst(i,j)] = distToPath([X(i,j);Y(i,j)], path);
+% end
+% end
+
+[dist dtst] = distToPath(X,Y,path);
+
+ndist = dist ./ max(max(dist));
+
+figure;imagesc(gridx,gridy,dtst + (0.75*ndist));hold on;
+plot(path(1,:)', path(2,:)', '-k', 'LineWidth',2);
\ No newline at end of file