Topotomo importing: initial function that creates the grain.proj structure

It still lacks the projection geometry

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

4_grains/gtGrainImportTopotomoScan.m

+function gr = gtGrainImportTopotomoScan(gr, p, det_ind, data_type)
+    if (~exist('p', 'var') || isempty(p))
+        p = gtLoadParameters();
+    end
+    if (~exist('data_type', 'var') || isempty(data_type))
+        data_type = 'single';
+    end
+
+    fprintf('Importing Topotomo acquisition into grain %d (phase: %d):\n', gr.id, gr.phaseid)
+    gr = gtCalculateGrain(gr, p, 'det_ind', det_ind);
+
+    fprintf(' - Loading blobs..')
+    c = tic();
+    phstep = gtAcqGetPhiStep(p, det_ind);
+    nbl = gtAcqTotNumberOfImages(p, det_ind);
+    nph = p.acq(det_ind).nproj_basetilt;
+    b_size = [p.acq(det_ind).xdet, p.acq(det_ind).ydet, nph];
+    roi_u = [1, p.acq(det_ind).xdet] + p.acq(det_ind).detroi_u_off;
+    roi_v = [1, p.acq(det_ind).ydet] + p.acq(det_ind).detroi_v_off;
+    roi_p = [1, nph] + p.acq(det_ind).range_basetilt(1) / phstep;
+
+    base_dir = p.acq(det_ind).collection_dir;
+    base_name = p.acq(det_ind).name;
+    % Creating blobs
+    bl = gtFwdSimBlobDefinition('blob_topo', nbl);
+    filen_name = fullfile(base_dir, [base_name, '0000.edf']);
+    info = edf_info(filen_name);
+    fprintf('\b\b (initialization done in %g s): ', toc(c))
+    c = tic();
+    for ii_b = 1:nbl
+        num_chars = fprintf('%03d/%03d', ii_b, nbl);
+
+        b_vol = zeros(b_size, data_type);
+        for ii_p = 1:nph
+            filen_name = fullfile(base_dir, [base_name, sprintf('%04d.edf', (ii_b - 1) * nph + ii_p)]);
+            b_vol(:, :, ii_p) = edf_read(filen_name, [], 'nodisp', info);
+        end
+        bl(ii_b).intm = b_vol;
+        bl(ii_b).mask = true(b_size);
+        bl(ii_b).bbuim = roi_u;
+        bl(ii_b).bbvim = roi_v;
+        bl(ii_b).bbpim = roi_p;
+        bl(ii_b).bbsize = b_size;
+        bl(ii_b).mbbu = bl(ii_b).bbuim;
+        bl(ii_b).mbbv = bl(ii_b).bbvim;
+        bl(ii_b).mbbp = bl(ii_b).bbpim;
+        bl(ii_b).mbbsize = bl(ii_b).bbsize;
+        bl(ii_b).intensity = sum(b_vol);
+
+        fprintf(repmat('\b', [1, num_chars]))
+        fprintf(repmat(' ', [1, num_chars]))
+        fprintf(repmat('\b', [1, num_chars]))
+    end
+    fprintf('\b\b: Done in %g seconds.\n', toc(c))
+
+    fprintf(' - Detecting phi index (phind)..')
+    c = tic();
+    lims_bt = p.acq(det_ind).range_basetilt;
+    within_lims = gr.allblobs(det_ind).phi > lims_bt(1) & gr.allblobs(det_ind).phi < lims_bt(2);
+    phinds = gr.allblobs(det_ind).phind(within_lims);
+    phinds_counts = histcounts(phinds, (1:5)-0.5);
+    [~, phind] = max(phinds_counts);
+
+    selectedph = gr.allblobs(det_ind).phind == phind;
+    fprintf('\b\b: Done in %g seconds.\n', toc(c))
+    fprintf('   + Chosen phind: %d (counts: [%s])\n', phind, sprintf(' %d', phinds_counts))
+
+    fprintf(' - Building "proj" structure..')
+    c = tic();
+    stack = arrayfun(@(x){sum(x.intm, 3)}, bl);
+    stack = cat(3, stack{:});
+
+    vol_size = [gr.proj(1).vol_size_y, gr.proj(1).vol_size_x, gr.proj(1).vol_size_z];
+    vol_size = gtGeoSam2Sam(vol_size, p.recgeo(1), p.recgeo(det_ind), 1);
+
+    proj = gtFwdSimProjDefinition('fwd_sim');
+
+    proj.stack = permute(stack, [1 3 2]);
+    proj.num_cols = size(proj.stack, 3);
+    proj.num_rows = size(proj.stack, 1);
+    proj.bl = bl;
+    proj.vol_size_x = vol_size(2);
+    proj.vol_size_y = vol_size(1);
+    proj.vol_size_z = vol_size(3);
+    proj.centerpix = gtGeoSam2Sam(gr.center, p.samgeo, p.recgeo(det_ind), 0);
+    proj.ondet = find(selectedph);
+    proj.included = (1:nbl)';
+    proj.selected = true(size(proj.ondet));
+
+%     sel = proj.ondet(proj.included(proj.selected));
+%     dvecsam = gr.allblobs(det_ind).dvecsam(sel, :);
+%     rot_s2l = gr.allblobs(det_ind).rot_s2l(:, :, sel);
+%     proj.geom = gtGeoProjForReconstruction(dvecsam, ...
+%         rot_s2l, gr.center, bbpos_det_grain, [], ...
+%         p.detgeo(det_ind), labgeo, samgeo, p.recgeo(det_ind), 'ASTRA_grain');
+
+    fprintf('\b\b: Done in %g seconds.\n', toc(c))
+
+    try
+        gr.proj(det_ind) = proj;
+    catch mexc
+        if gtCheckExceptionType(mexc, 'MATLAB:heterogeneousStrucAssignment')
+            fprintf(' - Updating legacy "proj" structure..')
+            c = tic();
+            % We now handle the fact that proj can be a legacy datastructure
+            base_proj = gtFwdSimProjDefinition('fwd_sim');
+            for ii_p = numel(gr.proj):-1:1
+                new_proj(ii_p) = gtAddMatFile(base_proj, gr.proj(ii_p), true);
+            end
+            gr.proj = new_proj;
+
+            gr.proj(det_ind) = proj;
+            fprintf('\b\b: Done in %g seconds.\n', toc(c))
+        end
+    end
+end