diff --git a/preRegistration/GetRandFrames.m b/preRegistration/GetRandFrames.m index dc7adfad..f9a8116e 100644 --- a/preRegistration/GetRandFrames.m +++ b/preRegistration/GetRandFrames.m @@ -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; @@ -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) = []; diff --git a/registration/registerBlocksAcrossPlanes.m b/registration/registerBlocksAcrossPlanes.m index e703f61f..1be6098c 100644 --- a/registration/registerBlocksAcrossPlanes.m +++ b/registration/registerBlocksAcrossPlanes.m @@ -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; @@ -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); @@ -273,40 +275,47 @@ 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 j0 - 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 @@ -314,15 +323,18 @@ 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