diff --git a/4_grains/GetRotAxeAngle.m b/4_grains/GetRotAxeAngle.m
deleted file mode 100644
index 01c5bb4b33f84588c946515e88df518703a4d101..0000000000000000000000000000000000000000
--- a/4_grains/GetRotAxeAngle.m
+++ /dev/null
@@ -1,70 +0,0 @@
-function [L0,Phi] = GetRotAxeAngle(T)
-% Function [L0,Phi] = GetRotAxeAngle(T)
-% for orthogonal matrix T(3,3) with det(T)=+1
-% return the orth of rotation L0(3,1) and
-% the angle of rotation Phi in [0,pi].
-% Author: Sergiy Iglin
-% e-mail: iglin@kpi.kharkov.ua
-% or: siglin@yandex.ru
-% personal page: http://iglin.exponenta.ru
-
-%========== Checking of datas ==================
-strerr='Is needed 1 input parameter: the orthogonal matrix T(3,3) with det(T)=+1!';
-tol=1e5*eps; % tolerance
-if nargin<1,
- error(strerr);
-end
-if ~isnumeric(T),
- error(strerr);
-end
-d=size(T);
-if ~(length(d)==2),
- error(strerr);
-end
-if ~all(d==3),
- error(strerr);
-end
-if max(max(abs(T'*T-eye(3))))>tol,
- error(strerr);
-end
-if abs(det(T)-1)>tol,
- error(strerr);
-end
-
-%============= Solution ================
-T1=T-eye(3);
-Lm=[cross(T1(:,1),T1(:,2)),cross(T1(:,2),T1(:,3)),...
- cross(T1(:,3),T1(:,1))];
-mLm=max(abs(Lm));
-[nmLm,inmLm]=max(mLm);
-L=Lm(:,inmLm);
-normL=norm(L);
-if normL<eps, % without rotation
- Phi=0;
- L0=zeros(3,1);
-else
- L0=L/normL; % orth
- if norm(L0-[1;0;0])<tol, % L0==i
- a=[0;1;0];
- b=T(:,2);
- else
- a=[1;0;0];
- b=T(:,1);
- end
- a1=a-(a'*L0)*L0;
- b1=cross(L0,a);
- c1=b-(a'*L0)*L0;
- r1=(a1.^2+b1.^2).^0.5;
- [r0,ir0]=max(r1);
- beta=atan2(b1(ir0),a1(ir0));
- ac=acos(c1(ir0)/r0);
- Phi=ac+beta;
- if norm(RotVecArAxe(a,L0,Phi)-b)>tol,
- Phi=-ac+beta;
- end
- if Phi>pi,
- L0=-L0;
- Phi=2*pi-Phi;
- end
-end
-return
\ No newline at end of file
diff --git a/4_grains/RotVecArAxe.m b/4_grains/RotVecArAxe.m
deleted file mode 100644
index f84f007067ce164a987cfdbaab2a13f4cf490aae..0000000000000000000000000000000000000000
--- a/4_grains/RotVecArAxe.m
+++ /dev/null
@@ -1,43 +0,0 @@
-function B = RotVecArAxe(A,L,Phi)
-% Function B = RotVecArAxe(L,Phi,A) rotate
-% the vector A(3,1) around the axe L(3,1)
-% into angle Phi radian counterclockwise.
-% Author: Sergiy Iglin
-% e-mail: iglin@kpi.kharkov.ua
-% or: siglin@yandex.ru
-% personal page: http://iglin.exponenta.ru
-
-%========== Checking of datas ==================
-strerr='Is needed 3 input parameters: vector A(3,1), axe L(3,1) and angle Phi!';
-if nargin<3,
- error(strerr);
-end
-if ~isnumeric(A)
- error(strerr);
-else
- A=A(:);
- if length(A)<3,
- error(strerr);
- end
- A=A(1:3);
-end
-if ~isnumeric(L)
- error(strerr);
-else
- L=L(:);
- if length(L)<3,
- error(strerr);
- end
- L=L(1:3);
- L0=L/norm(L); % orth
-end
-if ~isnumeric(Phi)
- error(strerr);
-else
- Phi=Phi(1);
-end
-
-%============= Solution ================
-cphi=cos(Phi);
-B=A*cphi+(A'*L0)*(1-cphi)*L0+cross(L0,A)*sin(Phi);
-return
\ No newline at end of file
diff --git a/4_grains/RotVecArAxed.m b/4_grains/RotVecArAxed.m
deleted file mode 100644
index b4139502aeec081784bb543c4def217009af71ab..0000000000000000000000000000000000000000
--- a/4_grains/RotVecArAxed.m
+++ /dev/null
@@ -1,43 +0,0 @@
-function B = RotVecArAxed(A,L,Phi)
-% Function B = RotVecArAxe(L,Phi,A) rotate
-% the vector A(3,1) around the axe L(3,1)
-% into angle Phi radian counterclockwise.
-% Author: Sergiy Iglin
-% e-mail: iglin@kpi.kharkov.ua
-% or: siglin@yandex.ru
-% personal page: http://iglin.exponenta.ru
-
-%========== Checking of datas ==================
-strerr='Is needed 3 input parameters: vector A(3,1), axe L(3,1) and angle Phi!';
-if nargin<3,
- error(strerr);
-end
-if ~isnumeric(A)
- error(strerr);
-else
- A=A(:);
- if length(A)<3,
- error(strerr);
- end
- A=A(1:3);
-end
-if ~isnumeric(L)
- error(strerr);
-else
- L=L(:);
- if length(L)<3,
- error(strerr);
- end
- L=L(1:3);
- L0=L/norm(L); % orth
-end
-if ~isnumeric(Phi)
- error(strerr);
-else
- Phi=Phi(1);
-end
-
-%============= Solution ================
-cphi=cosd(Phi);
-B=A*cphi+(A'*L0)*(1-cphi)*L0+cross(L0,A)*sind(Phi);
-return
\ No newline at end of file
diff --git a/4_grains/combine_grains.m b/4_grains/combine_grains.m
deleted file mode 100644
index 833b0510c796338cf60d657bba289603750b9893..0000000000000000000000000000000000000000
--- a/4_grains/combine_grains.m
+++ /dev/null
@@ -1,80 +0,0 @@
-function output = combine_grains(r_vectors,positions)
-
-%r_vectors = [grainid rx ry rz]
-%positions = [grainid grainx grainy grainz grainsize]
-%produced by gtReadGrains
-
-%avoid using this old version by mistake!
-check=inputwdefault('Should use version 2!!! Continue with this old version? [y/n]', 'n');
-if check=='n'
- return;
-end
-
-output=zeros(size(r_vectors,1),50);%assume no more than 50 grains to be combined...
-
-for i=1:size(r_vectors,1)
-
- dif = r_vectors(:,2:4)-repmat(r_vectors(i,2:4),size(r_vectors,1),1);
- dif = dif.^2;
- r_distance = sqrt(sum(dif,2));
-
- dif = positions(:,2:4)-repmat(positions(i,2:4),size(r_vectors,1),1);
- dif = dif.^2;
- p_distance = sqrt(sum(dif,2));
- %consider the average size of the two grains
- p_distance = p_distance./(0.5*(positions(:,5)+repmat(positions(i,5),size(r_vectors,1),1)));
-
- %close in orientation and distance
- dummy =[];
- %dummy = [dummy; find(r_distance < 0.03 & p_distance < 0.2)];
- %if super close in orientation, can relax distance requirement to allow
- %for irregular shapes
-dummy = [dummy; find(r_distance < 0.03 & p_distance < 0.4)];
- dummy = [dummy; find(p_distance < 0.2)];
-
- dummy = unique(dummy);
-
- output(i,1:(length(dummy))) = r_vectors(dummy,1)';
-
-end
-
-%crop output to only the required number of columns
-tmp = sum(output,1);
-dummy=find(tmp==0);
-output(:,dummy)=[];
-
-%get the unique rows - grains to be combined
-output = unique(output,'rows');
-%can still have multiple rows that refer to a single grain, eg ...;15 21
-%0;15 21 35;... in s5_dct5. Remove these...
-tmp = unique(output(:,1));
-output2=[];
-
-for i=1:length(tmp)
- dummy=find(output(:,1)==tmp(i));
- output2(i,:)=output(max(dummy),:);%take only the row with the most entries
-end
-output=output2;
-
-%also need to ensure that a grain only occurs once in the matrix
-maxid = max(output(:));%max grainid
-del_rows=[];
-for i=1:maxid;
- dummy=find(output==i);
-
- if length(dummy)>1
- [r,c]=ind2sub(size(output),dummy);%need to remove the entry in the first column
- dummy2=find(c==1);
- del_rows=[del_rows r(dummy2)];
- end
-
-end
-
-del_rows=unique(del_rows);
-output(del_rows,:)=[];
-
-
-
-
-%output is now a table of grains that should be grouped together...
-
diff --git a/4_grains/combine_grains2.m b/4_grains/combine_grains2.m
deleted file mode 100644
index fcf4547606373910be0daef5c5f816b0f87e6fdd..0000000000000000000000000000000000000000
--- a/4_grains/combine_grains2.m
+++ /dev/null
@@ -1,218 +0,0 @@
-function combines = combine_grains2(r_vectors,positions, bads)
-% <<<< NOT USED >>>>
-%r_vectors = [grainid rx ry rz]
-%positions = [grainid grainx grainy grainz grainsize]
-%produced by gtReadGrains
-%combine_grains2 - need to improve the handling of what combines with
-%what...
-%apparently need to use positions calculated from the reconstructed,
-%post-processed grains - the x/y/zcenters are rubbish for this, need
-%x/y/zcenter2s
-%check consistancy because of the danger of combining twins, or other close
-%grains
-
-%get rid of bads grains before starting, also those already merged
-
-if 1
-dummy=find(bads==1 | bads==2); %consider possibly bad grains
-
-r_vectors(dummy,:)=[];
-positions(dummy,:)=[];
-else
- disp('not dealing with bad grains!')
-end
-
-output=zeros(size(r_vectors,1),50);%assume no more than 50 grains to be combined...
-
-for i=1:size(r_vectors,1)
-
- dif = r_vectors(:,2:4)-repmat(r_vectors(i,2:4),size(r_vectors,1),1);
- dif = dif.^2;
- r_distance = sqrt(sum(dif,2));
-
- dif = positions(:,2:4)-repmat(positions(i,2:4),size(r_vectors,1),1);
- dif = dif.^2;
- p_distance = sqrt(sum(dif,2));
- %consider the average size of the two grains
- p_distance = p_distance./(0.5*(positions(:,5)+repmat(positions(i,5),size(r_vectors,1),1)));
-
- %close in orientation and distance
- dummy =[];
- %dummy = [dummy; find(r_distance < 0.03 & p_distance < 0.2)];
- %if super close in orientation, can relax distance requirement to allow
- %for irregular shapes
- dummy = [dummy; find(r_distance < 0.03 & p_distance < 0.4)];
- %can have found the wrong r space solution - will test consistancy later
- dummy = [dummy; find(p_distance < 0.2)];
-
- dummy = unique(dummy);
-
- output(i,1:(length(dummy))) = r_vectors(dummy,1)';
-
-end
-
-
-
-%crop output to only the required number of columns
-tmp = sum(output,1);
-dummy=find(tmp==0);
-output(:,dummy)=[];
-tmp = sum(output,2);
-dummy=find(tmp==0);
-output(dummy,:)=[];
-
-%get the unique rows - grains to be combined
-output = unique(output,'rows');
-
-
-%can still have multiple rows that refer to a single grain, eg ...;15 21
-%0;15 21 35;... in s5_dct5. Remove these...
-%want to keep the largest possible grain sets
-
-output2=zeros(length(output),50); %assume no more than 50 grains to be combined...
-count=1;
-test=0;
-test_length_old=size(output,1);
-
-while test==0
-
-for i=1:max(output(:)) % go though all grainids
-
-
- dummy=find(output==i);
- [a,b]=ind2sub(size(output),dummy); % a are the rows to combine
- tmp=output(a,:); % get those rows
- tmp=tmp(:);
- tmp=unique(tmp);
- tmp(find(tmp==0))=[];%remove zeros
- output2(count,1:length(tmp))=tmp; %and combined set of grains to output2
- count=count+1;
-
-end
-
-%remove empty rows and columns from output2
-tmp = sum(output2,1);
-dummy=find(tmp==0);
-output2(:,dummy)=[];
-tmp = sum(output2,2);
-dummy=find(tmp==0);
-output2(dummy,:)=[];
-
-test_length=size(output2,1);
-if test_length==test_length_old %no change in the sets
- test=1;
-else
- test=0;
- test_length_old=test_length;
- output=output2;
- output2=zeros(length(output),50);
- %and do it all again
-end
-
-end % end of while loop
-
-output2=unique(output2,'rows');
-%output2 is now a table of grains that should be grouped together...
-
-
-
-
-
-%Next step - where output2 suggests combining grains, test the consistancy
-%of the result. If putting the grain together results in a mess, should
-%not do it - either because it's a twin, or because two large grains with
-%com close together.
-
-combines=output2(find(output2(:,2)~=0),:);%rows containing more than a single grain
-
-results=[];
-for i=1:size(combines,1)
-
- plane_normal=[];
- hkl=[];
- list=0; %dummy entry, delete after
-
- for j=1:size(combines,2)
- if combines(i,j)~=0
- tmp=load(sprintf('4_grains/grain%d_/grain%d_.mat',combines(i,j),combines(i,j)));
- % tmp=load(sprintf('grain%d_/grain%d_.mat',combines(i,j),combines(i,j)));
- plane_normal=[plane_normal; tmp.pl(find(tmp.index),:)]; %collect data for rodrigues
- hkl=[hkl; tmp.hkl(find(tmp.index),:)]; %note - consider only those being used for reconstruction
- list((end+1):(end+length(find(tmp.index))))=j;
- end
- end
-
- list(1)=[];%remove the dummy 0 entry
-
- [a,good_list]=plot_rodrigues_consistancy(plane_normal,hkl,spacegroup,0);
- out=list(find(good_list==1));% the spots which are consistant
-
- for j=1:size(combines,2) %calc the fraction of each grain that remains good after the combination
- results(i,j)=length(find(out==j))/length(find(list==j));
- end
-
-end
-%results contains the fractions of spots which stay good in each grain when
-%grains are combined as suggested. Only do the combinations which do not
-%mess up consistancy
-consistancy_test=0.85; %if consistancy drops below 90% of what it was
-results(find(isnan(results)))=1; %set NaN to 1 to ignore
-results=results<consistancy_test; %keep those above threshold
-
-combines(find(results))=0; %remove the grains that should not be combined
-%combines is now the list of grains to combines, but with some rows with
-%only one, or no grainids.
-kill_list=[];
-for i=1:size(combines,1)
- dummy=find(combines(i,:)~=0);%how many non zero entries?
- if length(dummy)<=1
- kill_list=[kill_list i];
- end
-end
-combines(kill_list,:)=[];%remove the rows with less than two grainids
-
-
-%tidy up the combines table - remove zeros, sort
-for i=1:size(combines,1)
- combines(i,:)=sort(combines(i,:),'descend');
-end
-combines=unique(combines,'rows');%get unique rows
-tmp = sum(combines,1);%remove entirely zero columns
-dummy=find(tmp==0);
-combines(:,dummy)=[];
-%final sort to ascending order but with zeros to the end!
-for i=1:size(combines,1)
- tmp=length(find(combines(i,:)~=0)); %number of non-zeros in row
- combines(i,1:tmp)=sort(combines(i,1:tmp));
-end
-
-combines=sortrows(combines);%final sort by rows
-%can have lines like this:
-% 5 8 17 23 24 25 64 0 0
-% 5 8 17 24 25 64 0 0 0
-% 5 8 17 25 64 0 0 0 0
-% 5 8 23 24 25 64 0 0 0
-%get the longest list, kill the rest
-
-kill_a=[];
-for i=1:max(combines(:))
- a=find(combines(:,1)==i);
- if length(a)>1 %more than one row
- count=zeros(size(a));%count the number of entries
- for j=1:length(a)
- count(j)=length(find(combines(a(j),:)~=0));
- end
-% keep_a = a(find(count==max(count)));
- kill_a = [kill_a ; a(find(count~=max(count)))];
- end
-end
-combines(kill_a, :)=[];
-
-
-
-
-
-
-
-
-
diff --git a/4_grains/combine_grains2_twins.m b/4_grains/combine_grains2_twins.m
deleted file mode 100644
index 7f8515338d2897cb465ad76ac25d063c8745beb1..0000000000000000000000000000000000000000
--- a/4_grains/combine_grains2_twins.m
+++ /dev/null
@@ -1,217 +0,0 @@
-function combines = combine_grains2_twins(r_vectors,positions, bads, spacegroup)
-% <<<< NOT USED >>>>
-%r_vectors = [grainid rx ry rz]
-%positions = [grainid grainx grainy grainz grainsize]
-%produced by gtReadGrains
-%combine_grains2 - need to improve the handling of what combines with
-%what...
-%apparently need to use positions calculated from the reconstructed,
-%post-processed grains - the x/y/zcenters are rubbish for this, need
-%x/y/zcenter2s
-%check consistancy because of the danger of combining twins, or other close
-%grains
-
-%get rid of bads grains before starting, also those already merged
-
-if 1
-dummy=find(bads==1 | bads==2); %consider possibly bad grains
-
-r_vectors(dummy,:)=[];
-positions(dummy,:)=[];
-else
- disp('not dealing with bad grains!')
-end
-
-output=zeros(size(r_vectors,1),50);%assume no more than 50 grains to be combined...
-
-for i=1:size(r_vectors,1)
-
- dif = r_vectors(:,2:4)-repmat(r_vectors(i,2:4),size(r_vectors,1),1);
- dif = dif.^2;
- r_distance = sqrt(sum(dif,2));
-
- dif = positions(:,2:4)-repmat(positions(i,2:4),size(r_vectors,1),1);
- dif = dif.^2;
- p_distance = sqrt(sum(dif,2));
- %consider the average size of the two grains
- p_distance = p_distance./(0.5*(positions(:,5)+repmat(positions(i,5),size(r_vectors,1),1)));
-
- %close in orientation and distance
- dummy =[];
- %dummy = [dummy; find(r_distance < 0.03 & p_distance < 0.2)];
- %if super close in orientation, can relax distance requirement to allow
- %for irregular shapes
- dummy = [dummy; find(r_distance < 0.03 & p_distance < 0.4)];
- %can have found the wrong r space solution - will test consistancy later
- %dummy = [dummy; find(p_distance < 0.2)];
-
- dummy = unique(dummy);
-
- output(i,1:(length(dummy))) = r_vectors(dummy,1)';
-
-end
-
-%crop output to only the required number of columns
-tmp = sum(output,1);
-dummy=find(tmp==0);
-output(:,dummy)=[];
-tmp = sum(output,2);
-dummy=find(tmp==0);
-output(dummy,:)=[];
-
-%get the unique rows - grains to be combined
-output = unique(output,'rows');
-
-
-%can still have multiple rows that refer to a single grain, eg ...;15 21
-%0;15 21 35;... in s5_dct5. Remove these...
-%want to keep the largest possible grain sets
-
-output2=zeros(length(output),50); %assume no more than 50 grains to be combined...
-count=1;
-test=0;
-test_length_old=size(output,1);
-
-while test==0
-
-for i=1:max(output(:)) % go though all grainids
-
-
- dummy=find(output==i);
- [a,b]=ind2sub(size(output),dummy); % a are the rows to combine
- tmp=output(a,:); % get those rows
- tmp=tmp(:);
- tmp=unique(tmp);
- tmp(find(tmp==0))=[];%remove zeros
- output2(count,1:length(tmp))=tmp; %and combined set of grains to output2
- count=count+1;
-
-end
-
-%remove empty rows and columns from output2
-tmp = sum(output2,1);
-dummy=find(tmp==0);
-output2(:,dummy)=[];
-tmp = sum(output2,2);
-dummy=find(tmp==0);
-output2(dummy,:)=[];
-
-test_length=size(output2,1);
-if test_length==test_length_old %no change in the sets
- test=1;
-else
- test=0;
- test_length_old=test_length;
- output=output2;
- output2=zeros(length(output),50);
- %and do it all again
-end
-
-end % end of while loop
-
-output2=unique(output2,'rows');
-%output2 is now a table of grains that should be grouped together...
-
-
-
-
-
-%Next step - where output2 suggests combining grains, test the consistancy
-%of the result. If putting the grain together results in a mess, should
-%not do it - either because it's a twin, or because two large grains with
-%com close together.
-
-combines=output2(find(output2(:,2)~=0),:);%rows containing more than a single grain
-
-results=[];
-for i=1:size(combines,1)
-
- plane_normal=[];
- hkl=[];
- list=0; %dummy entry, delete after
-
- for j=1:size(combines,2)
- if combines(i,j)~=0
- tmp=load(sprintf('4_grains/grain%d_/grain%d_.mat',combines(i,j),combines(i,j)), 'pl', 'hkl', 'index');
- % tmp=load(sprintf('grain%d_/grain%d_.mat',combines(i,j),combines(i,j)));
- tmp.index=tmp.index(1:length(tmp.pl));
- plane_normal=[plane_normal; tmp.pl(find(tmp.index),:)]; %collect data for rodrigues
- hkl=[hkl; tmp.hkl(find(tmp.index),:)]; %note - consider only those being used for reconstruction
- list((end+1):(end+length(find(tmp.index))))=j;
- end
- end
-
- list(1)=[];%remove the dummy 0 entry
-
- [a,good_list]=plot_rodrigues_consistancy(plane_normal,hkl,spacegroup,0);
- out=list(find(good_list==1));% the spots which are consistant
-
- for j=1:size(combines,2) %calc the fraction of each grain that remains good after the combination
- results(i,j)=length(find(out==j))/length(find(list==j));
- end
-
-end
-%results contains the fractions of spots which stay good in each grain when
-%grains are combined as suggested. Only do the combinations which do not
-%mess up consistancy
-consistancy_test=0.75; %if consistancy drops below 90% of what it was
-results(find(isnan(results)))=1; %set NaN to 1 to ignore
-results=results<consistancy_test; %keep those above threshold
-
-combines(find(results))=0; %remove the grains that should not be combined
-%combines is now the list of grains to combines, but with some rows with
-%only one, or no grainids.
-kill_list=[];
-for i=1:size(combines,1)
- dummy=find(combines(i,:)~=0);%how many non zero entries?
- if length(dummy)<=1
- kill_list=[kill_list i];
- end
-end
-combines(kill_list,:)=[];%remove the rows with less than two grainids
-
-
-%tidy up the combines table - remove zeros, sort
-for i=1:size(combines,1)
- combines(i,:)=sort(combines(i,:),'descend');
-end
-combines=unique(combines,'rows');%get unique rows
-tmp = sum(combines,1);%remove entirely zero columns
-dummy=find(tmp==0);
-combines(:,dummy)=[];
-%final sort to ascending order but with zeros to the end!
-for i=1:size(combines,1)
- tmp=length(find(combines(i,:)~=0)); %number of non-zeros in row
- combines(i,1:tmp)=sort(combines(i,1:tmp));
-end
-
-combines=sortrows(combines);%final sort by rows
-%can have lines like this:
-% 5 8 17 23 24 25 64 0 0
-% 5 8 17 24 25 64 0 0 0
-% 5 8 17 25 64 0 0 0 0
-% 5 8 23 24 25 64 0 0 0
-%get the longest list, kill the rest
-
-kill_a=[];
-for i=1:max(combines(:))
- a=find(combines(:,1)==i);
- if length(a)>1 %more than one row
- count=zeros(size(a));%count the number of entries
- for j=1:length(a)
- count(j)=length(find(combines(a(j),:)~=0));
- end
-% keep_a = a(find(count==max(count)));
- kill_a = [kill_a ; a(find(count~=max(count)))];
- end
-end
-combines(kill_a, :)=[];
-
-
-
-
-
-
-
-
-
diff --git a/4_grains/gtCalcOmega.m b/4_grains/gtCalcOmega.m
deleted file mode 100644
index 09a6c97e1321fdd48d89e798932f6931f455fbc6..0000000000000000000000000000000000000000
--- a/4_grains/gtCalcOmega.m
+++ /dev/null
@@ -1,53 +0,0 @@
-function [Theta,Eta,Omega]=gtCalcOmega(normal, hkl)
-
-%calculates the Omega position at which diffraction will occur from a plane
-%defined by its plane normal and hkl type [h k l]
-
-%load data
-parameters=[];
-if isempty(parameters)
- load('parameters.mat');
-end
-
-%put plane normal in cylindrical polar coords
-length = sqrt(normal(1)^2 + normal(2)^2);
-alpha = atan2(normal(2), normal(1));
-vert = normal(3);
-
-%calc Theta for diffraction
-lambda = 12.398/parameters.acq.energy;
-if parameters.cryst.spacegroup == 225
-d=parameters.cryst.latticepar(1)/sqrt(hkl*hkl');
-else
- disp('only works for FCC at the mo...')
-end
-Theta = asin(lambda/(2*d));
-
-try
- %produce Omega angles
- temp(1) = acos((cos((pi/2)-Theta))/length);
- temp(2) = pi-temp(1);
- Omega = temp - alpha;
- Omega(3:4) = Omega(1:2)+pi;
-
- %put all in 0-2pi interval
- dummy = find(Omega<0);
- Omega(dummy) = Omega(dummy) + (2*pi);
- dummy = find(Omega>(2*pi));
- Omega(dummy) = Omega(dummy) - (2*pi);
-
- %calc Eta
-for i=1:size(Omega,2)
- alpha2 = alpha + Omega(i);
- y = length*sin(alpha2);
- Eta(i) = sin(y/vert);
- Theta(i) = Theta(1);
-end
-
-catch%if no solutions
- disp('no diffraction from this plane normal!')
-end
-
-
-
-
diff --git a/4_grains/gtCentroidOfRegion.m b/4_grains/gtCentroidOfRegion.m
deleted file mode 100644
index 8f44ac82ebf9b3c8bdd6103abd32eda458e30dd5..0000000000000000000000000000000000000000
--- a/4_grains/gtCentroidOfRegion.m
+++ /dev/null
@@ -1,10 +0,0 @@
-function c=gtCentroidOfRegion(m)
-% Gives the centroid x and y - first moment - of a matrix.
-% Minimum value: [1 1]
-% Maximum value: [size(m,1) size(m,2)]
-
-all=sum(sum(m));
-c(1)=sum(sum(m,2).*[1:size(m,1)]')/all;
-c(2)=sum(sum(m,1).*[1:size(m,2)])/all;
-
-end
diff --git a/4_grains/gtDoCombineGrains.m b/4_grains/gtDoCombineGrains.m
deleted file mode 100644
index 950d146992907d06818c8708ab65476c593cb1b8..0000000000000000000000000000000000000000
--- a/4_grains/gtDoCombineGrains.m
+++ /dev/null
@@ -1,167 +0,0 @@
-function gtDoCombineGrains(list)
-% <<<< NOT USED >>>>
-%list, generated by combine_grains.m, is a list of which grains in the
-%4_grains/ directory should be merged [g1 g2 0; g3 0 0; g4 g9 0...]
-%create the "new" grains in a new directory for safety
-
-parameters=[];
-load('parameters.mat');
-
-
-
-%warning('starting combine grains from grain 100')
-for i=1:size(list,1) %each of the new grains
-
-old_grain_number = list(i,1);
-new_grain_number = i;%+ 150%to distinguish from centre slice grains
-%make the new directory
-mkdir(sprintf('4_grains/new_grains/grain%d_',new_grain_number))
-
-
-dummy=find(list(i,2:end)~=0);
-if isempty(dummy) %if no combining of grains, just a renumbering...
-
- tmp = load(sprintf('4_grains/grain%d_/grain%d_.mat',list(i,1),list(i,1)));%load the old grain
- struct_ids=tmp.struct_ids;%need this for database update
- pair_vector=tmp.pair_vector;%need this for database update
-
- %copy the old thresholded volume file (grain%d_2.edf) to the new place
- try
- cmd=sprintf('cp 4_grains/grain%d_/grain%d_2.edf 4_grains/new_grains/grain%d_/grain%d_2.edf',...
- old_grain_number,old_grain_number,new_grain_number,new_grain_number);
- unix(cmd)
- end
- try
- %copy the old mat file to the new place
- cmd=sprintf('cp 4_grains/grain%d_/grain%d_.mat 4_grains/new_grains/grain%d_/grain%d_.mat',...
- old_grain_number,old_grain_number,new_grain_number,new_grain_number);
- unix(cmd)
- end
-
-else%if two or more grains to be combined
-
- %initialise variables fo the new grain
- zstart=[];zend=[];
- nx=[];x1=[];
- ny=[];y1=[];
- plane_normal=[];hkl=[];R_vector=[];
- stack=[];
- Omega=[];Theta=[];Eta=[];
- struct_ids=[];
- pair_vector=[];
- index=[];replaced=[];lambda=[];
- %some for calculations
- maxx=[];maxy=[];
-
- for j=1:size(list,2) %max number of old grains to combine
- if list(i,j)~=0 %ignoring zeros in list
- tmp = load(sprintf('4_grains/grain%d_/grain%d_.mat',list(i,j),list(i,j)));%load the old grain
- zstart = [zstart tmp.zstart];
- zend = [zend tmp.zend];
- x1 = [x1 tmp.x1];
- y1 = [y1 tmp.y1];
- maxx = [maxx (tmp.x1+tmp.nx)];
- maxy = [maxy (tmp.y1+tmp.ny)];
- stack = cat(3,stack,tmp.stack);
- plane_normal = [plane_normal; tmp.plane_normal];
- R_vector = [R_vector; tmp.R_vector];
- lambda = [lambda; tmp.lambda];
- hkl = [hkl; tmp.hkl];
- Omega = [Omega tmp.Omega];
- Theta = [Theta tmp.Theta];
- Eta = [Eta tmp.Eta];
- struct_ids = [struct_ids tmp.struct_ids];
-
- if isfield(tmp,'pair_vector')
- pair_vector = [pair_vector tmp.pair_vector];
- else
- addpair_vector = logical(ones(1,length(tmp.struct_ids)));
- pair_vector = [pair_vector tmp.pair_vector];
- end
-
- if isfield(tmp,'index') && length(tmp.index)==length(tmp.struct_ids)%some have wrong length index!
- addindex = reshape(tmp.index,1,length(tmp.struct_ids));%make into a row
- index = [index addindex];
- else
- addindex = logical(ones(1,length(tmp.struct_ids)));
- index = [index addindex];
- end
-
- if isfield(tmp,'replaced') && length(tmp.replaced)==length(tmp.struct_ids)
- addreplaced = reshape(tmp.replaced,1,length(tmp.struct_ids));%make into a row
- replaced = [replaced addreplaced];
- else
- addreplaced = logical(zeros(1,length(tmp.struct_ids)));
- replaced = [replaced addreplaced];
- end
-
- end
- end
-
-
- %sort where needed
- zstart = min(zstart);
- zend = max(zend);
- zcenter = ceil((zstart+zend)/2);
- x1 = min(x1);
- y1 = min(y1);
- maxx = max(maxx);
- maxy = max(maxy);
- nx = maxx-x1;
- ny = maxy-y1;
- xcenter = ceil((x1+nx)/2);
- ycenter = ceil((y1+ny)/2);
- [Omega,tmp]=sort(Omega);
- Theta = Theta(tmp);
- Eta = Eta(tmp);
- struct_ids = struct_ids(tmp);
- pair_vector=pair_vector(tmp);
- plane_normal = plane_normal(tmp,:);
- R_vector=mean(R_vector,1);
- lambda = mean(lambda,1);
- hkl = hkl(tmp,:);
- stack = stack(:,:,tmp);
- index = index(tmp);
- replaced = replaced(tmp);
-
- %may need also to remove duplicates....
- [struct_ids,tmp]=unique(struct_ids);
- pair_vector=pair_vector(tmp);
- Omega=Omega(tmp);
- Theta=Theta(tmp);
- Eta=Eta(tmp);
- index=index(tmp);
- replaced=replaced(tmp);
- plane_normal=plane_normal(tmp,:);
- hkl=hkl(tmp,:);
- stack=stack(:,:,tmp);
-
-
-
- %write the new grain%d_.mat
- name = sprintf('4_grains/new_grains/grain%d_/grain%d_.mat',new_grain_number,new_grain_number);
- save(name,'struct_ids','pair_vector','stack','Theta','Eta','Omega','plane_normal','hkl','R_vector',...
- 'zstart','zend','zcenter','x1','nx','xcenter','y1','ny','ycenter','index','replaced','lambda');
-
-
-end
-
- %update database with new grainids
- disp('update GrainID field in database')
- for k=1:length(struct_ids)
-
- if pair_vector(k)==2
- db_name=parameters.acq.pair_name;
- else
- db_name=parameters.acq.name;
- end
-
- mym(sprintf('update %sextspot set GrainID=%d where extspotID=%d',...
- db_name,new_grain_number,struct_ids(k)))
- mym(sprintf('update %sbb set GrainID=%d where extspotID=%d',...
- db_name,new_grain_number,struct_ids(k)))
- end
-
-end
-
-disp('Done!')
diff --git a/4_grains/gtDoCombineGrains2.m b/4_grains/gtDoCombineGrains2.m
deleted file mode 100644
index 01d0254e10d607e3b6574e511bf4b00ac427890f..0000000000000000000000000000000000000000
--- a/4_grains/gtDoCombineGrains2.m
+++ /dev/null
@@ -1,168 +0,0 @@
-function gtDoCombineGrains2(list, phaseid)
-% <<<< NOT USED >>>>
-%list, generated by combine_grains2.m, is a list of which grains in the
-%4_grains/ directory should be merged [g1 g2 0; g3 0 0; g4 g9 0...]
-
-%rather than moving everything around, rather just add data to the .mat
-%file of the first grain, and set the bad value of the others to 2. Make
-%sure that the bad value of this first grain is 0 (could have been 0.5
-%before)
-
-%avoid using this old version by mistake!
-check=inputwdefault('Should use version _360!!! continue with this old version? [y/n]', 'n');
-if check=='n'
- return;
-end
-
-load('parameters.mat');
-
-gtDBConnect();
-
-for i=1:size(list,1) %each of the new grains
-
- %initialise variables fo the new grain
- zstart=[];zend=[];
- nx=[];x1=[];
- ny=[];y1=[];
- plane_normal=[];hkl=[];R_vector=[];
- stack=[];
- Omega=[];Theta=[];Eta=[];
- struct_ids=[];
- pair_vector=[];
- index=[];replaced=[];lambda=[];
- %some for calculations
- maxx=[];maxy=[];
-
- for j=1:size(list,2) %max number of old grains to combine
- if list(i,j)~=0 %ignoring zeros in list
-
- tmp = load(sprintf('4_grains/grain%d_/grain%d_.mat',list(i,j),list(i,j)));%load the old grain
- zstart = [zstart tmp.zstart];
- zend = [zend tmp.zend];
- x1 = [x1 tmp.x1];
- y1 = [y1 tmp.y1];
- maxx = [maxx (tmp.x1+tmp.nx)];
- maxy = [maxy (tmp.y1+tmp.ny)];
- stack = cat(3,stack,tmp.stack);
- plane_normal = [plane_normal; tmp.plane_normal];
- R_vector = [R_vector; tmp.R_vector]; %in fact will recalculate this!
- lambda = [lambda; tmp.lambda];
- hkl = [hkl; tmp.hkl];
- Omega = [Omega tmp.Omega];
- Theta = [Theta tmp.Theta];
- Eta = [Eta tmp.Eta];
- struct_ids = [struct_ids tmp.struct_ids];
-
- if isfield(tmp,'pair_vector')
- pair_vector = [pair_vector tmp.pair_vector];
- else
- addpair_vector = logical(ones(1,length(tmp.struct_ids)));
- pair_vector = [pair_vector addpair_vector]%tmp.pair_vector];
- end
-
- if isfield(tmp,'index') && length(tmp.index)==length(tmp.struct_ids)%some have wrong length index!
- addindex = reshape(tmp.index,1,length(tmp.struct_ids));%make into a row
- index = [index addindex];
- else
- addindex = logical(ones(1,length(tmp.struct_ids)));
- index = [index addindex];
- end
-
- if isfield(tmp,'replaced') && length(tmp.replaced)==length(tmp.struct_ids)
- addreplaced = reshape(tmp.replaced,1,length(tmp.struct_ids));%make into a row
- replaced = [replaced addreplaced];
- else
- addreplaced = logical(zeros(1,length(tmp.struct_ids)));
- replaced = [replaced addreplaced];
- end
-
- end %ignoring zeros
- end %for each grain to be combined
-
-
- %sort where needed
- zstart = min(zstart);
- zend = max(zend);
- zcenter = ceil((zstart+zend)/2);
- x1 = min(x1);
- y1 = min(y1);
- maxx = max(maxx);
- maxy = max(maxy);
- nx = maxx-x1;
- ny = maxy-y1;
- xcenter = ceil((x1+nx)/2);
- ycenter = ceil((y1+ny)/2);
- [Omega,tmp]=sort(Omega);
- Theta = Theta(tmp);
- Eta = Eta(tmp);
- struct_ids = struct_ids(tmp);
- pair_vector=pair_vector(tmp);
- plane_normal = plane_normal(tmp,:);
- R_vector=mean(R_vector,1);
- lambda = mean(lambda,1);
- hkl = hkl(tmp,:);
- stack = stack(:,:,tmp);
- index = index(tmp);
- replaced = replaced(tmp);
-
- %may need also to remove duplicates....
- [struct_ids,tmp]=unique(struct_ids);
- pair_vector=pair_vector(tmp);
- Omega=Omega(tmp);
- Theta=Theta(tmp);
- Eta=Eta(tmp);
- index=index(tmp);
- replaced=replaced(tmp);
- plane_normal=plane_normal(tmp,:);
- hkl=hkl(tmp,:);
- stack=stack(:,:,tmp);
-
- %now have collected the new data, ready to write to the new .mat
-
- %first check orientation / consistancy
- [R_vector, good_list]=plot_rodrigues_consistancy(plane_normal(find(index),:),hkl(find(index),:),parameters.cryst(phaseid).spacegroup,0);
-
-
- %set other bad flags to 2
- for j=2:size(list,2)
- if list(i,j)~=0 %ignoring zeros in list
- bad=2;
- save(sprintf('4_grains/grain%d_/grain%d_.mat',list(i,j),list(i,j)),'bad','-append');
- end
- end
-
- %bad flag for the new grain = 0 (good)
- bad=0;
-
-
-
-
- %save combined info into the first grain mat file
- name = sprintf('4_grains/grain%d_/grain%d_.mat',list(i,1),list(i,1));
- save(name,'struct_ids','pair_vector','stack','Theta','Eta','Omega','plane_normal','hkl','R_vector',...
- 'zstart','zend','zcenter','x1','nx','xcenter','y1','ny','ycenter','index','replaced','lambda','bad','-append');
-
-
-
- %update database with new grainids
- disp('update GrainID field in database')
- for k=1:length(struct_ids)
-
- if pair_vector(k)==2
- db_name=parameters.acq.pair_name;
- else
- db_name=parameters.acq.name;
- end
-
- mym(sprintf('update %sextspot set GrainID=%d where extspotID=%d',...
- db_name,list(i,1),struct_ids(k)))
-
- end
-
- %write the stack (replaced) for the combined grain
- gtWriteStack_replaced(list(i,1));
-
-
-end
-
-disp('Done!')
diff --git a/4_grains/gtDoCombineGrains_360.m b/4_grains/gtDoCombineGrains_360.m
deleted file mode 100644
index a4a59b189427c612f78343bc1c1bbed95d288273..0000000000000000000000000000000000000000
--- a/4_grains/gtDoCombineGrains_360.m
+++ /dev/null
@@ -1,266 +0,0 @@
-function gtDoCombineGrains_360(list, phaseid)
-% <<<< NOT USED >>>>
-%list, generated by combine_grains2.m, is a list of which grains in the
-%4_grains/ directory should be merged [g1 g2 0; g3 0 0; g4 g9 0...]
-
-%rather than moving everything around, rather just add data to the .mat
-%file of the first grain, and set the bad value of the others to 2. Make
-%sure that the bad value of this first grain is 0 (could have been 0.5
-%before)
-
-%360 degree scan concept, no pair vectors
-
-load('parameters.mat');
-
-gtDBConnect()
-
-for i=1:size(list,1) %each of the new grains
-
- %initialise variables fo the new grain
- zstart=[];zend=[];
- nx=[];x1=[];
- ny=[];y1=[];
- xcenter=[]; ycenter=[]; zcenter=[];
- plane_normal=[];hkl=[];R_vector=[];
- stack=[];
- Omega=[];Theta=[];Eta=[];
- struct_ids=[];
- index=[];replaced=[];lambda=[];
- %some for calculations
- maxx=[];maxy=[];
-
- for j=1:size(list,2) %max number of old grains to combine
- if list(i,j)~=0 %ignoring zeros in list
- % modif sabine to keep original mat file
- cmd = sprintf('cp 4_grains/grain%d_/grain%d_.mat 4_grains/grain%d_/grain%d_.mat_original', list(i,j), list(i,j), list(i,j), list(i,j));
- % cmd = sprintf('cp grain%d_/grain%d_.mat grain%d_/grain%d_.mat_original', list(i,j), list(i,j), list(i,j), list(i,j));
-
- unix(cmd);
- % fin modif sabine
-
- tmp = load(sprintf('4_grains/grain%d_/grain%d_.mat',list(i,j),list(i,j)));%load the old grain
- % tmp = load(sprintf('grain%d_/grain%d_.mat',list(i,j),list(i,j)));%load the old grain
- zstart = [zstart tmp.zstart];
- zend = [zend tmp.zend];
- x1 = [x1 tmp.x1];
- y1 = [y1 tmp.y1];
- xcenter = [xcenter tmp.xcenter];
- ycenter = [ycenter tmp.ycenter];
- zcenter = [zcenter tmp.zcenter];
- maxx = [maxx (tmp.x1+tmp.nx)];
- maxy = [maxy (tmp.y1+tmp.ny)];
- stack = cat(3,stack,tmp.stack);
- plane_normal = [plane_normal; tmp.plane_normal];
- R_vector = [R_vector; tmp.R_vector]; %in fact will recalculate this!
- lambda = [lambda; tmp.lambda];
- hkl = [hkl; tmp.hkl];
- Omega = [Omega tmp.Omega];
- Theta = [Theta tmp.Theta];
- Eta = [Eta tmp.Eta];
- struct_ids = [struct_ids tmp.struct_ids];
-
-
- if isfield(tmp,'index') && length(tmp.index)==length(tmp.struct_ids)%some have wrong length index!
- addindex = reshape(tmp.index,1,length(tmp.struct_ids));%make into a row
- index = [index addindex];
- else
- addindex = logical(ones(1,length(tmp.struct_ids)));
- index = [index addindex];
- end
-
- if isfield(tmp,'replaced') && length(tmp.replaced)==length(tmp.struct_ids)
- addreplaced = reshape(tmp.replaced,1,length(tmp.struct_ids));%make into a row
- replaced = [replaced addreplaced];
- else
- addreplaced = logical(zeros(1,length(tmp.struct_ids)));
- replaced = [replaced addreplaced];
- end
-
- end %ignoring zeros
- end %for each grain to be combined
-
-
- %sort where needed
- zstart = min(zstart);
- zend = max(zend);
- %zcenter = ceil((zstart+zend)/2);
- zcenter=ceil(mean(zcenter));
- x1 = min(x1);
- y1 = min(y1);
- maxx = max(maxx);
- maxy = max(maxy);
- nx = maxx-x1;
- ny = maxy-y1;
-% xcenter = ceil((x1+nx)/2);
-% ycenter = ceil((y1+ny)/2);
- xcenter=ceil(mean(xcenter));
- ycenter=ceil(mean(ycenter));
- [Omega,tmp]=sort(Omega);
- Theta = Theta(tmp);
- Eta = Eta(tmp);
- struct_ids = struct_ids(tmp);
- plane_normal = plane_normal(tmp,:);
- R_vector=mean(R_vector,1);
-% lambda = mean(lambda,1); % modif sabine
- hkl = hkl(tmp,:);
- stack = stack(:,:,tmp);
- index = index(tmp);
- replaced = replaced(tmp);
-
- %may need also to remove duplicates....
- [struct_ids,tmp]=unique(struct_ids);
- Omega=Omega(tmp);
- Theta=Theta(tmp);
- Eta=Eta(tmp);
- index=index(tmp);
- replaced=replaced(tmp);
- plane_normal=plane_normal(tmp,:);
- hkl=hkl(tmp,:);
- stack=stack(:,:,tmp);
-
- % modif sabine
- if length(struct_ids)>40
- lambda = [0.1 0.05 0.02];
- elseif length(struct_ids)>30
- %lambda = [0.3 0.2 0.1]; %initial values
- lambda = [0.2 0.1 0.05]; %initial values - sheared difspots
- elseif length(struct_ids)>15
- %lambda = [0.5 0.3 0.1]; %initial values
- lambda = [0.3 0.2 0.1]; %initial values - sheared difspots
- else
- %lambda = [0.7 0.3 0.2]; %initial values
- lambda = [0.4 0.25 0.15]; %initial values - sheared difspots
- end
-
- % fin modif sabine
-
- if 1 %sheared difspots
- acq=parameters.acq;
-
- %don't know why this isn't work. use an average center for the moment
-% %If using sheared difspots, need to recalulate the grain centre, and then
-% %rebuild the difstack
-% %centre:
-% lines=[];% [point, direction]
-% %get all pair information available
-% for k=1:length(struct_ids)
-% mysqlcmd=sprintf('select samcentXA,samcentYA,samcentZA,ldirZ,ldirY,ldirZ from %s where difAID=%d or difBID=%d',...
-% acq.pair_tablename, struct_ids(k), struct_ids(k));
-% [a,b,c,d,e,f]=mym(mysqlcmd);
-% if ~isempty(a)
-% lines(end+1,:)=[a b c d e f];
-% end
-% end
-% keyboard
-% lines=unique(lines, 'rows');
-% grain.center=pofintersectexp(lines)'; % new center of mass
-% grain.xcenter= (acq.bb(3)/2)+grain.center(2);
-% grain.ycenter= (acq.bb(3)/2)+grain.center(1);
-% grain.zcenter= -grain.center(3)+acq.ydet/2-acq.bb(2);
-%
-
-
-
- %add other fields to grain to pass to dtShearDifspots
- grain.struct_ids=struct_ids;
- grain.theta=Theta;
- grain.eta=Eta;
- grain.xcenter=xcenter;
- grain.ycenter=ycenter;
- grain.zcenter=zcenter;
- %get Omega from difspot table to ensure no confusion between real/sheared
- %values
- for k=1:length(struct_ids)
- Omega(k)=mym(sprintf('select CentroidImage from %sdifspot where difspotid=%d', acq.name,struct_ids(k)));
- end
- Omega = Omega/(acq.nproj*2)*360;
- grain.omega=Omega;
- [stack, difomega]=gtShearDifspots(grain);
- %finally save the sheared values to the .mat file
- Omega=difomega;
- end
-
-
- %now have collected the new data, ready to write to the new .mat
-
- %first check orientation / consistancy
- [R_vector, good_list]=plot_rodrigues_consistancy(plane_normal(find(index),:),hkl(find(index),:),parameters.cryst(phaseid).spacegroup,0);
-
-
- %set other bad flags to 2
- for j=2:size(list,2)
- if list(i,j)~=0 %ignoring zeros in list
- bad=2;
- save(sprintf('4_grains/grain%d_/grain%d_.mat',list(i,j),list(i,j)),'bad','-append');
- end
- end
-
- %bad flag for the new grain = 0 (good)
- bad=0;
-
-
-
-
- %save combined info into the first grain mat file
-
- name = sprintf('4_grains/grain%d_/grain%d_.mat',list(i,1),list(i,1));
-% name = sprintf('grain%d_/grain%d_.mat',list(i,1),list(i,1));
-
-% modif sabine
-if length(struct_ids)>40
- lambda = [0.1 0.05 0.02];
- elseif length(struct_ids)>30
- %lambda = [0.3 0.2 0.1]; %initial values
- lambda = [0.2 0.1 0.05]; %initial values - sheared difspots
- elseif length(struct_ids)>15
- %lambda = [0.5 0.3 0.1]; %initial values
- lambda = [0.3 0.2 0.1]; %initial values - sheared difspots
- else
- %lambda = [0.7 0.3 0.2]; %initial values
- lambda = [0.4 0.25 0.15]; %initial values - sheared difspots
-end
-
-
-% fin modif sabine
-
-
-
- save(name,'struct_ids','stack','Theta','Eta','Omega','plane_normal','hkl','R_vector',...
- 'zstart','zend','zcenter','x1','nx','xcenter','y1','ny','ycenter','index','replaced','lambda','bad','-append');
-
-
- %update database with new grainids
- disp('update GrainID field in database')
- extspot_failed=0;
- for k=1:length(struct_ids)
-
- try
- mym(sprintf('update %sextspot set GrainID=%d where extspotID=%d',...
- parameters.acq.name,list(i,1),struct_ids(k)))
- catch
- extspot_failed=1;
- end
-
- mym(sprintf('update %sdifspot set GrainID=%d where difspotID=%d',...
- parameters.acq.name,list(i,1),struct_ids(k)))
- end
- if extspot_failed
- disp('no extspot table found for this sample - not updating!')
- end
-
- %write the stack (replaced) for the combined grain
-
- gtWriteStack_360(list(i,1));
-
-
- %repeat the autoselecting of which projections to use
- try
- %not many licences, and not critical
- gt_select_projections_auto(list(i,1), 0); %0 is plotflag
- end
- %rerun ART
- gtDoART(list(i,1));
-
-end
-
-disp('Done!')
diff --git a/4_grains/gtFindAngles.m b/4_grains/gtFindAngles.m
deleted file mode 100644
index 14d7265854213e6d042373e3994a3662c8191e9d..0000000000000000000000000000000000000000
--- a/4_grains/gtFindAngles.m
+++ /dev/null
@@ -1,44 +0,0 @@
-
-
-function [Theta,Eta] = gtFindAngles(struct_id, grainx, grainy, grainz, omega);
-parameters=[];
-%calculate Theta angle between grain position (x,y,z) and difspot of
-%struct_id
-%struct_id can be a vector...
-%omega should be input in degrees
-%output angles are in degrees
-
-omega=omega*pi/180;%convert to Omega
-
-for i=1:length(struct_id)
-
-%mysqlcmd=sprintf('select Xcentroid from %sbboxes inner join %sbb on bbID=bboxID where %sbb.extspotID=%d',...
-% parameters.acq.name,parameters.acq.name,parameters.acq.name,struct_id(i))
-%xdirect(i)=parameters.acq.bb(1)+mym(mysqlcmd);
-
-[bb,cent]=gtGetBBProps(struct_id(i));
-CentroidX=cent(1);
-xdirect(i)=parameters.acq.bb(1) + CentroidX;
-ydirect(i)=parameters.acq.bb(2)+grainz;
-
-mysqlcmd=sprintf('select CentroidX,CentroidY from %sdifspot where difspotID=%d',parameters.acq.name,struct_id(i));
-[xdiffr(i),ydiffr(i)]=mym(mysqlcmd);
-
-end
-
-
-r=sqrt((xdirect-xdiffr).^2+(ydirect-ydiffr).^2);
-
-grain2rot=sqrt((grainx-parameters.acq.bb(3)/2).^2+(grainy-parameters.acq.bb(3)/2).^2);
-angle=atan2(parameters.acq.bb(3)/2-grainy,grainx-parameters.acq.bb(3)/2);
-
-dz=grain2rot.*sin(-omega+angle)+parameters.acq.dist/parameters.acq.pixelsize;
-
-Theta=((atan(r./dz))/2)*180/pi;
-
-Eta=atan2((xdiffr-xdirect),(ydirect-ydiffr))*180/pi;
-
-dummy= find(Eta<0);
- Eta(dummy)=360+Eta(dummy);
-
-end
diff --git a/4_grains/gtFindAngles_pair.m b/4_grains/gtFindAngles_pair.m
deleted file mode 100644
index a48f68693c2a02ec12fdee5df0274dbd324dff7f..0000000000000000000000000000000000000000
--- a/4_grains/gtFindAngles_pair.m
+++ /dev/null
@@ -1,60 +0,0 @@
-
-
-function [Theta,Eta] = gtFindAngles_pair(struct_id, pair, grainx, grainy, grainz, omega, parameters);
-
-
-
-
-if isempty(parameters)
- load('parameters.mat');
-end
-%calculate Theta angle between grain position (x,y,z) and difspot of
-%struct_id
-%struct_id can be a vector...
-%omega should be input in degrees
-%output angles are in degrees
-%pair - 1 for this dataset, 2 if struct_id refers to the pair
-
-omega=omega*pi/180;%convert to radians
-
-for i=1:length(struct_id)
-
-%mysqlcmd=sprintf('select Xcentroid from %sbboxes inner join %sbb on bbID=bboxID where %sbb.extspotID=%d',...
-% parameters.acq.name,parameters.acq.name,parameters.acq.name,struct_id(i))
-%xdirect(i)=parameters.acq.bb(1)+mym(mysqlcmd);
-
-[bb,cent]=gtGetBBProps_pair(struct_id(i), pair(i));
-CentroidX=cent(1);
-xdirect(i)=parameters.acq.bb(1) + CentroidX;
-ydirect(i)=parameters.acq.bb(2)+grainz;
-
-if pair(i) == 2
- name = parameters.acq.pair_name;
-elseif pair(i) == 1
-name=parameters.acq.name;
-else
- error('need to specify which half of scan (pair=1 or 2)')
- return
-end
-
-mysqlcmd=sprintf('select CentroidX,CentroidY from %sdifspot where difspotID=%d',name,struct_id(i));
-[xdiffr(i),ydiffr(i)]=mym(mysqlcmd);
-
-end
-
-
-r=sqrt((xdirect-xdiffr).^2+(ydirect-ydiffr).^2);
-
-grain2rot=sqrt((grainx-parameters.acq.bb(3)/2).^2+(grainy-parameters.acq.bb(3)/2).^2);
-angle=atan2(parameters.acq.bb(3)/2-grainy,grainx-parameters.acq.bb(3)/2);
-
-dz=grain2rot.*sin(-omega+angle)+parameters.acq.dist/parameters.acq.pixelsize;
-
-Theta=((atan(r./dz))/2)*180/pi;
-
-Eta=atan2((xdiffr-xdirect),(ydirect-ydiffr))*180/pi;
-
-dummy= find(Eta<0);
- Eta(dummy)=360+Eta(dummy);
-
-end
diff --git a/4_grains/gtFindGroupInCand.m b/4_grains/gtFindGroupInCand.m
deleted file mode 100644
index 309e476970f3f73e2dd93b1c492dff34e6252505..0000000000000000000000000000000000000000
--- a/4_grains/gtFindGroupInCand.m
+++ /dev/null
@@ -1,71 +0,0 @@
-
-function [pgood,grainoutput]=gtFindGroupInCand(cand,tol,sample,ACM,spacegroup)
-% cand: all the lines that possibly belong to the same grain as the 1st line
-% (they are pair-consistent)
-% 1st line in "cand" is the actual line investigated
-
-% NOTE: p2 and p3 intersection might be far from the actual grain center if
-% they are near being parallel
-
-nof_cand=length(cand.id);
-
-p1=1;
-
-for p2=2:nof_cand % loop to find a 2nd line of the grain
- for p3=p2+1:nof_cand % loop to find a 3rd line of the grain
-
- pgood=false(nof_cand,1);
- pgood([p1,p2])=true;
-
- % Assume p2 belong to p1, try a 3rd one:
- % is p1,p2,p3 group-consistent?
-
- ok=gtTryAddLineToGroup(p3,pgood,cand,tol,sample,ACM);
-
- if ok
- pgood(p3)=true; % p1,p2,p3 are group consistent
-
- for p4=p3+1:nof_cand
- % having a group of p1,p2,p3 does p4 fit as well?
- ok=gtTryAddLineToGroup(p4,pgood,cand,tol,sample,ACM);
-
- if ok
- % p1,p2,p3,p4 are group-consistent.
- % We accept them as a good solution.
- % Which else could belong to this group?
- % Loop through all the remaining candidates and try whether they
- % fit. Keep the result, that is all the lines from this grain.
- pgood(p4)=true;
-
- for pn=p4+1:nof_cand
- ok=gtTryAddLineToGroup(pn,pgood,cand,tol,sample,ACM);
- if ok
- pgood(pn)=true;
- end
- end
-
- grainl=gtKeepLine(pgood,cand); % lines in the grain
-
- % Grain info:
- % modif sabine to run erik sample
- grainoutput=gtIndexCreateGrainOutputBuild(grainl); %,spacegroup);
-
- % grainoutput=gtIndexCreateGrainOutput(grainl,ACM);
- % end modif sabine
-
- return
-
- end
-
- end
- end
-
- end
-end
-
-% If no grain is found:
-pgood=false(nof_cand,1);
-pgood(1)=true;
-grainoutput=[];
-
-end
diff --git a/4_grains/gtGetBBProps_pair.m b/4_grains/gtGetBBProps_pair.m
deleted file mode 100644
index 5949b19deb98e95c54d7218ed9600590063fe659..0000000000000000000000000000000000000000
--- a/4_grains/gtGetBBProps_pair.m
+++ /dev/null
@@ -1,57 +0,0 @@
-function [bb,centroid]=gtGetBBProps_pair(struct_id, pair, bbtype)
-% usage: [bb,centroid]=gtGetBBProps(struct_id); % to get best bbox available
-% or [bb,centroid]=gtGetBBProps(struct_id,'search'); % to get search bbox
-% pair: 1 for this scan, 2 if the struct_id refers to the pair dataset
-parameters=[];
-
-if isempty(parameters)
- disp('Loading parameter file')
- load('parameters.mat');
-end
-
-
-if pair == 2
- name_root = parameters.acq.pair_name;
-elseif pair == 1
-name_root=parameters.acq.name;
-else
- Error('Need to specify "pair" as 1 or 2 - ak 1/2007')
-return
-end
-
-if exist('bbtype','var')
- % if specified, use the searchBB explicitly.
- disp('Using search bounding box')
-
- mysqlcmd=sprintf(...
- 'select searchbbID from %sextspot where extspotID=%d',...
- name_root,struct_id);
- searchbboxID=mym(mysqlcmd);
-
- mysqlcmd=sprintf(...
- 'select Xorigin, Yorigin, Xsize, Ysize,Xcentroid,Ycentroid from %sbboxes where bboxID=%d',...
- name_root,...
- searchbboxID);
-
-
-else
- % otherwise, use the bbview to return the best bounding box (and
- % centroid)
-
- mysqlcmd=sprintf(...
- 'select Xorigin,Yorigin,Xsize,Ysize,Xcentroid,Ycentroid from %sbb inner join %sbboxes on %sbb.bbID=%sbboxes.bboxID where %sbb.extspotID = %d', ...
- name_root,...
- name_root,...
- name_root,...
- name_root,...
- name_root,...
- struct_id);
-end
-
-[bb(1),bb(2),bb(3),bb(4),centroid(1),centroid(2)]=mym(mysqlcmd);
-
-
-
-
-%%get snake if possible, else use difspot silhouette
-%if mym(sprintf('select !isnull(bbID) from %sextspot where extspotID=%d',parameters.acq.name,struct_id)) == true;
\ No newline at end of file
diff --git a/4_grains/gtGetOmega_360.m b/4_grains/gtGetOmega_360.m
deleted file mode 100644
index a24a0d1472e03587a6fc0cebe07441096be14c77..0000000000000000000000000000000000000000
--- a/4_grains/gtGetOmega_360.m
+++ /dev/null
@@ -1,30 +0,0 @@
-function Omega = gtGetOmega_360(pair_ids)
-%get the omega angles corresponding to MaxImage of struct_ids vector
-%if there is a difA spot - omega is in range 0-180,
-%else if only a difB spot then omega will be in range 180-360
-%returns angle in DEGREES!
-
-%change to database query
-parameters=[];
-acq=parameters.acq;%to shorten sprintfs!
-Omega=[];
-
-
-for i=1:length(pair_ids)
-
- mysqlcmd = sprintf(['select ifnull(%sdifspot.MaxImage,%sdifspot.MaxImage+%d) '...
- 'from %s left join %sdifspot on difAID=%sdifspot.difspotID '...
- 'left join %sdifspot on difBID=%sdifspot.difspotID '...
- 'where %s.pairID=%d'],...
- acq.difA_name, acq.difB_name, acq.nproj,...
- acq.pair_tablename, acq.difA_name, acq.difA_name,...
- acq.difB_name, acq.difB_name,...
- acq.pair_tablename,pair_ids(i));
-
- Omega(i) = mym(mysqlcmd);
-
-end
-
-Omega = 180*(Omega./parameters.acq.nproj);
-
-end
\ No newline at end of file
diff --git a/4_grains/gtGetOmega_pair.m b/4_grains/gtGetOmega_pair.m
deleted file mode 100644
index d8798f4aeac78593a7f300c1cb08478fbddb55b6..0000000000000000000000000000000000000000
--- a/4_grains/gtGetOmega_pair.m
+++ /dev/null
@@ -1,31 +0,0 @@
-function Omega = gtGetOmega_pair(struct_ids, pair_vector, parameters)
-%get the omega angles corresponding to MaxImage of struct_ids vector
-%pair_vector - 1 or 2 for this half or other half of the scan
-%current directory.
-%returns angle in radians!
-
-
-
-%change to database query
-if isempty(parameters)
- load('parameters.mat');
-end
-
-Omega=[];
-for i=1:length(struct_ids)
-
- if pair_vector(i)==2
- Omega(i) = parameters.acq.nproj+mym(sprintf('select MaxImage from %sdifspot where difspotID = %d',parameters.acq.pair_name,struct_ids(i)));
- elseif pair_vector(i)==1
- Omega(i) = mym(sprintf('select MaxImage from %sdifspot where difspotID = %d',parameters.acq.name,struct_ids(i)));
- end
-
-end
-Omega = pi*(Omega./parameters.acq.nproj);
-warning('messing with omega!')
-%warning('Changed sign in gtGetOmega_pair!!! For ID11 data... ak')
-%Omega = -Omega;
-
-
-
-end
\ No newline at end of file
diff --git a/4_grains/gtGetSummedExtSpot.m b/4_grains/gtGetSummedExtSpot.m
deleted file mode 100644
index 40b057993d9bd7bec5adb422f2ac9017b954a873..0000000000000000000000000000000000000000
--- a/4_grains/gtGetSummedExtSpot.m
+++ /dev/null
@@ -1,90 +0,0 @@
-function im=gtGetSummedExtSpot(struct_id,varargin)
-% returns the summed full extinction spot image
-% conversion to database complete
-
-parameters=[];
-
-if isempty(parameters)
- disp('Loading parameter file')
- load('parameters.mat');
-end
-
-% connect to database
-connected=0;
-if ~isempty(varargin)
- connected=varargin{1};
-end
-if ~connected
-gtDBConnect
-end
-
-warning('contains proper bodged bb faffing, make nice after rerunning autofind')
-
-
-[bb(1),bb(2),bb(3),bb(4)]=mym(sprintf('select Xorigin,Yorigin,Xsize,Ysize from %sbb inner join %sbboxes on %sbb.bbID=%sbboxes.bboxID where %sbb.extspotID = %d', ...
- parameters.acq.name,...
- parameters.acq.name,...
- parameters.acq.name,...
- parameters.acq.name,...
- parameters.acq.name,...
- struct_id));
-
-%get snake if possible, else use difspot silhouette
-if mym(sprintf('select !isnull(bbID) from %sextspot where extspotID=%d',...
- parameters.acq.name,...
- struct_id)) == true;
-
-
-%put snake relative to boundingbox
-mysqlcmd=sprintf('select x,y from %ssnakes where snakeID=%d',parameters.acq.name,struct_id);
-
-[snakex,snakey]=mym(mysqlcmd);
-snakex = snakex - floor(bb(1));
-snakey = snakey - floor(bb(2));
-mask = poly2mask(snakex, snakey, bb(4), bb(3));
-
-
-
-else
- disp('Not using snake')
-%if we don't have a snake
-summeddif=gtGetSummedDifSpot(struct_id);
-
-%chqnge all this - save snake input in table during autofind, and reread
-%after. This recreation is stupuid
-summeddif_bw=summeddif>0.25*max(summeddif(:));
-summeddif_bw=bwmorph(summeddif_bw,'spur',inf);
-summeddif_bw=bwmorph(summeddif_bw,'erode',1);
-mask=bwmorph(summeddif_bw,'dilate',1);
-%%crop mask here
-tmp2=regionprops(double(mask),'BoundingBox');
-tmp2=tmp2.BoundingBox;
-
-tmp2=ceil(tmp2);%gt convention
-
-mask = gtCrop(mask,tmp2);
-bb(3:4)=tmp2(3:4);%change this!
-
-
-end
-
-%put BoundingBox into coordinates of the full image
-%not needed if using ext images, only if using the middle bit of fulls
-%bb(1) = ceil(bb(1) + parameters.acq.bb(1));
-%bb(2) = ceil(bb(2) + parameters.acq.bb(2));
-
-im=zeros(bb(4), bb(3));
-
-mysqlcmd=sprintf('select ExtStartImage,ExtEndImage from %sdifspot where difspotID=%d',parameters.acq.name,struct_id);
-[first,last]=mym(mysqlcmd);
-
-%disp([num2str(last-first+1) ' images in sum'])
-
-for i=first:last
- im = im + edf_read(sprintf('%s/1_preprocessing/ext/ext%04d.edf',parameters.acq.dir,i), bb);
-end
-
-%apply mask , either from snake or from difspot silhouette
-im = im.*mask;
-
-
diff --git a/4_grains/gtGetSummedExtSpot_pair.m b/4_grains/gtGetSummedExtSpot_pair.m
deleted file mode 100644
index 008069a7b365a8009ba4d75b8c0454f79a45a3de..0000000000000000000000000000000000000000
--- a/4_grains/gtGetSummedExtSpot_pair.m
+++ /dev/null
@@ -1,105 +0,0 @@
-function im=gtGetSummedExtSpot_pair(struct_id,pair,varargin)
-% returns the summed full extinction spot image
-% conversion to database complete
-%pair = 1 this dataset, pair = 2 pair dataset
-
-parameters=[];
-
-if isempty(parameters)
- disp('Loading parameter file')
- load('parameters.mat');
-end
-
-% connect to database
-connected=0;
-if ~isempty(varargin)
- connected=varargin{1};
-end
-if ~connected
-gtDBConnect
-end
-
-warning('contains proper bodged bb faffing, make nice after rerunning autofind')
-
-%get the name for this scan, or the pair dataset
-if pair == 2
- name=parameters.acq.pair_name;
-elseif pair ==1
- name=parameters.acq.name;
-end
-%path for reading images (assume pair directory is in same parent
-%directory)
-path=sprintf('../%s/',name);
-
-
-[bb(1),bb(2),bb(3),bb(4)]=mym(sprintf('select Xorigin,Yorigin,Xsize,Ysize from %sbb inner join %sbboxes on %sbb.bbID=%sbboxes.bboxID where %sbb.extspotID = %d', ...
- name, name, name, name, name, ...
- struct_id));
-
-%get snake if possible, else use difspot silhouette
-if mym(sprintf('select !isnull(bbID) from %sextspot where extspotID=%d',...
- name, struct_id)) == true;
-
-%put snake relative to boundingbox
-mysqlcmd=sprintf('select x,y from %ssnakes where snakeID=%d', name, struct_id);
-
-[snakex,snakey]=mym(mysqlcmd);
-snakex = snakex - floor(bb(1));
-snakey = snakey - floor(bb(2));
-mask = poly2mask(snakex, snakey, bb(4), bb(3));
-
-
-
-else
- disp('Not using snake')
-%if we don't have a snake
-summeddif=gtGetSummedDifSpot_pair(struct_id, pair);
-
-%change all this - save snake input in table during autofind, and reread
-%after. This recreation is stupuid
-summeddif_bw=summeddif>0.25*max(summeddif(:));
-summeddif_bw=bwmorph(summeddif_bw,'spur',inf);
-summeddif_bw=bwmorph(summeddif_bw,'erode',1);
-mask=bwmorph(summeddif_bw,'dilate',1);
-%%crop mask here
-tmp2=regionprops(double(mask),'BoundingBox');
-tmp2=tmp2.BoundingBox;
-
-tmp2=ceil(tmp2);%gt convention
-
-mask = gtCrop(mask,tmp2);
-bb(3:4)=tmp2(3:4);%change this!
-
-
-end
-
-%put BoundingBox into coordinates of the full image
-%not needed if using ext images, only if using the middle bit of fulls
-%bb(1) = ceil(bb(1) + parameters.acq.bb(1));
-%bb(2) = ceil(bb(2) + parameters.acq.bb(2));
-
-im=zeros(bb(4), bb(3));
-
-mysqlcmd=sprintf('select ExtStartImage,ExtEndImage from %sdifspot where difspotID=%d',name,struct_id);
-[first,last]=mym(mysqlcmd);
-
-%can have a problem because extspot bb extends beyond the edge of the ext image,
-%which has the size of parameters.acq.bb
-if bb(2)+bb(4) > parameters.acq.bb(4)
- bb(4) = parameters.acq.bb(4)-bb(2);
-end
-if bb(1)+bb(3) > parameters.acq.bb(3)
- bb(3) = parameters.acq.bb(3)-bb(1);
-end
-tmp_im=zeros(bb(4), bb(3)); % cropped image size
-
-for i=first:last
- tmp_im = tmp_im + edf_read(sprintf('%s/1_preprocessing/ext/ext%04d.edf',path,i), bb);
-end
-
-im=gtPlaceSubImage(tmp_im, im, 1, 1); % put the cropped image back into the correct size
-
-%apply mask , either from snake or from difspot silhouette
-im = im.*mask;
-
-
diff --git a/4_grains/gtIsPointInSample.m b/4_grains/gtIsPointInSample.m
deleted file mode 100644
index eef5275bc0e6a4a90732b853442f760d6b7ad3af..0000000000000000000000000000000000000000
--- a/4_grains/gtIsPointInSample.m
+++ /dev/null
@@ -1,25 +0,0 @@
-%
-% out=gtIsPointInSample(pcoord,sample_rad,sample_top,sample_bot,thr_geo)
-%
-% True, if the given in point is inside the illuminated cylindrical sample
-% volume. It works in LAB or Sample system.
-%
-% INPUT pcoord = [X,Y,Z] coordinates of the point
-% sample.rad = sample radius
-% sample.top = top of the sample
-% sample.bot = bottom of sample
-% tol_geo = geometrical tolerance (normally pos.)
-%
-% /all the input parameters must be in the same system and units/
-%
-% OUTPUT out = true if the given point is in the sample
-%
-
-function out=gtIsPointInSample(pcoord,sample,tol)
-
-out1=(pcoord(:,3) < sample.top+tol.geo) & (pcoord(:,3) > sample.bot-tol.geo);
-out2=(sqrt(pcoord(:,1).*pcoord(:,1)+pcoord(:,2).*pcoord(:,2)) < sample.rad+tol.geo);
-
-out=(out1&out2);
-
-end % of function
diff --git a/4_grains/gtMakeGrain.m b/4_grains/gtMakeGrain.m
deleted file mode 100644
index d2c49b80e7b24d13465e58115173321430bba08e..0000000000000000000000000000000000000000
--- a/4_grains/gtMakeGrain.m
+++ /dev/null
@@ -1,39 +0,0 @@
-%
-% Seperated from gtHandleSpot_360.
-
-function s = gtMakeGrain(sino, Omega, x_grain, y_grain, plot_flag)
-
- %back project good spots only
-
- binsino = sino>0.1*max(sino(:));
- img2=backpro(binsino,(pi/180)*Omega);
-
-
-
- %do imreconstruct from centre to threshold grain
- x_grain = ceil(x_grain); y_grain = ceil(y_grain);
- marker=zeros(size(img2));
- marker(y_grain, x_grain)=1;
- mask=double((img2>0.8*img2(y_grain, x_grain)));
- mask=imfill(mask,'holes');
-
- try
- grain=imreconstruct(marker,mask);
- if plot_flag
- figure
- h1=imshow(cat(3,img2,zeros(size(img2)),zeros(size(img2))),[]);
- hold on
- h2=imshow(cat(3,zeros(size(img2)),grain,zeros(size(img2))),[]);
- set(h2,'alphadata',0.5);
- end
- s=regionprops(grain,'Area','BoundingBox','Centroid');
-
- catch
- disp('Marker/mask problem - bombed out')
- accept=false;
- s=[];
- return
- end
-
-
-end
diff --git a/4_grains/gtMakeGrain3D.m b/4_grains/gtMakeGrain3D.m
deleted file mode 100644
index d2c49b80e7b24d13465e58115173321430bba08e..0000000000000000000000000000000000000000
--- a/4_grains/gtMakeGrain3D.m
+++ /dev/null
@@ -1,39 +0,0 @@
-%
-% Seperated from gtHandleSpot_360.
-
-function s = gtMakeGrain(sino, Omega, x_grain, y_grain, plot_flag)
-
- %back project good spots only
-
- binsino = sino>0.1*max(sino(:));
- img2=backpro(binsino,(pi/180)*Omega);
-
-
-
- %do imreconstruct from centre to threshold grain
- x_grain = ceil(x_grain); y_grain = ceil(y_grain);
- marker=zeros(size(img2));
- marker(y_grain, x_grain)=1;
- mask=double((img2>0.8*img2(y_grain, x_grain)));
- mask=imfill(mask,'holes');
-
- try
- grain=imreconstruct(marker,mask);
- if plot_flag
- figure
- h1=imshow(cat(3,img2,zeros(size(img2)),zeros(size(img2))),[]);
- hold on
- h2=imshow(cat(3,zeros(size(img2)),grain,zeros(size(img2))),[]);
- set(h2,'alphadata',0.5);
- end
- s=regionprops(grain,'Area','BoundingBox','Centroid');
-
- catch
- disp('Marker/mask problem - bombed out')
- accept=false;
- s=[];
- return
- end
-
-
-end
diff --git a/4_grains/gtMoveDifspotsIntoSingleTable.m b/4_grains/gtMoveDifspotsIntoSingleTable.m
deleted file mode 100644
index c834c1da69987005e1c8f0ac36bbb164572af988..0000000000000000000000000000000000000000
--- a/4_grains/gtMoveDifspotsIntoSingleTable.m
+++ /dev/null
@@ -1,95 +0,0 @@
-%
-% FUNCTION gtMoveDifspotsIntoSingleTable(parameters,deleteolddata)
-%
-% Copies 2 difspot tables into a single one.
-%
-% OPTIONAL INPUT
-% parameters: if empty, it will look for the parameters file, otherwise
-% define the following fields:
-% parameters.acq.nproj
-% parameters.acq.difA_name (without 'difspot' at the end)
-% parameters.acq.difB_name (without 'difspot' at the end)
-% deletetolddata: if true, old tables will be dropped afterwards.
-%
-
-function gtMoveDifspotsIntoSingleTable(parameters,deletetolddata)
-
-if ~exist('deletetolddata','var')
- deletetolddata=[];
-end
-
-if ~exist('parameters','var')
- parameters=[];
-end
-
-par_loaded=false;
-if isempty(parameters)
- par_loaded=true;
- load('parameters.mat');
-end
-
-nproj=parameters.acq.nproj;
-difnameA=[parameters.acq.difA_name 'difspot'];
-difnameAarch=[parameters.acq.difA_name 'difspot_xxx'];
-difnameB=[parameters.acq.difB_name 'difspot'];
-difnameBarch=[parameters.acq.difB_name 'difspot_xxx'];
-newtablename=parameters.acq.difA_name;
-
-
-gtDBConnect('admin')
-
-disp(' ')
-disp('Processing tables...')
-disp(' Values of CentroidImage will be taken from MaxImage.')
-disp(' ')
-
-mym(sprintf('DROP TABLE IF EXISTS %s',difnameAarch));
-mym(sprintf('DROP TABLE IF EXISTS %s',difnameBarch));
-
-mym(sprintf('RENAME TABLE %s TO %s',difnameA,difnameAarch))
-mym(sprintf('RENAME TABLE %s TO %s',difnameB,difnameBarch))
-
-gtDBCreateDifspotTable(newtablename)
-newtablename=[newtablename 'difspot'];
-
-disp(sprintf('Copying data from %s...',difnameA))
-mysqlcmd=sprintf(['INSERT INTO %s '...
- '(StartImage,EndImage,MaxImage,CentroidImage,ExtStartImage,ExtEndImage,'...
- 'Area,CentroidX,CentroidY,BoundingBoxXorigin,'...
- 'BoundingBoxYorigin,BoundingBoxXsize,BoundingBoxYsize,Integral) '...
- 'SELECT '...
- 'StartImage,EndImage,MaxImage,MaxImage,ExtStartImage,ExtEndImage,'...
- 'Area,CentroidX,CentroidY,BoundingBoxXorigin,'...
- 'BoundingBoxYorigin,BoundingBoxXsize,BoundingBoxYsize,Integral '...
- 'FROM %s'],newtablename,...
- difnameAarch);
-mym(mysqlcmd);
-
-disp(sprintf('Copying data from %s...',difnameB))
-mysqlcmd=sprintf(['INSERT INTO %s '...
- '(StartImage,EndImage,MaxImage,CentroidImage,ExtStartImage,ExtEndImage,'...
- 'Area,CentroidX,CentroidY,BoundingBoxXorigin,'...
- 'BoundingBoxYorigin,BoundingBoxXsize,BoundingBoxYsize,Integral) '...
- 'SELECT '...
- 'StartImage+%d,EndImage+%d,MaxImage+%d,MaxImage+%d,ExtStartImage+%d,ExtEndImage+%d,'...
- 'Area,CentroidX,CentroidY,BoundingBoxXorigin,'...
- 'BoundingBoxYorigin,BoundingBoxXsize,BoundingBoxYsize,Integral '...
- 'FROM %s'],newtablename,nproj,nproj,nproj,nproj,nproj,nproj,...
- difnameBarch);
-mym(mysqlcmd);
-
-if deletetolddata
- mym(sprintf('DROP TABLE %s',difnameAarch));
- mym(sprintf('DROP TABLE %s',difnameBarch));
- disp(sprintf('Table %s and %s have been dropped.',difnameA,difnameB))
- parameters.acq.difA_name=newtablename;
- parameters.acq.difB_name=[];
- if par_loaded
- save('parameters.mat',parameters)
- disp('Table names in the parameters file have been updated.')
- end
-end
-
-disp('Finished.')
-
-end % of function
diff --git a/4_grains/gtPairs2Grains360.m b/4_grains/gtPairs2Grains360.m
deleted file mode 100644
index 246438b8cd3c59e04f990f4e94d6ad9ab7e98556..0000000000000000000000000000000000000000
--- a/4_grains/gtPairs2Grains360.m
+++ /dev/null
@@ -1,354 +0,0 @@
-function [difAID, difBID, r_vector, grain_pos] = gtPairs2Grains360(x, phaseid)
-% <<<< NOT USED >>>>
-%for Andreas2_ multilayer dataset. Try indexing based on scattering
-%vectors only.
-%Simple approach for the small dataset
-
-dev_test=2;%what angular tolerence? (degrees)
-real_test=2;%what real space tolerance? (grainsize/real_test)
-plot_flag=1; %plot figures? - slow because of looking up spot images
-
-gtDBConnect
-parameters=[];
-load('parameters.mat');
-
-%first need to pick a large, paired 111 or 002 (or anything?)
-candidates = sfSelectPairedCandidates;
-%candidates is the pairID
-
-candidate=candidates(x);
-
-
-
-%get the data that will be used
-[pairIDs,plX,plY,plZ,thetatype]=mym('select pairID, plX, plY, plZ, thetatype from Andreas2_spotpairs where thetatype=1 or thetatype=2');
-[lorXA, lorYA, lorZA,lorXB, lorYB, lorZB, ldirX, ldirY, ldirZ]=mym('select samcentXA, samcentYA, samcentZA,samcentXB, samcentYB, samcentZB, ldirX, ldirY, ldirZ from Andreas2_spotpairs where thetatype=1 or thetatype=2');
-[difAIDs, difBIDs]=mym('select difAID, difBID from Andreas2_spotpairs where thetatype=1 or thetatype=2');
-
-linedata=[lorXA, lorYA, lorZA, ldirX, ldirY, ldirZ];
-keyboard
-
-plane_normals = zeros(length(plX), 3);
-plane_normals(:,1) = plX;
-plane_normals(:,2) = plY;
-plane_normals(:,3) = plZ;
-hkls=zeros(length(thetatype), 3);
-hkls(find(thetatype==1),:)=repmat([1 1 1],length(find(thetatype==1)),1);
-hkls(find(thetatype==2),:)=repmat([0 0 2],length(find(thetatype==2)),1);
-
-%deal with the first two diffraction rings only - low multiplicities means
-%less spots to consider - mixed blessing.
-allowed_angles11=sfGetAllowedAngles([1 1 1], [1 1 1], parameters.cryst(phaseid).crystal_system);
-allowed_angles12=sfGetAllowedAngles([1 1 1], [0 0 2], parameters.cryst(phaseid).crystal_system);
-allowed_angles22=sfGetAllowedAngles([0 0 2], [0 0 2], parameters.cryst(phaseid).crystal_system);
-allowed_angles11=repmat(allowed_angles11,length(plane_normals),1);
-allowed_angles12=repmat(allowed_angles12,length(plane_normals),1);
-allowed_angles22=repmat(allowed_angles22,length(plane_normals),1);
-
-%pick out the candidate
-testpair=find(pairIDs==candidate);
-
-subject_normal=plane_normals(testpair,:);
-subject_normal=repmat(subject_normal, size(plane_normals,1), 1);
-%interplanar angles in degrees
-angles=(180/pi)*acos(dot(subject_normal, plane_normals, 2));
-
-shortlist=[];
-
-if thetatype(testpair)==1
- %test against all 111-111 allowed angles
- dev_angles = abs(repmat(angles,1,size(allowed_angles11,2)) - allowed_angles11);
- dev_angles = min(dev_angles,[],2);
- shortlist=find(dev_angles<dev_test & thetatype==1);
- %plus do the 111-002 combinations
- dev_angles = abs(repmat(angles,1,size(allowed_angles12,2)) - allowed_angles12);
- dev_angles = min(dev_angles,[],2);
- shortlist = [shortlist; find(dev_angles<dev_test & thetatype==2)];
-elseif thetatype(testpair)==2
- %test against all 002-002 allowed angles
- dev_angles = abs(repmat(angles,1,size(allowed_angles22,2)) - allowed_angles22);
- dev_angles = min(dev_angles,[],2);
- shortlist=find(dev_angles<dev_test & thetatype==2);
- %plus do the 111-002 combinations
- dev_angles = abs(repmat(angles,1,size(allowed_angles12,2)) - allowed_angles12);
- dev_angles = min(dev_angles,[],2);
- shortlist = [shortlist; find(dev_angles<dev_test & thetatype==1)];
-end
-
-keyboard
-
-%This is a shortlist of possible plane normals according to scattering
-%vector direction.
-%...now, do real space
-
-%get difspot data
-[difAID, difBID]=mym(sprintf('select difAID, difBID from Andreas2_spotpairs where pairID=%d',candidate));
-[spotsizex(1), spotsizey(1)]=mym(sprintf('select boundingboxXsize, boundingboxYsize from Andreas2_difspot where difspotid=%d',difAID));
-[spotsizex(2), spotsizey(2)]=mym(sprintf('select boundingboxXsize, boundingboxYsize from Andreas2_difspot where difspotid=%d',difBID));
-spotsizex=mean(spotsizex);
-spotsizey=mean(spotsizey);
-spotsize=mean([spotsizex spotsizey]);
-
-distance=zeros(size(shortlist));
-points=zeros(length(shortlist),3);
-
-%quadratic... at^2 + bt +c = 0
-a= linedata(testpair,4)^2 + linedata(testpair,5)^2;
-b= -2*((linedata(testpair,1)*linedata(testpair,4))+(linedata(testpair,2)*linedata(testpair,5)));
-c= linedata(testpair,1)^2 + linedata(testpair,2)^2 - 200^2;
-
-t(1)= (-b+sqrt(b^2-(4*a*c)))/(2*a);
-t(2)= (-b-sqrt(b^2-(4*a*c)))/(2*a);
-tmin=min(t);
-tmax=max(t);
-
-
-
-if plot_flag
-%testpair line:
-point1=linedata(testpair,1:3);
-point1=point1-tmin*linedata(testpair,4:6);
-point2=point1-tmax*linedata(testpair,4:6);
-pdata=[point1; point2];
-figure
-hold on
-plot3(pdata(:,1), pdata(:,2), pdata(:,3),'r')
-end
-
-%pick out intersections from which a grain position can be refined
-grainpoints=[];
-for a=1:length(shortlist)
- point=pofintersectexp(linedata([testpair shortlist(a)], :))'; %best point between the two lines
- distance(a)=pointtolinedist(point, linedata(testpair, :));
- points(a,:)=point;
-
- %if lines are ~ parallel, this gets ugly
- if distance(a)>25
- plot3(point(1), point(2), point(3), 'go')
- elseif distance(a) >10
- plot3(point(1), point(2), point(3), 'ro')
- else
- if angles(shortlist(a))>5 %only plot if the angle is greater than 5 degrees
- %similarly, can only get grain position from lines more than a few
- %degrees apart.
- if plot_flag
- plot3(point(1), point(2), point(3), 'bo')
- end
- grainpoints=[grainpoints; point];%use to refine position
- end
- end
-end
-
-
-%then, refine grain position
-%need outlier rejecting least squares approach from grainpoints.
-%to work in a single variable, use t, where t is a coordinate along the path
-%of the original spot
-
-%hard code some sensible limits
-%tmin=3000
-%tmax=4000
-
-
-t_grain = fminbnd(@(t) sfFindGrain(grainpoints, t), tmin, tmax);
-
-x_grain=linedata(testpair,1)-t_grain*linedata(testpair,4);
-y_grain=linedata(testpair,2)-t_grain*linedata(testpair,5);
-z_grain=linedata(testpair,3)-t_grain*linedata(testpair,6);
-
-plot3(x_grain, y_grain, z_grain, 'xr', 'markersize', 20)
-
-%Finally (?) reject those lines which intersect too far from the grain
-%position
-dist = points-repmat([x_grain y_grain z_grain], size(points, 1), 1);
-dist = sqrt(sum((dist.*dist),2)); % distance of each point from grain
-shortlist2 = shortlist;
-shortlist2(find(dist>spotsize/real_test))=[];
-points2=points;
-points2(find(dist>spotsize/real_test), :)=[]; %points is now same length as shortlist2
-
-
-if plot_flag
- %look at these spots
- figure
- for a=1:length(shortlist2)
- im=gtGetSummedDifSpot(difAIDs(shortlist2(a)));
- subplot(2,length(shortlist2), a); imshow(im, [])
- im=gtGetSummedDifSpot(difBIDs(shortlist2(a)));
- subplot(2,length(shortlist2), a+length(shortlist2)); imshow(im, [])
- end
-end
-
-%run consistancy test of the final group of spots
-[r_vector, goods]=plot_rodrigues_consistancy(plane_normals(shortlist2,:), hkls(shortlist2,:), parameters.cryst(phaseid).spacegroup, plot_flag);
-points3=points2(goods,:);
-shortlist3=shortlist2(goods);
-
-if plot_flag
- figure
- for a=1:length(shortlist3)
- im=gtGetSummedDifSpot(difAIDs(shortlist3(a)));
- subplot(2,length(shortlist3), a); imshow(im, [])
- im=gtGetSummedDifSpot(difBIDs(shortlist3(a)));
- subplot(2,length(shortlist3), a+length(shortlist3)); imshow(im, [])
- end
-end
-
-
-%for output
-difAID=difAIDs(shortlist3);
-difBID=difBIDs(shortlist3);
-
-%should refine grain position again, having done the consistancy check
-t_grain = fminbnd(@(t) sfFindGrain(points3, t), tmin, tmax);
-
-x_grain=linedata(testpair,1)-t_grain*linedata(testpair,4);
-y_grain=linedata(testpair,2)-t_grain*linedata(testpair,5);
-z_grain=linedata(testpair,3)-t_grain*linedata(testpair,6);
-
-grain_pos=[x_grain y_grain z_grain];
-
-%then, call forward simulate to get anything missed, any unpaired spots,
-%and any higher order {hkl} pairs.
-%assemble data
-struct_ids=[difAID; difBID];
-[struct_ids_new] = gtForwardSimulate_pair2grain360(grain_pos, r_vector, struct_ids, parameters);
-
-%now, assemble all info for these new spots.
-struct_ids=[struct_ids; struct_ids_new'];
-grainTheta=[]; grainEta=[]; grainOmega=[]; GrainPlane_normals=[];
-
-for ll=1:length(struct_ids)
- %get info from pair table if spots are paired
-
- mysqlcmd=sprintf('select theta, eta, omega, plX, plY, plZ from Andreas2_spotpairs where difAID=%d',struct_ids(ll));
-
- [th,et,om,plX,plY,plZ]=mym(mysqlcmd);
-
- if ~isnan(th)
- if om<=180
- grainTheta(ll)=th; grainEta(ll)=et; grainOmega(ll)=om;
- grainPlane_normals(ll,1:3)=[plX plY plZ];
- else
- grainTheta(ll)=th; grainEta(ll)=180-et; grainOmega(ll)=om+180;
- if grainEta(ll)<0
- grainEta(ll)=grainEta(ll)+360;
- end
- grainPlane_normals(ll,1:3)=[-plX -plY -plZ];
- end
-
- else%if no info found in pair table, calculate it from grain position
- [th,et,om,plX,plY,plZ]=sfCalcAngles(struct_ids(ll), grain_pos);
- grainTheta(ll)=th; grainEta(ll)=et; grainOmega(ll)=om;
- grainPlane_normals(ll,1:3)=[plX plY plZ];
- end
-
-end
-
-%consistancy check again to be sure, shouldn't discard anything...
-%at this point, really need to have hkls assigned. For pairs, this in table,
-%but for others it needs to be done from scratch
-%note, gtIndexReflections can also return plane normals
-[tmp, grainHKLs] = gtIndexReflections(grainTheta, grainEta, grainOmega, parameters);
-
-
-keyboard
-%do we need a pseudo-twin ckeck somewhere?
-
-
-%remove everything claimed by a grain from the lists, and go again with the
-%next candidate.
-
-
-
-
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-function output = sfFindGrain(good_points, t_grain)
-%function for fminbnd to find grain position along original projection
-x_grain=linedata(testpair,1)-t_grain*linedata(testpair,4);
-y_grain=linedata(testpair,2)-t_grain*linedata(testpair,5);
-z_grain=linedata(testpair,3)-t_grain*linedata(testpair,6);
-
-dist = good_points-repmat([x_grain y_grain z_grain], size(good_points, 1), 1);
-dist = sqrt(sum((dist.*dist),2)); % distance of each point from grain
-%minimise sum of squares for the closest half of the points, to reduce
-%effect of outliers
-dist = sort(dist,1,'descend');
-tmp = ceil(size(dist,1)/2);
-dist(1:tmp)=[];
-output = dist'*dist;%sum of squares
-end
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- function cand = sfSelectPairedCandidates
- %returns all paired thetatype=1 pairids
- difA_name='Andreas2_';
- pairtablename='Andreas2_';
- cand=mym(sprintf([...
- 'select %sspotpairs.pairID from %sspotpairs inner join %sdifspot on difAID=difspotID '...
- 'where thetatype=1 or thetatype=2 order by boundingboxYsize desc'],...
- pairtablename,pairtablename,difA_name));
- end
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-function allowed_angles = sfGetAllowedAngles(hkl_1,hkl_2,crystal_system)
-%get the allowed interplanar angles
-hkl_m1 = gtCrystSignedHKLs( hkl_1, gtCrystGetSymmetryOperators(crystal_system) );
-hkl_m2 = gtCrystSignedHKLs( hkl_2, gtCrystGetSymmetryOperators(crystal_system) );
-allowed_angles = [];counter = 1;
-for l=1:size(hkl_m1,1)
- for jj=1:size(hkl_m2,1)
- tmp = acos ( dot(hkl_m1(l,:), hkl_m2(jj,:)) / (sqrt( (hkl_m1(l,:)*hkl_m1(l,:)')*(hkl_m2(jj,:)*hkl_m2(jj,:)') )) );
- tmp = tmp*180/pi;
- allowed_angles(counter) = tmp;
- counter = counter+1;
- end
-end
-allowed_angles=unique(allowed_angles);
-end
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-function [th, et, om, normalX normalY normalZ]= sfCalcAngles(struct_id, grainPos_sample)
-%get difspot pos
-
-[spotX,spotY,spotIm]=mym(sprintf('select CentroidX,CentroidY,MaxImage from Andreas2_difspot where difspotID=%d',struct_id));
-
-%calc omega
-om=180*(spotIm/parameters.acq.nproj);
-
-%put grain pos in setup (lab) reference
-grainPos_setup(1)= cosd(om)*grainPos_sample(1)+sind(om)*grainPos_sample(2) ;
-grainPos_setup(2)= -sind(om)*grainPos_sample(1)+cosd(om)*grainPos_sample(2) ;
-grainPos_setup(3)= grainPos_sample(3) ;
-
-%spot position in setup (lab) reference
-spotPos_setup=[(parameters.acq.dist/parameters.acq.pixelsize), spotX-parameters.acq.rotx, parameters.acq.bb(2)-spotY];
-
-%calc theta, eta
-deltaX = spotPos_setup(2)-grainPos_setup(2);
-deltaY = spotPos_setup(3)-grainPos_setup(3);%compared to gtMatchDifspots_v2b, sign of y has already been flipped
-grain_to_det = spotPos_setup(1)-grainPos_setup(1);
-
-th=acosd(grain_to_det / norm([grain_to_det deltaX deltaY]))/2;
-et=(180/pi)*atan2(deltaX, deltaY);
-if et<0
- et=360+et;
-end
-
-%calc plane normal (scattering vector)
-diff_vec_setup = [grain_to_det deltaX deltaY]/norm([grain_to_det deltaX deltaY]);
-pl_vec_setup = (diff_vec_setup-[1 0 0])/norm(diff_vec_setup-[1 0 0]) ;
-%transform to sample coordinates
-normalX = cosd(om)*pl_vec_setup(1)-sind(om)*pl_vec_setup(2) ;
-normalY = sind(om)*pl_vec_setup(1)+cosd(om)*pl_vec_setup(2) ;
-normalZ = pl_vec_setup(3) ;
-
-end
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-
-end
diff --git a/4_grains/gtPredictAngles.m b/4_grains/gtPredictAngles.m
deleted file mode 100644
index fc0b1710cb04abc9a6930103b63debb6fbe72bd3..0000000000000000000000000000000000000000
--- a/4_grains/gtPredictAngles.m
+++ /dev/null
@@ -1,76 +0,0 @@
-function [Theta,Eta,Omega]=gtPredictAngles(normal, hkl, parameters)
-
-%calculates the Omega position at which diffraction will occur from a plane
-%defined by its plane normal and hkl type [h k l]
-
-%30/10/2006 - modify to follow change in gtConsistancyCheck, and to produce
-%only the theta/eta/omega from the plane normal (not its inverse as well)
-
-%load data
-
-if isempty(parameters)
- load('parameters.mat');
-end
-
-%put plane normal in cylindrical polar coords
-length = sqrt(normal(1)^2 + normal(2)^2);
-alpha = atan2(normal(2), normal(1));
-vert = normal(3);
-
-%calc Theta for diffraction
-lambda = 12.398/parameters.acq.energy;
-
- %%%%%%%% cubic cases %%%%%%%%%%%%%
-if parameters.cryst.spacegroup == 225 || parameters.cryst.spacegroup == 227 || parameters.cryst.spacegroup == 229
- d=parameters.cryst.latticepar(1)/sqrt(hkl*hkl');
- % modif sabine
- %%%%%%%% hexagonal cases %%%%%%%%%%%%%
-elseif parameters.cryst.spacegroup == 663 || parameters.cryst.spacegroup == 194 || parameters.cryst.spacegroup == 167
- h=hkl(1);k=hkl(2) ;l=hkl(4);
- d=parameters.cryst.latticepar(1)/(4/3*(h^2+k^2+h*k)+(l*parameters.cryst.latticepar(1)/parameters.cryst.latticepar(3))^2)^0.5;
- % fin modif sabine
-else
- warning('unknown spacegroup!')
- return
-end
-Theta = asin(lambda/(2*d));
-
-%try
- %produce Omega angles - find four Omegas (hkl ->2, -h-k-l ->2)
- if ((cos((pi/2)-Theta))/length)<=1
- temp(1) = acos((cos((pi/2)-Theta))/length);
- temp(2) = -temp(1);
- Omega = alpha + pi + temp;
- else
- disp('no diffraction from this plane normal!')
- Theta=[];
- return
- end
-
- %calc Eta, make Theta to output
- %in our convention, Eta is measured clockwise from vertical on the
- %images.
-
- alpha2 = alpha - Omega;%plane normal angles at diffraction
-
- y = length*sin(alpha2);
- x = length*cos(alpha2);
-
- Eta = atan2(y,vert);%this was y, now +ve because of 'looking into beam" sense
- Theta(1:2) = Theta(1);%all Thetas are the same
-
- %put Omega in 0-2pi interval
- dummy = find(Omega<0);
- Omega(dummy) = Omega(dummy) + (2*pi);
- dummy = find(Omega>(2*pi));
- Omega(dummy) = Omega(dummy) - (2*pi);
-
- %put in degrees for human users
- Omega = Omega*180/pi;
- Theta = Theta*180/pi;
- Eta = Eta*180/pi;
-
-
-
-
-
diff --git a/4_grains/gtPredictAngles3D.m b/4_grains/gtPredictAngles3D.m
deleted file mode 100644
index c13c0fb9d5b594e20553075a9b69e62ed5b731cf..0000000000000000000000000000000000000000
--- a/4_grains/gtPredictAngles3D.m
+++ /dev/null
@@ -1,80 +0,0 @@
-function [Theta,Eta,Omega]=gtPredictAngles3D(normal, hkl, parameters)
-
-%calculates the Omega position at which diffraction will occur from a plane
-%defined by its plane normal and hkl type [h k l]
-
-%30/10/2006 - modify to follow change in gtConsistancyCheck, and to produce
-%only the theta/eta/omega from the plane normal (not its inverse as well)
-
-%load data
-
-if ~exist('parameters','var')
- load('parameters.mat');
-end
-
-%put plane normal in cylindrical polar coords
-length = sqrt(normal(1)^2 + normal(2)^2);
-alpha = atan2(normal(2), normal(1));
-vert = normal(3);
-
-%calc Theta for diffraction
-lambda = 12.398/parameters.acq.energy;
-
- %%%%%%%% cubic cases %%%%%%%%%%%%%
- %EDIT Erik begin
-if parameters.cryst.spacegroup == 225 || parameters.cryst.spacegroup == 227 || parameters.cryst.spacegroup == 229 || parameters.cryst.spacegroup == 221
- %EDIT Erik begin
- d=parameters.cryst.latticepar(1)/sqrt(hkl*hkl');
- % modif sabine
- %%%%%%%% hexagonal cases %%%%%%%%%%%%%
-elseif parameters.cryst.spacegroup == 663 || parameters.cryst.spacegroup == 194 || parameters.cryst.spacegroup == 167 || parameters.cryst.spacegroup == 152
- h=hkl(1);k=hkl(2) ;l=hkl(4);
- d=parameters.cryst.latticepar(1)/(4/3*(h^2+k^2+h*k)+(l*parameters.cryst.latticepar(1)/parameters.cryst.latticepar(3))^2)^0.5;
- % fin modif sabine
-else
- warning('unknown spacegroup ! ...aborting');
- return
-end
-
-Theta = asin(lambda/(2*d));
-
-%produce Omega angles - find four Omegas (hkl ->2, -h-k-l ->2)
-if ((cos((pi/2)-Theta))/length)<=1
- temp(1) = acos((cos((pi/2)-Theta))/length);
- temp(2) = -temp(1);
- Omega = alpha + pi + temp;
- %calc Eta, make Theta to output
- %in our new convention, Eta is measured clockwise from vertical on the
- %images when looking along the beam direction !
-
- alpha2 = alpha - Omega; %plane normal angles at diffraction
-
- y = length*sin(alpha2);
- x = length*cos(alpha2);
-
- Eta = atan2(-y,vert); %this was y, now -y because of 'looking along beam" sense
- Theta(1:2) = Theta(1); %all Thetas are the same
-
- %put Omega in 0-2pi interval
- dummy = find(Omega<0);
- Omega(dummy) = Omega(dummy) + (2*pi);
- dummy = find(Omega>(2*pi));
- Omega(dummy) = Omega(dummy) - (2*pi);
-
- %put in degrees for human users
- Omega = Omega*180/pi;
- Theta = Theta*180/pi;
- Eta = Eta*180/pi;
-else
- disp('no diffraction from this plane normal!')
- Theta=[];
- Eta=[];
- Omega=[];
-
-end
-
-
-
-
-
-
diff --git a/4_grains/gtPredictAngles_Si.m b/4_grains/gtPredictAngles_Si.m
deleted file mode 100644
index c93db5ba677bc3af4e4233a1b22bc60222f0ce7a..0000000000000000000000000000000000000000
--- a/4_grains/gtPredictAngles_Si.m
+++ /dev/null
@@ -1,74 +0,0 @@
-function [Theta,Eta,Omega]=gtPredictAngles_Si(normal, hkl)
-
-%calculates the Omega position at which diffraction will occur from a plane
-%defined by its plane normal and hkl type [h k l]
-
-%30/10/2006 - modify to follow change in gtConsistancyCheck, and to produce
-%only the theta/eta/omega from the plane normal (not its inverse as well)
-
-%_Si - simplified version - no parameters.mat
-%hard code stuff
-disp('hard coded variables - change in gtPredictAngles_Si')
-energy=40.2837
-latticepar=5.431
-
-%put plane normal in cylindrical polar coords
-length = sqrt(normal(1)^2 + normal(2)^2);
-alpha = atan2(normal(2), normal(1));
-vert = normal(3);
-
-%calc Theta for diffraction
-lambda = 12.398/energy;
-d=latticepar/sqrt(hkl*hkl');
-
-Theta = asin(lambda/(2*d));
-
-%try
- %produce Omega angles
- if ((cos((pi/2)-Theta))/length)<=1
- temp(1) = acos((cos((pi/2)-Theta))/length);
- temp(2) = -temp(1);
- Omega = alpha + pi + temp;
- else
- disp('no diffraction from this plane normal!')
- Theta=[];
- return
- end
-
- %calc Eta, make Theta to output
- %in our convention, Eta is measured clockwise from vertical on the
- %images.
-
- alpha2 = alpha - Omega;%plane normal angles at diffraction
-
- y = length*sin(alpha2);
- x = length*cos(alpha2);
-
- Eta = atan2(y,vert);%this was y, now +ve because of 'looking into beam" sense
- Theta(1:2) = Theta(1);%all Thetas are the same
-
-
-
- %put all in -(1/2)*pi - (3/2)*pi interval
- dummy = find(Omega<-pi/2);
- Omega(dummy) = Omega(dummy) + (2*pi);
- dummy = find(Omega>((1.5)*pi));
- Omega(dummy) = Omega(dummy) - (2*pi);
-
- %dummy = find(Eta<0);
- %Eta(dummy) = Eta(dummy) + (2*pi);
- %dummy = find(Eta>(2*pi));
- %Eta(dummy) = Eta(dummy) - (2*pi);
-
- %put in degrees for human users
- Omega = Omega*180/pi;
- Theta = Theta*180/pi;
- Eta = Eta*180/pi;
-
-%catch%if no solutions
-% disp('no diffraction from this plane normal!')
-%end
-
-
-
-
diff --git a/4_grains/gtPredictPositions.m b/4_grains/gtPredictPositions.m
deleted file mode 100644
index b68387d8ca6653a695c6eaab624d8dd8e1f49128..0000000000000000000000000000000000000000
--- a/4_grains/gtPredictPositions.m
+++ /dev/null
@@ -1,94 +0,0 @@
-%predict positions of spots in images
-%Marcelo, Andy, 30/10
-
-function [images,x_direct,y_direct,x_diffr,y_diffr] = gtPredictPositions(normal, hkl, grainid, plot_flag)
-%take a plane normal, its hkl type (to calc Theta), and a grainID (to calc
-%grain position in sample), and returns lists of max images, and x and y
-%positions of the ext and dif spots.
-%plot_flag=1 plots images
-
-%load data
-parameters=[];
-if isempty(parameters)
- load('parameters.mat');
-end
-
-grain_data = load(sprintf('4_grains/grain%d_/grain%d_.mat',grainid,grainid));
-grainx = grain_data.x1 + (grain_data.nx/2);
-grainy = grain_data.y1 + (grain_data.ny/2);
-zcenter = grain_data.zcenter;
-
-%use gtPredictAngles to get Theta, Eta, Omega - note gives a pair of angles
-[Theta,Eta,Omega] = gtPredictAngles(normal, hkl);
-%convert to radians
-Theta=Theta*pi/180;
-Eta=Eta*pi/180;
-Omega=Omega*pi/180;
-
-%get only Omegas within 0-180 degrees
-if parameters.acq.type=='180degree'
- %disp('180 degree data')
- dummy=find(Omega<0 | Omega>pi);
- Theta(dummy)=[];
- Eta(dummy)=[];
- Omega(dummy)=[];
-end
-
-%deal with grain position
-grain2rot=sqrt((grainx-parameters.acq.bb(3)/2).^2+(grainy-parameters.acq.bb(3)/2).^2);
-angle=atan2(parameters.acq.bb(3)/2-grainy,grainx-parameters.acq.bb(3)/2);
-
-%extinction spot position
-% y correction with the omega offset in the sample, that is a value to be added to X in the screen
-x_direct = parameters.acq.bb(1)+(parameters.acq.bb(3)/2)+(grain2rot.*cos(-Omega+angle));
-y_direct(1:length(Omega)) = parameters.acq.bb(2)+zcenter;
-
-%diffraction spot position
-% radius on the screen (in pixels) with the corrected distance from the sample to the
-% detector (y in sample coordinates) with the omega offset
-radius = tan(2*Theta).*((parameters.acq.dist/parameters.acq.pixelsize)+(grain2rot*sin(-Omega+angle)));
-x_diffr = x_direct + (radius.*sin(Eta)); % spot position (x position) in the screen
-y_diffr = y_direct - (radius.*cos(Eta)); % vertical coordinate on the screeen
-
-
-%get full image
- images = round((Omega/pi)*parameters.acq.proj);
-
- %plot results on full images
- if plot_flag
-
-for i=1:size(Omega,2)
-
- %convert Omega to image number
-
- start_image = images(i)-5;
- if start_image<0
- start_image=0;
- end
- end_image = images(i)+5;
- if end_image>5999
- end_image=5999;
- end
-
- start_image
- end_image
-
- full=zeros(parameters.acq.xdet,parameters.acq.ydet);
- for im=start_image:end_image
- full=full+edf_read(sprintf('1_preprocessed/full/full%04d.edf',im));
- end
- full=full./(end_image-start_image+1);
-
- figure
- imshow(full,[]);
- hold on;
- plot(x_direct(i),y_direct(i),'og')
- plot(x_diffr(i),y_diffr(i),'ro')
- drawnow
- disp('click image to close and see next spot (if exists...)')
- k=waitforbuttonpress;
-
-end
-
- end
-
diff --git a/4_grains/gtPredictPositions2.m b/4_grains/gtPredictPositions2.m
deleted file mode 100644
index a0871279d1ed7c280c03e6e699a8de4761bf2b4a..0000000000000000000000000000000000000000
--- a/4_grains/gtPredictPositions2.m
+++ /dev/null
@@ -1,104 +0,0 @@
-%predict positions of spots in images
-%version2 - designed to work in gtHandleSpot2 sequence, with no
-%grain%d_.mat yet written, all information passed as arguments
-
-%Marcelo, Andy, 30/10
-
-function [images,x_direct,y_direct,x_diffr,y_diffr] = gtPredictPositions2(normal, hkl, grainCentroid, plot_flag, parameters)
-%take a plane normal, its hkl type (to calc Theta), and a grainCentroid [x,y,z], and returns lists of max images, and x and y
-%positions of the ext and dif spots.
-%plot_flag=1 plots images
-
-%load data
-if isempty(parameters)
- load('parameters.mat');
-end
-grainx = grainCentroid(1);
-grainy =grainCentroid(2);
-zcenter = grainCentroid(3);
-
-%use gtPredictAngles to get Theta, Eta, Omega - note gives a pair of angles
-[Theta,Eta,Omega] = gtPredictAngles(normal, hkl);
-%convert to radians
-Theta=Theta*pi/180;
-Eta=Eta*pi/180;
-Omega=Omega*pi/180;
-
-%get only Omegas within 0-180 degrees for 180deg data
-if parameters.acq.type=='180degree'
- %disp('180 degree data')
- dummy=find(Omega<0 | Omega>pi);
- Theta(dummy)=[];
- Eta(dummy)=[];
- Omega(dummy)=[];
-end
-
-%deal with grain position
-grain2rot=sqrt((grainx-parameters.acq.bb(3)/2).^2+(grainy-parameters.acq.bb(3)/2).^2);
-angle=atan2(parameters.acq.bb(3)/2-grainy,grainx-parameters.acq.bb(3)/2);
-
-%extinction spot position
-% y correction with the omega offset in the sample, that is a value to be added to X in the screen
-x_direct = parameters.acq.bb(1)+(parameters.acq.bb(3)/2)+(grain2rot.*cos(-Omega+angle));
-y_direct(1:length(Omega)) = parameters.acq.bb(2)+zcenter;
-
-%diffraction spot position
-% radius on the screen (in pixels) with the corrected distance from the sample to the
-% detector (y in sample coordinates) with the omega offset
-radius = tan(2*Theta).*((parameters.acq.dist/parameters.acq.pixelsize)+(grain2rot*sin(-Omega+angle)));
-x_diffr = x_direct + (radius.*sin(Eta)); % spot position (x position) in the screen
-y_diffr = y_direct - (radius.*cos(Eta)); % vertical coordinate on the screeen
-
-%round to integer values
-x_direct = round(x_direct);
-y_direct = round(y_direct);
-x_diffr = round(x_diffr);
-y_diffr = round(y_diffr);
-
-%get full image
-images = round((Omega/pi)*parameters.acq.nproj);
-
-
-
-
- %~~~~~~~~plot results on full images~~~~~~~~~~~~~~~~~~~~~~
- if plot_flag
- for i=1:size(Omega,2)
-
- %convert Omega to image number
- if images(i) >= parameters.acq.nproj
- images(i)=images(i)-parameters.acq.nproj;
- path=sprintf('../%s/',parameters.acq.difB_name);
- else
- path=sprintf('../%s/',parameters.acq.difA_name);
- end
-
- start_image = images(i)-3;
- if start_image<0
- start_image=0;
- end
- end_image = images(i)+3;
- if end_image>parameters.acq.nproj-1;
- end_image=parameters.acq.nproj-1;
- end
-
- start_image
- end_image
-
- full=zeros(parameters.acq.xdet,parameters.acq.ydet);
- for im=start_image:end_image
- full=full+edf_read(sprintf('%s1_preprocessing/full/full%04d.edf',path,im));
- end
- full=full./(end_image-start_image+1);
-
- figure
- imshow(full,[]);
- hold on;
- plot(x_direct(i),y_direct(i),'og')
- plot(x_diffr(i),y_diffr(i),'ro')
- drawnow
- disp('click image to close and see next spot (if exists...)')
- k=waitforbuttonpress;
- end
- end
- %~~~~~~~~plot results on full images~~~~~~~~~~~~~~~~~~~~~~
diff --git a/4_grains/gtPredictPositions_360.m b/4_grains/gtPredictPositions_360.m
deleted file mode 100644
index 943915c38a212db6a9496a940dae5f8e7efb9af2..0000000000000000000000000000000000000000
--- a/4_grains/gtPredictPositions_360.m
+++ /dev/null
@@ -1,99 +0,0 @@
-%predict positions of spots in images
-%version2 - designed to work in gtHandleSpot2 sequence, with no
-%grain%d_.mat yet written, all information passed as arguments
-
-%Marcelo, Andy, 30/10
-
-function [images,x_direct,y_direct,x_diffr,y_diffr] = gtPredictPositions_360(normal, hkl, grainCentroid, plot_flag, parameters)
-%take a plane normal, its hkl type (to calc Theta), and a grainCentroid [x,y,z], and returns lists of max images, and x and y
-%positions of the ext and dif spots.
-%plot_flag=1 plots images
-
-%load data
-if isempty(parameters)
- load('parameters.mat');
-end
-grainx = grainCentroid(1);
-grainy =grainCentroid(2);
-zcenter = grainCentroid(3);
-
-%use gtPredictAngles to get Theta, Eta, Omega - note gives a pair of angles
-[Theta,Eta,Omega] = gtPredictAngles(normal, hkl, parameters);
-%convert to radians
-Theta=Theta*pi/180;
-Eta=Eta*pi/180;
-Omega=Omega*pi/180;
-
-%get only Omegas within 0-180 degrees for 180deg data
-if parameters.acq.type=='180degree'
- %disp('180 degree data')
- dummy=find(Omega<0 | Omega>pi);
- Theta(dummy)=[];
- Eta(dummy)=[];
- Omega(dummy)=[];
-end
-
-%deal with grain position
-grain2rot=sqrt((grainx-parameters.acq.bb(3)/2).^2+(grainy-parameters.acq.bb(3)/2).^2);
-angle=atan2(parameters.acq.bb(3)/2-grainy,grainx-parameters.acq.bb(3)/2);
-
-%extinction spot position
-% y correction with the omega offset in the sample, that is a value to be added to X in the screen
-x_direct = parameters.acq.bb(1)+(parameters.acq.bb(3)/2)+(grain2rot.*cos(-Omega+angle));
-y_direct(1:length(Omega)) = parameters.acq.bb(2)+zcenter;
-
-%diffraction spot position
-% radius on the screen (in pixels) with the corrected distance from the sample to the
-% detector (y in sample coordinates) with the omega offset
-radius = tan(2*Theta).*((parameters.acq.dist/parameters.acq.pixelsize)+(grain2rot*sin(-Omega+angle)));
-x_diffr = x_direct + (radius.*sin(Eta)); % spot position (x position) in the screen
-y_diffr = y_direct - (radius.*cos(Eta)); % vertical coordinate on the screeen
-
-%round to integer values
-x_direct = round(x_direct);
-y_direct = round(y_direct);
-x_diffr = round(x_diffr);
-y_diffr = round(y_diffr);
-
-%get full image
-images = round((Omega/pi)*parameters.acq.nproj);
-
- % modif sabine
- %plot_flag=1
-
- %~~~~~~~~plot results on full images~~~~~~~~~~~~~~~~~~~~~~
- if plot_flag
- for i=1:size(Omega,2)
-
- %convert Omega to image number
- start_image = images(i)-3;
- if start_image<0
- start_image=0;
- end
- end_image = images(i)+3;
- if end_image>(2*parameters.acq.nproj)-1;
- end_image=(2*parameters.acq.nproj)-1;
- end
-
- start_image
- end_image
-
- full=zeros(parameters.acq.xdet,parameters.acq.ydet);
- for im=start_image:end_image
- full=full+edf_read(sprintf('1_preprocessing/full/full%04d.edf',im));
- end
- full=full./(end_image-start_image+1);
-
- % ajout sabine%
- close all
- figure
- imshow(full,[]);
- hold on;
- plot(x_direct(i),y_direct(i),'og')
- plot(x_diffr(i),y_diffr(i),'ro')
- drawnow
- disp('click image to close and see next spot (if exists...)')
- k=waitforbuttonpress;
- end
- end
- %~~~~~~~~plot results on full images~~~~~~~~~~~~~~~~~~~~~~
diff --git a/4_grains/gtPredictPositions_3D.m b/4_grains/gtPredictPositions_3D.m
deleted file mode 100644
index 103c2533b4bb5b770feaa0b86b430117d24d4544..0000000000000000000000000000000000000000
--- a/4_grains/gtPredictPositions_3D.m
+++ /dev/null
@@ -1,277 +0,0 @@
-function [images,u_direct,v_direct,u_diffr,v_diffr,Omega,Eta,Theta] = gtPredictPositions_3D(phaseid, normal, hkl, grainCentroid, plot_flag, parameters, codeflag)
-% GTPREDICTPOSITIONS_3D predict positions of spots in images
-% [images,u_direct,v_direct,u_diffr,v_diffr,Omega,Eta,Theta] = gtPredictPositions_3D(phaseid, normal, hkl, grainCentroid, plot_flag, parameters, codeflag)
-% -----------------------------------------------------------------------------------------------------------------------------------------------
-%
-% if parameters.geo includes the neccessary fields for the arbitary
-% geometry, use the FED functions to predict the position of the
-% diffraction spot on the detector.
-% this may be a bit temporary until; FED has been updated
-%
-% take a plane normal, its hkl type (to calc Theta), and a grainCentroid [x,y,z], and returns lists of max images, and x and y
-% positions of the ext and dif spots.
-% plot_flag=1 plots images
-%
-% Version 005 01-08-2013 by LNervo
-% Replaced gtCalculateDist with gtCrystDSpacing and gtCrystTheta
-%
-% Modified on March 08 and March 14, 2012 by PReischig
-% Adapted input for gtCalculateDist.
-%
-% Modified on January 18, 2012 by LNervo
-% Cleaning and formatting
-%
-% Modified on June 28, 2011 by LNervo
-% Replace forward simulation part with the FED code
-% Add a flag to use the old code eve if the fed parameters are in
-% Add fed reordering of Omega
-%
-% Modified on March 08, 2009 by Marcelo, Andy, Wolfgang
-% Works with old coordinate system
-
-acq=parameters.acq;
-geo=parameters.geo;%obsolete...it does not exist anymore
-labgeo = parameters.labgeo;
-
-% do we have all the required fields?
-FEDgeofields{1}='detpos0';
-FEDgeofields{2}='detdiru0';
-FEDgeofields{3}='detdirv0';
-FEDgeofields{4}='beamdir0';
-FEDgeofields{5}='rotdir0';
-FEDgeofields{6}='detucen0';
-FEDgeofields{7}='detvcen0';
-
-% newcode
-if all(isfield(geo, FEDgeofields)) && codeflag
- % use the new functions, allowing arbitary geometry
- sam = gtSampleGeoInSampleSystem(parameters);
- % matrices for computing rotation tensors / arbitary geometry unpacking
- rotcomp = gtFedRotationMatrixComp(geo.rotdir0');
- beamdir0 = geo.beamdir0';
- % detector geometry
- det.uvcen = [geo.detucen0, geo.detvcen0]';
- det.Q = gtFedDetectorProjectionTensor(geo.detdiru0', geo.detdirv0',1,1); %det.Q=gtFedDetectorProjectionTensor(geo.detdiru0, geo.detdirv0)
- det.p = geo.detpos0'/acq.pixelsize; % in pixel
- det.n = cross(geo.detdiru0, geo.detdirv0)';
- % step size
- omstep = 180/acq.nproj;
-
- % given a hkl plane, calculate the interplanar distance.
- Theta = gtCrystTheta( gtCrystDSpacing(hkl, gtCrystHKL2CartesianMatrix(parameters.cryst(phaseid).latticepar) ), parameters.acq.energy);
- sinth = sind(Theta);
-
- % calc Omega angles
- [Omega,pllab,pls] = gtFedPredictOmega(normal',sinth,beamdir0,geo.rotdir0',rotcomp);
-
- if isempty(Omega)
- images = NaN;
- u_direct = NaN;
- v_direct = NaN;
- u_diffr = NaN;
- v_diffr = NaN;
- Theta = NaN;
- Eta = NaN;
- Omega = NaN;
- return;
- end
-
- if ~isempty(Omega)
- % initial S rotation tensor
- S1 = gtFedRotationTensor(Omega(1), rotcomp);
- S2 = gtFedRotationTensor(Omega(2), rotcomp);
- S3 = gtFedRotationTensor(Omega(3), rotcomp);
- S4 = gtFedRotationTensor(Omega(4), rotcomp);
-
- % Diffraction vector:
- d1 = gtFedPredictDiffVec(pllab(:,1), sinth, beamdir0);
- d2 = gtFedPredictDiffVec(pllab(:,2), sinth, beamdir0);
- d3 = gtFedPredictDiffVec(pllab(:,3), sinth, beamdir0);
- d4 = gtFedPredictDiffVec(pllab(:,4), sinth, beamdir0);
-
- d1 = d1/norm(d1); % d1'/norm(d1)
- d2 = d2/norm(d2);
- d3 = d3/norm(d3);
- d4 = d4/norm(d4);
-
- % Following the convention of the mathcing output, the omega value
- % smaller than 180deg (spot "a") is used in the pair data.
- % The two Friedel pairs are 1a-1b and 2a-2b.
-
- if Omega(1) < Omega(3)
- Omega1a = Omega(1);
- s1a = S1;
- s1b = S3;
- Omega1b = Omega(3);
- pl1a = pls(:,1)';
- pl1b = pls(:,3)';
- pllab1a = pllab(:,1)';
- pllab1b = pllab(:,3)';
- d1a = d1';
- d1b = d3';
- else
- Omega1a = Omega(3);
- s1a = S3;
- s1b = S1;
- Omega1b = Omega(1);
- pl1a = pls(:,3)';
- pl1b = pls(:,1)';
- pllab1a = pllab(:,3)';
- pllab1b = pllab(:,1)';
- d1a = d3';
- d1b = d1';
- end
-
- if Omega(2) < Omega(4)
- Omega2a = Omega(2);
- Omega2b = Omega(4);
- s2a = S2;
- s2b = S4;
- pl2a = pls(:,2)';
- pl2b = pls(:,4)';
- pllab2a = pllab(:,2)';
- pllab2b = pllab(:,4)';
- d2a = d2';
- d2b = d4';
- else
- Omega2a = Omega(4);
- Omega2b = Omega(2);
- s2a = S4;
- s2b = S2;
- pl2a = pls(:,4)';
- pl2b = pls(:,2)';
- pllab2a = pllab(:,4)';
- pllab2b = pllab(:,2)';
- d2a = d4';
- d2b = d2';
- end
-
- Omega=[Omega1a;Omega1b;Omega2a;Omega2b];
- Theta=Theta*ones(size(Omega));
-
- % Eta angles
- eta1a = gtGeoEtaFromDiffVec(d1a, labgeo); % takes row vector
- eta1b = gtGeoEtaFromDiffVec(d1b, labgeo);
- eta2a = gtGeoEtaFromDiffVec(d2a, labgeo);
- eta2b = gtGeoEtaFromDiffVec(d2b, labgeo);
-
- Eta=[eta1a;eta1b;eta2a;eta2b];
-
- % u,v coordinates on the detector
- uv1a=gtFedPredictUVW(s1a,d1a',grainCentroid',det.p,det.n,det.Q,det.uvcen,Omega1a,omstep);
- uv1b=gtFedPredictUVW(s1b,d1b',grainCentroid',det.p,det.n,det.Q,det.uvcen,Omega1b,omstep);
- uv2a=gtFedPredictUVW(s2a,d2a',grainCentroid',det.p,det.n,det.Q,det.uvcen,Omega2a,omstep);
- uv2b=gtFedPredictUVW(s2b,d2b',grainCentroid',det.p,det.n,det.Q,det.uvcen,Omega2b,omstep);
-
- images=round([uv1a(3);uv1b(3);uv2a(3);uv2b(3)]); %added round
- u_diffr=round([uv1a(1);uv1b(1);uv2a(1);uv2b(1)]);
- v_diffr=round([uv1a(2);uv1b(2);uv2a(2);uv2b(2)]);
-
- % calc u_dir and v_dir
- uv1a_direct=gtFedPredictUVW(s1a,beamdir0,grainCentroid',det.p,det.n,det.Q,det.uvcen,Omega1a,omstep);
- uv1b_direct=gtFedPredictUVW(s1b,beamdir0,grainCentroid',det.p,det.n,det.Q,det.uvcen,Omega1b,omstep);
- uv2a_direct=gtFedPredictUVW(s2a,beamdir0,grainCentroid',det.p,det.n,det.Q,det.uvcen,Omega2a,omstep);
- uv2b_direct=gtFedPredictUVW(s2b,beamdir0,grainCentroid',det.p,det.n,det.Q,det.uvcen,Omega2b,omstep);
-
- u_direct=round([uv1a_direct(1);uv1b_direct(1);uv2a_direct(1);uv2b_direct(1)]);
- v_direct=round([uv1a_direct(2);uv1b_direct(2);uv2a_direct(2);uv2b_direct(2)]);
- end % end if Omega
-
-else % end arbitrary geometry
-
- % use the standard geometry
- grainx = grainCentroid(1);
- grainy = grainCentroid(2);
- grainz = grainCentroid(3);
-
- % use gtPredictAngles to get Theta, Eta, Omega - gives a pair of angles
- [Theta,Eta,Omega] = gtPredictAngles3D(normal, hkl, parameters);
-
- if ~isempty(Theta)
- %get only Omegas within 0-180 degrees for 180deg data
- if strcmpi(parameters.acq.type,'180degree')
- %disp('180 degree data')
- dummy=find(Omega<0 | Omega>180);
- Theta(dummy)=[];
- Eta(dummy)=[];
- Omega(dummy)=[];
- end
-
- for i=1:length(Omega)
- % lab: grain position in Labsystem
- [labx,laby,labz]=gtTrSamToLab(grainx,grainy,grainz,Omega(i));
- lab=[labx; laby; labz];
- %diffracted beam direction (normalized)
- d=[cosd(2*Theta(i));-sind(2*Theta(i)).*sind(Eta(i));sind(2*Theta(i)).*cosd(Eta(i))];
-
- % intersection of diffraction beam with detector plane: labx+f*dx=detx
- f=(acq.dist/acq.pixelsize+parameters.match.corr_rot_to_det-labx)/d(1);
- % diffraction spot position
- pos=lab+f*d;
-
- u_diffr(i)=acq.rotu+pos(2); % old convention
- v_diffr(i)=acq.ydet/2-pos(3);
-
- u_direct(i)=acq.rotu+lab(2); % old convention
- v_direct(i)=acq.ydet/2-lab(3);
- end
-
- %round to integer values
- u_direct = round(u_direct);
- v_direct = round(v_direct);
- u_diffr = round(u_diffr);
- v_diffr = round(v_diffr);
- %get number of full image
- images = round((Omega/180)*acq.nproj);
-
- else
- images = NaN;
- u_direct = NaN;
- v_direct = NaN;
- u_diffr = NaN;
- v_diffr = NaN;
- Theta = NaN;
- Eta = NaN;
- Omega = NaN;
- end
-
-end % end standard geometry
-
-%~~~~~~~~plot results on full images~~~~~~~~~~~~~~~~~~~~~~
-if plot_flag
-
- for i=1:size(Omega,2)
-
- %convert Omega to image number
- start_image = images(i)-5;
- if start_image<0
- start_image=0;
- end
- end_image = images(i)+5;
- if end_image>(2*acq.nproj)-1;
- end_image=(2*acq.nproj)-1;
- end
- full=zeros(acq.ydet,acq.xdet);
- for im=start_image:end_image
- full=full+edf_read( fullfile('1_preprocessing','full', sprintf('full%04d.edf',im)) );
- end
- full=full./(end_image-start_image+1);
-
- % ajout sabine%
- close all
-
- figure();
- imshow((full),[-10 10]);
- text=sprintf('%d %d %d reflection',hkl(1),hkl(2), hkl(3));
- title(text);
- hold on;
- plot(u_direct(i),v_direct(i),'go');
- plot(u_diffr(i),v_diffr(i),'ro');
-
- disp('click image to close and see next spot (if exists...)')
- k=waitforbuttonpress;
- end
-end
-%~~~~~~~~plot results on full images~~~~~~~~~~~~~~~~~~~~~~
-
-end %end of function
diff --git a/4_grains/gtPredictPositions_Si.m b/4_grains/gtPredictPositions_Si.m
deleted file mode 100644
index 5110aa6e1edb0f009c4cc736633539750d2fd8b5..0000000000000000000000000000000000000000
--- a/4_grains/gtPredictPositions_Si.m
+++ /dev/null
@@ -1,32 +0,0 @@
-%predict positions of spots in images
-
-%_Si - simpler version. relative to "grain" at center of image and on
-%rotation axis, predict image coordinates from eta, theta
-
-%Marcelo, Andy, 30/10
-
-function [x_diffr,y_diffr] = gtPredictPositions_Si(Theta, Eta)
-
-disp('hard coded distance and pixel size in gtPredictPositions_Si')
-distance=11.452
-pixelsize=0.0024
-
-%convert to radians
-Theta=Theta*pi/180;
-Eta=Eta*pi/180;
-
-
-%diffraction spot position
-% radius on the screen (in pixels) with the corrected distance from the sample to the
-% detector (y in sample coordinates) with the omega offset
-radius = tan(2*Theta).*(distance/pixelsize);
-
-x_diffr = 1024 + (radius.*sin(Eta)); % spot position (x position) in the screen
-y_diffr = 1024 - (radius.*cos(Eta)); % vertical coordinate on the screeen
-disp('there might be a problem with this! but I dont know why...')
-
-
-%round to integer values
-x_diffr = round(x_diffr);
-y_diffr = round(y_diffr);
-
diff --git a/4_grains/gtPredictPositions_pair2grain.m b/4_grains/gtPredictPositions_pair2grain.m
deleted file mode 100644
index 83dd9f0b17b62a61283183df379c25b1a05f9ac4..0000000000000000000000000000000000000000
--- a/4_grains/gtPredictPositions_pair2grain.m
+++ /dev/null
@@ -1,119 +0,0 @@
-%predict positions of spots in images
-
-% gtPredictPositions_pair2grain
-%work with difspotpair grain indexing process. Different grain position
-%defination
-
-
-%Marcelo, Andy, 30/10
-
-function [images,x_direct,y_direct,x_diffr,y_diffr] = gtPredictPositions_pair2grain(normal, hkl, grainCentroid, plot_flag, parameters)
-%take a plane normal, its hkl type (to calc Theta), and a grainCentroid [x,y,z], and returns lists of max images, and x and y
-%positions of the ext and dif spots.
-%plot_flag=1 plots images
-
-%load data
-if isempty(parameters)
- load('parameters.mat');
-end
-
-
-%use gtPredictAngles to get Theta, Eta, Omega - note gives a pair of angles
-[Theta,Eta,Omega] = gtPredictAngles(normal, hkl);
-%returns in degree
-
-%get only Omegas within 0-180 degrees for 180deg data
-if parameters.acq.type=='180degree'
- %disp('180 degree data')
- dummy=find(Omega<0 | Omega>180);
- Theta(dummy)=[];
- Eta(dummy)=[];
- Omega(dummy)=[];
-end
-
-
-
-
-%deal with grain position
-[grainCentroid_om] = sfSample_to_setup_cor(grainCentroid, Omega');
-
-%extinction spot position
-x_direct = parameters.acq.rotx + grainCentroid_om(:,2); % y coord of grainpos corresponds to x on detector
-y_direct = parameters.acq.bb(2)-grainCentroid_om(:,3);
-%z in relative to the top to the sample bb, positive z moves spot towards top of image, hence negative
-
-%diffraction spot position
-radius = tand(2*Theta').*((parameters.acq.dist/parameters.acq.pixelsize)-grainCentroid_om(:,1));
-x_diffr = x_direct + (radius.*sind(Eta')); % spot position (x position) in the screen
-y_diffr = y_direct - (radius.*cosd(Eta')); % vertical coordinate on the screeen
-
-%round to integer values
-x_direct = round(x_direct);
-y_direct = round(y_direct);
-x_diffr = round(x_diffr);
-y_diffr = round(y_diffr);
-
-%get full image
-images = round((Omega'/180)*parameters.acq.nproj);
-
-
-
-
- %~~~~~~~~plot results on full images~~~~~~~~~~~~~~~~~~~~~~
- if plot_flag
- for i=1:size(Omega,2)
-
- %convert Omega to image number
- if images(i) >= parameters.acq.nproj
- images(i)=images(i)-parameters.acq.nproj;
- path=sprintf('../%s/',parameters.acq.difB_name);
- else
- path=sprintf('../%s/',parameters.acq.difA_name);
- end
-
- start_image = images(i)-3;
- if start_image<0
- start_image=0;
- end
- end_image = images(i)+3;
- if end_image>parameters.acq.nproj-1;
- end_image=parameters.acq.nproj-1;
- end
-
- start_image
- end_image
-
- full=zeros(parameters.acq.xdet,parameters.acq.ydet);
- for im=start_image:end_image
- full=full+edf_read(sprintf('%s1_preprocessing/full/full%04d.edf',path,im));
- end
- full=full./(end_image-start_image+1);
-
- figure
- imshow(full,[]);
- hold on;
- plot(x_direct(i),y_direct(i),'og')
- plot(x_diffr(i),y_diffr(i),'ro')
- drawnow
- disp('click image to close and see next spot (if exists...)')
- k=waitforbuttonpress;
- end
- end
- %~~~~~~~~plot results on full images~~~~~~~~~~~~~~~~~~~~~~
-
-
-
- function [insetup_vec]=sfSample_to_setup_cor(insample_vec,om)
- % om = omega angle in degrees
- %insetup_vec=([cosd(om) sind(om) 0; -sind(om) cosd(om) 0; 0 0
- %1]*insample_vec')' ;
- insetup_vec(:,1)= cosd(om)*insample_vec(1)+sind(om)*insample_vec(2) ;
- insetup_vec(:,2)= -sind(om)*insample_vec(1)+cosd(om)*insample_vec(2) ;
- insetup_vec(:,3)= insample_vec(3) ;
- end
-
-
-end
-
-
-
diff --git a/4_grains/gtPredictPositions_pair2grain360.m b/4_grains/gtPredictPositions_pair2grain360.m
deleted file mode 100644
index 405209885ffa4c00122ba9b1c217d6a27a3de57d..0000000000000000000000000000000000000000
--- a/4_grains/gtPredictPositions_pair2grain360.m
+++ /dev/null
@@ -1,77 +0,0 @@
-%predict positions of spots in images
-
-% gtPredictPositions_pair2grain
-%work with difspotpair grain indexing process. Different grain position
-%defination
-%360 no a / b half scans
-
-%Marcelo, Andy, 30/10
-
-function [images,x_direct,y_direct,x_diffr,y_diffr] = gtPredictPositions_pair2grain360(normal, hkl, grainCentroid, plot_flag, parameters)
-%take a plane normal, its hkl type (to calc Theta), and a grainCentroid [x,y,z], and returns lists of max images, and x and y
-%positions of the ext and dif spots.
-%plot_flag=1 plots images
-
-%load data
-if isempty(parameters)
- load('parameters.mat');
-end
-
-
-%use gtPredictAngles to get Theta, Eta, Omega - note gives a pair of angles
-[Theta,Eta,Omega] = gtPredictAngles(normal, hkl);
-%returns in degree
-
-%get only Omegas within 0-180 degrees for 180deg data
-if parameters.acq.type=='180degree'
- %disp('180 degree data')
- dummy=find(Omega<0 | Omega>180);
- Theta(dummy)=[];
- Eta(dummy)=[];
- Omega(dummy)=[];
-end
-
-
-
-
-%deal with grain position
-[grainCentroid_om] = sfSample_to_setup_cor(grainCentroid, Omega');
-
-%extinction spot position
-x_direct = parameters.acq.rotx + grainCentroid_om(:,2); % y coord of grainpos corresponds to x on detector
-y_direct = parameters.acq.bb(2)-grainCentroid_om(:,3);
-%z in relative to the top to the sample bb, positive z moves spot towards top of image, hence negative
-
-%diffraction spot position
-radius = tand(2*Theta').*((parameters.acq.dist/parameters.acq.pixelsize)-grainCentroid_om(:,1));
-x_diffr = x_direct + (radius.*sind(Eta')); % spot position (x position) in the screen
-y_diffr = y_direct - (radius.*cosd(Eta')); % vertical coordinate on the screeen
-
-%round to integer values
-x_direct = round(x_direct);
-y_direct = round(y_direct);
-x_diffr = round(x_diffr);
-y_diffr = round(y_diffr);
-
-%get full image
-images = round((Omega'/180)*parameters.acq.nproj);
-
-
-
-
-
-
- function [insetup_vec]=sfSample_to_setup_cor(insample_vec,om)
- % om = omega angle in degrees
- %insetup_vec=([cosd(om) sind(om) 0; -sind(om) cosd(om) 0; 0 0
- %1]*insample_vec')' ;
- insetup_vec(:,1)= cosd(om)*insample_vec(1)+sind(om)*insample_vec(2) ;
- insetup_vec(:,2)= -sind(om)*insample_vec(1)+cosd(om)*insample_vec(2) ;
- insetup_vec(:,3)= insample_vec(3) ;
- end
-
-
-end
-
-
-
diff --git a/4_grains/gtReadExtRoi_360.m b/4_grains/gtReadExtRoi_360.m
deleted file mode 100644
index 82d142146b9941e1057927ba713fcbad54205711..0000000000000000000000000000000000000000
--- a/4_grains/gtReadExtRoi_360.m
+++ /dev/null
@@ -1,76 +0,0 @@
-function [extstack,sino,centerslice,bb,index]=gtReadExtRoi_360(grain,parameters);
-
-np=length(grain.difspots);
-
-for j=1:np
- % get lab system grain COM Coordinates from grain structure (produced by PEter's
- % Indexing)
-% grain.center(1)=grain.ycenter-parameters.acq.bb(3)/2; % not very elegant: transform back for compatibility with calls of gtReadextroi outside of gtINWriteGrainFolder...
-% grain.center(2)=grain.xcenter-parameters.acq.bb(3)/2;
-% grain.center(3)=parameters.acq.ydet/2-parameters.acq.bb(2)-grain.zcenter;
-
- [xl,yl,zl]=gtTrSamToLab(grain.center(1),grain.center(2),grain.center(3),grain.omega(j));
- xim = yl+ parameters.acq.bb(3)/2; % grain COM coordinates in the direct beam images (assumes center of rotation in the center of the beam)
- yim = parameters.acq.ydet/2-zl-parameters.acq.bb(2);
-
-
- % get Boundingbox sizes and integration range from all difspots
- mymcmd=sprintf('select BoundingBoxXSize, BoundingBoxYSize, ExtStartImage, ExtEndImage, StartImage, EndImage, Integral from %sdifspot where difspotID=%d',parameters.acq.name,grain.difspots(j));
- [bbxs,bbys,startndx,endndx,difstartndx,difendndx,int]=mym(mymcmd);
- difbb(j,1)=max(1,round(xim-bbxs/2));
- difbb(j,2)=max(1,round(yim-bbys/2));
-
- bbxs=round(bbxs*1.1); %add 10% margin to bb
- bbys=round(bbys*1.1);
- integral(j)=int;
- ndx(j,:)=[startndx,endndx];
- difndx(j,:)=[difstartndx, difendndx];
- range(j)=endndx-startndx+1;
- difrange(j)=difendndx-difstartndx+1;
- bb(j,1)=max(1,round(xim-bbxs/2)); % horizontal start of bb in Ext image
- hend=min(parameters.acq.bb(3),bb(j,1)+bbxs-1); % limit horizontal extend to size of Ext image...
- bb(j,3)=hend-bb(j,1)+1;
-
- bb(j,4)=bbys;
-
-
-end
-
-% vertical bounds are identical for all extspots
-vertsize=mean(bb(:,4)); % first estimate of vertical size from average vertical size in difspots
-vertbeg=max(1,round(yim-vertsize/2));
-bb(:,2)=repmat(vertbeg,np,1);
-vertend=min(parameters.acq.bb(4),round(yim+vertsize/2));
-vertsize=vertend-vertbeg; % update vertical size of bb, taking into account finite size of ext image
-bb(:,4)=repmat(vertsize,np,1);
-centerslice=round(vertbeg+0.5*vertsize);
-
-
-
-% now read and sum under bb mask
-extstack=zeros(parameters.acq.bb(4),parameters.acq.bb(3),np);
-difstack=zeros(parameters.acq.bb(4),parameters.acq.bb(3),np);
-
-tmp=zeros(parameters.acq.bb(4),parameters.acq.bb(3));
-meanint=mean(integral);
-
-for j=1:np
- % get summed diffraction spot
- %dif=gtGetSummedDifspot(grain.difspots(j),parameters,1);
- %difstack(:,:,j)=gtPlaceSubImage(dif,tmp,difbb(j,1),difbb(j,2));
- %difstack(:,:,j)=meanint*difstack(:,:,j)/integral(j);
- im=zeros(bb(j,4),bb(j,3));
- for k=1:difrange(j)
- im = im+edf_read(sprintf('%s/1_preprocessing/ext/ext%04d.edf',parameters.acq.dir,difndx(j,1)+k-1),bb(j,:));
- end
-
- % place back into ext image and normalize with meanintensity of difspots
-
- extstack(:,:,j)=gtPlaceSubImage(im,tmp,bb(j,1),bb(j,2));
- extstack(:,:,j) = meanint*extstack(:,:,j)/integral(j);
-end
-
-index=ones(1,np); % think of some clever rejection of outliers here....
-
-sino=squeeze(extstack(centerslice,:,:));
-
diff --git a/4_grains/gtReadExtRoi_360wl.m b/4_grains/gtReadExtRoi_360wl.m
deleted file mode 100644
index c44f8175d3131c59cac6894466f6a9c0d8cc3b86..0000000000000000000000000000000000000000
--- a/4_grains/gtReadExtRoi_360wl.m
+++ /dev/null
@@ -1,95 +0,0 @@
-function [ext,sino,centerslice,bb,index]=gtReadExtRoi_360(grain,parameters);
-
-np=length(grain.difspots);
-vsizemax=0;
-hsizemax=0;
-avgint=grain.stat.intmean;
-
-for j=1:np
-
- % get lab system grain COM Coordinates
- [xl,yl,zl]=gtTrSamToLab(grain.center(1),grain.center(2),grain.center(3),grain.omega(j));
- uim = yl+ parameters.acq.bb(3)/2+0.5; % grain COM coordinates in the direct beam images (assumes looking into the beam, center of rotation in the center of the beam)
- vim = parameters.acq.ydet/2-zl-parameters.acq.bb(2)+0.5;
-
-
- % get Boundingbox sizes and integration range from all difspots
- mymcmd=sprintf('select BoundingBoxXSize, BoundingBoxYSize, ExtStartImage, ExtEndImage, StartImage, EndImage, Integral,Area from %sdifspot where difspotID=%d',parameters.acq.name,grain.difspots(j));
- [bbxs,bbys,startndx,endndx,difstartndx,difendndx,int,area]=mym(mymcmd);
-
-
- bbxs=round(bbxs*1.25); %add 25% margin to bb
- bbys=round(bbys*1.25);
-
- vsizemax=max(vsizemax,bbys);
- hsizemax=max(hsizemax,bbxs);
-
- integral(j)=int;
- ndx(j,:)=[startndx,endndx];
- difndx(j,:)=[difstartndx, difendndx];
- range(j)=endndx-startndx+1;
- difrange(j)=difendndx-difstartndx+1;
- bb(j,1)=max(1,round(uim-bbxs/2)); % horizontal start of bb in Ext image
- hend=min(parameters.acq.bb(3),bb(j,1)+bbxs-1); % limit horizontal extend to size of Ext image...
- bb(j,3)=hend-bb(j,1)+1;
-
- bb(j,4)=bbys;
- mym(gtDBUpdate(sprintf('%sdifspot', parameters.acq.name),'ExtBBXorigin',bb(j,1),'ExtBBYorigin',bb(j,2),'ExtBBXsize',bb(j,3),'ExtBBYsize',bb(j,4)));
-
-end
-
-% vertical bounds should be identical for all extspots
-vertsize=mean(bb(:,4)); % first estimate of vertical size from average vertical size in difspots
-vertbeg=max(1,round(vim-vertsize/2));
-bb(:,2)=repmat(vertbeg,np,1);
-vertend=min(parameters.acq.bb(4),round(vim+vertsize/2));
-vertsize=vertend-vertbeg+1; % update vertical size of bb, taking into account finite size of ext image
-bb(:,4)=repmat(vertsize,np,1);
-centerslice=round(vertbeg+0.5*vertsize);
-
-
-
-% now read and sum under bb mask
-extstack=zeros(parameters.acq.bb(4),parameters.acq.bb(3),np);
-%difstack=zeros(parameters.acq.bb(4),parameters.acq.bb(3),np);
-
-tmp=zeros(parameters.acq.bb(4),parameters.acq.bb(3));
-tmp2=zeros(vsizemax,hsizemax);
-ext=zeros(vsizemax,hsizemax,1,np);
-meanint=mean(integral);
-
-for j=1:np
- % get summed diffraction spot
- %dif=gtGetSummedDifspot(grain.difspots(j),parameters,1);
- %difstack(:,:,j)=gtPlaceSubImage(dif,tmp,difbb(j,1),difbb(j,2));
- %difstack(:,:,j)=meanint*difstack(:,:,j)/integral(j);
- im=zeros(bb(j,4),bb(j,3));
- for k=1:range(j)
- im = im+edf_read(sprintf('%s/1_preprocessing/ext/ext%04d.edf',parameters.acq.dir,ndx(j,1)+k-1),bb(j,:));
- end
-
-
- % place back into ext image and normalize with meanintensity of difspots
-
- im=im.*(im>0);
- im=avgint*im/sum(im(:));
-
-
- hstart=floor((hsizemax-size(im,2))/2+1);
- vstart=floor((vsizemax-size(im,1))/2+1);
- ext(:,:,1,j)=gtPlaceSubImage(im,tmp2,hstart,vstart);
-
-
- extstack(:,:,j)=gtPlaceSubImage(im,tmp,bb(j,1),bb(j,2));
-
- %extstack(:,:,j) = meanint*extstack(:,:,j)/sum(sum(extstack(:,:,j).*(extstack(:,:,j)>0)));
- name=sprintf('4_grains/grain%d_/extspot%d',grain.id,j);
- sdt_write(name,fliplr(extstack(:,:,j)),'float32'); % flipping because we look into the beam....
-
-end
-
-
-% only consider strong spots
-index=ones(np,1);
-sino=squeeze(extstack(centerslice,:,:));
-
diff --git a/4_grains/gtReadSpots.m b/4_grains/gtReadSpots.m
deleted file mode 100644
index a56da293492a9fc377ec422d129eab7920b8f51e..0000000000000000000000000000000000000000
--- a/4_grains/gtReadSpots.m
+++ /dev/null
@@ -1,79 +0,0 @@
-function [stack,sino,z_of_sino] = gtReadSpots(spot_id, struct_ids)
-% function [stack,sino,z_of_sino,varargout]=read_spots(struct_ids)
-% produce sinogram of all spots in struct_ids vector at the mid-height
-% of the seed spot (spot_id).
-% optional output argument orig holds original projection (Max)images for verification
-%could add orig, not currently supported
-
-%database enable!
-
-parameters=[];
-if isempty(parameters)
- load('parameters.mat');
-end
-
-
-
-
-stack=zeros(parameters.acq.bb(4),parameters.acq.bb(3),length(struct_ids));
-
-z_of_sino = mym(sprintf(...
- 'select ceil(Yorigin + (Ysize/2)) from %sbb inner join %sbboxes on %sbb.bbID = %sbboxes.bboxID where %sbb.extspotID = %d',...
- parameters.acq.name,...
- parameters.acq.name,...
- parameters.acq.name,...
- parameters.acq.name,...
- parameters.acq.name,...
- spot_id));
-
-%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-% consider 15 pixel shift - s5_dct5_ only
-%warning('15 pixel shift (s5_dct5_ only) active in gtReadSpots')
-%spot_Omega = (180/pi)*gtGetOmega(spot_id);
-%if spot_Omega>120
-%z_of_sino = z_of_sino-15;
-%end
-
-Omega = (180/pi)*gtGetOmega(struct_ids);%convert to degrees
-%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-%read images to produce sinogram at position z_of_sino
-for i=1:length(struct_ids)
-
- [BoundingBox(1),BoundingBox(2)] = mym(sprintf('select Xorigin,Yorigin from %sbb inner join %sbboxes on %sbb.bbID=%sbboxes.bboxID where %sbb.extspotID = %d',...
- parameters.acq.name,...
- parameters.acq.name,...
- parameters.acq.name,...
- parameters.acq.name,...
- parameters.acq.name,...
- struct_ids(i)));
-
- im = gtGetSummedExtSpot(struct_ids(i),1);% - 1 - already connected
- warning('bodge - change zero origin bounding boxes in gtReadSpots - fix in autofind')
- if BoundingBox(1)==0
- BoundingBox(1)=1;
- end
- if BoundingBox(2)==0
- BoundingBox(2)=1;
- end
-
- %place in stack
- stack(:,:,i) = gtPlaceSubImage(im, zeros(parameters.acq.bb(4), parameters.acq.bb(3)),ceil(BoundingBox(1)), ceil(BoundingBox(2)));
-
- %!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- %15 pixel offset - apply here
- % warning('15 pixel shift (s5_dct5_ only) active in gtReadSpots')
- % if(Omega(i)>120);%these need to be shifted by -15
- % stack(1:end-15,:,i) = stack(16:end,:,i);
- % stack(end-14:end,:,i) = zeros(15,parameters.acq.bb(3));%write zeros over the last 15 pixels
- % end
- %!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- stack(:,:,i)= 100*stack(:,:,i)/sum(sum(stack(:,:,i))); % think about some proper scaling here !
-
- %read sino from stack
- sino(:,i) = stack(z_of_sino, :, i);
-
-end
-
-end
diff --git a/4_grains/gtReadSpots_grain.m b/4_grains/gtReadSpots_grain.m
deleted file mode 100644
index d46fab4dddf3370c42dc55389cbc9c781b1d17d7..0000000000000000000000000000000000000000
--- a/4_grains/gtReadSpots_grain.m
+++ /dev/null
@@ -1,96 +0,0 @@
-function [stack,sino,z_of_sino] = gtReadSpots_grain(grainid)
-% function [stack,sino,z_of_sino,varargout]=read_spots(struct_ids)
-% produce sinogram of all spots in struct_ids vector at the mid-height
-% of the seed spot (spot_id).
-% optional output argument orig holds original projection (Max)images for verification
-%could add orig, not currently supported
-
-%can be used to "refresh" stack!
-
-%database enable!
-
-parameters=[];
-if isempty(parameters)
- load('parameters.mat');
-end
-
-grain_data=load(sprintf('4_grains/grain%d_/grain%d_.mat',grainid,grainid))
-struct_ids=grain_data.struct_ids;
-spot_id=struct_ids(1);%the first of the stack...
-
-
-stack=zeros(parameters.acq.bb(4),parameters.acq.bb(3),length(struct_ids));
-
-%if it doesn't already exist, add an index of replaced projections
-if isfield(grain_data,'replaced')
- replaced = logical(grain_data.replaced);
-else
- replaced = logical(zeros(length(struct_ids),1));
- save (sprintf('4_grains/grain%d_/grain%d_.mat',grainid,grainid), 'replaced', '-append')
-end
-
-
-z_of_sino = mym(sprintf(...
- 'select ceil(Yorigin + (Ysize/2)) from %sbb inner join %sbboxes on %sbb.bbID = %sbboxes.bboxID where %sbb.extspotID = %d',...
- parameters.acq.name,...
- parameters.acq.name,...
- parameters.acq.name,...
- parameters.acq.name,...
- parameters.acq.name,...
- spot_id));
-
-
-
-%read images to produce sinogram at position z_of_sino
-for i=1:length(struct_ids)
-
- if replaced(i)==0
- [BoundingBox(1),BoundingBox(2)] = mym(sprintf('select Xorigin,Yorigin from %sbb inner join %sbboxes on %sbb.bbID=%sbboxes.bboxID where %sbb.extspotID = %d',...
- parameters.acq.name,...
- parameters.acq.name,...
- parameters.acq.name,...
- parameters.acq.name,...
- parameters.acq.name,...
- struct_ids(i)));
-
- im = gtGetSummedExtSpot(struct_ids(i));
- %place in stack
- stack(:,:,i) = gtPlaceSubImage(im, zeros(parameters.acq.bb(4), parameters.acq.bb(3)),ceil(BoundingBox(1)), ceil(BoundingBox(2)));
-
- stack(:,:,i)= 100*stack(:,:,i)/sum(sum(stack(:,:,i))); % think about some proper scaling here !
-
- %read sino from stack
- sino(:,i) = stack(z_of_sino, :, i);
-
- elseif replaced(i)==1
- [BoundingBox(1),BoundingBox(2)] = mym(sprintf('select Xorigin,Yorigin from %sbboxes inner join %sextspot on SearchbbID=bboxID where %sbboxes.extspotID = %d',...
- parameters.acq.name,...
- parameters.acq.name,...
- parameters.acq.name,...
- struct_ids(i)));
-
- warning('bodge because some SearchBoundingBoxes have zero origins')
- if BoundingBox(1)==0
- BoundingBox(1)=1;
- end
- if BoundingBox(2)==0
- BoundingBox(2)=1;
- end
-
- im=edf_read(sprintf('2_difspot/difspot%05d.edf',struct_ids(i)));
- disp('Using transpose of summed dif image')
- im=transpose(im);
- disp('Transposing difspot')
-
- %place in stack
- stack(:,:,i) = gtPlaceSubImage(im, zeros(parameters.acq.bb(4), parameters.acq.bb(3)),ceil(BoundingBox(1)), ceil(BoundingBox(2)));
- stack(:,:,i)= 100*stack(:,:,i)/sum(sum(stack(:,:,i))); % think about some proper scaling here !
- %read sino from stack
- sino(:,i) = stack(z_of_sino,:,i);
-
-end
-
-end
-
-
-save(sprintf('4_grains/grain%d_/grain%d_.mat',grainid,grainid),'stack','-append')
\ No newline at end of file
diff --git a/4_grains/gtReadSpots_pair.m b/4_grains/gtReadSpots_pair.m
deleted file mode 100644
index 6ec3821a8c3997a473cd183244d6861e937d9ea8..0000000000000000000000000000000000000000
--- a/4_grains/gtReadSpots_pair.m
+++ /dev/null
@@ -1,88 +0,0 @@
-function [stack,sino,z_of_sino] = gtReadSpots_pair(spot_id, struct_ids,pair,pair_vector,replaced)
-% function [stack,sino,z_of_sino,varargout]=read_spots(struct_ids)
-% produce sinogram of all spots in struct_ids vector at the mid-height
-% of the seed spot (spot_id).
-%pair, pair_vector (1 for currect dataset, 2 for pair dataset)
-%replaced: 0=use extspots, 1=replace all with difspots (like WriteStack_replaced)
-
-%database enable!
-
-parameters=[];
-if isempty(parameters)
- load('parameters.mat');
-end
-
-
-
-stack=zeros(parameters.acq.bb(4),parameters.acq.bb(3),length(struct_ids));
-
-%get the z position for the sinogram
-if pair == 2
- name=parameters.acq.pair_name;
-elseif pair ==1
- name=parameters.acq.name;
-end
-z_of_sino = mym(sprintf(...
- 'select ceil(Yorigin + (Ysize/2)) from %sbb inner join %sbboxes on %sbb.bbID = %sbboxes.bboxID where %sbb.extspotID = %d',...
- name,name,name,name,name,...
- spot_id));
-
-
-Omega = (180/pi)*gtGetOmega_pair(struct_ids,pair_vector, parameters);%convert to degrees
-
-
-%read images to produce sinogram at position z_of_sino
-for i=1:length(struct_ids)
-
- if pair_vector(i) == 2
- name=parameters.acq.pair_name;
-elseif pair_vector(i) ==1
- name=parameters.acq.name;
- end
-
-
- if ~replaced
-
- [BoundingBox(1),BoundingBox(2)] = mym(sprintf('select Xorigin,Yorigin from %sbb inner join %sbboxes on %sbb.bbID=%sbboxes.bboxID where %sbb.extspotID = %d',...
- name,name,name,name,name,...
- struct_ids(i)));
-
- im = gtGetSummedExtSpot_pair(struct_ids(i), pair_vector(i) ,1);% - 1 - already connected
- warning('bodge - change zero origin bounding boxes in gtReadSpots - fix in autofind')
- if BoundingBox(1)==0
- BoundingBox(1)=1;
- end
- if BoundingBox(2)==0
- BoundingBox(2)=1;
- end
-
- %place in stack
- stack(:,:,i) = gtPlaceSubImage(im, zeros(parameters.acq.bb(4), parameters.acq.bb(3)),ceil(BoundingBox(1)), ceil(BoundingBox(2)));
- %scaling
- stack(:,:,i)= 100*stack(:,:,i)/sum(sum(stack(:,:,i))); % think about some proper scaling here !
-
- else %if replacing extspots with difspots
-
- BoundingBox=gtGetBBProps_pair(struct_ids(i), pair_vector(i));
- im = gtGetSummedDifSpot_pair(struct_ids(i), pair_vector(i));
- warning('bodge - change zero origin bounding boxes in gtReadSpots - fix in autofind')
- if BoundingBox(1)==0
- BoundingBox(1)=1;
- end
- if BoundingBox(2)==0
- BoundingBox(2)=1;
- end
-
- %place in stack
- stack(:,:,i) = gtPlaceSubImage(im, zeros(parameters.acq.bb(4), parameters.acq.bb(3)),ceil(BoundingBox(1)), ceil(BoundingBox(2)));
- %scaling
- stack(:,:,i)= 100*stack(:,:,i)/sum(sum(stack(:,:,i))); % think about some proper scaling here !
-
- end
-
- %read sino from stack
- sino(:,i) = stack(z_of_sino, :, i);
-
-end
-
-end
diff --git a/4_grains/gtRefMultiplicity.m b/4_grains/gtRefMultiplicity.m
deleted file mode 100644
index 8ac7b03ca5af2f9c0f14d8d1a50b4f1ea6af04f3..0000000000000000000000000000000000000000
--- a/4_grains/gtRefMultiplicity.m
+++ /dev/null
@@ -1,11 +0,0 @@
-function mult=gtRefMultiplicity(spacegroup)
-% <<<< NOT USED >>>>
-
-ref=parameters.cryst.hklsp;
-
-for i=1:size(ref,1)
- mult(i)=size(gtCrystSignedHKLs(ref(i,:), gtCrystGetSymmetryOperators([], spacegroup) ),1) ;
-end
-
-end
-
diff --git a/4_grains/gtThetaOfPairs.m b/4_grains/gtThetaOfPairs.m
deleted file mode 100644
index 1af56e28d5f7f874c7f81c929b5acd1074ebe76f..0000000000000000000000000000000000000000
--- a/4_grains/gtThetaOfPairs.m
+++ /dev/null
@@ -1,46 +0,0 @@
-%
-% function theta=gtThetaOfPairs(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 theta=gtThetaOfPairs(scXA,scYA,scXB,scYB,rottodet,tiltY,tiltZ,rotx,sorZ)
-
-% 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) ;
-%
-% theta=0.5*acosd(Dv1./sqrt(Dv1.^2+Dv2.^2+Dv3.^2)) ;
-
-Dv=gtDiffVecInLab(scXA,scYA,scXB,scYB,rottodet,tiltY,tiltZ,rotx,sorZ);
-
-theta=0.5*acosd(Dv(:,1)./sqrt(Dv(:,1).^2+Dv(:,2).^2+Dv(:,3).^2));
-
-end
diff --git a/4_grains/gt_simulate_fcc.m b/4_grains/gt_simulate_fcc.m
deleted file mode 100644
index 914c20884b47851af576609da32981175c35cc3a..0000000000000000000000000000000000000000
--- a/4_grains/gt_simulate_fcc.m
+++ /dev/null
@@ -1,75 +0,0 @@
-
-
-
-clear ttheta ind d;
-
-maxind=10;
-d0= 4.04;
-energy=15;
-minangle=70;
-maxangle=110;
-
-
-lambda=12.39/energy;
-
-i=1;
-
-for h=0:2:maxind
- for k=0:2:maxind
- for l=0:2:maxind
- hkl(i,:)=[h,k,l];
-
- d(i)=d0/sqrt(h.^2+k.^2+l.^2);
-
- x=lambda/2/d(i);
- if x<1
-
- ttheta(i)=2*asind(lambda/2/d(i));
- end
- i=i+1;
- end
- end
-end
-
-
-for h=1:2:maxind
- for k=1:2:maxind
- for l=1:2:maxind
- hkl(i,:)=[h,k,l];
-
- d(i)=d0/sqrt(h.^2+k.^2+l.^2);
-
- x=lambda/2/d(i);
- if x<1
-
- ttheta(i)=2*asind(lambda/2/d(i));
- end
- i=i+1;
- end
- end
-end
-
-
-
-
-ind=find(ttheta>minangle &ttheta<maxangle);
-
-titi=ttheta(ind);
-titi=unique(titi);
-bar(titi)
-shg
-
-keyboard
-hkls=hkl(ind,:)
-
-
-
-
-
-
-
-
-
-
-
-
\ No newline at end of file
diff --git a/4_grains/look_at_possibles.m b/4_grains/look_at_possibles.m
deleted file mode 100644
index 4f062b09c5ac437cb867032e017895da35770ad9..0000000000000000000000000000000000000000
--- a/4_grains/look_at_possibles.m
+++ /dev/null
@@ -1,61 +0,0 @@
-function spot_viewer(struct_ids, pairs, theta, eta, hkl, parameters)
-%show the relevant full images, marking the dif and ext spots
-%theta, eta, hkl just for comments, not used by the function
-
-
-%load data
-if isempty(parameters)
- load('parameters.mat');
-end
-
-for i=1:length(struct_ids)
-
- if pairs(i)==1
- name=parameters.acq.name
- else
- name=parameters.acq.pair_name
- end
-
- [im,x,y]=mym(sprintf('select MaxImage,CentroidX,CentroidY from %sdifspot where difspotID=%d',name,struct_ids(i)));
- if i==1
- full=edf_read(sprintf('../%s/1_preprocessing/full/full%04d.edf',name,im));
-
- imshow(full,[])
- hold on
- end
- switch hkl(i,1)
- case 1
- plot(x,y,'x')
- case 0
- plot(x,y,'o')
- case 2
- plot(x,y,'xr')
- case 3
- plot(x,y,'or')
- end
-
-
- %plot the extspot position that the difspot has been linked to by
- %autofind
-[bb,centroid]=gtGetBBProps_pair(struct_ids(i), pairs(i))
-
-plot(parameters.acq.bb(1)+bb(1),parameters.acq.bb(2)+bb(2),'ro')
-plot(parameters.acq.bb(1)+bb(1)+bb(4),parameters.acq.bb(2)+bb(2)+bb(4),'ro')
-
- switch hkl(i,1)
- case 1
- plot(parameters.acq.bb(1)+centroid(1),parameters.acq.bb(2)+centroid(2),'x')
- case 0
- plot(parameters.acq.bb(1)+centroid(1),parameters.acq.bb(2)+centroid(2),'o')
- case 2
- plot(parameters.acq.bb(1)+centroid(1),parameters.acq.bb(2)+centroid(2),'xr')
- case 3
- plot(parameters.acq.bb(1)+centroid(1),parameters.acq.bb(2)+centroid(2),'or')
- end
-
-
-disp(sprintf('should be theta=%f, eta=%f, (%d %d %d)',theta(i),eta(i),hkl(i,1:3)))
-
- % k=waitforbuttonpress;
- %close
-end
\ No newline at end of file
diff --git a/4_grains/look_at_possibles2.m b/4_grains/look_at_possibles2.m
deleted file mode 100644
index dfc9c819025b4ed8a4d651c788564b3304b10f73..0000000000000000000000000000000000000000
--- a/4_grains/look_at_possibles2.m
+++ /dev/null
@@ -1,38 +0,0 @@
-function look_at_possibles2(grainid)
-%version of look_at_possibles to look at the projections making up a grain,
-%showing the full images with the dif and ext spots marked
-
-
-%load data
-parameters=[];
-if isempty(parameters)
- load('parameters.mat');
-end
-
-grain_data=load(sprintf('4_grains/grain%d_/grain%d_.mat',grainid,grainid));
-struct_ids = grain_data.struct_ids;
-
-
-for i=1:length(struct_ids)
-
- [im1,im2,im,x,y]=mym(sprintf('select ExtStartImage,ExtEndImage,MaxImage,CentroidX,CentroidY from s5_dct5_difspot where difspotID=%d',struct_ids(i)));
-
- full=zeros(2048);
- for im=im1:im2
- full=full+edf_read(sprintf('1_preprocessed/full/full%04d.edf',im));
- end
-
-
- imshow(full,[])
- hold on
- plot(x,y,'x')
-
- %plot the extspot position that the difspot has been linked to by
- %autofind
- [x,y]=mym(sprintf('select CentroidX,CentroidY from s5_dct5_extspot where extspotID=%d',struct_ids(i)));
- plot(parameters.acq.bb(1)+x,parameters.acq.bb(2)+y,'ro')
-
-
- k=waitforbuttonpress;
- close
-end
\ No newline at end of file
diff --git a/4_grains/plot_rodrigues.m b/4_grains/plot_rodrigues.m
deleted file mode 100644
index 4d9098883011064634d1c853348fc0e52ae7df5d..0000000000000000000000000000000000000000
--- a/4_grains/plot_rodrigues.m
+++ /dev/null
@@ -1,116 +0,0 @@
-function max_coords=plot_rodrigues(plane_normal,hkl,spacegroup)
-
-%plane_normal = list of plane normals in sample coords
-%hkl a vector of reflection type eg [1 1 1; 0 0 2; 2 2 0 etc]
-
-save_plot=[];
-figure;
-hold on;
-
-for j=1:size(plane_normal,1)
-
-%________calculate h vector___________
-%from eq (3.7) in Poulsen
-%B is identity matrix in cubic material
-%gtnewGetReflections returns matrix of all reflections related by symmetry
-%to hkl, using Riso GetCubicSymOp.m
-
-hkl_reflections = gtCrystSignedHKLs(hkl(j,:), gtCrystGetSymmetryOperators([],spacegroup));
-
-for i=1:length(hkl_reflections)
-h = hkl_reflections(i,:);
-h = h./(sqrt(h*h'));
-
-%________calculate y vector___________
-% y is the normalised scattering vector in sample system
-% == plane normal
-y = plane_normal(j,:);
-
-
-%_____calculate rodrigues space vector_____
-%formulae (3.16), (3.17) in Poulsen
-
-r_zero = cross(h,y)./(1+(dot(h,y)));
-r_direction = (h + y)./(1+(dot(h,y)));
-
-%max line length needed ~1.5, say +-1.5 for saftey
-t=(-1.5:0.01:1.5)';
-r_plot = repmat(r_zero,length(t),1) + (t*r_direction);
-
-
-
-%remove everything outside of a fundamental region (assume cube)
-lim=sqrt(2)-1;
-for i=1:3
-dummy=find(r_plot(:,i)<-lim | r_plot(:,i)>lim);
-r_plot(dummy,:)=[];
-end
-
-%collect all plot data for analysis
-save_plot=[save_plot; r_plot];
-
-
-
-%____plot______
- if hkl(j,:)*hkl(j,:)'==[1 1 1]*[1 1 1]' %(111) - red
- plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'r');
-
- elseif hkl(j,:)*hkl(j,:)'==[0 0 2]*[0 0 2]' %(002) - blue
- plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'b');
-
- elseif hkl(j,:)*hkl(j,:)'==[2 2 0]*[2 2 0]' %(220) - green
- plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'g');
-
- elseif hkl(j,:)*hkl(j,:)'==[3 1 1]*[3 1 1]' %(311) - black
- plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'k');
-
- end
-
-
-end%loop through all hkl variants
-
-%uncomment this to see line by line in rodrigues space.
-%k=waitforbuttonpress %can watch as each is plotted in turn
-
-end%loop through all normals
-
-axis equal
-
-
-%find intersections
-
-inc=0.04;%scan space in these steps (size of the sampling element)
-max_score=0;
-max_coords=[0 0 0];
-score=[];
-count=1;
-for i=-lim:inc:lim
-for j=-lim:inc:lim
-for k=-lim:inc:lim
-
-dummy = find(save_plot(:,1)>i & save_plot(:,1)<(i+inc) &save_plot(:,2)>j & save_plot(:,2)<(j+inc) &save_plot(:,3)>k & save_plot(:,3)<(k+inc));
-score(count)=length(dummy);
-if (score(count)>max_score)
- max_score=score;
- max_coords=[i j k];
-end
-
-count=count+1;
-end
-end
-end
-
-max_score;
-max_coords;
-element = inc*[0 0 0; 0 0 1; 0 1 0; 0 1 1; 1 0 0; 1 0 1; 1 1 0; 1 1 1];
-max_element = repmat(max_coords,8,1)+element;
-max_coords = max_coords+0.5*[inc inc inc];%this is the approximate solution
-
-plot3(max_element(:,1),max_element(:,2),max_element(:,3),'go')
-plot3(max_coords(:,1),max_coords(:,2),max_coords(:,3),'co')
-
-
-
-
-
-
diff --git a/4_grains/plot_rodrigues_consistancy.m b/4_grains/plot_rodrigues_consistancy.m
deleted file mode 100644
index 143c426f54fd5a17fbb900cae503dd07b46d0c3e..0000000000000000000000000000000000000000
--- a/4_grains/plot_rodrigues_consistancy.m
+++ /dev/null
@@ -1,213 +0,0 @@
-function [pof, good_plnorms, dsqr]=plot_rodrigues_consistancy(plane_normal,hkl,spacegroup,plot_flag)
-%plane_normal = list of plane normals in sample coords
-%hkl a vector of reflection type eg [1 1 1; 0 0 2; 2 2 0 etc]
-%plot_flag: 1 to plot figure, 0 to not plot anything
-%add a consistancy check - which plane normals don't join at the
-%intersection?
-
-inc=0.1; % scan space in these steps (size of the sampling element) for coarse solution
-inc_ext=1.1*inc; % extended space to determine precise point of intersection; (should be larger than 0.87*inc)
-consistancy_test=0.2; %select those lines within distance consistancy_test*inc of the intersection (was 0.2)
-
-save_plot=[];
-
-if plot_flag
-figure;
-hold on;
-end
-
-lines.cor=[]; %[...; r_zero r_direction; ...] equations
-lines.plnorm=[]; %[...; j; ...] which plane normal
-
-for j=1:size(plane_normal,1) % loop of plane_normals
-
-%________calculate h vector___________
-%from eq (3.7) in Poulsen
-%B is identity matrix in cubic material
-%gtnewGetReflections returns matrix of all reflections related by symmetry
-%to hkl, using Riso GetCubicSymOp.m
-
-hkl_reflections = gtCrystSignedHKLs(hkl(j,:), gtCrystGetSymmetryOperators([], spacegroup) );
-if size(hkl_reflections,1)~=1%skip if hkl = [0 0 0]
-
-for i=1:length(hkl_reflections) % loop of hkl reflections
-h = hkl_reflections(i,:);
-h = h./(sqrt(h*h'));
-
-%________calculate y vector___________
-% y is the normalised scattering vector in sample system
-% == plane normal
-y = plane_normal(j,:);
-
-
-%_____calculate rodrigues space vector_____
-%formulae (3.16), (3.17) in Poulsen
-
-r_zero = cross(h,y)./(1+(dot(h,y)));
-r_direction = (h + y)./(1+(dot(h,y)));
-
-%collect plot data
-lines.cor=[lines.cor; r_zero r_direction];
-lines.plnorm=[lines.plnorm; j];
-
-%max line length needed ~1.5, say +-1.5 for saftey
-t=(-1.5:0.01:1.5)';
-r_plot = repmat(r_zero,length(t),1) + (t*r_direction);
-
-%remove everything outside of a fundamental region (assume cube)
-lim=sqrt(2)-1;
-for i=1:3
-dummy=find(r_plot(:,i)<-lim | r_plot(:,i)>lim);
-r_plot(dummy,:)=[];
-end
-
-%remove the corners of the truncated cube
-lim=1/sqrt(3);
-%calculate distance of each point along each <111> directions
-vector=[1 1 1]./sqrt(3);%normalise
-vector=repmat(vector,size(r_plot,1),1);
-tmp = sum( ((abs(r_plot)).*vector),2); %projections along <111>
-dummy=find(tmp>lim);
-length(dummy);
-r_plot(dummy,:)=[];
-
-
-%collect all plot data for analysis
-save_plot=[save_plot; r_plot];
-
-
-%____plot______
-if plot_flag
- if all(hkl(j,:)==[1 1 1]) %(111) - red
- plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'r');
- elseif all(hkl(j,:)==[1 1 0]) %(110) - cyan - BCC
- plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'c');
- elseif all(hkl(j,:)==[0 0 2]) %(002) - blue
- plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'b');
- elseif all(hkl(j,:)==[2 1 1]) %(110) - magenta - BCC
- plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'m');
- elseif all(hkl(j,:)==[2 2 0]) %(220) - green
- plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'g');
- elseif all(hkl(j,:)==[3 1 1]) %(311) - black
- plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'k');
- end
-end
-
-end%loop through all hkl variants
-end%if hkl is okay (skip hkl = [0 0 0])
-
-%uncomment this to see line by line in rodrigues space.
-%k=waitforbuttonpress %can watch as each is plotted in turn
-
-end%loop through all normals
-
-if plot_flag
-axis equal
-end
-
-
-%find intersections
-lim=sqrt(2)-1;
-max_score=0;
-max_coords=[0 0 0];
-score=[];
-count=1;
-for i=-lim:inc:lim
-for j=-lim:inc:lim
-for k=-lim:inc:lim
-
-dummy = find(save_plot(:,1)>i & save_plot(:,1)<(i+inc) &save_plot(:,2)>j & save_plot(:,2)<(j+inc) &save_plot(:,3)>k & save_plot(:,3)<(k+inc));
-score(count)=length(dummy);
-if (score(count)>max_score)
- max_score=score;
- max_coords=[i j k];
-end
-
-count=count+1;
-end
-end
-end
-
-element = inc*[0 0 0; 0 0 1; 0 1 0; 0 1 1; 1 0 0; 1 0 1; 1 1 0; 1 1 1];
-max_element = repmat(max_coords,8,1)+element;
-max_coords = max_coords+0.5*[inc inc inc]; %this is the approximate solution
-
-
-
-
-if plot_flag
- plot3(max_element(:,1),max_element(:,2),max_element(:,3),'go')
- plot3(max_coords(:,1),max_coords(:,2),max_coords(:,3),'co')
-end
-
-%Peter's bit - more precise handling of intersections
-% Determine precise point of intersection inside the sampling element
-
-%normalise the r_direction part of the line equations
-%note that it is necessary to normalise vectors for the distance calculations to
-%work!
-lines.cor=normlinedir(lines.cor);
-
-line_cand.cor=[];%line equations which pass close to max_coords
-line_cand.plnorm=[];%to which plane normal (struct_id entry) does line belong - to ensure only one line from each normal
-
-% pick lines which pass close to the sampling element: candidate lines
-distances=pointtolinedist(max_coords,lines.cor);
-dummy=find(distances< inc_ext);
-line_cand.cor = lines.cor(dummy,:);
-line_cand.plnorm = lines.plnorm(dummy);
-
-
-% If more than one line from a given plane is present in the set, then choose the one that fits the best.
-pof_prel=pofintersectexp(line_cand.cor)'; %initial least squares best intersection
-
-% Alternative (cleaner!) loop to find the closest line from a given plane
-% to the point of intersection, if more than one is present
- distances=pointtolinedist(pof_prel, line_cand.cor);%distance from each candidate line to point
-dummy2=[];%list of lines to delete
-for k=1:size(plane_normal,1)
- dummy=find(plane_normal==k);
- if length(dummy)>2%if more than one line, delete all but the closest
- %- make a list of those to delete
- dummy2 = [dummy2 find(distances~=min(distances(dummy)) & line_cand.plnorm ==k)];
- end
-end
-line_cand.cor(dummy2,:)=[];%do the deletion
-line_cand.plnorm(dummy2,:)=[];%do the deletion
-
-%line_cand now contains information only about lines which pass near the
-%approx intersection, and only the best line arising from each plane
-%normal
-
-%repeat loop, to try and choose more precisely
-pof_prel=pofintersectexp(line_cand.cor)'; %least squares best intersection
-
-%find distances of remaining lines to point
-distances=pointtolinedist(pof_prel, line_cand.cor);%distance from each candidate line to point
-test_distance = consistancy_test*inc;
-dummy2=find(distances>test_distance);
-
-%%display the distances and which plane normals they come from
-%distances
-%line_cand.plnorm
-
-line_cand.cor(dummy2,:)=[];%do the deletion
-line_cand.plnorm(dummy2,:)=[];%do the deletion
-
-
-%%[line_cand.cor,line_cand.plnorm]
-
-%final least squares fit with good lines only, display results
-pof=pofintersectexp(line_cand.cor)';
-dsqr=pointtolinedistsqrsum(max_coords, line_cand.cor);
-
-%consistancy check part
-good_plnorms = 1:size(plane_normal,1);%index of plane normals
-good_plnorms = ismember(good_plnorms,line_cand.plnorm);
-
-
-if plot_flag
- plot3(pof(1),pof(2),pof(3),'mo','markersize',25);
- plot3(pof(1),pof(2),pof(3),'m*','markersize',25);
-end
-
diff --git a/4_grains/plot_rodrigues_consistancy_grain.m b/4_grains/plot_rodrigues_consistancy_grain.m
deleted file mode 100644
index 2ad69d9d5886719bd5e303ba35e9e84692d77e3f..0000000000000000000000000000000000000000
--- a/4_grains/plot_rodrigues_consistancy_grain.m
+++ /dev/null
@@ -1,263 +0,0 @@
-function [read_out, pof, good_plnorms, dsqr]=plot_rodrigues_consistancy_grain(grainids, spacegroup, plot_flag)
-
-%plane_normal = list of plane normals in sample coords
-%hkl a vector of reflection type eg [1 1 1; 0 0 2; 2 2 0 etc]
-%plot_flag: 1 to plot figure, 0 to not plot anything
-%add a consistancy check - which plane normals don't join at the
-%intersection?
-
-%grain - tool for looking at grains/twins post reconstruction. Given one or more
-%grainids, do the consistancy check.
-
-inc=0.03; % scan space in these steps (size of the sampling element) for coarse solution
-inc_ext=1.1*inc; % extended space to determine precise point of intersection; (should be larger than 0.87*inc)
-consistancy_test=0.2; %select those lines within distance consistancy_test*inc of the intersection (was 0.2)
-
-save_plot=[];
-
-if isempty(plot_flag)
-plot_flag=1;%always plot
-end
-
-
-%read in data
-plane_normal=[];
-hkl=[];
-read_out=[];
-for i=grainids
- tmp=load(sprintf('4_grains/grain%d_/grain%d_.mat', i, i), 'plane_normal', 'pl', 'hkl', 'index');
-
- tmp.index=tmp.index(1:length(tmp.hkl));
- if isfield(tmp, 'pl')
- plane_normal=[plane_normal; tmp.pl(find(tmp.index), :)];%try using only selected reflections
- else
- plane_normal=[plane_normal; tmp.plane_normal(find(tmp.index), :)];%try using only selected reflections
- end
- hkl=[hkl; tmp.hkl(find(tmp.index), :)];%try using only selected reflections
- read_out=[read_out length(find(tmp.index))];
-end
-
-
-[pof, good_plnorms, dsqr]=plot_rodrigues_consistancy(plane_normal,hkl,spacegroup,plot_flag);
-
-
-good_plnorms2=good_plnorms;
-
-for i=1:length(read_out)
- tmp=good_plnorms2(1:read_out(i));
- good_plnorms2(1:read_out(i))=[];
- read_out(i)=length(find(tmp))/length(tmp);
-end
-
-
-
-
-% lines.cor=[]; %[...; r_zero r_direction; ...] equations
-% lines.plnorm=[]; %[...; j; ...] which plane normal
-%
-% for j=1:size(plane_normal,1) % loop of plane_normals
-%
-% %________calculate h vector___________
-% %from eq (3.7) in Poulsen
-% %B is identity matrix in cubic material
-% %gtnewGetReflections returns matrix of all reflections related by symmetry
-% %to hkl, using Riso GetCubicSymOp.m
-%
-% hkl_reflections = gtCrystSignedHKLs(hkl(j,:), gtCrystGetSymmetryOperators([], spacegroup));
-% if size(hkl_reflections,1)~=1%skip if hkl = [0 0 0]
-%
-% for i=1:length(hkl_reflections) % loop of hkl reflections
-% h = hkl_reflections(i,:);
-% h = h./(sqrt(h*h'));
-%
-% %________calculate y vector___________
-% % y is the normalised scattering vector in sample system
-% % == plane normal
-% y = plane_normal(j,:);
-%
-%
-% %_____calculate rodrigues space vector_____
-% %formulae (3.16), (3.17) in Poulsen
-%
-% r_zero = cross(h,y)./(1+(dot(h,y)));
-% r_direction = (h + y)./(1+(dot(h,y)));
-%
-% %collect plot data
-% lines.cor=[lines.cor; r_zero r_direction];
-% lines.plnorm=[lines.plnorm; j];
-%
-% %max line length needed ~1.5, say +-1.5 for saftey
-% t=(-1.5:0.01:1.5)';
-% r_plot = repmat(r_zero,length(t),1) + (t*r_direction);
-%
-% %remove everything outside of a fundamental region (assume cube)
-% lim=sqrt(2)-1;
-% for i=1:3
-% dummy=find(r_plot(:,i)<-lim | r_plot(:,i)>lim);
-% r_plot(dummy,:)=[];
-% end
-%
-% %remove the corners of the truncated cube
-% lim=1/sqrt(3);
-% %calculate distance of each point along each <111> directions
-% vector=[1 1 1]./sqrt(3);%normalise
-% vector=repmat(vector,size(r_plot,1),1);
-% tmp = sum( ((abs(r_plot)).*vector),2); %projections along <111>
-% dummy=find(tmp>lim);
-% length(dummy);
-% r_plot(dummy,:)=[];
-%
-%
-% %collect all plot data for analysis
-% save_plot=[save_plot; r_plot];
-%
-%
-% %____plot______
-% if plot_flag
-% if all(hkl(j,:)==[1 1 1]) %(111) - red
-% plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'r');
-% elseif all(hkl(j,:)==[1 1 0]) %(110) - cyan - BCC
-% plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'c');
-% elseif all(hkl(j,:)==[0 0 2]) %(002) - blue
-% plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'b');
-% elseif all(hkl(j,:)==[2 1 1]) %(110) - magenta - BCC
-% plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'m');
-% elseif all(hkl(j,:)==[2 2 0]) %(220) - green
-% plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'g');
-% elseif all(hkl(j,:)==[3 1 1]) %(311) - black
-% plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'k');
-% end
-% end
-%
-% end%loop through all hkl variants
-% end%if hkl is okay (skip hkl = [0 0 0])
-%
-% %uncomment this to see line by line in rodrigues space.
-% %k=waitforbuttonpress %can watch as each is plotted in turn
-%
-% end%loop through all normals
-%
-% if plot_flag
-% axis equal
-% end
-%
-%
-% %find intersections
-% lim=sqrt(2)-1;
-% max_score=0;
-% max_coords=[0 0 0];
-% score=[];
-% count=1;
-% for i=-lim:inc:lim
-% for j=-lim:inc:lim
-% for k=-lim:inc:lim
-%
-% dummy = find(save_plot(:,1)>i & save_plot(:,1)<(i+inc) &save_plot(:,2)>j & save_plot(:,2)<(j+inc) &save_plot(:,3)>k & save_plot(:,3)<(k+inc));
-% score(count)=length(dummy);
-% if (score(count)>max_score)
-% max_score=score;
-% max_coords=[i j k];
-% end
-%
-% count=count+1;
-% end
-% end
-% end
-%
-% max_score;
-% max_coords;
-% element = inc*[0 0 0; 0 0 1; 0 1 0; 0 1 1; 1 0 0; 1 0 1; 1 1 0; 1 1 1];
-% max_element = repmat(max_coords,8,1)+element;
-% max_coords = max_coords+0.5*[inc inc inc]; %this is the approximate solution
-%
-%
-%
-%
-% if plot_flag
-% plot3(max_element(:,1),max_element(:,2),max_element(:,3),'go')
-% plot3(max_coords(:,1),max_coords(:,2),max_coords(:,3),'co')
-% end
-%
-% %Peter's bit - more precise handling of intersections
-% % Determine precise point of intersection inside the sampling element
-%
-% %normalise the r_direction part of the line equations
-% %note that it is necessary to normalise vectors for the distance calculations to
-% %work!
-% lines.cor=normlinedir(lines.cor);
-%
-% line_cand.cor=[];%line equations which pass close to max_coords
-% line_cand.plnorm=[];%to which plane normal (struct_id entry) does line belong - to ensure only one line from each normal
-%
-% % pick lines which pass close to the sampling element: candidate lines
-% distances=pointtolinedist(max_coords,lines.cor);
-% dummy=find(distances< inc_ext);
-% line_cand.cor = lines.cor(dummy,:);
-% line_cand.plnorm = lines.plnorm(dummy);
-%
-%
-% % If more than one line from a given plane is present in the set, then choose the one that fits the best.
-% pof_prel=pofintersectexp(line_cand.cor)'; %initial least squares best intersection
-%
-% % Alternative (cleaner!) loop to find the closest line from a given plane
-% % to the point of intersection, if more than one is present
-% distances=pointtolinedist(pof_prel, line_cand.cor);%distance from each candidate line to point
-% dummy2=[];%list of lines to delete
-% for k=1:size(plane_normal,1)
-% dummy=find(plane_normal==k);
-% if length(dummy)>2%if more than one line, delete all but the closest
-% %- make a list of those to delete
-% dummy2 = [dummy2 find(distances~=min(distances(dummy)) & line_cand.plnorm ==k)];
-% end
-% end
-% line_cand.cor(dummy2,:)=[];%do the deletion
-% line_cand.plnorm(dummy2,:)=[];%do the deletion
-%
-% %line_cand now contains information only about lines which pass near the
-% %approx intersection, and only the best line arising from each plane
-% %normal
-%
-% %repeat loop, to try and choose more precisely
-% pof_prel=pofintersectexp(line_cand.cor)'; %least squares best intersection
-%
-% %find distances of remaining lines to point
-% distances=pointtolinedist(pof_prel, line_cand.cor);%distance from each candidate line to point
-% test_distance = consistancy_test*inc;
-% dummy2=find(distances>test_distance);
-%
-% %%display the distances and which plane normals they come from
-% %distances
-% %line_cand.plnorm
-%
-% line_cand.cor(dummy2,:)=[];%do the deletion
-% line_cand.plnorm(dummy2,:)=[];%do the deletion
-%
-%
-% %%[line_cand.cor,line_cand.plnorm]
-%
-% %final least squares fit with good lines only, display results
-% max_coords;
-% pof_prel;
-% pof=pofintersectexp(line_cand.cor)'
-% dsqr=pointtolinedistsqrsum(max_coords, line_cand.cor)
-%
-% %consistancy check part
-% good_plnorms = 1:size(plane_normal,1);%index of plane normals
-% good_plnorms = ismember(good_plnorms,line_cand.plnorm);
-%
-% good_plnorms2=good_plnorms;
-%
-% for i=1:length(read_out)
-% tmp=good_plnorms2(1:read_out(i));
-% good_plnorms2(1:read_out(i))=[];
-% read_out(i)=length(find(tmp))/length(tmp);
-% end
-%
-%
-%
-%
-% if plot_flag
-% plot3(pof(1),pof(2),pof(3),'mo','markersize',25);
-% plot3(pof(1),pof(2),pof(3),'m*','markersize',25);
-% end
-%
diff --git a/4_grains/plot_rodrigues_consistancy_sortgrains.m b/4_grains/plot_rodrigues_consistancy_sortgrains.m
deleted file mode 100644
index 09e4ad268eca2dff0857ee058dbf4767d4e1fc19..0000000000000000000000000000000000000000
--- a/4_grains/plot_rodrigues_consistancy_sortgrains.m
+++ /dev/null
@@ -1,403 +0,0 @@
-function [pof, good_plnorms, dsqr]=plot_rodrigues_consistancy_sortgrains(plane_normal,hkl,spacegroup,latticepar,plot_flag)
-% [pof, good_plnorms, dsqr]=plot_rodrigues_consistancy_sortgrains(plane_normal,hkl,spacegroup,latticepar,plot_flag)
-% -----------------------------------------------------------------------------------------------------------------
-% Peter's modification of plot_rodrigues_consistancy for sorting grains;
-%
-% INPUT:
-% plane_normal = list of plane normals in sample coords
-% hkl = a vector of reflection type eg [1 1 1; 0 0 2; 2 2 0 etc...]
-% for hexagonal this this be a 4 index notation
-% spacegroup = parameters.cryst(phaseid).spacegroup
-% latticepar = lattice parameters
-% plot_flag = true to plot figure, false to not plot anything
-%
-% Version 001 19-10-2007 by AKing
-% add consistancy check
-% Version 002 19-10-2012 by LNervo
-% add spacegroup argument
-
-save_plot=[];
-
-if plot_flag
- figure;
- hold on;
-end
-
-lines.cor=[];
-lines.plnorm=[];
-consistancy_test=0.2; %select those lines within distance consistancy_test*inc of the intersection (was 0.2)
-
-switch spacegroup
- case {225, 229} % ie cubic lattice
- inc=0.1; % scan space in these steps (size of the sampling element)
- inc_ext=1.1*inc; % extended space to determine precise point of intersection; (should be larger than 0.87*inc)
- lim=sqrt(2)-1;
- limz=sqrt(2)-1;
- for j=1:size(plane_normal,1) % loop of plane_normals
-
- %________calculate h vector___________
- %from eq (3.7) in Poulsen
- %B is identity matrix in cubic material
- %gtnewGetReflections returns matrix of all reflections related by symmetry
- %to hkl, using Riso GetCubicSymOp.m
-
- hkl_reflections = gtCrystSignedHKLs(hkl(j,:), gtCrystGetSymmetryOperators([], spacegroup) );
- if size(hkl_reflections,1)~=1%skip if hkl = [0 0 0]
- for i=1:size(hkl_reflections,1) % loop of hkl reflections
- h = hkl_reflections(i,:);
- h = h./(sqrt(h*h'));
- %________calculate y vector___________
- % y is the normalised scattering vector in sample system
- % == plane normal
- y = plane_normal(j,:);
-
- %_____calculate rodrigues space vector_____
- %formulae (3.16), (3.17) in Poulsen
-
- r_zero = cross(h,y)./(1+(dot(h,y)));
- r_direction = (h + y)./(1+(dot(h,y)));
-
- lines.cor=[lines.cor; r_zero r_direction];
- lines.plnorm=[lines.plnorm; j];
-
- %max line length needed ~1.5, say +-1.5 for saftey
- t=(-1.5:0.01:1.5)';
- r_plot = repmat(r_zero,length(t),1) + (t*r_direction);
-
- %remove everything outside of a fundamental region (assume cube)
- lim=sqrt(2)-1;
- for i=1:3
- dummy=find(r_plot(:,i)<-lim | r_plot(:,i)>lim);
- r_plot(dummy,:)=[];
- end
-
- %remove the corners of this region also
- lim=1/sqrt(3);
- %calculate distance of each point along each <111> directions
- vector=[1 1 1]./sqrt(3);%normalise
- vector=repmat(vector,size(r_plot,1),1);
- tmp = sum( ((abs(r_plot)).*vector),2); %projections along <111>
- dummy=find(tmp>lim);
- r_plot(dummy,:)=[];
-
-
-
- %collect all plot data for analysis
- save_plot=[save_plot; r_plot];
-
-
-
- %____plot______
- if plot_flag
- if all(hkl(j,:)==[1 1 1]) %(111) - red
- plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'r');
-
- elseif all(hkl(j,:)==[0 0 2]) %(002) - blue
- plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'b');
-
- elseif all(hkl(j,:)==[2 2 0]) %(220) - green
- plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'g');
-
- elseif all(hkl(j,:)==[3 1 1]) %(311) - black
- plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'k');
-
- else %all others in pink
- plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'m');
-
- end
- end
-
- end%loop through all hkl variants
- end%if hkl is okay (skip hkl = [0 0 0])
- %uncomment this to see line by line in rodrigues space.
- %k=waitforbuttonpress %can watch as each is plotted in turn
-
- end%loop through all normals
- %%%%%%%end case 225/ 229
-
-
- case 663
-
- inc=0.08;% scan space in these steps (size of the sampling element)
- inc_ext=1.1*inc; % extended space to determine precise point of intersection; (should be larger than 0.87*inc)
- lim=1;
- limz=2-sqrt(3);
-
- plane_normal_cart=plane_normal;
- for j=1:size(plane_normal_cart,1) % loop of plane_normals
- %________calculate h vector___________
- %from eq (3.7) in Poulsen
- %B is identity matrix in cubic material
- %gtnewGetReflections returns matrix of all reflections related by symmetry
- %to hkl, using Riso GetCubicSymOp.m
-
- hkil_reflections_hex = gtCrystSignedHKLs(hkl(j,:), gtCrystGetSymmetryOperators([], spacegroup) );
-
- if size(hkil_reflections_hex,1)~=1%skip if hkl = [0 0 0]
- for i=1:size(hkil_reflections_hex,1) % loop of hkil reflections
-
- % passage de l espace direct hexagonal en coordonnees cartesiennes
- h(1)= hkil_reflections_hex(i,1) + 0.5 * hkil_reflections_hex(i,2);
- h(2)= 3^0.5/2 *hkil_reflections_hex(i,2);
- h(3)= hkil_reflections_hex(i,4);
- % normalisation of the cartesian coordinates space
- % allow for c/a ratio in hexagonal unit cell
- h(1)=h(1)*2/(sqrt(3)*latticepar(1));
- h(2)=h(2)*2/(sqrt(3)*latticepar(1));
- h(3)=h(3)/latticepar(3);
- % normalise to unit vector
- h = h./(sqrt(h*h'));
-
- %________calculate y vector___________
- % y is the normalised scattering vector in cartesian system
- % == plane normal
- y = plane_normal_cart(j,:);
- y = y./(sqrt(y*y'));
-
- %_____calculate rodrigues space vector_____
- %formulae (3.16), (3.17) in Poulsen
-
- r_zero = cross(h,y)./(1+(dot(h,y)));
- r_direction = (h + y)./(1+(dot(h,y)));
-
- lines.cor=[lines.cor; r_zero r_direction];
- lines.plnorm=[lines.plnorm; j];
-
-
- %max line length needed ~1.5, say +-1.5 for saftey
- t=(-3:0.01:3)';
- r_plot = repmat(r_zero,length(t),1) + (t*r_direction);
-
- %remove everything outside of a fundamental region (assume cf Heinz_Neumann acta. crist.)
- % base and top of the prism
- dummy=find(abs(r_plot(:,1))>=1 | abs(r_plot(:,2))>=1 | abs(r_plot(:,3))>=(2-sqrt(3)) );
- r_plot(dummy,:)=[];
- % remove facette of the faces
- dummy=find((r_plot(:,1)+sqrt(3)*r_plot(:,2)>=2) | (r_plot(:,1)-sqrt(3)*r_plot(:,2))>=2);
- r_plot(dummy,:)=[];
- dummy=find((-r_plot(:,1)+sqrt(3)*r_plot(:,2)>=2) | (-r_plot(:,1)-sqrt(3)*r_plot(:,2))>=2);
- r_plot(dummy,:)=[];
- dummy=find((r_plot(:,2)+sqrt(3)*r_plot(:,1))>=2 | (r_plot(:,2)-sqrt(3)*r_plot(:,1))>=2);
- r_plot(dummy,:)=[];
- dummy=find((-r_plot(:,2)+sqrt(3)*r_plot(:,1))>=2 | (-r_plot(:,2)-sqrt(3)*r_plot(:,1))>=2);
- r_plot(dummy,:)=[];
- %collect all plot data for analysis
- save_plot=[save_plot; r_plot];
-
-
-
- %____plot______ % adapted to the snow case
- if plot_flag
- % keyboard
- %plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'m');
- hkil_hex=hkl;
- if all(hkil_hex(j,:)==[0 0 0 2])
-
- plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'r');
- %figure(2)
- %quiver3(0,0,0,h(1),h(2),h(3),'r')
-
- elseif all(hkil_hex(j,:)==[1 1 -2 0])
-
- plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'b');
- %figure(2)
- %quiver3(0,0,0,h(1),h(2),h(3),'b')
-
- elseif all(hkil_hex(j,:)==[1 -1 0 0])
-
- plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'g');
- %figure(2)
- %quiver3(0,0,0,h(1),h(2),h(3),'g')
- elseif all(hkil_hex(j,:)==[1 -1 0 1])
-
- plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'k');
- %figure(2)
- %quiver3(0,0,0,h(1),h(2),h(3),'k')
- elseif all(hkil_hex(j,:)==[1 1 -2 2])
-
- plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'c');
- %figure(2)
- %quiver3(0,0,0,h(1),h(2),h(3),'c')
-
- % elseif all(hkil_hex(j,:)==[1 1 -2 1])
- % plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'k');
- %
- % elseif all(hkil_hex(j,:)==[1 -1 0 2])
- % plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'m');
- %
- % elseif all(hkil_hex(j,:)==[1 1 -2 3])
- % plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'c');
- %
- % elseif all(hkil_hex(j,:)==[1 -1 0 3])
- % plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'y');
- %
- % elseif all(hkil_hex(j,:)==[1 1 -2 4])
- % plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'r:');
- %
- % elseif all(hkil_hex(j,:)==[1 -1 0 4])
- % plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'b:');
- %
- % elseif all(hkil_hex(j,:)==[0 0 0 4])
- % plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'g:');
- %
- % elseif all(hkil_hex(j,:)==[-2 0 2 1])
- % plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'m:');
- %
- % elseif all(hkil_hex(j,:)==[-2 0 2 3])
- % plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'c:');
- %
- % elseif all(hkil_hex(j,:)==[-3 1 2 1])
- % plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'y:');
- %
- % elseif all(hkil_hex(j,:)==[-3 1 2 0])
- % plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'b-');
- % else
- % plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'m');
- %
- end
- % legend(' 0 0 0 -2', '1 1 -2 0', '1 -1 0 0', '1 1 -2 1', '1 -1 0 2','1 1 -2 3', '1 1 -2 4','1 -1 0 4', ' 0 0 0 4','-2 0 2 1','-2 0 2 3','-3 1 2 3','-3 1 2 0');
- end
-
- end %loop through all hkl variants
- end %if hkl is okay (skip hkl = [0 0 0])
- % uncomment this to see line by line in rodrigues space.
- % k=waitforbuttonpress %can watch as each is plotted in turn
-
- end%loop through all normals
-
-end % end switch space group
-
-
-
-if plot_flag
- axis equal
-end
-lines.cor;
-lines.plnorm;
-size(lines.cor)
-
-
-
-%find intersections
-tic
-max_score=0;
-max_coords=[0 0 0];
-score=[];
-count=1;
-for i=-lim:inc:lim
- for j=-lim:inc:lim
- for k=-limz:inc:limz
-
- dummy = find(save_plot(:,1)>i & save_plot(:,1)<(i+inc) &save_plot(:,2)>j & save_plot(:,2)<(j+inc) &save_plot(:,3)>k & save_plot(:,3)<(k+inc));
- score(count)=length(dummy);
- if (score(count)>max_score)
- max_score=score;
- max_coords=[i j k];
- end
-
- count=count+1;
- end
- end
-end
-
-element = inc*[0 0 0; 0 0 1; 0 1 0; 0 1 1; 1 0 0; 1 0 1; 1 1 0; 1 1 1];
-max_element = repmat(max_coords,8,1)+element;
-max_coords = max_coords+0.5*[inc inc inc]; %this is the approximate solution
-
-if plot_flag
- plot3(max_element(:,1),max_element(:,2),max_element(:,3),'go')
- plot3(max_coords(:,1),max_coords(:,2),max_coords(:,3),'co')
-end
-
-
-% Determine precise point of intersection inside the sampling element
-
-lines.cor=normlinedir(lines.cor);
-
-line_cand.cor=[];%line equations which pass close to max_coords
-line_cand.plnorm=[];%to which plane normal (struct_id entry) does line belong - to ensure only one line from each normal
-
-%for j=1:size(lines.cor,1) % picks lines which pass close to the sampling element: candidate lines
-% if (pointtolinedist(max_coords, lines.cor(j,:)) < inc_ext)
-% line_cand.cor=[line_cand.cor; lines.cor(j,:)];
-% line_cand.plnorm=[line_cand.plnorm; lines.plnorm(j)];
-% end
-%end
-distances=pointtolinedist(max_coords,lines.cor);
-dummy=find(distances< inc_ext);
-line_cand.cor = lines.cor(dummy,:);
-line_cand.plnorm = lines.plnorm(dummy);
-
-
-% If more than one line from a given plane is present in the set, then
-% choose the one that fits the best.
-pof_prel=pofintersectexp(line_cand.cor)';
-
-% for k=1:size(plane_normal,1) % delete redundant lines from line_cand
-% dis=1e+38;
-% m=1;
-% while (m <= size(line_cand.cor,1))
-% if (line_cand.plnorm(m) == k)
-% if pointtolinedist(pof_prel, line_cand.cor(m,:)) < dis
-% dis=pointtolinedist(pof_prel, line_cand.cor(m,:));
-% dis_line.cor=line_cand.cor(m,:);
-% end
-% line_cand.cor(m,:)=[];
-% line_cand.plnorm(m)=[];
-% m=m-1;
-% end
-% m=m+1;
-% end
-% if dis ~= 1e+38
-% line_cand.cor=[line_cand.cor; dis_line.cor(1:6)];
-% line_cand.plnorm=[line_cand.plnorm; k];
-% end
-% end
-
-% Alternative (cleaner!) loop to find the closest line from a given plane
-% to the point of intersection, if more than one is present
-distances=pointtolinedist(pof_prel, line_cand.cor);%distance from each candidate line to point
-dummy2=[];%list of lines to delete
-for k=1:size(plane_normal,1)
- dummy=find(plane_normal==k);
- if length(dummy)>2%if more than one line, delete all but the closest
- %- make a list of those to delete
- dummy2 = [dummy2 find(distances~=min(distances(dummy)) & line_cand.plnorm ==k)];
- end
-end
-line_cand.cor(dummy2,:)=[];%do the deletion
-line_cand.plnorm(dummy2,:)=[];%do the deletion
-
-
-%line_cand now contains information only about lines which pass near the
-%approx intersection, and only the best line arising from each plane
-%normal
-
-%repeat loop, to try and choose more precisely
-pof_prel=pofintersectexp(line_cand.cor)'; %least squares best intersection
-
-%find distances of remaining lines to point
-distances=pointtolinedist(pof_prel, line_cand.cor);%distance from each candidate line to point
-test_distance = consistancy_test*inc;
-dummy2=find(distances>test_distance);
-
-%display the distances and which plane normals they come from
-
-line_cand.cor(dummy2,:)=[];%do the deletion
-line_cand.plnorm(dummy2,:)=[];%do the deletion
-
-
-%final least squares fit with good lines only, display results
-pof=pofintersectexp(line_cand.cor)';
-dsqr=pointtolinedistsqrsum(max_coords, line_cand.cor);
-
-%consistancy check part
-good_plnorms = 1:size(plane_normal,1);%index of plane normals
-good_plnorms = ismember(good_plnorms,line_cand.plnorm);
-
-
-if plot_flag
- plot3(pof(1),pof(2),pof(3),'mo','markersize',25);
- plot3(pof(1),pof(2),pof(3),'m*','markersize',25);
-end
-toc
diff --git a/4_grains/plot_rodrigues_indexed.m b/4_grains/plot_rodrigues_indexed.m
deleted file mode 100644
index d3497ba294ea93081877969a1272defbe7799f86..0000000000000000000000000000000000000000
--- a/4_grains/plot_rodrigues_indexed.m
+++ /dev/null
@@ -1,111 +0,0 @@
-function max_coords=plot_rodrigues_indexed(plane_normal,hkl)
-
-%plane_normal = list of plane normals in sample coords
-% for caase where reflections have been indexed - ie (0 2 -2)
-% rather than {2 2 0}
-%hkl a vector of reflections [1 -1 1; 0 0 2; -2 2 0 etc]
-
-save_plot=[];
-figure;
-hold on;
-
-for j=1:size(plane_normal,1)
-
-%________calculate h vector___________
-%from eq (3.7) in Poulsen
-%B is identity matrix in cubic material
-%gtnewGetReflections returns matrix of all reflections related by symmetry
-%to hkl, using Riso GetCubicSymOp.m
-
-
-h = hkl(j,:);
-h = h./(sqrt(h*h'));
-
-%________calculate y vector___________
-% y is the normalised scattering vector in sample system
-% == plane normal
-y = plane_normal(j,:);
-
-
-%_____calculate rodrigues space vector_____
-%formulae (3.16), (3.17) in Poulsen
-
-r_zero = cross(h,y)./(1+(dot(h,y)));
-r_direction = (h + y)./(1+(dot(h,y)));
-
-%max line length needed ~1.5, say +-1.5 for saftey
-t=(-1.5:0.01:1.5)';
-r_plot = repmat(r_zero,length(t),1) + (t*r_direction);
-
-%remove everything outside of a fundamental region (assume cube)
-lim=sqrt(2)-1;
-for i=1:3
-dummy=find(r_plot(:,i)<-lim | r_plot(:,i)>lim);
-r_plot(dummy,:)=[];
-end
-
-%collect all plot data for analysis
-save_plot=[save_plot; r_plot];
-
-
-
-%____plot______
- if hkl(j,:)*hkl(j,:)'==[1 1 1]*[1 1 1]' %(111) - red
- plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'r');
-
- elseif hkl(j,:)*hkl(j,:)'==[0 0 2]*[0 0 2]' %(002) - blue
- plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'b');
-
- elseif hkl(j,:)*hkl(j,:)'==[2 2 0]*[2 2 0]' %(220) - green
- plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'g');
-
- elseif hkl(j,:)*hkl(j,:)'==[3 1 1]*[3 1 1]' %(311) - black
- plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'k');
-
- end
-
-%uncomment this to see line by line in rodrigues space.
-%k=waitforbuttonpress %can watch as each is plotted in turn
-
-end%loop through all normals
-
-axis equal
-
-
-%find intersections
-
-inc=0.04;%scan space in these steps (size of the sampling element)
-max_score=0;
-max_coords=[0 0 0];
-score=[];
-count=1;
-for i=-lim:inc:lim
-for j=-lim:inc:lim
-for k=-lim:inc:lim
-
-dummy = find(save_plot(:,1)>i & save_plot(:,1)<(i+inc) &save_plot(:,2)>j & save_plot(:,2)<(j+inc) &save_plot(:,3)>k & save_plot(:,3)<(k+inc));
-score(count)=length(dummy);
-if (score(count)>max_score)
- max_score=score;
- max_coords=[i j k];
-end
-
-count=count+1;
-end
-end
-end
-
-max_score;
-max_coords;
-element = inc*[0 0 0; 0 0 1; 0 1 0; 0 1 1; 1 0 0; 1 0 1; 1 1 0; 1 1 1];
-max_element = repmat(max_coords,8,1)+element;
-max_coords = max_coords+0.5*[inc inc inc];%this is the approximate solution
-
-plot3(max_element(:,1),max_element(:,2),max_element(:,3),'go')
-plot3(max_coords(:,1),max_coords(:,2),max_coords(:,3),'co')
-
-
-
-
-
-
diff --git a/4_grains/plot_rodrigues_precise.m b/4_grains/plot_rodrigues_precise.m
deleted file mode 100644
index 8f2dc6ecfb81b1180e1ceae92e6eefcb59113e64..0000000000000000000000000000000000000000
--- a/4_grains/plot_rodrigues_precise.m
+++ /dev/null
@@ -1,205 +0,0 @@
-function [pof, dsqr]=plot_rodrigues_precise(plane_normal,hkl,spacegroup,plot_flag)
-% <<<< NOT USED >>>>
-%plane_normal = list of plane normals in sample coords
-%hkl a vector of reflection type eg [1 1 1; 0 0 2; 2 2 0 etc]
-%plot_flag: 1 to plot figure, 0 to not plot anything
-
-inc=0.1; % scan space in these steps (size of the sampling element)
-inc_ext=1.1*inc; % extended space to determine precise point of intersection; (should be larger than 0.87*inc)
-
-save_plot=[];
-
-if plot_flag
-figure;
-hold on;
-end
-
-lines.cor=[];
-lines.plnorm=[];
-
-for j=1:size(plane_normal,1) % loop of plane_normals
-
-%________calculate h vector___________
-%from eq (3.7) in Poulsen
-%B is identity matrix in cubic material
-%gtnewGetReflections returns matrix of all reflections related by symmetry
-%to hkl, using Riso GetCubicSymOp.m
-
-hkl_reflections = gtCrystSignedHKLs(hkl(j,:), gtCrystGetSymmetryOperators([],spacegroup));
-
-if size(hkl_reflections,1)~=1%skip if hkl = [0 0 0]
-for i=1:size(hkl_reflections,1) % loop of hkl reflections
-h = hkl_reflections(i,:);
-h = h./(sqrt(h*h'));
-
-%________calculate y vector___________
-% y is the normalised scattering vector in sample system
-% == plane normal
-y = plane_normal(j,:);
-
-
-%_____calculate rodrigues space vector_____
-%formulae (3.16), (3.17) in Poulsen
-
-r_zero = cross(h,y)./(1+(dot(h,y)));
-r_direction = (h + y)./(1+(dot(h,y)));
-
-
-
-
-
-
-lines.cor=[lines.cor; r_zero r_direction];
-lines.plnorm=[lines.plnorm; j];
-
-
-
-
-
-%max line length needed ~1.5, say +-1.5 for saftey
-t=(-1.5:0.01:1.5)';
-r_plot = repmat(r_zero,length(t),1) + (t*r_direction);
-
-%remove everything outside of a fundamental region (assume cube)
-lim=sqrt(2)-1;
-for i=1:3
-dummy=find(r_plot(:,i)<-lim | r_plot(:,i)>lim);
-r_plot(dummy,:)=[];
-end
-
-%remove the corners of this region also
-lim=1/sqrt(3);
-%calculate distance of each point along each <111> directions
-vector=[1 1 1]./sqrt(3);%normalise
-vector=repmat(vector,size(r_plot,1),1);
-tmp = sum( ((abs(r_plot)).*vector),2); %projections along <111>
-dummy=find(tmp>lim);
-r_plot(dummy,:)=[];
-
-
-
-%collect all plot data for analysis
-save_plot=[save_plot; r_plot];
-
-
-
-%____plot______
-if plot_flag
- if all(hkl(j,:)==[1 1 1]) %(111) - red
- plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'r');
-
- elseif all(hkl(j,:)==[0 0 2]) %(002) - blue
- plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'b');
-
- elseif all(hkl(j,:)==[2 2 0]) %(220) - green
- plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'g');
-
- elseif all(hkl(j,:)==[3 1 1]) %(311) - black
- plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'k');
-
- else %all others in pink
- plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'m');
-
- end
-end
-
-end%loop through all hkl variants
-end%if hkl is okay (skip hkl = [0 0 0])
-%uncomment this to see line by line in rodrigues space.
-%k=waitforbuttonpress %can watch as each is plotted in turn
-
-end%loop through all normals
-
-if plot_flag
-axis equal
-end
-
-%find intersections
-
-max_score=0;
-max_coords=[0 0 0];
-score=[];
-count=1;
-for i=-lim:inc:lim
-for j=-lim:inc:lim
-for k=-lim:inc:lim
-
-dummy = find(save_plot(:,1)>i & save_plot(:,1)<(i+inc) &save_plot(:,2)>j & save_plot(:,2)<(j+inc) &save_plot(:,3)>k & save_plot(:,3)<(k+inc));
-score(count)=length(dummy);
-if (score(count)>max_score)
- max_score=score;
- max_coords=[i j k];
-end
-
-count=count+1;
-end
-end
-end
-
-element = inc*[0 0 0; 0 0 1; 0 1 0; 0 1 1; 1 0 0; 1 0 1; 1 1 0; 1 1 1];
-max_element = repmat(max_coords,8,1)+element;
-max_coords = max_coords+0.5*[inc inc inc]; %this is the approximate solution
-
-if plot_flag
- plot3(max_element(:,1),max_element(:,2),max_element(:,3),'go')
- plot3(max_coords(:,1),max_coords(:,2),max_coords(:,3),'co')
-end
-
-
-
-
-
-
-
-
-
-
-% Determine precise point of intersection inside the sampling element
-
-lines.cor=normlinedir(lines.cor);
-
-line_cand.cor=[];%line equations which pass close to max_coords
-line_cand.plnorm=[];%to which plane normal (struct_id entry) does line belong - to ensure only one line from each normal
-
-for j=1:size(lines.cor,1) % picks lines which pass close to the sampling element: candidate lines
- if (pointtolinedist(max_coords, lines.cor(j,:)) < inc_ext)
- line_cand.cor=[line_cand.cor; lines.cor(j,:)];
- line_cand.plnorm=[line_cand.plnorm; lines.plnorm(j)];
- end
-end
-
-%normlinedir(line_cand.cor);
-
- % If more than one line from a given plane is present in the set, then choose the one that fits the best.
-
-pof_prel=pofintersectexp(line_cand.cor)';
-
-for k=1:size(plane_normal,1) % delete redundant lines from line_cand
- dis=1e+38;
- m=1;
- while (m <= size(line_cand.cor,1))
- if (line_cand.plnorm(m) == k)
- if pointtolinedist(pof_prel, line_cand.cor(m,:)) < dis
- dis=pointtolinedist(pof_prel, line_cand.cor(m,:));
- dis_line.cor=line_cand.cor(m,:);
- end
- line_cand.cor(m,:)=[];
- line_cand.plnorm(m)=[];
- m=m-1;
- end
- m=m+1;
- end
- if dis ~= 1e+38
- line_cand.cor=[line_cand.cor; dis_line.cor(1:6)];
- line_cand.plnorm=[line_cand.plnorm; k];
- end
-end
-
-pof=pofintersectexp(line_cand.cor)';
-dsqr=pointtolinedistsqrsum(max_coords, line_cand.cor);
-
-if plot_flag
- plot3(pof(1),pof(2),pof(3),'mo','markersize',25);
- plot3(pof(1),pof(2),pof(3),'m*','markersize',25);
-end
-
diff --git a/zUtil_Strain/gtKeepLine.m b/zUtil_Strain/gtKeepLine.m
deleted file mode 100644
index dff6ca51cb913ed04d777aba6cfde67e21040778..0000000000000000000000000000000000000000
--- a/zUtil_Strain/gtKeepLine.m
+++ /dev/null
@@ -1,26 +0,0 @@
-% to be deleted
-% new macro: gtIndexSelectLines
-
-function linekept=gtKeepLine(keeplist,line)
-
- linekept.id=line.id(keeplist);
- linekept.pairid=line.pairid(keeplist);
- linekept.aid=line.aid(keeplist);
- linekept.bid=line.bid(keeplist);
- linekept.theta=line.theta(keeplist);
- linekept.eta=line.eta(keeplist);
- linekept.omega=line.omega(keeplist);
- linekept.int=line.int(keeplist);
- linekept.bbxs=line.bbxs(keeplist);
- linekept.bbys=line.bbys(keeplist);
-
- linekept.ca=line.ca(keeplist,:);
- linekept.cb=line.cb(keeplist,:);
- linekept.dir=line.dir(keeplist,:);
- linekept.pl=line.pl(keeplist,:);
-
- linekept.thetatype=line.thetatype(keeplist);
- linekept.hkl=line.hkl(keeplist,:);
- %linekept.lRv=line.lRv(:,:,keeplist);
-
-end