gtDefSyntheticGrainCreate: simplified grain creation

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

zUtil_Deformation/gtDefSyntheticGrainCreate.m

@@ -53,27 +53,12 @@ parameters.labgeo = labgeo;
 %%% Load grain parameters
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-grsize = fedpars.grainsize; % x, y, z number of voxels in grain
-gvsize = fedpars.gvsize;    % x, y, z number of voxels in an element
-
-% Grain volume
-grvol = gtFedTestGenerateGrainVolume(grsize);
-fedpars.grainvol = grvol;
-
-% Determine number of volume elements along X, Y, Z, based on grain volume
-% and element size:
-[grenv(1), grenv(2), grenv(3)] = size(grvol);
-ngv = ceil(grenv ./ gvsize);
-fedpars.ngv = ngv;
-
-% Size of grain envelope
-grenv = ngv .* gvsize;
-fedpars.grainenv = grenv;
-
-% Padded grain volume
-grvol_pad = false(grenv);
-grvol_pad(1:size(grvol, 1), 1:size(grvol, 2), 1:size(grvol, 3)) = grvol;
-grvol     = grvol_pad;
+if (~exist('dmvol', 'var') || isempty(dmvol))
+    grsize = fedpars.grainsize; % x, y, z number of voxels in grain
+else
+    grsize = [size(dmvol, 1), size(dmvol, 2), size(dmvol, 3)];
+end
+ngv = grsize;
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%% Volume sampling and volume elements
@@ -82,8 +67,8 @@ grvol     = grvol_pad;
 fprintf('\b\b: Done in %g seconds.\nAllocating memory for elements..', toc(c))
 c = tic();
 
-nv    = ngv(1)*ngv(2)*ngv(3);          % total number of vol elements
-svvol = gvsize(1)*gvsize(2)*gvsize(3); % size of sampling volume
+% total number of vol elements
+nv = ngv(1) * ngv(2) * ngv(3);
 
 nd = sum(fedpars.dcomps);
 
@@ -157,22 +142,18 @@ end
 fprintf('\b\b: Done in %g seconds.\nComputing voxel centers..', toc(c))
 c = tic();
 
+fedpars.vox000sam = grain.center - (ngv / 2 + 0.5) .* recgeo(1).voxsize;
+
 % Elements' positions in the sample volume in mm!
 [gv.cs, gv.d] = compute_voxel_centers(gv, dmvol, fedpars, recgeo(1));
-
-if (isfield(fedpars, 'is_real_data') && fedpars.is_real_data)
-    gc = grain.center;
-else
-    % Grain center position in sample in mm!
-    gc = fedpars.vox000sam + (grsize/2 + [0.5, 0.5, 0.5]) .* recgeo(1).voxsize;
-    grain.center = gc;
-    fedpars.graincentsam = gc;
-end
+size(gv.cs)
+size(gv.d)
+size(gv.used_ind)
 
 fprintf('\b\b: Done in %g seconds.\n', toc(c))
 
-if ~all(gradok)
-    disp('Gradients are higher than limit for components:')
+if (~all(gradok))
+    disp('Gradients are higher than the limits for components:')
     disp(find(~gradok)')
 end
 
@@ -187,12 +168,10 @@ disp(min(dmvol(:)))
 
 % Cleaning grain struct
 if (isfield(grain, 'id'))
-    grain = struct('id', grain.id, 'phaseid', grain.phaseid, ...
-        'center', grain.center, 'R_vector', grain.R_vector);
-else
-    grain = struct('id', 1, 'phaseid', grain.phaseid, ...
-        'center', grain.center, 'R_vector', grain.R_vector);
+    grain.id = 1;
 end
+grain = struct('id', grain.id, 'phaseid', grain.phaseid, ...
+    'center', grain.center, 'R_vector', grain.R_vector);
 
 grain = gtCalculateGrain(grain, parameters);
 
@@ -280,6 +259,8 @@ function [gvcs, gvd] = compute_voxel_centers(gv, dmvol, fedpars, recgeo)
     nv = size(gv.ind, 2);
     nv_ones = ones(nv, 1);
 
+    base_gvcs = fedpars.vox000sam(nv_ones, :)' + gv.ind .* recgeo.voxsize(nv_ones, :)';
+
     if (isfield(fedpars, 'volume_super_sampling') ...
             && (fedpars.volume_super_sampling > 1))
 
@@ -307,9 +288,7 @@ function [gvcs, gvd] = compute_voxel_centers(gv, dmvol, fedpars, recgeo)
 
                     offset_subvoxel = [ss_x, ss_y, ss_z];
 
-                    gvcs(:, :, num_subvoxels-counter) = ...
-                        fedpars.vox000sam(nv_ones, :)' ...
-                        + gv.ind .* recgeo.voxsize(nv_ones, :)' ...
+                    gvcs(:, :, num_subvoxels-counter) = base_gvcs ...
                         + offset_subvoxel(nv_ones, :)';
 
                     offset_pos = offset_subvoxel ./ recgeo.voxsize;
@@ -370,7 +349,7 @@ function [gvcs, gvd] = compute_voxel_centers(gv, dmvol, fedpars, recgeo)
             end
         end
     else
-        gvcs = fedpars.vox000sam(nv_ones, :)' + gv.ind .* recgeo.voxsize(nv_ones, :)';
+        gvcs = base_gvcs;
         gvd = gv.d;
     end
 end