fixing a coordinate bug (detposcorner > pixels) in gtForwardSimulation_v2 -...

fixing a coordinate bug (detposcorner > pixels) in gtForwardSimulation_v2  - load tables at beginning

git-svn-id: https://svn.code.sf.net/p/dct/code/trunk@554 4c865b51-4357-4376-afb4-474e03ccb993

2_difspot/gtDBBlob2SpotTable.m

@@ -158,14 +158,17 @@ for difspotID = first:last
     % Write edf file on disk
     if flag_writeedf
         
-        % Read full edf-s, create summed image
-        im = zeros(bb(5),bb(4));
-        
-        for fullID = bb(3):bb(3)+bb(6)-1
-            im = im + edf_read(sprintf('1_preprocessing/full/full%04d.edf', fullID),...
-                bb([1 2 4 5]));
+        if any(strcmp(seg.difspotmask,{'none', 'blob2D'}))
+            % Read full edf-s, create summed image
+            im = zeros(bb(5),bb(4));
+
+            for fullID = bb(3):bb(3)+bb(6)-1
+                im = im + edf_read(sprintf('1_preprocessing/full/full%04d.edf', fullID),...
+                    bb([1 2 4 5]));
+            end
         end
         
+        
         % Check minimum area condition
         if isempty(thr_area)
             mask2D = [];

2_difspot/gtPostSegmentation.m

@@ -74,8 +74,7 @@ if test~=0
         overwrite_flag=true;
         gtDBCreateDifspotTable(parameters.acq.name, overwrite_flag);
     else
-        disp('Keeping difspot data.  Quitting gtPostSegmentation')
-        return
+        disp('Keeping difspot data')
     end
 end
         

5_reconstruction/gtForwardSimulate_v2.m

@@ -67,12 +67,10 @@ phases_num    = length(parameters.cryst);    % number of phases in the material
 
 cryst         = parameters.cryst(phaseID);
 grainfilename = fullfile(acq.dir, '4_grains', sprintf('phase_%02d', phaseID), 'index.mat');
+fprintf('\nLoading indexing results...\n');
 grains        = load(grainfilename);
 grain_num     = length(grains.grain);
 
-name          = acq.name;
-difspottable  = [name,'difspot'];
-
 
 %%% this section should be replaced with a function dealing with different
 % versions of parameters names:  fed, v1, v2...
@@ -83,14 +81,53 @@ detdirv0     = labgeo.detdirv;
 rotdir       = labgeo.rotdir;
 beamdir      = labgeo.beamdir;
 detrefpos    = labgeo.detrefpos / acq.pixelsize;   % detector reference position in pixels
-detposcorner = labgeo.detorig;  % the (0,0) corner of the detector in Lab system (edge of detector pixel)
+detposcorner = labgeo.detorig / acq.pixelsize;  % the (0,0) corner of the detector in Lab system (edge of detector pixel)
 
 rotcomp      = gtMathsRotationMatrixComp(rotdir,'row');
 omstep       = 180 / acq.nproj;
 imagerange   = round(fsim.omegarange/omstep);   % number of images summed up when searching raw data for a missing spot
 
-for n = first : last
 
