gtForwardSimulation/gtCalculateGrain: fixed grain center (pixels vs. mms)

Added new field centerpix for the purpouse!

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

4_grains/gtCalculateGrain.m

@@ -69,6 +69,7 @@ end
 % Rodrigues vector (grain orientation)
 R_vec = grain.R_vector;
 phase = grain.phaseid;
+gcsam = grain.center';
 
 cryst  = parameters.cryst(phase);
 labgeo = parameters.labgeo;
@@ -88,10 +89,7 @@ if (~isfield(labgeo, 'detnorm'))
     labgeo.detnorm = gtMathsCross(labgeo.detdiru, labgeo.detdirv);
 end
 
-pixelsize = mean([labgeo.pixelsizeu labgeo.pixelsizev]);
-
 uvorig  = [labgeo.detrefu, labgeo.detrefv]';
-csam    = (grain.center .* pixelsize)';
 omstep  = 180/parameters.acq.nproj;
 
 %%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -106,7 +104,6 @@ dsp_ref = gtCrystDSpacing(double(cryst.hkl(:, cryst.thetatypesp)), Bmat);
 
 % Plane normals in Cartesian CRYSTAL coordinates
 pl_cry = gtCrystHKL2Cartesian(double(cryst.hklsp), Bmat)';
-num_plane_normals = size(pl_cry, 1);
 
 % Compute all orientation matrices g (Reminder: Vc = g . Vs)
 g = gtMathsRod2OriMat(R_vec');
@@ -127,6 +124,7 @@ if (app.usestrain)
     [pl_samd, drel] = gtStrainPlaneNormals(pl_sam', defT); % unit vector
 else
     pl_samd = pl_sam';
+    num_plane_normals = size(pl_sam, 1);
     drel = ones(1, num_plane_normals);
 end
 
@@ -244,7 +242,7 @@ eta = gtGeoEtaFromDiffVec(dveclab, labgeo);
 dvecsam = gtGeoLab2Sam(dveclab, om, labgeo, samgeo, true);
 
 % u,v,w coordinates on the detector
-csam_v = csam(:, ones(1, size(rot, 3)));
+csam_v = gcsam(:, ones(1, size(rot, 3)));
 uvw = gtFedPredictUVWMultiple(rot, dveclab', csam_v, labgeo.detrefpos', ...
     labgeo.detnorm', Qdet, uvorig, om', omstep)';
 

5_reconstruction/gtForwardSimulate_v2.m

@@ -143,9 +143,6 @@ for n = first : last
 
     gr = grain{n}; % output from INDEXTER
 
-    % graincenter contains COM coordinates in pixels
-    graincenter_pix = gr.center ./ pixelsize;
-
     %%% Predict spot positions
     gr = gtCalculateGrain(gr, parameters);
     u  = gr.allblobs.uv(:, 1);
@@ -403,12 +400,15 @@ for n = first : last
     %%% Now fill in the information for the output
     %  remove reflections for which no segmented spot has been found
 
-    out.selected  = selected;  % vector of logicals with pre-selected spots active
-    out.included  = included; % This is needed by Orientation Reconstruction
     out.id        = gr.id;
     out.phaseid   = phaseID;
     out.R_vector  = gr.R_vector;
-    out.center    = graincenter_pix;
+    out.center    = gr.center;
+    % We should use gtGeoSam2Sam with samgeo and recgeo here!
+    out.centerpix = gr.center ./ pixelsize;
+
+    out.selected  = selected;  % vector of logicals with pre-selected spots active
+    out.included  = included; % This is needed by Orientation Reconstruction
     out.ondet     = onDet;
     out.allblobs  = gr.allblobs;
     out.full      = full;
@@ -449,9 +449,9 @@ for n = first : last
     out.vol             = [];
     out.segbb           = zeros(1,6);
     out.threshold       = 0;
-    out.shift(1) = round((acq.bb(3)-out.proj.vol_size_x)/2 + out.center(1)) + 1;
-    out.shift(2) = round((acq.bb(3)-out.proj.vol_size_y)/2 + out.center(2)) + 1;
-    out.shift(3) = round((acq.bb(4)-out.proj.vol_size_z)/2 + out.center(3)) + 1;
+    out.shift(1) = round((acq.bb(3)-out.proj.vol_size_x)/2 + out.centerpix(1)) + 1;
+    out.shift(2) = round((acq.bb(3)-out.proj.vol_size_y)/2 + out.centerpix(2)) + 1;
+    out.shift(3) = round((acq.bb(4)-out.proj.vol_size_z)/2 + out.centerpix(3)) + 1;
 
     if (fsim.save_grain)
         fprintf(' (Done).\n - Saving ouput..');
@@ -724,7 +724,7 @@ function candidateIDs = findCandidates(spotsProps, reflPos, segmentedSpots)
     cand_area = segmentedSpots.BoundingBoxXsize(image_ok).*segmentedSpots.BoundingBoxYsize(image_ok);
     area_ok   = find(cand_area >= 0.5 * spotsProps.XsizeMin .* spotsProps.YsizeMin);
     image_area_ok  = image_ok(area_ok);
-    
+
     Xdist = abs(segmentedSpots.BoundingBoxXorigin(image_area_ok) ...
         + 0.5 * segmentedSpots.BoundingBoxXsize(image_area_ok) - reflPos.u);
     Ydist = abs(segmentedSpots.BoundingBoxYorigin(image_area_ok) ...