diff --git a/1_preprocessing/gtPreprocessing.m b/1_preprocessing/gtPreprocessing.m
index 2b6dc5893fd4153f5b73a3fd67b02bd81a8e37fb..75b512e477e40de115ef0b9a4cd58ca8ea8b3923 100755
--- a/1_preprocessing/gtPreprocessing.m
+++ b/1_preprocessing/gtPreprocessing.m
@@ -299,6 +299,8 @@ parameters.acq.rotx=parameters.acq.rotu;
if isfield(parameters.acq, 'no_direct_beam') && strcmpi(parameters.acq.no_direct_beam,'true')
parameters.prep.udrift=zeros(1, totproj);
parameters.prep.vdrift=zeros(1, totproj);
+ % hardcoded sep 27 lnervo
+ parameters.prep.udriftabs=zeros(1,totproj);
else
detcenter_abs=parameters.acq.xdet/2+0.5;
bboxcenter_abs=(2*parameters.acq.bb(1)+parameters.acq.bb(3)-1)/2;
@@ -357,7 +359,6 @@ disp('Using OAR parameters:')
OAR_parameters.walltime=sprintf('%02d:%02d:00', OARh, OARm);
OAR_parameters.mem=2000;
-disp(sprintf('Number of OAR jobs: %d', OARjobs));
disp('OAR parameters:')
disp(OAR_parameters)
% could adjust this calculation
diff --git a/1_preprocessing/gtSetupCalibrationQuali.m b/1_preprocessing/gtSetupCalibrationQuali.m
index e1adfd3cb05d427791af810d240eb1c12f9660bc..16d95c3823077153d58703803f3e6b0a5b196704 100755
--- a/1_preprocessing/gtSetupCalibrationQuali.m
+++ b/1_preprocessing/gtSetupCalibrationQuali.m
@@ -52,7 +52,7 @@ if strcmpi(parameters.calib.type,'180degree') || strcmpi(parameters.calib.type,'
disp(' ')
disp('Creating dark.edf image of calibration scan')
parameters.calib.name=parameters.calib.dir(find(parameters.calib.dir=='/',1,'last')+1:end);
- xmlfname_calib=sprintf('%s/%s.xml',parameters.calib.dir,parameters.calib.name);
+ xmlfname_calib=sprintf('%s/%s.xml',parameters.calib.dir,parameters.calib.name)
if ~exist(xmlfname_calib,'file')
error('Cannot proceed - missing xml file of calibration scan.')
end
diff --git a/3_match_extspot/gtCalculateTwotheta.m b/3_match_extspot/gtCalculateTwotheta.m
index e639ca4fd20964c779d31d4ff449c91d1f1ba7b3..efaa9929927db44098249207c2dec0b40ad59e75 100644
--- a/3_match_extspot/gtCalculateTwotheta.m
+++ b/3_match_extspot/gtCalculateTwotheta.m
@@ -48,7 +48,8 @@ energy=acq.energy;
lambda=gtConvEnergyToWavelength(energy);%mod LNervo
latticepar=acq.latticepar;
-if isfield(acq,'detanglemin') && isfield(acq,'detanglemax')
+if (isfield(acq,'detanglemin') && isfield(acq,'detanglemax')) && ...
+ (~isempty(acq.detanglemin) && ~isempty(acq.detanglemax))
minangle=acq.detanglemin;
maxangle=acq.detanglemax;
else
@@ -279,8 +280,7 @@ if nargout==2
results.reflections=hkl; %from file
results.tt=tt; %from file
results.d0=d0; %from file
- results.twoth=twoth';
- results.dsp=dsp';
+
end
diff --git a/zUtil_OAR/OAR_make.m b/zUtil_OAR/OAR_make.m
index 561cd63d29c2dca114c081410b515fc6013fca95..6f41bf0725c6224b93bc7fe861340f8d4bd39229 100644
--- a/zUtil_OAR/OAR_make.m
+++ b/zUtil_OAR/OAR_make.m
@@ -21,25 +21,27 @@ function OAR_make(executable,first,last,njobs,otherparameters,submitflag,varargi
%
% Version 008 22-09-2011 by LNervo
-% Adjust part of oardirectives (=varargin)
+% Adjust part of oardirectives (=varargin). Passed by pairs
%
-% Available directives:
+% Available directives (default values are shown)
+%
+% NAME VALUE OAR SYNTAX
% interactive false -I / --interactive
-% queue='default';% -q <queue> / --queue=<queue>
-% name=executable;% -n <text> / --name=<text>
-% notify='';% --notify "mail:name\@domain.com"
-% resubmit=-1;% --resumbit=<jobid>
-% xml=false;% -X / --xml
-% stdout='';% -O <outfile> / --stdout=<outfile>
-% stderr='';% -E <errfile> / --stderr=<errfile>
-% version=false;% -V / --version
-% node=1;
-% cpu=1;
-% core=1;
-% gpu='NO';% -p "gpu='NO'"
-% host='';% -p "host like '<name>%'"
-% mem=2000;
-% walltime='2:00:00';
+% queue 'default' -q <queue> / --queue=<queue>
+% name executable -n <text> / --name=<text>
+% notify '' --notify "mail:name\@domain.com"
+% resubmit -1 --resumbit=<jobid>
+% xml false -X / --xml
+% stdout '' -O <outfile> / --stdout=<outfile>
+% stderr '' -E <errfile> / --stderr=<errfile>
+% version false -V / --version
+% node 1 -l /node=<nnode>
+% cpu 1 -l /cpu=<ncpu>
+% core 1 -l /core=<ncore>
+% gpu 'NO' -p gpu = '<gpuvalue>'
+% host '' -p host like '<name>%'
+% mem 2000 -p mem > <memory>
+% walltime '2:00:00' -l walltime=<walltime>
%
% Version 007 29-04-2011 by LNervo:
% Change name of structure parameters into oarparameters
@@ -51,6 +53,8 @@ function OAR_make(executable,first,last,njobs,otherparameters,submitflag,varargi
% clean the screen each time
clc
debug=false;
+
+% default OAR directives
app.interactive=false;% -I / --interactive
app.queue='default';% -q <queue> / --queue=<queue>
app.name=executable;% -n <text> / --name=<text>
@@ -69,6 +73,7 @@ app.mem=2000;
app.walltime='2:00:00';
app=parse_pv_pairs(app,varargin);
+
oarpars=app;
%%%%%%%%%%%%%%
@@ -357,7 +362,7 @@ if strcmpi(matlab,'false')
cmd=sprintf('more %s | grep -m 2 %s ',mpath,fname);
[s,help]=system(cmd);
if isempty(help)
- help=sprintf('# This file has been generated by using MATLAB function OAR_make_v6,\nrunning the module ''%s''\n',executable);
+ help=sprintf('# This file has been generated by using MATLAB function OAR_make_v7,\nrunning the module ''%s''\n',executable);
else
fprintf(fid,'# %s',help);
end
diff --git a/zUtil_Optimization/addMatFile.m b/zUtil_Optimization/addMatFile.m
new file mode 100644
index 0000000000000000000000000000000000000000..1e35ce1f11af3a3546b08087e021f14a5094bb8d
--- /dev/null
+++ b/zUtil_Optimization/addMatFile.m
@@ -0,0 +1,32 @@
+function updateVar=addMatFile(mat1,var,update)
+% updateVar=addMatFile(mat1,var,update)
+
+if ~exist('update','var') || isempty(update)
+ update=false;
+end
+
+names=fieldnames(var);
+
+for i=1:length(names)
+ if ~isfield(mat1,char(names(i)))
+ mat1.(char(names(i)))=var.(char(names(i)));
+ else
+ if ~update
+ variable=char(names(i))
+ oldvalue=mat1.(char(names(i)))
+ newvalue=var.(char(names(i)))
+
+ check=inputwdefault('Do you want to replace it with the new value?','n')
+ if strcmpi(check,'y')
+ mat1.(char(names(i)))=var.(char(names(i)));
+ end
+
+ else
+ mat1.(char(names(i)))=var.(char(names(i)));
+ end
+ end
+end
+
+updateVar=mat1;
+end
+
\ No newline at end of file
diff --git a/zUtil_Optimization/calculateUVWfromGeo.m b/zUtil_Optimization/calculateUVWfromGeo.m
new file mode 100644
index 0000000000000000000000000000000000000000..7339e77c64b3e6bce0e9a64678e73dfe24b2e076
--- /dev/null
+++ b/zUtil_Optimization/calculateUVWfromGeo.m
@@ -0,0 +1,41 @@
+%% Compute uv coordinates and w(omega) values from grain given an arbitrary
+%% geometry
+function fed=calculateUVWfromGeo(geometry,grain)
+
+detdiru = geometry.detdiru0';
+detdirv = geometry.detdirv0';
+Qdet = gtFedDetectorProjectionTensor(detdiru,detdirv,1,1);
+tn = cross(detdiru,detdirv);
+detpos = geometry.detpos0';
+uvorig = [geometry.detucen0, geometry.detvcen0]';
+omstep=180/geometry.nproj;
+
+csam=grain.center';
+
+% grain.allblobs.[variable]
+srot=grain.srot;
+dvec=grain.dvec;
+omega=grain.omega;
+
+fed.uv = [];
+fed.uvw = [];
+
+n=size(omega,1)/4;
+for j=1:n % must be from 0 to n-1
+
+ % from grainFed=gtFedGenerateRandomGRain2(grain{grainid},full)
+ S=srot(12*j-11:12*j,:);
+ D=dvec(4*j-3:4*j,:);
+ O=omega(4*j-3:4*j); % omega values ordered like small, small+180, big, big+180 for the 4
+
+ uv1a=gtFedPredictUVW(S(1:3,:),D(1,:)',csam,detpos,tn,Qdet,uvorig,O(1),omstep);
+ uv1b=gtFedPredictUVW(S(4:6,:),D(2,:)',csam,detpos,tn,Qdet,uvorig,O(2),omstep);
+ uv2a=gtFedPredictUVW(S(7:9,:),D(3,:)',csam,detpos,tn,Qdet,uvorig,O(3),omstep);
+ uv2b=gtFedPredictUVW(S(10:12,:),D(4,:)',csam,detpos,tn,Qdet,uvorig,O(4),omstep);
+
+ fed.uv = [fed.uv;uv1a(1),uv1a(2);uv1b(1),uv1b(2);uv2a(1),uv2a(2);uv2b(1),uv2b(2)];
+ fed.uvw = [fed.uvw;uv1a(1),uv1a(2),uv1a(3);uv1b(1),uv1b(2),uv1b(3);uv2a(1),uv2a(2),uv2a(3);uv2b(1),uv2b(2),uv2b(3)];
+
+end % size /4
+
+end % end function
\ No newline at end of file
diff --git a/zUtil_Optimization/checkSpot.m b/zUtil_Optimization/checkSpot.m
new file mode 100644
index 0000000000000000000000000000000000000000..cb0997655d2edcb82549420f559f69547ca64210
--- /dev/null
+++ b/zUtil_Optimization/checkSpot.m
@@ -0,0 +1,11 @@
+function checkSpot(parameters,bbox)
+
+for i=1:length(bbox.expected_id)
+ vol=gtDBBrowseDiffractionVolume(parameters.acq.name,bbox.expected_id(i));
+ figure;imshow(sum(vol,3),[0 2000])
+ bbox.all(i,:)
+ title(sprintf('#%d id %d x,y [%f %f]',i,bbox.expected_id(i),bbox.xcen(i),bbox.ycen(i)))
+ waitforbuttonpress
+end
+
+end
\ No newline at end of file
diff --git a/zUtil_Optimization/findSpot.m b/zUtil_Optimization/findSpot.m
new file mode 100644
index 0000000000000000000000000000000000000000..59bd33ee64feb3ab35ebed656277790615f98fc9
--- /dev/null
+++ b/zUtil_Optimization/findSpot.m
@@ -0,0 +1,29 @@
+function [spot,range,box,omegamean]=findSpot(parameters,imagemin,imagemax,ymin,ymax,xmin,xmax,number)
+
+dir=parameters.acq.dir;
+cd(sprintf('%s/1_preprocessing/full',dir));
+
+sum=zeros(parameters.acq.xdet);
+for i=imagemin:imagemax
+ im=edf_read(sprintf('full%04d.edf',i));
+ sum=sum+im;
+end
+
+spot=sum(ymin:ymax,xmin:xmax);
+box=[ymin ymax xmin xmax];
+range=[imagemin imagemax];
+
+cd(dir);
+name=sprintf('spot%d.edf',number);
+if exist(name,'file')
+ fprintf('The file %s already exist. ',name)
+ check=inputwdefault('Overwrite it?','y');
+ if strcmpi(check,'y')
+ edf_write(spot,name);
+ end
+else
+ edf_write(spot,name);
+end
+omegamean=mean([imagemin imagemax])*180/parameters.acq.nproj;
+
+end
\ No newline at end of file
diff --git a/zUtil_Optimization/findSpotsWoIndexter.m b/zUtil_Optimization/findSpotsWoIndexter.m
new file mode 100644
index 0000000000000000000000000000000000000000..56ba27943abbb8cb728a8675c232a464e8fb2dde
--- /dev/null
+++ b/zUtil_Optimization/findSpotsWoIndexter.m
@@ -0,0 +1,106 @@
+function bbox=findSpotsWoIndexter(parameters,uvfed,omegafed,bbsizemean,areamin,areamax)
+% bbox=findSpotsWoIndexter(parameters,uvfed,omegafed,bbsizemean,areamin,areamax)
+
+if nargin < 6
+ disp('Usage: bbox=findSpotsWoIndexter(parameters,uvfed,omegafed,bbsizemean,areamin,areamax)')
+ return
+end
+
+if isfield(parameters,'acq')
+ name=parameters.acq.name;
+ omstep=180/parameters.acq.nproj;
+else
+ name=parameters.name;
+ omstep=180/parameters.nproj;
+end
+bbsize=bbsizemean; % pixels
+
+expected_id=[];
+image=[];
+omega=[];
+x0=[];
+y0=[];
+xsize=[];
+ysize=[];
+xcen=[];
+ycen=[];
+feduv=[];
+area=[];
+ind=[];
+dbuv=[];
+
+dif_Xsize=bbsize;
+dif_Ysize=bbsize;
+deltaImage=10;
+
+
+gtDBConnect
+
+for i=1:length(uvfed)
+
+ myquery=sprintf(['select difspotID from %sdifspot where '...
+ '(%d between CentroidImage-%d and CentroidImage+%d) '...
+ 'and (abs(%sdifspot.CentroidX-%f)/%d)<0.5 '...
+ 'and (abs(%sdifspot.CentroidY-%f)/%d)<0.5 '...
+ 'and (Area between %d and %d) '...
+ 'order by ((abs(BoundingBoxXsize-%d)/%d)+(abs(BoundingBoxYsize-%d)/%d)) limit 1'],...
+ name,...
+ round(omegafed(i)/omstep),deltaImage,deltaImage,...
+ name, uvfed(i,1), dif_Xsize,...
+ name, uvfed(i,2), dif_Ysize,...
+ areamin,areamax,...
+ dif_Xsize, dif_Xsize, dif_Ysize, dif_Ysize);
+
+ [difspotID] = mym(myquery);
+
+ if ~isempty(difspotID)
+ expected_id=[expected_id; difspotID];
+ secquery=sprintf('select BoundingBoxXOrigin, BoundingBoxYOrigin, BoundingBoxXSize, BoundingBoxYSize, CentroidX, CentroidY, CentroidImage, Area from %sdifspot where difspotID=%d',...
+ name, difspotID);
+
+ [BoundingBoxXorigin, BoundingBoxYorigin, BoundingBoxXsize, BoundingBoxYsize, CentroidX, CentroidY, CentroidImage, Area]=mym(secquery);
+
+ % disp(sprintf('CentroidImage between CentroidImage-10 and CentroidImage+10 = %d',CentroidImage))
+ % disp(sprintf('abs(CentroidX-uvfed(i,1))/dif_Xsize = %f',abs(CentroidX-uvfed(i,1))/dif_Xsize))
+ % disp(sprintf('abs(CentroidY-uvfed(i,2))/dif_Ysize = %f',abs(CentroidY-uvfed(i,2))/dif_Ysize))
+ % disp(sprintf('Area between %d and %d = %d',areamin,areamax,Area))
+
+
+
+ image=[image; CentroidImage];
+ omega=[omega; CentroidImage*omstep];
+ x0=[x0; BoundingBoxXorigin];
+ y0=[y0; BoundingBoxYorigin];
+ xsize=[xsize; BoundingBoxXsize];
+ ysize=[ysize; BoundingBoxYsize];
+ xcen=[xcen; CentroidX];
+ ycen=[ycen; CentroidY];
+ uvdb=[CentroidX CentroidY];
+ dbuv=[dbuv;uvdb];
+ area=[area; Area];
+ ind=[ind;i];
+ feduv=[feduv;uvfed(i,:)];
+ end
+%
+end
+
+bbox.expected_id=expected_id;
+bbox.x0=x0;
+bbox.y0=y0;
+bbox.xsize=xsize;
+bbox.ysize=ysize;
+bbox.xcen=xcen;
+bbox.ycen=ycen;
+bbox.image=image;
+bbox.omega=omega;
+bbox.area=area;
+bbox.dbuv=dbuv;
+bbox.all=[expected_id x0 y0 xcen ycen xsize ysize omega];
+bbox.title='expected_id x0 y0 xcen ycen xsize ysize omega';
+bbox.query={'CentroidImage',[-deltaImage deltaImage],'Area',[areamin areamax],...
+ 'CentroidX-u', [dif_Xsize 0.5],...
+ 'CentroidY-v', [dif_Ysize 0.5]};
+bbox.ind=ind;
+bbox.feduv=feduv;
+
+end
\ No newline at end of file
diff --git a/zUtil_Optimization/gtCalculateDist.m b/zUtil_Optimization/gtCalculateDist.m
new file mode 100644
index 0000000000000000000000000000000000000000..e1663c20ea2d9d5081657df371e95b0cf471c734
--- /dev/null
+++ b/zUtil_Optimization/gtCalculateDist.m
@@ -0,0 +1,87 @@
+function [dsp,crystal_system,tt]=gtCalculateDist(hkl,parameters)
+% [dsp,crystal_system]=gtCalculateDist(hkl,parameters)
+% hkl <double 1x3>
+
+if nargin < 2
+ disp('Usage: [dsp,crystal]=gtCalculateDist(hkl,parameters)')
+ return
+end
+
+latticepar=parameters.acq.latticepar;
+spacegroup=parameters.acq.spacegroup;
+
+a=latticepar(1);b=latticepar(2);c=latticepar(3);
+alpha=latticepar(4);beta=latticepar(5);gamma=latticepar(6);
+
+if size(hkl,2)==4 %hexagonal materials
+ h=hkl(1);k=hkl(2);l=hkl(4);
+else
+ h=hkl(1);k=hkl(2);l=hkl(3);
+end
+
+if ~isfield(parameters.acq,'hermann_mauguin')
+ hermann=readSpaceGroup(spacegroup);
+end
+
+% calculate twotheta (twoth) for different spacegroups
+if between(spacegroup,1,2)
+ crystal='triclinic';
+ Y=h^2/a^2*(sind(alpha))^2 + k^2/b^2*(sind(beta))^2 + l^2/c^2*(sind(gamma))^2;
+ Z=2*k*l/b/c*(cosd(beta)*cosd(gamma)-cosd(alpha)) + ...
+ 2*l*h/c/a*(cosd(gamma)*cosd(alpha)-cosd(beta)) + ...
+ 2*h*k/a/b*(cosd(alpha)*cosd(beta)-cosd(gamma));
+ invdsp2=1/(1-(cosd(alpha))^2-(cosd(beta))^2-(cosd(gamma))^2+2*cosd(alpha)*cosd(beta)*cosd(gamma)) * ( Y + Z );
+end
+if between(spacegroup,3,15)
+ crystal='monoclinic';
+ invdsp2=h^2/(a*sind(beta))^2 + k^2/b^2 + l^2/(c*sind(beta))^2 - (2*h*l*cosd(beta)) / (a*c*(sind(beta))^2);
+end
+if between(spacegroup,16,74)
+ crystal='orthorhombic';
+ invdsp2=h^2/a^2+k^2/b^2+l^2/c^2;
+end
+if between(spacegroup,75,142)
+ crystal='tetragonal';
+ invdsp2=(h^2+k^2)/a^2+l^2/c^2;
+end
+if between(spacegroup,143,167)
+ crystal='trigonal';
+ if ~isempty(strfind(hermann,'P')) % hexagonal system
+ invdsp2=4/(3*a^2) * (h^2 + k^2 + h*k) + l^2/c^2;
+ elseif ~isempty(strfind(hermann,'R')) % rhombohedral system
+ Y1=h^2+k^2+h*k;
+ Y2=2*(h*k+h*l+k*l);
+ Z1=(sind(alpha))^2;
+ Z2=(cosd(alpha))^2-cosd(alpha);
+ W=1+2*(cosd(alpha))^3-3*(cosd(alpha))^2;
+ invdsp2=1/a^2*(Y1*Z1+Y2*Z2)/W;
+ end
+end
+if between(spacegroup,168,194) || spacegroup==663 % snow case
+ crystal='hexagonal';
+ invdsp2=4/(3*a^2) * (h^2 + k^2 + h*k) + l^2/c^2;
+end
+if between(spacegroup,195,230)
+ crystal='cubic';
+ invdsp2=(h^2+k^2+l^2)/a^2;
+end
+dsp=1/(invdsp2)^0.5;
+lambda=gtConvEnergyToWavelength(parameters.acq.energy);
+twoth=2*asind(lambda/(2*dsp));
+
+if nargout == 2
+ crystal_system=crystal;
+end
+
+if nargout == 3
+ crystal_system=crystal;
+ tt=twoth;
+end
+
+if ~isfield(parameters.acq,'crystal_system')
+ parameters.acq.crystal_system=crystal;
+ save parameters parameters
+end
+
+
+end
\ No newline at end of file
diff --git a/zUtil_Optimization/gtCalculateUVWfromGeo.m b/zUtil_Optimization/gtCalculateUVWfromGeo.m
new file mode 100644
index 0000000000000000000000000000000000000000..7063046746e9071e6fed716e66ae25573c9e354b
--- /dev/null
+++ b/zUtil_Optimization/gtCalculateUVWfromGeo.m
@@ -0,0 +1,38 @@
+%% Compute uv coordinates and w(omega) values from grain given an arbitrary
+%% geometry
+function out=gtCalculateUVWfromGeo(geometry,grainFed)
+
+detdiru = geometry.detdiru0';
+detdirv = geometry.detdirv0';
+Qdet = gtFedDetectorProjectionTensor(detdiru,detdirv,1,1);
+tn = cross(detdiru,detdirv);
+detpos = geometry.detpos0';
+uvorig = [geometry.detucen0, geometry.detvcen0]';
+omstep=180/geometry.nproj;
+
+csam=grainFed.center';
+
+fed.uv = [];
+fed.uvw = [];
+
+n=size(grainFed.allblobs.omega,1)/4;
+for j=1:n % must be from 0 to n-1
+
+ % from grainFed=gtFedGenerateRandomGRain2(grain{grainid},full)
+ S=grainFed.allblobs.srot(12*j-11:12*j,:);
+ D=grainFed.allblobs.dvec(4*j-3:4*j,:);
+ O=grainFed.allblobs.omega(4*j-3:4*j); % omega values ordered like small, small+180, big, big+180 for the 4
+
+ uv1a=gtFedPredictUVW(S(1:3,:),D(1,:)',csam,detpos,tn,Qdet,uvorig,O(1),omstep);
+ uv1b=gtFedPredictUVW(S(4:6,:),D(2,:)',csam,detpos,tn,Qdet,uvorig,O(2),omstep);
+ uv2a=gtFedPredictUVW(S(7:9,:),D(3,:)',csam,detpos,tn,Qdet,uvorig,O(3),omstep);
+ uv2b=gtFedPredictUVW(S(10:12,:),D(4,:)',csam,detpos,tn,Qdet,uvorig,O(4),omstep);
+
+ fed.uv = [fed.uv;uv1a(1),uv1a(2);uv1b(1),uv1b(2);uv2a(1),uv2a(2);uv2b(1),uv2b(2)];
+ fed.uvw = [fed.uvw;uv1a(1),uv1a(2),uv1a(3);uv1b(1),uv1b(2),uv1b(3);uv2a(1),uv2a(2),uv2a(3);uv2b(1),uv2b(2),uv2b(3)];
+
+end % size /4
+
+out=fed;
+
+end % end function
\ No newline at end of file
diff --git a/zUtil_Optimization/gtCheckParameters.m b/zUtil_Optimization/gtCheckParameters.m
new file mode 100644
index 0000000000000000000000000000000000000000..8fb0ce4f66443c4fe7697077aaf748e3ab05e1a7
--- /dev/null
+++ b/zUtil_Optimization/gtCheckParameters.m
@@ -0,0 +1,185 @@
+function gtCheckParameters(parameters)
+
+field{1,1}='detpos0'; field{1,2}='Detector position';
+field{2,1}='detdiru0'; field{2,2}='Detector u-direction';
+field{3,1}='detdirv0'; field{3,2}='Detector v-direction';
+field{4,1}='beamdir0'; field{4,2}='Beam direction';
+field{5,1}='rotdir0'; field{5,2}='Rotation axis direction';
+field{6,1}='detucen0'; field{6,2}='Detector u-centre';
+field{7,1}='detvcen0'; field{7,2}='Detector v-centre';
+field{8,1}='detdir0'; field{8,2}='Detector direction';
+field{9,1}='hermann_mauguin'; field{9,2}='Hermann Mauguin short symbol';
+field{10,1}='xop'; field{10,2}='Output from xop/Diamond';
+field{11,1}='crystal_system'; field{11,2}='crystal system';
+field{12,1}='detanglemin'; field{12,2}='Detector minimum 2Theta';
+field{13,1}='detanglemax'; field{13,2}='Detector maximum 2Theta';
+field{14,1}='no_direct_beam'; field{14,2}='No beam into the detector';
+field{15,1}='rotation_axis'; field{15,2}='rotation axis';
+field{16,1}='dist'; field{16,2}='distance sample-detector';
+field{16,1}='name'; field{16,2}='name of the dataset';
+
+if ~exist('paraemters','var')
+ load('parameters.mat');
+end
+
+names=fieldnames(parameters.acq);
+
+for i=1:length(names)
+ list{i,1}=names{i};
+end
+
+if all(isfield(parameters.acq,field(:,1)))
+ % check if values are correct
+
+ correct=true;
+ %%% detpos %%%
+ if norm(parameters.acq.detpos0)~=parameters.acq.dist
+ disp('The norm of detpos0 and dist are different. They may be corrupted.')
+ correct=false;
+ end
+
+ %%% vertical detector %%%
+ if strcmpi(parameters.acq.no_direct_beam,'f') || ...
+ strcmpi(parameters.acq.no_direct_beam,'false') || ...
+ strcmpi(parameters.acq.no_direct_beam,'no')
+ parameters.acq.no_direct_beam=false;
+ save parameters parameters
+ end
+ if strcmpi(parameters.acq.no_direct_beam,'t') || ...
+ strcmpi(parameters.acq.no_direct_beam,'true') || ...
+ strcmpi(parameters.acq.no_direct_beam,'yes')
+ parameters.acq.no_direct_beam=true;
+ save parameters parameters
+ end
+ if parameters.acq.no_direct_beam || strcmpi(parameters.acq.rotation_axis,'horizontal')
+ if strcmpi(parameters.acq.rotation_axis,'vertical')
+ disp('The rotation axis should not be vertical in this case.')
+ correct=false;
+ end
+ if parameters.acq.rotdir0(3)~=0
+ disp('The rotation axis should be the horizontal plane. Check the third component.')
+ correct=false;
+ end
+ if (parameters.acq.no_direct_beam && strcmpi(parameters.acq.rotation_axis,'vertical')) || ...
+ (~parameters.acq.no_direct_beam && strcmpi(parameters.acq.rotation_axis,'horizontal'))
+ disp('Check the rotation axis and the direct beam, which may be corrupted.')
+ correct=false;
+ end
+
+ end
+
+ %%% beam %%%
+ if parameters.acq.beamdir0~=[1 0 0]
+ disp('The beam direction may be corrupted.')
+ parameters.acq.beamdir0=[1 0 0];
+ save parameters parameters
+ disp('beam direction corrected')
+ end
+
+ %%% detdir0 %%%
+ if parameters.acq.detdir0~=cross(parameters.acq.detdiru0,parameters.acq.detdirv0)
+ disp('The detector direction may be corrupted. Check the u,v directions too.')
+ parameters.acq.detdir0=cross(parameters.acq.detdiru0,parameters.acq.detdirv0);
+ disp('detector direction corrected.')
+ end
+
+ %%% angles %%%
+ [min,max]=gtLimitsTwoTheta(0,1)
+ if parameters.acq.detanglemin~=min
+ disp('detanglemin may be corrupted.')
+ [min,max]=gtLimitsTwoTheta(1,1);
+ min
+ disp('detanglemin corrected.')
+ end
+ if parameters.acq.detanglemax~=max
+ disp('detanglemax may be corrupted.')
+ [min,max]=gtLimitsTwoTheta(1,1);
+ max
+ disp('detanglemax corrected.')
+ end
+
+ %%% latticepar %%%
+ if strcmpi(parameters.acq.crystal_system,'cubic')
+ lat=parameters.acq.latticepar;
+ if ~isequal(lat(1),lat(2),lat(3)) || ~isequal(lat(4),lat(5),lat(6))
+ disp('Lattice parameters must be [a a a alpha alpha alpha]')
+ correct=false;
+ end
+ end
+ if strcmpi(parameters.acq.crystal_system,'hexagonal')
+ lat=parameters.acq.latticepar;
+ if ~isequal(lat(1),lat(2)) || ~isequal(lat(4),lat(5))
+ disp('Lattice parameters must be [a a c alpha alpha gamma]')
+ correct=false;
+ end
+ end
+
+ %%% name %%%
+ count=0;
+ for i=1:length(names)
+ if ischar(parameters.acq.(char(names(i))))
+ ind=strfind(char(parameters.acq.(char(names(i)))),parameters.acq.name);
+ if ~isempty(ind)
+ count=count+1;
+ end
+ end
+ end
+ if count~=6
+ disp('The name of one or more directory or table is corrupted.')
+ correct=false;
+ end
+
+ if ~correct
+ disp('Check the values of the parameters...')
+ parameters.acq=gtModifyStructure(parameters.acq, list,'Display parameters.acq:');
+ save parameters parameters
+ %gtCheckParameters
+ else
+ disp('Parameters are corrected.')
+ end
+
+else
+ f=isfield(parameters.acq,field(:,1));
+ field(:,2)=num2cell(f);
+ for i=1:length(field)
+ if isfield(parameters.acq,char(field(i,1)))
+ name=parameters.acq.(char(field(i,1)));
+ if isvector(name)
+ field(i,3)={name};
+ elseif isfloat(name)
+ field(i,3)=num2cell(name);
+ else
+ field(i,3)=cellstr(name);
+ end
+ else
+ field{i,3}='';
+ fprintf('The variable %s is missing...\n',char(field(i,1)))
+ missing.(char(field(i,1)))='';
+ end
+ end
+ disp(field)
+ disp(' ')
+ disp(missing)
+
+ newnames=fieldnames(missing);
+
+ for l=1:length(newnames)
+ parameters.acq.(char(newnames(l)))='';
+ newlist{l,1}=newnames{l};
+ for k=1:size(field,2)
+ if strcmpi(newnames{l},field{k,2})
+ newlist{l,2}=field{k,2};
+ else
+ newlist{l,2}='';
+ end
+ end
+ end
+
+ disp('Add the values of the missing parameters...')
+ parameters.acq = gtModifyStructure(parameters.acq, newlist,'Missing parameters.acq:');
+ save parameters parameters
+ %gtCheckParameters
+
+end
+
+end
\ No newline at end of file
diff --git a/zUtil_Optimization/gtFindUVW3.m b/zUtil_Optimization/gtFindUVW3.m
new file mode 100644
index 0000000000000000000000000000000000000000..71988b7bfa979233e6ba46120714cb2288e0c328
--- /dev/null
+++ b/zUtil_Optimization/gtFindUVW3.m
@@ -0,0 +1,85 @@
+function [F,M,x] = gtFindUVW3(x,xdata,csam,uvorig)
+% detdiru (3x1)
+% detdirv (3x1)
+% detpos (3x1)
+% uvorig (2x1)
+% x (1x9)
+%
+% xdata (Nx4,3)
+% srot (Nx3,3)
+% dvec (N,3)
+
+if size(x,1) ~= 1 && size(x,2) ~= 9
+ disp('Please check the size of x')
+ return
+end
+
+%detdiru=[x(1);sqrt(1-x(1)^2-x(2)^2);x(2)];
+%detdirv=[x(3);x(4);sqrt(1-x(3)^2-x(4)^2)];
+%detpos=[x(5);x(6);x(7)];
+
+detpos=[x(7);x(8);x(9)];
+
+% Gram - Schmidt process for orthonormalization
+
+x1=[x(1);x(2);x(3)];
+x2=[x(4);x(5);x(6)];
+x3=cross(x1,x2);
+norm12=dot(x1,x2); % alpha
+norm23=dot(x2,x3); % beta
+norm13=dot(x1,x3); % gamma
+
+%disp('Gram-Schmidt process')
+y=orthonormGS(x1,x2,x3);
+x(1:6)=[y(1,:),y(2,:)];
+% new detector reference system
+M=y;
+
+
+detdiru=y(1,:)';
+normu=norm(detdiru);
+detdirv=y(2,:)';
+normv=norm(detdirv);
+tn=y(3,:)';
+normt=norm(tn);
+
+% normalize
+detdiru=detdiru/norm(detdiru);
+detdirv=detdirv/norm(detdirv);
+
+% split the xdata into two matrices
+n3=size(xdata,1)/4*3;
+srot=xdata(1:n3,:);
+dvec=xdata(n3+1:end,:);
+
+detscaleu=1;
+detscalev=1;
+
+% Q - projection matrix from LAB to DETECTOR coordinates %2x3
+Qdet = [ detdiru'*detscaleu; detdirv'*detscalev];
+
+F=[];
+n=size(dvec,1);
+for j=1:n
+ % rotation tensor
+ S = srot(3*j-2:3*j,:); %3x3
+ % diffraction vector
+ D = dvec(j,:)'; %3x1 unit vector
+
+ %grain coord. in LAB system
+ clab = S*csam; %3x1
+ % diffraction factor
+ dfac = tn'*(detpos-clab)/(tn'*D); %1x1
+
+ % diffraction vector in the u,v plane from uvorig
+ duvpl = clab + dfac*D - detpos; %3x1
+
+ uv = uvorig + Qdet*duvpl; %2x1
+
+ F = [F; uv(1),uv(2)]; % Nx2
+
+end
+
+
+
+end % end function
\ No newline at end of file
diff --git a/zUtil_Optimization/gtForwardSimulate2_ln.m b/zUtil_Optimization/gtForwardSimulate2_ln.m
new file mode 100644
index 0000000000000000000000000000000000000000..8782fc02a4f304dfea4ad0e6fc4cf87517207117
--- /dev/null
+++ b/zUtil_Optimization/gtForwardSimulate2_ln.m
@@ -0,0 +1,830 @@
+function [grain,o] = gtForwardSimulate2_ln(grain,parameters,varargin)
+%% GTFORWARDSIMULATE2_LN Forward simulation
+% [grain,o] = gtForwardSimulate2_ln(grain,parameters,varargin)
+%
+% should be integrated into gtINDEXTER later
+%
+% - searches database for missing reflections ...
+% - identify new pairs
+% - calculate theta,eta,omega,hkl,pn, et... for each of the new spot
+% - update grain structure
+% - update spotpair table
+% - update difspot table
+%
+% WL 04/2009: overworked previous version gtForwardSimulate_360
+%
+% INPUT:
+% grain{grainid}
+% parameters
+%
+%
+% OPTIONAL INPUT:
+% update_db (<logical 1x1> {0})
+% show_results (<logical 1x1> {0})
+% connectdb (<logical 1x1> {1})
+%
+%
+% OUTPUT:
+% grain{grainid} updated with additional fields:
+% - newdifspots difspotID of added spots
+% - newpairs pairID of added pairs
+% - completeness ratio of theoretically expected over actually found reflections
+% - missing_refl contains information concerning those reflections
+% which could not be found in the database
+%
+% o (output) has fields:
+% - image : image number
+% - omega : omega value in degrees
+% - udirect : u (horiz.) coordinate of extspot
+% - vdirect : v (vertical) coordinate of extspot
+% - udiffr : u (horiz.) coordinate of difspot
+% - vdiffr : v (vertical) coordinate of difspot
+% - omega : omega in degrees
+% - eta : eta in degrees
+% - theta : theta in degrees
+% - h,k,l : (signed) hkl indices of the reflection
+% - flag : NaN: reflection geometrically not accessible
+% -1: no spot matching the search criteria has been found, although expected (overlaps, segmentation problem, etc...)
+% 0: spot position outside sensitive area of the detector
+% 1: new spot to be included in calculation of completeness
+% 2: spot pair which has already been assigned to the grain by gtINDEXTER
+% 4: spotpair which has been matched, but not assigned to the grain during Indexing
+% 5: problem
+% - difspotID : difspotID...
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+if nargin < 2
+ disp('Usage: [grain,o] = gtForwardSimulate2_ln(grain,parameters,varargin)')
+ grain=[];
+ o=[];
+ return
+end
+app.connectdb=1;
+app.update_db=0;
+app.show_results=0;
+app=parse_pv_pairs(app,varargin);
+
+
+
+% new fields
+grain.missing_refl=[]; % will be used to keep track of reflections which could not be found in the database (e.g. segmentation problem, weak intensity, overlaps...)
+grain.newdifspots=[]; % difspotID of spots which have been added to the grain
+grain.nof_newpairs=[]; % pairID of pairs identified in newdifspots
+grain.completeness=[]; % ratio of found over expected reflections
+
+% connect to the database
+if app.connectdb
+ gtDBConnect
+end
+
+
+%% read parameters
+
+acq=parameters.acq;
+name=acq.name;
+spacegroup=acq.spacegroup;
+difspottable=[name,'difspot'];
+db=parameters.database;
+latticepar=acq.latticepar;
+
+%get the mean difspot, X and Y sizes for the current spots in the grain
+dif_Xsize=grain.stat.bbxsmean;
+dif_Xsize_ratio=grain.stat.bbxsrat*1.3;
+dif_Ysize=grain.stat.bbysmean;
+dif_Ysize_ratio=grain.stat.bbysrat*1.3;
+
+dif_Xsearch=2*grain.stat.bbxsstd;
+dif_Ysearch=4*grain.stat.bbysstd;
+
+
+% Deal with the problem that old versions of gtMatchDifspots did not add fields
+% etaB and omegaB to the spotpair_table
+[a,b,c,d,e,f]=mym(sprintf('show columns from %s like "etaB"',acq.pair_tablename));
+if isempty(a)
+ %add fields etaB and omegaB to old tables
+ mym(sprintf('ALTER TABLE %s.%s ADD COLUMN etaB FLOAT AFTER omega',db.name,acq.pair_tablename));
+ mym(sprintf('ALTER TABLE %s.%s ADD COLUMN omegaB FLOAT AFTER etaB',db.name,acq.pair_tablename));
+ disp('adding etaB and omegaB fields to the spotpair table')
+ % here we need to run through the table and fill in these two fields -
+ % maybe we should add also the hkl and plane normal information
+ % this still needs to be done, sorry.
+end
+
+if ~isfield(grain.stat,'areamean')
+ % we need to add this value in gtINDEXTER
+ grain.stat.areamean=0.7*grain.stat.bbxsmean*grain.stat.bbysmean; % 0.7 is a rough estimate of the spot area / BB-area ratio
+end
+dif_int=grain.stat.intmean;
+dif_area=grain.stat.areamean;
+dif_area_ratio=dif_Xsize_ratio*dif_Ysize_ratio;
+
+
+%% get all possible hkls
+
+[hkltypes_used,info]=gtCrystSpacegroupReflections;
+
+disp('Simulated HKL reflections:')
+disp(hkltypes_used)
+disp(['hkltypes_used: ' num2str(size(hkltypes_used))])
+
+
+% 4 columns needed for hexagonal
+if strcmpi(acq.crystal_system,'hexagonal')
+ hkltypes_used(:,4) = hkltypes_used(:,3);
+ hkltypes_used(:,3) = -hkltypes_used(:,1)-hkltypes_used(:,2);
+end
+
+hklsp = [];
+hkl = [];
+thtype = [];
+
+for i = 1:size(hkltypes_used,1)
+ allhkls_i = gtCrystSignedHKLs(hkltypes_used(i,:),spacegroup);
+ nhkls = size(allhkls_i,1)/2;
+ hklsp = [hklsp; allhkls_i(nhkls+1:end,:)];
+ hkl = [hkl; repmat(hkltypes_used(i,:),nhkls,1)];
+ thtype = [thtype; repmat(i,nhkls,1)];
+end
+
+%disp(['hkl: ' num2str(size(hkl))])
+%disp(['allhkls_i: ' num2str(size(allhkls_i))])
+disp(['hklsp: ' num2str(size(hklsp))])
+
+disp('Getting cartesian coordinates from lattice coordinates...')
+for i=1:size(hklsp,1)
+ normalised_hkls(i,:)=gtHKL2Cart_ln(hklsp(i,:),spacegroup); % modified by LNervo
+end
+
+% produce all possible plane normals using R_vector
+g = Rod2g(grain.R_vector);
+% warning('using g, rather than inv(g), in eq. 3.6 in Poulsen - Rod2g or coords problem?')
+all_normals = (g * normalised_hkls')';
+all_hkls=hklsp;
+%disp(['all_normals: ' num2str(size(all_normals))])
+
+grain.hklsp=hklsp;
+grain.all_normals=all_normals;
+
+%% Forward simulation
+% now do forward simuation - where do we expect dif and ext spots ?
+disp('Starting forward simulation...')
+
+j=1;
+for i=1:size(all_normals,1)
+
+ [image,udirect,vdirect,udiffr,vdiffr,omega,eta,theta] = gtPredictPositions_3D_ln(all_normals(i,:), all_hkls(i,:), grain.center, 0, parameters, 1);
+
+ len=length(image);
+ h=repmat(all_hkls(i,1),len,1);
+ k=repmat(all_hkls(i,2),len,1);
+
+ if size(all_hkls,2)==4
+ i_hkl=repmat(all_hkls(i,3),len,1);
+ l=repmat(all_hkls(i,4),len,1);
+ else
+ l=repmat(all_hkls(i,3),len,1);
+ i_hkl=NaN;
+ i_hkl=repmat(i_hkl,len,1);
+ end
+
+
+ if isnan(image)
+ o.flag(j)=NaN;
+ o.h(j)=h;
+ o.k(j)=k;
+ o.i(j)=i_hkl;
+ o.l(j)=l;
+ o.ind(j)=i;
+ disp(sprintf('Reflection %s is not accessible',num2str(all_hkls(i,:))))
+% if size(all_hkls,2)==3
+% disp(sprintf('Reflection %d %d %d is not accessible',all_hkls(i,1),all_hkls(i,2),all_hkls(i,3)));
+% else
+% disp(sprintf('Reflection %d %d %d %d is not accessible',all_hkls(i,1),all_hkls(i,2),all_hkls(i,3),all_hkls(i,4)));
+% end
+ j=j+1;
+ else
+ for ind=1:len
+ o.image(j) = image(ind);
+ o.udirect(j) = udirect(ind);
+ o.vdirect(j) = vdirect(ind);
+ o.udiffr(j) = udiffr(ind);
+ o.vdiffr(j) = vdiffr(ind);
+ o.omega(j) = omega(ind);
+ o.eta(j) = eta(ind);
+ o.theta(j) = theta(ind);
+ o.h(j) = h(ind);
+ o.k(j) = k(ind);
+ o.i(j) = i_hkl(ind);
+ o.l(j) = l(ind);
+ o.flag(j) = 1;
+ o.difspotID(j) = 0;
+ o.ind(j) = i;
+ j=j+1;
+ end
+ end
+
+end
+
+nop=grain.nof_pairs;
+grain.hklA=zeros(nop,4);
+
+%% Interogate the database(s) to find predicted spots
+
+
+% variables expected and found will be used for calculation of completeness
+expected=0;
+found=0;
+radius1=sqrt((dif_Xsize/2)^2+(dif_Ysize/2)^2); % half diagonal of average difspot BoundingBox
+
+for i=1:length(o.image)
+ radius2=sqrt((o.udiffr(i)-acq.xdet/2)^2+(o.vdiffr(i)-acq.ydet/2)^2); % distance between center and spot
+ radius=radius1+radius2; % approximate (maximum) radius of the difspot
+ % don't look at spots which potentially touch the blind area or direct beam area of the detector
+ if ( o.udiffr(i) < dif_Xsize/2 || o.udiffr(i) > acq.xdet-dif_Xsize/2 ...
+ || o.vdiffr(i) < dif_Ysize/2 || o.vdiffr(i) > acq.ydet-dif_Ysize/2 ...
+ || radius > acq.maxradius ...
+ || ( o.udiffr(i)+dif_Xsize/2>parameters.seg.bb(1) && o.udiffr(i)-dif_Xsize/2 < parameters.seg.bb(1)+parameters.seg.bb(3) ...
+ && o.vdiffr(i)+dif_Ysize/2 > parameters.seg.bb(2) && o.vdiffr(i)-dif_Ysize/2 < parameters.seg.bb(2)+parameters.seg.bb(4) ) ...
+ )
+ o.flag(i)=0; % don't consider these reflections
+ if isnan(o.i(i))
+ disp(sprintf('%d %d %d spot close to border or outside of active area - skipping...',o.h(i),o.k(i),o.l(i)));
+ else
+ disp(sprintf('%d %d %d %d spot close to border or outside of active area - skipping...',o.h(i),o.k(i),o.i(i),o.l(i)));
+ end
+ continue % skip this reflection and go to next one
+ else
+ expected=expected+1;
+ % use selection criteria similar to match difspots (match)
+ match=parameters.match;
+ thr_max_offset=min(20,round(match.thr_max_offset/(abs(sind(o.eta(i)))))); % allow for a 1/sin(eta) dependency on omega search range - limited to a max offset of 20 images
+ firstquery=sprintf(['select difspotID, CentroidX, CentroidY, CentroidImage, BoundingBoxXsize, BoundingBoxYsize, Area, grainID from %sdifspot where '...
+ '(%d between CentroidImage-%d and CentroidImage+%d) '...
+ 'and (abs(%sdifspot.CentroidX-%d)/%d)<0.5 '...
+ 'and (abs(%sdifspot.CentroidY-%d)/%d)<0.5 '...
+ 'order by ((abs(BoundingBoxXsize-%d)/%d)+(abs(BoundingBoxYsize-%d)/%d)) limit 1'],...
+ name,...
+ o.image(i), thr_max_offset, thr_max_offset,...
+ name, o.udiffr(i), dif_Xsize,...
+ name, o.vdiffr(i), dif_Ysize,...
+ dif_Xsize, dif_Xsize, dif_Ysize, dif_Ysize); % use relaxed conditions for the COM positions
+
+ [candidates, CentroidX, CentroidY, CentroidImage, BoundingBoxXsize, BoundingBoxYsize, Area, grainID] = mym(firstquery);
+
+ if isempty(candidates) % segmentation problem ?
+
+ o.flag(i)=-1;
+ if isnan(o.i(i)) % get the display right for hkil notation in hexagonal, trigonal, etc...
+ disp(sprintf('Can not find %d %d %d reflection (segmentation problem ?) skipping...',o.h(i),o.k(i),o.l(i)));
+ grain.missing_refl=[grain.missing_refl;o.h(i),o.k(i),o.l(i)];
+
+ else
+ disp(sprintf('Can not find %d %d %d %d reflection (segmentation problem ?) skipping...',o.h(i),o.k(i),o.i(i),o.l(i)));
+ grain.missing_refl=[grain.missing_refl;o.h(i),o.k(i),o.i(i),o.l(i)];
+ end
+
+ if app.show_results
+ disp(firstquery)
+
+ %convert Omega to image number
+ start_image = o.image(i)-5;
+ if start_image<0
+ start_image=0;
+ end
+ end_image = o.image(i)+5;
+ if end_image>(2*acq.nproj)-1;
+ end_image=(2*acq.nproj)-1;
+ end
+
+ start_image
+ end_image
+
+ full=zeros(acq.ydet,acq.xdet);
+ 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%
+
+ figure(1)
+ imshow((full),[-100 100]);
+ impixelinfo
+ text=sprintf('%d %d %d reflection',o.h(i),o.k(i), o.l(i));
+ title(text)
+ hold on;
+ plot(o.udirect(i),o.vdirect(i),'go')
+ plot(o.udiffr(i),o.vdiffr(i),'ro')
+ disp('click image to close and check next reflection')
+ %keyboard
+ waitforbuttonpress;
+ end
+
+ else
+ found=found+1;
+
+ % we assume that part of the intensity stems from our grain and increment the completeness counter...
+ % now refine the search criteria
+ query=sprintf(['select difspotID from %sdifspot where '...
+ '(%d between CentroidImage-%d and CentroidImage+%d) '...
+ 'and (abs(%sdifspot.CentroidX-%d)/%d)<%f '...
+ 'and (abs(%sdifspot.CentroidY-%d)/%d)<%f '...
+ 'and (BoundingBoxXsize between %f and %f) and (BoundingBoxYsize between %f and %f) '...
+ 'and (Area between %f and %f) '...
+ 'and (grainID=0 or grainID=%d) '... % only look at spots which have not already been assigned
+ 'order by ((abs(BoundingBoxXsize-%d)/%d)+(abs(BoundingBoxYsize-%d)/%d)) limit 1'],... % this could be improved according to Peters goodness criteria...
+ name,...
+ o.image(i), thr_max_offset, thr_max_offset,...
+ name, o.udiffr(i), dif_Xsize, dif_Xsearch/dif_Xsize,...
+ name, o.vdiffr(i), dif_Ysize, dif_Xsearch/dif_Xsize,... % use the same (more relaxed) X condition for the Y COM position
+ dif_Xsize/dif_Xsize_ratio, dif_Xsize*dif_Xsize_ratio, dif_Ysize/dif_Ysize_ratio, dif_Ysize*dif_Ysize_ratio,...
+ dif_area/dif_area_ratio, dif_area*dif_area_ratio,...
+ grain.id,...
+ dif_Xsize, dif_Xsize, dif_Ysize, dif_Ysize);
+
+ [difspotID] = mym(query);
+
+ % classify spots with flag:
+ % NaN: not accessible reflection
+ % -1 : no intensity found at expected position (loose criteria)
+ % 0 : there is some intensity - but the spot does not match strictsearch criteria
+ % 1 : spot matching the strict criteria
+ % 2 : spot which already had been found by gtINDEXTER
+
+ if isempty(difspotID)
+ disp(sprintf('Found some intensity for the %d %d %d reflection, but the spot does not match the strict search criteria - excluding this spot...)...',o.h(i),o.k(i),o.l(i)));
+
+ o.flag(i)=0;
+
+ if app.show_results
+ sprintf(['CentroidImage: %d is expected between %d and %d\n'...
+ 'CentroidX: %f is expected between %f and %f \n '...
+ 'CentroidY: %f is expected between %f and %f \n' ...
+ 'BoundingBoxXsize: %d is expected between %f and %f \n' ...
+ 'BoundingBoxYsize: %d is expected between %f and %f \n' ...
+ 'Area: %f: is expected between %f and %f \n' ...
+ 'grainID: %d should be 0 or %d'],...
+ CentroidImage, CentroidImage-thr_max_offset, CentroidImage+thr_max_offset,...
+ CentroidX,o.udiffr(i)-dif_Xsearch,o.udiffr(i)+dif_Xsearch,...
+ CentroidY,o.vdiffr(i)-dif_Ysearch,o.vdiffr(i)+dif_Ysearch,...
+ BoundingBoxXsize,dif_Xsize/dif_Xsize_ratio, dif_Xsize*dif_Xsize_ratio,...
+ BoundingBoxYsize,dif_Ysize/dif_Ysize_ratio, dif_Ysize*dif_Ysize_ratio,...
+ Area, dif_area/dif_area_ratio, dif_area*dif_area_ratio,...
+ grainID,grainID)
+
+
+ canspot=gtGetSummedDifSpot(candidates,parameters);
+ figure(2);imshow(canspot,[])
+ %convert Omega to image number
+ start_image = o.image(i)-5;
+ if start_image<0
+ start_image=0;
+ end
+ end_image = o.image(i)+5;
+ if end_image>(2*acq.nproj)-1;
+ end_image=(2*acq.nproj)-1;
+ end
+
+ start_image
+ end_image
+
+ full=zeros(acq.ydet,acq.xdet);
+ 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%
+
+ figure(1)
+ imshow((full),[-100 100]);
+ impixelinfo
+ text=sprintf('%d %d %d reflection',o.h(i),o.k(i), o.l(i));
+ title(text)
+ hold on;
+ plot(o.udirect(i),o.vdirect(i),'go')
+ plot(o.udiffr(i),o.vdiffr(i),'ro')
+ disp('click image to close and check next reflection')
+ %keyboard
+ k=waitforbuttonpress;
+ end
+ % elseif length(difspotID>1)
+ % select the most probable spot...to be implemented
+ else
+ disp(sprintf('Found spot matching the strict criteria for the %d %d %d reflection',o.h(i),o.k(i),o.l(i)));
+ o.flag(i)=1;
+ o.difspotID(i)=difspotID; % this spot has been detected
+ index=find(grain.difspots==difspotID); % do we already have a difspot with this ID number?
+
+ if ~isempty(index)
+ o.flag(i)=2; % this spot had already been paired and indexed by gtINDEXTER
+ if index>nop % this spot is in the second half of the scan - we can update the hkl information of the pair spot in the first half...
+ index=index-nop;
+ grain.hklA(index,:)=-[o.h(i), o.k(i), o.i(i), o.l(i)]; % the pair spot has negative hkil indices
+ else
+ grain.hklA(index,:)=[o.h(i), o.k(i), o.i(i), o.l(i)];
+ end
+
+ if app.update_db
+ mysqlcmd=dbUpdate(difspottable,'difspotID',difspotID,'GrainID',grain.id);
+ mym(mysqlcmd);
+ end
+ end
+
+ end % if a difspot has been found at the expected position
+ end % if isempty candidates
+ end % if there could be a spot on the detector
+end % loop over predicted spot positions
+
+
+grain.completeness=found/expected;
+
+%% Are there any Friedel pairs in the remaining list of detected spots ?
+%
+% Pairs will be processed as such (calculation of
+% eta, theta, omega from the Friedel pair...
+difspotID=[];
+j=1;
+for i=1:length(o.flag)
+ if o.flag(i)==1
+ ind(j)=i; % ind indicates lines in output structure containing unpaired spots
+ omega(j)=o.omega(i);
+ eta(j)=o.eta(i);
+ theta(j)=o.theta(i);
+ hkl(j,1:4)=[o.h(i),o.k(i),o.i(i),o.l(i)];
+ difspotID(j)=o.difspotID(i);
+ j=j+1;
+ end
+end
+
+
+
+grain.newdifspots=difspotID;
+
+grain.hklB=zeros(nop,4);
+grain.plA=zeros(nop,3);
+grain.plB=zeros(nop,3);
+%nop=grain.nof_pairs;
+grain.difspotA=zeros(nop,1);
+grain.difspotB=zeros(nop,1);
+grain.omegaA=zeros(nop,1);
+grain.omegaB=zeros(nop,1);
+grain.etaA=zeros(nop,1);
+grain.etaB=zeros(nop,1);
+grain.thetaA=zeros(nop,1);
+grain.thetaB=zeros(nop,1);
+
+% define new variables
+grain.difspotA(1:nop)=grain.difspots(1:nop);
+grain.difspotB(1:nop)=grain.difspots(nop+1:end);
+grain.omegaA(1:nop)=grain.omega(1:nop);
+grain.omegaB(1:nop)=grain.omega(1:nop)+180;
+grain.etaA(1:nop)=grain.eta(1:nop);
+grain.etaB(1:nop)=180-grain.eta(1:nop);
+grain.thetaA(1:nop)=grain.theta(1:nop);
+grain.thetaB(1:nop)=grain.theta(1:nop);
+
+grain.hklB(1:nop,:)=-grain.hklA(1:nop,:);
+grain.plA(1:nop,:)=grain.pl(1:nop,:);
+grain.plB(1:nop,:)=-grain.plA(1:nop,:);
+% Have any of these spots already been matched, but not been assigned to the grain ?
+
+
+if size(grain.hkl,2)==3
+ grain.hkl(:,4)=grain.hkl(:,3);
+ grain.hkl(:,3)=NaN;
+end
+
+for j=1:length(difspotID)
+ [pid, gid, difaid, difbid, omega_p, eta_p, theta_p]=mym(sprintf('select pairid, grainid, difaid, difbid, omega, eta, theta from %s where difaid=%d or difbid=%d',acq.pair_tablename, difspotID(j), difspotID(j)));
+ if ~isempty(pid) && isnan(gid) %if there is a pair, we can update the grainid to this grain to stop it being considered again
+
+ if length(pid)>=2
+ disp('hm - found two or more entries for the same difspot in pairtable...will take the first one - database corrupted ?!?');
+ pid=pid(1);
+ gid=gid(1);
+ difaid=difaid(1);
+ difbid=difbid(1);
+ omega_p=omega_p(1);
+ eta_p=eta_p(1);
+ theta_p=theta_p(1);
+ end
+
+ if difaid==difspotID(j)
+ pairspotID=difbid;
+ else
+ pairspotID=difaid;
+ end
+ [tmp,index]=find(o.difspotID==pairspotID);
+ if isempty(index)
+ disp(sprintf('Matching thinks that spot %d has a pair, but ForwardSimulation can not find it, we will exclude this spot...',difspotID(j)));
+ o.flag(ind(j))=5; % 5 for problem
+ else
+ nop=nop+1;
+ o.flag(ind(j))=4;
+ o.flag(index)=4;
+ h=o.h(ind(j));
+ k=o.k(ind(j));
+ i=o.i(ind(j));
+ l=o.l(ind(j));
+
+ mysqlcmd=sprintf('select plX, plY, plZ from %s where difAID=%d',acq.pair_tablename, difaid);
+ [plX, plY, plZ]=mym(mysqlcmd);
+
+
+ grain.pairid=[grain.pairid,pid];
+ grain.difspotA=[grain.difspotA;difaid];
+ grain.difspotB=[grain.difspotB;difbid];
+ grain.omegaA=[grain.omegaA;omega_p];
+ grain.omegaB=[grain.omegaB;omega_p+180];
+ grain.etaA=[grain.etaA;eta_p];
+ grain.etaB=[grain.etaB;180-eta_p];
+ grain.thetaA=[grain.thetaA;theta_p];
+ grain.thetaB=[grain.thetaB;theta_p];
+ grain.hklA=[grain.hkl;[h k i l]];
+ grain.hklB=[-grain.hkl;-[h k i l]];
+ grain.plA=[grain.pl;[plX plY plZ]];
+ grain.plB=[-grain.pl;-[plX plY plZ]];
+
+ if app.update_db
+ %update grainID in spotpair table
+ mysqlcmd=dbUpdate(acq.pair_tablename,'pairid',pid,'GrainID',grain.id);
+ mym(mysqlcmd);
+ end
+ end
+ end
+end
+
+hkl=[];
+
+j=1;
+
+% Looking for pairs in the remaining list of unpaired spots
+ind=[];
+for i=1:length(o.flag)
+ if o.flag(i)==1
+ ind(j)=i; % ind contains line-ids of remaining flag "1" spots in output structure
+ omega(j)=o.omega(i);
+ eta(j)=o.eta(i);
+ theta(j)=o.theta(i);
+ hkl(j,1:4)=[o.h(i),o.k(i),o.i(i),o.l(i)];
+ difspotID(j)=o.difspotID(i);
+ j=j+1;
+ end
+end
+
+%% try to find the 180 degree offset (omega-pair) spot of the spot
+
+newpairs=[];
+j=1;
+
+for i=1:length(ind)
+
+
+ index1=[];
+ index2=[];
+ index=[];
+
+ [dummy,dummy2,index1]=intersect(-hkl(i,:),hkl,'rows'); % find the first instance of -h-k-l in the list of spots
+ if ~isempty(index1)
+ hkl(index1,:)=[0 0 0 0]; % a first -h-k-l has been found - set it temporarily to 0 0 0
+ [dummy,dummy2,index2]=intersect(-hkl(i,:),hkl,'rows'); % is there a second -h -k -l ?
+ hkl(index1,:)=-hkl(i,:); % set back the first one to its original values
+ if ~isempty(index2)
+ index=[index1,index2]; % these are the possible pair candidates - now investigate if valid
+ end
+
+ pair_omega=mod(omega(i)+180,360); % calculate the 180 degree offset omega value
+ delta_omega=abs(omega(index)-pair_omega);
+ goodone=find(delta_omega==min(delta_omega),1); % select the -h-k-l reflection fulfilling omega2=omega1+180
+
+ thr_max_offset=min(20,2*round(match.thr_max_offset/(abs(sind(eta(i))))))*180/acq.nproj;
+
+ if delta_omega(goodone)<thr_max_offset
+ newpairs(j,1)=i;
+ newpairs(j,2)=index(goodone);
+
+ j=j+1;
+ else
+ if isnan(o.i(1))
+ disp(sprintf('Could not find the omega pair spot of %d %d %d', hkl(i,1),hkl(i,2),hkl(i,4)));
+ else
+ disp(sprintf('Could not find the omega pair spot of %d %d %d %d', hkl(i,1),hkl(i,2),hkl(i,3),hkl(i,4)));
+ end
+ end %if delta_omega
+ end
+end
+
+newpairs=sort(newpairs,2);
+newpairs=unique(newpairs,'rows');
+
+if ~isempty(newpairs)
+ newhkl=hkl(newpairs(:,1),:);
+else
+ newhkl=[];
+end
+
+%% get eta, theta, omega from these pairs
+
+% Theoretical Bragg angles from energy and crystal structure
+% modified by LNervo
+%[tmp,results]=gtnew_calculate_twotheta() ;
+%comp_2thetas=results.centre; % coloumn vector of all valid 2theta angles in degrees
+info
+comp_2thetas=info.tt; % coloumn vector of all valid 2theta angles in degrees
+thr_ang=match.thr_ang;
+
+
+corr_rot_to_det=match.corr_rot_to_det; % correction of rot_to_det in pixels (0)
+rot_to_det=acq.dist/acq.pixelsize+corr_rot_to_det;
+tiltY=match.tiltY;
+tiltZ=match.tiltZ;
+sorZ = floor(acq.ydet/2) + 0.5; % between two central pixels (if even)
+
+
+
+
+% now fill in the information into the grain structure
+% first spot should correspond to spot in first half of scan
+grain.nof_newpairs=size(newpairs,1);
+
+
+for i=1:size(newpairs,1)
+ if omega(newpairs(i,1))>omega(newpairs(i,2))
+ newpairs(i,:)=newpairs(i,[2,1]);
+ newhkl(i,:)=-newhkl(i,:); % check this !!!!
+ end
+
+ nop=nop+1;
+ difAID=difspotID(newpairs(i,1));
+ difBID=difspotID(newpairs(i,2));
+ query=sprintf('select CentroidX,CentroidY,CentroidImage,BoundingBoxXsize,BoundingBoxYsize,Integral from %sdifspot where DifspotID=%d',name,difAID);
+ [sp1.centX,sp1.centY,sp1.cenIM,sp1.bbXsize,sp1.bbYsize,sp1.integral]=mym(query);
+ query=sprintf('select CentroidX,CentroidY,CentroidImage,BoundingBoxXsize,BoundingBoxYsize,Integral from %sdifspot where DifspotID=%d',name,difBID);
+ [sp2.centX,sp2.centY,sp2.cenIM,sp2.bbXsize,sp2.bbYsize,sp2.integral]=mym(query);
+
+ theta= gtThetaOfPairs(sp1.centX,sp1.centY,sp2.centX,sp2.centY,rot_to_det,tiltY,tiltZ,acq.rotx,sorZ) ; % check with Peter if tilts are taken correctly into account
+ % Check theta angle and plane family for each candidate
+ [theta_OK,angle_diff,thetatype] = sfCheck_thetas(sp1.centX,sp1.centY,sp2.centX,sp2.centY,...
+ rot_to_det,tiltY,tiltZ,acq.rotx,sorZ,comp_2thetas,2*thr_ang); % increased the angle tolerance to 2*thr_ang here
+
+
+ eta= gtEtaOfPairs(sp1.centX,sp1.centY,sp2.centX,sp2.centY,tiltY,tiltZ,acq.rotx,sorZ) ;
+ etaB=mod(180-eta,360); % eta of spotB
+ omega1=sp1.cenIM/acq.nproj*180; %nproj is number of images over 180 degrees
+ omega2=sp2.cenIM/acq.nproj*180;
+ omegamean=mean([omega1,omega2-180]);
+ omegaB=omegamean+180;
+ sp1.id=difAID;
+ sp2.id=difBID;
+ % Diffraction vector in LAB coordinates:
+ diff_vec_setup= gtDiffVecInLab(sp1.centX,sp1.centY,sp2.centX,sp2.centY,rot_to_det,tiltY,tiltZ,acq.rotx,sorZ) ;
+
+ % Plane normal in LAB coordinates:
+ pl_vec_setup= gtDiffPlaneNormInLab(sp1.centX,sp1.centY,sp2.centX,sp2.centY,rot_to_det,tiltY,tiltZ,acq.rotx,sorZ);
+
+ % Diffraction vector and plane normal in SAMPLE coordinates:
+ [diff_vec_sam(1),diff_vec_sam(2),diff_vec_sam(3)] = gtTrLabToSam(diff_vec_setup(1),diff_vec_setup(2),diff_vec_setup(3),omegamean) ;
+ [pl_vec_sam(1),pl_vec_sam(2),pl_vec_sam(3)] = gtTrLabToSam(pl_vec_setup(1),pl_vec_setup(2),pl_vec_setup(3),omegamean) ;
+
+ % Spot centers from SCINTILLATOR PLANE into LAB coordinates (in 3D):
+ [lab_cent1(1),lab_cent1(2),lab_cent1(3)] = gtTrScToLab(sp1.centX,sp1.centY,rot_to_det,tiltY,tiltZ,acq.rotx,sorZ) ;
+ [lab_cent2(1),lab_cent2(2),lab_cent2(3)] = gtTrScToLab(sp2.centX,sp2.centY,rot_to_det,tiltY,tiltZ,acq.rotx,sorZ) ;
+
+ % Spot centers from LAB into SAMPLE coordinates:
+ [samcentA(1),samcentA(2),samcentA(3)] = gtTrLabToSam(lab_cent1(1),lab_cent1(2),lab_cent1(3),omega1) ;
+ [samcentB(1),samcentB(2),samcentB(3)] = gtTrLabToSam(lab_cent2(1),lab_cent2(2),lab_cent2(3),omega2) ;
+
+ % Average intensity and bounding box sizes from the two spots:
+ av_intint=(sp1.integral+sp2.integral)/2 ;
+ av_bbXsize=(sp1.bbXsize+sp2.bbXsize)/2 ;
+ av_bbYsize=(sp1.bbYsize+sp2.bbYsize)/2 ;
+
+ if app.update_db
+ % Variables in table
+ % samcent: difspot centroid in the sample system
+ % lines connecing the two spots in the sample system
+ % lor: line origin (centroid of spot A);
+ % ldir: unit vector of line direction
+ % pl: plane normal in sample system
+
+ mysqlcmd=sprintf(['INSERT INTO %s (difAID, difBID, theta, eta, omega, etaB, omegaB, '...
+ 'avintint, avbbXsize, avbbYsize, '...
+ 'samcentXA, samcentYA, samcentZA, samcentXB, samcentYB, samcentZB, '...
+ 'ldirX, ldirY, ldirZ, plX, plY, plZ, thetatype) ',...
+ 'VALUES (%d,%d,%0.20g,%0.20g,%0.20g,%0.20g,%0.20g,%0.20g,%0.20g,%0.20g,'...
+ '%0.20g,%0.20g,%0.20g,%0.20g,%0.20g,%0.20g,%0.20g,%0.20g,%0.20g,%0.20g,%0.20g,%0.20g,%d)'],...
+ acq.pair_tablename,sp1.id,sp2.id,theta,eta,omegamean,etaB,omegaB,av_intint,av_bbXsize,av_bbYsize,...
+ samcentA(1),samcentA(2),samcentA(3),samcentB(1),samcentB(2),samcentB(3),...
+ diff_vec_sam(1),diff_vec_sam(2),diff_vec_sam(3),pl_vec_sam(1),pl_vec_sam(2),pl_vec_sam(3),thetatype) ;
+ mym(mysqlcmd);
+ pairid=mym(sprintf('select pairid from %s where difaid=%d',acq.pair_tablename,sp1.id));
+ grain.pairid=[grain.pairid,pairid(1)];
+ end
+
+
+ grain.difspotA(nop)=difspotID(newpairs(i,1));
+ grain.difspotB(nop)=difspotID(newpairs(i,2));
+ grain.omegaA(nop)=omegamean;
+ grain.omegaB(nop)=omegamean+180;
+ grain.etaA(nop)=eta;
+ grain.etaB(nop)=etaB;
+ grain.thetaA(nop)=theta;
+ grain.thetaB(nop)=theta;
+ grain.plA(nop,:)=pl_vec_setup;
+ grain.plB(nop,:)=-pl_vec_setup;
+ grain.hklA(nop,:)=newhkl(i,:);
+ grain.hklB(nop,:)=-newhkl(i,:);
+
+ % and fill in the output structure
+ ind1=ind(newpairs(i,1));
+ ind2=ind(newpairs(i,2));
+ o.omega(ind1)=omegamean;
+ o.omega(ind2)=omegamean+180;
+ o.theta(ind1)=theta;
+ o.theta(ind2)=theta;
+ o.eta(ind1)=eta;
+ o.eta(ind2)=180-eta;
+ o.flag(ind1)=4;
+ o.flag(ind2)=4;
+end
+
+
+%%
+% try to calculate in the missing information for the remaining, unpaired
+% spots based on the measured spot positions and known grain position
+
+
+ind=[];
+j=1;
+for i=1:length(o.image)
+ if o.flag(i)==1
+ ind(j)=i; % ind now contains the remaining unpaired spots
+ j=j+1;
+ end
+end
+
+grain.nof_singles=length(ind);
+
+% allocate vectors containing single spots
+grain.difspotC=zeros(length(ind),1);
+grain.omegaC=zeros(length(ind),1);
+grain.etaC=zeros(length(ind),1);
+grain.thetaC=zeros(length(ind),1);
+grain.hklC=zeros(length(ind),4);
+grain.plC=zeros(length(ind),3);
+
+for j=1:length(ind)
+
+ [grainx,grainy,grainz]=gtTrSamToLab(grain.center(1),grain.center(2),grain.center(3),o.omega(ind(j)));
+
+ % calculate the distances in pixels with respect to the laboratory system
+ dx= acq.dist/acq.pixelsize-grainx;
+ dy= o.udiffr(ind(j))-acq.rotx-grainy;
+ dz= acq.ydet/2-o.vdiffr(ind(j))-grainz;
+
+ grain.omegaC(j) = o.omega(ind(j));
+ grain.etaC(j)= atan2(dy,dz)/pi*180;
+ grain.thetaC(j) = atand(sqrt(dy.^2+dz.^2)/dx)/2;
+ grain.difspotC(j)=o.difspotID(ind(j));
+
+ kh=[dx dy dz];
+ kh=kh/norm(kh); % normalized difracted beam direction
+ k0=[1 0 0]; % normalized incoming beam direction
+ G= (kh-k0)/norm(kh-k0); % normalized Laboratory system scattering vector
+ %refl = find(abs(tt.centre-Theta(ndx+j)*2)==min(abs(tt.centre-Theta(ndx+j)*2))); % find out what type of reflection this Theta corresponds to
+ grain.hklC(j,:) = [o.h(ind(j)), o.k(ind(j)), o.i(ind(j)), o.l(ind(j))];
+ [pn(1),pn(2),pn(3)]= gtTrLabToSam(G(1),G(2),G(3),o.omega(ind(j))); % normalized plane_normal in sample coordinate system
+ grain.plC(j,:)=[pn(1),pn(2),pn(3)];
+
+end
+
+
+if any(isnan(o.i))
+ grain.hklA=grain.hklA(:,[1 2 4]);
+ grain.hklB=grain.hklB(:,[1 2 4]);
+ grain.hklC=grain.hklC(:,[1 2 4]);
+end
+
+disp(sprintf('\n\n Forward Simulation has detected %d new spotpairs and %d single spots\n\n',grain.nof_newpairs,grain.nof_singles));
+
+
+%%
+ function [theta_OK,angle_diff,thetatype]=sfCheck_thetas(scXA,scYA,scXB,scYB,rottodet,tiltY,tiltZ,rotx,sorZ,comp_2thetas,thr_ang)
+ % Checks whether the position of two spots gives a diffraction vector which
+ % corresponds to a valid 2theta angle within the given threshold.
+
+ angle=2*gtThetaOfPairs(scXA,scYA,scXB,scYB,rottodet,tiltY,tiltZ,rotx,sorZ);
+
+ diff=angle-comp_2thetas;
+
+ if all(abs(diff) > thr_ang) % checks if all values are non-zero
+ theta_OK=false;
+ else
+ theta_OK=true ;
+ end
+ [mini,thetatype]=min(abs(diff));
+ angle_diff=diff(thetatype);
+ end
+
+end
\ No newline at end of file
diff --git a/zUtil_Optimization/gtFullOptimizer3.m b/zUtil_Optimization/gtFullOptimizer3.m
new file mode 100644
index 0000000000000000000000000000000000000000..7b3a49eaf1c9d256658e7c1ed33cc98579f53a45
--- /dev/null
+++ b/zUtil_Optimization/gtFullOptimizer3.m
@@ -0,0 +1,72 @@
+function gtFullOptimizer3(grainid, parameters, difspots)
+% gtFullOptimizer3(grainid, parameters, difspots)
+% grainid = grain{grainid}
+% difspots = list of spots id
+
+% save some info in the opt structure
+opt.id=grainid.id;
+opt.workdir=pwd;
+
+if ~isfield(grainid,'difspots')
+ grainid.difspots=difspots;
+end
+
+ind=strfind(pwd,'/');
+opt.dataset=opt.workdir(ind(end)+1:end);
+
+opt.grain='grain.mat';
+opt.parfile='parameters.mat';
+opt.gcenter=grainid.center;
+opt.gR_vector=grainid.R_vector;
+
+% Create the full summed difspots image
+if isfield(grainid,'difspots')
+ full=gtPlaceDifspotinFull_ln(grainid.difspots,parameters,'save',1);
+else
+ full=gtPlaceDifspotinFull_ln(difspots,parameters);
+end
+
+% generate grain fed
+name=sprintf('gr%d_fed.mat',grainid.id);
+if ~exist(name,'file')
+ gr=gtFedGenerateRandomGrain2(grainid,full,parameters,'save',1);
+else
+ check=inputwdefault('Do you want to load the existing fed grain?','y');
+ if strcmpi(check,'y')
+ gr=load(name);
+ gr=gr.grain;
+ else
+ gr=gtFedGenerateRandomGrain2(grainid,full,parameters);
+ end
+end
+
+% load original geometry
+[geometry,pippo]=gtOptimizeDetector2(gr,parameters); % was gtOptimizeDetector
+opt.fed.geo=geometry;
+opt.fed.uv=gr.allblobs.uv;
+% removes NaN
+ind=find(isnan(opt.fed.uv(:,1))); % size -> size/2
+opt.ind.nan=ind;
+opt.fed.uv(ind,:)=[];
+opt.db.uv=[pippo.bbox.xcen pippo.bbox.ycen];
+% removes zeros from database
+ind2=find(opt.db.uv(:,1)==0);
+ind3=find(opt.db.uv(:,1)~=0);
+opt.db.uv(ind2,:)=[];
+opt.ind.zeros=ind2;
+opt.ind.nonzeros=ind3;
+
+% fit the spots positions taken from the database with the current geometry as guess
+disp('*** Fit ***')
+fitgeo=gtTuneDetector3(opt.fed.geo,gr,opt.db.uv,ind3);
+opt.fit.geo=fitgeo;
+opt.fit.uv=fitgeo.Fnew(:,1:2);
+
+
+name=sprintf('opt%d.mat',grainid.id);
+save(name,'opt');
+
+tmp=gtPlotUVW(full,opt.fit.uv,'r',opt.fed.geo,[0 200],'After optimization : ');
+
+
+end % end function
diff --git a/zUtil_Optimization/gtGetGeometry.m b/zUtil_Optimization/gtGetGeometry.m
new file mode 100644
index 0000000000000000000000000000000000000000..fc23452db4e2e6ddefde37844b8dec455b4ee9f4
--- /dev/null
+++ b/zUtil_Optimization/gtGetGeometry.m
@@ -0,0 +1,54 @@
+function geo = gtGetGeometry(parameters, fed)
+%% FUNCTION GTGETGEOMETRY Get geometry from parameters
+% geo = gtGetGeometry(parameters, fed)
+%
+% fed = flag for fed analysis <logical 1x1> {0}
+%
+% Version 002 21-09-2011 LNervo
+
+if ~exist('fed','var') || isempty(fed)
+ fed=false;
+end
+if fed
+geometry{1}='detpos0';
+geometry{2}='detdiru0';
+geometry{3}='detdirv0';
+geometry{4}='beamdir0';
+geometry{5}='rotdir0';
+geometry{6}='detucen0';
+geometry{7}='detvcen0';
+
+if ~fed
+ geometry{8}='detdir0';
+ geometry{9}='pixelsize';
+ geometry{10}='nproj';
+ geometry{11}='rotation_axis';
+end
+
+if ~exist('parameters','var') || isempty(parameters)
+ load parameters
+end
+
+if isfield(parameters,'acq')
+ acq=parameters.acq;
+else
+ acq=parameters;%25/08/2011
+end
+
+if ~all(isfield(acq,geometry))
+ disp('Some parameters are missing. Running gtCheckParameters...')
+ gtCheckParameters
+ return
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% save the current geometry
+for i=1:length(geometry)
+ var=char(geometry(i));
+ if fed
+ var=var(1:end-1);
+ end
+ geo.(var)=acq.(char(geometry(i)));
+end
+
+end
\ No newline at end of file
diff --git a/zUtil_Optimization/gtLimitsTwoTheta.m b/zUtil_Optimization/gtLimitsTwoTheta.m
new file mode 100644
index 0000000000000000000000000000000000000000..9930db66cb033d58b66b05427e99a491adf5c157
--- /dev/null
+++ b/zUtil_Optimization/gtLimitsTwoTheta.m
@@ -0,0 +1,147 @@
+function [minangle,maxangle]=gtLimitsTwoTheta(updatepar,calc)
+%% GTLIMITSTWOTHETA Calculate TwoTheta maximum and minimum
+% [minangle,maxangle]=gtLimitsTwoTheta(updatepar)
+%
+% Load parameters.mat and save (if updatepar=1) new fields in parameters.acq:
+%
+% detanglemin = minimum twotheta seen from the detector
+% detanglemax = maximum twotheta seen from the detector
+%
+% maxradius = maximum radius from the center of the detector
+% detdir0 = direction normal to the detector (been the detector)
+%
+% Possible detector configuration:
+% 1) inline detector and centered w.r.t. beam direction - detpos0(x,0,0)
+% 2) inline detector with z-offset - detpos0(x,0,z)
+% 3) vertical detector centered w.r.t. the sample - detpos0(x,0,z)
+%
+% FUNCTIONS USED:
+%
+% gtConvEnergyToWavelength(energy) : lambda
+% -> energy <double 1x1> [keV]
+% lambda <double 1x1> [Angstrom (Ã…)]
+%
+
+if exist('parameters.mat','file')
+ load parameters;
+else
+ error('Cannot find the parameters.mat file!')
+end
+
+if ~exist('updatepar','var') || isempty(updatepar)
+ updatepar=0;
+end
+if ~exist('calc','var') || isempty(calc)
+ calc=1;
+end
+
+acq=parameters.acq;
+spacegroup=acq.spacegroup;
+energy=acq.energy;
+lambda=gtConvEnergyToWavelength(energy);
+latticepar=acq.latticepar;
+detpos=acq.detpos0
+
+% Calculate minimum and maximum twotheta angle of the experiment
+
+upd=calc;
+if isfield(parameters.acq,'detanglemin') && isfield(parameters.acq,'detanglemax') && ...
+ ~isempty(parameters.acq.detanglemin) && ~isempty(parameters.acq.detanglemax) ...
+ && parameters.acq.detanglemax~=0
+
+ minangle=parameters.acq.detanglemin;
+ maxangle=parameters.acq.detanglemax;
+ upd=false;
+else
+ upd=true;
+end
+
+
+if upd
+ if ~isfield(parameters.acq,'maxradius')
+ maxradius=sqrt((acq.ydet/2)^2+(acq.xdet/2)^2);% in pixel
+ parameters.acq.maxradius=maxradius;
+ if updatepar
+ save parameters parameters
+ end
+ else
+ maxradius=acq.maxradius;
+ end
+
+ % dist in mm
+ % maxradius in pixel
+ % detpos0 in mm
+
+ if isfield(parameters.acq,'detdir0')
+ detdir=parameters.acq.detdir0;
+ else
+ detdir=cross(parameters.acq.detdiru0,parameters.acq.detdirv0);
+ parameters.acq.detdir0=detdir;
+ if updatepar
+ save parameters parameters
+ end
+ end
+
+ x=[1 0 0];
+ y=[0 1 0];
+ z=[0 0 1];
+ if abs(dot(detdir,x))==1 % detector inline normal // axis X
+ % suppose detpos(2)=0
+ % detector must be perpendicular to the plane x-z
+ % case inline detector and centered w.r.t. beam direction
+ disp('Inline detector')
+ R=maxradius*acq.pixelsize;
+ if detpos(3)==0
+ minangle=0;
+ else
+ % case inline detector with z-dir offset
+ minangle=atand( ((detpos(3) - R)/detpos(1)) );
+ end
+ maxangle=atand( ((detpos(3) + R)/detpos(1)) );
+
+ elseif abs(dot(detdir,z))==1 % detector vertical normal // axis Z
+ disp('Vertical detector')
+ % case vertical detector offset by xd w.r.t. z-axis (LAB ref.)
+ % detector must be perpendicular to the plane x-z
+ % suppose detpos(2)=0
+ % D distance sample center - detector center
+ %D=sqrt(detpos(1)+detpos(3)); % ?????????
+ R=maxradius*acq.pixelsize;
+ maxangle=atand( abs(detpos(3))/(detpos(1)-R) );
+ minangle=atand( abs(detpos(3))/(detpos(1)+R) );
+ if minangle < 0
+ minangle = abs(minangle)+90;
+ end
+ if maxangle < 0
+ maxangle = abs(maxangle)+90;
+ end
+
+ elseif abs(dot(detdir,y))==1 % detector vertical but rotated normal // axis Y
+ disp('Vertical detector - rotated 90 deg')
+ % case vertical detector
+ % suppose detpos(3)=0
+ % D distance sample center - detector center
+ R=maxradius*acq.pixelsize;
+ %D=sqrt(detpos(1)*detpos(1)+detpos(2)*detpos(2));
+ maxangle=atand( abs(detpos(2))/(detpos(1)-R) );
+ minangle=atand( abs(detpos(2))/(detpos(1)+R) );
+ if maxangle < 0
+ maxangle = abs(maxangle)+90;
+ end
+ if minangle < 0
+ minangle = abs(minangle)+90;
+ end
+ else
+ error('Check if detector direction is unusual. Normally the detector is inline or vertical.')
+ end
+
+ parameters.acq.detanglemin=minangle;
+ parameters.acq.detanglemax=maxangle;
+ if updatepar
+ save parameters parameters
+ end
+end
+
+end
+
+%%%
\ No newline at end of file
diff --git a/zUtil_Optimization/gtOptimizeDetector2.m b/zUtil_Optimization/gtOptimizeDetector2.m
new file mode 100644
index 0000000000000000000000000000000000000000..9c1657b69e4e5d25a41dfe791d2ecdbb10956e2d
--- /dev/null
+++ b/zUtil_Optimization/gtOptimizeDetector2.m
@@ -0,0 +1,333 @@
+function [geo,newgr]=gtOptimizeDetector2(grainFed,parameters)
+%% newgr=gtOptimizeDetector2(grainFed,parameters)
+%
+% Get geometry from parameters, get difspots positions and bbox from
+% database, predict positions from fed code and
+%
+% INPUT:
+% grainFed = output from gtFedGenerateRandomGrain2
+% parameters
+%
+% OUTPUT:
+% geo = geometry from the parameters file
+% newgr = structure with geometry (current geometry),
+% bbox (values read from database)
+% fed (values predicted with fed code)
+% table (comparison between fed and database
+% values)
+
+geo=gtGetGeometry(parameters);
+
+%% Current geometry
+
+disp('Getting the detector geometry...')
+
+detdiru = geo.detdiru0';%
+detdirv = geo.detdirv0';%
+Qdet = gtFedDetectorProjectionTensor(detdiru,detdirv,1,1);%
+tn = cross(detdiru,detdirv);%
+detpos = (geo.detpos0./geo.pixelsize)'; % pixels
+omstep = 180/geo.nproj;
+uvorig = [geo.detucen0, geo.detvcen0]';
+
+disp('...done!')
+
+%% Getting info from grain
+
+if isempty(grainFed)
+ disp('Calculating output from gtFedGenerateRandomGrain2...')
+ grain=load('grain.mat');
+ full=gtPlaceDifspotinFull_ln(grain{grainFed.id}.difspots,parameters);
+ gr=gtFedGenerateRandomGrain2(grain{grainFed.id},full);
+ grainFed=gr;
+end
+
+grainID=grainFed.id;
+csam = grainFed.center';
+R_vector=grainFed.R_vector;
+difspots=grainFed.difspots;
+
+all=grainFed.allblobs;
+fedtable=[];
+
+
+% NaN rows
+indn=find(isnan(all.uvw(:,1)));
+
+hklsp=all.hklsp;hklsp(indn,:)=[];
+srot=all.srot;
+for i=1:size(indn)
+ srot(indn*3-2:indn*3,:)=[];
+end
+gfomega=all.omega;gfomega(indn,:)=[];
+dvec=all.dvec;dvec(indn,:)=[];
+theta=all.theta;theta(indn,:)=[];
+eta=all.eta;eta(indn,:)=[];
+uvwfed=all.uvw;uvwfed(indn,:)=[];
+uvfed=all.uv;uvfed(indn,:)=[];
+images=round((gfomega/180)*geo.nproj);
+
+
+fed.grainID=grainID;
+fed.center=csam;
+fed.R_vector=R_vector;
+fed.hklsp=hklsp;
+fed.omega=gfomega;
+fed.theta=theta;
+fed.eta=eta;
+fed.uvwfed=uvwfed;
+fed.uvfed=uvfed;
+fed.difspots=difspots';
+
+%% Get info from the database
+
+match=parameters.match;
+name=parameters.acq.name;
+
+%get the mean difspot, X and Y sizes for the current spots in the grain
+dif_Xsize=grainFed.stat.bbxsmean;
+dif_Xsize_ratio=grainFed.stat.bbxsrat*1.3;
+dif_Ysize=grainFed.stat.bbysmean;
+dif_Ysize_ratio=grainFed.stat.bbysrat*1.3;
+
+dif_Xsearch=2*grainFed.stat.bbxsstd;
+dif_Ysearch=4*grainFed.stat.bbysstd;
+
+if ~isfield(grainFed.stat,'areamean')
+ % we need to add this value in gtINDEXTER
+ grainFed.stat.areamean=0.7*grainFed.stat.bbxsmean*grainFed.stat.bbysmean; % 0.7 is a rough estimate of the spot area / BB-area ratio
+end
+
+dif_int=grainFed.stat.intmean;
+dif_area=grainFed.stat.areamean;
+dif_area_ratio=dif_Xsize_ratio*dif_Ysize_ratio;
+
+
+
+
+gtDBConnect
+
+grain_id=[];
+expected_id=[];
+found_id=[];
+image=[];
+omega=[];
+x0=[];
+y0=[];
+xsize=[];
+ysize=[];
+xcen=[];
+ycen=[];
+bbxcen=[];
+bbycen=[];
+area=[];
+
+
+for i=1:length(gfomega)
+
+ thr_max_offset=min(20,round(match.thr_max_offset/(abs(sind(eta(i)))))); % allow for a 1/sin(eta) dependency on omega search range - limited to a max offset of 20 images
+ mysqlcmd=sprintf(['select grainID, difspotID, CentroidX, CentroidY, '...
+ 'BoundingBoxXorigin, BoundingBoxYorigin, BoundingBoxXsize, '...
+ 'BoundingBoxYsize, CentroidImage, Area from %sdifspot where '...
+ '(%d between CentroidImage-%d and CentroidImage+%d) '...
+ 'and (abs(%sdifspot.CentroidX-%d)/%d)<0.5 '...
+ 'and (abs(%sdifspot.CentroidY-%d)/%d)<0.5 '...
+ 'order by ((abs(BoundingBoxXsize-%d)/%d)+(abs(BoundingBoxYsize-%d)/%d)) limit 1'],...
+ name,...
+ round(gfomega(i)/omstep), thr_max_offset, thr_max_offset,...
+ name, uvfed(i,1), dif_Xsize,...
+ name, uvfed(i,2), dif_Ysize,...
+ dif_Xsize, dif_Xsize, dif_Ysize, dif_Ysize);
+
+
+ [grainID, expected, CentroidX, CentroidY, ...
+ BoundingBoxXorigin, BoundingBoxYorigin, ...
+ BoundingBoxXsize, BoundingBoxYsize, ...
+ CentroidImage, Area] = mym(mysqlcmd);
+
+ if isempty(expected)
+ grain_id=[grain_id; 0];
+ expected_id=[expected_id; 0];
+ image=[image; 0];
+ omega=[omega; 0];
+ x0=[x0; 0];
+ y0=[y0; 0];
+ xsize=[xsize; 0];
+ ysize=[ysize; 0];
+ xcen=[xcen; 0];
+ ycen=[ycen; 0];
+ bbxcen=[bbxcen; 0];
+ bbycen=[bbycen; 0];
+ area=[area; 0];
+ found_id=[found_id; 0];
+ else
+ grain_id=[grain_id; grainID];
+ expected_id=[expected_id; expected];
+ fprintf('expected_id: %d \n',expected);
+ image=[image; CentroidImage];
+ omega=[omega; CentroidImage*omstep];
+ x0=[x0; BoundingBoxXorigin];
+ y0=[y0; BoundingBoxYorigin];
+ xsize=[xsize; BoundingBoxXsize];
+ ysize=[ysize; BoundingBoxYsize];
+ xcen=[xcen; CentroidX];
+ ycen=[ycen; CentroidY];
+ bbxcen=[bbxcen; BoundingBoxXorigin+BoundingBoxXsize/2];
+ bbycen=[bbycen; BoundingBoxYorigin+BoundingBoxYsize/2];
+ area=[area; Area];
+
+ query=sprintf(['select difspotID from %sdifspot where '...
+ '(%d between CentroidImage-%d and CentroidImage+%d) '...
+ 'and (abs(%sdifspot.CentroidX-%d)/%d)<%f '...
+ 'and (abs(%sdifspot.CentroidY-%d)/%d)<%f '...
+ 'and (BoundingBoxXsize between %f and %f) and (BoundingBoxYsize between %f and %f) '...
+ 'and (Area between %f and %f) '...
+ 'order by ((abs(BoundingBoxXsize-%d)/%d)+(abs(BoundingBoxYsize-%d)/%d)) limit 1'],... % this could be improved according to Peters goodness criteria...
+ name,...
+ round(gfomega(i)/omstep), thr_max_offset, thr_max_offset,...
+ name, uvfed(i,1), dif_Xsize, dif_Xsearch/dif_Xsize,...
+ name, uvfed(i,2), dif_Ysize, dif_Ysearch/dif_Ysize,... % use the same (more relaxed) X condition for the Y COM position
+ dif_Xsize/dif_Xsize_ratio, dif_Xsize*dif_Xsize_ratio, dif_Ysize/dif_Ysize_ratio, dif_Ysize*dif_Ysize_ratio,...
+ dif_area/dif_area_ratio, dif_area*dif_area_ratio,...
+ dif_Xsize, dif_Xsize, dif_Ysize, dif_Ysize);
+
+ [found] = mym(query);
+
+ if ~isempty(found)
+ found_id=[found_id; found];
+ fprintf('found_id: %d \n',found);
+ else
+ found_id=[found_id; 0];
+ end
+
+ end % end expected
+
+end %loop on omega
+
+% save info into parameters
+
+bbox.grain_id=grain_id;
+bbox.expected_id=expected_id;
+bbox.found_id=found_id;
+bbox.x0=x0;
+bbox.y0=y0;
+bbox.xsize=xsize;
+bbox.ysize=ysize;
+bbox.xcen=xcen;
+bbox.ycen=ycen;
+bbox.bbxcen=bbxcen;
+bbox.bbycen=bbycen;
+bbox.image=image;
+bbox.omega=omega;
+bbox.area=area;
+
+fed.omstep=omstep;
+
+%% Find the differences between found and expected
+
+diff=findDifspots(bbox.expected_id,bbox.found_id,name);
+diff.query='expected_id and found_id';
+bbox.diff=diff;
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% Compute uv coordinates and w(omega) values
+
+% fed.UV = [];
+% fed.UVW = [];
+%
+% n=size(gfomega,1)/4;
+% for j=1:n % must be from 0 to n-1
+%
+% % from grainFed=gtFedGenerateRandomGRain2(grain{grainid},full)
+% S=all.srot(12*j-11:12*j,:);
+% D=all.dvec(4*j-3:4*j,:);
+% O=all.omega(4*j-3:4*j); % omega values ordered like small, small+180, big, big+180 for the 4
+%
+% uv1a=gtFedPredictUVW(S(1:3,:),D(1,:)',csam,detpos,tn,Qdet,uvorig,O(1),omstep)
+% uv1b=gtFedPredictUVW(S(4:6,:),D(2,:)',csam,detpos,tn,Qdet,uvorig,O(2),omstep)
+% uv2a=gtFedPredictUVW(S(7:9,:),D(3,:)',csam,detpos,tn,Qdet,uvorig,O(3),omstep)
+% uv2b=gtFedPredictUVW(S(10:12,:),D(4,:)',csam,detpos,tn,Qdet,uvorig,O(4),omstep)
+%
+% fed.UV = [fed.UV;uv1a(1),uv1a(2);uv1b(1),uv1b(2);uv2a(1),uv2a(2);uv2b(1),uv2b(2)];
+% fed.UVW = [fed.UVW;uv1a(1),uv1a(2),uv1a(3);uv1b(1),uv1b(2),uv1b(3);uv2a(1),uv2a(2),uv2a(3);uv2b(1),uv2b(2),uv2b(3)];
+% end % size /4
+
+newgr.fed=fed;
+newgr.geometry=geo;
+
+
+%% Try to match the difspots and the hkl values
+
+if ~isempty(bbox.found_id)
+ diff2=findDifspots(bbox.found_id,fed.difspots,name);
+diff2.query='found_id and difspots';
+bbox.diff2=diff2;
+end
+
+newgr.bbox=bbox;
+
+size(fed.omega)
+size(bbox.omega)
+size(fed.uvfed)
+size(bbox.xcen)
+size(bbox.ycen)
+size(bbox.expected_id)
+size(bbox.found_id)
+size(fed.hklsp)
+
+if size(fed.omega,1)~=size(bbox.omega,1)
+ return
+end
+% table to compare fed results with database
+table=int32([fed.omega bbox.omega ...
+ fed.uvfed(:,1) bbox.xcen ...
+ fed.uvfed(:,2) bbox.ycen ...
+ bbox.expected_id bbox.found_id ...
+ fed.hklsp]);
+newgr.table=table;
+nonzero=find(table(:,2)~=0);
+newgr.ind=nonzero;
+
+%newgr.ind=find(bbox.expected_id~=0);
+
+end % end function
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+
+function list=findDifspots(list1,list2,name)
+
+
+[rows,ia,ib,neq,dif]=findDifferentRowIntoMatrix(list1,list2);
+
+dif(find(dif(:,1)==0),:)=[];
+list.ind=dif(:,2)';
+list.ID=dif(:,1)';
+
+for j=1:size(dif(:,1),1)
+ mysqlcmd=sprintf(['select grainID, CentroidX, CentroidY, CentroidImage, Integral, '...
+ 'BoundingBoxXorigin, BoundingBoxYorigin, BoundingBoxXsize, BoundingBoxYsize, Area from %sdifspot where difspotID=%d'],...
+ name, dif(j,1));
+
+ [grainID, CentroidX, CentroidY, CentroidImage, Integral, BoundingBoxXorigin, BoundingBoxYorigin, BoundingBoxXsize, BoundingBoxYsize, Area]=mym(mysqlcmd);
+
+ if ~isempty(BoundingBoxXorigin)
+ list.x0(j)=BoundingBoxXorigin;
+ list.y0(j)=BoundingBoxYorigin;
+ list.xsize(j)=BoundingBoxXsize;
+ list.ysize(j)=BoundingBoxYsize;
+ list.xcen(j)=CentroidX;
+ list.ycen(j)=CentroidY;
+ list.image(j)=CentroidImage;
+ list.int(j)=Integral;
+ list.area(j)=Area;
+ list.grainid(j)=grainID;
+ end
+end
+
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
diff --git a/zUtil_Optimization/gtOptimizeVertical.m b/zUtil_Optimization/gtOptimizeVertical.m
new file mode 100644
index 0000000000000000000000000000000000000000..37ffe174729f3045fd8aeec64db5e46885d7ce75
--- /dev/null
+++ b/zUtil_Optimization/gtOptimizeVertical.m
@@ -0,0 +1,79 @@
+function [newgeo,bbox,full_all,gr]=gtOptimizeVertical(grain,parameters,bbxsizemean,areamin,areamax,color,greys,save)
+% [newgeo,bbox,full_all,gr]=gtOptimizeVertical(grain,parameters,bbxsizemean,areamin,areamax,[color],[greys],[save])
+%
+% INPUT:
+% grain = grain{grainid}
+% parameters = parameters file
+% bbxsizemean = discrepancy between fed u,v and spot centroid u,v
+% | difspot.CentroidX - u | < 0.5 * bbxsizemean
+% | difspot.CentroidY - v | < 0.5 * bbxsizemean
+% areamin/areamax = Area < areamax & Area > areamin
+% color = color of predicted u,v points
+% greys = grey scale
+% save = true if saving
+%
+% OUTPUT:
+% newgeo = optimized geometry
+% bbox = structure with information from the database and from fed
+% full_all = image with spots
+% gr = grain updated structure
+
+if nargin < 5
+ disp('Usage: gtOptimizeVertical(grain,parameters,bbxsizemean,areamin,areamax,[color],[greys],[save])')
+ return
+end
+if isempty(color)
+ color='b';
+end
+if isempty(greys)
+ greys=[];
+end
+if ~exist('save','var') || isempty(save)
+ save=false;
+end
+
+
+gr=gtFedGenerateRandomGrain2(grain,zeros(parameters.acq.xdet),parameters,'save',save);
+uv=gr.allblobs.uv;
+omega=gr.allblobs.omega;
+
+% find spots' indexes intersecting the detector
+[uvdet,inddet]=findOnDet(uv);
+omegadet=omega(inddet);
+length(uv)
+
+% search in the database
+bbox=findSpotsWoIndexter(parameters,uvdet,omegadet,bbxsizemean,areamin,areamax);
+length(bbox.feduv)
+
+% place spots in the image (figure not shown)
+full_all=gtPlaceDifspotinFull_ln(bbox.expected_id,parameters,'save',save,'connected',true);
+
+% plot the spots using the current geometry
+tmp=gtPlotUVW(full_all,bbox.feduv,'b',parameters.acq,greys,sprintf('Before optimization : |x-u| < %2.0f : area [%d %d]',bbxsizemean*0.5,areamin,areamax),false);
+
+% get the indexes of spots found in the database from the full uv list
+[r,ia,ib,d]=findRowIntoMatrix(uv,bbox.feduv);
+[r2,ia2,ib2,d2]=findRowIntoMatrix(uv(ia,:),bbox.feduv);
+ind=sortrows([ib2 ia]);
+indsorted=ind(:,2);
+bbox.indsorted=indsorted;
+
+% optimize detector
+newgeo=gtTuneDetector3(parameters.acq,gr,bbox.dbuv,indsorted);
+
+disp('detdiru0: ')
+disp(newgeo.detdiru0)
+disp('detdirv0: ')
+disp(newgeo.detdirv0)
+
+% plot the spots positions with the optimized geometry
+tmp2=gtPlotUVW(full_all,newgeo.Fnew(:,1:2),color,newgeo,greys,sprintf('After optimization : blue [bbox.dbuv] %s [bbox.feduv]',color),true);
+hold on
+for i=1:length(bbox.dbuv)
+ plot(bbox.dbuv(i,1),bbox.dbuv(i,2),'g.','MarkerSize',20)
+end
+
+
+end
+
diff --git a/zUtil_Optimization/gtPlotUVW.m b/zUtil_Optimization/gtPlotUVW.m
new file mode 100644
index 0000000000000000000000000000000000000000..7f12b1d64ede1e113e8e6c8c56e636a74b5a633d
--- /dev/null
+++ b/zUtil_Optimization/gtPlotUVW.m
@@ -0,0 +1,141 @@
+function noiseuv=gtPlotUVW(im,uv,color,geo,scale,titlename,debug,noise,numnoise)
+% noiseuv=gtPlotUVW(im,uv,color,geo,titlename,debug,noise,numnoise)
+%
+% INPUT:
+% im = output from gtPlaceDifspotsinfull <double 2048x2048>
+% uv = u,v list to plot <double Nx2>
+% color = color to be used in the plot <char 1x1>
+% geo = output from gtTuneDetector3
+% scale = gray scaled from -100 to 400
+% titlename = title name <string 1x1>
+% debug = print comments or info <logical 1x1> {0}
+% noise = number of pixels noise <int8 1x1 {2}>
+% numnoise = percentage of noised spots <double 1x1 {}>
+%
+% OUTPUT:
+% noiseuv = list u,v noised
+
+figure
+if ~isempty(im)
+ if ~isempty(scale)
+ imshow(im,scale)
+ else
+ imshow(im,[])
+ end
+ hold on
+end
+
+if ~exist('titlename','var') || isempty(titlename)
+ titlename='';
+end
+
+if ~exist('debug','var') || isempty(debug)
+ debug=false;
+end
+if ~exist('noise','var') || isempty(noise)
+ noise=2;
+end
+if ~exist('numnoise','var') || isempty(numnoise)
+ numnoise=[];
+end
+
+N=size(uv,1);
+
+detusize=2048;
+detvsize=2048;
+uvorig=[1024.5;1024.5];
+beamdir=geo.beamdir0';
+detdiru=geo.detdiru0';
+detdirv=geo.detdirv0';
+detdir=cross(detdiru,detdirv);
+detpos=geo.detpos0; % mm
+
+if isfield(geo,'Xzero')
+ olddiru=geo.Xzero(1:3);
+ olddirv=geo.Xzero(4:6);
+ oldpos=geo.Xzero(7:9).*geo.pixelsize; % mm
+else
+ olddiru=detdiru;
+ olddirv=detdirv;
+ oldpos=geo.detpos0;%.*geo.pixelsize
+end
+
+
+% Image frame
+plot([0 detusize+1],[0 0],'k')
+plot([0 detusize+1],[detvsize+1 detvsize+1],'k')
+plot([0 0],[0 detvsize+1],'k')
+plot([detusize+1 detusize+1],[0 detvsize+1],'k')
+
+%% Midlines
+plot([0 detusize+1],[detvsize/2+0.5 detvsize/2+0.5],'-.y','Linewidth',1)
+plot([detusize/2+0.5 detusize/2+0.5],[0 detvsize+1],'-.y','Linewidth',1)
+
+Qdet = gtFedDetectorProjectionTensor(detdiru,detdirv,1,1);
+
+%% Beam direction in image
+beamuv = Qdet*beamdir;
+% Arrow in figure indicating beam direction
+
+if (dot(beamdir,detdir)-1)<0.001 % // 1
+ plot(uvorig(1),uvorig(2),'.w','MarkerSize',20)
+elseif (dot(beamdir,detdir)+1)<0.001 % // -1
+ plot(uvorig(1),uvorig(2),'xw','MarkerSize',20)
+else
+ quiver(uvorig(1),uvorig(2),beamuv(1),beamuv(2),1000,'-w','Linewidth',2)
+end
+impixelinfo
+
+%% origin u,v
+plot(0,0,'.y','MarkerSize',20)
+plot(0,0,'oy','MarkerSize',10)
+
+% unit vectors u,v
+quiver(0,0,0,1,150,'-y','Linewidth',2)
+quiver(0,0,1,0,150,'-y','Linewidth',2)
+
+%% center displacement
+shift=oldpos-detpos;
+newcenu=uvorig(1)+dot(shift,olddiru)./geo.pixelsize;
+newcenv=uvorig(2)+dot(shift,olddirv)./geo.pixelsize;
+plot(newcenu,newcenv,'+y','MarkerSize',20)
+
+
+for i=1:N
+ plot(uv(i,1),uv(i,2),['o' color],'MarkerSize',10)
+end
+
+% title
+title(titlename)
+
+
+%% noise
+
+if ~isempty(numnoise)
+
+ rnoise=(rand(N,size(uv,2))-0.5)*noise*2;
+ ind=randperm(N);
+ color=rand(1,3);
+ newuv(ind,:)=rnoise(ind,:)+uv(ind,:);
+ plot(newuv(ind,1),newuv(ind,2),['o' color],'MarkerSize',10)
+ noiseuv=newuv;
+
+else
+ noiseuv=uv;
+end
+
+
+if debug
+ disp('detpos:')
+ disp(detpos)
+ disp('oldpos:')
+ disp(oldpos)
+ disp('shift:')
+ disp(shift)
+ disp('beamuv:')
+ disp(beamuv)
+end
+
+end
+
+
\ No newline at end of file
diff --git a/zUtil_Optimization/gtTuneDetector3.m b/zUtil_Optimization/gtTuneDetector3.m
new file mode 100644
index 0000000000000000000000000000000000000000..caa805c30b4dc831fece6aee3b8fe8818d1ce32a
--- /dev/null
+++ b/zUtil_Optimization/gtTuneDetector3.m
@@ -0,0 +1,200 @@
+function newgeo=gtTuneDetector3(geometry,grainFed,uvdb,induv,debug,varargin)
+%% newgeo=gtTuneDetector3(geometry,grainFed,uvdb,induv,debug,varargin)
+% INPUT:
+% grainFed = output from gtFedGenerateRandomGrain2
+% geometry = output from gtOptimizeDetector
+% uvdb = uv from database
+% induv = indexes of uvdb to be used from the full list of
+% grainFed.allblobs.uv
+% debug = print comments
+% varargin = optional arguments for options
+% options = optimset(varargin);
+% DEFAULT: options =
+% optimset('TolFun',1e-2,'MaxFunEvals',3000,'MaxIter',700)
+%
+% OUTPUT:
+% newgeo = new geometry
+%
+%
+% OPTIMSET PARAMETERS for MATLAB
+% Display - Level of display [ off | iter | notify | final ]
+% MaxFunEvals - Maximum number of function evaluations allowed
+% [ positive integer ]
+% MaxIter - Maximum number of iterations allowed [ positive scalar ]
+% TolFun - Termination tolerance on the function value [ positive scalar ]
+% TolX - Termination tolerance on X [ positive scalar ]
+% FunValCheck - Check for invalid values, such as NaN or complex, from
+% user-supplied functions [ {off} | on ]
+% OutputFcn - Name(s) of output function [ {[]} | function ]
+% All output functions are called by the solver after each
+% iteration.
+% PlotFcns - Name(s) of plot function [ {[]} | function ]
+% Function(s) used to plot various quantities in every iteration
+%
+
+tic
+disp('Start gtTuneDetector3...')
+
+if ~exist('debug','var') || isempty(debug)
+ debug=false;
+end
+
+%% get info from newgr and grainFed
+newgeo=geometry;
+
+pixelsize = geometry.pixelsize;
+detdiru = geometry.detdiru0;
+detdirv = geometry.detdirv0;
+detpos = geometry.detpos0
+check=inputwdefault('Is the detpos correct (unit length must be mm)? y/n','y');
+if strcmpi(check,'n')
+ detpos=input('Insert the right value of detpos between []: ');
+else
+
+end
+detpos=detpos./pixelsize;
+uvorig = [geometry.detucen0, geometry.detvcen0]';
+omstep = 180/geometry.nproj;
+pixelsize = geometry.pixelsize;
+
+if ~isfield(geometry,'detdir0')
+ detdir = cross(detdiru,detdirv);
+else
+ detdir = geometry.detdir0;
+end
+
+if dot(detdiru,detdirv)>1e-4
+ disp('u,v are not perpendicular. Quitting...')
+ return
+end
+
+% grain.allblobs.[variable]
+srot = grainFed.allblobs.srot;
+dvec = grainFed.allblobs.dvec;
+omega = grainFed.allblobs.omega;
+csam = grainFed.center';
+
+uvwfed = grainFed.allblobs.uvw;
+
+%% method minimum squares
+
+% parameters for lsqcurvefit
+% ux uy uz vx vy vz dx dy dz
+X0=[detdiru detdirv detpos];
+
+% avoid NAN difspot entries
+% % ind=find(~isnan(uvwfed(:,1))); % empty u coord.
+% % dvec=dvec(ind,:);
+% % omega=omega(ind,:);
+% % uvwfed=uvwfed(ind,:);
+
+% apply induv to the lists with size (Nx3,1)
+dvec=dvec(induv,:);
+omega=omega(induv,:);
+
+pippo=[];
+for i=1:size(dvec,1)
+ j=induv(i);% was ind(i)
+ pippo=[pippo;srot(j*3-2:j*3,:)];
+end
+srot=pippo;
+
+%srot=srot(induv,:);
+
+xdata=[srot;dvec];
+ydata=uvdb; % uv coordinates found in the database (Nx2)
+
+
+if debug
+ size(srot)
+ size(dvec)
+ size(omega)
+ size(xdata)
+ size(ydata)
+end
+
+newgeo.uvwfed=uvwfed;
+newgeo.Xzero=X0;
+
+% calculate the old uvw coordinates from the current geometry
+[Fzero,Mzero] = gtFindUVW3(X0,xdata,csam,uvorig);
+Fzero(:,3)=mod(omega,360)/omstep;
+newgeo.Fzero=Fzero;
+% current detector reference system
+newgeo.Mzero=Mzero;
+
+% comparison between f(x,xdata) and ydata
+% min sum [(f(x,xdata,par)-ydata)^2]
+% x
+%
+% input: X0 = the current geometry parameters 1x12
+% xdata = (Nx4,3) list from [srot;dvec]
+% ydata = [bbox.xcen,bbox.ycen]
+% par = list of constants in the function (csam, detpos)
+%
+% output: X = the optimized geometry parameters from lsq
+
+
+% change some options for the fit
+if ~exist('varargin','var')|| isempty(varargin)
+ options = optimset('TolFun',1e-2,'MaxFunEvals',3000,'MaxIter',700);
+else
+ options = optimset(varargin);
+end
+
+disp('Start fitting...')
+f = @(x,xdata) gtFindUVW3(x,xdata,csam,uvorig);
+UB=[];
+LB=[];
+[X,a,b,c,d,e,f] = lsqcurvefit(f,X0,xdata,ydata,[],[],options);
+
+% calculate the new uvw (F) values with the optimized geometry (X)
+[F,M,Xnew] = gtFindUVW3(X,xdata,csam,uvorig);
+F(:,3)=mod(omega,360)/omstep;
+newgeo.X = X;
+newgeo.F = F;
+newgeo.M = M;
+% recalculate uvw (Fnew) with the orthonormal axes (Xnew)
+
+[Fnew,Mnew] = gtFindUVW3(Xnew,xdata,csam,uvorig);
+Fnew(:,3)=mod(omega,360)/omstep;
+
+newgeo.Xnew=Xnew;
+newgeo.Fnew=Fnew;
+% new detector reference system
+newgeo.Mnew=Mnew;
+
+% calculate the dot products
+newgeo.dotuv=dot(Mnew(1,:),Mnew(2,:));
+newgeo.dotut=dot(Mnew(1,:),Mnew(3,:));
+newgeo.dotvt=dot(Mnew(2,:),Mnew(3,:));
+
+% get the new geometry
+newgeo.detdiru0=Mnew(1,:);
+newgeo.detdirv0=Mnew(2,:);
+newpos=Xnew(7:9).*pixelsize
+check=inputwdefault('Must the detpos be reconverted into mm? y/n','n');
+if strcmpi(check,'y')
+ newpos=newpos.*pixelsize;
+end
+newgeo.detdir0=Mnew(3,:);
+
+% new center
+oldpos=X0(7:9).*pixelsize;
+newgeo.detpos0=newpos;
+olddiru=X0(1:3);
+olddirv=X0(4:6);
+shift=oldpos-newpos;
+newcenu=uvorig(1)+dot(shift,olddiru)./pixelsize;
+newcenv=uvorig(2)+dot(shift,olddirv)./pixelsize;
+
+newgeo.newdetucen0=newcenu;
+newgeo.newdetvcen0=newcenv;
+newgeo.dist=norm(newpos);
+
+toc
+
+
+
+
+end % end function
\ No newline at end of file