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

Cleaning. Deleted old pair matching functions in folder 3_pairmatching. 

git-svn-id: https://svn.code.sf.net/p/dct/code/trunk@508 4c865b51-4357-4376-afb4-474e03ccb993
parent f63ce796
No related branches found
No related tags found
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 0 additions and 3016 deletions
3_pairmatching/gtMATCHCorrectGeo.m deleted 100644 → 0
+ 0
− 812
View file @ f63ce796
%
% FUNCTION out=gtMATCHCorrectGeo(varargin)
%
% Gives geometrical correction for the rot. axis to detector distance,
% and detector tilts (Y and Z) according to data in a spotpair or calibration table.
%
% X-ray wavelengths are calculated from the Bragg angle of each pair and
% the input lattice parameters. Their disgttribution is optimized by least squares
% to find the three unknown parameters.
% By default (CostMode 0), the deviation of lambda-s from the input value
% that's measured by Bond's method is optimized.
% The other option (CostMode 1) is to optimize the std of the lambda-s
% ignoring the input lamda value. However, this is not recommended as it
% may give wrong results that look nice.
%
% If the optimization does not work out, plot the cost function and find a
% minima manually.
%
%
% INPUT: out=gtCorrectGeo(varargin)
%
% Parameter names and default values:
%
% OptPar = 'all' - parameters to be optimized
% all - both tilts and distance (recommended)
% tilts - only the two tilts, distance is fixed
% dist - only the distance
% PlaneFams = [] - which plane families to use ([1,2,5,6,etc.])
% ThetaBounds = [] - boundaries of restricted theta region to use
% (manual rejection of outliers)
% ShowResults = true - show graphical output
% Parameters = [] - all parameters (otherwise parameter file is used)
% GuessTiltY = 0 - guess value for tiltY
% GuessTiltZ = 0 - guess value for tiltZ
% GuessDist = 1 - guess value for distance factor (multiplicator of input value)
% BoundsTilts = 5 - boundary of search region for tilts (+-this value in deg)
% BoundsDist = 0.3 - boundary of search region for distance factor (1+-this value)
% CostMode = 0 - type of cost function to apply
% 0 - lambda deviation from input value (recommended)
% 1 - lambda standard deviation (doesn't use input value)
% Attention! This may give wrong results that look nice!
% TolX = 1e-6 - optimization end criterion (accuracy) on parameters
% TolFun = 1e-9 - optimization end criterion (differential) on cost function
% Calibration = true - true - uses calibration table
% false - uses normal spotpair table
% ShowCostFun = false - calculates and shows the cost function on the
% region [BoundsTilts,BoundsTilts,BoundDist]
% IncTilts = 0.2 - increment of tilts for showing the cost function
% IncDist = 0.02 - increment of distance factor for showing the cost function
%
%
% OUTPUT: fields in the out variable:
%
% costfunvals = cost function values as a volume, if it was calculated
% tiltY = tilt around axis Y in Lab system (in degrees)
% tiltZ = tilt around axis Z in Lab system (in degrees)
% rottodet_ = detector distance in pixel/mm
% rottodet_corr_ = correction for the distance
% (precise value = measured + correction)
% par = parameters used in the optimization
% optimout = output of the matlab optimization routine
%
% Peter Reischig, ESRF, 08/2009
%
%%
function out=gtMATCHCorrectGeo(varargin)
par.OptPar = 'all'; % parameters to be optimized ('tilts'/'dist'/'all')
par.PlaneFams = []; % which plane families to use ([1,2,5,6,etc.])
par.ThetaBounds = []; % boundaries of restricted theta region to use
par.ShowResults = true; % show graphical output
par.Parameters = []; % parameters like in parameter file
par.GuessTiltY = 0; % guess value for tiltY
par.GuessTiltZ = 0; % guess value for tiltZ
par.GuessDist = 1; % guess value for distance factor
par.BoundsTilts = 5; % boundary of search region for tilts (+-this value in deg)
par.BoundsDist = 0.3; % boundary of search region for distance factor (1+-this value)
par.CostMode = 0; % type of cost function to apply
par.TolX = 1e-6; % optimization end criterion (accuracy) on parameters
par.TolFun = 1e-9; % optimization end criterion (differential) on cost function
par.Calibration = true; % true - uses calibration table; false - uses normal spotpair table
par.ShowCostFun = false; % calculates and shows the cost function on the region [BoundsTilts,BoundsTilts,BoundDist]
par.IncTilts = 0.2; % increment of tilts for showing the cost function
par.IncDist = 0.02; % increment of distance factor for showing the cost function
par=parse_pv_pairs(par,varargin);
if isempty(par.Parameters)
load('parameters.mat');
else
parameters = par.Parameters;
end
if isempty(par.ThetaBounds) && isfield('parameters','match_calib') %,'theta_bounds')
par.ThetaBounds = parameters.match_calib.theta_bounds;
end
disp('Parameters used:')
disp(par)
tic
measured_rottodet_mm = parameters.acq.dist
measured_rottodet_pix = parameters.acq.dist/parameters.acq.pixelsize
rotx = parameters.acq.rotx
% Set sample coordinate system
% sorX = on the rotation axis by definition
% sorY = on the rotation axis by definition
% sorZ = floor(parameters.acq.ydet/2) % at the geometrical center of the image
if isfield(parameters.acq,'sampleorigZ')
sorZ = parameters.acq.sampleorigZ;
else
sample = gtSampleGeoInSampleSystem(parameters);
sorZ = sample.sorZim;
end
if isfield(parameters.acq, 'difA_name')
difspottable=[parameters.acq.difA_name 'difspot'] ; % difspot table A set
else
difspottable=[parameters.acq.name 'difspot'] ;
end
if par.Calibration
pairtable=parameters.acq.calib_tablename;
else
pairtable=parameters.acq.pair_tablename;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% SET input for calculation
lambda=gtConvEnergyToWavelength(parameters.acq.energy) % in angstroms
reflections = parameters.cryst.hklsp; % type of reflections
d0 = parameters.cryst.dsp; % dspacing calculated from hkl
gtDBConnect();
nof_pairs= mym(sprintf('select count(*) from %s',pairtable)) ;
disp('Reading data from database...')
if isempty(par.PlaneFams)
[centXimA,centYimA,th,eta,om,thetatype]=mym(sprintf(['SELECT '...
'%s.CentroidX, %s.CentroidY, %s.theta, %s.eta, %s.omega, %s.thetatype '...
'FROM (%s INNER JOIN %s ON %s.difspotID=%s.difAID) ORDER BY %s.pairID'], ...
difspottable,difspottable,pairtable,pairtable,pairtable,pairtable,difspottable,...
pairtable,difspottable,pairtable,pairtable)) ;
[centXimB,centYimB]=mym(sprintf(['SELECT %s.CentroidX, %s.CentroidY ' ...
'FROM (%s INNER JOIN %s ON %s.difspotID=%s.difBID) ORDER BY %s.pairID'], ...
difspottable,difspottable,difspottable,...
pairtable,difspottable,pairtable,pairtable)) ;
else
centXimA = [] ;
centYimA = [] ;
centXimB = [] ;
centYimB = [] ;
th = [] ;
eta = [] ;
om = [] ;
thetatype = [] ;
for ii=1:length(par.PlaneFams)
[centXimA_r,centYimA_r,th_r,eta_r,om_r,thetatype_r]=mym(sprintf(['SELECT '...
'%s.CentroidX, %s.CentroidY, %s.theta, %s.eta, %s.omega, %s.thetatype '...
'FROM (%s INNER JOIN %s ON %s.difspotID=%s.difAID) ' ...
'WHERE %s.thetatype=%d ORDER BY %s.pairID'], ...
difspottable,difspottable,pairtable,pairtable,pairtable,pairtable,difspottable,...
pairtable,difspottable,pairtable,pairtable,par.PlaneFams(ii),pairtable)) ;
[centXimB_r,centYimB_r]=mym(sprintf(['SELECT %s.CentroidX, %s.CentroidY ' ...
' FROM (%s INNER JOIN %s ON %s.difspotID=%s.difBID)' ...
' WHERE %s.thetatype=%d ORDER BY %s.pairID'], ...
difspottable,difspottable,difspottable,...
pairtable,difspottable,pairtable,pairtable,par.PlaneFams(ii),pairtable)) ;
centXimA = [centXimA; centXimA_r];
centYimA = [centYimA; centYimA_r];
centXimB = [centXimB; centXimB_r];
centYimB = [centYimB; centYimB_r];
th = [th; th_r];
eta = [eta; eta_r];
om = [om; om_r];
thetatype = [thetatype; thetatype_r];
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Apply theta boundaries
if ~isempty(par.ThetaBounds)
todel=[];
for i=1:length(thetatype)
if th(i)<par.ThetaBounds(thetatype(i),1) || th(i)>par.ThetaBounds(thetatype(i),2)
todel=[todel; i];
end
end
centXimA(todel) = [];
centYimA(todel) = [];
centXimB(todel) = [];
centYimB(todel) = [];
th(todel) = [];
eta(todel) = [];
om(todel) = [];
thetatype(todel) = [];
nof_pairs = length(thetatype);
end
disp('Number of pairs used:')
disp(nof_pairs)
disp('Processing data...')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Display the cost function
if par.ShowCostFun
sizetilt = floor(par.BoundsTilts/par.IncTilts*2) + 1;
sizedist = floor(par.BoundsDist/par.IncDist*2) + 1;
ts = -par.BoundsTilts + (0:sizetilt-1)*par.IncTilts;
ds = 1 -par.BoundsDist + (0:sizedist-1)*par.IncDist;
disp(' ')
disp('Calculating the cost function over the specified range... ')
disp(' tilts values:')
disp(ts)
disp(' ')
disp(' distance factor values:')
disp(ds)
disp(' ')
vs = zeros(sizetilt,sizetilt,sizedist);
for i = 1:sizetilt % tiltY
for j= 1:sizetilt % tiltZ
for k= 1:sizedist % dist factor
var(1) = ts(i);
var(2) = ts(j);
var(3) = ds(k);
cvec = sfCostAll(var,d0,centXimA,centYimA,centXimB,centYimB,...
rotx,sorZ,measured_rottodet_pix,lambda,par.CostMode);
vs(i,j,k) = sqrt(sum(cvec.^2));
end
end
disp(sprintf(' %f percent is done.',i/sizetilt*100))
end
disp(' ')
disp('Coordinates in output volume: 1 - tiltY; 2 - tiltZ; 3 - dist')
disp(' ')
disp('Minimum location of cost function:')
[minval,minloc]=min(vs(:));
[mini,minj,mink]=ind2sub(size(vs),minloc);
disp([mini,minj,mink])
min_tiltY = ts(mini);
min_tiltZ = ts(minj);
min_dist = ds(mink);
disp(sprintf(' tiltY = %f',min_tiltY))
disp(sprintf(' tiltZ = %f',min_tiltZ))
disp(sprintf(' dist = %f',min_dist))
gtVolView(vs)
set(gcf,'name','Cost function')
out.costfunvals = vs;
else
out.costfunvals = [];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Parameters to optimize
% Optimize for the two tilts
if strcmpi(par.OptPar,'tilt')
[opt,fval,res,exitflag,optimout] = lsqnonlin(@(var) sfCostTilts(var,d0,...
centXimA,centYimA,centXimB,centYimB,...
rotx,sorZ,measured_rottodet_pix,lambda,par.CostMode),...
[ par.GuessTiltY, par.GuessTiltZ ],...
[-par.BoundsTilts, -par.BoundsTilts ],...
[ par.BoundsTilts, par.BoundsTilts ],...
optimset('TolX',par.TolX,'TolFun',par.TolFun,'Display','iter'));
if any(exitflag==[1,2,3])
disp(' ')
disp('An optimimum for correction has been found.')
disp(' ')
else
disp(' ')
disp('No optimized values for correction were found.')
disp(' ')
opt.optimout = optimout;
return
end
tiltY = opt(1) ; % tilt around Y in degrees
tiltZ = opt(2) ; % tilt around Z in degrees
rottodet_mm = measured_rottodet_mm;
% Optimize for distance (rottodet)
elseif strcmpi(par.OptPar,'dist')
tiltY = parameters.match.tiltY;
tiltZ = parameters.match.tiltZ;
[opt,fval,exitflag,optimout] = fminbnd(@(var) sfCostDist(var,d0,...
centXimA,centYimA,centXimB,centYimB,...
rotx,sorZ,measured_rottodet_pix,tiltY,tiltZ,lambda,par.CostMode),...
1-par.BoundsDist, 1+par.BoundsDist,...
optimset('TolX',par.TolX,'Display','iter'));
if exitflag==1
disp(' ')
disp('An optimimum for correction has been found.')
disp(' ')
else
disp(' ')
disp('No optimized value for correction was found.')
disp(' ')
opt.optimout = optimout;
return
end
rottodet_mm = opt*measured_rottodet_pix*parameters.acq.pixelsize;
% Optimize both for the two tilts and the distance
elseif strcmpi(par.OptPar,'all')
[opt,fval,res,exitflag,optimout] = lsqnonlin(@(var) sfCostAll(var,d0,...
centXimA,centYimA,centXimB,centYimB,...
rotx,sorZ,measured_rottodet_pix,lambda,par.CostMode),...
[ par.GuessTiltY, par.GuessTiltZ, par.GuessDist],...
[-par.BoundsTilts, -par.BoundsTilts, 1-par.BoundsDist],...
[ par.BoundsTilts, par.BoundsTilts, 1+par.BoundsDist],...
optimset('TolX',par.TolX,'TolFun',par.TolFun,'Display','iter'));
if any(exitflag==[1,2,3])
disp(' ')
disp('An optimimum for correction has been found.')
disp(' ')
else
disp(' ')
disp('No optimized values for correction were found.')
disp(' ')
opt.optimout = optimout;
return
end
tiltY=opt(1) ;
tiltZ=opt(2) ;
rottodet_mm=opt(3)*measured_rottodet_pix*parameters.acq.pixelsize;
else
error('Optimization parameter can be: ''all''/''tilt''/''dist'' ')
end
rottodet_pix=rottodet_mm/parameters.acq.pixelsize ;
corr_rottodet_pix=rottodet_pix-measured_rottodet_pix ;
corr_rottodet_mm=corr_rottodet_pix*parameters.acq.pixelsize ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Calculation of optimized values
for i=1:length(centXimA)
dv = gtMatchDiffVecInLab(centXimA(i),centYimA(i),...
centXimB(i),centYimB(i),rottodet_pix,tiltY,tiltZ,rotx,sorZ);
th_op(i,1) = gtMatchThetaOfPair(dv);
eta_op(i,1) = gtMatchEtaOfPoint(dv);
end
lambda_op = d0.*sind(th_op)*2;
for i=1:max(thetatype)
nof_thetatype(i) = sum(thetatype==i);
avgthetas_op(i) = mean(th_op(thetatype==i));
end
% psi_op = d0.*sind(th_op) ;
% std_psi_op = std(psi_op) ;
%
% d_op = lambda./(2.*psi_op) ; % lattice parameter
% std_d_op = std(d_op) ;
% opt_d = lambda/2/mean(psi_op) ;
%
% for i=1:size(reflections,1)
% w(i)=1./norm(reflections(i,:)) ; % plane geometry (in angstroms)
% end
%
% opt_thetas=asind(mean(psi_op)./w) ; % optimized diffraction angle for each plane family
%
% mean_psi_op=mean(psi_op);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Calculation of original values from the database for comparison
% psi_0 = v.*sind(th) ;
% std_psi_0 = std(psi_0) ;
% d_0 = lambda./(2.*psi_0) ; % original lattice parameter
% std_d_0 = std(d_0) ;
lambda_0 = d0.*sind(th)*2;
for i=1:max(thetatype)
avgthetas_0(i) = mean(th(thetatype==i));
end
disp('Getting twotheta values from parameters file...')
valid_2thetas = parameters.cryst.twoth ; % coloumn vector of all valid 2theta angles in degrees
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Create output
out.tiltY = tiltY;
out.tiltZ = tiltZ;
out.rottodet_pix = rottodet_pix;
out.rottodet_corr_pix = corr_rottodet_pix;
out.rottodet_mm = rottodet_mm;
out.rottodet_corr_mm = corr_rottodet_mm;
out.par = par;
out.optimout = optimout;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Display results
disp(' ')
disp(' ')
disp('Number of pairs used for correction:'), disp(length(thetatype)) ;
disp(' Number of pairs per plane-family:'), disp(nof_thetatype);
% disp('Wavelength in angstroms:')
% disp(lambda)
disp('Final value of the objective function:')
disp(fval)
% disp('Mean of psi (dimensionless) [input; database; optimized]:')
% disp([lambda/2/latticepar;mean(psi_0);mean_psi_op])
%
% disp('Std of psi:')
% disp([0;std_psi_0;std_psi_op])
%
% disp('Mean of lattice parameter (angstrom) [input; database; optimized]:')
% disp([latticepar;mean(d_0);opt_d])
%
% disp('Std of lattice parameter:')
% disp([0;std_d_0;std_d_op])
disp('Mean of lambda [input; database; optimized]:')
disp([lambda; mean(lambda_0); mean(lambda_op)])
disp('Std of lambda:')
disp([0; std(lambda_0); std(lambda_op)])
disp('Mean diffraction angles (deg) [computed; database; optimized]:')
disp(valid_2thetas/2)
disp(avgthetas_0)
disp(avgthetas_op)
disp('Rotation axis to detector distance (mm) [input; optimized; difference]:')
disp([measured_rottodet_mm; rottodet_mm; corr_rottodet_mm])
disp('Rotation axis to detector distance (pixel) [input; optimized; difference]:')
disp([measured_rottodet_pix; rottodet_pix; corr_rottodet_pix])
disp('Tilt around axis Y in lab system (deg):')
disp(tiltY)
disp('Tilt around axis Z in lab system (deg):')
disp(tiltZ)
disp(' ')
toc
disp(' ')
% disp('Input psi (dimensionless):')
% disp(lambda/2/latticepar)
% disp('Mean of optimized psi:')
% disp(mean_psi_op)
% disp('Mean of originally measured psi:')
% disp(mean(psi_0))
% disp('Std of originally measured psi (lambda/2d):')
% disp(std_psi_0)
% disp('Std of optimized psi:')
% disp(std_psi_op)
% disp('Input lattice parameter (angstrom):')
% disp(latticepar)
% disp('Mean of originally measured lattice parameter:')
% disp(mean(d_0))
% disp('Mean of optimized lattice parameter:')
% disp(opt_d)
% disp('Std of originally measured lattice parameter:')
% disp(std_d_0)
% disp('Std of optimized lattice parameter:')
% disp(std_d_op)
% disp('Diffraction angles computed from input energy and lattice parameter (deg):')
% disp(valid_2thetas/2)
% disp('Mean of originally measured diffraction angles:')
% disp(thetas_0)
% disp('Mean of optimized diffraction angles:')
% disp(opt_thetas)
% disp('Input rotation axis to detector distance (mm):')
% disp(measured_rottodet_mm)
% disp('Optimized rotation axis to detector distance (mm):')
% disp(rottodet_mm)
% if strcmpi(par.OptPar,'dist')
% disp('Optimization range (mm):')
% disp([lim_rottodet1 lim_rottodet2]*measured_rottodet_mm)
% end
% disp('Correction of rotation axis to detector distance (mm):')
% disp(corr_rottodet_mm)
% disp('Input rotation axis to detector distance (pixel):')
% disp(measured_rottodet_pix)
% disp('Optimized rotation axis to detector distance (pixel):')
% disp(rottodet_pix)
% disp('Correction of rotation axis to detector distance (pixel):')
% disp(corr_rottodet_pix)
% if strcmpi(par.OptPar,'dist')
% disp('Optimization range (pixel):')
% disp([lim_rottodet1 lim_rottodet2]*measured_rottodet_pix)
% end
if par.ShowResults
figure('name','Histogram of theta values')
% max value in the histogram be > 30 - to make it nice
for i=0:9
bins_th = floor(min(th_op)*1000)/1000 : 0.0005+0.0005*i : ceil(max(th_op)*1000)/1000;
hist_th = hist(th_op, bins_th);
if max(hist_th)>30
break
end
end
hist(th_op, bins_th);
h = findobj(gca,'Type','patch');
set(h,'FaceColor','g','EdgeColor','g')
hold on ;
hist(th,bins_th) ;
xlabel('theta (deg)') ;
title('BLUE: original values from database; GREEN: optimized') ;
figure('name','theta - eta OPTIMIZED')
plot(th_op,eta_op,'.g','MarkerSize',6);
hold on ;
for i=1:length(avgthetas_op)
if ~isempty(find(par.PlaneFams==i, 1))
plot([avgthetas_op(i) avgthetas_op(i)],[0 360],'-.k','LineWidth',1);
end
end
xlabel('theta (deg)') ;
ylabel('eta (deg)') ;
figure('name','theta - eta ORIGINAL')
hold on ;
if ~isempty(par.ThetaBounds) % draw rectangles under the restricted theta ranges
for i=1:size(par.ThetaBounds,1)
if par.ThetaBounds(i,2) > par.ThetaBounds(i,1)
rectangle('Position',[par.ThetaBounds(i,1),0,...
par.ThetaBounds(i,2)-par.ThetaBounds(i,1),360],...
'FaceColor','y','EdgeColor','none')
end
end
end
for i=1:length(avgthetas_0)
if ~isempty(find(par.PlaneFams==i, 1))
plot([avgthetas_0(i) avgthetas_0(i)],[0 360],'-.k','LineWidth',1);
end
end
plot(th,eta,'.r','MarkerSize',6);
xlabel('theta (deg)')
ylabel('eta (deg)')
title('in YELLOW: restricted theta range used for optimization')
figure('name','Histogram of lambda values')
% max value in the histogram be > 30 - to make it nice
for i=0:9
bins_l = floor(min(lambda_op)*10000)/10000 : 0.00002+0.00002*i : ceil(max(lambda_op)*10000)/10000;
hist_l = hist(lambda_op, bins_l);
if max(hist_l)>30
break
end
end
hist(lambda_op,bins_l);
h = findobj(gca,'Type','patch');
set(h,'FaceColor','g','EdgeColor','g')
hold on ;
hist(lambda_0,bins_l) ;
xlabel('lambda (angstrom)') ;
title('BLUE: original lambda values from database; GREEN: optimized') ;
% figure('name','Histogram of lattice parameters')
% hist(d_op,length(thetatype)/8);
% h = findobj(gca,'Type','patch');
% set(h,'FaceColor','g','EdgeColor','g')
% hold on ;
% hist(d_0,nof_pairs/8) ;
% xlabel('lattice parameter (angstrom)') ;
% title('BLUE: original values (as in database); GREEN: optimized') ;
% figure('name','Histogram of psi values')
% hist(psi_op,length(thetatype)/8);
% h = findobj(gca,'Type','patch');
% set(h,'FaceColor','g','EdgeColor','g')
% hold on ;
% hist(psi_0,nof_pairs/8) ;
% xlabel('psi (dimensionless)') ;
% title('BLUE: original psi (=\lambda/2d) values (as in database); GREEN: optimized') ;
% figure
% plot(eta_op,d_op,'.g','MarkerSize',6);
% hold on ;
% plot([0 360],[opt_d opt_d],'-.k','LineWidth',1);
% xlabel('eta (deg)') ;
% ylabel('d (Angstroms)') ;
% title('d - eta OPTIMIZED') ;
%
% figure
% plot(eta,d_0,'.r','MarkerSize',6);
% hold on ;
% plot([0 360],[mean(d_0) mean(d_0)],'-.k','LineWidth',1);
% xlabel('eta (deg)') ;
% ylabel('d (Angstroms)') ;
% title('d - eta ORIGINAL') ;
% figure
% plot(eta_op,psi_op,'.g','MarkerSize',6);
% hold on ;
% plot([0 360],[mean(psi_op) mean(psi_op)],'-.k','LineWidth',1);
% xlabel('eta (deg)') ;
% ylabel('psi') ;
% title('psi (=\lambda/2d)- eta OPTIMIZED') ;
%
% figure
% plot(eta,psi_0,'.r','MarkerSize',6);
% hold on ;
% plot([0 360],[mean(psi_0) mean(psi_0)],'k','LineWidth',1,'LineStyle','-.');
% xlabel('eta (deg)') ;
% ylabel('psi') ;
% title('psi (=\lambda/2d) - eta ORIGINAL') ;
end % of show results
end % of function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Sub-functions in use
function cost = sfCostDist(var,d0,centUA,centVA,centUB,centVB,rotU,sorZ,...
measured_rottodet_pix,tiltY,tiltZ,lambda_input,costmode)
% Input vectors:
% var = rottodet_pix
% diffracton vectors
dv = gtMatchDiffVecInLab(centUA,centVA,centUB,centVB,...
measured_rottodet_pix*var,tiltY,tiltZ,rotU,sorZ);
% Bragg angles
theta = gtMatchThetaOfPair(dv);
% lambda-s from Friedel pairs should ideally be equal to lambda_input (Bragg's law)
lambda = d0.*sind(theta)*2;
% Value to be minimized (a scalar required for fminbnd)
switch costmode
case 0 % Deviation of lambda from input value is to be minimized
cost = sqrt(sum((lambda-lambda_input).^2));
case 1 % Standard deviation of lambda is to be minimized (doesn't use input value)
cost = std(lambda);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function cost = sfCostTilts(var,d0,centUA,centVA,centUB,centVB,rotU,sorZ,...
measured_rottodet_pix,lambda,costmode)
% Input vectors:
% var(1) = tiltY
% var(2) = tiltZ
% diffracton vectors
dv = gtMatchDiffVecInLab(centUA,centVA,centUB,centVB,...
measured_rottodet_pix,var(1),var(2),rotU,sorZ);
cost = sfCostFunction(dv,d0,lambda,costmode);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function cost = sfCostAll(var,d0,centUA,centVA,centUB,centVB,rotU,sorZ,...
measured_rottodet_pix,lambda,costmode)
% Input vectors:
% var(1) = tiltY
% var(2) = tiltZ
% var(3) = rotation axis to detector distance (same unit as inp-s)
% diffracton vectors
dv = gtMatchDiffVecInLab(centUA,centVA,centUB,centVB,...
measured_rottodet_pix*var(3),var(1),var(2),rotU,sorZ);
cost = sfCostFunction(dv,d0,lambda,costmode);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% cost function for 'tilts' and 'all'
function cost=sfCostFunction(dv,d0,lambda_input,costmode)
% Bragg angles
theta = gtMatchThetaOfPair(dv);
% lambda-s from Friedel pairs should ideally be equal to lambda_input (Bragg's law)
lambda = d0.*sind(theta)*2;
% Value to be minimized (a vector required for lsqnonlin)
switch costmode
case 0 % Deviation of lambda from input value is to be minimized
cost = lambda-lambda_input;
case 1 % Standard deviation of lambda is to be minimized (doesn't use input value)
cost = lambda-mean(lambda);
end
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function cost=sfCostFunction(dv,v,lambda,latticepar,costmode)
%
% % Bragg angles
% theta = gtMatchThetaOfPair(dv);
%
% % input psi value
% psi_input = lambda/(2*latticepar);
%
% % psi-s from Friedel pairs should ideally be equal to psi_input (Bragg's law)
% psi = v.*sind(theta);
%
%
% % Value to be minimized
%
% switch costmode
%
% case 0 % Deviation of psi from input value is to be minimized
% cost = psi-psi_input;
%
% case 1 % Deviation of lattice parameter from input value is to be minimized
% a = lambda/2./psi; % lattice parameters
% %cost = sqrt(sum((a-latticepar).^2)) ;
% cost = a-latticepar;
%
% case 2 % Standard deviation of psi is to be minimized
% cost = psi-mean(psi);
%
% case 3 % Standard deviation of lattice parameter is to be minimized
% a = lambda/2./psi; % lattice parameters
% cost = a-mean(a);
%
% end
end
\ No newline at end of file
3_pairmatching/gtMATCHDifspots.m deleted 100644 → 0
+ 0
− 831
View file @ f63ce796
This diff is collapsed.
3_pairmatching/gtMATCHEvaluate.m deleted 100644 → 0
+ 0
− 479
View file @ f63ce796
%
% FUNCTION gtMATCHEvaluate(varargin)
%
% Gives tables and figures as feedback on the matching of diffraction spots
% and matching tolerances by showing data of specified or all pairs in the
% spotpair table.
% Does not modify any data.
%
%
% USAGE
%
% gtMATCHEvaluate
% gtMATCHEvaluate([],[],'calibration',1)
% gtMATCHEvaluate([],[],'calibration',1,'table',1)
% gtMATCHEvaluate([99,100,102],parameters,'pairfigure',1,'pairdata',1,...
% 'pairdist',0,'tolhist',0,'dopause','button')
%
%
% OPTIONAL INPUT
% {default values in brackets}
%
% PairID: vector of pairID-s to be considered; if empty, all are considered {[]}
% Parameters: all parameters like in the parameters file {[]}
% PairFigure: pairs are shown graphically {false}
% PairData: table of spot properties of individual pairs {false}
% Table: table of spot properties of all specified pairs {false}
% TolHist: histogram of spot properties of pairs {true}
% PairDist: pair distribution vs omega, eta, and eta vs theta {true}
% DoPause: pause in seconds after each dispaled pair or 'button' {0}
%
%
% OUTPUT figures and tables for evaluation and comparison
%
%
% Peter Reischig, ESRF, 06/2009
%
function gtMATCHEvaluate(varargin)
% if ~exist('pairID','var')
% pairID=[];
% end
%
% if ~exist('parameters','var')
% parameters=[];
% end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Set parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
par.PairID = [];
par.Parameters = [];
par.DoPause = [];
par.PairFigure = false;
par.PairData = false;
par.Table = false;
par.TolHist = false; %true;
par.PairDist = true;
par.Calibration = false;
par=parse_pv_pairs(par,varargin);
if isempty(par.Parameters)
load('parameters.mat');
else
parameters = par.Parameters ;
end
if par.Calibration
% Thresholds for image position
thr_ang = parameters.match_calib.thr_ang; % 2theta angle
thr_ang_scale = parameters.match_calib.thr_ang_scale; % addition to thr_ang scaled wtih 2theta
thr_max_offset = parameters.match_calib.thr_max_offset; % MaxImage offset
thr_ext_offset = parameters.match_calib.thr_ext_offset; % ExtStartImage and ExtEndImage offset
%thr_genim_offset = parameters.match_calib.thr_genim_offset; % at least one of the three images should be offset as maximum this value
corr_rot_to_det = parameters.match_calib.corr_rot_to_det; % correction of rot_to_det in pixels
% Thresholds
thr_intint = parameters.match_calib.thr_intint; % integrated intensity
thr_area = parameters.match_calib.thr_area; % area of spots inside the bounding box
thr_bbsize = parameters.match_calib.thr_bbsize; % bounding box size
thr_goodness = parameters.match_calib.thr_goodness;
minsizeX = parameters.match_calib.minsizeX;
minsizeY = parameters.match_calib.minsizeY;
else
% Thresholds for image position
thr_ang = parameters.match.thr_ang; % 2theta angle
thr_ang_scale = parameters.match.thr_ang_scale; % addition to thr_ang scaled wtih 2theta
thr_max_offset = parameters.match.thr_max_offset; % MaxImage offset
thr_ext_offset = parameters.match.thr_ext_offset; % ExtStartImage and ExtEndImage offset
%thr_genim_offset = parameters.match.thr_genim_offset; % at least one of the three images should be offset as maximum this value
corr_rot_to_det = parameters.match.corr_rot_to_det; % correction of rot_to_det in pixels
% Thresholds
thr_intint = parameters.match.thr_intint; % integrated intensity
thr_area = parameters.match.thr_area; % area of spots inside the bounding box
thr_bbsize = parameters.match.thr_bbsize; % bounding box size
thr_goodness = parameters.match.thr_goodness;
minsizeX = parameters.match.minsizeX;
minsizeY = parameters.match.minsizeY;
end
sample_radius = parameters.acq.bb(3)/2;
sample_top = parameters.acq.bb(2);
sample_bot = parameters.acq.bb(2)+parameters.acq.bb(4);
rot_to_det = parameters.acq.dist/parameters.acq.pixelsize + corr_rot_to_det;
% Origin of SAMPLE coordinates
% sorX = on the rotation axis by definition
% sorY = on the rotation axis by definition
% sorZ = defined in the parameters file
if isfield(parameters.acq, 'difA_name')
difspottable=[parameters.acq.difA_name 'difspot'];
else
difspottable=[parameters.acq.name 'difspot'];
end
if par.Calibration
pairtable = parameters.acq.calib_tablename ;
else
pairtable = parameters.acq.pair_tablename ;
end
gtDBConnect
% Number of preojections per 180 degrees
nproj=parameters.acq.nproj;
% Theoretical 2theta values
comp_2thetas=parameters.cryst.twoth; % coloumn vector of all valid 2theta angles in degrees
% Pair ID-s to be used for evaluation
if isempty(par.PairID)
par.PairID=mym(sprintf('Select pairID from %s',pairtable));
end
cmp=[];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Loop through pairs; collect data; print tables
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if par.PairData
disp(' ')
disp(' COMPARISON OF SPECIFIED PAIRS ')
disp(' ')
end
if par.PairData || par.PairFigure || par.TolHist || par.Table
for i=1:length(par.PairID)
% Load pairs to be used
[actID,sp1.id,sp2.id,goodness,theta,eta,omega]=mym(sprintf('select pairID,difAID,difBID,goodness,theta,eta,omega from %s where pairID=%d', pairtable, par.PairID(i)));
% Data of spot A
[sp1.id,sp1.maxim,sp1.extstim,sp1.extendim,sp1.area,sp1.centU,sp1.centV,sp1.bbU,sp1.bbV,sp1.bbUsize,sp1.bbVsize,sp1.integral]=...
mym(sprintf(['select difspotID,CentroidImage,ExtStartImage,ExtEndImage,Area,CentroidX,CentroidY,BoundingBoxXorigin,BoundingBoxYorigin,'...
'BoundingBoxXsize,BoundingBoxYsize,Integral from %s where difspotID=%d'],difspottable,sp1.id)) ;
% Data of spot B
[sp2.id,sp2.maxim,sp2.extstim,sp2.extendim,sp2.area,sp2.centU,sp2.centV,sp2.bbU,sp2.bbV,sp2.bbUsize,sp2.bbVsize,sp2.integral]=...
mym(sprintf(['select difspotID,CentroidImage,ExtStartImage,ExtEndImage,Area,CentroidX,CentroidY,BoundingBoxXorigin,BoundingBoxYorigin,'...
'BoundingBoxXsize,BoundingBoxYsize,Integral from %s where difspotID=%d'],difspottable,sp2.id)) ;
% Collect data of all
add_cmp.pairID=actID;
add_cmp.id1=sp1.id;
add_cmp.id2=sp2.id;
add_cmp.goodness=goodness;
add_cmp.maxim=sp2.maxim-sp1.maxim-nproj;
add_cmp.extstim=sp2.extstim-sp1.extstim-nproj;
add_cmp.extendim=sp2.extendim-sp1.extendim-nproj;
add_cmp.integral=sp2.integral/sp1.integral*100;
add_cmp.area=sp2.area/sp1.area*100;
add_cmp.bbUsize=sp2.bbUsize/sp1.bbUsize*100;
add_cmp.bbVsize=sp2.bbVsize/sp1.bbVsize*100;
cmp=[cmp, add_cmp];
% Projected bounding box of irradiated sample volume on opposite image
% 180deg offset
pr_bb = gtMatchProjectedSampleVol(sp1.centU,sp1.centV,sample_radius,sample_top,sample_bot,rot_to_det,parameters.acq.rotx) ;
% Print table with data of actual pair
if par.PairData
disp(' ')
disp(' difspotID CentroidIm ExtStartIm ExtEndIm Integral Area BBoxXsize BBoxYsize')
disp('+----------+ +----------+ +----------+ +----------+ +----------+ +----------+ +----------+ +----------+')
disp(sprintf(' %10d %10.3f %10d %10d %10d %10d %10d %10d', ...
sp1.id,sp1.maxim,sp1.extstim,sp1.extendim,sp1.integral,sp1.area,sp1.bbUsize,sp1.bbVsize)) ;
disp(sprintf(' %10d %10.3f %10d %10d %10d %10d %10d %10d', ...
sp2.id,sp2.maxim,sp2.extstim,sp2.extendim,sp2.integral,sp2.area,sp2.bbUsize,sp2.bbVsize)) ;
disp(sprintf(' Pair %5d %8.3fim %8dim %8dim %9.2f%% %9.2f%% %9.2f%% %9.2f%% ', ...
actID,add_cmp.maxim,add_cmp.extstim,add_cmp.extendim,add_cmp.integral,add_cmp.area,add_cmp.bbUsize,add_cmp.bbVsize));
end
% Draw pairfigure
if par.PairFigure
if par.Calibration
gtMatchShowPairFigure(1,parameters.acq,parameters.match_calib,sp1,sp2,pr_bb,comp_2thetas,0,[],[],[],omega,theta,eta,actID)
else
gtMatchShowPairFigure(1,parameters.acq,parameters.match,sp1,sp2,pr_bb,comp_2thetas,0,[],[],[],omega,theta,eta,actID)
end
end
if strcmpi('button',par.DoPause)
waitforbuttonpress
else
pause(par.DoPause)
end
end % of loop
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Print table of required pairs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if par.Table
disp(' ')
disp(' ')
disp(' TABLE OF ALL SPECIFIED PAIRS ')
disp(' ')
disp(' difspotIDA difspotIDB Pair ID CentroidIm ExtStartIm ExtEndIm Integral Area BBoxXsize BBoxYsize')
disp('+----------+ +----------+ +----------+ +----------+ +----------+ +----------+ +----------+ +----------+ +----------+ +----------+')
for i=1:length(cmp)
if ~isempty(cmp(i).pairID)
disp(sprintf(' %10d %10d %10d %8.3fim %8dim %8dim %9.2f%% %9.2f%% %9.2f%% %9.2f%% ', ...
cmp(i).id1,cmp(i).id2,cmp(i).pairID,cmp(i).maxim,cmp(i).extstim,cmp(i).extendim,cmp(i).integral,cmp(i).area,cmp(i).bbUsize,cmp(i).bbVsize));
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Draw histograms
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if par.TolHist
% Collect data for histogram display
for i=1:length(cmp)
hist_maxim(i)=cmp(i).maxim;
hist_extstim(i)=cmp(i).extstim;
hist_extendim(i)=cmp(i).extendim;
hist_integral(i)=cmp(i).integral;
hist_area(i)=cmp(i).area;
hist_bbUsize(i)=cmp(i).bbUsize;
hist_bbVsize(i)=cmp(i).bbVsize;
hist_goodness(i)=cmp(i).goodness;
end
% Determine bins
areabins=(floor(1/thr_area*10)/10-0.1:round((thr_area-1)*10)/100:ceil(thr_area*10)/10+0.1)*100;
bbbins=(floor(1/thr_bbsize*10)/10-0.1:round((thr_bbsize-1)*10)/100:ceil(thr_bbsize*10)/10+0.1)*100;
intbins=(floor(1/thr_intint*10)/10-0.1:round((thr_intint-1)*10)/100:ceil(thr_intint*10)/10+0.1)*100;
centbins=-ceil(thr_max_offset)-1:0.1:ceil(thr_max_offset)+1;
extbins=-ceil(thr_ext_offset)-1:ceil(thr_ext_offset)+1;
% Draw figures
figure('name','Difference in Centroid Images B-A')
hold on
xlim([min(centbins) max(centbins)])
hist(hist_maxim,centbins)
hist_maxim_max=max(hist(hist_maxim,centbins));
plot([-thr_max_offset -thr_max_offset],[0 hist_maxim_max],'r--')
plot([thr_max_offset thr_max_offset],[0 hist_maxim_max],'r--')
figure('name','Difference in Extinction Start Images; B-A')
hold on
xlim([min(extbins) max(extbins)])
hist(hist_extstim,extbins)
hist_extstim_max=max(hist(hist_extstim,extbins));
plot([-thr_ext_offset -thr_ext_offset],[0 hist_extstim_max],'r--')
plot([thr_ext_offset thr_ext_offset],[0 hist_extstim_max],'r--')
figure('name','Difference in Extinction End Images; B-A')
hold on
xlim([min(extbins) max(extbins)])
hist(hist_extendim,extbins)
hist_extendim_max=max(hist(hist_extendim,extbins));
plot([-thr_ext_offset -thr_ext_offset],[0 hist_extendim_max],'r--')
plot([thr_ext_offset thr_ext_offset],[0 hist_extendim_max],'r--')
figure('name','Ratio of Integrated Intensities; B/A %')
hold on
xlim([min(intbins) max(intbins)])
hist(hist_integral,intbins)
hist_integral_max=max(hist(hist_integral,intbins));
plot([100/thr_intint 100/thr_intint],[0 hist_integral_max],'r--')
plot([100*thr_intint 100*thr_intint],[0 hist_integral_max],'r--')
figure('name','Ratio of Areas; B/A %')
hold on
xlim([min(areabins) max(areabins)])
hist(hist_area,areabins)
hist_area_max=max(hist(hist_area,areabins));
plot([100/thr_area 100/thr_area],[0 hist_area_max],'r--')
plot([100*thr_area 100*thr_area],[0 hist_area_max],'r--')
figure('name','Ratio of BBox X size; B/A %')
hold on
xlim([min(bbbins) max(bbbins)])
hist(hist_bbUsize,bbbins)
hist_bbUsize_max=max(hist(hist_bbUsize,bbbins));
plot([100/thr_bbsize 100/thr_bbsize],[0 hist_bbUsize_max],'r--')
plot([100*thr_bbsize 100*thr_bbsize],[0 hist_bbUsize_max],'r--')
figure('name','Ratio of BBox Y size; B/A %')
hold on
xlim([min(bbbins) max(bbbins)])
hist(hist_bbVsize,bbbins)
hist_bbVsize_max=max(hist(hist_bbVsize,bbbins));
plot([100/thr_bbsize 100/thr_bbsize],[0 hist_bbVsize_max],'r--')
plot([100*thr_bbsize 100*thr_bbsize],[0 hist_bbVsize_max],'r--')
figure('name','Histogram of goodness of match')
hold on
hist(hist_goodness,max(length(hist_goodness)/20,10))
hist_goodness_max=max(hist(hist_goodness,max(length(hist_goodness)/20,10)));
if thr_goodness ~= Inf
plot([thr_goodness thr_goodness],[0 hist_goodness_max],'r--')
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Draw other figures
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if par.PairDist
% Bin width for matching percentage
binwidth=5;
if par.Calibration
% Use matching constraint for selection
if ~isempty(parameters.match_calib.addconstr)
addconstr=['and ' parameters.match_calib.addconstr];
else
addconstr=[];
end
else
% Use matching constraint for selection
if ~isempty(parameters.match.addconstr)
addconstr=['and ' parameters.match.addconstr];
else
addconstr=[];
end
end
% All centroid image (~omega) values
[centimA]=mym(sprintf('SELECT CentroidImage FROM %s WHERE CentroidImage<=%d AND BoundingBoxXsize>%0.20g AND BoundingBoxYsize>%0.20g %s',difspottable,nproj,minsizeX,minsizeY,addconstr))*180/nproj;
[centimB]=mym(sprintf('SELECT CentroidImage FROM %s WHERE CentroidImage>%d AND BoundingBoxXsize>%0.20g AND BoundingBoxYsize>%0.20g %s',difspottable,nproj,minsizeX,minsizeY,addconstr))*180/nproj;
% Distribution in the histograms
edgesA=0:binwidth:180;
edgesB=180:binwidth:360;
centimA=histc(centimA,edgesA);
centimB=histc(centimB,edgesB);
centom=min(centimA,centimB);
% All angles
[theta,eta,omega]=mym(sprintf('SELECT theta,eta,omega from %s',pairtable));
match_omega=histc(omega,edgesA)';
if size(match_omega,1)==1
match_omega=match_omega';
end
edges_eta=0:binwidth:360;
match_eta=histc(eta,edges_eta);
matchpc=match_omega./centom*100;
% Draw figures
figure('name','Theta vs omega')
hold on
if isfield('parameters','match_calib') %,'theta_bounds')
thbs=parameters.match_calib.theta_bounds;
if ~isempty(thbs)
for i=1:size(thbs,1)
if thbs(i,2) > thbs(i,1)
rectangle('Position',[0,thbs(i,1),180,thbs(i,2)-thbs(i,1)],...
'FaceColor','y','EdgeColor','none')
end
end
end
end
plot(omega,theta,'b.')
for i=1:length(comp_2thetas)
plot([0 380],[comp_2thetas(i) comp_2thetas(i)]/2,'g-.')
end
axis([0 180 min(comp_2thetas)/2-1 max(comp_2thetas)/2+1])
figure('name','Eta vs theta')
hold on
if isfield('parameters','match_calib') %,'theta_bounds')
thbs=parameters.match_calib.theta_bounds;
if ~isempty(thbs)
for i=1:size(thbs,1)
if thbs(i,2) > thbs(i,1)
rectangle('Position',[thbs(i,1),0,thbs(i,2)-thbs(i,1),360],...
'FaceColor','y','EdgeColor','none')
end
end
end
end
plot(theta,eta,'b.')
for i=1:length(comp_2thetas)
plot([comp_2thetas(i) comp_2thetas(i)]/2,[0 380],'g-.')
bndlow = (comp_2thetas(i)-thr_ang-thr_ang_scale*comp_2thetas(i))/2;
bndhigh = (comp_2thetas(i)+thr_ang+thr_ang_scale*comp_2thetas(i))/2;
plot([bndlow bndlow],[0 380],'r--')
plot([bndhigh bndhigh],[0 380],'r--')
end
axis([min(comp_2thetas)/2-1 max(comp_2thetas)/2+1 0 380])
figure('name','No. of pairs vs theta')
edges_theta=round((min(comp_2thetas)/2-0.1)*100)/100:0.002:max(comp_2thetas)/2+0.1;
theta_h=histc(theta,edges_theta);
bar(edges_theta,theta_h)
h=findobj(gca,'Type','patch');
set(h,'EdgeColor','none')
hold on
for i=1:length(comp_2thetas)
plot([comp_2thetas(i) comp_2thetas(i)]/2,[0 max(theta_h)],'g-.')
bndlow = (comp_2thetas(i)-thr_ang-thr_ang_scale*comp_2thetas(i))/2;
bndhigh = (comp_2thetas(i)+thr_ang+thr_ang_scale*comp_2thetas(i))/2;
plot([bndlow bndlow],[0 max(theta_h)],'r--')
plot([bndhigh bndhigh],[0 max(theta_h)],'r--')
end
xlim([min(edges_theta) max(edges_theta)])
figure('name','No. of pairs found vs eta')
bar((edges_eta+binwidth/2)',match_eta)
h=findobj(gca,'Type','patch');
set(h,'EdgeColor','none')
xlim([0,360])
figure('name','Matching percentage vs omega')
bar((edgesA+binwidth/2)',matchpc)
xlim([0,180])
h=findobj(gca,'Type','patch');
set(h,'EdgeColor','none')
set(h,'FaceColor','r')
figure('name','No. of spots and no. of pairs vs omega')
bar((edgesA+binwidth/2)',[centom,match_omega])
h=findobj(gca,'Type','patch');
set(h,'EdgeColor','none')
xlim([0,180])
end
end % of function
3_pairmatching/gtMatchDefaultParameters.m deleted 100644 → 0
+ 0
− 31
View file @ f63ce796
%
% FUNCTION par_match=gtMatchDefaultParameters
%
% Stores and provides the default parameters for difspot matching.
%
%
function par_match=gtMatchDefaultParameters()
par_match.thr_ang = 0.2;
par_match.thr_ang_scale = 0;
par_match.thr_max_offset = 2;
par_match.thr_ext_offset = 10;
par_match.thr_genim_offset = 2;
par_match.corr_rot_to_det = 0;
par_match.thr_for_corr = 0.2;
par_match.thr_goodness = Inf;
par_match.thr_intint = 6;
par_match.thr_area = 1.6;
par_match.thr_bbsize = 1.5;
par_match.minsizeX = 5;
par_match.minsizeY = 5;
par_match.tiltY = 0;
par_match.tiltZ = 0;
par_match.centmode = 0;
par_match.addconstr = [];
end
\ No newline at end of file
3_pairmatching/gtMatchDefaultParametersCalib.m deleted 100644 → 0
+ 0
− 33
View file @ f63ce796
%
% FUNCTION par_match=gtMatchDefaultParametersCalib
%
% Stores and provides the default calibration parameters for difspot
% matching.
%
function par_match=gtMatchDefaultParametersCalib()
par_match.thr_ang = 45;
par_match.thr_ang_scale = 0;
par_match.thr_max_offset = 2;
par_match.thr_ext_offset = 5;
par_match.thr_genim_offset = 2;
par_match.corr_rot_to_det = 0;
par_match.thr_for_corr = 0;
par_match.thr_goodness = 0.2;
par_match.thr_intint = 3;
par_match.thr_area = 1.2;
par_match.thr_bbsize = 1.2;
par_match.minsizeX = 10;
par_match.minsizeY = 10;
par_match.tiltY = 0;
par_match.tiltZ = 0;
par_match.centmode = 0;
par_match.addconstr = [];
par_match.theta_bounds = [];
end
\ No newline at end of file
3_pairmatching/gtMatchDiffVecInLab.m deleted 100644 → 0
+ 0
− 38
View file @ f63ce796
%
% FUNCTION dv=gtMatchDiffVecInLab(scUA,scVA,scUB,scVB,rottodet,tiltY,tiltZ,rotU,sorZ)
%
% Given the Scintillator coordinates of Friedel pairs of diffraction spots
% A and B, computes the LAB coordinates of the diffraction vector they
% correspond to, according to tilts and rotation axis to detector
% (scintillator) distance.
%
% INPUT scUA,scVA - coloumn vectors of Scintillator coordinates of point A (in pixels or mm-s)
% scUB,scVB - coloumn vectors of Scintillator coordinates of point B (in pixels or mm-s)
% rottodet - rotation axis to detector (scintillator) distance at
% rotU,sorZ (same units as sc-s)
% tiltY - tilt around axis Y in Lab system (in degrees)
% tiltZ - tilt around axis Z in Lab system (in degrees)
% rotU - horizontal location of the rotation axis in the images
% sorZ - origin Z of Sample system
%
% OUTPUT dv - coloumn vector of diffraction vector LAB coordinates
%
function dv=gtMatchDiffVecInLab(scUA,scVA,scUB,scVB,rottodet,tiltY,tiltZ,rotU,sorZ)
[labXA,labYA,labZA] = gtMatchSc2Lab(scUA,scVA,rottodet,tiltY,tiltZ,rotU,sorZ);
[labXB,labYB,labZB] = gtMatchSc2Lab(scUB,scVB,rottodet,tiltY,tiltZ,rotU,sorZ);
% Flip X and Y coordinates of B (to get the real position of spot B, mirror
% on the rotation axis while staying in the same LABA reference). Diffraction
% vector is parallel to the path between A and B mirrored:
dv = [labXA,labYA,labZA] - [-labXB,-labYB,labZB];
% Normalize:
dnorm = sqrt(dv(:,1).^2+dv(:,2).^2+dv(:,3).^2) ;
dv(:,1) = dv(:,1)./dnorm ;
dv(:,2) = dv(:,2)./dnorm ;
dv(:,3) = dv(:,3)./dnorm ;
end
3_pairmatching/gtMatchEtaOfPoint.m deleted 100644 → 0
+ 0
− 56
View file @ f63ce796
%
% FUNCTION etas=gtMatchEtaOfPoint(a,b)
%
% Determines the corresponding eta angles of points on the detector or from
% diffraction vectors.
%
% INPUT
% Number of input arguments:
% 1 -> vectors (n,3) of LAB coordinates X,Y,Z of diffraction vectors
% 2 -> LAB coordinates Y and Z of points
%
% OUTPUT etas(n,1) - coloumn vector of eta angles
%
function etas=gtMatchEtaOfPoint(varargin)
if nargin == 1 % takes a diffraction vector
inp = varargin{1};
y= inp(:,2);
z= inp(:,3);
elseif nargin == 2 % takes two coordinates on detector
y= varargin{1};
z= varargin{2};
else
error('Number of input arguments must be either 1 or 2!')
end
etas=zeros(length(y),1);
for i=1:length(y)
if y(i) > 0 % 0 < eta < 180deg
if z(i) > 0 % 0 < eta < 90deg
etas(i)=atand(y(i)/z(i)) ;
elseif z(i) == 0
etas(i)=90;
else % 90deg < eta < 180deg
etas(i)=atand(y(i)/z(i))+180 ;
end
elseif y(i) == 0
if z(i) > 0
etas(i)=0 ;
elseif z(i) == 0
etas(i)=0;
else
etas(i)=180 ;
end
else % 180deg < eta < 360deg
if z(i) < 0 % 180deg < eta < 270deg
etas(i)=atand(y(i)/z(i))+180 ;
elseif z(i) == 0
etas(i)=270;
else % 270deg < eta < 360deg
etas(i)=atand(y(i)/z(i))+360 ;
end
end
end
3_pairmatching/gtMatchLab2Sam.m deleted 100644 → 0
+ 0
− 32
View file @ f63ce796
%
% FUNCTION [samX,samY,samZ]=gtMatchLab2Sam(labX,labY,labZ,om)
%
% LAB -> SAMPLE coordinate transformation (omega rotation)
%
% Given the LAB coordinates of a point, transforms them into Sample
% coordinates according to omega (where omega is normally positive between
% 0 and 360 degrees)
%
% LAB system is right-handed and defined as:
% X - beam direction, origin on the rotation axis, positive towards the camera
% Y - lateral, origin on the rotation axis
% Z - vertical, origin defined in the parameters file (usually at the centre of
% the image), positive upwards
%
% SAMPLE system is fixed to the sample. Coordinates are defined as LAB
% coordinates at omega=0.
%
% INPUT labX,labY,labZ - coloumn vectors of LAB coordinates (arb. units)
% om - coloumn vector of omega angles (in degrees)
%
% OUTPUT samX,samY,samZ - coloumn vectors of SAMPLE coordinates
%
function [samX,samY,samZ]=gtMatchLab2Sam(labX,labY,labZ,om)
samX= cosd(om).*labX - sind(om).*labY ;
samY= sind(om).*labX + cosd(om).*labY ;
samZ= labZ ;
end
3_pairmatching/gtMatchPlaneNormInLab.m deleted 100644 → 0
+ 0
− 37
View file @ f63ce796
%
% FUNCTION plnorm=gtMatchPlaneNormInLab(scXA,scYA,scXB,scYB,rottodet,tiltY,tiltZ,rotx,sorZ)
%
% Given the Scintillator coordinates of pairs of diffraction spots A and B,
% computes the LAB coordinates of the diffracting plane normals they correspond
% to, according to tilts and rotation axis to detector (scintillator) distance.
%
% INPUT scXA = coloumn vector of X coordinates on the Scintillator of point A (pixels or mm-s)
% scYA = coloumn vector of Y coordinates on the Scintillator of point A (pixels or mm-s)
% scXB = coloumn vector of X coordinates on the Scintillator of point B (pixels or mm-s)
% scYB = coloumn vector of Y coordinates on the Scintillator of point B (pixels or mm-s)
% rottodet = rotation axis to detector (scintillator) distance at rotx,sorZ (pixels or mm-s)
% /scX,scY,rottodet have to be in the same units; e.g. pixels or mm-s/
% tiltY = tilt around axis Y in Lab system (in degrees)
% tiltZ = tilt around axis Z in Lab system (in degrees)
% rotx = location of rotation axis in the images (X coord. in pixels)
% sorZ = Z=0 location in the images (Y coord. in pixels at the rotation axis)
%
% OUTPUT plnorm = coloumn vector of plane normals of size(plnorm)=[n,3]
%
%
function plnorm=gtMatchPlaneNormInLab(dv)
% plane normal:
plv = dv-repmat([1 0 0],size(dv,1),1); % plane vector
%plnorm = plv/norm(plv);
plvnorm=sqrt(plv(:,1).^2+plv(:,2).^2+plv(:,3).^2) ; % norm of plane vectors
plnorm(:,1)= plv(:,1)./plvnorm ; % normalized plane vectors
plnorm(:,2)= plv(:,2)./plvnorm ;
plnorm(:,3)= plv(:,3)./plvnorm ;
end
3_pairmatching/gtMatchProjectedSampleVol.m deleted 100644 → 0
+ 0
− 112
View file @ f63ce796
%
% FUNCTION bb = gtMatchProjectedSampleVol(spotU, spotV, sample_radius,
% sample_top, sample_bot, rot_to_det, rotU)
%
% Given the location of a dffraction spotA, computes the bounding box of
% the projection of a cylindrical sample as visible on the detector plane
% exactly 180 degrees offset.
%
% INPUT:
%
% Input is given in IMAGE coordinates in pixels.
% Origin: spotU=0 anywhere, e.g. at the left edge of the image; positive to the right in the image
% spotV=0 anywhere, e.g. at the top of the image, positive downwards
%
%
% spotU,spotV - spotA centroid U coordinate in its own image
% sample_radius - radius of the cylindrical sample
% sample_top - V coordinate in pixels of the top of the illuminated part
% sample_bot - V coordinate in pixels of the bottom of the illuminated part
% (sample_top < sample_bot)
% rot_to_det - rotation axis to detector distance in pixels
% rotU - position U of the rotation axis in the image in pixels
%
% OUTPUT:
%
% Bounding-box of projection, upper left and lower right corners in the offset image
%
% bb = [up_leftU; up_leftV; bo_rightU; bo_rightV]
% i.e. [left_edgeU; top_edgeV; right_edgeU; bottom_edgeV]
%
%
% Peter Reischig, ESRF, 11/2006
%
% /Verified/
%%
%don't understand why, but seems to not give enough vertical range for a
%wide but flat sample volume.
function bb=gtMatchProjectedSampleVol(spotU, spotV, sample_radius, sample_top, sample_bot, rot_to_det, rotU)
spotU=spotU-rotU ;
a=sqrt(rot_to_det^2+spotU^2) ;
beta=asind(sample_radius/a) ;
if spotU >=0
alpha=asind(rot_to_det/a) ;
else
alpha=180-asind(rot_to_det/a) ;
end
% upper-left U (left edge)
bb(1,1)=rotU - spotU+2*rot_to_det*cotd(alpha+beta) ;
% upper-left V (top edge, lowest V)
bb(2,1)=spotV+(sample_top-spotV)*(2*rot_to_det/(rot_to_det+sign(sample_top-spotV)*sample_radius)) ;
% bottom right U (right edge)
bb(3,1)=rotU - spotU+2*rot_to_det*cotd(alpha-beta) ;
% bottom-right V (bottom edge, highest V)
bb(4,1)=spotV+(sample_bot-spotV)*(2*rot_to_det/(rot_to_det-sign(sample_bot-spotV)*sample_radius)) ;
end
% keyboard
%
% %define sample "corners"
% corners(1,:) = [sample_radius sample_bot rot_to_det-sample_radius];
% corners(2,:) = [-sample_radius sample_bot rot_to_det-sample_radius];
% corners(3,:) = [sample_radius sample_top rot_to_det-sample_radius];
% corners(4,:) = [-sample_radius sample_top rot_to_det-sample_radius];
% corners(5,:) = [sample_radius sample_bot rot_to_det+sample_radius];
% corners(6,:) = [-sample_radius sample_bot rot_to_det+sample_radius];
% corners(7,:) = [sample_radius sample_top rot_to_det+sample_radius];
% corners(8,:) = [-sample_radius sample_top rot_to_det+sample_radius];
%
% %done above
% %spotU=spotU-rotU;
%
% %relative distance from spot to corner
% deltaX=corners(:,1)-spotU;
% deltaY=corners(:,2)-spotV;
% %relative distance from spot to projection of corner
% deltaX=deltaX.*(2*rot_to_det./corners(:,3));
% deltaY=deltaY.*(2*rot_to_det./corners(:,3));
% %take extremes
% deltaXmin=min(deltaX);
% deltaXmax=max(deltaX);
% deltaYmin=min(deltaY);
% deltaYmax=max(deltaY);
% %assemble output
% % leftX (min X)
% bb(1)=spotU+deltaXmin;
% %allow for rotated detector
% bb(1)=rotU-bb(1);
% % upY (min Y)
% bb(2)=spotV+deltaYmin;
% % rightX (max X)
% bb(3)=spotU+deltaXmax;
% %allow for rotated detector
% bb(3)=rotU-bb(3);
% % botY (max Y)
% bb(4)=spotV+deltaYmax;
3_pairmatching/gtMatchSc2Lab.m deleted 100644 → 0
+ 0
− 48
View file @ f63ce796
%
% FUNCTION [labX,labY,labZ]=gtMatchSc2Lab(scU,scV,rottodet,tiltY,tiltZ,rotU,sorZ))
%
% SCINTILLATOR -> LAB coordinate transformation
%
% Given the Scintillator coordinates of a point, transforms them into LAB
% coordinates according to tilts and rotation axis to detector (scintillator)
% distance.
%
% SCINTILLATOR coordinates are defined as the image coordinates (after
% undistortion correction). Images are seen looking into the beam.
% U = horizontal, positive to the right, origin is at the rotation axis location
% (at pixel rotU)
% V = downwards, origin is at pixel sorZ (origin Z of SAMPLE system)
%
% LAB system is right-handed and defined as:
% X - beam direction, origin on the rotation axis, positive towards the camera
% Y - lateral, origin on the rotation axis
% Z - vertical, origin defined in the parameters file (usually at the centre of
% the image), positive upwards
%
% INPUT scU,scV - coloumn vectors of Scintillator coordinates (in pixels or mm-s)
% rottodet - rotation axis to detector (scintillator) distance at
% rotU,sorZ (same units as scU,scV)
% tiltY - tilt around axis Y in Lab system (in degrees)
% tiltZ - tilt around axis Z in Lab system (in degrees)
% rotU - horizontal location of the rotation axis in the images
% sorZ - origin Z of Sample system
%
% OUTPUT labX,labY,labZ - coloumn vectors the LAB coordinates
%
function [labX,labY,labZ]=gtMatchSc2Lab(scU,scV,rottodet,tiltY,tiltZ,rotU,sorZ)
% Transformation matrix:
TR=[-sind(tiltZ)*cosd(tiltY) -sind(tiltY);...
cosd(tiltZ) 0 ;...
sind(tiltZ)*sind(tiltY) -cosd(tiltY)];
% Vector in LAB coordinates:
labvec = TR*[scU'-rotU; scV'-sorZ];
labX = labvec(1,:)' + rottodet;
labY = labvec(2,:)';
labZ = labvec(3,:)';
end
3_pairmatching/gtMatchShowPairFigure.m deleted 100644 → 0
+ 0
− 244
View file @ f63ce796
%
% FUNCTION gtMatchShowPairFigure(fig,par_acq,par_match,spA,spB,pr_bb,...
% comp_2thetas,theta_band,spC,fullimnoA,fullimnoB,omega,theta,eta,pairID)
%
% Plots a figure of matched or unmatched Friedel pairs.
% Used by gtMatchDifspots, gtShowDifspots, gtShowPairs, etc.
%
% INPUT
% fig: figure number for display (e.g. gcf); if 0, only fullimage is
% updated
% parameters: parameters loaded from parameters file
% spA,spB: data of spot A and its pair B (e.g. from database) with the
% fields:
% spA.id: difspotID
% spA.maxim: # of maximum or centroid image
% spA.bbX,spA.bbY,spA.bbXsize,spA.bbYsize: bounding box
% spA.centX,spA.centY: centroid of spot
% if spB.id is empty, spot B is not displayed
% OPTIONAL INPUT
% pr_bb: projected sample bounding box; if 0, it's calculated from
% parameters; not displayed if empty or not specified
% comp_2thetas: vector of computed 2theta-s
% theta_band: if true, tolerance bands are displayed for first pair
% spC: any other spot candidates with a structure like spB
% fullimnoA,fullimnoB: # of full images to be displayed; if not
% specified, spA and spB.maxim (or spA.maxim+nproj)
% will be displayed
% omega,theta,eta: to be displayed above the full images
% pairID displayed in figure if specified
%
%
%
% by Peter Reischig ESRF,10/2007
%
function gtMatchShowPairFigure(fig,par_acq,par_match,spA,spB,pr_bb,...
comp_2thetas,theta_band,spC,fullimnoA,fullimnoB,omega,theta,eta,pairID)
parameters.acq=par_acq;
parameters.match=par_match;
if isfield(spA,'bbU')
spA.bbX=spA.bbU;
spA.bbY=spA.bbV;
spA.bbXsize=spA.bbUsize;
spA.bbYsize=spA.bbVsize;
spA.centX=spA.centU;
spA.centY=spA.centV;
spB.bbX=spB.bbU;
spB.bbY=spB.bbV;
spB.bbXsize=spB.bbUsize;
spB.bbYsize=spB.bbVsize;
spB.centX=spB.centU;
spB.centY=spB.centV;
end
text_dist=20;
rot_to_det=parameters.acq.dist/parameters.acq.pixelsize+parameters.match.corr_rot_to_det;
% Projected sample bb
if ~exist('pr_bb','var')
pr_bb=[];
end
if pr_bb==0
sample_radius=parameters.acq.bb(3)/2;
sample_top=parameters.acq.bb(2);
sample_bot=parameters.acq.bb(2)+parameters.acq.bb(4);
pr_bb=projected_sample(spA.centX,spA.centY,sample_radius,sample_top,sample_bot,rot_to_det,parameters.acq.rotx);
end
% 2thetas and candidate spots
if ~exist('comp_2thetas','var')
comp_2thetas=[];
end
if ~exist('theta_band','var')
theta_band=false;
end
if ~exist('spC','var')
spC=[];
end
if ~isfield(spC,'id')
spC.id=[];
end
% Full images to display
if ~exist('fullimnoA','var')
fullimnoA=[];
end
if isempty(fullimnoA)
fullimnoA=spA.maxim;
end
if ~exist('fullimnoB','var')
fullimnoB=[];
end
if isempty(fullimnoB)
fullimnoB=spB.maxim;
end
if isempty(fullimnoB)
fullimnoB=mod(fullimnoA+parameters.acq.nproj,2*parameters.acq.nproj);
end
fullimnoA=round(fullimnoA);
fullimnoB=round(fullimnoB);
% Angles
if ~exist('theta','var')
theta=[];
end
if ~exist('eta','var')
eta=[];
end
if ~exist('omega','var')
omega=[];
end
% PairID
if ~exist('pairID','var')
pairID=[];
end
fullimages_directory=sprintf('%s/1_preprocessing/full',parameters.acq.dir) ;
fullimB=edf_read([fullimages_directory, sprintf('/full%0.4d.edf',fullimnoB)]);
spotimA=edf_read([fullimages_directory, sprintf('/full%0.4d.edf',fullimnoA)],...
[spA.bbX spA.bbY spA.bbXsize spA.bbYsize]);
%spotimA=gtCrop(fullimA,[spA.bbX spA.bbY spA.bbXsize spA.bbYsize]);
%directim=gtCrop(fullimA,parameters.acq.bb);
clims=[-max(spotimA(:)) max(spotimA(:))];
% Spot images
if fig==0
figure(gcf) % only updates full image (search pair of same spot )
else
figure(fig) % updates the whole figure with new spot
clf(fig)
set(gcf,'name','Difspot Pairs')
set(gcf,'Units','normalized','Position',[0 0 1 1]);
subp_totheight=max([0.2 1/round(length(spC.id)/2+1)]);
subp_height=subp_totheight-0.02;
subp_width=0.19;
% Image of specified spot
subplot('position',[0 1-subp_totheight subp_width subp_height])
difspim=gtGetSummedDifSpot(spA.id,parameters,'connected',1);
imshow(difspim,clims)
hold on
rectangle('Position',[0.5,0.5,spA.bbXsize,spA.bbYsize],'EdgeColor','r');
plot(spA.centX-spA.bbX+1,spA.centY-spA.bbY+1,'r+');
title(sprintf('DifspotA: %d', spA.id));
axis([0,spA.bbXsize+1,0,spA.bbYsize+1])
% Image of pair if exists
if ~isempty(spB.id)
subplot('position',[subp_width 1-subp_totheight subp_width subp_height])
difspim=gtGetSummedDifSpot(spB.id,parameters,'connected',1);
imshow(fliplr(difspim),clims)
hold on
rectangle('Position',[0.5,0.5,spB.bbXsize,spB.bbYsize],'EdgeColor','g');
plot(2*(spB.bbXsize/2+0.5)-(spB.centX-spB.bbX+1),spB.centY-spB.bbY+1,'g+'); % spot id flipped
title(sprintf('DifspotB: %d', spB.id));
axis([0,spB.bbXsize+1,0,spB.bbYsize+1])
end
% Images of candidate spots
for spi=1:length(spC.id)
subplot('position',[mod(spi+1,2)*subp_width floor((-spi-2)/2)*subp_totheight+1 subp_width subp_height])
difspim=gtGetSummedDifSpot(spC.id(spi),parameters,'connected',1);
imshow(fliplr(difspim),clims)
hold on
plot(2*(spC.bbXsize(spi)/2+0.5)-(spC.centX(spi)-spC.bbX(spi)+1),spC.centY(spi)-spC.bbY(spi)+1,'b+'); % spot flipped
title(sprintf('DifspotB: %d',spC.id(spi)));
axis([0,spC.bbXsize(spi),0,spC.bbYsize(spi)])
end
end
% Display full images
subplot('Position',[0.4 0 0.6 1])
%totim=gtPlaceSubImage(directim,fliplr(fullimB),parameters.acq.bb(1),parameters.acq.bb(2));
totim=fliplr(fullimB);
totim=gtPlaceSubImage(spotimA,totim,spA.bbX,spA.bbY);
imshow(totim,clims);
hold on
title(sprintf(['PairID: %d DifspotA: %d DifspotB: %d #FullimA: %d #FullimB: %d Omega: %5.2fdeg'...
' Theta: %5.2fdeg Eta: %5.2fdeg \n\n (shown is real position of DifspotA and FullimB flipped)'],...
pairID,spA.id,spB.id,fullimnoA,fullimnoB,omega,theta,eta)) ;
rectangle('Position',[parameters.acq.bb(1)-2,parameters.acq.bb(2)-2,parameters.acq.bb(3)+4,parameters.acq.bb(4)+4],'EdgeColor','c');
rectangle('Position',[spA.bbX-2,spA.bbY-2,spA.bbXsize+4,spA.bbYsize+4],'EdgeColor','r') ;
text(spA.bbX+spA.bbXsize+text_dist,spA.bbY+spA.bbYsize+text_dist,int2str(spA.id),'Color','r');
if ~isempty(pr_bb)
rectangle('Position',[2*parameters.acq.rotx-pr_bb(1)-(pr_bb(3)-pr_bb(1)),pr_bb(2),(pr_bb(3)-pr_bb(1)),(pr_bb(4)-pr_bb(2))],'EdgeColor','b');
end
% Display pair if exists
if ~isempty(spB.id)
rectangle('Position',[2*parameters.acq.rotx-spB.bbX-spB.bbXsize-2,spB.bbY-2,spB.bbXsize+4,spB.bbYsize+4],'EdgeColor','g') ;
plot([spA.centX, 2*parameters.acq.rotx-spB.centX],[spA.centY, spB.centY],'c') ;
text(2*parameters.acq.rotx-spB.bbX,spB.bbY+spB.bbYsize+text_dist,int2str(spB.id),'color','g');
end
% Candidate spot centroids in full image
for k=1:length(spC.id)
plot(2*parameters.acq.rotx-spC.centX(k),spC.centY(k),'g+') ;
text(2*parameters.acq.rotx-spC.centX(k)+text_dist,spC.centY(k)+text_dist,int2str(spC.id(k)),'color','g');
end
%Circles along 2theta cones
t=(0:1:360)' ;
if ~isempty(comp_2thetas)
circRthetas=2*rot_to_det*tand(comp_2thetas) ;
for k=1:length(circRthetas)
circthetasX=spA.centX+cosd(t)*circRthetas(k) ;
circthetasY=spA.centY+sind(t)*circRthetas(k) ;
plot(circthetasX,circthetasY,'b') ;
end
if theta_band
circRthetas=2*rot_to_det*tand(comp_2thetas+parameters.match.thr_ang) ;
for k=1:length(circRthetas)
circthetasX=spA.centX+cosd(t)*circRthetas(k) ;
circthetasY=spA.centY+sind(t)*circRthetas(k) ;
plot(circthetasX,circthetasY,'y') ;
end
circRthetas=2*rot_to_det*tand(comp_2thetas-parameters.match.thr_ang) ;
for k=1:length(circRthetas)
circthetasX=spA.centX+cosd(t)*circRthetas(k) ;
circthetasY=spA.centY+sind(t)*circRthetas(k) ;
plot(circthetasX,circthetasY,'y') ;
end
end
end
drawnow
end % of function
3_pairmatching/gtMatchThetaOfPair.m deleted 100644 → 0
+ 0
− 16
View file @ f63ce796
%
% FUNCTION thetas=gtMatchThetaOfPair(dv)
%
% Determines the corresponding Bragg angles from diffraction vector LAB
% coordinates.
%
% INPUT dv(n,3) - X,Y,Z coordinates of the diffraction vectors
%
% OUTPUT thetas(n,1) - coloumn vector of theta angles
%
function thetas=gtMatchThetaOfPair(dv)
thetas=0.5*acosd(dv(:,1)./sqrt(dv(:,1).^2+dv(:,2).^2+dv(:,3).^2));
end
3_pairmatching/gtSetupPairMatching.m deleted 100644 → 0
+ 0
− 247
View file @ f63ce796
function gtSetupPairMatching()
% GTSETUPPAIRMATCHING Helper function to wrap around Peter's matching functions
% gtSetupPairMatching
% -------------------
% Modified on January 21, 2012 by Lnervos
parameters=[];
load('parameters.mat');
% close figures before starting
close all
% make a small figure, top left of screen, with the required buttons
dim.ScreenSize=get(0, 'ScreenSize');
dim.FigureSize=[440 300];
dim.Margin=25;
% now can create the figure
fh=2; % figure 1 is used by display of match difspots
figure(fh);
set(fh, 'Visible', 'on', 'Position', [dim.Margin dim.ScreenSize(4)-dim.FigureSize(2)-dim.Margin dim.FigureSize(1) dim.FigureSize(2)],...
'menubar', 'none', 'name', 'Control pair matching');
% building the GUI
% HEADER TEXT AT THE TOP
uicontrol(fh, 'Style', 'text', 'String', {'Use these buttons to run pair matching', 'Look at the desktop / terminal for output'}, ...
'Position', [1 dim.FigureSize(2)-70 dim.FigureSize(1) 60]);
% CALIBRATION COLUMN (L)
uicontrol(fh, 'Style', 'text', 'String', 'Calibration - match good spots only to refine geometry', ...
'Position', [10 dim.FigureSize(2)-130 200 50]);
RunCalibButton=uicontrol(fh, 'Style', 'pushbutton', 'String', 'Do Calibration Matching', ...
'Position', [10 dim.FigureSize(2)-170 200 30], 'callback', @sfRunCalibButtonFunction);
ModCalibParamsButton=uicontrol(fh, 'Style', 'pushbutton', 'String', 'Change Calib Parameters', ...
'Position', [10 dim.FigureSize(2)-210 200 30], 'callback', @sfModCalibParamsButtonFunction);
CorGeoCalibButton=uicontrol(fh, 'Style', 'pushbutton', 'String', 'Refine geometry from Calib', ...
'Position', [10 dim.FigureSize(2)-250 200 30], 'callback', @sfCorGeoCalibButtonFunction);
% REGULAR COLUMN (R)
uicontrol(fh, 'Style', 'text', 'String', {'Matching - match all spots', 'Set MaxTime to Inf to run whole dataset'}, ...
'Position', [220 dim.FigureSize(2)-130 210 50]);
RunNormalButton=uicontrol(fh, 'Style', 'pushbutton', 'String', 'Do Matching', ...
'Position', [220 dim.FigureSize(2)-170 210 30], 'callback', @sfRunNormalButtonFunction);
ModNormalParamsButton=uicontrol(fh, 'Style', 'pushbutton', 'String', 'Change Matching Parameters', ...
'Position', [220 dim.FigureSize(2)-210 210 30], 'callback', @sfModNormalParamsButtonFunction);
CorGeoNormalButton=uicontrol(fh, 'Style', 'pushbutton', 'String', 'Refine geometry from Matching', ...
'Position', [220 dim.FigureSize(2)-250 210 30], 'callback', @sfCorGeoNormalButtonFunction);
% QUIT BUTTON
QuitButton=uicontrol(fh, 'Style', 'pushbutton', 'String', 'Finish', ...
'Position', [(dim.FigureSize(1)/2)-50 dim.FigureSize(2)-290 100 30], 'callback', {@sfQuitButtonFunction, fh});
% wait until finished - thanks Greg
quit=0;
while ~quit
drawnow
end
% quit function
close(fh)
quit=0;
end % end of function
%%
% SUB-sFUNCTIONS
function sfRunCalibButtonFunction(varargin)
% default parameters
par.ShowFindings = false;
par.DoPause = 0;
par.Calibration = true;
par.MaxTime = Inf;
% list/captions
list{1,1}='ShowFindings'; list{1,2}={'Display results for each pair'};
list{2,1}='DoPause'; list{2,2}={'Time to pause on pair, or "button" for mouse click'};
par=gtModifyStructure(par, list);
gtMATCHDifspots('ShowFindings', par.ShowFindings, 'DoPause', par.DoPause, 'Calibration', par.Calibration)
end
function sfModCalibParamsButtonFunction(varargin)
parameters=[];
load('parameters.mat');
list{1,1} = 'thr_ang'; list{1,2} = 'Angular threshold when searching';
list{2,1} = 'thr_ang_scale'; list{2,2} = '???';
list{3,1} = 'thr_max_offset'; list{3,2} = 'Centroid image offset threshold (#images)';
list{4,1} = 'thr_ext_offset'; list{4,2} = 'ExtStart / ExtEnd offset threshold (#images)';
list{5,1} = 'thr_genim_offset'; list{5,2} = 'Threshold for at least one of Centroid/ExtStart/ExtEnd';
list{6,1} = 'corr_rot_to_det'; list{6,2} = 'Correction to rot to det distance (pixels)';
list{7,1} = 'thr_for_corr'; list{7,2} = 'Threshold goodness for using image correction';
list{8,1} = 'thr_goodness'; list{8,2} = 'Acceptable limit: 0:perfect, <0.3:good, >0.7:bad';
list{9,1} = 'thr_intint'; list{9,2} = 'Integrated intensity threshold ratio';
list{10,1} = 'thr_area'; list{10,2} = 'Bounding box area ratio';
list{11,1} = 'thr_bbsize'; list{11,2} = 'Bounding box size ratio';
list{12,1} = 'minsizeX'; list{12,2} = 'Minimum bounding box X size (pixels)';
list{13,1} = 'minsizeY'; list{13,2} = 'Minimum bounding box Y size (pixels)';
list{14,1} = 'tiltY'; list{14,2} = 'Detector tilt around Y-axis (degrees)' ;
list{15,1} = 'tiltZ'; list{15,2} = 'Detector tilt around Z-axis (degrees)';
list{16,1} = 'centmode'; list{16,2} = '0:cent of mass 1:spot correl 2:im correl';
parameters.match_calib=gtModifyStructure(parameters.match_calib, list);
save parameters parameters
end
function sfCorGeoCalibButtonFunction(varargin)
% default parameters
par.OptPar = 'all'; % parameters to be optimized ('tilts'/'dist'/'all')
par.PlaneFams = []; % which plane families to use ([1,2,5,6,etc.])
par.ThetaBounds = []; % boundaries of restricted theta region to use
par.ShowResults = true; % show graphical output
par.Parameters = []; % parameters like in parameter file
par.GuessTiltY = 0; % guess value for tiltY
par.GuessTiltZ = 0; % guess value for tiltZ
par.GuessDist = 1; % guess value for distance factor
par.BoundsTilts = 5; % boundary of search region for tilts (+-this value in deg)
par.BoundsDist = 0.3; % boundary of search region for distance factor (1+-this value)
par.CostMode = 0; % type of cost function to apply
par.TolX = 1e-6; % optimization end criterion (accuracy) on parameters
par.TolFun = 1e-9; % optimization end criterion (differential) on cost function
par.Calibration = true; % true - uses calibration table; false - uses normal spotpair table
par.ShowCostFun = false; % calculates and shows the cost function on the region [BoundsTilts,BoundsTilts,BoundDist]
par.IncTilts = 0.2; % increment of tilts for showing the cost function
par.IncDist = 0.02; % increment of distance factor for showing the cost function
list{1,1}='OptPar'; list{1,2} = 'parameters to be optimized (tilt / dist /all)';
list{2,1}='PlaneFams'; list{2,2} = 'which plane families to use ([1,2,5,6,etc.])';
list{3,1}='ThetaBounds'; list{3,2} = 'boundaries of restricted theta region to use';
list{4,1}='ShowResults'; list{4,2} = 'show graphical output';
list{5,1}='GuessTiltY'; list{5,2} = 'guess value for tiltY';
list{6,1}='GuessTiltZ'; list{6,2} = 'guess value for tiltZ';
list{7,1}='GuessDist'; list{7,2} = 'guess value for distance factor';
list{8,1}='BoundsTilts'; list{8,2} = 'boundary of search region for tilts (+-this value in deg)';
list{9,1}='BoundsDist'; list{9,2} = 'boundary of search region for distance factor (1+-this value)';
list{10,1}='CostMode'; list{10,2} = 'type of cost function to apply';
list{11,1}='TolX'; list{11,2} = 'optimization end criterion (accuracy) on parameters';
list{12,1}='TolFun'; list{12,2} = 'optimization end criterion (differential) on cost function ';
list{13,1}='ShowCostFun'; list{13,2} = 'calculates and shows the cost function on the region [BoundsTilts,BoundsTilts,BoundDist]';
list{14,1}='IncTilts'; list{14,2} = 'increment of tilts for showing the cost function';
list{15,1}='IncDist'; list{15,2} = 'increment of distance factor for showing the cost function';
par=gtModifyStructure(par, list);
gtMATCHCorrectGeo('OptPar', par.OptPar, 'PlaneFams', par.PlaneFams, 'ThetaBounds', par.ThetaBounds, 'ShowResults', par.ShowResults,...
'GuessTiltY', par.GuessTiltY, 'GuessTiltZ', par.GuessTiltZ, 'GuessDist', par.GuessDist, 'BoundsTilts', par.BoundsTilts, ...
'BoundsDist', par.BoundsDist, 'CostMode', par.CostMode, 'TolX', par.TolX, 'TolFun', par.TolFun, 'Calibration', par.Calibration, ...
'ShowCostFun', par.ShowCostFun, 'IncTilts', par.IncTilts, 'IncDist', par.IncDist);
end
function sfRunNormalButtonFunction(varargin)
% default parameters
par.ShowFindings = false;
par.DoPause = 0;
par.Calibration = false;
par.MaxTime = Inf;
par.UpdateTable = true;
% list/captions
list{1,1}='ShowFindings'; list{1,2}={'Display results for each pair'};
list{2,1}='DoPause'; list{2,2}={'Time to pause on pair, or "button" for mouse click'};
list{3,1}='MaxTime'; list{3,2}={'Max time to spend matching in secs, or Inf'};
list{4,1}='UpdateTable'; list{4,2}={'Update pair table - normally true'};
par=gtModifyStructure(par, list);
gtMATCHDifspots('ShowFindings', par.ShowFindings, 'DoPause', par.DoPause, 'MaxTime', par.MaxTime);
end
function sfModNormalParamsButtonFunction(varargin)
parameters=[];
load('parameters.mat');
list{1,1} = 'thr_ang'; list{1,2} = 'Angular threshold when searching';
list{2,1} = 'thr_ang_scale'; list{2,2} = '???';
list{3,1} = 'thr_max_offset'; list{3,2} = 'Centroid image offset threshold (#images)';
list{4,1} = 'thr_ext_offset'; list{4,2} = 'ExtStart / ExtEnd offset threshold (#images)';
list{5,1} = 'thr_genim_offset'; list{5,2} = 'Threshold for at least one of Centroid/ExtStart/ExtEnd';
list{6,1} = 'corr_rot_to_det'; list{6,2} = 'Correction to rot to det distance (pixels)';
list{7,1} = 'thr_for_corr'; list{7,2} = 'Threshold goodness for using image correction';
list{8,1} = 'thr_goodness'; list{8,2} = 'Acceptable limit: 0:perfect, <0.3:good, >0.7:bad';
list{9,1} = 'thr_intint'; list{9,2} = 'Integrated intensity threshold ratio';
list{10,1} = 'thr_area'; list{10,2} = 'Bounding box area ratio';
list{11,1} = 'thr_bbsize'; list{11,2} = 'Bounding box size ratio';
list{12,1} = 'minsizeX'; list{12,2} = 'Minimum bounding box X size (pixels)';
list{13,1} = 'minsizeY'; list{13,2} = 'Minimum bounding box Y size (pixels)';
list{14,1} = 'tiltY'; list{14,2} = 'Detector tilt around Y-axis (degrees)' ;
list{15,1} = 'tiltZ'; list{15,2} = 'Detector tilt around Z-axis (degrees)';
list{16,1} = 'centmode'; list{16,2} = '0:cent of mass 1:spot correl 2:im correl';
parameters.match=gtModifyStructure(parameters.match, list);
save parameters parameters
end
function sfCorGeoNormalButtonFunction(varargin)
% default parameters
par.OptPar = 'all'; % parameters to be optimized ('tilts'/'dist'/'all')
par.PlaneFams = []; % which plane families to use ([1,2,5,6,etc.])
par.ThetaBounds = []; % boundaries of restricted theta region to use
par.ShowResults = true; % show graphical output
par.Parameters = []; % parameters like in parameter file
par.GuessTiltY = 0; % guess value for tiltY
par.GuessTiltZ = 0; % guess value for tiltZ
par.GuessDist = 1; % guess value for distance factor
par.BoundsTilts = 5; % boundary of search region for tilts (+-this value in deg)
par.BoundsDist = 0.3; % boundary of search region for distance factor (1+-this value)
par.CostMode = 0; % type of cost function to apply
par.TolX = 1e-6; % optimization end criterion (accuracy) on parameters
par.TolFun = 1e-9; % optimization end criterion (differential) on cost function
par.Calibration = false; % true - uses calibration table; false - uses normal spotpair table
par.ShowCostFun = false; % calculates and shows the cost function on the region [BoundsTilts,BoundsTilts,BoundDist]
par.IncTilts = 0.2; % increment of tilts for showing the cost function
par.IncDist = 0.02; % increment of distance factor for showing the cost function
list{1,1}='OptPar'; list{1,2} = 'parameters to be optimized (tilt / dist /all)';
list{2,1}='PlaneFams'; list{2,2} = 'which plane families to use ([1,2,5,6,etc.])';
list{3,1}='ThetaBounds'; list{3,2} = 'boundaries of restricted theta region to use';
list{4,1}='ShowResults'; list{4,2} = 'show graphical output';
list{5,1}='GuessTiltY'; list{5,2} = 'guess value for tiltY';
list{6,1}='GuessTiltZ'; list{6,2} = 'guess value for tiltZ';
list{7,1}='GuessDist'; list{7,2} = 'guess value for distance factor';
list{8,1}='BoundsTilts'; list{8,2} = 'boundary of search region for tilts (+-this value in deg)';
list{9,1}='BoundsDist'; list{9,2} = 'boundary of search region for distance factor (1+-this value)';
list{10,1}='CostMode'; list{10,2} = 'type of cost function to apply';
list{11,1}='TolX'; list{11,2} = 'optimization end criterion (accuracy) on parameters';
list{12,1}='TolFun'; list{12,2} = 'optimization end criterion (differential) on cost function ';
list{13,1}='ShowCostFun'; list{13,2} = 'calculates and shows the cost function on the region [BoundsTilts,BoundsTilts,BoundDist]';
list{14,1}='IncTilts'; list{14,2} = 'increment of tilts for showing the cost function';
list{15,1}='IncDist'; list{15,2} = 'increment of distance factor for showing the cost function';
par=gtModifyStructure(par, list);
gtMATCHCorrectGeo('OptPar', par.OptPar, 'PlaneFams', par.PlaneFams, 'ThetaBounds', par.ThetaBounds, 'ShowResults', par.ShowResults,...
'GuessTiltY', par.GuessTiltY, 'GuessTiltZ', par.GuessTiltZ, 'GuessDist', par.GuessDist, 'BoundsTilts', par.BoundsTilts, ...
'BoundsDist', par.BoundsDist, 'CostMode', par.CostMode, 'TolX', par.TolX, 'TolFun', par.TolFun, 'Calibration', par.Calibration, ...
'ShowCostFun', par.ShowCostFun, 'IncTilts', par.IncTilts, 'IncDist', par.IncDist);
end
function sfQuitButtonFunction(~, ~, ~)
quit=1;
assignin('caller', 'quit', quit)
end
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment