diff --git a/4_grains/GtThreshold.m b/4_grains/GtThreshold.m
index a0d3f4c2ff1fdff97d9bc611d9ccf998f1d4bd45..fd719792d07d7dc586bf7e7fe1fbbe7ac1e31534 100644
--- a/4_grains/GtThreshold.m
+++ b/4_grains/GtThreshold.m
@@ -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
diff --git a/5_reconstruction/gtRecDefaultParameters.m b/5_reconstruction/gtRecDefaultParameters.m
index f3f73e775fe9901190ef38d7a99cb275803ba6c8..b44b07d5576848a4da5802f09e4bdc27772fff8e 100644
--- a/5_reconstruction/gtRecDefaultParameters.m
+++ b/5_reconstruction/gtRecDefaultParameters.m
@@ -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
diff --git a/zUtil_Parameters/build_list_v2.m b/zUtil_Parameters/build_list_v2.m
index 9367d0da11b9e61be393d8f0d5954c65133a5207..5c3a375043c5bbd94f31bce2f5fbacd02d21067a 100644
--- a/zUtil_Parameters/build_list_v2.m
+++ b/zUtil_Parameters/build_list_v2.m
@@ -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