blobTrackopencv.m

% Syntax 
% blobTrackopencv('video',g,'maxFrame',1);
%
%  test1 = blobTrackopencv('video',g,'show',0,'mask',mask,'NumWorkers',12);
% 
%  % checks directly the results
%  test1.montagePatch('sidelength',128,'frames',1:1000:30000);
function [track,o] = blobTrackopencv(varargin);
o =mbparse.output();

p = inputParser();
p.addParamValue('video',[]);
p.addParamValue('videoFile',[]);
p.addParamValue('frames',[]);
p.addParamValue('show',false);
p.addParamValue('nfig','new');

p.addParamValue('useMask',true);
p.addParamValue('useBackground',false);

p.addParamValue('updateBackground',false);
p.addParamValue('updateBackgroundInterval',20);
p.addParamValue('updateBackgroundExtent',100);
p.addParamValue('useEnhancement',false);
p.addParamValue('NumWorkers',1);
p.addParamValue('opencv',false);

p.addParamValue('firstFrame',[]);
p.addParamValue('maxFrame',[]);

p.addParamValue('mask',[]);

p.addParamValue('binarizationThreshold',0.9);
p.addParamValue('erosion',8);

p.addParamValue('intermediateSaveInterval',[]);
p.addParamValue('intermediateSaveFile','blobTracktemp.mat');



p.parse(varargin{:});
q = p.Results;


if ~isempty(q.frames)
    if length(q.frames)==1
       o.echo('Evaluation of a single frame (%d)',q.frames(1)); 
       q.firstFrame = q.frames(1);
       q.maxFrame   = q.frames(1);
        
    end   
end





if ~isempty(q.videoFile)
    % passed as string
    videoFile = q.videoFile;
    video     = [];
    
    
else
    % passed as video
    videoFile = q.video.fileLocation;
    video     = q.video;
end

capTest = cv.VideoCapture(videoFile);
N = capTest.get('FrameCount');

u = capTest.read();
 
 if isempty(u)
    o.echo('Empty frame');
    return
 end
 
delete(capTest);

convertedMask = cast(q.mask,class(u));

cap = cv.VideoCapture(videoFile);

% N = video.NumberOfFrames();

if isempty(q.maxFrame)
    maxNum = N;
else
    maxNum = q.maxFrame;
end


centers = zeros(N,2);
areas   = zeros(N,1);
failedFrames = zeros(N,1);
counter = 1;

if q.show
thisFigure = mbgraph.figure(q.nfig);
end

eroder = strel('diamond',q.erosion);
eroderMatrix = eroder.getnhood;


maxIntensity = 255;

if isempty(q.firstFrame)
    firstFrame = 1;
else
    firstFrame = q.firstFrame;
    cap.set('PosFrames', q.firstFrame);
end

o.echo('Reading  frames %d through %d',{firstFrame,maxNum});

check  = mbparse.multicore.checkPool('matlab_workers',q.NumWorkers);


tic;


if check.NumWorkersFinal==1
    for i=firstFrame:maxNum;
        
        u = cap.read();
        
        if isempty(u)
            o.echo('Empty frame %d',i);
             failedFrames(i) = 1;
            continue;
        end
        u = rgb2gray(u);
        
        % inverse
        uinv = imcomplement(u);
        uinv = uinv.*convertedMask;
        
        % erosion
        ueroded  = cv.erode(uinv,'Element',eroderMatrix);
        
        % binarization
        % ubinarized =mbimg.binarize(mbimg.normalize(ueroded),q.binarizationThreshold);
        ubinarized = cv.threshold(ueroded,q.binarizationThreshold*maxIntensity,'Method','Binary');
        
        
        
        % segmentation
        
        [myCenter,usegmented] =  extracMouseLocal(ubinarized);
        
        if isempty(usegmented)
            failedFrames(i) = 1;
            continue;
        end
        
        centers(i,:) = myCenter;
        
        
        if q.show
            sm = 2;
            sn = 2;
            sc = 1;
            subplot(sm,sn,sc);sc = sc+1;
            imshow(u,[]);
            h = mbimg.showTile('centerhv',myCenter,'width',60,'height',60);
            title('estimated location on original');
            
            subplot(sm,sn,sc);sc = sc+1;
            imshow(ueroded,[]);
            title('Erosion');
            
            subplot(sm,sn,sc);sc = sc+1;
            imshow(ubinarized,[]);
            title('binarization');
            
            subplot(sm,sn,sc);sc = sc+1;
            imshow(usegmented,[]);
            title('segmented');
            
        end
        
    end
    
else
    
    disp('Task parallelization.');
    myFrames = firstFrame:maxNum;
    frameGroups = mbparse.multicore.distributions.array2consecutiveArrays(myFrames,check.NumWorkersFinal);
    
    for i=1:check.NumWorkersFinal
        %capArray{i} = cv.VideoCapture(videoFile);
        firstFrames(i) = frameGroups{i}(1);
        %capArray{i}.set('PosFrames',frameGroups{i}(1));
    end
    
    thresholdOutOfStructure = q.binarizationThreshold;
    
    spmd 
        %testParallelSimple(labindex,frameGroups{labindex},capArray(labindex),....
         %   convertedMask,thresholdOutOfStructure);
        centersLab = testParallel(labindex,frameGroups{labindex},firstFrames,videoFile,....
              convertedMask,thresholdOutOfStructure,eroderMatrix,maxIntensity);
    end
    
    % reconstitute matrix
    for i=1:length(centersLab)
        myChunk = centersLab{i};
        centers(frameGroups{i},:) = myChunk(frameGroups{i},:) ;
    end
    
