Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Fixed bug in GetRandFrames... #167

Open
wants to merge 4 commits into
base: master
Choose a base branch
from
Open

Fixed bug in GetRandFrames... #167

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
8c0a035
Fixed bug in GetRandFrames when using a specific experiment to get ta…
Mar 19, 2019
acacef8
Merge branch 'master' into alignAcrossPlanes
May 17, 2019
131a004
Fixes bug during interpolation of planes (didn't work of number of fr…
May 17, 2019
2fb6320
Fixed bug so that frame alignment works without specifying the specif…
May 23, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
125 changes: 50 additions & 75 deletions preRegistration/GetRandFrames.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,10 +12,6 @@
targetImage = getOr(ops, {'targetImage'}, []);

ntifs = sum(cellfun(@(x) numel(x), fs)); % total number of tiff files
nfmax = max(1, round(ops.NimgFirstRegistration/ntifs)); % number of tiffs to take per file
if nfmax>=2000
nfmax = 1999;
end

% size of images
nbytes = fs{1}(1).bytes;
Expand All @@ -29,96 +25,75 @@
indx = 0;
IMG = zeros(Ly, Lx, nplanes, ops.NimgFirstRegistration, 'single');

% grab random frames from target experiment
if ~isempty(targetImage)
k = targetImage(1);
% compute number of frames in tiff
if abs(nbytes - fs{k}(1).bytes)>1e3
nbytes = fs{k}(1).bytes;
nFr = nFrames(fs{k}(1).name);
end
startTiff = floor(length(fs{k})/2 - ...
if ~isempty(targetImage) % grab frames around half time during target experiment
expInds = targetImage(1);
startTiff = floor(length(fs{expInds})/2 - ...
nplanes*nchannels*ops.NimgFirstRegistration/2/nFr);
iplane = mod(((startTiff-1)*nFr)/nchannels-1, nplanes) + 1;
j = startTiff;
while indx < ops.NimgFirstRegistration
offset = 0;
if j == 1
offset = nchannels*nplanes;
tiffInds{expInds} = startTiff:length(fs{expInds});
iplaneK = 1 - (startTiff-1)*nFr; % to calculate index of first frame belonging to plane 1 in startTiff
nfmax = floor(nFr/nplanes/nchannels)-1;
else
expInds = 1:length(ops.SubDirs);
tiffInds = cell(1, length(ops.SubDirs));
for k = 1:length(ops.SubDirs)
tiffInds{k} = 1:length(fs{k});
end
iplaneK = 1;
nfmax = max(1, round(ops.NimgFirstRegistration/ntifs)); % number of tiffs to take per file
if nfmax>=2000
nfmax = 1999;
end
end

for k = expInds
iplane0 = iplaneK;
nchannels = getOr(ops, {'nchannels'}, 1);
if ~isempty(ops.expts)
if ismember(ops.expts(k), getOr(ops, 'expred', []))
nchannels = ops.nchannels_red;
end
end

for j = tiffInds{k}
% compute number of frames in tiff if size different from previous
if abs(nbytes - fs{k}(j).bytes)>1e3
nbytes = fs{k}(j).bytes;
nFr = nFrames(fs{k}(j).name);
end
numFr = min(ops.NimgFirstRegistration-indx, nFr/nchannels/nplanes);
if j==1
offset = nchannels*nplanes;
else
offset = 0;
end
% check if there are enough frames in the tiff
if nFr<(nchannels*nplanes*nfmax+offset)
continue;
end

iplane0 = mod(iplane0-1, nplanes) + 1;
% load red channel if red_align == 1
if red_align
ichanset = [offset + nchannels*(nplanes-iplane) + [rchannel;...
nchannels*nplanes*numFr]; nchannels];
ichanset = [offset + nchannels*(iplane0-1) + [rchannel;...
nchannels*nplanes*nfmax]; nchannels];
else
ichanset = [offset + nchannels*(nplanes-iplane) + [ichannel;...
nchannels*nplanes*numFr]; nchannels];
ichanset = [offset + nchannels*(iplane0-1) + [ichannel;...
nchannels*nplanes*nfmax]; nchannels];
end
iplane = mod(iplane + nFr/nchannels - 1, nplanes) + 1;
data = loadFramesBuff(fs{k}(j).name, ichanset(1), ichanset(2), ichanset(3));
iplane0 = iplane0 - nFr/nchannels;
data = loadFramesBuff(fs{k}(j).name, ichanset(1),ichanset(2), ichanset(3));
data = reshape(data, Ly, Lx, nplanes, []);

IMG(:,:,:,indx+(1:size(data,4))) = data;
indx = indx + size(data,4);
j = j + 1;
end
% grab frames from all files
else
for k = 1:length(ops.SubDirs)
iplane0 = 1;
nchannels = ops.nchannels;
if ~isempty(ops.expts)
if ismember(ops.expts(k), getOr(ops, 'expred', []))
nchannels = ops.nchannels_red;
end
end

for j = 1:length(fs{k})
% compute number of frames in tiff if size different from previous
if abs(nbytes - fs{k}(j).bytes)>1e3
nbytes = fs{k}(j).bytes;
nFr = nFrames(fs{k}(j).name);
end
if j==1
offset = nchannels*nplanes;
else
offset = 0;
end
% check if there are enough frames in the tiff
if nFr<(nchannels*nplanes*nfmax+offset)
continue;
end

iplane0 = mod(iplane0-1, nplanes) + 1;
% load red channel if red_align == 1
if red_align
ichanset = [offset + nchannels*(iplane0-1) + [rchannel;...
nchannels*nplanes*nfmax]; nchannels];
else
ichanset = [offset + nchannels*(iplane0-1) + [ichannel;...
nchannels*nplanes*nfmax]; nchannels];
end
iplane0 = iplane0 - nFr/nchannels;
data = loadFramesBuff(fs{k}(j).name, ichanset(1),ichanset(2), ichanset(3));
data = reshape(data, Ly, Lx, nplanes, []);

IMG(:,:,:,indx+(1:size(data,4))) = data;
indx = indx + size(data,4);

if indx>ops.NimgFirstRegistration
break;
end
end

if indx>ops.NimgFirstRegistration
break;
end
end

if indx>ops.NimgFirstRegistration
break;
end
end

IMG(:,:,:,(1+indx):end) = [];
48 changes: 30 additions & 18 deletions registration/registerBlocksAcrossPlanes.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -226,14 +226,15 @@
end

% align frames and write movie
t2 = 0;
t2 = 0; % number of all frames of already aligned experiments (divisible by nplanes)
% if two consecutive files have as many bytes, they have as many frames
nbytes = 0;
xyValid = true(ops.Ly, ops.Lx);
for k = 1:length(fs)
dataPrev = zeros(ops.Ly, ops.Lx, nplanes, 'int16');
dataPrev = zeros(ops.Ly, ops.Lx, nplanes, 'int16'); % for last nplanes frames of previous tiff file
iplane0 = 1:1:ops.nplanes;
t1 = 0;
t1 = 0; % number of all frames of already aligned tiffs in current experiment
tPlanes = zeros(1,indInterpolate(end)) + t2/nplanes; % index of imaging cycle of last aligned frame for each plane
for j = 1:length(fs{k})
if abs(nbytes - fs{k}(j).bytes)>1e3
nbytes = fs{k}(j).bytes;
Expand All @@ -257,6 +258,7 @@
end

% align frames according to determined registration offsets
% (without interpolating across planes)
[dreg, xyValid] = BlockRegMovie(data, ops, ...
dsall(t1+t2+(1:size(data,3)),:,:), xyValid);

Expand All @@ -273,56 +275,66 @@
if interpolateAcrossPlanes && ~isempty(RegFileBinLocation)
% use frames of plane that match best to target image of
% current plane and average across similar planes
dataNext = [];
for i = indInterpolate
ind1 = iplane0(planesToInterpolate(i)) : nplanes : size(data,3);
sh = shifts(ceil((t1+t2)/nplanes + (1:length(ind1))));
uniqueShifts = unique(sh);
dataNext = []; % first nplanes frames of following tiff if needed for interpolation
for i = indInterpolate % go through all planes (target images) that occur throughout all experiments
ind1 = find(iplane0 == planesToInterpolate(i)); % find indices of all frames belonging to current target frame
ind1 = ind1 : nplanes : size(dreg,3);
sh = shifts(tPlanes(i) + (1:length(ind1))); % find shifts in z for this target image in current tiff
tPlanes(i) = tPlanes(i) + length(ind1);
% sh = shifts(ceil((t1+t2)/nplanes + (1:length(ind1))));
uniqueShifts = unique(sh); % find planes covering all shifts
dwrite = zeros(size(dreg,1),size(dreg,2),length(ind1));
for s = uniqueShifts
planesInv = find(~isnan(profiles(:,i-s)))';
planesInv = find(~isnan(profiles(:,i-s)))'; % find planes used for interpolation
diffs = planesInv - i;
for pl = 1:length(planesInv)
ind2 = ind1 + diffs(pl); % smallest possible value: -nplanes+1,
ind2 = ind1 + diffs(pl); % find frames of current interpolation plane in tiff
% smallest possible value: -nplanes+1,
% largest possible values: size(dreg,3)+nplanes
ind2(sh~=s) = [];
ind3 = sh==s; % of those frames only use those where shift matches current shift level
ind2(~ind3) = [];
ind3 = find(ind3);
weight = profiles(planesInv(pl),i-s);
if ind2(1)<1 % frame is part of previously loaded tiff file (can happen if frames per tiff are not mupltiple of planes)
dwrite(:,:,1) = dwrite(:,:,1) + ...
weight .* double(dataPrev(:,:,end-ind2(1)));
weight .* double(dataPrev(:,:,end+ind2(1)));
ind2(1) = [];
ind3(1) = [];
end
if ind2(end)>size(dreg,3) % frame is part of upcoming tiff file
if j<length(fs{k})
if isempty(dataNext) % load first few frames of next tiff and register them
ichanset = [ichannel; nplanes; ops.nchannels];
data = loadFramesBuff(fs{k}(j+1).name, ...
dataN = loadFramesBuff(fs{k}(j+1).name, ...
ichanset(1), ichanset(2), ichanset(3));
if BiDiPhase
yrange = 2:2:Ly;
if BiDiPhase>0
data(yrange, (1+BiDiPhase):Lx,:) = data(yrange, 1:(Lx-BiDiPhase),:);
dataN(yrange, (1+BiDiPhase):Lx,:) = dataN(yrange, 1:(Lx-BiDiPhase),:);
else
data(yrange, 1:Lx+BiDiPhase,:) = data(yrange, 1-BiDiPhase:Lx,:);
dataN(yrange, 1:Lx+BiDiPhase,:) = dataN(yrange, 1-BiDiPhase:Lx,:);
end
end
dataNext = register_movie(data, ops, ...
dataNext = register_movie(dataN, ops, ...
dsall((t1+t2) + size(dreg,3) ...
+ (1:nplanes),:));
end
dwrite(:,:,end) = dwrite(:,:,end) + weight .*...
double(dataNext(:,:,ind2(end)-size(dreg,3)));
end
ind2(end) = [];
ind3(end) = [];
end
dwrite(:,:,ceil(ind2/nplanes)) = dwrite(:,:,ceil(ind2/nplanes)) + ...
dwrite(:,:,ind3) = dwrite(:,:,ind3) + ...
weight .* double(dreg(:,:,ind2));
% dwrite(:,:,ceil(ind2/nplanes)) = dwrite(:,:,ceil(ind2/nplanes)) + ...
% weight .* double(dreg(:,:,ind2));
end
end
fwrite(fidIntpol{ops.planesToProcess(planesToInterpolate(i))}, dwrite, class(data));
end
end
dataPrev = dreg(:,:,end-nplanes+1:end);
dataPrev = dreg(:,:,max(1,end-nplanes+1):end);
t1 = t1 + size(dreg,3);
iplane0 = iplane0 - nFr/ops.nchannels;
end
Expand Down