+fprintf('\nLoading reflections from database...\n');
+difspottable = [parameters.acq.name,'difspot'];
+
+% Load difspot data in 'allspots'  /  spotpair data in 'allpairs'
+mymcmd = sprintf(['SELECT difspotID, CentroidX, CentroidY, CentroidImage, '...
+    'BoundingBoxXsize, BoundingBoxYsize, BoundingBoxXorigin, BoundingBoxYorigin '...
+    'FROM %s'],difspottable);
+
+[difspotIDs, CentroidXs, CentroidYs, CentroidImages, BoundingBoxXsizes, ...
+    BoundingBoxYsizes, BoundingBoxXorigins, BoundingBoxYorigins] = mym(mymcmd);
+
+numspots           = max(difspotIDs);
+difspotID          = NaN(numspots, 1);
+CentroidX          = NaN(numspots, 1);
+CentroidY          = NaN(numspots, 1);
+CentroidImage      = NaN(numspots, 1);
+BoundingBoxXsize   = NaN(numspots, 1);
+BoundingBoxYsize   = NaN(numspots, 1);
+BoundingBoxXorigin = NaN(numspots, 1);
+BoundingBoxYorigin = NaN(numspots, 1);
+
+for ii = 1: length(difspotIDs)
+    jj = difspotIDs(ii);
+    difspotID(jj)          = jj; 
+    CentroidX(jj)          = CentroidXs(ii);
+    CentroidY(jj)          = CentroidYs(ii);
+    CentroidImage(jj)      = CentroidImages(ii);
+    BoundingBoxXsize(jj)   = BoundingBoxXsizes(ii);
+    BoundingBoxYsize(jj)   = BoundingBoxYsizes(ii);
+    BoundingBoxXorigin(jj) = BoundingBoxXorigins(ii);
+    BoundingBoxYorigin(jj) = BoundingBoxYorigins(ii);
+end
+
+    
+mymcmd    = sprintf('select ldirX, ldirY, ldirZ, omega, omegaB, difAID, difBID from %s ',acq.pair_tablename);
+[allpairs.ldirX, allpairs.ldirY, allpairs.ldirZ, allpairs.omegaA, allpairs.omegaB, allpairs.difAID, allpairs.difBID] = mym(mymcmd);
+
+for n = first : last
+    fprintf('\nRunning forward simulation for grain %d...\n',n);
+    
     gr = grains.grain{n};     %  output from INDEXTER
 
     graincenter = gr.center / acq.pixelsize;    % graincenter contains COM coordinates in pixels 
@@ -149,22 +186,24 @@ for n = first : last
     bb        = zeros(numspots, 4);                % difspot BoundingBox from database
     status    = cell(numspots, 1);                 % string corresponding to flag (summarizing result of forward simulation)
     spotid    = zeros(numspots, 1);                % Mysql difspottable difspotID
-    check     = zeros(numspots, 5);
+    check     = zeros(numspots, 2);
     hklsp     = zeros(numspots,size(gr.hkl,1));    % signed hk(i)l indices of reflection
     omega     = zeros(numspots, 1);                % rotation angle in degrees
     proj_geom = zeros(numspots, 12);               % Vector used for ROI grain reconstruction in ASTRA (grain shifted to center of roi volume)
     proj_geom_full = zeros(numspots, 12);          % Vector describing full projection geometry (full images, grain at nominal position in sample volume)
     proj_geom_abs  = zeros(numspots, 12);          % Vector describung absorption image projection geometry
 
-    r         = zeros(numspots, 3);                % direction of diffracted beam
+    r         = zeros(numspots, 3);                % direction of diffracted beam used in ASTRA
     r_th      = zeros(numspots, 3);                % theoretically predicted direction of diffracted beam (based on R_vector)
     r_ex      = zeros(numspots, 3);                % estimated direction of diffracted beam, based on grain & difspot COM positions
     U         = zeros(numspots, 1);                % COM u positions of difspots from difspottable
     V         = zeros(numspots, 1);                % COM v positions of difspots from difspottable
     full      = zeros(acq.xdet, acq.ydet);
 
-    %%% Main loop:  examine all hklsp and search corresponding spots / intensity on the detector...
-
+    
+    
+    
+    %%% Main loop:  examine all hklsp and search corresponding spots / intensity on the detector..
     for ii = 1 : numspots
         currentSpotID = ondet(ii);
         pl      = gr.allblobs.pl(currentSpotID, :);
@@ -180,27 +219,19 @@ for n = first : last
         % improved...
         thr_max_offset = round(min(fsim.MaxOmegaOffset, (fsim.omegarange/(abs(sind(eta)))))/omstep);
 
