Commit 26b2b8cd authored by Zheheng Liu's avatar Zheheng Liu
Browse files

some changes on gtPreprocessing and twin cluster reconstruction.


Signed-off-by: Zheheng Liu's avatarZheheng Liu <zheheng.liu@esrf.fr>
parent 22241f7d
......@@ -134,9 +134,14 @@ save(parameters_name,'parameters');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
darkname = fullfile(parameters.acq.dir,'0_rawdata',parameters.acq.name,'dark.edf');
darkendname = fullfile(parameters.acq.dir,'0_rawdata',parameters.acq.name,'darkend0000.edf');
[dark,info] = edf_read(darkendname);
if ~exist(darkname,'file')
disp('Creating dark.edf image of scan.')
[dark,info] = edf_read(darkendname);
for ii_dark = 1:(parameters.acq.ndark-1)
darkendname = fullfile(parameters.acq.dir,'0_rawdata',parameters.acq.name,sprintf('darkend04%.edf',ii_dark));
darkend = edf_read(darkendname, [], [], info);
dark = dark + darkend;
end
dark = dark/parameters.acq.ndark;
info.datatype = 'float32';
edf_write(dark,darkname,info);
......
function [int_profile, refs] = gtGrainComputeIncomingBeamIntensity(gr, parameters, det_index, refs, radius)
function [int_profile, refs] = gtGrainComputeIncomingBeamIntensity(gr, parameters, det_index, refs, radius, int_margin)
% FUNCTION [int_profile, refs] = gtGrainComputeIncomingBeamIntensity(gr, parameters, det_index, refs, radius)
if (parameters.acq.interlaced_turns > 0)
......@@ -20,6 +20,14 @@ function [int_profile, refs] = gtGrainComputeIncomingBeamIntensity(gr, parameter
radius = 3;
end
if (~exist('int_margin', 'var') || isempty(int_margin))
try
int_margin = load('mean_beam_intensities_in_margins_of_flat_field.mat');
int_margin = int_margin.mean_beam_intensities_in_margins_of_flat_field;
catch
int_margin = 1;
end
end
omegas = gr.allblobs.omega(gr_included);
gc_det_pos = gtGeoGrainCenterSam2Det(gr.center, omegas, parameters);
......@@ -40,6 +48,12 @@ function [int_profile, refs] = gtGrainComputeIncomingBeamIntensity(gr, parameter
int_profile = interp3(refs, det_positions_tot(:, 1, :), det_positions_tot(:, 2, :), det_positions_tot(:, 3, :), 'linear', 0);
int_profile = sum(bsxfun(@times, int_profile, weights), 3);
int_nums = img_nums * (numel(int_margin) - 1) / parameters.acq.nproj / 2 + 1;
int_nums(int_nums < 1) = 1;
int_nums(int_nums > numel(int_margin)) = numel(int_margin);
int_margin = [int_margin, int_margin(end)];
int_time = (1 - int_nums + floor(int_nums)) .* int_margin(floor(int_nums)) + (int_nums - floor(int_nums)) .* int_margin(floor(int_nums) + 1);
int_profile = int_profile .* int_time(:);
end
function refs = load_reference_imgs(parameters)
......
......@@ -79,13 +79,17 @@ function grn = gtGrainMerge(gr1, gr2, p, mode, flag_merge_bl)
gr2_bl_to_merge = gr2.proj(ii_d).bl(merge_blob_ind_on_gr2_bl_list);
if ~isempty(gr1_bl_to_merge) && ~isempty(gr2_bl_to_merge)
tmp_bl_merge = gtGrainMergeProjBl([gr1_bl_to_merge, gr2_bl_to_merge]);
elseif ~isempty(gr1_bl_to_merge)
tmp_bl_merge = gr1_bl_to_merge;
elseif ~isempty(gr2_bl_to_merge)
tmp_bl_merge = gr2_bl_to_merge;
elseif flag_merge_bl == 1
if ~isempty(gr1_bl_to_merge)
tmp_bl_merge = gr1_bl_to_merge;
elseif ~isempty(gr2_bl_to_merge)
tmp_bl_merge = gr2_bl_to_merge;
else
error('gtGrainMerge:missing_blob', ...
'Both gr1.proj(%d).bl(%d) and gr2.proj(%d).bl(%d) are empty', ii_d, merge_blob_ind_on_gr1_bl_list, ii_d, merge_blob_ind_on_gr2_bl_list)
end
else
error('gtGrainMerge:missing_blob', ...
'Both gr1.proj(%d).bl(%d) and gr2.proj(%d).bl(%d) are empty', ii_d, merge_blob_ind_on_gr1_bl_list, ii_d, merge_blob_ind_on_gr2_bl_list)
tmp_bl_merge = grn.proj(ii_d).bl(ii_merge_bl);
end
grn.proj(ii_d).bl(ii_merge_bl) = tmp_bl_merge;
end
......
......@@ -39,38 +39,19 @@ function bl_merge = gtGrainMergeProjBl(bl_array, difspotID_array)
if exist('difspotID_array','var') && numel(difspotID_array(:)) == numel(bl_array(:))
[~, tmp_ind] = unique(difspotID_array(:));
bl_array = bl_array(tmp_ind);
% Merge intensity
intensity_array = cat(1, bl_array(:).intensity);
tmp_intensity = sum(intensity_array);
intensity_array = intensity_array/tmp_intensity;
% Merge mask and intm
tmp_mask = false(bl_merge.bbsize);
tmp_intm = single(zeros(bl_merge.bbsize));
for i_bl = 1:numel(bl_array(:))
tmp_vol_ind = sfVolInd(bound(:, :, i_bl), bound_merge, bl_merge.bbsize);
tmp_mask(tmp_vol_ind) = tmp_mask(tmp_vol_ind) + bl_array(i_bl).mask(:);
tmp_intm(tmp_vol_ind) = tmp_intm(tmp_vol_ind) + intensity_array(i_bl) * bl_array(i_bl).intm(:);
end
else
% Merge mask, intm and intensity.
intensity_array = cat(1, bl_array(:).intensity);
tmp_mask = false(bl_merge.bbsize);
tmp_intm = single(zeros(bl_merge.bbsize));
for i_bl = 1:numel(bl_array(:))
tmp_vol_ind = sfVolInd(bound(:, :, i_bl), bound_merge, bl_merge.bbsize);
% Ignore the bl whose mask overlaps the ones of previous bls
if any(bl_array(i_bl).mask(:) & tmp_mask(tmp_vol_ind))
intensity_array(i_bl) = 0;
else
tmp_mask(tmp_vol_ind) = tmp_mask(tmp_vol_ind) + bl_array(i_bl).mask(:);
tmp_intm(tmp_vol_ind) = tmp_intm(tmp_vol_ind) + intensity_array(i_bl) * bl_array(i_bl).intm(:);
end
end
tmp_intensity = sum(intensity_array(:));
tmp_intm = tmp_intm/tmp_intensity;
end
% Merge mask, intm and intensity.
tmp_mask = false(bl_merge.bbsize);
tmp_intm = single(zeros(bl_merge.bbsize));
for i_bl = 1:numel(bl_array(:))
tmp_vol_ind = sfVolInd(bound(:, :, i_bl), bound_merge, bl_merge.bbsize);
% Remove the overlaps
tmp_not_overlapped_vol = (bl_array(i_bl).mask(:) & ~tmp_mask(tmp_vol_ind)) .* bl_array(i_bl).intm(:) .* bl_array(i_bl).intensity;
tmp_mask(tmp_vol_ind) = tmp_mask(tmp_vol_ind) | bl_array(i_bl).mask(:);
tmp_intm(tmp_vol_ind) = tmp_intm(tmp_vol_ind) + tmp_not_overlapped_vol;
end
tmp_intensity = sum(tmp_intm(:));
tmp_intm = tmp_intm/tmp_intensity;
bl_merge.intensity = tmp_intensity;
bl_merge.mask = tmp_mask;
......
......@@ -15,5 +15,5 @@ function [int_fact,phi1, phi2, range] = gtCrystGeometricFactor(thetatype, eta, r
[~, ~, phi1] = gtMathsCordLength(R1, r, k(1));
[~, ~, phi2] = gtMathsCordLength(R2, r, k(3));
range = abs(phi1 - phi2);
int_fact = range./min(range);
int_fact = range/min(range);
end
......@@ -447,7 +447,7 @@ classdef GtPhase < handle
'bb_r_space', bb_r_space, 'bb_o_space', bb_o_space);
end
function cluster = makeCluster(included_ids, ref_r_vec, bb_r_space, bb_o_space, cluster_type, active)
function cluster = makeCluster(included_ids, ref_r_vec, bb_r_space, bb_o_space, cluster_type, twin_paths, active)
% FUNCTION cluster = GtPhase.makeCluster(included_ids, ref_r_vec, bb_r_space, bb_o_space, cluster_type, active)
% Cluster types:
% 0 - Regular cluster with grains that are close both in real
......@@ -465,7 +465,7 @@ classdef GtPhase < handle
cluster = struct( ...
'included_ids', included_ids, 'ref_r_vec', ref_r_vec, ...
'bb_r_space', bb_r_space, 'bb_o_space', bb_o_space, ...
'cluster_type', cluster_type, 'active', logical(active));
'cluster_type', cluster_type, 'twin_paths', twin_paths, 'active', logical(active));
end
function obj = loadobj(obj)
......
......@@ -22,6 +22,7 @@ classdef Gt6DReconstructionAlgorithmFactory < handle
% 'lorentz' : num_interp is rescaled by the Lorentz factor
num_interp_type = 'lorentz';
mean_beam_intensities_in_margins_of_flat_field = 1;
parameters;
end
......@@ -101,6 +102,12 @@ classdef Gt6DReconstructionAlgorithmFactory < handle
% use_predicted_scatter_ints to remove problems due to overlaps
self.use_matrix_row_rescaling = true;
self.use_predicted_scatter_ints = max(self.use_predicted_scatter_ints, 1);
try
tmp = load('mean_beam_intensities_in_margins_of_flat_field.mat');
self.mean_beam_intensities_in_margins_of_flat_field = tmp.mean_beam_intensities_in_margins_of_flat_field;
catch
self.mean_beam_intensities_in_margins_of_flat_field = 1;
end
blobs = reshape(blobs, [], 1);
......@@ -1339,7 +1346,7 @@ classdef Gt6DReconstructionAlgorithmFactory < handle
if (self.use_predicted_scatter_ints > 1)
try
fprintf(' - Applying incoming beam intensity corrections\n')
[beam_intensities, refs] = gtGrainComputeIncomingBeamIntensity(ref_or, self.parameters, det_ind);
[beam_intensities, refs] = gtGrainComputeIncomingBeamIntensity(ref_or, self.parameters, det_ind, [], [], self.mean_beam_intensities_in_margins_of_flat_field);
beam_in_ints_avg = beam_intensities(sel);
beam_in_ints_avg = beam_in_ints_avg / mean(beam_in_ints_avg);
catch mexc
......@@ -1378,7 +1385,7 @@ classdef Gt6DReconstructionAlgorithmFactory < handle
or.allblobs = ors_allbls(ii_o);
if (numel(beam_in_ints_avg) > 1)
try
[beam_in_int_or, refs] = gtGrainComputeIncomingBeamIntensity(or, self.parameters, det_ind, refs);
[beam_in_int_or, refs] = gtGrainComputeIncomingBeamIntensity(or, self.parameters, det_ind, refs, [], self.mean_beam_intensities_in_margins_of_flat_field);
beam_in_int_or = beam_in_int_or(sel);
beam_in_int_or = beam_in_int_or / mean(beam_in_int_or);
catch mexc
......
Markdown is supported
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