end
toc



o.assignprop('centers',fliplr(centers));
o.assignprop('areas',areas);
o.assignprop('failedFrames',failedFrames);

if ~isempty(video)
    
    
    track = mbvid.videotrack.patchtrack('video',video,'data',fliplr(centers));
else
    o.echo('No video available! Use mbvid.videotrack.patchtrack on your own!');
    track = [];
end


o.ok = true;

delete(cap);


function [Centroid,bwmouse] =  extracMouseLocal(bw)

cc = bwconncomp(bw);


% check the size of secon componet
indexSecondBest = 1;

if isempty(cc.PixelIdxList)
    Centroid = [0,0];
    bwmouse  = [];
    return
    
end

for i=1:length(cc.PixelIdxList)
   myArea(i) = length(cc.PixelIdxList{i}); 
end

[nestArea,orderedIndices] = sort(myArea,'descend');

indexSecondBest = orderedIndices(1);


mousePixelIndices = cc.PixelIdxList{indexSecondBest};
N = length(cc.PixelIdxList{indexSecondBest});

io = bw(cc.PixelIdxList{indexSecondBest});
bwmouse = false(size(bw));
bwmouse(mousePixelIndices) = true;

s  = regionprops(bwmouse, 'centroid');
Centroid = s.Centroid;
% o.assignprop('Centroid',s.Centroid);
% o.assignprop('Area',N);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% PArallelization
%
function testParallelSimple(myIndex,myFrames,cap,convertedMask,binarizationThreshold);
o = mbparse.output();
o.echo('I am Matlab worker %d, starting with frame %d to %d',{myIndex,myFrames(1),myFrames(end)});
framesInTask = myFrames;

for i=framesInTask;
    
    u = cap.read();
    
end

function centers = testParallel(myIndex,myFrames,firstFrames,videoFile,convertedMask,binarizationThreshold,eroderMatrix,maxIntensity);

o = mbparse.output();
o.leading = sprintf('worker %d',myIndex);
o.echo(' frames %d through %d',{myFrames(1),myFrames(end)});
framesInTask = myFrames;

cap = cv.VideoCapture(videoFile);

%o.echo(' reading first frame: %d ',cap.get('PosFrames'));
cap.set('PosFrames',firstFrames(myIndex));
%o.echo(' reading again first frame: %d ',cap.get('PosFrames'));

if cap.get('PosFrames') ~= firstFrames(myIndex)
    o.echo(' First frame actually read in this worker: %d (should be %d)',{cap.get('PosFrames'),firstFrames(myIndex)});
 
    o.echo(' This is a vicious error: empty pointer in worker %d',myIndex);
    o.echo(' ... switching to manual upwind of cv.VideoCapture (pointers in this system do not accept set("PosFrames") ');
    delete(cap);
    cap = cv.VideoCapture(videoFile);
    for i=1:firstFrames(myIndex)-1;
        buffer = cap.read();
    end
    
end

if cap.get('PosFrames') ~= firstFrames(myIndex)
    o.echo(' First frame actually read in this worker: %d (should be %d)',{cap.get('PosFrames'),firstFrames(myIndex)});
    o.echo(' FATAL ERROR:  worker %d could not recover indetermination',myIndex);
    centers = NaN*ones(max(framesInTask),2);
    delete(cap);
    return
end

centers = zeros(max(framesInTask),2);

counter = 0;
for i=framesInTask;
    counter = counter+1;
    u = cap.read();

    if isempty(u)
        o.echo(' Empty frame (counter: %d)',counter);
        return
    end 
     
u = rgb2gray(u);
%     
%     % inverse
     uinv = imcomplement(u);
     uinv = uinv.*convertedMask;
%     
%     % erosion
     ueroded  = cv.erode(uinv,'Element',eroderMatrix);
%     
%     % binarization
%     % ubinarized =mbimg.binarize(mbimg.normalize(ueroded),q.binarizationThreshold);
     ubinarized = cv.threshold(ueroded,binarizationThreshold*maxIntensity,'Method','Binary');
%    
%     % segmentation
%     
     [myCenter,usegmented] =  extracMouseLocal(ubinarized); 
%     
     if isempty(usegmented)
         failedFrames(i) = 1;
         continue;
     end
%     
      centers(i,:) = myCenter;
    
 
end
o.echo('   Matlab worker %d ended!',{myIndex});
delete(cap);



function centers = singleTask(framesInTask,cap,convertedMask,eroderMatrix,maxIntensity)

for i=framesInTask;
    
    u = cap.read();
    u = rgb2gray(u);
    
    % inverse
    uinv = imcomplement(u);
    uinv = uinv.*convertedMask;
    
    % erosion
    ueroded  = cv.erode(uinv,'Element',eroderMatrix);
    
    % binarization
    % ubinarized =mbimg.binarize(mbimg.normalize(ueroded),q.binarizationThreshold);
    ubinarized = cv.threshold(ueroded,q.binarizationThreshold*maxIntensity,'Method','Binary');
   
    % segmentation
    
    [myCenter,usegmented] =  extracMouseLocal(ubinarized); 
    
    if isempty(usegmented)
        failedFrames(i) = 1;
        continue;
    end
    
    centers(i,:) = myCenter;
    
 
end