From 2f3ce06395a0a9f654d1481a595cd1f2faf3d96c Mon Sep 17 00:00:00 2001
From: Peter Reischig <peter.reischig@esrf.fr>
Date: Tue, 12 Jun 2012 15:52:03 +0000
Subject: [PATCH] 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
---
4_grains/gt2LinesAngle.m | 13 -
4_grains/gt2LinesDist.m | 27 --
4_grains/gt2LinesIntersection.m | 50 ---
4_grains/gtAngConsCosMat.m | 12 -
4_grains/gtAngConsCosVec.m | 14 -
4_grains/gtAngConsMat.m | 17 --
4_grains/gtAngularCons.m | 43 ---
4_grains/gtAngularConsCos.m | 17 --
4_grains/gtCheckGroupCons.m | 16 -
4_grains/gtCheckPairCons.m | 67 ----
4_grains/gtDiffPlaneNormInLab.m | 42 ---
4_grains/gtDiffVecInLab.m | 45 ---
4_grains/gtEtaOfPairs.m | 71 -----
4_grains/gtEtaOfPoint.m | 41 ---
4_grains/gtEvaluateDifspotMatch.m | 335 --------------------
4_grains/gtMatchDifspots.m | 474 -----------------------------
4_grains/gtMatchDifspots_snow.m | 489 ------------------------------
17 files changed, 1773 deletions(-)
delete mode 100644 4_grains/gt2LinesAngle.m
delete mode 100644 4_grains/gt2LinesDist.m
delete mode 100644 4_grains/gt2LinesIntersection.m
delete mode 100644 4_grains/gtAngConsCosMat.m
delete mode 100644 4_grains/gtAngConsCosVec.m
delete mode 100644 4_grains/gtAngConsMat.m
delete mode 100644 4_grains/gtAngularCons.m
delete mode 100644 4_grains/gtAngularConsCos.m
delete mode 100644 4_grains/gtCheckGroupCons.m
delete mode 100644 4_grains/gtCheckPairCons.m
delete mode 100644 4_grains/gtDiffPlaneNormInLab.m
delete mode 100644 4_grains/gtDiffVecInLab.m
delete mode 100644 4_grains/gtEtaOfPairs.m
delete mode 100644 4_grains/gtEtaOfPoint.m
delete mode 100644 4_grains/gtEvaluateDifspotMatch.m
delete mode 100644 4_grains/gtMatchDifspots.m
delete mode 100644 4_grains/gtMatchDifspots_snow.m
diff --git a/4_grains/gt2LinesAngle.m b/4_grains/gt2LinesAngle.m
deleted file mode 100644
index 95a814d5..00000000
--- a/4_grains/gt2LinesAngle.m
+++ /dev/null
@@ -1,13 +0,0 @@
-%
-% 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
diff --git a/4_grains/gt2LinesDist.m b/4_grains/gt2LinesDist.m
deleted file mode 100644
index e46a75e2..00000000
--- a/4_grains/gt2LinesDist.m
+++ /dev/null
@@ -1,27 +0,0 @@
-% 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
diff --git a/4_grains/gt2LinesIntersection.m b/4_grains/gt2LinesIntersection.m
deleted file mode 100644
index d21fc8f2..00000000
--- a/4_grains/gt2LinesIntersection.m
+++ /dev/null
@@ -1,50 +0,0 @@
-%
-% 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
diff --git a/4_grains/gtAngConsCosMat.m b/4_grains/gtAngConsCosMat.m
deleted file mode 100644
index 99d4d27c..00000000
--- a/4_grains/gtAngConsCosMat.m
+++ /dev/null
@@ -1,12 +0,0 @@
-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
diff --git a/4_grains/gtAngConsCosVec.m b/4_grains/gtAngConsCosVec.m
deleted file mode 100644
index e9d7e725..00000000
--- a/4_grains/gtAngConsCosVec.m
+++ /dev/null
@@ -1,14 +0,0 @@
-
-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
diff --git a/4_grains/gtAngConsMat.m b/4_grains/gtAngConsMat.m
deleted file mode 100644
index d94853a7..00000000
--- a/4_grains/gtAngConsMat.m
+++ /dev/null
@@ -1,17 +0,0 @@
-
-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
diff --git a/4_grains/gtAngularCons.m b/4_grains/gtAngularCons.m
deleted file mode 100644
index 28b16c06..00000000
--- a/4_grains/gtAngularCons.m
+++ /dev/null
@@ -1,43 +0,0 @@
-
-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
diff --git a/4_grains/gtAngularConsCos.m b/4_grains/gtAngularConsCos.m
deleted file mode 100644
index 86bf05e0..00000000
--- a/4_grains/gtAngularConsCos.m
+++ /dev/null
@@ -1,17 +0,0 @@
-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
diff --git a/4_grains/gtCheckGroupCons.m b/4_grains/gtCheckGroupCons.m
deleted file mode 100644
index 254e5416..00000000
--- a/4_grains/gtCheckGroupCons.m
+++ /dev/null
@@ -1,16 +0,0 @@
-
-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
diff --git a/4_grains/gtCheckPairCons.m b/4_grains/gtCheckPairCons.m
deleted file mode 100644
index 82f3e008..00000000
--- a/4_grains/gtCheckPairCons.m
+++ /dev/null
@@ -1,67 +0,0 @@
-%
-% 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
-
diff --git a/4_grains/gtDiffPlaneNormInLab.m b/4_grains/gtDiffPlaneNormInLab.m
deleted file mode 100644
index 3f3eaeeb..00000000
--- a/4_grains/gtDiffPlaneNormInLab.m
+++ /dev/null
@@ -1,42 +0,0 @@
-%
-% 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
diff --git a/4_grains/gtDiffVecInLab.m b/4_grains/gtDiffVecInLab.m
deleted file mode 100644
index 9598c06f..00000000
--- a/4_grains/gtDiffVecInLab.m
+++ /dev/null
@@ -1,45 +0,0 @@
-%
-% 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
diff --git a/4_grains/gtEtaOfPairs.m b/4_grains/gtEtaOfPairs.m
deleted file mode 100644
index f47dd458..00000000
--- a/4_grains/gtEtaOfPairs.m
+++ /dev/null
@@ -1,71 +0,0 @@
-%
-% 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
diff --git a/4_grains/gtEtaOfPoint.m b/4_grains/gtEtaOfPoint.m
deleted file mode 100644
index ebd81dd0..00000000
--- a/4_grains/gtEtaOfPoint.m
+++ /dev/null
@@ -1,41 +0,0 @@
-% 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
diff --git a/4_grains/gtEvaluateDifspotMatch.m b/4_grains/gtEvaluateDifspotMatch.m
deleted file mode 100644
index 8c468a35..00000000
--- a/4_grains/gtEvaluateDifspotMatch.m
+++ /dev/null
@@ -1,335 +0,0 @@
-%
-% 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
-
diff --git a/4_grains/gtMatchDifspots.m b/4_grains/gtMatchDifspots.m
deleted file mode 100644
index 4f6bdc97..00000000
--- a/4_grains/gtMatchDifspots.m
+++ /dev/null
@@ -1,474 +0,0 @@
-%
-% 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
diff --git a/4_grains/gtMatchDifspots_snow.m b/4_grains/gtMatchDifspots_snow.m
deleted file mode 100644
index 9468f93d..00000000
--- a/4_grains/gtMatchDifspots_snow.m
+++ /dev/null
@@ -1,489 +0,0 @@
-%
-% 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
--
GitLab