-        % loose search criteria - is there any intensity at the predicted spot
-        % position ?
-        firstquery = sprintf(...
-            ['SELECT difspotID, CentroidX, CentroidY, CentroidImage, ' ...
-                'BoundingBoxXsize, BoundingBoxYsize, BoundingBoxXorigin, ' ...
-                'BoundingBoxYorigin' ...
-             '    FROM %sdifspot' ...
-             '    WHERE (%d between CentroidImage-%d and CentroidImage+%d) '...
-                'and (abs(CentroidX-%d) / BoundingBoxXsize) < 0.5 '...
-                'and (abs(CentroidY-%d) / BoundingBoxYsize) < 0.5 '...
-             '    ORDER BY ((abs(BoundingBoxXsize-%d)/%d)+(abs(BoundingBoxYsize-%d)/%d)) limit 1'],...
-            name, ...
-            image_n, thr_max_offset, thr_max_offset, ...
-            u(currentSpotID), ...
-            v(currentSpotID), ...
-            dif_Xsize, dif_Xsize, dif_Ysize, dif_Ysize);
-
-        [difspotID, CentroidX, CentroidY, CentroidImage, BoundingBoxXsize, ...
-            BoundingBoxYsize, BoundingBoxXorigin, BoundingBoxYorigin] = mym(firstquery);
-
-        if isempty(difspotID)
+        % Look for difspots appying loose search criteria - is there any intensity at the predicted spot
+        % position ?  Choose the one which is closest to the expected size
+        
+        
+        
+        image_ok      = find((CentroidImage >= image_n - thr_max_offset) & (CentroidImage <= image_n + thr_max_offset));
+        intensity_ok  = find(abs(CentroidX(image_ok) - u(currentSpotID)) < max(dif_Xsize/2, BoundingBoxXsize(image_ok)/2 - dif_Xsize/2) & ...
+                abs(CentroidY(image_ok) - v(currentSpotID)) < max(dif_Ysize, BoundingBoxYsize(image_ok)/2 - dif_Ysize/2));
+            
+        candidateIDs   = image_ok(intensity_ok);
+       
+        
+        if isempty(candidateIDs)
             % no candidates have been found: segmentation / overlap / intensity / ... problem ?
 
             flag(ii) = 1; % No intensity at expected spot position
@@ -244,58 +275,52 @@ for n = first : last
                 end % app.assemble_figure
             end % app.check_spot
 
-        else % isempty(difspotID)
-            U(ii)    = CentroidX;
-            V(ii)    = CentroidY;
-            Z(ii)    = CentroidImage;
-            bb(ii, :) = [BoundingBoxXorigin, BoundingBoxYorigin, ...
-                         BoundingBoxXsize, BoundingBoxYsize];
-            spotid(ii) = difspotID;
+        else % isempty(candidateIDs)
+            % in case there are several candidates, take the one which has
+            % the Area close to the expected Area.  This selection "rule" could be
+            % improved...
+
+            [area, area_index]   = sort(abs(BoundingBoxXsize(candidateIDs) .* BoundingBoxYsize(candidateIDs) - dif_Xsize*dif_Ysize));
+            [image, image_index] = sort(abs(CentroidImage(candidateIDs) - image_n));
+            candidateID  = candidateIDs(image_index(1));
+            
+            U(ii)  = CentroidX(candidateID);
+            V(ii)  = CentroidY(candidateID);
+            Z      = CentroidImage(candidateID);
+            bb(ii, :)  = [BoundingBoxXorigin(candidateID), BoundingBoxYorigin(candidateID), ...
+                            BoundingBoxXsize(candidateID), BoundingBoxYsize(candidateID)];
+            spotid(ii) = candidateID;
+            
             flag(ii) = 3; % Found a segmented difspot at expected position - we will check if it does match strict criteria
-            spot = gtGetSummedDifSpot(difspotID, parameters, 1);
+            spot = gtGetSummedDifSpot(candidateID, parameters, 1);
             spot = spot / sum(spot(:)) * dif_Xsize * dif_Xsize * dif_Ysize;   % here we force the intensity to be the same in all spots -
 
             if (fsim.assemble_figure)
                 full = gtPlaceSubImage2(spot, full, bb(ii, 1), bb(ii, 2), 'summed');
             end
