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

Grain Segmentation: made it possible to reach iterative thresholding from multiple codepaths 


Initially introduced by Laura

Signed-off-by: default avatarNicola Vigano <nicola.vigano@esrf.fr>
parent 2160f8a9
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
4_grains/GtThreshold.m
+ 142
− 112
View file @ f5b2b080
......@@ -18,10 +18,12 @@ classdef GtThreshold < handle
obj.param = load('parameters.mat');
obj.param = obj.param.parameters;
[obj.param.rec, rej_params] = parse_pv_pairs(obj.param.rec, varargin);
if (~isfield(obj.param.rec, 'do_morph_recon'))
obj.param.rec.do_morph_recon = false;
par_thr = obj.getThrParams();
[par_thr, rej_params] = parse_pv_pairs(par_thr, varargin);
if (~isfield(par_thr, 'do_morph_recon'))
par_thr.do_morph_recon = false;
end
obj.param.rec.thresholding = par_thr;
obj.conf = parse_pv_pairs(obj.conf, rej_params);
end
......@@ -35,7 +37,7 @@ classdef GtThreshold < handle
end
defaults = [];
defaults.percentile = obj.param.rec.percentile2;
defaults.percentile = obj.param.rec.thresholding.percentile;
defaults.phases = {};
for phaseID = 1:length(sample.phases)
......@@ -55,7 +57,7 @@ classdef GtThreshold < handle
[defaults, sample] = obj.generateDefaults(filetable);
if (exist('percentile', 'var'))
obj.param.rec.percentile2 = percentile;
obj.param.rec.thresholding.percentile = percentile;
parameters = obj.param;
save('parameters.mat', 'parameters', '-v7.3');
......@@ -101,9 +103,9 @@ classdef GtThreshold < handle
grain = obj.loadGrain(phaseID, grainID);
% Actually threshold now!
if (isfield(obj.param.rec, 'ANDY') && obj.param.rec.ANDY)
if (obj.checkIsIterative())
% iterative search
grain = obj.thresholdAutoSingleGrainANDY(grain);
grain = obj.thresholdAutoSingleGrainIter(grain);
else
% normal way
grain = obj.thresholdAutoSingleGrain(grain);
......@@ -131,7 +133,7 @@ classdef GtThreshold < handle
gauge.delete();
disp('Now writing to sample.mat')
samplefile = fullfile(obj.param.acq.dir, '4_grains', 'sample.mat');
samplefile = fullfile(obj.param.acq.dir, '4_grains', 'sample.mat');
sample.mergeToFile(filetable, samplefile, 'sample');
if (nargout >= 1)
......@@ -149,12 +151,19 @@ classdef GtThreshold < handle
if (exist('center', 'var'))
do_morpho = true;
else
center = zeros(1, 3);
center = [];
do_morpho = false;
warning('THRESHOLD:less_capability', ...
'Morphological reconstruction is disabled')
end
grain = obj.segmentAndDoMorph(grain, threshold, do_morpho, center);
if (obj.checkIsIterative())
% iterative search
grain = obj.thresholdAutoSingleGrainIter(grain, threshold, center);
else
% normal way
grain = obj.segmentAndDoMorph(grain, threshold, do_morpho, center);
end
if isempty(grain.seg)
fprintf('\n\nWarning! Empty volume in grain %d.\n\n', grainID);
......@@ -176,7 +185,13 @@ classdef GtThreshold < handle
fprintf('Segmenting Grain %d, from phase %d..', grainID, phaseID);
grain = obj.loadGrain(phaseID, grainID);
grain = obj.thresholdAutoSingleGrain(grain);
if (obj.checkIsIterative())
% iterative search
grain = obj.thresholdAutoSingleGrainIter(grain);
else
% normal way
grain = obj.thresholdAutoSingleGrain(grain);
end
if isempty(grain.seg)
fprintf('\n\nWarning! Empty volume in grain %d.\n\n',grainID);
......@@ -198,10 +213,15 @@ classdef GtThreshold < handle
% GTTHRESHOLD/volumeThreshold
% segvol = volumeThreshold(obj, volume, threshold, center)
if (exist('center', 'var') && all(center > 0))
if (isempty(volume))
gtError('GtThreshold:wrong_argument', ...
'Volume to segment is empty!')
end
if (exist('center', 'var') && ~isempty(center) && all(center > 0))
do_morpho = true;
else
center = zeros(1, 3);
center = [];
do_morpho = false;
warning('THRESHOLD:less_capability', ...
'Morphological reconstruction is disabled')
......@@ -232,6 +252,11 @@ classdef GtThreshold < handle
grain.vol = grain_details.ODF6D.intensity;
grain.shift = grain_details.ODF6D.shift;
end
if (~isfield(grain, 'vol') || isempty(grain.vol))
gtError('GtThreshold:wrong_argument', ...
'Volume to segment is empty!')
end
end
function saveGrain(obj, phaseID, grainID, grain)
......@@ -258,6 +283,25 @@ classdef GtThreshold < handle
end
end
function is_iterative = checkIsIterative(obj)
is_iterative = isfield(obj.param.rec, 'thresholding') ...
&& (obj.param.rec.thresholding.iter_num > 0);
end
function thr_params = getThrParams(obj)
% GTTHRESHOLD/getThrParams
% FUNCTION that deals with parameters versions
if (isfield(obj.param.rec, 'thresholding'))
thr_params = obj.param.rec.thresholding;
else
warning('GtThreshold:old_params', ...
'Old "rec" parameters. You should produce a new version.')
thr_params = obj.param.rec;
thr_params.percentile = thr_params.percentile2;
end
end
function grain = thresholdAutoSingleGrain(obj, grain)
% GTTHRESHOLD/THRESHOLDAUTOSINGLEGRAIN Segment 3D grain volume
%
......@@ -265,7 +309,7 @@ classdef GtThreshold < handle
%
% INPUTS: grain.vol graylevel reconstruction of grain volume
% rec contents of parameters.rec structure:
% .percentile2 adaptive threshold see below
% .percentile adaptive threshold see below
% .do_morph_recon remove not connected areas by morphological
% reconstruction (<logical> {true})
% OUTPUTS:
......@@ -280,13 +324,10 @@ classdef GtThreshold < handle
% We want a threshold value below the peak corresponding to the grain, but
% above the background peak.
if (~isfield(grain, 'vol') || isempty(grain.vol))
gtError('GtThreshold:wrong_argument', ...
'Volume to segment is empty!')
end
vol = grain.vol;
par_thr = obj.getThrParams();
% calulate intensity in the central region of the grain
centerx = round(size(vol, 1) * 0.45) : round(size(vol, 1) * 0.55);
centery = round(size(vol, 2) * 0.45) : round(size(vol, 2) * 0.55);
......@@ -300,22 +341,91 @@ classdef GtThreshold < handle
stdvol = std( vol(indexes) );
meanvol = mean( vol(indexes) );
%if (stdvol > meanvol)
% threshold_val = obj.param.rec.percentile1/100 * meanvol;
%else
% threshold_val = meanvol - obj.param.rec.percentile2/100 * stdvol;
%end
% The above cases are not continuous! Use this instead:
threshold_val = meanvol - obj.param.rec.percentile2/100 * stdvol;
center(1) = round(size(vol, 1) / 2);
center(2) = round(size(vol, 2) / 2);
center(3) = round(size(vol, 3) / 2);
threshold_val = meanvol - par_thr.percentile/100 * stdvol;
center = round([size(vol, 1), size(vol, 2), size(vol, 3)] / 2);
grain = obj.segmentAndDoMorph( grain, threshold_val, ...
obj.param.rec.do_morph_recon, ...
par_thr.do_morph_recon, ...
center);
end
function grain = thresholdAutoSingleGrainIter(obj, grain, threshold_val, center)
% GTTHRESHOLD/THRESHOLDAUTOSINGLEGRAINITER Segment 3D grain volume
% iteratively
%
% grain = thresh.thresholdAutoSingleGrainIter(grain, threshold_val, center)
%
% INPUTS: grain.vol graylevel reconstruction of grain volume
%
% OUTPUTS:
% grain.thresh Threshold value
% grain.seg Segmented grain volume
% grain.Area Area of grain
% grain.segbb BoundingBox (gt convention) for insertion of seg
% volume into sample volume
%
% explanation - we want to end up with a grain that has a volume
% which is consistent with the difspot sizes. Adjust the threshold
% until we achieve this.
% We need other fields:
% parameters.rec.thresholding.iter_factor
% parameters.rec.thresholding.num_iter
vol = grain.vol;
par_thr = obj.getThrParams();
% slow and stupid for the moment
% get bounding box sizes from pair table
query = sprintf([ ...
'SELECT avbbXsize, avbbYsize ' ...
' FROM %s ' ...
' WHERE grainid = %d' ], ...
obj.param.acq.pair_tablename, grain.id);
[bbx, bby] = mym(query);
bbx = median(bbx);
bby = median(bby);
target = bbx * bbx * bby * par_thr.iter_factor;
if (~exist('threshold_val', 'var') || isempty(threshold_val))
% calulate intensity in the central region of the grain
centerx = round(size(vol, 1) * 0.45) : round(size(vol, 1) * 0.55);
centery = round(size(vol, 2) * 0.45) : round(size(vol, 2) * 0.55);
centerz = round(size(vol, 3) * 0.45) : round(size(vol, 3) * 0.55);
% use this as a first guess
threshold_val = gtImgMeanValue(vol(centerx, centery, centerz)) ;
end
if (~exist('center', 'var') || isempty(center) || any(center <= 0))
% centre point for morphological
center = round([size(vol, 1), size(vol, 2), size(vol, 3)] / 2);
end
happy = false; %:-(
count = 0;
while (~happy && count < par_thr.iter_num)
% segment, test size
grain = obj.segmentAndDoMorph( grain, threshold_val, ...
par_thr.do_morph_recon, center);
test = sum(grain.seg(:));
delta = 100 * (test - target) / target;
if (abs(delta) > 1)
delta_sign = sign(delta);
delta_mag = min(threshold_val * abs(delta) / 100, 0.2);
threshold_val = threshold_val + (delta_sign * delta_mag);
else
happy = true; %:-)
end
count = count + 1;
end
fprintf('finished grain %d after %d iterations, error %0.2f\n', ...
grain.id, count, delta)
end
function grain = segmentAndDoMorph(obj, grain, threshold, do_morpho, center)
% GTTHRESHOLD/SEGMENTANDDOMORPH
......@@ -352,15 +462,10 @@ classdef GtThreshold < handle
function segvol = segmentVol(~, vol, threshold, do_morpho, errorMsg, center)
% GTTHRESHOLD/SEGMENTVOL
if (isempty(vol))
gtError('GtThreshold:wrong_argument', ...
'Volume to segment is empty!')
end
segvol = vol > threshold;
% morphological reconstruction removes not connected areas
if (do_morpho && all(center > 0))
if (do_morpho && ~isempty(center))
try
marker = false(size(vol));
marker(center(1), center(2), center(3)) = true;
......@@ -371,80 +476,5 @@ classdef GtThreshold < handle
end
end
end
function grain = thresholdAutoSingleGrainANDY(obj, grain)
% GTTHRESHOLD/THRESHOLDAUTOSINGLEGRAINANDY Segment 3D grain volume
%
% Test version, Andy
%
% function grain = thresh.thresholdAutoSingleGrainANDY(grain)
%
% INPUTS: grain.vol graylevel reconstruction of grain volume
%
% OUTPUTS:
% grain.thresh Threshold value
% grain.seg Segmented grain volume
% grain.Area Area of grain
% grain.segbb BoundingBox (gt convention) for insertion of seg
% volume into sample volume
%
% explanation - we want to end up with a grain that has a volume
% which is consistent with the difspot sizes. Adjust the threshold
% until we achieve this. Probably need another field in
% parameters.rec.ANDY_f and a constant.
if (~isfield(grain, 'vol') || isempty(grain.vol))
gtError('GtThreshold:wrong_argument', ...
'Volume to segment is empty!')
end
vol = grain.vol;
% slow and stupid for the moment
% get bounding box sizes from pair table
query = sprintf([ ...
'SELECT avbbXsize, avbbYsize ' ...
'FROM %s ' ...
'WHERE grainid = %d' ], ...
obj.param.acq.pair_tablename, grain.id);
[bbx, bby] = mym(query);
bbx = median(bbx);
bby = median(bby);
target = bbx * bbx * bby * obj.param.rec.ANDY_f ;
% calulate intensity in the central region of the grain
centerx = round(size(vol, 1) * 0.45) : round(size(vol, 1) * 0.55);
centery = round(size(vol, 2) * 0.45) : round(size(vol, 2) * 0.55);
centerz = round(size(vol, 3) * 0.45) : round(size(vol, 3) * 0.55);
% use this as a first guess
threshold_val = gtImgMeanValue(vol(centerx, centery, centerz)) ;
% centre point for morphological
center(1) = round(size(vol, 1) / 2);
center(2) = round(size(vol, 2) / 2);
center(3) = round(size(vol, 3) / 2);
happy = false; %:-(
count = 0;
while (~happy && count<obj.param.rec.ANDY_nitr)
% segment, test size
grain = obj.segmentAndDoMorph( grain, threshold_val, ...
obj.param.rec.do_morph_recon, center);
test = sum(grain.seg(:));
delta = 100 * (test-target)/target;
if (abs(delta) > 1)
delta_sign = sign(delta);
delta_mag = min(threshold_val*abs(delta)/100, 0.2);
threshold_val = threshold_val + (delta_sign * delta_mag);
else
happy = true; %:-)
end
count = count+1 ;
end
fprintf('finished grain %d after %d iterations, error %0.2f\n', grain.id, count, delta)
end
end
end
5_reconstruction/gtRecDefaultParameters.m
+ 3
− 3
View file @ f5b2b080
......@@ -19,7 +19,7 @@ function par_rec = gtRecDefaultParameters(varargin)
par_rec.grains = gtRecGrainsDefaultParameters(conf.grains_algo);
% Grain segmentation related values
par_rec.percentile1 = 25;
par_rec.percentile2 = 20;
par_rec.do_morph_recon = true;
par_rec.thresholding = struct( ...
'percentile', 20, 'do_morph_recon', true, ...
'num_iter', 0, 'iter_factor', 1);
end
zUtil_Parameters/build_list_v2.m
+ 12
− 7
View file @ f5b2b080
......@@ -605,17 +605,12 @@ list.fsim(end+1, :) = {'oversizeVol', ...
% reconstruction
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
list.rec = cell(0, 4);
list.rec(end+1, :) = {'percentile1', ...
'Alternative adaptive threshold: percentile/100 * mean below grain average value', 'double', 1};
list.rec(end+1, :) = {'percentile2', ...
'Adaptive threshold: percentile/100 * std below grain average value', 'double', 1};
list.rec(end+1, :) = {'do_morph_recon', ...
'Perform morphological reconstruction during segmentation', 'logical', 1};
list.rec(end+1, :) = {'absorption', ...
'Absorption volume reconstruction information', 'struct', 0};
list.rec(end+1, :) = {'grains', ...
'Grains volume reconstruction information', 'struct', 0};
list.rec(end+1, :) = {'thresholding', ...
'Grains volume segmentation information', 'struct', 0};
list.rec__absorption = cell(0, 4);
list.rec__absorption(end+1, :) = {'algorithm', ...
......@@ -637,6 +632,16 @@ list.rec__grains(end+1, :) = {'list', ...
list.rec__grains(end+1, :) = {'options', ...
'Reconstruction algorithm''s options', 'struct', 0};
list.rec__thresholding = cell(0, 4);
list.rec__thresholding(end+1, :) = {'percentile', ...
'Threshold: percentile/100 * std', 'double', 1};
list.rec__thresholding(end+1, :) = {'do_morph_recon', ...
'Perform morphological reconstruction during segmentation', 'logical', 1};
list.rec__thresholding(end+1, :) = {'num_iter', ...
'Number of iterations used in thresholding grains'' volume', 'double', 1};
list.rec__thresholding(end+1, :) = {'iter_factor', ...
'Weight factor for the iterative segmentation', 'double', 1};
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% fed
......
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