Phantom generation: small cleanup

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

7_fed2/gtFedTestLoadData_v2.m

@@ -180,62 +180,48 @@ for ii = 1:nbl
     bl(ii).thetatype = grain.allblobs.thetatype(lbl(ii));
 
     % Initial strained signed plane normal in SAMPLE reference
-    bl(ii).plsam  = grain.allblobs.pl(lbl(ii), :);
+    bl(ii).plsam = grain.allblobs.pl(lbl(ii), :);
 
     % Initial strained non-signed plane normal in SAMPLE reference
     bl(ii).plorig = grain.allblobs.plorig(lbl(ii), :);
 
     % Initial strained signed plane normal in the diffraction position in
     % LAB reference
-    bl(ii).pl0    = grain.allblobs.pllab(lbl(ii), :);
+    bl(ii).pl0 = grain.allblobs.pllab(lbl(ii), :);
 
     % Initial strained diffraction angles and rotation matrix
     bl(ii).th0    = grain.allblobs.theta(lbl(ii));     % initial theta
-    bl(ii).sinth0 = sind(bl(ii).th0);                  % initial sin(theta)
+    bl(ii).sinth0 = grain.allblobs.sintheta(lbl(ii));  % initial sin(theta)
     bl(ii).eta0   = grain.allblobs.eta(lbl(ii));       % initial eta
     bl(ii).om0    = grain.allblobs.omega(lbl(ii));     % initial omega
     bl(ii).omind  = grain.allblobs.omind(lbl(ii));     % initial omega index (1..4)
-    bl(ii).S0     = gtFedRotationTensor(bl(ii).om0, lab.rotcomp); % initial rotation tensor
+    bl(ii).S0     = grain.allblobs.srot(:, :, lbl(ii)); % initial rotation tensor
 
     % Initial strained diffraction vector in LAB reference
-    bl(ii).dvec0   = grain.allblobs.dvec(lbl(ii), :);
+    bl(ii).dvec0 = grain.allblobs.dvec(lbl(ii), :);
 
     % Initial strained diffraction vector in SAMPLE reference
     bl(ii).dvecsam = grain.allblobs.dvecsam(lbl(ii), :);
 
+    gc_lab = gtGeoSam2Lab(gc', bl(ii).S0', parameters.labgeo, parameters.samgeo, false);
 