-
-            image_ok = between(CentroidImage, image_n - thr_max_offset, image_n + thr_max_offset);
-            centx_ok = abs(CentroidX - u(currentSpotID)) < dif_Offset;
-            centy_ok = abs(CentroidY - v(currentSpotID)) < dif_Offset;
-            xsize_ok = between(BoundingBoxXsize, dif_XsizeMin, dif_XsizeMax);
-            ysize_ok = between(BoundingBoxYsize, dif_YsizeMin, dif_YsizeMax);
-
-            if (fsim.verbose)
-                % useful feedback to find out which of the more strict search criteria fail...
-                sprintf('%d < %d < %d : %d',image_n - thr_max_offset, CentroidImage, image_n + thr_max_offset, image_ok)
-                sprintf('abs(%d - %d) < %d : %d', CentroidX, u(currentSpotID), dif_Offset, centx_ok)
-                sprintf('abs(%d - %d) < %d : %d', CentroidY, v(currentSpotID), dif_Offset, centy_ok)
-                sprintf('%d  < %d < %d : %d', dif_XsizeMin, BoundingBoxXsize, dif_XsizeMax, xsize_ok)
-                sprintf('%d  < %d < %d : %d', dif_YsizeMin, BoundingBoxYsize, dif_YsizeMax, ysize_ok)
-                inputwdefault('continue','yes')
-            end
-
-            check(ii,:) = [image_ok, centx_ok, centy_ok, xsize_ok, ysize_ok];
-
+             
+            cent_ok    = sqrt((U(ii) - u(currentSpotID)).^2 + (V(ii) - v(currentSpotID)).^2) <= dif_Offset;
+            size_ok    = all([between(bb(ii,3), dif_XsizeMin, dif_XsizeMax), between(bb(ii,4), dif_YsizeMin, dif_YsizeMax)]);
+            
+            
+            check(ii,:)  = [cent_ok, size_ok];
+            
+       
             if all(check(ii,:)), flag(ii) = 4; end % otherwise we keep flag == 3
 
-            if find(gr.difspotidA == difspotID)
+            if find(gr.difspotidA == candidateID)
                 % If available, we want to use the angles as determined from the
                 % Friedel pairs - interrogate spotpair table
-                flag(ii) = 5; % This spot is part of a pair and has been used by INDEXTER
-                query    = sprintf('select ldirX, ldirY, ldirZ, omega from %s where difAID=%d', ...
-                                    parameters.acq.pair_tablename, difspotID);
-                [ldirX, ldirY, ldirZ, omegaA] = mym(query);
-                omega(ii) = omegaA;
-                r(ii,:) = [ldirX, ldirY, ldirZ];
+                flag(ii) = 5; % This spot is part of a pair and has been used by INDEXTER              
+                index = find(allpairs.difAID == candidateID);
+                omega(ii) = allpairs.omegaA(index);
+                r(ii,:) = [allpairs.ldirX(index), allpairs.ldirY(index), allpairs.ldirZ(index)];
                 Rottens = gtMathsRotationTensor(-omega(ii), rotcomp);
-            elseif find(gr.difspotidB == difspotID)
+            elseif find(gr.difspotidB == candidateID)
                 flag(ii) = 5; % This spot is part of a pair and has been used by INDEXTER
-                query    = sprintf('select ldirX, ldirY, ldirZ, omegaB from %s where difBID=%d', ...
-                                    parameters.acq.pair_tablename, difspotID);
-                [ldirX, ldirY, ldirZ, omegaB] = mym(query);
-                omega(ii) = omegaB; % in this case we have to reverse the line direction
-                r(ii,:)   = -[ldirX, ldirY, ldirZ];
+                index = find(allpairs.difBID == candidateID);
+                omega(ii) = allpairs.omegaB(index);  % in this case we have to reverse the line direction
+                r(ii,:)   = -[allpairs.ldirX(index), allpairs.ldirY(index), allpairs.ldirZ(index)];
                 Rottens = gtMathsRotationTensor(-omega(ii), rotcomp);
             else
                 % we found a segmentedd spot -  can use the theoretically
@@ -303,7 +328,7 @@ for n = first : last
                 % COM of grain and spot) for backprojection
 
                 % theoretically prediceted omega angle
-                Rottens    = gtMathsRotationTensor(-om_th, rotcomp);
+                Rottens    = gtMathsRotationTensor(-om_th, rotcomp);  
                 % theoretically predicted beam direction
                 r_th(ii,:) = beamdir*Rottens + 2*sind(theta)*pl;
 
@@ -317,6 +342,7 @@ for n = first : last
                 end
             end
 
