Skip to content
Snippets Groups Projects
Commit 2f3ce063 authored by Peter Reischig's avatar Peter Reischig Committed by Nicola Vigano
Browse files

Cleaning. Deleted old functions in folder 4_grains and adjusted a few accordingly. 

git-svn-id: https://svn.code.sf.net/p/dct/code/trunk@510 4c865b51-4357-4376-afb4-474e03ccb993
parent 75a2c583
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 0 additions and 1773 deletions
4_grains/gt2LinesAngle.m deleted 100644 → 0
+ 0
− 13
View file @ 75a2c583
%
% Calculates the angle (in degrees) of two straight lines (l1 and l2) in 3D.
%
% INPUT: Two straight lines l1 and l2 = [ point_X point_Y point_Z dir_X dir_Y dir_Z ]
%
% OUTPUT: Angle (in degrees) between l1 and l2 = out
%
function out = gt2LinesAngle(l1,l2)
out = acosd(dot(l2(4:6),l1(4:6))/norm(l2(4:6))/norm(l1(4:6))) ;
end
4_grains/gt2LinesDist.m deleted 100644 → 0
+ 0
− 27
View file @ 75a2c583
% Gives the spatial distance between a given base line and other 3D lines pair-wise
%
% INPUT: bl= base line [ point_X point_Y point_Z dir_X dir_Y dir_Z ]
% ol= other lines [ point_X point_Y point_Z dir_X dir_Y dir_Z ]
%
% OUTPUT: Spatial distance between bl and ol-s
%
function out = gt2LinesDist(bl,ol)
nol=size(ol,1);
blr=repmat(bl,nol,1);
a=cross(blr(:,4:6),ol(:,4:6),2);
out = abs(dot(ol(:,1:3)-blr(:,1:3),a,2)./sqrt(sum(a.*a,2))) ;
end
% function out = gt2LinesDist(l1,l2)
%
% a=cross(l1(4:6),l2(4:6)) ;
%
% out = abs(dot(l2(1:3)-l1(1:3),a)/norm(a)) ;
%
% end
4_grains/gt2LinesIntersection.m deleted 100644 → 0
+ 0
− 50
View file @ 75a2c583
%
% Gives the closest point between a given base line and other 3D lines pair-wise
% (their 'point of intersection').
%
% INPUT: bl[point_X point_Y point_Z norm_dir_X norm_dir_Y norm_dir_Z] =
% equation of the base line, direction vector normalized
% ol[point_X point_Y point_Z norm_dir_X norm_dir_Y norm_dir_Z] =
% equation of the ohter 3D lines, direction vector normalized
%
% OUTPUT: out[ X;Y;Z ] = closest point to pairs (`intersection point`)
%
function out = gt2LinesIntersection(bl,ol)
nol=size(ol,1);
A11=bl(4:6)*bl(4:6)';
A12=-ol(:,4:6)*bl(4:6)';
A21=A12;
A22=sum(ol(:,4:6).*ol(:,4:6),2);
c0=ol(:,1:3)-repmat(bl(1:3),nol,1);
c11=c0*bl(4:6)';
c21=sum(-c0.*ol(:,4:6),2);
out=zeros(nol,3);
for i=1:nol
if bl==ol(i,:)
out(i,1:3)=bl(1:3);
continue
end
A=[A11 A12(i); A21(i) A22(i)];
c(1,1)=c11(i);
c(2,1)=c21(i);
par=A\c;
out1=(bl(1:3)+par(1)*bl(4:6))' ;
out2=(ol(i,1:3)+par(2)*ol(i,4:6))' ;
out(i,1:3) = (out1+out2)/2 ;
end
end
4_grains/gtAngConsCosMat.m deleted 100644 → 0
+ 0
− 12
View file @ 75a2c583
function ACM=gtAngConsCosMat(spacegroup)
no_ref=size(parameters.cryst.hklsp,1);
for i=1:no_ref
for j=i:no_ref
ACM{i,j}=gtAngularConsCos(spacegroup,i,j);
ACM{j,i}=ACM{i,j} ;
end
end
end
4_grains/gtAngConsCosVec.m deleted 100644 → 0
+ 0
− 14
View file @ 75a2c583
function ACM_vec=gtAngConsCosVec(spacegroup)
ACM=gtAngConsMat(spacegroup);
ACM_vec=[];
for i=1:size(ACM,1)
for j=i:size(ACM,2)
ACM_vec=[ACM_vec; ACM{i,j}] ;
end
end
end
4_grains/gtAngConsMat.m deleted 100644 → 0
+ 0
− 17
View file @ 75a2c583
function [ACM,ACMU]=gtAngConsMat(cryst)
no_ref=size(cryst.hklsp,1);
%ACM=zeros(no_ref,no_ref) ;
for i=1:no_ref
for j=i:no_ref
[uniangle,angle_mult]=gtAngularCons(cryst,i,j);
ACM{i,j}=uniangle;
ACM{j,i}=uniangle;
ACMU{i,j}=angle_mult; % multiplicity
ACMU{j,i}=angle_mult;
end
end
end
4_grains/gtAngularCons.m deleted 100644 → 0
+ 0
− 43
View file @ 75a2c583
function [uniangle,angle_mult]=gtAngularCons(cryst,f1,f2)
ref=cryst.hklsp;
spacegroup=cryst.spacegroup;
angle=[];
%%% if cubic (221, 225, 229) %%%
if spacegroup==225 || spacegroup==229 || spacegroup==221 || spacegroup==227
v1=gtGetReflections(ref(f1,:)) ;
v2=gtGetReflections(ref(f2,:)) ;
%%% if hexagonal (194, 663) %%%
elseif spacegroup==194 || spacegroup==663 || spacegroup==167 || spacegroup==152
v1=gtGetReflections_sab(ref(f1,:)) ;
v2=gtGetReflections_sab(ref(f2,:)) ;
v1=gtHKL2Cart(v1);
v2=gtHKL2Cart(v2);
%%% otherwise error %%%
else
error('Currently only spacegroups 225, 229, 194, 663 (167?) supported')
end
for i=1:size(v1,1)
for j=1:size(v2,1)
if v1(i,:)==v2(j,:)
angle=[angle; 0];
elseif v1(i,:)==-v2(j,:)
angle=[angle; 0];
else
angle=[angle; acosd(abs((v1(i,:)*v2(j,:)')/norm(v1(i,:))/norm(v2(j,:))))] ;
end
end
end
uniangle=unique(angle);
for i=1:length(uniangle)
angle_mult(i)=length(find(uniangle(i)==angle));
end
end
4_grains/gtAngularConsCos.m deleted 100644 → 0
+ 0
− 17
View file @ 75a2c583
function cos=gtAngularConsCos(spacegroup,f1,f2)
ref=parameters.cryst.hklsp;
cos=[];
v1=gtGetReflections(ref(f1,:)) ;
v2=gtGetReflections(ref(f2,:)) ;
for i=1:size(v1,1)
for j=1:size(v2,1)
cos=[cos; abs(v1(i,:)*v2(j,:)')/norm(v1)/norm(v2)] ;
end
end
cos=unique(cos);
cos=[1; cos];
end
4_grains/gtCheckGroupCons.m deleted 100644 → 0
+ 0
− 16
View file @ 75a2c583
function groupcons=gtCheckGroupCons(tryline,goodl,goodpofinters,tol,sample,ACM)
nof_goods=length(goodl.id);
paircons=gtCheckPairCons(tryline,goodl,tol,sample,ACM);
groupcons=false;
if length(paircons)==nof_goods
if pointtolinedist(goodpofinters,[tryline.ca tryline.dir])<2*tol.dist;
groupcons=true;
end
end
end
4_grains/gtCheckPairCons.m deleted 100644 → 0
+ 0
− 67
View file @ 75a2c583
%
% Checks pair consistency in a set of lines.
%
% act: the actual line (structure)
% cand: candidates (structure) to be pair consistent with it
% tol: tolerances (structure)
% sample: sample geometry parameters (structure)
% ACM: angular consistency matrix
function [paircons,progress,cand,pofinters,angle]=gtCheckPairCons(act,cand,tol,sample,ACM)
paircons=(1:length(cand.id))'; % initial indices of consistent lines
progress=length(cand.id);
% integrated intensities close enough?
cons=(act.int/tol.int<cand.int) & (act.int*tol.int>cand.int);
cand=gtKeepLine(cons,cand);
paircons(~cons)=[];
progress=[progress; length(cand.id)];
% bbox x size close enough?
cons=(act.bbxs/tol.bbxs<cand.bbxs) & (act.bbxs*tol.bbxs>cand.bbxs);
cand=gtKeepLine(cons,cand);
paircons(~cons)=[];
progress=[progress; length(cand.id)];
% bbox y size close enough?
cons=(act.bbys/tol.bbys<cand.bbys) & (act.bbys*tol.bbys>cand.bbys);
cand=gtKeepLine(cons,cand);
paircons(~cons)=[];
progress=[progress; length(cand.id)];
% spatial distance small enough?
cons=gt2LinesDist([act.ca act.dir],[cand.ca cand.dir])<tol.dist;
cand=gtKeepLine(cons,cand);
paircons(~cons)=[];
progress=[progress; length(cand.id)];
% is intersection in the sample volume?%
% NOTE: intersection of 2 lines might actually be far from the grain center
% if they are near being parallel. This could be taken into account here...
% On the other hand, there should be no parallel lines in a grain set.
pofinters=gt2LinesIntersection([act.ca act.dir],[cand.ca cand.dir]);
cons=gtIsPointInSample(pofinters,sample,tol);
cand=gtKeepLine(cons,cand);
pofinters(~cons,:)=[];
paircons(~cons)=[];
progress=[progress; length(cand.id)];
% plane angles consistent?
cons=false(length(cand.id),1);
angle=zeros(length(cand.id),1);
for i=1:length(cand.id)
ACV=ACM{act.thetatype,cand.thetatype(i)};
ang=gt2PlanesAngle(act.pl,cand.pl(i,:));
angle(i)=ang;
cons(i)=min(abs(ACV-ang))<tol.ang;
end
cand=gtKeepLine(cons,cand);
pofinters(~cons,:)=[];
angle(~cons)=[];
paircons(~cons)=[];
progress=[progress; length(cand.id)];
end
4_grains/gtDiffPlaneNormInLab.m deleted 100644 → 0
+ 0
− 42
View file @ 75a2c583
%
% FUNCTION plnorm=gtDiffPlaneNormInLab(scXA,scYA,scXB,scYB,rottodet,tiltY,tiltZ,rotx,sorZ)
%
% Given the Scintillator coordinates of pairs of diffraction spots,
% computes the diffracting plane normals they correspond to, according to
% tilts and rotation axis to detector (scintillator) distance.
%
% LAB system is right-handed and defined as:
% origin X = on the rotation axis, positive towards the camera
% origin Y = on the rotation axis
% origin Z = floor(parameters.acq.ydet/2) in the middle of the central
% pixel row of the image, positive upwards
%
% INPUT scXA = coloumn vector of X coordinates on the Scintillator of point A (pixels or mm-s)
% scYA = coloumn vector of Y coordinates on the Scintillator of point A (pixels or mm-s)
% scXB = coloumn vector of X coordinates on the Scintillator of point B (pixels or mm-s)
% scYB = coloumn vector of Y coordinates on the Scintillator of point B (pixels or mm-s)
% rottodet = rotation axis to detector (scintillator) distance at rotx,sorZ (pixels or mm-s)
% /scX,scY,rottodet have to be in the same units; e.g. pixels or mm-s/
% tiltY = tilt around axis Y in Lab system (in degrees)
% tiltZ = tilt around axis Z in Lab system (in degrees)
% rotx = location of rotation axis in the images (X coord. in pixels)
% sorZ = Z=0 location in the images (Y coord. in pixels at the rotation axis)
%
% OUTPUT plnorm = coloumn vector of plane normals of size(plnorm)=[n,3]
%
%
function plnorm=gtDiffPlaneNormInLab(scXA,scYA,scXB,scYB,rottodet,tiltY,tiltZ,rotx,sorZ)
dv= gtDiffVecInLab(scXA,scYA,scXB,scYB,rottodet,tiltY,tiltZ,rotx,sorZ); % diffraction vector
plv= dv-repmat([1 0 0],size(dv,1),1); % plane vector
plvnorm=sqrt(plv(:,1).^2+plv(:,2).^2+plv(:,3).^2) ; % norm of plane vectors
plnorm(:,1)= plv(:,1)./plvnorm ; % normalized plane vectors
plnorm(:,2)= plv(:,2)./plvnorm ;
plnorm(:,3)= plv(:,3)./plvnorm ;
end
4_grains/gtDiffVecInLab.m deleted 100644 → 0
+ 0
− 45
View file @ 75a2c583
%
% function Dv=gtDiffVecInLab(scXA,scYA,scXB,scYB,rottodet,tiltY,tiltZ,rotx,sorZ)
%
% Given the Scintillator coordinates of a pair of diffraction spots,
% computes the real theta angle they correspond to, according to tilts and
% rotation axis to detector (scintillator) distance.
%
% LAB system is right-handed and defined as:
% origin X = on the rotation axis, positive towards the camera
% origin Y = on the rotation axis
% origin Z = floor(parameters.acq.ydet/2) in the middle of the central
% pixel row of the image, positive upwards
%
% INPUT scXA = coloumn vector of X coordinates on the Scintillator of point A (pixels or mm-s)
% scYA = coloumn vector of Y coordinates on the Scintillator of point A (pixels or mm-s)
% scXB = coloumn vector of X coordinates on the Scintillator of point B (pixels or mm-s)
% scYB = coloumn vector of Y coordinates on the Scintillator of point B (pixels or mm-s)
% rottodet = rotation axis to detector (scintillator) distance at rotx,sorZ (pixels or mm-s)
% /scX,scY,rottodet have to be in the same units; e.g. pixels or mm-s/
% tiltY = tilt around axis Y in Lab system (in degrees)
% tiltZ = tilt around axis Z in Lab system (in degrees)
% rotx = location of rotation axis in the images (X coord. in pixels)
% sorZ = Z=0 location in the images (Y coord. in pixels at the rotation axis)
%
% e.g. ([234;334;434],[1444;1644;1944],[134;234;534],[555;545;535],
% 3600,0.123,0.421,1026,1024)
%
% OUTPUT theta = coloumn vector of real theta angles
%
% e.g. [7.3948; 7.1822; 6.9281]
%
%
function Dv=gtDiffVecInLab(scXA,scYA,scXB,scYB,rottodet,tiltY,tiltZ,rotx,sorZ)
Dv(:,1)=2*rottodet-sind(tiltZ).*(scXA+scXB-2*rotx)-sind(tiltY).*cosd(tiltZ).*(scYA+scYB-2*sorZ) ;
Dv(:,2)=cosd(tiltZ).*(scXA+scXB-2*rotx)-sind(tiltY).*sind(tiltZ).*(scYA+scYB-2*sorZ) ;
Dv(:,3)=-cosd(tiltY).*(scYA-scYB) ;
Dnorm=sqrt(Dv(:,1).^2+Dv(:,2).^2+Dv(:,3).^2) ;
Dv(:,1)=Dv(:,1)./Dnorm ;
Dv(:,2)=Dv(:,2)./Dnorm ;
Dv(:,3)=Dv(:,3)./Dnorm ;
end
4_grains/gtEtaOfPairs.m deleted 100644 → 0
+ 0
− 71
View file @ 75a2c583
%
% function etas=gtEtaOfPairs(scXA,scYA,scXB,scYB,tiltY,tiltZ,rotx,sorZ)
%
% Given the Scintillator coordinates of a pair of diffraction spots,
% computes the real eta angle they correspond to, according to tilts.
%
% LAB system is right-handed and defined as:
% origin X = on the rotation axis, positive towards the camera
% origin Y = on the rotation axis
% origin Z = floor(parameters.acq.ydet/2) in the middle of the central
% pixel row of the image, positive upwards
%
% INPUT scXA = coloumn vector of X coordinates on the Scintillator of point A (pixels or mm-s)
% scYA = coloumn vector of Y coordinates on the Scintillator of point A (pixels or mm-s)
% scXB = coloumn vector of X coordinates on the Scintillator of point B (pixels or mm-s)
% scYB = coloumn vector of Y coordinates on the Scintillator of point B (pixels or mm-s)
% /scX,scY have to be in the same units; e.g. pixels or mm-s/
% tiltY = tilt around axis Y in Lab system (in degrees)
% tiltZ = tilt around axis Z in Lab system (in degrees)
% rotx = location of rotation axis in the images (X coord. in pixels)
% sorZ = Z=0 location in the images (Y coord. in pixels at the rotation axis)
%
% e.g. ([234;334;434],[1444;1644;1944],[134;234;534],[555;545;535],
% 0.123,0.421,1026,1024)
%
% OUTPUT etas = coloumn vector of real eta angles
%
% e.g. [22.5448; 165.1822; 345.9281]
%
%
function etas=gtEtaOfPairs(scXA,scYA,scXB,scYB,tiltY,tiltZ,rotx,sorZ)
% Diffraction vector coordinates
%Dv1=2*rottodet-sind(tiltZ).*(scXA+scXB-2*rotx)-sind(tiltY).*cosd(tiltZ).*(scYA+scYB-2*sorZ) ;
Dv2=cosd(tiltZ).*(scXA+scXB-2*rotx)-sind(tiltY).*sind(tiltZ).*(scYA+scYB-2*sorZ) ;
Dv3=-cosd(tiltY).*(scYA-scYB) ;
etas=gtEtaOfPoint(Dv2,Dv3);
% etas=zeros(length(Dv2),1);
%
% for i=1:length(Dv2)
% if Dv2(i) > 0 % 0 < eta < 180deg
% if Dv3(i) > 0 % 0 < eta < 90deg
% etas(i)=atand(Dv2(i)/Dv3(i)) ;
% elseif Dv3(i) == 0
% etas(i)=90;
% else % 90deg < eta < 180deg
% etas(i)=atand(Dv2(i)/Dv3(i))+180 ;
% end
% elseif Dv2(i) == 0
% if Dv3(i) > 0
% etas(i)=0 ;
% elseif Dv3(i) == 0
% etas(i)=0;
% else
% etas(i)=180 ;
% end
% else % 180deg < eta < 360deg
% if Dv3(i) < 0 % 180deg < eta < 270deg
% etas(i)=atand(Dv2(i)/Dv3(i))+180 ;
% elseif Dv3(i) == 0
% etas(i)=270;
% else % 270deg < eta < 360deg
% etas(i)=atand(Dv2(i)/Dv3(i))+360 ;
% end
% end
% end
end
4_grains/gtEtaOfPoint.m deleted 100644 → 0
+ 0
− 41
View file @ 75a2c583
% Eta angle from lab coordinates (e.g. diff vectors).
%
% Lab coordinates (right-handed):
% x: along the beam
% y: horizontal (+x in images)
% z: vertical upwards (-y in images)
function etas=gtEtaOfPoint(y,z)
etas=zeros(length(y),1);
for i=1:length(y)
if y(i) > 0 % 0 < eta < 180deg
if z(i) > 0 % 0 < eta < 90deg
etas(i)=atand(y(i)/z(i)) ;
elseif z(i) == 0
etas(i)=90;
else % 90deg < eta < 180deg
etas(i)=atand(y(i)/z(i))+180 ;
end
elseif y(i) == 0
if z(i) > 0
etas(i)=0 ;
elseif z(i) == 0
etas(i)=0;
else
etas(i)=180 ;
end
else % 180deg < eta < 360deg
if z(i) < 0 % 180deg < eta < 270deg
etas(i)=atand(y(i)/z(i))+180 ;
elseif z(i) == 0
etas(i)=270;
else % 270deg < eta < 360deg
etas(i)=atand(y(i)/z(i))+360 ;
end
end
end
4_grains/gtEvaluateDifspotMatch.m deleted 100644 → 0
+ 0
− 335
View file @ 75a2c583
%
% FUNCTION gtEvaluateDifspotMatch(whichID,varargin)
%
% Does not modify any data, only gives feedback on how good
% the matching of diffraction spots is by showing pairs from database
% according to the given difspotID-s from the first set.
% Run from the main folder of the scan (it uses the parameters file).
%
% USAGE e.g.
% gtEvaluateDifspotMatch([],'pair_dist')
% gtEvaluateDifspotMatch([200:300],'pairfigure','pairdata','dopause','button')
% gtEvaluateDifspotMatch([1000:4000],'table','tol_hist')
%
% INPUT
% whichID: vector of difspotID-s from the first set to be shown; if
% empty, all are shown
% 'pairfigure': pairs found are shown graphically
% 'pairdata': table of spot properties for comparison
% 'table': table for comparison (feedback to check tolerances)
% 'tol_hist': histogram of properties of pair members (to check tolerances)
% 'pair_dist': pair distribution vs omega, eta, and eta vs theta
% 'dopause': pause in sec or 'button' after each pair. Use it in pair:
% 'dopause','button' or e.g. 'dopause',2 .
%
% OUTPUT figures and tables for comparison
%
%
% Peter Reischig, ESRF, 01/2007
%
function gtEvaluateDifspotMatch(whichID,varargin)
load('parameters.mat');
%% Set parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
par.dopause=[];
par.pairfigure=false;
par.pairdata=false;
par.table=false;
par.tol_hist=false;
par.pair_dist=false;
par=sfParseFunctionInput(par,varargin);
% Thresholds for image position
thr_ang=parameters.match.thr_ang; % angular threshold; absolute value in degrees
thr_max_offset=parameters.match.thr_max_offset; % MaxImage offset threshold; absolut value in #images
thr_ext_offset=parameters.match.thr_ext_offset; % ExtStartImage and ExtEndImage offset threshold; absolut value in #images
thr_genim_offset=parameters.match.thr_genim_offset; % at least one of the three images should be offset as maximum this value in #images
corr_rot_to_det=parameters.match.corr_rot_to_det; % correction of rot_to_det in pixels
% Thresholds ( for search within the limits: basevalue/thr & & basevalue*thr )
thr_intint=parameters.match.thr_intint; % integrated intensity
thr_area=parameters.match.thr_area; % area of spots inside the bounding box
thr_bbsize=parameters.match.thr_bbsize; % bounding box size
sample_radius=parameters.acq.bb(3)/2;
sample_top=parameters.acq.bb(2);
sample_bot=parameters.acq.bb(2)+parameters.acq.bb(4);
rot_to_det=parameters.acq.dist/parameters.acq.pixelsize + corr_rot_to_det;
% Set sample coordinate system
% sorX = on the rotation axis by definition
% sorY = on the rotation axis by definition
%sorZ=sample_top; % at the top of the sample bounding box (acquisition); lowest number in pixels
%fullimages_directory=sprintf('%s/1_preprocessing/full',parameters.acq.dir) ;
pairtable=parameters.acq.pair_tablename ;
if isfield(parameters.acq, 'difA_name')
difspottable=[parameters.acq.difA_name 'difspot'];
else
difspottable=[parameters.acq.name 'difspot'];
end
nproj=parameters.acq.nproj; % number of preojections per 180 degrees
gtDBConnect
%close all
comp_2thetas=parameters.cryst.twoth; % coloumn vector of all valid 2theta angles in degrees
if ~exist('whichID','var')
whichID=[];
end
if isempty(whichID)
whichID=mym(sprintf('Select difspotID from %s',difspottable));
end
cmp=[];
%% Main loop for spotsA
if par.pairdata || par.pairfigure || par.tol_hist || par.table
for i=1:length(whichID)
[sp1.id,sp1.maxim,sp1.extstim,sp1.extendim,sp1.area,sp1.centX,sp1.centY,sp1.bbX,sp1.bbY,sp1.bbXsize,sp1.bbYsize,sp1.integral]=...
mym(sprintf(['select difspotID,CentroidImage,ExtStartImage,ExtEndImage,Area,CentroidX,CentroidY,BoundingBoxXorigin,BoundingBoxYorigin,'...
'BoundingBoxXsize,BoundingBoxYsize,Integral from %s where difspotID=%d'],difspottable,whichID(i))) ;
[pairID,sp1.pair,theta,eta,omega]=mym(sprintf('select pairID,difBID,theta,eta,omega from %s where difAID=%d', pairtable, whichID(i)));
if isempty(sp1.pair)
if par.pairdata
disp(sprintf('Difspot of ID %d is not matched',whichID(i)))
end
continue
end
[sp2.id,sp2.maxim,sp2.extstim,sp2.extendim,sp2.area,sp2.centX,sp2.centY,sp2.bbX,sp2.bbY,sp2.bbXsize,sp2.bbYsize,sp2.integral]=...
mym(sprintf(['select difspotID,CentroidImage,ExtStartImage,ExtEndImage,Area,CentroidX,CentroidY,BoundingBoxXorigin,BoundingBoxYorigin,'...
'BoundingBoxXsize,BoundingBoxYsize,Integral from %s where difspotID=%d'],difspottable,sp1.pair)) ;
add_cmp.pairID=pairID;
add_cmp.id1=sp1.id;
add_cmp.id2=sp2.id;
add_cmp.maxim=sp2.maxim-sp1.maxim;
add_cmp.extstim=sp2.extstim-sp1.extstim;
add_cmp.extendim=sp2.extendim-sp1.extendim;
add_cmp.integral=sp2.integral/sp1.integral*100;
add_cmp.area=sp2.area/sp1.area*100;
add_cmp.bbXsize=sp2.bbXsize/sp1.bbXsize*100;
add_cmp.bbYsize=sp2.bbYsize/sp1.bbYsize*100;
cmp=[cmp, add_cmp];
pr_bb = projected_sample(sp1.centX,sp1.centY,sample_radius,sample_top,sample_bot,rot_to_det,parameters.acq.rotx) ;
if par.pairdata
disp(' ')
disp(' difspotID CentroidIm ExtStartIm ExtEndIm Integral Area BBoxXsize BBoxYsize')
disp('+----------+ +----------+ +----------+ +----------+ +----------+ +----------+ +----------+ +----------+')
disp(sprintf(' %10d %10.3f %10d %10d %10d %10d %10d %10d', ...
sp1.id,sp1.maxim,sp1.extstim,sp1.extendim,sp1.integral,sp1.area,sp1.bbXsize,sp1.bbYsize)) ;
disp(sprintf(' %10d %10.3f %10d %10d %10d %10d %10d %10d', ...
sp2.id,sp2.maxim,sp2.extstim,sp2.extendim,sp2.integral,sp2.area,sp2.bbXsize,sp2.bbYsize)) ;
disp(sprintf(' Pair %5d %8.3fim %8dim %8dim %9.2f%% %9.2f%% %9.2f%% %9.2f%% ', ...
pairID,add_cmp.maxim,add_cmp.extstim,add_cmp.extendim,add_cmp.integral,add_cmp.area,add_cmp.bbXsize,add_cmp.bbYsize));
end
if par.pairfigure
gtShowPairFigure(1,parameters,sp1,sp2,pr_bb,comp_2thetas,0,[],[],[],omega,theta,eta)
end
if strcmpi('button',par.dopause)
waitforbuttonpress
else
pause(par.dopause)
end
end % end of main loop for spots1
end
if par.table
disp(' ')
disp(' TABLE FOR COMPARISON ')
disp(' ')
disp(' difspotIDA difspotIDB Pair ID CentroidIm ExtStartIm ExtEndIm Integral Area BBoxXsize BBoxYsize')
disp('+----------+ +----------+ +----------+ +----------+ +----------+ +----------+ +----------+ +----------+ +----------+ +----------+')
for i=1:length(cmp)
if ~isempty(cmp(i).pairID)
disp(sprintf(' %10d %10d %10d %8.3fim %8dim %8dim %9.2f%% %9.2f%% %9.2f%% %9.2f%% ', ...
cmp(i).id1,cmp(i).id2,cmp(i).pairID,cmp(i).maxim,cmp(i).extstim,cmp(i).extendim,cmp(i).integral,cmp(i).area,cmp(i).bbXsize,cmp(i).bbYsize));
end
end
end
if par.tol_hist
for i=1:length(cmp)
hist_maxim(i)=cmp(i).maxim;
hist_extstim(i)=cmp(i).extstim;
hist_extendim(i)=cmp(i).extendim;
hist_integral(i)=cmp(i).integral;
hist_area(i)=cmp(i).area;
hist_bbXsize(i)=cmp(i).bbXsize;
hist_bbYsize(i)=cmp(i).bbYsize;
end
bins=max([10 ceil(length(cmp)/20)]);
centbins=nproj-ceil(thr_max_offset):0.1:nproj+ceil(thr_max_offset);
extbins=(nproj-ceil(thr_ext_offset)):(nproj+ceil(thr_ext_offset));
figure('name','Difference in Centroid Images B-A')
hist(hist_maxim,centbins)
figure('name','Difference in Extinction Start Images; B-A')
hist(hist_extstim,extbins)
figure('name','Difference in Extinction End Images; B-A')
hist(hist_extendim,extbins)
figure('name','Ratio of Integrated Intensities; B/A %')
hist(hist_integral,bins)
figure('name','Ratio of Areas; B/A %')
hist(hist_area,bins)
figure('name','Ratio of BBox X size; B/A %')
hist(hist_bbXsize,bins)
figure('name','Ratio of BBox Y size; B/A %')
hist(hist_bbYsize,bins)
end
if par.pair_dist
% Matching percentage
binwidth=5;
[centimA]=mym(sprintf('SELECT CentroidImage from %s where CentroidImage<=%d',difspottable,nproj))*180/nproj;
[centimB]=mym(sprintf('SELECT CentroidImage from %s where CentroidImage>%d',difspottable,nproj))*180/nproj;
edgesA=0:binwidth:180;
edgesB=180:binwidth:360;
centimA=histc(centimA,edgesA);
centimB=histc(centimB,edgesB);
centom=min(centimA,centimB);
[theta,eta,omega]=mym(sprintf('SELECT theta,eta,omega from %s',pairtable));
match_omega=histc(omega,edgesA);
edges_eta=0:binwidth:360;
match_eta=histc(eta,edges_eta);
matchpc=match_omega./centom*100;
figure('name','Eta vs omega')
plot(omega,eta,'b.')
axis([0 180 0 380])
figure('name','Theta vs omega')
plot(omega,theta,'b.')
hold on
for i=1:length(comp_2thetas)
plot([0 380],[comp_2thetas(i) comp_2thetas(i)]/2,'g-.')
end
axis([0 180 min(comp_2thetas)/2-1 max(comp_2thetas)/2+1])
figure('name','Eta vs theta')
plot(theta,eta,'b.')
hold on
for i=1:length(comp_2thetas)
plot([comp_2thetas(i) comp_2thetas(i)]/2,[0 380],'g-.')
end
axis([min(comp_2thetas)/2-1 max(comp_2thetas)/2+1 0 380])
figure('name','No. of pairs vs theta')
edges_theta=round((min(comp_2thetas)/2-0.1:180/nproj/5:max(comp_2thetas)/2+0.1)*100)/100;
theta_h=histc(theta,edges_theta);
bar(edges_theta,theta_h)
h=findobj(gca,'Type','patch');
set(h,'EdgeColor','none')
hold on
for i=1:length(comp_2thetas)
plot([comp_2thetas(i) comp_2thetas(i)]/2,[0 max(theta_h)],'g-.')
end
xlim([min(edges_theta) max(edges_theta)])
figure('name','No. of pairs found vs eta')
bar((edges_eta+binwidth/2)',match_eta)
h=findobj(gca,'Type','patch');
set(h,'EdgeColor','none')
xlim([0,360])
figure('name','Matching percentage vs omega')
bar((edgesA+binwidth/2)',matchpc)
xlim([0,180])
h=findobj(gca,'Type','patch');
set(h,'EdgeColor','none')
set(h,'FaceColor','r')
figure('name','No. of spots and no. of pairs vs omega')
bar((edgesA+binwidth/2)',[centom,match_omega])
h=findobj(gca,'Type','patch');
set(h,'EdgeColor','none')
xlim([0,180])
end
end % of function
%% Parse input arguments
function par=sfParseFunctionInput(par,addpar)
if length(addpar)<=0 % just return the defaults
return
end
if ~isstruct(par)
error 'No structure for defaults was supplied'
end
propnames = fieldnames(par);
lpropnames = lower(propnames);
for i=1:length(lpropnames)
if islogical(par.(propnames{i}))
j=1;
while j<=length(addpar)
if strmatch(lpropnames{i},lower(addpar{j}),'exact')
par.(propnames{i})=~par.(propnames{i});
addpar(j)=[];
else
j=j+1;
end
end
end
end
n=length(addpar)/2;
if n~=floor(n)
error 'Property/value pairs must come in PAIRS.'
end
for i=1:n
p_i=lower(addpar{2*i-1});
v_i=addpar{2*i};
ind=strmatch(p_i,lpropnames,'exact'); % index of matching property
if isempty(ind)
ind=find(strncmp(p_i,lpropnames,length(p_i)));
if isempty(ind)
error(['No matching property found for: ',addpar{2*i-1}])
elseif length(ind)>1
error(['Ambiguous property name: ',addpar{2*i-1}])
end
end
par.(propnames{ind})=v_i;
end
end % of sfParseFunctionInput
4_grains/gtMatchDifspots.m deleted 100644 → 0
+ 0
− 474
View file @ 75a2c583
%
% FUNCTION gtMatchDifspots(par_match,savepar,show_findings,dopause)
%
% (Used to be gtMatchDifspots_v3)
% Matches difspots of a 360 degree scan, from a single difspot table.
% It can take into account tilts and correlate difspots (check if it
% gives any improvement) for better angular precision.
% Speed shouldn't depend much on the number of difspots, since it works
% on a fraction of the spots (creates and uses a smaller spotpairtable,
% named 'processtable').
%
% USAGE e.g. gtMatchDifspots
% gtMatchDifspots(par_match)
% gtMatchDifspots(par_match,0,0,'button')
%
% OPTIONAL INPUT
%
% par_match: matching parameters; if not specified it looks into
% the parameters file
% savepar: saves the last used match parameters in the parameters file
% show_findings: if true, the pairs found and candidates are shown in figures
% dopause: pause time in sec after a pair is found; if set as 'button',
% it waits for button press
%
% OUTPUT loads data in spotpairs table given in the parameters file
%
%
% MATCHING TOLERANCES read from parameters.match (default values in brackets)
%
% centmode: centroid of difspots based on:
% 0 or []: uses center of mass from difspot table (default)
% 1: difspots from blobs or edf-s are correlated
% 2: difspots - from summed fulls - are correlated
% addconstr: additional constraints on which spots to consider in the database
%
% Scintillator tilts:
% tiltY: tilt around lab Y (horizontal) axis in degrees (0)
% tiltZ: tilt around lab Z (vertical) axis in degreees (0)
%
% Given spotA, which spotB-s to consider.
%
% Diffraction vector:
% thr_ang: theta angular deviation; absolute value in degrees (0.2)
% corr_rot_to_det: correction of rotation-axis - detector dist. in pixels (0)
%
% Thresholds for full image positions:
% thr_max_offset: CentroidImage offset threshold; absolut value in #images (2)
% thr_ext_offset: ExtStartImage and ExtEndImage offset threshold; absolut value in #images (10)
% thr_genim_offset: at least one of the three images should be offset as maximum this value in #images (1)
%
% Thresholds for difspot images (for search between the limits: basevalue/thr & & basevalue*thr)
% thr_intint: integrated intensity (3)
% thr_area: area of spots inside the bounding box (1.2)
% thr_bbsize: bounding box size (1.2)
%
%
%
% by Peter Reischig, ESRF, 11/2007
%
%%
function gtMatchDifspots(par_match,savepar,show_findings,dopause)
parameters=[];
% modif sabine
load('parameters.mat');
tic
if ~exist('par_match','var')
par_match=[];
end
if ~exist('savepar','var')
savepar=false;
end
if ~exist('show_findings','var')
show_findings=false;
end
if ~exist('dopause','var')
dopause=0;
end
% Default parameters
if isempty(par_match) && ~isfield(parameters,'match')
disp('Using default parameters.')
parameters.match.thr_ang=0.2;
parameters.match.thr_max_offset=2;
parameters.match.thr_ext_offset=10;
parameters.match.thr_genim_offset=1;
parameters.match.corr_rot_to_det=0;
parameters.match.thr_intint=3;
parameters.match.thr_area=1.2;
parameters.match.thr_bbsize=1.2;
parameters.match.tiltY=0;
parameters.match.tiltZ=0;
parameters.match.centmode=0;
parameters.match.addconstr=[];
elseif ~isempty(par_match)
parameters.match=par_match;
end
disp(parameters.match)
if savepar
save([ parameters.acq.dir '/parameters.mat'],'parameters');
disp(' ')
disp('Actual matching parameters are saved in the parameters file.')
end
%% Set up search parameters
% Tolerances for image position
thr_ang=parameters.match.thr_ang; % angular threshold; absolute value in degrees (0.2)
thr_max_offset=parameters.match.thr_max_offset; % CentroidImage offset threshold; absolut value in #images (2)
thr_ext_offset=parameters.match.thr_ext_offset; % ExtStartImage and ExtEndImage offset threshold; absolut value in #images (10)
thr_genim_offset=parameters.match.thr_genim_offset; % at least one of the three images should be offset as maximum this value in #images (1)
corr_rot_to_det=parameters.match.corr_rot_to_det; % correction of rot_to_det in pixels (0)
% Tolerances (for search within the limits: basevalue/thr & & basevalue*thr)
thr_intint=parameters.match.thr_intint; % integrated intensity (3)
thr_area=parameters.match.thr_area; % area of spots inside the bounding box (1.2)
thr_bbsize=parameters.match.thr_bbsize; % bounding box size (1.2)
% Other options
centmode=parameters.match.centmode;
tiltY=parameters.match.tiltY;
tiltZ=parameters.match.tiltZ;
addconstr=parameters.match.addconstr;
if ~isempty(parameters.match.addconstr)
addconstr=[' and ' addconstr];
end
sample_radius=parameters.acq.bb(3)/2;
% Sample geometry in imagesm:
sample_top=parameters.acq.bb(2);
sample_bot=parameters.acq.bb(2)+parameters.acq.bb(4);
rot_to_det=parameters.acq.dist/parameters.acq.pixelsize+corr_rot_to_det;
% Set sample coordinate system
% sorX = on the rotation axis by definition
% sorY = on the rotation axis by definition
sample=gtSampleGeoInSampleSystem(parameters);
sorZ=sample.sorZim; % floor(parameters.acq.ydet/2);
if isfield(parameters.acq, 'difA_name')
difspottable=[parameters.acq.difA_name 'difspot'];
makeprocesstable=[parameters.acq.difA_name 'process_'];
processtable=[parameters.acq.difA_name 'process_difspot'];
else
difspottable=[parameters.acq.name 'difspot'];
makeprocesstable=[parameters.acq.name 'process_'];
processtable=[parameters.acq.name 'process_difspot'];
end
pairtable=parameters.acq.pair_tablename ;
nproj=parameters.acq.nproj; % number of preojections per 180 degrees
nofspotspertable=2000; % no.of spots in the processtable (2000 seems good)
comp_2thetas=parameters.cryst.twoth; % coloumn vector of all valid 2theta angles in degrees
disp(' ')
%close all
gtDBCreateDifspotTable(makeprocesstable,1)
gtDBCreateSpotPairTable(pairtable,1)
nof_spots1=mym(sprintf('select count(*) from %s where CentroidImage between %0.20g and %0.20g %s',difspottable,0,nproj,addconstr))
nof_spots2=mym(sprintf('select count(*) from %s where CentroidImage between %0.20g and %0.20g %s',difspottable,nproj,2*nproj,addconstr))
[spotsA,centim]=mym(sprintf(['select difspotID,CentroidImage from %s where CentroidImage<=%d %s',...
' order by CentroidImage asc, Integral desc'],difspottable,nproj+thr_max_offset,addconstr));
max_difspotID=mym(sprintf('select max(difspotID) from %s',difspottable));
pairedspots=false(max_difspotID,1);
%if isempty(startat)
mym(sprintf('truncate %s',pairtable)) ;
disp('Data in spotpairs table have been deleted, if there were any.') ;
istart=1;
%else
% error('Update this part if needed at all...')
% pairedA=mym(sprintf('select difAID from %s',pairtable));
% pairedB=mym(sprintf('select difBID from %s',pairtable));
% pairedspots(pairedA)=true;
% pairedspots(pairedB)=true;
% [tmp,istart]=min(abs(values_orderby-startat));
%end
pairsfound=0;
show_bands=true;
disp('Processing data...')
disp('Progress:')
disp(' #Spots SpotID A Integral Pairs found Matching');
%% Main loop for spots from first half
for i=istart:length(spotsA)
sp1.id=spotsA(i);
if pairedspots(sp1.id)
continue
end
% create new processtable
if mod(i-istart,nofspotspertable)==0
centim_begin=centim(i)+nproj-thr_max_offset;
centim_end=centim(min([i+nofspotspertable length(centim)]))+nproj+thr_max_offset;
mym(sprintf('TRUNCATE %s',processtable))
mysqlcmd=sprintf(['INSERT INTO %s SELECT difspotID,StartImage,EndImage,MaxImage,ExtStartImage,ExtEndImage, Area, CentroidX, CentroidY,'...
'CentroidImage, BoundingBoxXorigin,BoundingBoxYorigin,BoundingBoxXsize,BoundingBoxYsize,Integral FROM %s WHERE CentroidImage BETWEEN %0.20g AND %0.20g']...
,processtable,difspottable,centim_begin,centim_end);
mym(mysqlcmd);
end
[sp1.maxim,sp1.extstim,sp1.extendim,sp1.area,sp1.centX,sp1.centY,sp1.bbX,sp1.bbY,sp1.bbXsize,sp1.bbYsize,sp1.integral]=...
mym(sprintf(['select CentroidImage,ExtStartImage,ExtEndImage,Area,CentroidX,CentroidY,BoundingBoxXorigin,BoundingBoxYorigin,'...
'BoundingBoxXsize,BoundingBoxYsize,Integral from %s where difspotID=%d %s'],difspottable,sp1.id,addconstr)) ;
pr_bb = projected_sample(sp1.centX,sp1.centY,sample_radius,sample_top,sample_bot,rot_to_det,parameters.acq.rotx) ;
% Set limits for search between lim1..lim2
% CentroidImage, ExtStart, Extend limits for pair
[maxim1,maxim2]=sfCalcBoundaries(sp1.maxim,thr_max_offset,nproj);
[extstim1,extstim2]=sfCalcBoundaries(sp1.extstim,thr_ext_offset,nproj);
[extendim1,extendim2]=sfCalcBoundaries(sp1.extendim,thr_ext_offset,nproj);
% At least one of them should also fulfill thr_genim_offset
[maxim1_gen,maxim2_gen]=sfCalcBoundaries(sp1.maxim,thr_genim_offset,nproj);
[extstim1_gen,extstim2_gen]=sfCalcBoundaries(sp1.extstim,thr_genim_offset,nproj);
[extendim1_gen,extendim2_gen]=sfCalcBoundaries(sp1.extendim,thr_genim_offset,nproj);
% Select candidate spots2 from database by all the criteria
[sps2.id, sps2.maxim, sps2.extstim, sps2.extendim, sps2.area, sps2.centX, sps2.centY, sps2.bbX, sps2.bbY, sps2.bbXsize, sps2.bbYsize, sps2.integral]= ...
mym(sprintf(['SELECT difspotID, CentroidImage, ExtStartImage, ExtEndImage, Area, CentroidX, CentroidY, BoundingBoxXorigin, BoundingBoxYorigin, BoundingBoxXsize, ', ...
'BoundingBoxYsize, Integral FROM %s WHERE ', ...
'(CentroidX between %0.20g and %0.20g) AND (CentroidY between %0.20g and %0.20g) ', ...
'AND (BoundingBoxXsize between %0.20g and %0.20g) AND (BoundingBoxYsize between %0.20g and %0.20g) ', ...
'AND (Integral between %0.20g and %0.20g) ', ...
'AND (Area between %0.20g and %0.20g) ', ...
'AND (CentroidImage between %0.20g and %0.20g) ', ...
'AND (ExtStartImage between %d and %d) ', ...
'AND (ExtEndImage between %d and %d) ', ...
'AND ((CentroidImage between %0.20g and %0.20g) OR ', ...
'(ExtStartImage between %d and %d) OR ', ...
'(ExtEndImage between %d and %d)) ', ...
' %s'], ...
processtable,pr_bb(1),pr_bb(3),pr_bb(2),pr_bb(4), ...
sp1.bbXsize/thr_bbsize, sp1.bbXsize*thr_bbsize, ...
sp1.bbYsize/thr_bbsize, sp1.bbYsize*thr_bbsize, ...
sp1.integral/thr_intint, sp1.integral*thr_intint, ...
sp1.area/thr_area, sp1.area*thr_area, ...
maxim1, maxim2, ...
extstim1, extstim2, ...
extendim1, extendim2, ...
maxim1_gen, maxim2_gen, ...
extstim1_gen, extstim2_gen, ...
extendim1_gen, extendim2_gen, ...
addconstr)) ;
% Delete candidates already paired
sps2.maxim(pairedspots(sps2.id))=[];
sps2.extstim(pairedspots(sps2.id))=[];
sps2.extendim(pairedspots(sps2.id))=[];
sps2.area(pairedspots(sps2.id))=[];
sps2.centX(pairedspots(sps2.id))=[];
sps2.centY(pairedspots(sps2.id))=[];
sps2.bbX(pairedspots(sps2.id))=[];
sps2.bbY(pairedspots(sps2.id))=[];
sps2.bbXsize(pairedspots(sps2.id))=[];
sps2.bbYsize(pairedspots(sps2.id))=[];
sps2.integral(pairedspots(sps2.id))=[];
sps2.id(pairedspots(sps2.id))=[];
%% Theta angle check and selection of best candidate by area, integral, bbXsize, bbYsize, deviation from theoretical theta
sqrsum=Inf ;
for j=1:length(sps2.id)
[theta_OK,angle_diff,thetatype] = sfCheck_thetas(sp1.centX,sp1.centY,sps2.centX(j),sps2.centY(j),...
rot_to_det,tiltY,tiltZ,parameters.acq.rotx,sorZ,comp_2thetas,thr_ang);
if theta_OK==true
sqrsumnew = ((sps2.area(j)-sp1.area)/sp1.area)^2 + ((sps2.bbXsize(j)-sp1.bbXsize)/sp1.bbXsize)^2 + ...
((sps2.bbYsize(j)-sp1.bbYsize)/sp1.bbYsize)^2 + ((sps2.integral(j)-sp1.integral)/sp1.integral)^2 + ...
(min(abs(angle_diff))/thr_ang)^2 + ...
((sps2.extstim(j)-sp1.extstim)/thr_ext_offset)^2 + ((sps2.extendim(j)-sp1.extendim)/thr_ext_offset)^2 ;
if (sqrsumnew < sqrsum)
sqrsum=sqrsumnew ;
pairj=j ;
thetatype_out=thetatype ;
end
end
end
%% Being left with the right pair spot, compute output parameters and loading database
if sqrsum < Inf
pairsfound=pairsfound+1 ;
% Save pair data sp2
sp2.id=sps2.id(pairj);
sp2.maxim=sps2.maxim(pairj);
sp2.extstim=sps2.extstim(pairj);
sp2.extendim=sps2.extendim(pairj);
sp2.area=sps2.area(pairj);
sp2.centX=sps2.centX(pairj);
sp2.centY=sps2.centY(pairj);
sp2.integral=sps2.integral(pairj);
sp2.bbXsize=sps2.bbXsize(pairj);
sp2.bbYsize=sps2.bbYsize(pairj);
sp2.bbX=sps2.bbX(pairj);
sp2.bbY=sps2.bbY(pairj);
% Leave only the wrong candidates in sps2
sps2.maxim(pairj)=[];
sps2.extstim(pairj)=[];
sps2.extendim(pairj)=[];
sps2.area(pairj)=[];
sps2.centX(pairj)=[];
sps2.centY(pairj)=[];
sps2.bbX(pairj)=[];
sps2.bbY(pairj)=[];
sps2.bbXsize(pairj)=[];
sps2.bbYsize(pairj)=[];
sps2.integral(pairj)=[];
sps2.id(pairj)=[];
pairedspots(sp1.id)=true;
pairedspots(sp2.id)=true;
switch centmode
case {0, []}
case 1 % correlate difspots (edf-s or blobs)
[corrx,corry]=gtCorrelateDifspots(sp1.id,sp2.id,parameters,centmode,1);
sp2.centX=sp2.centX-corrx; % the negative of corr is applied, since
% corr referred to the flipped image
sp2.centY=sp2.centY+corry;
case 2 % correlate summed fulls over the difspot bb
[corrx,corry]=gtCorrelateDifspots(sp1.id,sp2.id,parameters,centmode,1);
sp2.centX=sp2.centX-corrx; % the negative of corr is applied, since
% corr referred to the flipped image
sp2.centY=sp2.centY+corry;
end
theta= gtThetaOfPairs(sp1.centX,sp1.centY,sp2.centX,sp2.centY,rot_to_det,tiltY,tiltZ,parameters.acq.rotx,sorZ) ;
eta= gtEtaOfPairs(sp1.centX,sp1.centY,sp2.centX,sp2.centY,tiltY,tiltZ,parameters.acq.rotx,sorZ) ;
diff_vec_setup= gtDiffVecInLab(sp1.centX,sp1.centY,sp2.centX,sp2.centY,rot_to_det,tiltY,tiltZ,parameters.acq.rotx,sorZ) ;
pl_vec_setup= gtDiffPlaneNormInLab(sp1.centX,sp1.centY,sp2.centX,sp2.centY,rot_to_det,tiltY,tiltZ,parameters.acq.rotx,sorZ);
omega1=180/nproj*sp1.maxim;
omega2=180/nproj*sp2.maxim; % (sp2.maxim should always be larger)
omega=(omega1+omega2-180)/2; % the average of the two
[diff_vec_sam(1),diff_vec_sam(2),diff_vec_sam(3)] = gtTrLabToSam(diff_vec_setup(1),diff_vec_setup(2),diff_vec_setup(3),omega) ;
[pl_vec_sam(1),pl_vec_sam(2),pl_vec_sam(3)] = gtTrLabToSam(pl_vec_setup(1),pl_vec_setup(2),pl_vec_setup(3),omega) ;
[lab_cent1(1),lab_cent1(2),lab_cent1(3)] = gtTrScToLab(sp1.centX,sp1.centY,rot_to_det,tiltY,tiltZ,parameters.acq.rotx,sorZ) ;
[lab_cent2(1),lab_cent2(2),lab_cent2(3)] = gtTrScToLab(sp2.centX,sp2.centY,rot_to_det,tiltY,tiltZ,parameters.acq.rotx,sorZ) ;
[samcentA(1),samcentA(2),samcentA(3)] = gtTrLabToSam(lab_cent1(1),lab_cent1(2),lab_cent1(3),omega1) ;
[samcentB(1),samcentB(2),samcentB(3)] = gtTrLabToSam(lab_cent2(1),lab_cent2(2),lab_cent2(3),omega2) ;
av_intint=(sp1.integral+sp2.integral)/2 ;
av_bbXsize=(sp1.bbXsize+sp2.bbXsize)/2 ;
av_bbYsize=(sp1.bbYsize+sp2.bbYsize)/2 ;
% Variables in table
% samcent: difspot centroid in the sample system
% lines connecing the two spots in the sample system
% lor: line origin (centroid of spot A);
% ldir: unit vector of line direction
% pl: plane normal in sample system
mysqlcmd=sprintf(['INSERT INTO %s (difAID, difBID, theta, eta, omega, avintint, avbbXsize, avbbYsize, '...
'samcentXA, samcentYA, samcentZA, samcentXB, samcentYB, samcentZB, '...
'ldirX, ldirY, ldirZ, plX, plY, plZ, thetatype) ',...
'VALUES (%d,%d,%0.20g,%0.20g,%0.20g,%0.20g,%0.20g,%0.20g,%0.20g,%0.20g,%0.20g,%0.20g,%0.20g,%0.20g,%0.20g,%0.20g,%0.20g,%0.20g,%0.20g,%0.20g,%d)'],...
pairtable,sp1.id,sp2.id,theta,eta,omega,av_intint,av_bbXsize,av_bbYsize,...
samcentA(1),samcentA(2),samcentA(3),samcentB(1),samcentB(2),samcentB(3),...
diff_vec_sam(1),diff_vec_sam(2),diff_vec_sam(3),pl_vec_sam(1),pl_vec_sam(2),pl_vec_sam(3),thetatype_out) ;
mym(mysqlcmd);
%% Display findings
if show_findings
gtShowPairFigure(1,parameters,sp1,sp2,pr_bb,comp_2thetas,show_bands,...
sps2,[],[],omega,theta,eta)
show_bands=false;
if strcmpi('button',dopause)
waitforbuttonpress
else
pause(dopause)
end
end
end % pair found
if mod(i,100)==0
disp(sprintf('%10d %10d %10.1f %10d %3.2f%%',i,sp1.id,sp1.integral,pairsfound,pairsfound/(i-istart+1)*100)) ;
end
end % end of main loop for spotsA
disp('Total number of diffraction spots in first half:'), disp(nof_spots1) ;
disp('Total number of diffraction spots in second half:'), disp(nof_spots2) ;
disp('Number of pairs found:'), disp(pairsfound) ;
disp('Total number of pairs in spotpairtable:')
disp(mym(sprintf('Select count(pairID) from %s',pairtable)));
disp(' ');
toc
disp(' ')
gtEvaluateDifspotMatch([],'pair_dist')
end % of function
%% Sub-functions used
% Given an image number, returns the boundaries of two ranges (as image numbers)
% with an offset of 180degree +-tol offset.
function [b1,b2]=sfCalcBoundaries(imno1,tol,nproj)
b1=imno1+nproj-tol;
b2=imno1+nproj+tol;
end
function [theta_OK,diff,thetatype]=sfCheck_thetas(scXA,scYA,scXB,scYB,rottodet,tiltY,tiltZ,rotx,sorZ,comp_2thetas,thr_ang)
% Checks whether the position of two spots gives a diffraction vector which
% corresponds to a valid 2theta angle within the given threshold.
angle=2*gtThetaOfPairs(scXA,scYA,scXB,scYB,rottodet,tiltY,tiltZ,rotx,sorZ);
diff=comp_2thetas-angle;
if all(abs(comp_2thetas-angle) > thr_ang) % checks if all values are non-zero
theta_OK=false ;
thetatype=[] ;
else
theta_OK=true ;
[mini,thetatype]=min(abs(comp_2thetas-angle));
end
end
4_grains/gtMatchDifspots_snow.m deleted 100644 → 0
+ 0
− 489
View file @ 75a2c583
%
% FUNCTION gtMatchDifspots(par_match,centmode,tiltY,tiltZ,addconstr,
% savepar,show_findings,dopause)
%
% (Used to be gtMatchDifspots_v3)
% Matches difspots of a 360 degree scan, from a single difspot table.
% It can take into account tilts and correlate difspots (check if it
% gives any improvement) for better angular precision.
% Speed shouldn't depend much on the number of difspots, since it works
% on a fraction of the spots (creates and uses a smaller spotpairtable,
% named 'processtable').
%
% USAGE e.g. gtMatchDifspots
% gtMatchDifspots(par_match)
% gtMatchDifspots(par_match,1,0,0,'BoundingBoxXsize>10 and BoundingBoxYsize>10')
% gtMatchDifspots(par_match,[],0.1,0.2,[],[],1,'button')
%
% OPTIONAL INPUT
%
% par_match: matching parameters; if not specified it looks into
% the parameters file
% centmode: centroid of difspots based on:
% 0 or []: uses center of mass from difspot table
% 1: difspots from blobs or edf-s are correlated
% 2: difspots - from summed fulls - are correlated
% tiltY: tilt around horizontal axis (in degrees)
% tiltZ: tilt around vertical axis (in degreees)
% addconstr: additional constraints on which spots to consider in the database
% savepar: saves the last used match parameters in the parameters file
% show_findings: if true, the pairs found and candidates are shown in figures
% dopause: pause time in sec after a pair is found; if set as 'button',
% it waits for button press
%
% OUTPUT loads data in spotpairs table given in the parameters file
%
%
% MATCH TOLERANCES parameters.match... and default values in brackets
%
% Given spotA, which spotB-s to consider.
%
% Diffraction vector:
% thr_ang: theta angular deviation; absolute value in degrees (0.2)
% corr_rot_to_det: correction of rotation-axis - detector dist. in pixels (0)
%
% Thresholds for full image positions:
% thr_max_offset: CentroidImage offset threshold; absolut value in #images (2)
% thr_ext_offset: ExtStartImage and ExtEndImage offset threshold; absolut value in #images (10)
% thr_genim_offset: at least one of the three images should be offset as maximum this value in #images (1)
%
% Thresholds for difspot images (for search between the limits: basevalue/thr & & basevalue*thr)
% thr_intint: integrated intensity (3)
% thr_area: area of spots inside the bounding box (1.2)
% thr_bbsize: bounding box size (1.2)
%
%
%
%
% by Peter Reischig, ESRF, 11/2007
%
%%
function gtMatchDifspots(par_match,centmode,tiltY,tiltZ,addconstr,savepar,show_findings,dopause)
parameters=[];
load('parameters.mat');
tic
if ~exist('par_match','var')
par_match=[];
end
if ~exist('centmode','var')
centmode=[];
end
if isempty(centmode)
centmode=0;
end
if ~exist('tiltY','var')
tiltY=0;
end
if ~exist('tiltZ','var')
tiltZ=0;
end
if ~exist('addconstr','var')
addconstr=[];
elseif isempty(addconstr)
addconstr=[];
else
addconstr=[' and ' addconstr];
end
if ~exist('savepar','var')
savepar=false;
end
if ~exist('show_findings','var')
show_findings=false;
end
if ~exist('startat','var')
startat=[];
end
%% Default parameters
if isempty(par_match) && ~isfield(parameters.match,'thr_ang') % modif sabine ||
% replaced by && otherwise it doesnot read the match parameters you add
% in the parameters file
disp('Using default parameters.')
parameters.match.thr_ang=0.2;
parameters.match.thr_max_offset=2;
parameters.match.thr_ext_offset=10;
parameters.match.thr_genim_offset=1;
parameters.match.corr_rot_to_det=0;
parameters.match.thr_intint=3;
parameters.match.thr_area=1.2;
parameters.match.thr_bbsize=1.2;
elseif ~isempty(par_match)
parameters.match=par_match;
end
parameters.match.centmode=centmode;
parameters.match.tiltY=tiltY;
parameters.match.tiltZ=tiltZ;
parameters.match.addconstr=addconstr;
disp(parameters.match)
if savepar
save('parameters.mat','parameters')
disp(' ')
disp('Actual matching parameters are saved in the parameters file.')
end
%% Set up search parameters
% Tolerances for image position
thr_ang=parameters.match.thr_ang; % angular threshold; absolute value in degrees (0.2)
thr_max_offset=parameters.match.thr_max_offset; % CentroidImage offset threshold; absolut value in #images (2)
thr_ext_offset=parameters.match.thr_ext_offset; % ExtStartImage and ExtEndImage offset threshold; absolut value in #images (10)
thr_genim_offset=parameters.match.thr_genim_offset; % at least one of the three images should be offset as maximum this value in #images (1)
corr_rot_to_det=parameters.match.corr_rot_to_det; % correction of rot_to_det in pixels (0)
% Tolerances (for search within the limits: basevalue/thr & & basevalue*thr)
thr_intint=parameters.match.thr_intint; % integrated intensity (3)
thr_area=parameters.match.thr_area; % area of spots inside the bounding box (1.2)
thr_bbsize=parameters.match.thr_bbsize; % bounding box size (1.2)
sample_radius=parameters.acq.bb(3)/2;
% Sample geometry in imagesm:
sample_top=parameters.acq.bb(2);
sample_bot=parameters.acq.bb(2)+parameters.acq.bb(4);
rot_to_det=parameters.acq.dist/parameters.acq.pixelsize+corr_rot_to_det;
% Set sample coordinate system
% sorX = on the rotation axis by definition
% sorY = on the rotation axis by definition
sample=gtSampleGeoInSampleSystem(parameters);
sorZ=sample.sorZim; % floor(parameters.acq.ydet/2);
if isfield(parameters.acq, 'difA_name')
difspottable=[parameters.acq.difA_name 'difspot'];
makeprocesstable=[parameters.acq.difA_name 'process_'];
processtable=[parameters.acq.difA_name 'process_difspot'];
else
difspottable=[parameters.acq.name 'difspot'];
makeprocesstable=[parameters.acq.name 'process_'];
processtable=[parameters.acq.name 'process_difspot'];
end
pairtable=parameters.acq.pair_tablename ;
nproj=parameters.acq.nproj; % number of preojections per 180 degrees
nofspotspertable=2000; % no.of spots in the processtable (2000 seems good)
comp_2thetas=parameters.cryst.twoth; % coloumn vector of all valid 2theta angles in degrees
disp(' ')
%close all
gtDBCreateDifspotTable(makeprocesstable,1)
gtDBCreateSpotPairTable(pairtable,1)
nof_spots1=mym(sprintf('select count(*) from %s where CentroidImage between %0.20g and %0.20g %s and Area>300',difspottable,0,nproj,addconstr));
nof_spots2=mym(sprintf('select count(*) from %s where CentroidImage between %0.20g and %0.20g %s and Area>300',difspottable,nproj,2*nproj,addconstr));
[spotsA,centim]=mym(sprintf(['select difspotID,CentroidImage from %s where CentroidImage<=%d %s and Area>300',...
' order by CentroidImage asc, Integral desc'],difspottable,nproj+thr_max_offset,addconstr));
max_difspotID=mym(sprintf('select max(difspotID) from %s',difspottable));
pairedspots=false(max_difspotID,1);
if isempty(startat)
mym(sprintf('truncate %s',pairtable)) ;
disp('Data in spotpairs table have been deleted, if there were any.') ;
istart=1;
else
error('Update this part if needed at all...')
% pairedA=mym(sprintf('select difAID from %s',pairtable));
% pairedB=mym(sprintf('select difBID from %s',pairtable));
% pairedspots(pairedA)=true;
% pairedspots(pairedB)=true;
% [tmp,istart]=min(abs(values_orderby-startat));
end
pairsfound=0;
show_bands=true;
disp('Processing data...')
disp('Progress:')
disp(' #Spots SpotID A Integral Pairs found Matching');
%% Main loop for spots from first half
for i=istart:length(spotsA)
sp1.id=spotsA(i);
if pairedspots(sp1.id)
continue
end
% create new processtable
if mod(i-istart,nofspotspertable)==0
centim_begin=centim(i)+nproj-thr_max_offset;
centim_end=centim(min([i+nofspotspertable length(centim)]))+nproj+thr_max_offset;
mym(sprintf('TRUNCATE %s',processtable))
mysqlcmd=sprintf(['INSERT INTO %s SELECT difspotID,StartImage,EndImage,MaxImage,ExtStartImage,ExtEndImage, Area, CentroidX, CentroidY,'...
'CentroidImage, BoundingBoxXorigin,BoundingBoxYorigin,BoundingBoxXsize,BoundingBoxYsize,Integral FROM %s WHERE CentroidImage BETWEEN %0.20g AND %0.20g']...
,processtable,difspottable,centim_begin,centim_end);
mym(mysqlcmd);
end
[sp1.maxim,sp1.extstim,sp1.extendim,sp1.area,sp1.centX,sp1.centY,sp1.bbX,sp1.bbY,sp1.bbXsize,sp1.bbYsize,sp1.integral]=...
mym(sprintf(['select CentroidImage,ExtStartImage,ExtEndImage,Area,CentroidX,CentroidY,BoundingBoxXorigin,BoundingBoxYorigin,'...
'BoundingBoxXsize,BoundingBoxYsize,Integral from %s where difspotID=%d %s'],difspottable,sp1.id,addconstr)) ;
pr_bb = projected_sample(sp1.centX,sp1.centY,sample_radius,sample_top,sample_bot,rot_to_det,parameters.acq.rotx) ;
% Set limits for search between lim1..lim2
% CentroidImage, ExtStart, Extend limits for pair
[maxim1,maxim2]=sfCalcBoundaries(sp1.maxim,thr_max_offset,nproj);
[extstim1,extstim2]=sfCalcBoundaries(sp1.extstim,thr_ext_offset,nproj);
[extendim1,extendim2]=sfCalcBoundaries(sp1.extendim,thr_ext_offset,nproj);
% At least one of them should also fulfill thr_genim_offset
[maxim1_gen,maxim2_gen]=sfCalcBoundaries(sp1.maxim,thr_genim_offset,nproj);
[extstim1_gen,extstim2_gen]=sfCalcBoundaries(sp1.extstim,thr_genim_offset,nproj);
[extendim1_gen,extendim2_gen]=sfCalcBoundaries(sp1.extendim,thr_genim_offset,nproj);
% Select candidate spots2 from database by all the criteria
[sps2.id, sps2.maxim, sps2.extstim, sps2.extendim, sps2.area, sps2.centX, sps2.centY, sps2.bbX, sps2.bbY, sps2.bbXsize, sps2.bbYsize, sps2.integral]= ...
mym(sprintf(['SELECT difspotID, CentroidImage, ExtStartImage, ExtEndImage, Area, CentroidX, CentroidY, BoundingBoxXorigin, BoundingBoxYorigin, BoundingBoxXsize, ', ...
'BoundingBoxYsize, Integral FROM %s WHERE ', ...
'(CentroidX between %0.20g and %0.20g) AND (CentroidY between %0.20g and %0.20g) ', ...
'AND (BoundingBoxXsize between %0.20g and %0.20g) AND (BoundingBoxYsize between %0.20g and %0.20g) ', ...
'AND (Integral between %0.20g and %0.20g) ', ...
'AND (Area between %0.20g and %0.20g) ', ...
'AND (CentroidImage between %0.20g and %0.20g) ', ...
'AND (ExtStartImage between %d and %d) ', ...
'AND (ExtEndImage between %d and %d) ', ...
'AND ((CentroidImage between %0.20g and %0.20g) OR ', ...
'(ExtStartImage between %d and %d) OR ', ...
'(ExtEndImage between %d and %d)) ', ...
' %s'], ...
processtable,pr_bb(1),pr_bb(3),pr_bb(2),pr_bb(4), ...
sp1.bbXsize/thr_bbsize, sp1.bbXsize*thr_bbsize, ...
sp1.bbYsize/thr_bbsize, sp1.bbYsize*thr_bbsize, ...
sp1.integral/thr_intint, sp1.integral*thr_intint, ...
sp1.area/thr_area, sp1.area*thr_area, ...
maxim1, maxim2, ...
extstim1, extstim2, ...
extendim1, extendim2, ...
maxim1_gen, maxim2_gen, ...
extstim1_gen, extstim2_gen, ...
extendim1_gen, extendim2_gen, ...
addconstr)) ;
% Delete candidates already paired
sps2.maxim(pairedspots(sps2.id))=[];
sps2.extstim(pairedspots(sps2.id))=[];
sps2.extendim(pairedspots(sps2.id))=[];
sps2.area(pairedspots(sps2.id))=[];
sps2.centX(pairedspots(sps2.id))=[];
sps2.centY(pairedspots(sps2.id))=[];
sps2.bbX(pairedspots(sps2.id))=[];
sps2.bbY(pairedspots(sps2.id))=[];
sps2.bbXsize(pairedspots(sps2.id))=[];
sps2.bbYsize(pairedspots(sps2.id))=[];
sps2.integral(pairedspots(sps2.id))=[];
sps2.id(pairedspots(sps2.id))=[];
%% Theta angle check and selection of best candidate by area, integral, bbXsize, bbYsize, deviation from theoretical theta
sqrsum=Inf ;
for j=1:length(sps2.id)
[theta_OK,angle_diff,thetatype] = sfCheck_thetas(sp1.centX,sp1.centY,sps2.centX(j),sps2.centY(j),...
rot_to_det,tiltY,tiltZ,parameters.acq.rotx,sorZ,comp_2thetas,thr_ang);
if theta_OK==true
sqrsumnew = ((sps2.area(j)-sp1.area)/sp1.area)^2 + ((sps2.bbXsize(j)-sp1.bbXsize)/sp1.bbXsize)^2 + ...
((sps2.bbYsize(j)-sp1.bbYsize)/sp1.bbYsize)^2 + ((sps2.integral(j)-sp1.integral)/sp1.integral)^2 + ...
(min(abs(angle_diff))/thr_ang)^2 + ...
((sps2.extstim(j)-sp1.extstim)/thr_ext_offset)^2 + ((sps2.extendim(j)-sp1.extendim)/thr_ext_offset)^2 ;
if (sqrsumnew < sqrsum)
sqrsum=sqrsumnew ;
pairj=j ;
thetatype_out=thetatype ;
end
end
end
%% Being left with the right pair spot, compute output parameters and loading database
if sqrsum < Inf
pairsfound=pairsfound+1 ;
% Save pair data sp2
sp2.id=sps2.id(pairj);
sp2.maxim=sps2.maxim(pairj);
sp2.extstim=sps2.extstim(pairj);
sp2.extendim=sps2.extendim(pairj);
sp2.area=sps2.area(pairj);
sp2.centX=sps2.centX(pairj);
sp2.centY=sps2.centY(pairj);
sp2.integral=sps2.integral(pairj);
sp2.bbXsize=sps2.bbXsize(pairj);
sp2.bbYsize=sps2.bbYsize(pairj);
sp2.bbX=sps2.bbX(pairj);
sp2.bbY=sps2.bbY(pairj);
% Leave only the wrong candidates in sps2
sps2.maxim(pairj)=[];
sps2.extstim(pairj)=[];
sps2.extendim(pairj)=[];
sps2.area(pairj)=[];
sps2.centX(pairj)=[];
sps2.centY(pairj)=[];
sps2.bbX(pairj)=[];
sps2.bbY(pairj)=[];
sps2.bbXsize(pairj)=[];
sps2.bbYsize(pairj)=[];
sps2.integral(pairj)=[];
sps2.id(pairj)=[];
pairedspots(sp1.id)=true;
pairedspots(sp2.id)=true;
switch centmode
case {0, []}
case 1 % correlate difspots (edf-s or blobs)
[corrx,corry]=gtCorrelateDifspots(sp1.id,sp2.id,parameters,centmode,1);
sp2.centX=sp2.centX-corrx; % the negative of corr is applied, since
% corr referred to the flipped image
sp2.centY=sp2.centY+corry;
case 2 % correlate summed fulls over the difspot bb
[corrx,corry]=gtCorrelateDifspots(sp1.id,sp2.id,parameters,centmode,1);
sp2.centX=sp2.centX-corrx; % the negative of corr is applied, since
% corr referred to the flipped image
sp2.centY=sp2.centY+corry;
end
theta= gtThetaOfPairs(sp1.centX,sp1.centY,sp2.centX,sp2.centY,rot_to_det,tiltY,tiltZ,parameters.acq.rotx,sorZ) ;
eta= gtEtaOfPairs(sp1.centX,sp1.centY,sp2.centX,sp2.centY,tiltY,tiltZ,parameters.acq.rotx,sorZ) ;
diff_vec_setup= gtDiffVecInLab(sp1.centX,sp1.centY,sp2.centX,sp2.centY,rot_to_det,tiltY,tiltZ,parameters.acq.rotx,sorZ) ;
pl_vec_setup= gtDiffPlaneNormInLab(sp1.centX,sp1.centY,sp2.centX,sp2.centY,rot_to_det,tiltY,tiltZ,parameters.acq.rotx,sorZ);
omega1=180/nproj*sp1.maxim;
omega2=180/nproj*sp2.maxim; % (sp2.maxim should always be larger)
omega=(omega1+omega2-180)/2; % the average of the two
[diff_vec_sam(1),diff_vec_sam(2),diff_vec_sam(3)] = gtTrLabToSam(diff_vec_setup(1),diff_vec_setup(2),diff_vec_setup(3),omega) ;
[pl_vec_sam(1),pl_vec_sam(2),pl_vec_sam(3)] = gtTrLabToSam(pl_vec_setup(1),pl_vec_setup(2),pl_vec_setup(3),omega) ;
[lab_cent1(1),lab_cent1(2),lab_cent1(3)] = gtTrScToLab(sp1.centX,sp1.centY,rot_to_det,tiltY,tiltZ,parameters.acq.rotx,sorZ) ;
[lab_cent2(1),lab_cent2(2),lab_cent2(3)] = gtTrScToLab(sp2.centX,sp2.centY,rot_to_det,tiltY,tiltZ,parameters.acq.rotx,sorZ) ;
[samcentA(1),samcentA(2),samcentA(3)] = gtTrLabToSam(lab_cent1(1),lab_cent1(2),lab_cent1(3),omega1) ;
[samcentB(1),samcentB(2),samcentB(3)] = gtTrLabToSam(lab_cent2(1),lab_cent2(2),lab_cent2(3),omega2) ;
av_intint=(sp1.integral+sp2.integral)/2 ;
av_bbXsize=(sp1.bbXsize+sp2.bbXsize)/2 ;
av_bbYsize=(sp1.bbYsize+sp2.bbYsize)/2 ;
% Variables in table
% samcent: difspot centroid in the sample system
% lines connecing the two spots in the sample system
% lor: line origin (centroid of spot A);
% ldir: unit vector of line direction
% pl: plane normal in sample system
mysqlcmd=sprintf(['INSERT INTO %s (difAID, difBID, theta, eta, omega, avintint, avbbXsize, avbbYsize, '...
'samcentXA, samcentYA, samcentZA, samcentXB, samcentYB, samcentZB, '...
'ldirX, ldirY, ldirZ, plX, plY, plZ, thetatype) ',...
'VALUES (%d,%d,%0.20g,%0.20g,%0.20g,%0.20g,%0.20g,%0.20g,%0.20g,%0.20g,%0.20g,%0.20g,%0.20g,%0.20g,%0.20g,%0.20g,%0.20g,%0.20g,%0.20g,%0.20g,%d)'],...
pairtable,sp1.id,sp2.id,theta,eta,omega,av_intint,av_bbXsize,av_bbYsize,...
samcentA(1),samcentA(2),samcentA(3),samcentB(1),samcentB(2),samcentB(3),...
diff_vec_sam(1),diff_vec_sam(2),diff_vec_sam(3),pl_vec_sam(1),pl_vec_sam(2),pl_vec_sam(3),thetatype_out) ;
mym(mysqlcmd);
%% Display findings
if show_findings
gtShowPairFigure(1,parameters,sp1,sp2,pr_bb,comp_2thetas,show_bands,...
sps2,[],[],omega,theta,eta)
show_bands=false;
if strcmpi('button',dopause)
waitforbuttonpress
else
pause(dopause)
end
end
end % pair found
if mod(i,100)==0
disp(sprintf('%10d %10d %10.1f %10d %3.2f%%',i,sp1.id,sp1.integral,pairsfound,pairsfound/(i-istart+1)*100)) ;
end
end % end of main loop for spotsA
disp('Total number of diffraction spots in first half:'), disp(nof_spots1) ;
disp('Total number of diffraction spots in second half:'), disp(nof_spots2) ;
disp('Number of pairs found:'), disp(pairsfound) ;
disp('Total number of pairs in spotpairtable:')
disp(mym(sprintf('Select count(pairID) from %s',pairtable)));
disp(' ');
toc
disp(' ')
gtEvaluateDifspotMatch([],'pair_dist')
end % of function
%% Sub-functions used
% Given an image number, returns the boundaries of two ranges (as image numbers)
% with an offset of 180degree +-tol offset.
function [b1,b2]=sfCalcBoundaries(imno1,tol,nproj)
b1=imno1+nproj-tol;
b2=imno1+nproj+tol;
end
function [theta_OK,diff,thetatype]=sfCheck_thetas(scXA,scYA,scXB,scYB,rottodet,tiltY,tiltZ,rotx,sorZ,comp_2thetas,thr_ang)
% Checks whether the position of two spots gives a diffraction vector which
% corresponds to a valid 2theta angle within the given threshold.
angle=2*gtThetaOfPairs(scXA,scYA,scXB,scYB,rottodet,tiltY,tiltZ,rotx,sorZ);
diff=comp_2thetas-angle;
if all(abs(comp_2thetas-angle) > thr_ang) % checks if all values are non-zero
theta_OK=false ;
thetatype=[] ;
else
theta_OK=true ;
[mini,thetatype]=min(abs(comp_2thetas-angle));
end
end
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment