Newer
Older
function [samp_ors, estim_space_bbox_pix, estim_orient_bbox] = gt6DCreateProjDataFromTwinnedGrainFwdProj(parent_id, variants, phase_id, varargin)
Nicola Vigano
committed
% FUNCTION samp_ors = gt6DCreateProjDataFromTwinnedGrainFwdProj(parent_id, variants, phase_id, varargin)
% INPUT:
% parent_id: Parent grain ID (it could be the parent grain structure)
% variants: an array of structures defined like the following:
% struct('type', {'gr_id' | 'gr_str' | 'var_num'}, 'data', {...} )
% OUTPUT:
% samp_ors: is an array of grain structures containing the parent and its
% twins
if (~exist('phase_id', 'var'))
phase_id = 1;
end
conf = struct( ...
'verbose', false, ...
'det_ind', 1, ...
'min_eta', 15, ...
'rspace_oversize', 1.1, ...
'check_spots', false, ...
'stack_oversize', 1.4, ...
'save', false );
conf = parse_pv_pairs(conf, varargin);
p = gtLoadParameters();
symm = gtCrystGetSymmetryOperators(p.cryst(phase_id).crystal_system);
Nicola Vigano
committed
if (numel(parent_id) > 1)
error('gt6DCreateProjDataFromTwinnedGrainFwdProj:wrong_argument', ...
'Only one grain should be passed as parent')
Nicola Vigano
committed
end
if (~isstruct(parent_id))
grs(1) = gtLoadGrain(phase_id, parent_id);
grs(1) = gtCalculateGrain(grs(1), p);
Nicola Vigano
committed
grs = parent_id;
Nicola Vigano
committed
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
if (isfield(grs(1), 'g_twins') && ~isempty(grs(1).twin_vars))
theo_fwdproj_variants = cat(1, grs(1).twin_vars(:).R_vector);
else
theo_fwdproj_variants = [];
end
variants_num = numel(variants);
fprintf('Loading variants: ')
for ii_var = variants_num:-1:1
num_chars = fprintf('%02d (type: %s)', variants_num-ii_var+1, variants(ii_var).type);
switch (variants(ii_var).type)
case {'gr_id', 'gr_str'}
try
if (strcmpi(variants(ii_var).type, 'gr_id'))
gr_var_id = variants(ii_var).data;
grs(ii_var+1) = gtLoadGrain(phase_id, variants(ii_var).data);
else
gr_var_id = variants(ii_var).data.id;
grs(ii_var+1) = variants(ii_var).data;
end
grs(ii_var+1) = gtCalculateGrain(grs(ii_var+1), p);
catch mexc
gtPrintException(mexc, 'Reforward-projecting grain because of:')
grs(ii_var+1) = gtForwardSimulate_v2(gr_var_id, gr_var_id, ...
pwd, phase_id, 'save_grain', false, 'display_figure', false);
end
if (~isempty(theo_fwdproj_variants))
dists = theo_fwdproj_variants - repmat(grs(ii_var+1).R_vector, [size(theo_fwdproj_variants, 1), 1]);
dists = 2 * atand(sqrt(sum(dists .^ 2, 2)));
indx = find(dists < 1);
if (~isempty(indx))
fprintf('\n - Corresponding to precomputed variant: %d\n', indx);
fprintf(['Loading variants: ' repmat(' ', [1 num_chars])]);
end
end
case 'var_num'
grs(ii_var+1) = twin_fwd_sim_to_gr(grs(1), grs(1).twin_vars(variants(ii_var).data));
otherwise
error('gt6DCreateProjDataFromTwinnedGrainTheoFwdProj:wrong_argument', ...
'Variant type: %s cannot be handled', variants(ii_var).type)
end
fprintf(repmat('\b', [1 num_chars]));
fprintf(repmat(' ', [1 num_chars]));
fprintf(repmat('\b', [1 num_chars]));
end
fprintf('Done.\n')
grs_list = cat(2, grs(:).id);
num_grains = numel(grs_list);
if (conf.verbose)
fprintf('R_vectors:\n')
cat(1, grs(:).R_vector)
fprintf('Reciprocal disorientations:\n')
for ii_g_1 = 1:num_grains
grs_to_check = (ii_g_1+1):num_grains;
if (numel(grs_to_check) > 0)
fprintf(' - Disorientations from grain: %d\n', grs(ii_g_1).id)
for ii_g_2 = grs_to_check
dis_angle = gtDisorientation(grs(ii_g_1).R_vector', grs(ii_g_2).R_vector', symm);
fprintf(' + Grain %d: %f\n', grs(ii_g_2).id, dis_angle);
end
end
end
end
refgr = grs(1);
grs(2:end) = reorder_reflections(grs(2:end), refgr);
refgr_proj = refgr.proj(conf.det_ind);
ref_ondet = refgr_proj.ondet;
ref_included = refgr_proj.included;
ref_selected = refgr_proj.selected;
Nicola Vigano
committed
if (conf.verbose && false)
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
try
refgr = gtComputeIncomingBeamIntensity(refgr, p);
beam_ints = refgr.beam_intensity(ref_selected);
beam_ints = beam_ints / mean(beam_ints);
catch mexc
gtPrintException(mexc, 'Skipping beam intensity renormalization')
beam_ints = 1;
end
blob_ints = refgr.intensity(ref_selected);
vis_spots = ref_ondet(ref_included(ref_selected));
thetatypes = refgr.allblobs.thetatype(vis_spots);
if (isfield(p.cryst(phase_id), 'int_exp'))
predicted_ints = p.cryst(phase_id).int_exp(thetatypes)';
f = figure();
ax = axes('parent', f);
plot(ax, blob_ints)
hold(ax, 'on')
plot(ax, predicted_ints / mean(predicted_ints) * mean(blob_ints), 'y')
predicted_ints = predicted_ints .* refgr.allblobs.lorentzfac(vis_spots);
plot(ax, predicted_ints / mean(predicted_ints) * mean(blob_ints), 'r')
predicted_ints = predicted_ints .* beam_ints;
plot(ax, predicted_ints / mean(predicted_ints) * mean(blob_ints), 'g')
end
predicted_ints = p.cryst(phase_id).int(thetatypes)';
f = figure();
ax = axes('parent', f);
plot(ax, blob_ints)
hold(ax, 'on')
plot(ax, predicted_ints / mean(predicted_ints) * mean(blob_ints), 'y')
predicted_ints = predicted_ints .* refgr.allblobs.lorentzfac(vis_spots);
plot(ax, predicted_ints / mean(predicted_ints) * mean(blob_ints), 'r')
predicted_ints = predicted_ints .* beam_ints;
plot(ax, predicted_ints / mean(predicted_ints) * mean(blob_ints), 'g')
fprintf('%d) %d, %g\n', ...
[(1:numel(blob_ints))', refgr.allblobs.thetatype(vis_spots), ...
blob_ints - predicted_ints / mean(predicted_ints) * mean(blob_ints)]')
hold(ax, 'off')
end
if (conf.verbose)
vis_sel = ref_ondet(ref_included);
% vis_sel = ref_ondet;
fprintf('%02d) w1 %6.2f, w2 %6.2f, n1 %6.2f, n2 %6.2f, w_inds (%d<->%d), hkls [%2d %2d %2d]<->[%2d %2d %2d] <- diff: %7.2f (!! %d, L %g) [bad n1 %d n2 %d]\n', ...
[(1:numel(vis_sel))', ...
grs(1).allblobs.omega(vis_sel), ...
grs(2).allblobs.omega(vis_sel), ...
grs(1).allblobs.eta(vis_sel), ...
grs(2).allblobs.eta(vis_sel), ...
grs(1).allblobs.omind(vis_sel), ...
grs(2).allblobs.omind(vis_sel), ...
grs(1).allblobs.hklsp(vis_sel, [1 2 4]), ...
grs(2).allblobs.hklsp(vis_sel, [1 2 4]), ...
abs(mod(grs(1).allblobs.omega(vis_sel) - grs(2).allblobs.omega(vis_sel) + 180, 360) - 180), ...
abs(mod(grs(1).allblobs.omega(vis_sel) - grs(2).allblobs.omega(vis_sel) + 180, 360) - 180) < grs(1).allblobs.lorentzfac(vis_sel), ...
grs(1).allblobs.lorentzfac(vis_sel), ...
acosd(abs(cosd(grs(1).allblobs.eta(vis_sel)))) < conf.min_eta, ...
acosd(abs(cosd(grs(2).allblobs.eta(vis_sel)))) < conf.min_eta, ...
]');
end
vol_size = [ ...
refgr_proj.vol_size_y, ...
refgr_proj.vol_size_x, ...
refgr_proj.vol_size_z ];
vol_half_size = vol_size / 2;
space_bbox = [ ...
refgr_proj.centerpix - vol_half_size * conf.rspace_oversize, ...
refgr_proj.centerpix + vol_half_size * conf.rspace_oversize ];
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
estim_space_bbox_pix = [floor(space_bbox(1:3)), ceil(space_bbox(4:6))];
bbox_size_pix = estim_space_bbox_pix(4:6) - estim_space_bbox_pix(1:3);
estim_space_bbox_mm = [ ...
gtGeoSam2Sam(estim_space_bbox_pix(1:3), p.recgeo, p.samgeo, false, false), ...
gtGeoSam2Sam(estim_space_bbox_pix(4:6), p.recgeo, p.samgeo, false, false) ];
bbox_size_mm = estim_space_bbox_mm(4:6) - estim_space_bbox_mm(1:3);
img_sizes = zeros(num_grains, 2);
for ii_g = num_grains:-1:1
sampler = GtOrientationSampling(p, grs(ii_g), 'verbose', conf.verbose);
r_vecs = sampler.guess_ODF_BB()';
estim_orient_bbox = [min(r_vecs, [], 1), max(r_vecs, [], 1)];
bbox_size_rod = estim_orient_bbox(4:6) - estim_orient_bbox(1:3);
% oversizing the orienation a bit
delta_bbox_size_rod = bbox_size_rod * 0.05;
estim_orient_bbox = estim_orient_bbox + [-delta_bbox_size_rod, delta_bbox_size_rod];
bbox_size_rod = estim_orient_bbox(4:6) - estim_orient_bbox(1:3);
bbox_size_deg = 2 * atand(bbox_size_rod);
if (conf.verbose)
fprintf('\n');
fprintf('Estimated spatial voxel BBox: [%3d, %3d, %3d] -> [%3d, %3d, %3d]\n', estim_space_bbox_pix);
fprintf(' BBox size: %3d, %3d, %3d (%f, %f, %f mm)\n', bbox_size_pix, bbox_size_mm);
fprintf(' Estimated orientation BBox: [%3.3f, %3.3f, %3.3f] -> [%3.3f, %3.3f, %3.3f]\n', estim_orient_bbox);
fprintf(' BBox size: %3.3f, %3.3f, %3.3f (deg)\n', bbox_size_deg);
fprintf('\n');
end
refor = struct( ...
'id', grs(ii_g).id, 'phaseid', grs(ii_g).phaseid, ...
'center', estim_space_bbox_mm(1:3) + bbox_size_mm / 2, ...
'R_vector', grs(ii_g).R_vector );
refor = gtCalculateGrain(refor, p);
refor.bb_ors = get_6D_space_corners(estim_space_bbox_mm, ...
bbox_size_mm, estim_orient_bbox, bbox_size_rod, refgr, p);
img_sizes(ii_g, :) = [size(grs(ii_g).proj.stack, 1), size(grs(ii_g).proj.stack, 3), ];
samp_ors(ii_g) = refor;
end
samp_ors = reorder_reflections(samp_ors, refgr);
[stackUSize, stackVSize] = get_stack_UV_size(cat(1, img_sizes), p, conf);
bool_ref_inc = false(size(refgr.allblobs.omega));
bool_ref_inc(ref_ondet(ref_included)) = true;
bool_ref_sel = false(size(refgr.allblobs.omega));
bool_ref_sel(ref_ondet(ref_included(ref_selected))) = true;
for ii_g = 1:num_grains
fprintf('Loading raw images for grain %d: ', ii_g)
if (ii_g > 1)
twingr = grs(ii_g);
twingr_proj = twingr.proj(conf.det_ind);
twin_ondet = twingr_proj.ondet;
twin_inc = twingr_proj.included;
twin_sel = twingr_proj.selected;
Nicola Vigano
committed
bool_twin_ondet = false(size(refgr.allblobs.omega));
bool_twin_ondet(twin_ondet) = true;
bool_twin_inc = false(size(refgr.allblobs.omega));
bool_twin_inc(twin_ondet(twin_inc)) = true;
bool_twin_sel = false(size(refgr.allblobs.omega));
bool_twin_sel(twin_ondet(twin_inc(twin_sel))) = true;
Nicola Vigano
committed
% Indices of reflections that exist on the detector for the
% twin, and were included for the parent
bool_test_inc = bool_ref_inc & bool_twin_ondet;
inc_ref_omegas = samp_ors(1).allblobs.omega(bool_test_inc);
inc_ref_Lfac = samp_ors(1).allblobs.lorentzfac(bool_test_inc);
inc_ref_etas = samp_ors(1).allblobs.eta(bool_test_inc);
inc_twin_omegas = samp_ors(ii_g).allblobs.omega(bool_test_inc);
inc_twin_etas = samp_ors(ii_g).allblobs.eta(bool_test_inc);
w_diff = abs(mod(inc_ref_omegas - inc_twin_omegas + 180, 360) - 180);
n_diff = abs(mod(inc_ref_etas - inc_twin_etas + 180, 360) - 180);
test_shared = (w_diff < inc_ref_Lfac) & (n_diff < 5);
Nicola Vigano
committed
% let's add the shared as included
shared_parent = false(size(refgr.allblobs.omega));
Nicola Vigano
committed
bool_test_inc_indx = find(bool_test_inc);
shared_parent(bool_test_inc_indx(test_shared)) = true;
Nicola Vigano
committed
bool_twin_inc_redundant = bool_twin_inc & shared_parent;
% So if it is shared, we use the parent's blob
bool_twin_inc = bool_twin_inc | shared_parent;
% If the parent enables it, we enable it as well
bool_twin_sel = bool_twin_sel | (shared_parent & bool_ref_sel);
Nicola Vigano
committed
% We now find the blobs in the old twin bl structure that are
% shared with the parent
shared_parent_inc_redundant = bool_twin_inc_redundant(twin_ondet(twin_inc));
twin_inc = find(bool_twin_inc(twin_ondet));
twin_sel = bool_twin_sel(twin_ondet(twin_inc));
shared_parent_pos = bool_ref_inc & shared_parent;
shared_parent_pos = shared_parent_pos(ref_ondet(ref_included));
% We cut it down to the size of the included
shared_parent = shared_parent(twin_ondet(twin_inc));
local_ondet = twin_ondet;
local_included = twin_inc;
local_selected = twin_sel;
else
local_ondet = ref_ondet;
local_included = ref_included;
local_selected = ref_selected;
end
num_blobs = numel(local_included);
if (ii_g == 1)
blobs = grs(ii_g).proj.bl;
else
blobs = gtFwdSimBlobDefinition('blob', num_blobs);
Nicola Vigano
committed
blobs(~shared_parent) = grs(ii_g).proj.bl(~shared_parent_inc_redundant);
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
% Using the same because it gives the same size/uvw-limits
blobs(shared_parent) = samp_ors(1).proj.bl(shared_parent_pos);
end
for ii_b = 1:num_blobs
num_chars = fprintf('%02d/%02d', ii_b, num_blobs);
blobs(ii_b) = include_blob(blobs(ii_b), stackUSize, stackVSize);
fprintf(repmat('\b', [1 num_chars]));
fprintf(repmat(' ', [1 num_chars]));
fprintf(repmat('\b', [1 num_chars]));
end
fprintf('Done.\n')
spots = arrayfun(@(x){sum(x.intm, 3)}, blobs);
spots = permute(cat(3, spots{:}), [1 3 2]);
proj.centerpix = estim_space_bbox_pix(1:3) + bbox_size_pix / 2;
proj.bl = blobs;
proj.stack = spots;
vol_size = bbox_size_pix + mean(bbox_size_pix) * 0.3;
proj.vol_size_x = vol_size(2);
proj.vol_size_y = vol_size(1);
proj.vol_size_z = vol_size(3);
proj.ondet = local_ondet;
proj.included = local_included;
proj.selected = local_selected;
if (ii_g > 1)
proj.shared_bls_parent = zeros(size(shared_parent));
proj.shared_bls_parent(shared_parent) = find(shared_parent_pos);
samp_ors(ii_g).proj(conf.det_ind).shared_bls_twins(shared_parent_pos) = true;
else
proj.shared_bls_twins = false(num_blobs, 1);
end
samp_ors(ii_g).proj(conf.det_ind) = proj;
Nicola Vigano
committed
if (conf.check_spots && ii_g > 1)
samp_ors(ii_g).proj.selected = gtGuiGrainMontage(samp_ors(ii_g), [], true);
end
end
if (conf.verbose)
f = figure();
ax = axes('parent', f);
hold(ax, 'on');
for ii_g = 1:num_grains
gt6DPlotOrientationBBox(ax, cat(1, samp_ors(ii_g).bb_ors(:).R_vector));
scatter3(ax, grs(ii_g).R_vector(1), grs(ii_g).R_vector(2), grs(ii_g).R_vector(3), 30);
end
hold(ax, 'off');
drawnow();
end
if (conf.save)
str_ids = sprintf('_%04d', grs_list);
grain_filename = fullfile(p.acq.dir, '4_grains', ...
sprintf('phase_%02d', phase_id), ...
sprintf('grain_cluster%s.mat', str_ids));
fprintf('Saving the cluster file "%s"..', grain_filename)
save(grain_filename, 'samp_ors', '-v7.3');
fprintf('\b\b: Done.\n')
fprintf('Saving to sample.mat..')
sample = GtSample.loadFromFile();
sample.phases{phase_id}.clusters(end+1) = ...
GtPhase.makeCluster(grs_list, cat(1, samp_ors(:).R_vector), estim_space_bbox_pix, [], 1);
sample.saveToFile();
fprintf('\b\b: Done.\n')
end
end
function orients = get_6D_space_corners(estim_space_bbox_mm, bbox_size_mm, estim_orient_bbox, bbox_size_rod, refgr, p)
orients = cell(2, 2, 2, 2, 2, 2);
for ii_x1 = 0:1
base_x1 = estim_space_bbox_mm(1:3) + [ii_x1 * bbox_size_mm(1), 0, 0];
for ii_x2 = 0:1
base_x2 = base_x1 + [0, ii_x2 * bbox_size_mm(2), 0];
for ii_x3 = 0:1
base_x3 = base_x2 + [0, 0, ii_x3 * bbox_size_mm(3)];
for ii_r1 = 0:1
base_r1 = estim_orient_bbox(1:3) + [ii_r1 * bbox_size_rod(1), 0, 0];
for ii_r2 = 0:1
base_r2 = base_r1 + [0, ii_r2 * bbox_size_rod(2), 0];
for ii_r3 = 0:1
base_r3 = base_r2 + [0, 0, ii_r3 * bbox_size_rod(3)];
% fprintf(' center: %f, %f, %f, r_vec: %f, %f, %f\n', base_x3, base_r3)
orients{ii_r3+1, ii_r2+1, ii_r1+1, ii_x3+1, ii_x2+1, ii_x1+1} = struct( ...
'id', refgr.id, 'phaseid', refgr.phaseid, ...
'center', base_x3, 'R_vector', base_r3);
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
end
end
end
end
end
end
orients = gtCalculateGrain_p(orients, p);
orients = [orients{:}];
orients = reorder_reflections(orients, refgr);
end
function orients = reorder_reflections(orients, refgr)
refgr_omind = refgr.allblobs.omind;
refgr_omind_ind = [ ...
find(refgr_omind == 1); find(refgr_omind == 2); ...
find(refgr_omind == 3); find(refgr_omind == 4) ];
for ii_g = 1:numel(orients)
orients(ii_g) = gtGrainAllblobsFilterOrder(orients(ii_g), refgr_omind_ind);
end
end
function [u_size, v_size] = get_stack_UV_size(img_sizes, p, conf)
max_img_sizes = [max(img_sizes(:, 1)), max(img_sizes(:, 2))];
stack_imgs_oversize = min(p.fsim.oversize, conf.stack_oversize);
u_size = round(max_img_sizes(1) * stack_imgs_oversize);
v_size = round(max_img_sizes(2) * stack_imgs_oversize);
if (conf.verbose)
fprintf('\n');
fprintf(' Maximum images size: [%3d, %3d]\n', max_img_sizes);
fprintf('Stack images size (oversize: %1.2f): [%3d, %3d]\n', stack_imgs_oversize, u_size, v_size);
fprintf('\n');
end
end
function blob = include_blob(blob, stackUSize, stackVSize)
new_bbsize = [stackUSize, stackVSize, blob.bbsize(3)];
shift = floor((new_bbsize - blob.bbsize) / 2);
new_bbuim = blob.bbuim - shift(1);
new_bbuim(2) = new_bbuim(1) + new_bbsize(1) - 1;
new_bbvim = blob.bbvim - shift(2);
new_bbvim(2) = new_bbvim(1) + new_bbsize(2) - 1;
blob.bbuim = new_bbuim;
blob.bbvim = new_bbvim;
blob.bbsize = new_bbsize;
blob.intm = gtPlaceSubVolume(zeros(new_bbsize, 'like', blob.intm), blob.intm, shift);
blob.mask = logical(gtPlaceSubVolume(zeros(new_bbsize, 'uint8'), blob.mask, shift));
end
Nicola Vigano
committed
function gr = twin_fwd_sim_to_gr(ref_gr, twin)
gr = ref_gr;
gr.R_vector = twin.R_vector;
gr.allblobs = twin.allblobs;
gr.proj = twin.proj;
% And the useless stuff
gr.check = twin.fwd_sim.check;
gr.flag = twin.fwd_sim.flag;
gr.spotid = twin.fwd_sim.spotid;
% only included spots info
gr.cands = twin.fwd_sim.cands;
gr.likelihood = twin.fwd_sim.likelihood;
gr.difspotID = twin.fwd_sim.difspotID;
gr.uv_exp = twin.fwd_sim.uvw(:, 1:2);
gr.om_exp = twin.fwd_sim.uvw(:, 3);
gr.bb_exp = twin.fwd_sim.bb;
gr.intensity = twin.fwd_sim.intensity;
end