+            
             hoffset = round((hsize - bb(ii,3)) / 2);
             voffset = round((vsize - bb(ii,4)) / 2);
 
@@ -348,7 +374,7 @@ for n = first : last
             detposFull = detposcorner + detrefu*detdiru0 + detrefv*detdirv0;
             proj_geom_full(ii, :)   = proj_geom(ii,:);
             proj_geom_full(ii, 4:6) = (detposFull * Rottens);
-        end % if isempty(difspotID)
+        end % if isempty(candidateID)
     end % loop over predicted spot positions
 
     found_reflections  = length(find(flag > 1));
@@ -374,8 +400,13 @@ for n = first : last
     %%% Now fill in the information for the output
     %  remove reflections for which no segmented spot has been found
 
-    ind           = find(flag >= 3); % indices of spots to be included in output structure
-
+    if fsim.include_cand
+        ind       = find(flag >= 3); % indices of spots to be included in output structure
+    else
+        ind       = find(flag >= 4);
+    end
+    
+   
     out.selected  = flag(ind) >= 4;  % vector of logicals with pre-selected spots active
     out.id        = gr.id;
     out.phaseid   = phaseID;
@@ -406,7 +437,9 @@ for n = first : last
     out.proj.vol_size_y = hsize;
     out.proj.vol_size_z = vsize;
     out.proj.num_iter   = parameters.rec.num_iter;
-
+    out.vol             = [];
+    out.segbb           = zeros(6,1);
+    out.threshold       = 0;
     out.shift(1) = round((acq.bb(3)-hsize)/2 + graincenter(1)) + 1;
     out.shift(2) = round((acq.bb(3)-hsize)/2 + graincenter(2)) + 1;
     out.shift(3) = round((acq.bb(4)-vsize)/2 + graincenter(3)) + 1;
@@ -441,9 +474,13 @@ function [spot, m] = sfGetRoi(start_image, end_image, acq, bb)
     if all([bb(1) > 0, bb(2) > 0, bb(1)+bb(3) < acq.xdet, bb(2)+bb(4) < acq.ydet])
         fullImgsDir = fullfile(acq.dir, '1_preprocessing', 'full');
         indexes = start_image : end_image;
+         filename = fullfile(fullImgsDir, sprintf('full%04d.edf', indexes(1)));
+         info     = edf_info(filename);
+         
+         
         for img_n = 1:nimages
             fullImgName = fullfile(fullImgsDir, sprintf('full%04d.edf', indexes(img_n)));
-            spot(:, :, img_n) = edf_read(fullImgName, bb);
+            spot(:, :, img_n) = edf_read(fullImgName, bb, [], info);
             m(img_n) = gtImgMeanValue2(spot(centery, centerx, img_n));
         end
     end

5_reconstruction/gtFsimDefaultParameters.m

@@ -5,12 +5,12 @@ fsim.omegarange      = 0.5;       % look in  +/- 1 degrees from predicted image
 fsim.MaxOmegaOffset  = 1;         % Maximum omega offset
 fsim.bb_size         = [];        % default dimensions of ROI, read from the FULL images - if empty the zeropadded projection size will be used
 fsim.assemble_figure = true;      % assemble the results of forward simulation into a full image saved in the grain%d_.mat file
-fsim.display_figure  = true;      %
+fsim.display_figure  = false;     %
 fsim.Grayscale       = [-20 200]; % default colorlimits when visualizing diffraction spots
 fsim.Rdist_factor    = 3;         % allow for multiple times Rdiststd in calulcation of lateral spot offset
 fsim.bbsize_factor   = 2;         % allow for +/- times std difference in horizontal and vertical bb size
 fsim.oversize        = 1.5;       % projections are zeropadded before reconstruction
 fsim.use_th          = true;      % use theoretically predicted diffraction angles for new detected spots
-fsim.verbose         = false;     % 
-
+fsim.verbose         = false;     % tell which search criteria has failed
+fsim.include_cand    = false;     % include candidate spots (flag = 3  :  wrong size) into difspot stack  (deselected)
 end

5_reconstruction/gtSetupReconstruction.m

