diff --git a/2_difspot/gtDBBlob2SpotTable.m b/2_difspot/gtDBBlob2SpotTable.m index 8d01655a5af5d8575811a2461c02a9a37cb56235..a8e061ce788ae3d6ee401a1fed85544b8c4e164e 100644 --- a/2_difspot/gtDBBlob2SpotTable.m +++ b/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 = []; diff --git a/2_difspot/gtPostSegmentation.m b/2_difspot/gtPostSegmentation.m index bfd99c4c5689285671512472ca7e76cac1fb2e1f..a0397ed27563d6222b9fc39f0d15d9cc51a8137f 100644 --- a/2_difspot/gtPostSegmentation.m +++ b/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 diff --git a/5_reconstruction/gtForwardSimulate_v2.m b/5_reconstruction/gtForwardSimulate_v2.m index 7e6d71394fa69870c69bb38c25d78b66c29a82ad..2b37fbcb3d2908d33d6e64bfd942cf343e630ebc 100644 --- a/5_reconstruction/gtForwardSimulate_v2.m +++ b/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 diff --git a/5_reconstruction/gtFsimDefaultParameters.m b/5_reconstruction/gtFsimDefaultParameters.m index 996a74ab0023a8ac966503726341a8b19aa547c1..4b282d35888e1bc1a75d975768a117e926fdd39d 100644 --- a/5_reconstruction/gtFsimDefaultParameters.m +++ b/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 diff --git a/5_reconstruction/gtSetupReconstruction.m b/5_reconstruction/gtSetupReconstruction.m index bc3c34d7870ffebfb488a23adbc98ce4cc05f21a..b590be24b4e7e033f1df18d5ff3a5a221bafc203 100644 --- a/5_reconstruction/gtSetupReconstruction.m +++ b/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)); diff --git a/zUtil_DataStructures/gtCreateGtSample.m b/zUtil_DataStructures/gtCreateGtSample.m index d1acd8088ca54160e2181bd26e2eb494844f4184..962f6f5bca1d815172b2f3e73e3aba0e8a23cd32 100644 --- a/zUtil_DataStructures/gtCreateGtSample.m +++ b/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); diff --git a/zUtil_Imaging/gtCheckSegmentation.m b/zUtil_Imaging/gtCheckSegmentation.m index a3c79dce92cc9182fee99a5cd5a3be20f12038a5..4af4cf39ed5c8285e881863efbc4d6a84f013c85 100644 --- a/zUtil_Imaging/gtCheckSegmentation.m +++ b/zUtil_Imaging/gtCheckSegmentation.m @@ -1,7 +1,7 @@ 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) diff --git a/zUtil_Parameters/build_list_v2.m b/zUtil_Parameters/build_list_v2.m index a88bdbf2ed27aa8db602e5e1fd8853b67763b007..e0595a21515f288fc36c75c49277a8100c6fdde6 100644 --- a/zUtil_Parameters/build_list_v2.m +++ b/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