-    % Calculate initial mean diffraction parameters of blob based on grain center
-
-    % Initial centroid u in image
-    bl(ii).u0im = gtFedPredictUVW(bl(ii).S0, bl(ii).dvec0', gc, ...
+    % Initial centroid in image
+    bl(ii).u0im = gtFedPredictUVW(eye(3), bl(ii).dvec0', gc_lab', ...
         lab.detrefpos, lab.detnorm, lab.Qdet, lab.detrefuv, bl(ii).om0, ...
         lab.omstep)';
 
     if isempty(bl(ii).u0im)
         blobsdiffaway = [blobsdiffaway, ii];
     else
-        bl(ii).dfac0  = gtFedPredictDiffFac(bl(ii).S0, bl(ii).dvec0', ...
-            gc,lab.detrefpos, lab.detnorm);
-
-        % Mean U tensor of blob is calculated from mean parameters
-        %bl(ii).U0   = gtFedDerUmatrix(bl(ii).plsam',bl(ii).om0,bl(ii).S0,lab.rotcomp,gc,...
-        %    bl(ii).sinth0,bl(ii).dvec0',bl(ii).dfac0,lab.detrefpos,lab.detnorm,lab.Qdet,lab.beamdir,lab.omstep,fedpars.beamchroma);
-
-        % For the moment, use Umean instead of Umax to estimate extreme peak
-        % shifts:
-
-        % u limit for all elements (u displacement at extreme strain combination);
-        %bl(ii).ulim   = ceil(abs(bl(ii).U0)*fedpars.dlim')';
-
-        % Padding around blob volume according to ulim
-        %bl(ii).addbl  = bl(ii).ulim + fedpars.blobsizeadd;
-        bl(ii).addbl  = fedpars.blobsizeadd;
+        % Padding around blob volume
+        bl(ii).addbl = fedpars.blobsizeadd;
 
         % Check if blob hits the detector area
         uok = (bl(ii).u0im(1) < parameters.acq.xdet) & (bl(ii).u0im(1) > 1);
         vok = (bl(ii).u0im(2) < parameters.acq.ydet) & (bl(ii).u0im(2) > 1);
-        if uok && vok
-          blobondet = [blobondet, ii];
+
+        if (uok && vok)
+            blobondet = [blobondet, ii];
         end
     end
 end
@@ -299,18 +285,11 @@ for ii = 1:nbl
                      round( (bl(ii).bbwim(1) + bl(ii).bbwim(2))/2 ) ];
 
     % Projection parameters
-    bl(ii).proj_geom = gtGeoProjForReconstruction(bl(ii).dvecsam, ...
-                       bl(ii).om0, gc', ...
-                       [bl(ii).bbuim(1) bl(ii).bbvim(1) bl(ii).bbsize(1:2)], ...
-                       [], parameters.labgeo, parameters.samgeo, ...
-                       parameters.recgeo, 'ASTRA_grain');
-
-
-    %bl(ii).int0 = zeros(bl(ii).bbsize);
-
-    %bl(ii).int  = zeros(bl(ii).bbsize);
-    %bl(ii).intd = zeros(bl(ii).bbsize);
-    %bl(ii).intm = zeros(bl(ii).bbsize);
+    bl(ii).proj_geom = gtGeoProjForReconstruction( ...
+        bl(ii).dvecsam, bl(ii).om0, gc', ...
+        [bl(ii).bbuim(1) bl(ii).bbvim(1) bl(ii).bbsize(1:2)], ...
+        [], parameters.labgeo, parameters.samgeo, parameters.recgeo, ...
+        'ASTRA_grain');
 
     %%% Needed for FED
 %     bl(ii).int  = zeros(bl(ii).bbsize, fedpars.bltype);
@@ -320,12 +299,6 @@ for ii = 1:nbl
 
     % u in blob = u in image minus u at blob origin
     bl(ii).u0bl = bl(ii).u0im - bl(ii).bbor;
-
-    % To save memory vs. simlicity
-    %    bl(i).int0 = single(zeros(bl(i).bbsize));
-    %    bl(i).int  = single(zeros(bl(i).bbsize));
-    %    bl(i).intd = single(zeros(bl(i).bbsize));
-    %    bl(i).intm = single(zeros(bl(i).bbsize));
 end
 
 toc
@@ -339,12 +312,7 @@ tic
 
 gvcs = gv.cs;
 
-% !!! parfor
 for ii = nbl:-1:1
-    %gvb(ii).uc0im   = zeros(size(gvcs));
-    %gvb(ii).ucim    = zeros(size(gvcs));
-    %gvb(ii).ucbl    = zeros(size(gvcs));
-
     gvb(ii).uc0bl   = zeros(size(gvcs));
     gvb(ii).uc      = zeros(size(gvcs));
     gvb(ii).ucd     = zeros(size(gvcs));
@@ -353,11 +321,6 @@ for ii = nbl:-1:1
     gvb(ii).ucbl8ini= zeros(1, size(gvcs,2)*8);
     gvb(ii).ucbl8ok = 0;
 
-    %gvb(ii).uc0bl   = zeros(size(gvcs),'single');
-    %gvb(ii).uc      = zeros(size(gvcs),'single');
-    %gvb(ii).ucblmod = zeros(size(gvcs),'single');
-    %gvb(ii).ucbl8   = zeros([size(gvcs,1), size(gvcs,2)*8],'single');
-
     % Expand rotation tensors
     subs = {(1:3)', (1:3)', ones(1,nv)};
     rot  = bl(ii).S0(subs{:});
@@ -376,16 +339,6 @@ for ii = nbl:-1:1
     % Initial u of gv centre relative to blob origin
     bbor = bl(ii).bbor';
     gvb(ii).uc0bl = uc0im - bbor(:,ones(1,nv));
-
-    % Actual u of gv centre relative to blob origin
-    %gvb(ii).ucbl = gvb(ii).uc0bl;
-
-    % Actual peak shift of gv centre
-    %gvb(ii).uc = zeros(3,nv);
-
-    % Actual plane normal
-    %plsam = bl(ii).plsam';
-    %gv.plsam(:,:,ii) = plsam(:,ones(1,nv));
 end
 
 toc
@@ -487,43 +440,6 @@ for ii = 1:nbl
     fedpars.blobnextdet(ii) = abs(det(bl(ii).U0*bl(nii).U0'));
 end
 
-% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-% %%% Compute initial diffraction patterns
-% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%
-% disp('Computing initial blob intensities...')
-% gv.d(:) = 0;
-% [gvb, bl] = gtFedFwdProjExact(gv, gvb, bl, fedpars);
-%
-% for ii = 1:nbl
-%     %bl(ii).int0 = bl(ii).int;
-%     gvb(ii).ucd(:) = 0;
-% end
-%
-% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-% %%% Compute strain derivatives
-% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%
-% disp('Computing uvw derivatives...')
-%
-% % only adds gvb.U field
-% gvb = gtFedDerUmatrixExact2(gv, gvb, bl, fedpars);
-% 
-% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-% %%% Average projection matrices
-% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%
-% % Mean U matrix of blob
-% for ii = 1:nbl
-%     bl(ii).U0 = reshape(mean(gvb(ii).U, 2), 3, []);
-% end
-%
-% fedpars.U0  = vertcat(bl(:).U0);
-% fedpars.U0w = fedpars.U0(3:3:end,:);
-%
-% fedpars.T0  = pinv(fedpars.U0);
-% fedpars.T0w = pinv(fedpars.U0w);
-
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%% Transform into single
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -547,6 +463,4 @@ end
 disp('True blob dimensions:')
 disp(vertcat(bl(:).mbbsize))
 
-% beep, pause(1), beep, pause(1), beep
-
 end % of main function