Modified geometry fitting algorithm; now non-linearized rotations.

Signed-off-by: Peter Reischig <peter.reischig@esrf.fr>

3_pairmatchingGUI/gtMatchFitParameters.m

 function handles = gtMatchFitParameters(handles)
 
 % Does the detector and lattice parameter fitting simultaneously for all
-% the specified parameters via the reflections marked by the user. It uses 
+% the specified parameters using the reflections marked by the user. It uses 
 % Matlab's built-in 'lsqnonlin' algorithm.
-
+%
+% The initial labgeo parameters used in the fitting are loaded from the
+% parameters file into handles.fitting.labgeobase when the GUI is launched.
+% These values are then corrected by finding the best fitting parameters.
+% Some fitting parameters are directly listed in labgeobase and some are 
+% rotations relative to labgeobase. As an initial guess for the fitting, 
+% the latest corrections are applied to labgeobase.
+% To change labgeobase to the latest (saved) version, the GUI needs 
+% restarting. In this case the rotations will be reset to 0.
 
 % Check available Optimization Toolbox License
 ok = license('checkout','Optimization_Toolbox');
@@ -138,13 +146,13 @@ allboundshigh(1:3) = handles.DetParsCorrAct(1:3) + ...
                      handles.labgeo.detsizeu*...
                      handles.labgeo.pixelsizeu;
 
-% Detector vector directions;
-allboundslow([4,6,7])  = -0.5;
-allboundshigh([4,6,7]) =  0.5;
+% Detector vector rotations in degrees;
+allboundslow([4,6,7])  = -30;
+allboundshigh([4,6,7]) =  30;
 
 % Detector vectors U-V angle in degrees
-allboundslow(5)  = 45;
-allboundshigh(5) = 135;
+allboundslow(5)  = 60;
+allboundshigh(5) = 120;
 
 % Pixel sizes
 allboundslow(8)  = handles.DetParsCorrAct(8)*0.5;
@@ -175,8 +183,10 @@ disp('PARAMETERS FITTING')
     @(newvars) sfDeviation(newvars,newvarsind,oldvars,labgeobase,cryst,...
     lambdamode,lambda_true,pairs),...
     varsguess, boundslow, boundshigh,...
-    optimset('TolX',handles.fitting.tolx, 'TolFun',handles.fitting.tolfun,...
-    'TypicalX',handles.fitting.typicalx(newvarsind),'FunValCheck','off', 'Display','iter'));
+    optimset('TolX',handles.fitting.tolx, 'TolFun',handles.fitting.tolfun, ...
+    'DiffMinChange',handles.fitting.derstep, 'DiffMaxChange',handles.fitting.derstep, ...
+    'FunValCheck','off', 'Display','iter'));
+    %'TypicalX',handles.fitting.typicalx(newvarsind),'FunValCheck','off', 'Display','iter'));
 
 % if any(exitflag==[1,2,3,4])
 %     disp(' ')

3_pairmatchingGUI/gtMatchInitialize.m

@@ -165,14 +165,14 @@ handles.fitting.labgeobase.rotdir  = handles.parameters.labgeo.rotdir;
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 % Row names in table
-rownames{1} = 'Position X';
-rownames{2} = 'Position Y';
-rownames{3} = 'Position Z';
-rownames{4} = 'Tilt in plane';
+rownames{1} = 'Position X (mm)';
+rownames{2} = 'Position Y (mm)';
+rownames{3} = 'Position Z (mm)';
+rownames{4} = 'Tilt in plane (deg)';
 rownames{5} = 'Angle U-V (deg)';
-rownames{6} = 'Tilt around U';
-rownames{7} = 'Tilt around V';
-rownames{8} = 'Pixel size mean';
+rownames{6} = 'Tilt around U (deg)';
+rownames{7} = 'Tilt around V (deg)';
+rownames{8} = 'Pixel size mean (um)';
 rownames{9} = 'Pixel size ratio';
 
 set(handles.DetPars_Table,'RowName',rownames);
@@ -182,9 +182,9 @@ set(handles.DetPars_Table,'RowName',rownames);
 handles.DetParsCorrSaved    = zeros(1,9);
 handles.DetParsCorrSaved([1,2,3]) = handles.labgeo.detrefpos(1:3);
 handles.DetParsCorrSaved([4,6,7]) = 0;
