Random fixes and added extent creation code to work with latest tracker...

Random fixes and added extent creation code to work with latest tracker version/segmentation writer.

matlab/GoTrackMovies.m

@@ -100,7 +100,10 @@ for i=startidx:endidx
         results(i) = trackingStruct;
         save(fullfile(outputroot, [trackingStruct.fullname '_a' mkalphastr(trackingStruct.alpha) '.mat']), 'trackingStruct');
 
-        writeSegData(fullfile(outputroot, segdatname), results(i).rgDetect, results(i).rgCC);
+        iminfo = imfinfo(fullfile(trackingStruct.moviepath, trackingStruct.moviename));
+        imSize = [iminfo(1).Height iminfo(1).Width];
+        extents = calcTraceExtents(trackingStruct, imSize);
+        writeSegData(fullfile(outputroot, segdatname), results(i).rgDetect, results(i).rgCC, extents);
         %writeTraceData([outdatname '_trace.dat'], [jPath';iPath']);
         writeMotionBlurData(fullfile(outputroot, mbdatname), results(i).rgMotionBlur);
 

matlab/MakeBDNFMovieList.m

+function movielist = MakeBDNFMovieList(inputroot)
+    movielist = [];
+    
+    dlist=dir(fullfile(inputroot, 'HD*'));
+    dirnums = parseDirList(dlist, 'HD', 'BDNF MIDDLE');
+    [srt, dorder] = sort(dirnums);
+    for i=1:length(dlist)
+        if ( dirnums(dorder(i)) < 0 )
+            continue;
+        end
+        
+        ROOT = dlist(dorder(i)).name;
+        
+        flist=dir(fullfile(inputroot, ROOT, '*BDNF*.tif'));
+        fnums = parseDirList(flist, 'cell ', 'middle BDNF.tif');
+        [srt, forder] = sort(fnums);
+        for j=1:length(flist)
+            if ( fnums(forder(j)) < 0 )
+                continue;
+            end
+            
+            cellname = ['cell ' num2str(fnums(forder(j))) ' middle'];
+            moviename = flist(forder(j)).name;
+            moviepath = fullfile(inputroot, ROOT);
+            moviefilename = fullfile(moviepath, moviename);
+            tracefilename = fullfile(inputroot, 'OutputKymos', [ROOT '_' cellname '_path.mat']);
+            kymofilename = fullfile(inputroot, 'OutputKymos', [ROOT '_' cellname '_kymo.tif']);
+            if ( exist(moviefilename, 'file') && exist(tracefilename, 'file') && exist(kymofilename, 'file') )
+                movielist(end+1,1).root = [ROOT '_'];
+                movielist(end,1).cellname = cellname;
+                movielist(end,1).moviename = moviename;
+                movielist(end,1).moviepath = moviepath;
+                movielist(end,1).tracefilename = tracefilename;
+                movielist(end,1).kymofilename = kymofilename;
+                movielist(end,1).datasetnum = dirnums(dorder(i));
+            end
+        end
+    end
+end
+
+function nums = parseDirList(list, prefix, postfix)
+    nums = -ones(length(list),1);
+    for i=1:length(list)
+        stidx = 1;
+        endidx = length(list(i).name);
+        if ( exist('prefix','var') && ~isempty(prefix) )
+            stidx = strfind(lower(list(i).name), lower(prefix));
+            
+            if ( stidx ~= 1 )
+                continue;
+            end
+            
+            stidx = stidx + length(prefix);
+        end
+        
+        if ( exist('postfix','var') && ~isempty(postfix) )
+            endidx = strfind(lower(list(i).name), lower(postfix));
+            endidx = endidx - 1;
+        end
+        
+        if ( isempty(stidx) || isempty(endidx) || (stidx == endidx) )
+            continue;
+        end
+        
+        tkstr = strtrim(list(i).name(stidx:endidx));
+        tkidx = strfind(tkstr, ' ');
+        tkidx = [1 tkidx length(tkstr)];
+        
+        for j=1:length(tkidx)-1
+            num = str2num(tkstr(tkidx(j):tkidx(j+1)));
+            if ( ~isempty(num) )
+                nums(i) = num;
+                break;
+            end
+        end
+        
+    end
+end
+
+    
\ No newline at end of file

