function M = find_xyz_roi( s, volume )
% find_xyz_roi finds the ids and x, y and z-coordinates of each roi in a
% volume
% INPUT: Svoboda data structure s and a volume number (>=2)
% OUTPUT: matrix M (4, # cells in volume), where first row is the cell id,
% and the second, third and fourth the x, y and z-coordinate
[ncellsvolume, ntime] = size(s.timeSeriesArrayHash.value{volume}.valueMatrix);
%% check number of cells
ncellsok = 1;
nctot = 0;
nplanes = length(s.timeSeriesArrayHash.value{volume}.imagingPlane);
for np = 1:nplanes
nctot = nctot+length(s.timeSeriesArrayHash.value{volume}.imagingPlane{np}.ids);
end
if ~(nctot == ncellsvolume)
disp('Warning: Number of recorded cells not equal to sum of cells in planes')
ncellsok = 0;
end
M = nan*ones(4,ncellsvolume);
nn=0;
for np = 1:nplanes
%% find cell ids for this plane
cellidsthisplane = s.timeSeriesArrayHash.value{volume}.imagingPlane{np}.ids;
ncellsthisplane = length(cellidsthisplane);
if ncellsok == 1
% Extra check: are cell ids from different locations (value matrix and imaging planes) the same?
if ~ (sum(s.timeSeriesArrayHash.descrHash{volume}.value{1}(np).roiIds == cellidsthisplane) == ncellsthisplane)
disp('Something went wrong, cell ids not correct')
keyboard
end
else
% Since # cells in value matrix is not the same as # cells in
% imaging planes, check for each cell id in imaging planes whether it was recorded
cellokvec = zeros(1,ncellsthisplane);
for nc = 1:ncellsthisplane
if ~isempty(find(s.timeSeriesArrayHash.value{volume}.ids == cellidsthisplane(nc)))
cellokvec(nc) = 1;
end
end
ncellsthisplane = sum(cellokvec);
if ncellsthisplane>0
cellidsthisplane = cellidsthisplane(cellokvec);
else
disp(['Plane ' num2str(np) ' has no recorded cells'])
cellidsthisplane = [];
end
end
M(1, nn+1:nn+ncellsthisplane) = cellidsthisplane;
%% z-coordinate (depth):
% planes are 15 micrometers separated
% each volume contains 3 planes
% so volume 2, plane 1 correspondse to '0', plane 2 to '15', etc
M(4, nn+1:nn+ncellsthisplane) = ((np-1)+(volume-2)*3)*15;
%% find x and y coordinate: take average over pixel locations
for nc = 1:ncellsthisplane
M(2:3, nn+nc) = mean(s.timeSeriesArrayHash.descrHash{volume}.value{1}(np).rois(nc).cornersXY,2);
end
nn = nn+ncellsthisplane;
end
end