@@ -49,7 +49,12 @@ disp(['Setting up Reconstruction for phase ' num2str(phaseID) ' (' parameters.cr
 
 samplefile = fullfile(parameters.acq.dir, '4_grains', 'sample.mat');
 if exist(samplefile, 'file')
+    check=inputwdefault('4_grains/sample.mat already exists - do you want to reinitialize it ? [y/n]','n');
+    if strcmpi(check,'y')
+        sample = GtSample();
+    else
        sample = GtSample.loadFromFile(samplefile);
+    end
 else
        sample = GtSample();
 end
@@ -87,7 +92,7 @@ last  = length(grain);
 
 
 % get reconstruction parameters
-check = inputwdefault('Use default parameters? [y/n]', 'y');
+check = inputwdefault('Use default parameters for grain reconstruction ? [y/n]', 'y');
 if ~strcmpi(check, 'y')
     list = build_list_v2();
     parameters.rec = gtModifyStructure(parameters.rec, list.rec(:,1:2));

zUtil_DataStructures/gtCreateGtSample.m

@@ -33,26 +33,38 @@ for phaseid = 1 : phaseNum
 
     phaseDir = sprintf('phase_%02d', phaseid);
     fprintf('Phase %02d: ', phaseid);
+    
+    %% Check if we have already reconstructed volumes or not
+    grainFile = sprintf('grain_%04d.mat', 1);
+    grainPath = fullfile('4_grains', phaseDir, grainFile);
+    grainInfo = load(grainPath);
+    
+    if all(isfield(grainInfo,{'vol', 'segbb', 'threshold'}))
+       variables = {'shift', 'vol', 'segbb', 'selected', 'threshold', 'completeness', 'difspotID'};
+       reconstructed  = 1;
+    else
+       variables = {'shift', 'selected', 'completeness', 'difspotID'};
+       reconstructed = 0;
+    end
+    
     for grainid = 1 : grainsNum
         fprintf('%04d/%04d', grainid, grainsNum);
 
         grainFile = sprintf('grain_%04d.mat', grainid);
         grainPath = fullfile('4_grains', phaseDir, grainFile);
-        grainInfo = load(grainPath, 'shift', 'vol', 'segbb', 'selected', 'threshold', 'completeness', 'difspotID');
+        grainInfo = load(grainPath, variables{:});
 
         phase.setR_vector(grainid, grain{grainid}.R_vector);
         phase.setCenter(grainid, grain{grainid}.center);
-        if (~isfield(grainInfo, 'segbb'))
-            phase.setBoundingBox(grainid, [grainInfo.shift, size(grainInfo.vol)]);
-        else    
+        if reconstructed
             phase.setBoundingBox(grainid, grainInfo.segbb);
+            phase.setThreshold(grainid, grainInfo.threshold);
         end
 
         selected = zeros(length(grainInfo.difspotID),1);
         selected(grainInfo.selected) = 1;
         phase.setSelectedDiffspots(grainid, selected);
-
-        phase.setThreshold(grainid, grainInfo.threshold);
+        
         phase.setCompleteness(grainid, grainInfo.completeness);
         phase.setDifspotID(grainid, grainInfo.difspotID);
         

zUtil_Imaging/gtCheckSegmentation.m

 function gtCheckSegmentation(imnum,parameters,clims)
  
    querie=sprintf('select XIndex,YIndex from %sdifblobdetails where ZIndex=%d',parameters.acq.name,imnum);
-   full=gtGetFullImage(imnum,clims);
+   full=gtGetFullImage(imnum,parameters);
    spots=zeros(parameters.acq.ydet,parameters.acq.xdet);
    [x,y]=mym(querie);
    for i=1:length(x)

zUtil_Parameters/build_list_v2.m

@@ -737,6 +737,9 @@ list.fsim{11,1} = 'use_th';                            list.fsim{11,2} = 'Use th
 list.fsim{12,1} = 'verbose';                           list.fsim{12,2} = 'Display the search criteria for each of the spots';
                                                        list.fsim{12,3} = 'logical';
                                                        list.fsim{12,4} = 2;
+list.fsim{13,1} = 'include_cand';                      list.fsim{13,2} = 'Include difspot candidates into the stack  - automatically disabled';
+                                                       list.fsim{13,3} = 'logical';
+                                                       list.fsim{13,4} = 2;
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % reconstruction