matlab/NewGoSegmentIt.m

+%function [rgDetect, rgMotionBlur, imKymo, LKymo, iPath, jPath, totalDetections] = GoSegmentIt(movieroot, moviefilename, tracefilename, kymofilename, alpha, extMaxH)
+function segResultStruct = NewGoSegmentIt(inResultStruct, tracefilename, kymofilename)
+
+segResultStruct = inResultStruct;
+
+rgDetect = [];
+rgCC = [];
+rgMotionBlur = [];
+LKymo = [];
+iPath = [];
+jPath = [];
+locNodeLUT = 0;
+globNodeLUT = [];
+
+alpha = inResultStruct.alpha;
+extMaxH = inResultStruct.extH;
+tracerange = inResultStruct.tracerange;
+
+bCheckMotionBlur = 1;
+% if ( nargout < 2 )
+%     bCheckMotionBlur = 0;
+% end 
+
+Movie_Name = fullfile(inResultStruct.moviepath, inResultStruct.moviename);
+load(tracefilename);
+imKymo=imread(kymofilename);
+
+if ( isempty(tracerange) )
+    tracerange = [1 length(iPath)];
+    segResultStruct.tracerange = tracerange;
+end
+
+segResultStruct.imKymo = imKymo;
+segResultStruct.LKymo = LKymo;
+segResultStruct.iPath = iPath;
+segResultStruct.jPath = jPath;
+
+nKymo=0*imKymo;
+hold off; imagesc(imKymo);hold on;colormap gray
+RADIUS_DILATE=12;
+RADIUS=10;
+info=imfinfo(Movie_Name);
+NFRAMES = length(info);
+
+sed = strel('disk',RADIUS_DILATE);
+se = strel('disk',RADIUS);
+
+rcFG={};
+rcPKS={};
+
+nbdsize = 5;
+nbdoffset = floor(nbdsize / 2);
+nbdsize = 2*nbdoffset+1;
+rgDetect = [];
+
+totalDetections = 0;
+
+for t=1: NFRAMES
+    im=(imread(Movie_Name,t));
+    im=mat2gray(im);
+%
+    h = fspecial('gaussian', 7, 5);
+    imfilt = imfilter(im, h, 'symmetric');
+    for j=1:50
+        imfilt = imfilter(imfilt, h, 'symmetric');
+    end
+    imhfreq = max((im - imfilt), zeros(size(im)));
+    medim = medfilt2(imhfreq,[4,4]);
+    %im = imhfreq;
+    
+    bw=logical(0*im);
+    idx=sub2ind(size(im),round(iPath(tracerange(1):tracerange(2))), round(jPath(tracerange(1):tracerange(2))));
+    bw(idx)=1;
+    bwd=imdilate(bw,sed);
+    bw=imdilate(bw,se);
+
+%
+    %level=alpha*graythresh(im(bwd));
+    level=alpha*graythresh(medim(bwd));
+    level = min(level,0.95);
+    %bwi=im2bw(im,level);
+    bwi=im2bw(medim,level);
+    bwi(~bw)=0;
+
+    %bwMax=imextendedmax(im,extMaxH,4)& bwi;
+    bwMax=imextendedmax(imhfreq,extMaxH,4)& bwi;
+    bwMax=bwmorph(bwMax,'shrink',Inf);
+
+    [rfg cfg]=find(bwi);
+    rcFG{t}=[rfg cfg];
+    %[rpk cpk]=find(bwMax & bwi);
+    %rcPKS{t}=[rpk cpk];
+
+    [L num]=bwlabel(bwi,4);
+    hold off;imagesc(imhfreq);hold on
+    %hold off;imagesc(im);hold on
+    [r c] = find(bwi);
+    plot(c,r,'.g');
+    [r c] = find(bwMax);
+    plot(c,r,'.r');
+
+% 
+% 
+%     bwp=bwboundaries(bwi);
+%     hold off;imagesc(im);hold on
+% %     if 1==t brighten(0.33),end
+%     for i=1:length(bwp)
+%         plot(bwp{i}(:,2),bwp{i}(:,1),'-r')
+%     end
+%     [r c]=find(bwi);
+%     plot(c,r,'.g')
+% 
+%     % plot all pts
+%     for n=1:num
+%         [r c]=find(L==n);
+%         for i=1:length(r)
+%             pKymo=LKymo(r(i),c(i));
+%             if 0==pKymo
+%                 dd=(jPath-c(i)).^2+(iPath-r(i)).^2;
+%                 [m pKymo]=min(dd);
+%             end
+%             plot(t,pKymo,'.b')
+%         end
+%     end
+    % plot maxima
+    j = 1;
+    rcPKS{t} = [];
+    for n=1:num
+        [r c]=find(L==n);
+        if length(r)>10,ccc='.m';
+        elseif length(r)>6,ccc='.r';
+        elseif length(r)>3,ccc='.g';
+        else ccc='.y';
+        end
+        ccc='.r';
+        for i=1:length(r)
+            if ~bwMax (r(i),c(i)),continue,end
+            pKymo=LKymo(r(i),c(i));
+            if 0==pKymo
+                dd=(jPath-c(i)).^2+(iPath-r(i)).^2;
+                [m pKymo]=min(dd);
+            end
+            plot(t,pKymo,ccc)
+            
+            yl = max(r(i)-nbdoffset,1);
+            xl = max(c(i)-nbdoffset,1);
+            yh = min(size(im,1),r(i)+nbdoffset);
+            xh = min(size(im,2),c(i)+nbdoffset);
+
+            yoffmin = yl - (r(i)-nbdoffset);
+            yoffmax = yh - (r(i)+nbdoffset);
+            xoffmin = xl - (c(i)-nbdoffset);
+            xoffmax = xh - (c(i)+nbdoffset);
+
+            peaknbd = -ones(2*nbdoffset+1, 2*nbdoffset+1);
+            peaknbd((yoffmin+1):(nbdsize+yoffmax),(xoffmin+1):(nbdsize+xoffmax)) = im(yl:yh,xl:xh);
+
+            rgDetect{t}(j,:) = [r(i) c(i) peaknbd(:)'];
+            rcPKS{t}(j,:) = [r(i) c(i)];
+            rgCC{t}{j} = [r c];
+            j = j + 1;
+            
+            totalDetections = totalDetections + 1;
+        end
+%         pix=find(L==n);
+%         if  ~any(bwMax(pix))
+%             rmean=round(mean(r));
+%             cmean=round(mean(c));
+%             bwMax (rmean,cmean)=1;
+%             pKymo=LKymo(rmean,cmean);
+%             if 0==pKymo
+%                 dd=(jPath-cmean).^2+(iPath-rmean).^2;
+%                 [m pKymo]=min(dd);
+%             end
+%             plot(t,pKymo,'xg')
+%         end
+    end
+
+    if ( t < NFRAMES )
+        locNodeLUT(t+1) = j+locNodeLUT(t)-1;
+    end
+
+    drawnow
+end
+
+if ( ~bCheckMotionBlur )
+    rgMotionBlur = [];
+    return;
+end
+
+lenpix = size(imKymo,1);
+if ( (totalDetections / (lenpix*NFRAMES)) > 0.20 )
+    disp(['Too many detections ' inResultStruct.moviename]);
+    rgDetect = [];
+    rgMotionBlur = [];
+    return;
+end
+
+t = 1;
+idx = 1;
+for i=1:totalDetections
+    globNodeLUT(i,:) = [t,idx];
+    idx = idx + 1;
+    if ( idx > size(rgDetect{t},1) )
+        idx = 1;
+        t = t + 1;
+    end
+end
+
+VMAX=20;
+rgMotionBlur = cell(1,length(rgDetect)-1);
+for t=1: NFRAMES-1
+    
+    for nxtfrm = 1:min(1,NFRAMES-t)
+        bwtFG=logical(0*im);
+        idx=sub2ind(size(im),rcFG{t}(:,1),rcFG{t}(:,2));
+        bwtFG(idx)=1;
+
+        bwt1FG=logical(0*im);
+        idx=sub2ind(size(im),rcFG{t+nxtfrm}(:,1),rcFG{t+nxtfrm}(:,2));
+        bwt1FG(idx)=1;
+
+        [L num]=bwlabel(bwtFG,8);
+        [L1 num1]=bwlabel(bwt1FG,8);
+%         figure();im=(imread(Movie_Name,t));im(im>stats(2))=stats(2);im=mat2gray(im);imshow(im);hold on;
+%         for i=1:size(rcPKS{t},1)
+%             [r c] = find(L==L(rcPKS{t}(i,1),rcPKS{t}(i,2)));
+%     %         imtst = zeros(size(im));
+%     %         imtst(tstpts) = 1;
+%     %         btst = bwboundaries(imtst);
+%     %         elfit = EllipseDirectFit([btst{1}(:,2),btst{1}(:,1)]);
+%             plot(c,r,'.r');
+%         end
+        %hold off;
+
+        stats = regionprops(L, 'MajorAxisLength', 'MinorAxisLength', 'Orientation', 'Centroid');
+    
+        for i=1:size(rcPKS{t},1)
+            curlabel = L(rcPKS{t}(i,1),rcPKS{t}(i,2));
+            [r c]=find(L==curlabel);
+            
+            if ( nxtfrm == 1 )
+                rgMotionBlur{t}{i} = [];
+            end
+            
+%             curvdir = [cosd(stats(curlabel).Orientation), -sind(stats(curlabel).Orientation)];
+%             plot([stats(curlabel).Centroid(1)-0.5*stats(curlabel).MajorAxisLength*curvdir(1) stats(curlabel).Centroid(1)+0.5*stats(curlabel).MajorAxisLength*curvdir(1)], [stats(curlabel).Centroid(2)-0.5*stats(curlabel).MajorAxisLength*curvdir(2) stats(curlabel).Centroid(2)+0.5*stats(curlabel).MajorAxisLength*curvdir(2)],'-b');
+            
+            for j=1:size(rcPKS{t+nxtfrm},1)
+                vdir = (rcPKS{t+nxtfrm}(j,:)-rcPKS{t}(i,:)) / nxtfrm;
+                vij=norm(vdir);
+                
+%                 if ( vij > 0 )
+%                     vdir = vdir / vij;
+%                     ndir = [-vdir(2),vdir(1)];
+%                 else
+%                     vdir = [1,0];
+%                     ndir = [0,1];
+%                 end
+%                 
+%                 dottanlengths = dot(ones(size(r,1),1)*vdir, ([r,c] - ones(size(r,1),1)*rcPKS{t}(i,:)),2);
+%                 dotnormlengths = dot(ones(size(r,1),1)*ndir, ([r,c] - ones(size(r,1),1)*rcPKS{t}(i,:)),2);
+%                 tlen = (max(dottanlengths) - min(dottanlengths));
+%                 nlen = (max(dotnormlengths) - min(dotnormlengths));
+%                 
+%                 vapprx = (tlen - nlen) / (150 / 650);
+                
+                if vij>VMAX,continue,end
+                
+                [r1 c1]=find(L1==L1(rcPKS{t+nxtfrm}(j,1),rcPKS{t+nxtfrm}(j,2)));
+                dmin=Inf;
+                for k=1:length(r1)
+                    d=sqrt( (r1(k)-r).^2 + (c1(k)-c).^2);
+                    dd=min(d);
+                    if dd<dmin
+                        dmin=dd;
+                    end
+                end
+
+                p1=LKymo(rcPKS{t}(i,1),rcPKS{t}(i,2));
+                if 0==p1
+                    dd=(jPath-rcPKS{t}(i,2)).^2+(iPath-rcPKS{t}(i,1)).^2;
+                    [m p1]=min(dd);
+                end
+
+                p2=LKymo(rcPKS{t+nxtfrm}(j,1),rcPKS{t+nxtfrm}(j,2));
+                if 0==p2
+                    dd=(jPath-rcPKS{t+nxtfrm}(j,2)).^2+(iPath-rcPKS{t+nxtfrm}(j,1)).^2;
+                    [m p2]=min(dd);
+                end
+                if dmin< (nxtfrm-(150/650))*vij+1
+%                     if ( (vij >= 1) && abs(dot(curvdir, vdir)) < 0.5 )
+%                         plot([t t+nxtfrm],[p1 p2],'-r')
+%                     else
+                        plot([t t+nxtfrm],[p1 p2],'-y')
+%                     end
+                    %rgMotionBlur{t}{i}(1:2,end+1) = [j;1];
+                else
+                    plot([t t+nxtfrm],[p1 p2],'-b')
+                    %rgMotionBlur{t}{i}(1:2,end+1) = [j;0];
+                end
+                %rgMotionBlur{t}{i}(:,end+1) = [locNodeLUT(t+nxtfrm)+j;dmin;curvdir';abs(dot(curvdir, vdir))];
+                rgMotionBlur{t}{i}(1:2,end+1) = [locNodeLUT(t+nxtfrm)+j;(dmin/nxtfrm)];
+
+                plot(t,p1,'.r')
+                plot(t+nxtfrm,p2,'.r')
+            end
+        end
+    
+    end
+    
+drawnow
+end
+
+segResultStruct.rgDetect = rgDetect;
+segResultStruct.rgCC = rgCC;
+segResultStruct.locNodeLUT = locNodeLUT;
+segResultStruct.globNodeLUT = globNodeLUT;
+segResultStruct.rgMotionBlur = rgMotionBlur;
+
+end
\ No newline at end of file

matlab/calcTraceExtents.m

+function extents = calcTraceExtents(trackingStruct, imSize)
+    radius = 10;
+    extents = round([min(trackingStruct.iPath)-radius max(trackingStruct.iPath)+radius; min(trackingStruct.jPath)-radius max(trackingStruct.jPath)+radius]);
+    
+    if ( extents(1,1) < 1 )
+        extents(1,1) = 1;
+    end
+    
+    if ( extents(2,1) < 1 )
+        extents(2,1) = 1;
+    end
+    
+    if ( extents(1,2) > imSize(1) )
+        extents(1,2) = imSize(1);
+    end
+    
+    if ( extents(2,2) > imSize(2) )
+        extents(2,2) = imSize(2);
+    end
+end
\ No newline at end of file

matlab/mkalphastr.m

+function str = mkalphastr(alpha)
+    str = num2str(alpha, '%3.2f');
+end
\ No newline at end of file

src/tracker.cpp

@@ -280,6 +280,8 @@ int main(int argc, char* argv[])
 			BuildBestPaths(bestOutEdges, t, iLookback);
 		}
 
+		printf("t = %d, %d detections\n", t, rgDetectLengths[t]);
+
 		for ( int destPtIdx=0; destPtIdx < rgDetectLengths[t+1]; ++destPtIdx)
 		{
 			int nextGIdx = GetGlobalIdx(t+1, destPtIdx);