-handles.DetParsCorrSaved(5) = acosd(handles.labgeo.detdiru*handles.labgeo.detdirv'); %handles.labgeo.uvangle;
-handles.DetParsCorrSaved(8) = 0.5*(handles.labgeo.pixelsizeu+handles.labgeo.pixelsizev); %handles.labgeo.pixelsizemean;
-handles.DetParsCorrSaved(9) = handles.labgeo.pixelsizeu/handles.labgeo.pixelsizev; %handles.labgeo.pixelsizeratio;
+handles.DetParsCorrSaved(5) = acosd(handles.labgeo.detdiru * handles.labgeo.detdirv');
+handles.DetParsCorrSaved(8) = 1000*0.5*(handles.labgeo.pixelsizeu + handles.labgeo.pixelsizev);
+handles.DetParsCorrSaved(9) = handles.labgeo.pixelsizeu / handles.labgeo.pixelsizev;
 
 % Actual and last fitted values
 handles.DetParsCorrLast = handles.DetParsCorrSaved;
@@ -207,14 +207,14 @@ set(handles.DetPars_Table,'Data',dettable);
 % Detector position
 handles.DetParsIncrement([1,2,3]) = handles.labgeo.detsizeu*handles.labgeo.pixelsizeu*0.01;
 
-% Tilts vectors increment
-handles.DetParsIncrement([4,6,7]) = 0.05; 
+% Tilts increment
+handles.DetParsIncrement([4,6,7]) = 0.5; 
 
 % UV angle in degrees
 handles.DetParsIncrement(5) = 0.5;
 
 % Pixel sizes
-handles.DetParsIncrement(8) = handles.labgeo.pixelsizeu*0.1;
+handles.DetParsIncrement(8) = handles.DetParsCorrSaved(8)*0.1;
 handles.DetParsIncrement(9) = 0.1;
 
 
@@ -267,21 +267,24 @@ handles.LattParsIncrement(4:6) = 0.5;
 % 15 = 'beta';
 % 16 = 'gamma';
 
-% Typical values for fitting to estimate gradients 
-handles.fitting.typicalx([1,2,3])    = handles.labgeo.pixelsizeu; % detector position
-handles.fitting.typicalx([4,6,7])    = 0.001;  % detector vector orientations
-handles.fitting.typicalx(5)          = 0.1;    % detector U-V angle in degrees
-handles.fitting.typicalx(8)          = handles.labgeo.pixelsizeu*0.001;  % mean pixel size
-handles.fitting.typicalx(9)          = 0.001;  % pixel size ratio
-handles.fitting.typicalx([10,11,12]) = 0.001;  % unit cell size in Angstrom
-handles.fitting.typicalx([13,14,15]) = 0.1;    % unit cell angles in degrees
+% Typical values of fitted parameters to scale finite differences 
+% handles.fitting.typicalx([1,2,3])    = 1;     % detector position in mm
+% handles.fitting.typicalx([4,6,7])    = 0.1;   % detector vector rotations in degrees
+% handles.fitting.typicalx(5)          = 0.1;   % detector U-V angle in degrees
+% handles.fitting.typicalx(8)          = 1;     % mean pixel size in microns
+% handles.fitting.typicalx(9)          = 1;     % pixel size ratio
+% handles.fitting.typicalx([10,11,12]) = 1;     % unit cell size in Angstrom
+% handles.fitting.typicalx([13,14,15]) = 10;    % unit cell angles in degrees
 
-% Termination tolerance on fitted parameters (??)
-handles.fitting.tolx = 1e-7;
+% Termination tolerance on fitted parameters
+handles.fitting.tolx = 1e-6;
 
 % Termination tolerance on the fitting function value
 handles.fitting.tolfun = 1e-10;
 
+% Step size of finite differences to calculate derivatives
+handles.fitting.derstep = 1e-7;
+
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% Define blank plot handles 
@@ -317,11 +320,24 @@ handles = gtMatchLoadDifspotData(handles);
 %% Theoretical Bragg angles - load crystallographic data
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
+% X-ray beam wavelength in Angstroms 
+lambda_true = gtConvEnergyToWavelength(handles.parameters.acq.energy);
+
 % {hkl}-s, and theoretical Bragg angles from energy and crystal structure
 for ip = 1:handles.nof_phases
-    handles.CompPhaseThetas{ip}  = double(parameters.cryst(ip).theta);
+    
+    %handles.CompPhaseThetas{ip}  = double(parameters.cryst(ip).theta);
     handles.CompPhaseHKLs{ip}    = double(parameters.cryst(ip).hkl);
     handles.CompIntensities{ip}  = double(parameters.cryst(ip).int);
+    
+    % 'B matrix': HLK -> Cartesian transformation
+    Bmat = gtCrystHKL2CartesianMatrix(handles.LattParsCorrAct(ip,:));
+    % d-spacings
+    dsp_comp = gtCrystDSpacing(handles.CompPhaseHKLs{ip}, Bmat);
+    % Bragg angle theta; avoid using the saved theta values in
+    % parameters.cryst, so that the angles surely correspond to the energy
+    handles.CompPhaseThetas{ip} = asind(0.5*lambda_true./dsp_comp);
+
     handles.CompThetaIndices{ip} = 1:length(handles.CompPhaseThetas{ip});
     handles.CompThetasTicked{ip} = true(1,length(handles.CompPhaseThetas{ip}));
     
@@ -334,16 +350,22 @@ for ip = 1:handles.nof_phases
   
 end
 
-% Tolerance for distinguishing unique thetas in degrees
-handles.ThetaUniqueTol = 1e-6;
-
+% Tolerance for distinguishing unique {hkl} families in degrees
+if isfield(parameters.match,'uniquetheta')%
+    handles.ThetaUniqueTol = parameters.match.uniquetheta;
+else
+    handles.ThetaUniqueTol = 0.1;
+end
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% Theta filter limits
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-handles.ThetaLimits = 0.5*[handles.labgeo.detanglemin, ...
-                           handles.labgeo.detanglemax];
+if ~isempty(parameters.match.thetalimits)
+    handles.ThetaLimits = parameters.match.thetalimits;
+else
+    handles.ThetaLimits = 0.5*[handles.labgeo.detanglemin, ...
+                               handles.labgeo.detanglemax];
+end
 handles.ThetaLimLastLow  = handles.ThetaLimits(1);
 handles.ThetaLimLastHigh = handles.ThetaLimits(2);
 

3_pairmatchingGUI/gtMatchRecalculateLabgeo.m

@@ -3,24 +3,27 @@ function labgeo = gtMatchRecalculateLabgeo(newvars,newvarsind,oldvars,...
 
 % Recalculates labgeo parameters based on new and old values relative to
 % labgeobase. New values come from the user or from fitting.
+% It uses all the 'newvars' that are specified. For those parameters which 
+% have no new value specified, it uses the 'oldvars'. 'Labgeo' is then 
+% recalculated according to this complete list relative to 'labgeobase'.
 %
 % List of parameters in oldvars: 
 %
-% 1  = 'Position X';
-% 2  = 'Position Y';
-% 3  = 'Position Z';
-% 4  = 'Tilt in plane';
-% 5  = 'Angle U-V';
-% 6  = 'Tilt around U';
-% 7  = 'Tilt around V';
-% 8  = 'Pixel size mean';
-% 9  = 'Pixel size ratio';
-% 10 = 'lengtha';
-% 12 = 'lengthb';
-% 13 = 'lengthc';
-% 14 = 'alpha';
-% 15 = 'beta';
-% 16 = 'gamma';
+% 1  = Position X (mm)
+% 2  = Position Y (mm)
+% 3  = Position Z (mm)
+% 4  = Tilt in plane (deg)
+% 5  = Angle U-V     (deg)
+% 6  = Tilt around U (deg)
+% 7  = Tilt around V (deg)
+% 8  = Pixel size mean  (um)
+% 9  = Pixel size ratio (1)
+% 10 = lengtha (A)
+% 12 = lengthb (A)
+% 13 = lengthc (A)
+% 14 = alpha   (deg)
+% 15 = beta    (deg)
+% 16 = gamma   (deg)
 
 
 
@@ -32,24 +35,28 @@ actvars = oldvars;
 actvars(newvarsind) = newvars;
 
 % Detector position
-labgeo.detrefpos  = actvars(1:3);
+labgeo.detrefpos(1,:) = actvars(1:3);
 
 % Pixel sizes
-labgeo.pixelsizeu     = actvars(8)*sqrt(actvars(9));
-labgeo.pixelsizev     = actvars(8)/sqrt(actvars(9));
+labgeo.pixelsizeu = 0.001*actvars(8)*sqrt(actvars(9));
+labgeo.pixelsizev = 0.001*actvars(8)/sqrt(actvars(9));
 
-% Tilt in plane (increment U vector appr. along V)
-detdiru = labgeobase.detdiru + cross(labgeobase.detnorm,...
-                               labgeobase.detdiru)*actvars(4);
+% Tilts around U,V and detector normal
+rmcU = gtMathsRotationMatrixComp(labgeobase.detdiru, 'row');
+rmcV = gtMathsRotationMatrixComp(labgeobase.detdirv, 'row');
+rmcN = gtMathsRotationMatrixComp(labgeobase.detnorm, 'row');
 
-% Tilt around V - right handed rotation
-detdiru = detdiru - labgeobase.detnorm*actvars(6);
-detdiru = detdiru/sqrt(detdiru*detdiru');
-labgeo.detdiru = detdiru;
+rotU = gtMathsRotationTensor(actvars(6), rmcU);
+rotV = gtMathsRotationTensor(actvars(7), rmcV);
+rotN = gtMathsRotationTensor(actvars(4), rmcN);
+
+% Complete rotation matrix 
+rot = rotU*rotV*rotN;
+
+% Apply rotation matrix to row vectors
+detdiru = labgeobase.detdiru*rot;
+detdirv = labgeobase.detdirv*rot;
 
-% Tilt around U - right handed rotation
-detdirv = labgeobase.detdirv + labgeobase.detnorm*actvars(7);
-detdirv = detdirv/sqrt(detdirv*detdirv');
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% Recalculate V according to angle between U and V 
@@ -58,22 +65,24 @@ detdirv = detdirv/sqrt(detdirv*detdirv');
 % Actual U-V angle
 actangle = acosd(detdiru*detdirv');
 
-% Difference vector normalised
-diffvec = detdirv - detdiru;
-diffvec = diffvec/sqrt(diffvec*diffvec');
+% Update detector normal
+detnorm = cross(detdiru, detdirv);
+detnorm = detnorm/sqrt(detnorm*detnorm');
 
-% Correction vector along diffvec
-corrvec = diffvec*sind(actvars(5)-actangle)/...
-          cosd(actvars(5)-actangle/2);
+% Rotation matrix to adjust angle between U-V
+rmcN = gtMathsRotationMatrixComp(detnorm, 'row');
+rotN = gtMathsRotationTensor(0.5*(actvars(5) - actangle), rmcN);
 
-% Apply correction
-detdirv = detdirv + corrvec;
+% Apply correction opposite ways 
+detdiru = detdiru*rotN';
+detdirv = detdirv*rotN;
 
-% Renormalise V vector
+% Renormalise U and V vector
+labgeo.detdiru = detdiru/sqrt(detdiru*detdiru');
 labgeo.detdirv = detdirv/sqrt(detdirv*detdirv');
 
 % Recalculate angle for checking (should be = actvars(5))
-%labgeo.uvangle = acosd(labgeo.detdiru*labgeo.detdirv');
+%uvangle = acosd(labgeo.detdiru*labgeo.detdirv')
 
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -81,8 +90,7 @@ labgeo.detdirv = detdirv/sqrt(detdirv*detdirv');
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 % Update detector normal
-detnorm = cross(labgeo.detdiru, labgeo.detdirv);
-labgeo.detnorm = detnorm/sqrt(detnorm*detnorm');
+labgeo.detnorm = detnorm;
 
 % Update detector origin: pixel (0,0) coordinates in LAB ref.
 labgeo.detorig = labgeo.detrefpos - ...
@@ -90,3 +98,4 @@ labgeo.detorig = labgeo.detrefpos - ...
     labgeo.detdirv*labgeobase.detrefv*labgeo.pixelsizev;
 
 
+end % of function
\ No newline at end of file