Skip to content
Snippets Groups Projects
Commit dd05096d authored by Nicola Vigano's avatar Nicola Vigano
Browse files

Segmentation/doublethr: tried to cope with strange behavior and reformatted code a little bit 


Signed-off-by: default avatarNicola Vigano <nicola.vigano@esrf.fr>

git-svn-id: https://svn.code.sf.net/p/dct/code/trunk@652 4c865b51-4357-4376-afb4-474e03ccb993
parent 221fddf7
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
2_difspot/gtSegmentDiffractionBlobs_doublethr.m
+ 146
− 209
View file @ dd05096d
......@@ -68,21 +68,20 @@ function gtSegmentDiffractionBlobs_doublethr(~, ~, workingdirectory)
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% setup
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if isdeployed
global GT_DB
global GT_MATLAB_HOME
load('workspaceGlobal.mat');
% first=str2double(first);
% last=str2double(last);
verbose = false; % don't show all the messages - as this makes the .out file very large
% don't show all the messages - as this makes the .out file very large
out = GtConditionalOutput(false);
disp('running compiled, so will show minimal messages only')
else
verbose = true;
out = GtConditionalOutput(true);
disp('displaying the info messages')
end % special case for running interactively in current directory
......@@ -118,19 +117,18 @@ thr_grow_rat = seg.thr_grow_rat; % relative threshold for growing seeds
% thr_grow is used for growing the seeds, and linearly depends on thr_seed.
% This linearity is constrained by the upper and lower limits for thr_grow:
if isfield(parameters.seg,'thr_grow_low')
if isfield(parameters.seg, 'thr_grow_low')
thr_grow_low = parameters.seg.thr_grow_low;
else
thr_grow_low = [];
end
if isfield(parameters.seg,'thr_grow_high')
if isfield(parameters.seg, 'thr_grow_high')
thr_grow_high = parameters.seg.thr_grow_high;
else
thr_grow_high = [];
end
if isfield(seg, 'extendblobinc')
extendblobinc = seg.extendblobinc; % when increasing bounding box use these incs.
else
......@@ -143,22 +141,18 @@ if isfield(seg, 'background_subtract') && seg.background_subtract==true
background_subtract = true;
end
xdet = parameters.acq.xdet;
ydet = parameters.acq.ydet;
% if the median background is not subtracted from the fulls as they are
% created, then we need a record of what the median background of each full
% is, to be used by sfGetVol (which only reads part images and therefore
% can't calculate the median itself.
fullmedianvals = NaN(totproj+1,1); % if totproj-7199, need 7200
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% connect to database
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
gtDBConnect();
table_difblob = sprintf('%sdifblob',parameters.acq.name);
......@@ -173,26 +167,23 @@ table_seeds = sprintf('%sseeds',parameters.acq.name);
% for controlling parallel jobs
table_seeds_tmp = sprintf('%sseeds_tmp',parameters.acq.name);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% wait until gtSeedSegmentation has finished
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wait = 1;
while wait
test = mym(sprintf('select count(val) from %sfullmedianvals', ...
test = mym(sprintf('SELECT count(val) FROM %sfullmedianvals', ...
parameters.acq.name));
if test > totproj % test will reach (for eg) 7200, totproj will be 7199
wait = 0; % go go go!
else
pause(30)
disp(sprintf('waiting for SeedSegmentation to finish... done %d of %d',...
test, totproj+1))
fprintf('Waiting for SeedSegmentation to finish: done %d of %d\n', ...
test, totproj+1)
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% get seed data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
......@@ -207,43 +198,37 @@ seedlist = [a, b, c, d, e]; % don't need bad at this point
% main loop over seeds
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
edfHeader = [];
edf_info = [];
for ss = 1:length(seedlist)
% this seed data
seedID = seedlist(ss, 1);
seedx = seedlist(ss, 2);
seedy = seedlist(ss, 3);
seedz = seedlist(ss, 4);
seedInt = seedlist(ss, 5);
seedX = seedlist(ss, 2);
seedY = seedlist(ss, 3);
seedZ = seedlist(ss, 4);
seedInt = seedlist(ss, 5); % Intensity
% tests - has this seed been taken? file table style
try
test = mym(sprintf('insert into %s (seedID) values (%d)',...
test = mym(sprintf('INSERT INTO %s (seedID) VALUES (%d)',...
table_seeds_tmp, seedID));
catch mexc
% Let's just silence the error!
test = false;
end
if test ~= 1
if verbose
disp('another job has already taken this seed - skipping')
end
if (test ~= 1)
out.fprintf('another job has already taken this seed - skipping\n')
continue
end
% tests - is this seed now bad?
test = mym(sprintf('select bad from %s where seedID=%d', table_seeds,...
test = mym(sprintf('SELECT bad FROM %s WHERE seedID=%d', table_seeds,...
seedID));
if test == 1
if verbose
disp('this seed already bad - skipping')
end
if (test == 1)
out.fprintf('this seed already bad - skipping\n');
continue
end
......@@ -260,13 +245,13 @@ for ss = 1:length(seedlist)
% continue
% end
if verbose, disp('seed ok'); end
out.fprintf('seed ok\n');
% if seed ok, grow seed to blob
test = [1 1 1 1 1 1];
% Apply lower and upper limit for thr_grow to grow seed
thr_grow = thr_grow_rat*seedInt;
thr_grow = thr_grow_rat * seedInt;
if ~isempty(thr_grow_low) && (thr_grow < thr_grow_low)
thr_grow = thr_grow_low;
......@@ -276,88 +261,72 @@ for ss = 1:length(seedlist)
thr_grow = thr_grow_high;
end
% for the moment use a dummy bb +-15 around seed pixel
% use extendblobinc for this
ebi = ceil(extendblobinc/2);
% find min. and max. x,y,z values for the origin of the temporary seed
% bounding box
tx1 = max(seedx-ebi(1), 1);
ty1 = max(seedy-ebi(2), 1);
tx1 = max(seedX-ebi(1), 1);
ty1 = max(seedY-ebi(2), 1);
tx2 = min(seedx+ebi(1), parameters.acq.xdet);
ty2 = min(seedy+ebi(2), parameters.acq.ydet);
tx2 = min(seedX+ebi(1), parameters.acq.xdet);
ty2 = min(seedY+ebi(2), parameters.acq.ydet);
tz1 = max(seedz-ebi(3), 0);
tz2 = min(seedz+ebi(3), totproj);
tz1 = max(seedZ-ebi(3), 0);
tz2 = min(seedZ+ebi(3), totproj);
% temporary bounding box around seed
bb = [tx1 ty1 tz1 tx2-tx1+1 ty2-ty1+1 tz2-tz1+1];
while any(test)
% extend the bb
bb = sfExtendBB(test, bb, extendblobinc, totproj, xdet, ydet);
% get the volume - trying to save header decoding
if (isfield(edfHeader, 'headerlength'))
[greyvol, fullmedianvals] = sfGetVol(bb, table_fullmedianvals, ...
fullmedianvals, background_subtract, seg.bbox, edfHeader);
else
[greyvol, fullmedianvals, edfHeader] = sfGetVol(bb, ...
table_fullmedianvals, fullmedianvals, background_subtract, ...
seg.bbox);
end
[greyvol, fullmedianvals, edf_info] = sfGetVol(bb, table_fullmedianvals, ...
fullmedianvals, background_subtract, seg.bbox, edf_info);
% label volume
threshvol = greyvol > thr_grow;
threshvol = bwlabeln(threshvol);
% pick out label containing seed
bloblabel = threshvol(seedy-bb(2)+1 ,seedx-bb(1)+1 ,seedz-bb(3)+1);
bloblabel = threshvol(seedY-bb(2)+1, seedX-bb(1)+1, seedZ-bb(3)+1);
threshvol = (threshvol == bloblabel);
% apply the label as a mask
rawgreyvol = greyvol; % save the un-masked greyvol
greyvol = greyvol.*threshvol;
greyvol = greyvol .* threshvol;
% test - connected to a higher seed value?
% ... set all lower connected seeds to bad
test_int = sfTestInt(greyvol, seedInt);
if test_int
if verbose
disp('blob contains a value higher than the seed - bad')
end
if (test_int)
out.fprintf('blob contains a value higher than the seed - bad\n');
end
% test - max size exceeded?
% ... set all connected lower seeds to bad - any lower seed will give a
% lower th2 and hence larger blob
test_size = sfTestBlobSize(threshvol, maxblobsize);
if test_size
if verbose
disp('blob blob exceeds the maximum size - bad')
end
if (test_size)
out.fprintf('blob blob exceeds the maximum size - bad\n')
end
if test_int || test_size
if (test_int || test_size)
% blob is bad - set contain, lower value seeds to zero
sfBadBlob(greyvol, bb, seedInt, thr_seed, table_seeds)
sfBadBlob(greyvol, bb, seedInt, thr_seed, table_seeds);
break
end
% test - does blob reach edges of sub volume?
% ...do the while loop again
test = sfTestVol(threshvol, bb, totproj, xdet, ydet);
if ~any(test)
if verbose
disp('blob fits within volume...')
end
out.fprintf('blob fits within volume...\n');
if isfield(parameters.seg, 'no_edges') && parameters.seg.no_edges
% adjust the mask to avoid edges in the projection when
......@@ -366,26 +335,18 @@ for ss = 1:length(seedlist)
end
sfWriteBlob(greyvol, bb, minblobsize, table_difblob, ...
table_difblobdetails, seedID)
table_difblobdetails, seedID);
else
if verbose
disp('increasing subvolume around blob')
end
out.fprintf('increasing subvolume around blob\n');
end
end % - while growing blob
end % - of seedlist loop
end % - of main function
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% SUB-FUNCTIONS
% SUB-FUNCTIONS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% write a good blob to the database
......@@ -432,7 +393,7 @@ if all(bb(3:4) > minblobsize(1:2)) % size test
'zIndex', zz, ...
'graylevel', vals);
end
% could write difspot.edfs here, without ever passing
% through difblobs in the database. Would need to add
% functionality from gtDBBlob2SpotTable_WriteDifspots
......@@ -445,29 +406,31 @@ function test = sfTestInt(greyvol, seedInt)
% trivial test of blob max intensity
% working with double precision numbers - build in a tiny safety margin for
% rounding
test = 0;
if max(greyvol(:)) > (seedInt+0.00000000001)
test = 1;
end
test = (max(greyvol(:)) > (seedInt+0.00000000001));
end % of sub-function
function test = sfTestBlobSize(vol, maxblobsize)
% find the extent of the thresholded blob, check against max size
profilex = sum(sum(vol,3), 1)';
profiley = sum(sum(vol,3), 2);
profilez = squeeze(sum(sum(vol,1), 2));
sizex = find(profilex, 1, 'last') - find(profilex, 1, 'first');
sizey = find(profiley, 1, 'last') - find(profiley, 1, 'first');
sizez = find(profilez, 1, 'last') - find(profilez, 1, 'first');
if any([sizex sizey sizez] > maxblobsize(1:3))
test = 1;
else
test = 0;
end
profilex = sum(sum(abs(vol), 3), 1)';
profiley = sum(sum(abs(vol), 3), 2);
profilez = squeeze(sum(sum(abs(vol), 1), 2));
sizex = find(profilex, 1, 'last') - find(profilex, 1, 'first');
sizey = find(profiley, 1, 'last') - find(profiley, 1, 'first');
sizez = find(profilez, 1, 'last') - find(profilez, 1, 'first');
test = any([sizex sizey sizez] > maxblobsize(1:3));
% Additional test for consistency check! should not be needed.
[sizex, sizey, sizez] = size(vol);
if ((~test) && any([sizex sizey sizez] > (2 * maxblobsize(1:3))) )
test = true;
% If this happens, it's not a good thing: the growth of the blob is gone
% beyond the limits and it is now including useless things
warning('SEGMENTATION:logic_fail', ...
'Blob BB extended with zeros beyond the double of blob size');
end
end % of sub-function
......@@ -476,32 +439,31 @@ function sfBadBlob(greyvol, bb, seedInt, thr_seed, table_seeds)
% set these to bad=1 in the table
% voxels at the edge of the subvol could be indentified wrongly as seeds.
% mostly copied from seed_threshold
for z = 1:size(greyvol, 3)
slice = greyvol(:, :, z);
if (max(slice(:)) > 0)
labels = bwlabel(slice>thr_seed);
for ii = 1:max(labels)
spot = find(labels == ii);
[val, index] = max(slice(spot));
[yy, xx] = ind2sub(size(slice), spot(index));
if (xx == 1) || (yy == 1) || (xx == size(slice, 2)) ...
|| (yy == size(slice, 1)) || (val > seedInt)
% not a real seed / intensity higher than original seed
continue
else
% set to bad in the seeds table
mym(sprintf('update %s set bad=1 where x=%d and y=%d and z=%d', ...
table_seeds, xx+bb(1)-1, yy+bb(2)-1, z+bb(3)-1))
end
end % loop through labels
end % skip empty slices
end % z loop
for zz = 1:size(greyvol, 3)
slice = greyvol(:, :, zz);
if (max(slice(:)) > 0)
labels = bwlabel(slice > thr_seed);
for ii = 1:max(labels)
spot = find(labels == ii);
[val, index] = max(slice(spot));
[yy, xx] = ind2sub(size(slice), spot(index));
if (xx == 1) || (yy == 1) || (xx == size(slice, 2)) ...
|| (yy == size(slice, 1)) || (val > seedInt)
% not a real seed / intensity higher than original seed
continue
else
% set to bad in the seeds table
mym(sprintf('UPDATE %s SET bad=1 WHERE x=%d AND y=%d AND z=%d', ...
table_seeds, xx+bb(1)-1, yy+bb(2)-1, zz+bb(3)-1))
end
end % loop through labels
end % skip empty slices
end % z loop
end % % of sub-function
......@@ -509,91 +471,67 @@ function test = sfTestVol(vol, bb, totproj, xdet, ydet)
% does a thresholded blob extend to the edges of a volume? And which
% edges? Account for the limits of the dataset...
test = zeros(1,6); % ~same format as bb
if bb(1) > 1 % low x
if ~isempty(find(vol(:,1,:), 1, 'first'))
test(1) = 1;
end
end
if bb(2) > 1 % low y
if ~isempty(find(vol(1,:,:), 1, 'first'))
test(2) = 1;
end
end
if bb(3) > 0 % low z (omega)
if ~isempty(find(vol(:,:,1), 1, 'first'))
test(3) = 1;
end
end
if bb(1)+bb(4)-1 < xdet % high x
if ~isempty(find(vol(:,end,:), 1, 'first'))
test(4) = 1;
end
end
if bb(2)+bb(5)-1 < ydet % high y
if ~isempty(find(vol(end,:,:), 1, 'first'))
test(5) = 1;
end
end
if bb(3)+bb(6)-1 < (totproj) % high z (omega)
if ~isempty(find(vol(:,:,end), 1, 'first'))
test(6) = 1;
end
end
test = zeros(1, 6); % ~same format as bb
% low x
test(1) = ( (bb(1) > 1) && ~isempty(find(vol(:, 1, :), 1, 'first')) );
% low y
test(2) = ( (bb(2) > 1) && ~isempty(find(vol(1, :, :), 1, 'first')) );
% low z (omega)
test(3) = ( (bb(3) > 0) && ~isempty(find(vol(:, :, 1), 1, 'first')) );
% high x
test(4) = ( ((bb(1)+bb(4)-1) < xdet) && ~isempty(find(vol(:, end, :), 1, 'first')) );
% high y
test(5) = ( ((bb(2)+bb(5)-1) < ydet) && ~isempty(find(vol(end, :, :), 1, 'first')) );
% high z (omega)
test(6) = ( ((bb(3)+bb(6)-1) < totproj) && ~isempty(find(vol(:, :, end), 1, 'first')) );
end % of sub-function
function [vol, fullmedianvals, edfHeader] = sfGetVol(bb, table_fullmedianvals, ...
fullmedianvals, background_subtract, segbb, edfHeader)
function [vol, fullmedianvals, info] = sfGetVol(bb, table_fullmedianvals, ...
fullmedianvals, background_subtract, segbb, info)
%get the full volume of the given bounding box
% median background info in table
% define a volume
vol = zeros([bb(5) bb(4) bb(6)]);
% read the images
for z = bb(3):(bb(3)+bb(6)-1)
filename = sprintf('1_preprocessing/full/full%04d.edf', z);
if (~exist('edfHeader', 'var'))
edfHeader = edf_info(filename);
end
% if we don't have the median value we need, calculate it
if background_subtract
if isnan(fullmedianvals(z+1))
bg = mym(sprintf('SELECT val FROM %s WHERE ndx=%d', table_fullmedianvals, z));
fullmedianvals(z+1) = bg;
% define a volume
vol = zeros([bb(5) bb(4) bb(6)]);
% read the images
for z = bb(3):(bb(3)+bb(6)-1)
% if we don't have the median value we need, calculate it
if background_subtract
if isnan(fullmedianvals(z + 1))
bg = mym(sprintf('SELECT val FROM %s WHERE ndx = %d', ...
table_fullmedianvals, z));
fullmedianvals(z + 1) = bg;
else
bg = fullmedianvals(z + 1);
end
else
bg = fullmedianvals(z+1);
bg = 0;
end
else
bg = 0;
% read into the volume, subtracting the median background
filename = sprintf('1_preprocessing/full/full%04d.edf', z);
[img, info] = edf_read(filename, bb([1 2 4 5]), false, info);
vol(:, :, z-bb(3)+1) = img - bg;
end
% read into the volume, subtracting the median background
vol(:,:,z-bb(3)+1) = edf_read(filename, bb([1 2 4 5]), false, edfHeader) - bg;
end
% set anything inside the seg bounding box to zeros
% put seg.bbox relative to the volume
% work out which pixels to set to zero
mins = segbb(1:2) - bb(1:2) + 1;
maxs = mins + segbb(3:4) - 1;
% Cut negative (unphysical) values (should also solve some unlikely issues)
vol(vol < 0) = 0;
% fix limits
mins = max(mins, 1);
maxs = min(maxs, bb(4:5));
% set anything inside the seg bounding box to zeros
% put seg.bbox relative to the volume
% work out which pixels to set to zero
mins = segbb(1:2) - bb(1:2) + 1;
maxs = mins + segbb(3:4) - 1;
% if we have to do something about it
vol(mins(2):maxs(2), mins(1):maxs(1), :) = 0;
% fix limits
mins = max(mins, 1);
maxs = min(maxs, bb(4:5));
% if we have to do something about it
vol(mins(2):maxs(2), mins(1):maxs(1), :) = 0;
end % of sub-function
......@@ -644,19 +582,18 @@ function greyvol = sfNoEdges(rawgreyvol, threshvol)
% the outline of the summed spot is true to th2 (intensities may be
% different) but there will be no lines due to summing of subsequent images
% omega range
profilez = squeeze(sum(sum(threshvol,1), 2));
startz = find(profilez, 1, 'first');
endz = find(profilez, 1, 'last');
% xy mask
xymask = (sum(threshvol, 3)) > 0;
profilez = squeeze(sum(sum(threshvol, 1), 2));
startz = find(profilez, 1, 'first');
endz = find(profilez, 1, 'last');
% new mask
for z = startz:endz
threshvol(:,:,z) = xymask;
end
% xy mask
xymask = (sum(threshvol, 3)) > 0;
% apply the new mask
greyvol = rawgreyvol .* threshvol;
% new mask
for z = startz:endz
threshvol(:, :, z) = xymask;
end
% apply the new mask
greyvol = rawgreyvol .* threshvol;
end % of sub-function
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment