diff --git a/3_pairmatching/gtMATCHCorrectGeo.m b/3_pairmatching/gtMATCHCorrectGeo.m
deleted file mode 100644
index 1cec278bad9cddd03d93755f88e44279e992c8f9..0000000000000000000000000000000000000000
--- a/3_pairmatching/gtMATCHCorrectGeo.m
+++ /dev/null
@@ -1,812 +0,0 @@
-%
-% 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
diff --git a/3_pairmatching/gtMATCHDifspots.m b/3_pairmatching/gtMATCHDifspots.m
deleted file mode 100644
index 8b42531c586009f86964b7947c245c7b82251453..0000000000000000000000000000000000000000
--- a/3_pairmatching/gtMATCHDifspots.m
+++ /dev/null
@@ -1,831 +0,0 @@
-%
-% FUNCTION gtMATCHDifspots(varargin)
-%
-% Finds Friedel pairs by matching diffraction spots in a difspot table from
-% a 360 degree scan from. Outputs pair data in the spotpair table. Goodness
-% of match and correlated difspot ID-s are written in an output file.
-% Matching and its goodness is based on metadata of spots, and image correlation
-% in difficult cases. It can take into account detector tilts for better angular
-% precision.
-% In calibration mode, only good matches are considered and the cristallographic
-% criterion is ignored (no restriction on Bragg angles). This allows for setup
-% calibration (finding exact tilts and detector distance).
-%
-% OPTIONAL INPUT
-% {default values in brackets}
-%
-% Calibration: activates calibration mode {false}
-% UpdateTable: if true, the Spotpair table will be updated {true}
-% ParMatch: matching parameters; if not specified it looks into
-% the parameters file (parameters.match or parameters.match_calib)
-% SavePar: saves the last used matching parameters in the parameters
-% file {false}
-% ShowFindings: if true, the pairs found and candidates are shown in
-% figures {false}
-% DoPause: pause time in sec after a pair is found; if set as 'button',
-% it waits for button press {0}
-%
-%
-% OUTPUT
-%
-% Writes pair data in the spotpairs table (or calibration table in calibration mode)
-% given in the parameters file.
-% Saves a file 'match_outputX.mat' containing:
-% corredpairs: correlated (dubious) matches; all the pairs which required image
-% correlation for selecting the best candidate
-% data order: [pairID, spotA.ID, spotB.ID, ignored candidates]
-% goodness: goodness of fit value of each pair based on their metadata
-%
-%
-% MATCHING PARAMETERS AND TOLERANCES
-%
-% Parameters used for matching should be contained in the field
-% parameters.match or parameters.match_calib {default values in brackets}.
-%
-% centmode: centroid of difspots based on:
-% 0 or []: uses center of mass from difspot table (default)
-% 1: difspots from blobs or edf-s are correlated
-% 2: difspots - from summed fulls - are correlated
-% addconstr: additional constraints on which spots to consider in the database
-%
-% tiltY: scintillator tilt around lab Y (horizontal) axis in degrees {0}
-% tiltZ: scintillator tilt around lab Z (vertical) axis in degreees {0}
-% corr_rot_to_det: correction for rotation-axis - detector dist. in pixels {0}
-%
-% Thresholds for the diffraction vector:
-% thr_ang: 2theta angular deviation; in degrees >0
-% thr_ang_scale: addition to thr_ang linearly scaled with 2theta; degree/degree >=0
-%
-% Thresholds for full image position offsets:
-% thr_max_offset: CentroidImage offset threshold; in #images >=0
-% thr_ext_offset: ExtStartImage and ExtEndImage offset threshold; in #images >=0
-% thr_genim_offset: at least one of the three images should be offset as
-% maximum this value in #images >=0
-%
-% Thresholds for difspot images (search between limits: basevalue/thr & & basevalue*thr):
-% thr_intint: integrated intensity
-% thr_area: area of spots inside the bounding box
-% thr_bbsize: bounding box size
-%
-% Minimum difspot size to be considered:
-% minsizeX: minimum X size in pixels
-% minsizeY: minimum Y size in pixels
-%
-% Thresholds for goodness of match:
-% thr_goodness: maximum accepted goodness of match based on metadata >=0
-% (0 - perfect match, <0.3 - good match, >0.7 - bad match)
-% thr_for_corr: if the difference in "goodness of match" is smaller than this value
-% between difspot candidates of a pair, the best candidate will be
-% selected based on image correlation
-%
-%
-% by Peter Reischig, ESRF, 06/2009
-%
-%
-% Modification Andy 2011 - add MaxTime as a variable, to specify an limit
-% to how long the function will run before returning. This is to allow it
-% to be run interactively on a few spots from a large dataset, looking at
-% the results after say 30 seconds, without having to hit ctrl-c to stop
-% it.
-
-
-%%
-
-
-function gtMATCHDifspots(varargin)
-
-% Number of spots handled in one bunch in processtable
-nofspotspertable = 500; % (500 seems efficient)
-
-
-disp(' ')
-disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
-disp('%% gtMatchDifspots - search for Friedel pairs %%')
-disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
-
-
-%load('parameters.mat');
-tic
-
-% Default parameters
-%par.ParMatch = [];
-par.Parameters = [];
-par.SavePar = false;
-par.ShowFindings = false;
-par.UpdateTable = true;
-par.DoPause = 0;
-par.Calibration = false;
-par.MaxTime = Inf;
-
-par=parse_pv_pairs(par,varargin);
-warning('off','Images:initSize:adjustingMag');
-
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%% Set up match parameters
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-if isempty(par.Parameters)
- load('parameters.mat');;
- par.Parameters = parameters;
-else
- parameters = par.Parameters;
-end
-
-give_warning = false;
-
-if par.Calibration
-
- disp(' running in CALIBRATION MODE ')
-
- if isfield(parameters,'match_calib')
- if ~isfield(parameters.match_calib,'thr_for_corr')
- tmp = gtMatchDefaultParametersCalib;
- parameters.match_calib.thr_for_corr = tmp.thr_for_corr;
- disp(' ')
- disp('Using default "thr_for_corr" which was not found in the parameters file.')
- give_warning = true;
- end
- if ~isfield(parameters.match_calib,'thr_goodness')
- tmp = gtMatchDefaultParametersCalib;
- parameters.match_calib.thr_goodness = tmp.thr_goodness;
- disp(' ')
- disp('Using default "thr_goodness" which was not found in the parameters file.')
- give_warning = true;
- end
- par_match = parameters.match_calib;
- else
- disp('Using default calibration matching parameters.')
- parameters.match_calib = gtMatchDefaultParametersCalib;
- par_match = gtMatchDefaultParametersCalib; % par_match is what is actually gonna be used
- give_warning = true;
- end
-
-else
-
- if isfield(parameters,'match')
- if ~isfield(parameters.match,'thr_for_corr')
- tmp = gtMatchDefaultParameters;
- parameters.match.thr_for_corr = tmp.thr_for_corr;
- disp(' ')
- disp('Using default "thr_for_corr" which was not found in the parameters file.')
- give_warning = true;
- end
- if ~isfield(parameters.match,'thr_goodness')
- tmp = gtMatchDefaultParameters;
- parameters.match.thr_goodness = tmp.thr_goodness;
- disp(' ')
- disp('Using default "thr_goodness" which was not found in the parameters file.')
- give_warning = true;
- end
- par_match = parameters.match;
- else
- disp('Using default matching parameters.')
- parameters.match = gtMatchDefaultParameters;
- par_match = gtMatchDefaultParameters; % par_match is what is actually gonna be used
- give_warning = true;
- end
-
-end
-
-
-if par.SavePar
- save('parameters.mat','parameters')
- disp(' ')
- disp('Actual matching parameters are saved in the parameters file.')
-elseif give_warning
- disp(' ')
- disp('WARNING! Actual matching parameters will NOT be saved in the parameters file.')
-end
-
-
-disp(' ')
-disp('Parameters used:')
-disp(par_match)
-
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%% Short names of tolerances
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-% Tolerances for image position
-thr_ang = par_match.thr_ang;
-thr_ang_scale = par_match.thr_ang_scale;
-thr_max_offset = par_match.thr_max_offset;
-thr_ext_offset = par_match.thr_ext_offset;
-thr_genim_offset = par_match.thr_genim_offset;
-corr_rot_to_det = par_match.corr_rot_to_det;
-
-% Tolerances (for search within the limits: basevalue/thr & & basevalue*thr)
-thr_intint = par_match.thr_intint;
-thr_area = par_match.thr_area;
-thr_bbsize = par_match.thr_bbsize;
-
-% Other options
-centmode = par_match.centmode;
-tiltY = par_match.tiltY;
-tiltZ = par_match.tiltZ;
-thr_for_corr = par_match.thr_for_corr;
-thr_goodness = par_match.thr_goodness;
-minsizeX = par_match.minsizeX;
-minsizeY = par_match.minsizeY;
-
-addconstr = par_match.addconstr;
-if ~isempty(par_match.addconstr)
- addconstr=[' AND ' addconstr];
-end
-
-% Sample geometry in IMAGE coordinates in pixels
-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;
-
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%% Definition of SAMPLE coordinate system
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-
-% Define the Z origin of the SAMPLE coordinate system
-if isfield(parameters.acq,'sampleorigZ')
- sorZ = parameters.acq.sampleorigZ;
-else
- % Set origin of sample coordinate system
- % sorX = on the rotation axis by definition
- % sorY = on the rotation axis by definition
- sample = gtSampleGeoInSampleSystem(parameters);
- sorZ = sample.sorZim;
- parameters.acq.sampleorigZ = sorZ;
- if par.SavePar
- save('parameters.mat','parameters')
- disp(' ')
- disp(sprintf('Using default origin for SAMPLE Z coordinates: %f',sorZ))
- disp('Value has been added to the parameters file.')
- else
- disp(sprintf('Using default origin for SAMPLE Z coordinates: %f',sorZ))
- disp('Value has NOT been added to the parameters file.')
- give_warning = true;
- end
-end
-
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%% Setup tables and load data
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-% Set table names
-if isfield(parameters.acq, 'difA_name')
- difspottable=[parameters.acq.difA_name 'difspot'];
- makeprocesstable=[parameters.acq.difA_name 'process_'];
- processtable=[parameters.acq.difA_name 'process_difspot'];
-else
- difspottable=[parameters.acq.name 'difspot'];
- makeprocesstable=[parameters.acq.name 'process_'];
- processtable=[parameters.acq.name 'process_difspot'];
-end
-
-if par.Calibration
- pairtable = parameters.acq.calib_tablename ;
-else
- pairtable = parameters.acq.pair_tablename ;
-end
-
-% Number of projections per 180 degrees
-nproj = parameters.acq.nproj;
-
-% Theoretical Bragg angles from energy and crystal structure
-disp('Loading theoretical Bragg angles from energy and crystal structure')
-comp_2thetas = parameters.cryst.twoth; % coloumn vector of all valid 2theta angles in degrees
-
-disp(' ')
-
-% Recreate blank tables (overwrite)
-gtDBCreateDifspotTable(makeprocesstable,1)
-
-% Number of spots to be considered ...
-% ... in first half
-nof_spots1 = mym(sprintf(['SELECT count(*) FROM %s WHERE CentroidImage BETWEEN '...
- '%0.20g AND %0.20g AND BoundingBoxXsize>%0.20g AND BoundingBoxYsize>%0.20g %s'],...
- difspottable,0,nproj,minsizeX,minsizeY,addconstr));
-% ... in second half
-nof_spots2 = mym(sprintf(['SELECT count(*) FROM %s WHERE CentroidImage BETWEEN '...
- '%0.20g AND %0.20g AND BoundingBoxXsize>%0.20g AND BoundingBoxYsize>%0.20g %s'],...
- difspottable,nproj,2*nproj,minsizeX,minsizeY,addconstr));
-
-disp(' ')
-disp('Number of spots considered in the first half (A):'), disp(nof_spots1) ;
-disp('Number of spots considered in the second half (B):'), disp(nof_spots2) ;
-
-% Spot A-s from first half of scan (centim -> omega)
-[spotsA,centim] = mym(sprintf(['SELECT difspotID,CentroidImage FROM %s WHERE CentroidImage<=%d'...
- ' AND BoundingBoxXsize>%0.20g AND BoundingBoxYsize>%0.20g %s ORDER',...
- ' BY CentroidImage asc, Integral desc'],difspottable,nproj+thr_max_offset,...
- minsizeX,minsizeY,addconstr));
-
-max_difspotID = mym(sprintf('SELECT max(difspotID) FROM %s',difspottable));
-pairedspots = false(max_difspotID,1); % will indicate if a spot has been paired
-
-% Update pairtable
-if par.UpdateTable
- mym(sprintf('truncate %s',pairtable)) ;
- disp(' ')
- disp('Data in spotpairs table have been deleted, if there were any.') ;
-end
-
-istart=1; % this is obsolete
-
-% Andy - start timing
-c0=clock;
-c0s=c0(4)*3600 + c0(5)*60 + c0(6);
-
-pairsfound = 0;
-corredpairs = [];
-goodness_out = [];
-show_bands = true;
-
-disp(' ')
-disp('Processing data...')
-disp(' ')
-disp(' #Spots A Omega #Pairs found Matching');
-
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%% Main search loop through spots from first half
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-for i=istart:length(spotsA)
-
- sp1.id=spotsA(i); % id of spot A
-
- if pairedspots(sp1.id) % spot has been paired
- continue
- end
-
- % Create new processtable for the next group of spots
- if mod(i-istart,nofspotspertable) == 0
- centim_begin = centim(i)+nproj-thr_max_offset;
- centim_end = centim(min([i+nofspotspertable length(centim)]))+nproj+thr_max_offset;
-
- mym(sprintf('TRUNCATE %s',processtable))
- mysqlcmd=sprintf(['INSERT INTO %s (difspotID,StartImage,EndImage,MaxImage,'...
- 'ExtStartImage,ExtEndImage, Area, CentroidX, CentroidY,CentroidImage,'...
- 'BoundingBoxXorigin,BoundingBoxYorigin,BoundingBoxXsize,BoundingBoxYsize,'...
- 'Integral)'...
- ' SELECT difspotID,StartImage,EndImage,MaxImage,'...
- 'ExtStartImage,ExtEndImage, Area, CentroidX, CentroidY,CentroidImage,'...
- 'BoundingBoxXorigin,BoundingBoxYorigin,BoundingBoxXsize,BoundingBoxYsize,'...
- 'Integral FROM %s WHERE CentroidImage BETWEEN %0.20g AND %0.20g AND '...
- 'BoundingBoxXsize>%0.20g AND BoundingBoxYsize>%0.20g'],...
- processtable,difspottable,centim_begin,centim_end,minsizeX,minsizeY);
- mym(mysqlcmd);
- end
-
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %% Preselection of candidate spots from second half
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
- % Load difspot properties of Spot A
- [sp1.maxim,sp1.extstim,sp1.extendim,sp1.area,sp1.centX,sp1.centY,...
- sp1.bbX,sp1.bbY,sp1.bbXsize,sp1.bbYsize,sp1.integral] =...
- mym(sprintf(['SELECT CentroidImage,ExtStartImage,ExtEndImage,Area,CentroidX,CentroidY,'...
- 'BoundingBoxXorigin,BoundingBoxYorigin,BoundingBoxXsize,BoundingBoxYsize,Integral'...
- ' FROM %s WHERE difspotID=%d %s'],difspottable,sp1.id,addconstr)) ;
-
- % Boundaries of the sample projected on the opposite detector plane from
- % the location of Spot A (this defines the boundaries of the search area).
- % Input is given in IMAGE coordinates in pixels.
- pr_bb = gtMatchProjectedSampleVol(sp1.centX,sp1.centY,sample_radius,sample_top,...
- sample_bot,rot_to_det,parameters.acq.rotx) ;
-
- % Set limits for search between lim1..lim2
- % CentroidImage, ExtStart, Extend limits for pair
- [maxim1,maxim2] = sfCalcBoundaries(sp1.maxim,thr_max_offset,nproj);
- [extstim1,extstim2] = sfCalcBoundaries(sp1.extstim,thr_ext_offset,nproj);
- [extendim1,extendim2] = sfCalcBoundaries(sp1.extendim,thr_ext_offset,nproj);
-
- % At least one of them should also fulfill thr_genim_offset
- [maxim1_gen,maxim2_gen] = sfCalcBoundaries(sp1.maxim,thr_genim_offset,nproj);
- [extstim1_gen,extstim2_gen] = sfCalcBoundaries(sp1.extstim,thr_genim_offset,nproj);
- [extendim1_gen,extendim2_gen] = sfCalcBoundaries(sp1.extendim,thr_genim_offset,nproj);
-
-
- % Select candidate spots2 from database by all the criteria
- [sps2.id, sps2.maxim, sps2.extstim, sps2.extendim, sps2.area,...
- sps2.centX, sps2.centY, sps2.bbX, sps2.bbY, sps2.bbXsize, sps2.bbYsize, sps2.integral] = ...
- mym(sprintf(['SELECT difspotID, CentroidImage, ExtStartImage, ExtEndImage,', ...
- 'Area, CentroidX, CentroidY, BoundingBoxXorigin, BoundingBoxYorigin, BoundingBoxXsize, ', ...
- 'BoundingBoxYsize, Integral FROM %s WHERE ', ...
- '(CentroidX BETWEEN %0.20g AND %0.20g) AND (CentroidY BETWEEN %0.20g AND %0.20g) ', ...
- 'AND (BoundingBoxXsize BETWEEN %0.20g AND %0.20g) AND (BoundingBoxYsize BETWEEN %0.20g AND %0.20g) ', ...
- 'AND (Integral BETWEEN %0.20g AND %0.20g) ', ...
- 'AND (Area BETWEEN %0.20g AND %0.20g) ', ...
- 'AND (CentroidImage BETWEEN %0.20g AND %0.20g) ', ...
- 'AND (ExtStartImage BETWEEN %d AND %d) ', ...
- 'AND (ExtEndImage BETWEEN %d AND %d) ', ...
- 'AND ((CentroidImage BETWEEN %0.20g AND %0.20g) OR ', ...
- '(ExtStartImage BETWEEN %d AND %d) OR ', ...
- '(ExtEndImage BETWEEN %d AND %d)) ', ...
- ' %s'], ...
- processtable,pr_bb(1),pr_bb(3),pr_bb(2),pr_bb(4), ...
- sp1.bbXsize/thr_bbsize, sp1.bbXsize*thr_bbsize, ...
- sp1.bbYsize/thr_bbsize, sp1.bbYsize*thr_bbsize, ...
- sp1.integral/thr_intint, sp1.integral*thr_intint, ...
- sp1.area/thr_area, sp1.area*thr_area, ...
- maxim1, maxim2, ...
- extstim1, extstim2, ...
- extendim1, extendim2, ...
- maxim1_gen, maxim2_gen, ...
- extstim1_gen, extstim2_gen, ...
- extendim1_gen, extendim2_gen, ...
- addconstr)) ;
-
- % Delete candidates that are already paired
- sps2.maxim(pairedspots(sps2.id)) =[];
- sps2.extstim(pairedspots(sps2.id)) =[];
- sps2.extendim(pairedspots(sps2.id)) =[];
- sps2.area(pairedspots(sps2.id)) =[];
- sps2.centX(pairedspots(sps2.id)) =[];
- sps2.centY(pairedspots(sps2.id)) =[];
- sps2.bbX(pairedspots(sps2.id)) =[];
- sps2.bbY(pairedspots(sps2.id)) =[];
- sps2.bbXsize(pairedspots(sps2.id)) =[];
- sps2.bbYsize(pairedspots(sps2.id)) =[];
- sps2.integral(pairedspots(sps2.id)) =[];
- sps2.id(pairedspots(sps2.id)) =[];
-
-
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %% Check theta angles; goodness of match
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
- goodness=Inf(length(sps2.id),1); % will conatin goodness of match values
- thetatypes=NaN(length(sps2.id),1); % will contain the plane family
- pairj=[]; % will contain the final pair index, if anything is found
-
- for j=1:length(sps2.id)
-
- % Check theta angle and plane family for each candidate
- [theta_OK,angle_diff,thetatype_j] = sfCheckTheta(sp1.centX,sp1.centY,sps2.centX(j),sps2.centY(j),...
- rot_to_det,tiltY,tiltZ,parameters.acq.rotx,sorZ,comp_2thetas,thr_ang,thr_ang_scale);
-
- if theta_OK==true
-
- % Relative deviations of spot properties ...
- reldevs = [max(sps2.area(j)/sp1.area, sp1.area/sps2.area(j))-1; ...
- max(sps2.integral(j)/sp1.integral, sp1.integral/sps2.integral(j))-1; ...
- max(sps2.bbXsize(j)/sp1.bbXsize, sp1.bbXsize/sps2.bbXsize(j))-1; ...
- max(sps2.bbYsize(j)/sp1.bbYsize, sp1.bbYsize/sps2.bbYsize(j))-1; ...
- abs(sps2.extstim(j)-sp1.extstim-nproj); ...
- abs(sps2.extendim(j)-sp1.extendim-nproj); ...
- abs(sps2.maxim(j)-sp1.maxim-nproj); ...
- abs(angle_diff)];
-
- % ... weighted by corresponding tolerances ...
- sqrsumnew = reldevs./[thr_area; thr_intint; thr_bbsize; thr_bbsize; ...
- thr_ext_offset; thr_ext_offset; thr_max_offset; ...
- thr_ang];
-
- % .. give the goodness of match for each candidate:
- goodness(j) = sqrt(sum(sqrsumnew.^2)); % goodness of match
-
- thetatypes(j) = thetatype_j; % plane family
-
- end
-
- end
-
-
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %% Selection of best candidate
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
- % Number of candidates found
-
- if length(sps2.id)==1 && theta_OK % only one spot meets search criteria in table
- if goodness < thr_goodness
- pairj=1; % if goodness is smaller than limit -> accept match
- thetatype=thetatypes; % corresponding plane family
- else
- pairj=[]; % otherwise ignore finding
- end
- elseif length(sps2.id)>1 % more than one spot meet search criteria
-
- [bestgoodness,pairj]=min(goodness); % best goodness of match
-
- if bestgoodness < thr_goodness % at least one valid pair has been found
-
- thetatype=thetatypes(pairj);
-
- % If there are other candidates that match spotA almost as good as the
- % best one, all those will be correlated with spotA, and best
- % matching will be kept as its pair
- tobecorred=false(length(sps2.id),1);
- for j=1:length(sps2.id)
- % if goodness of match is close to best
- if goodness(j)<bestgoodness+thr_for_corr && goodness(j) < thr_goodness
- tobecorred(j)=true;
- end
- end
-
- % Correlate difspot candidates with original if properties are too close
- if sum(tobecorred)>1
- % Create blank image for corrrelation (max size needed)
- pad0=zeros(max([sps2.bbYsize(tobecorred); sp1.bbYsize]), ...
- max([sps2.bbXsize(tobecorred); sp1.bbXsize]));
-
- % Image of spotA
- imA=gtGetSummedDifSpot(sp1.id,parameters,'connected',1); % difspot image A
- impadA=pad0;
- impadA(1:size(imA,1),1:size(imA,2))=imA; % fill in image
-
- corrval=zeros(length(sps2.id),1);
-
- for j=1:length(sps2.id)
- if tobecorred(j)
-
- % Image of spotB
- imB=fliplr(gtGetSummedDifSpot(sps2.id(j),parameters,'connected',1));
- impadB=pad0;
- impadB(1:size(imB,1),1:size(imB,2))=imB;
-
- corrval(j)=max2(abs(normxcorr2(impadA,impadB))); % normalized correlation value
-
- % To avoid licence shortage, the function is renamed ...
- %corrval(j)=max2(abs(gtMatch_normxcorr2(impadA,impadB))); % normalized correlation value
-
- end
- end
-
- [maxcorrval,pairj]=max(corrval); % pairj indicates index of best fitting spot
- thetatype=thetatypes(pairj); % plane family
-
- % Correlated pairs are added to corredpairs:
- % data order: [pairID, spotA.ID, spotB.ID, other candidates]
- tobecorred(pairj)=false;
- corredpairs{length(corredpairs)+1}=[pairsfound+1, sp1.id, sps2.id(pairj) sps2.id(tobecorred)'];
-
- end
-
- else
- pairj=[]; % none of the candidates makes a valid theta angle
- end
-
- end
-
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %% Compute output parameters; load database
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
- if ~isempty(pairj) % a pair has been found
- pairsfound=pairsfound+1 ;
- goodness_out=[goodness_out; goodness(pairj)];
-
- % Put pair data in sp2
- sp2.id=sps2.id(pairj);
- sp2.maxim=sps2.maxim(pairj);
- sp2.extstim=sps2.extstim(pairj);
- sp2.extendim=sps2.extendim(pairj);
- sp2.area=sps2.area(pairj);
- sp2.centX=sps2.centX(pairj);
- sp2.centY=sps2.centY(pairj);
- sp2.integral=sps2.integral(pairj);
- sp2.bbXsize=sps2.bbXsize(pairj);
- sp2.bbYsize=sps2.bbYsize(pairj);
- sp2.bbX=sps2.bbX(pairj);
- sp2.bbY=sps2.bbY(pairj);
-
- % Leave only the wrong candidates in sps2
- sps2.maxim(pairj)=[];
- sps2.extstim(pairj)=[];
- sps2.extendim(pairj)=[];
- sps2.area(pairj)=[];
- sps2.centX(pairj)=[];
- sps2.centY(pairj)=[];
- sps2.bbX(pairj)=[];
- sps2.bbY(pairj)=[];
- sps2.bbXsize(pairj)=[];
- sps2.bbYsize(pairj)=[];
- sps2.integral(pairj)=[];
- sps2.id(pairj)=[];
-
- % Mark the two spots as paired
- pairedspots(sp1.id)=true;
- pairedspots(sp2.id)=true;
-
- % Centre of mass positions or difspot correlation used for the calculation
- % of pair properties
- switch centmode
- case {0, []} % use center of mass
-
- case 1 % correlate difspots (edf-s or blobs)
- [corrx,corry]=gtCorrelateDifspots(sp1.id,sp2.id,parameters,centmode,1);
- sp2.centX=sp2.centX-corrx; % the negative of corr is applied, since
- % corr referred to the flipped image
- sp2.centY=sp2.centY+corry;
-
- case 2 % correlate summed fulls under the difspot bb
- [corrx,corry]=gtCorrelateDifspots(sp1.id,sp2.id,parameters,centmode,2);
- sp2.centX=sp2.centX-corrx; % the negative of corr is applied, since
- % corr referred to the flipped image
- sp2.centY=sp2.centY+corry;
- end
-
- % Diffraction vector in LAB coordinates:
- dv_lab = gtMatchDiffVecInLab(sp1.centX,sp1.centY,sp2.centX,sp2.centY,...
- rot_to_det,tiltY,tiltZ,parameters.acq.rotx,sorZ) ;
-
-
- % Theta and eta from the pair:
- theta = gtMatchThetaOfPair(dv_lab) ;
- eta = gtMatchEtaOfPoint(dv_lab) ;
- etaB = mod(180-eta,360); % eta of spotB
-
-
- % Plane normal in LAB coordinates:
- pln_lab = gtMatchPlaneNormInLab(dv_lab);
-
- % Omega correction:
- omega1 = 180/nproj*sp1.maxim;
- omega2 = 180/nproj*sp2.maxim; % (sp2.maxim should always be larger)
- omegaA = (omega1+omega2-180)/2; % the average of the two
- omegaB = omegaA+180; % omega of spotB
-
-
- % Diffraction vector and plane normal in SAMPLE coordinates:
- [dv_sam(1), dv_sam(2), dv_sam(3)] = gtMatchLab2Sam(dv_lab(1), dv_lab(2), dv_lab(3), omegaA) ;
- [pln_sam(1), pln_sam(2), pln_sam(3)] = gtMatchLab2Sam(pln_lab(1), pln_lab(2), pln_lab(3), omegaA) ;
-
- % Spot centers from SCINTILLATOR PLANE into LAB coordinates (in 3D):
- [labcentA(1),labcentA(2),labcentA(3)] = gtMatchSc2Lab(sp1.centX,sp1.centY,rot_to_det,tiltY,tiltZ,parameters.acq.rotx,sorZ) ;
- [labcentB(1),labcentB(2),labcentB(3)] = gtMatchSc2Lab(sp2.centX,sp2.centY,rot_to_det,tiltY,tiltZ,parameters.acq.rotx,sorZ) ;
-
- % Spot centers from LAB into SAMPLE coordinates (with corrected omega-s !!):
- [samcentA(1),samcentA(2),samcentA(3)] = gtMatchLab2Sam(labcentA(1), labcentA(2), labcentA(3), omegaA) ;
- [samcentB(1),samcentB(2),samcentB(3)] = gtMatchLab2Sam(labcentB(1), labcentB(2), labcentB(3), omegaB) ;
-
- % 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 par.UpdateTable
- % 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, goodness, '...
- '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,%0.20g,%d)'],...
- pairtable,sp1.id,sp2.id,goodness(pairj),theta,eta,omegaA,etaB,omegaB,av_intint,av_bbXsize,av_bbYsize,...
- samcentA(1),samcentA(2),samcentA(3),samcentB(1),samcentB(2),samcentB(3),...
- dv_sam(1),dv_sam(2),dv_sam(3),pln_sam(1),pln_sam(2),pln_sam(3),thetatype) ;
- mym(mysqlcmd);
- end
-
-
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %% Display findings; create output file; graphical output
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
- % Show images of paired spots
- if par.ShowFindings
-
- if par.Calibration
- gtMatchShowPairFigure(1,parameters.acq,parameters.match_calib,sp1,sp2,...
- pr_bb,comp_2thetas,show_bands,sps2,[],[],omegaA,theta,eta,pairsfound)
- else
- gtMatchShowPairFigure(1,parameters.acq,parameters.match,sp1,sp2,...
- pr_bb,comp_2thetas,show_bands,sps2,[],[],omegaA,theta,eta,pairsfound)
- end
-
- show_bands=false;
- if strcmpi('button',par.DoPause)
-
- waitforbuttonpress
- else
- pause(par.DoPause)
- end
- end
-
- end % of loop pair found
-
- % Display progress
- if mod(i,100)==0 && pairsfound>0
- disp(sprintf('%10d %10.1f %10d %3.2f%%',i,omegaA,pairsfound,pairsfound/(i-istart+1)*100)) ;
-
- % timing check
- c1=clock;
- c1s=c1(4)*3600 + c1(5)*60 + c1(6);
- cdif=c1s-c0s;
- if cdif>par.MaxTime
- disp(' ');disp('Reached specified MaxTime'); pause(1); disp(' ')
- break
- end
- end
-
-end % end of main loop for spotsA
-
-if pairsfound==0
- disp(' ')
- disp('There has been no pair found. Try to reset matching parameters.')
- disp(' ')
- return
-end
-
-disp(sprintf('%10d %10.1f %10d %3.2f%%',i,omegaA,pairsfound,pairsfound/(i-istart+1)*100)) ;
-
-disp(' ')
-disp(' ')
-disp('Number of spots considered in the first half (A):'), disp(nof_spots1) ;
-disp('Number of spots considered in the second half (B):'), disp(nof_spots2) ;
-disp('Number of pairs found:'), disp(pairsfound) ;
-disp(' ')
-
-if isempty(corredpairs)
- disp('No dubious matches have been found.')
-else
- disp('Dubious matches have been found and saved in output file.')
-end
-
-% Save correlated (dubious) pairs in output file
-fname=gtLastFileName('match_output','new');
-goodness=goodness_out;
-save(fname,'corredpairs','goodness')
-disp(' ')
-disp(sprintf('Output has been saved in %s',fname))
-
-disp(' ');
-toc
-disp(' ')
-
-if ~par.SavePar && give_warning
- disp('WARNING! Actual matching parameters have NOT been saved in parameters file.')
- disp(' ')
-end
-
-if ~par.UpdateTable
- disp('WARNING! Spotpair table has NOT been updated.')
- disp(' ')
-else
- if par.Calibration
- gtMATCHEvaluate('Parameters',parameters,'Calibration',1) % graphical evaluation of matching
- else
- gtMATCHEvaluate('Parameters',parameters) % graphical evaluation of matching
- end
-end
-
-
-end % of function
-
-
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%% Sub-functions used
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-
-function [b1,b2]=sfCalcBoundaries(imno1,tol,nproj)
-% Given an image number, returns the boundaries of two ranges (as image numbers)
-% with an offset of 180degree +-tol offset.
-
- b1=imno1+nproj-tol;
- b2=imno1+nproj+tol;
-
-end
-
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-function [theta_OK,angle_diff,thetatype]=sfCheckTheta(scXA,scYA,scXB,scYB,...
- rottodet,tiltY,tiltZ,rotx,sorZ,comp_2thetas,thr_ang,thr_ang_scale)
-% Checks whether the position of two spots gives a diffraction vector which
-% corresponds to a valid 2theta angle within the given threshold.
-
- % diffraction vector:
- dv = gtMatchDiffVecInLab(scXA,scYA,scXB,scYB,rottodet,tiltY,tiltZ,rotx,sorZ);
-
- % two theta:
- angle = 2*gtMatchThetaOfPair(dv);
-
- % discrepancy in 2theta
- diff = angle-comp_2thetas;
-
- if any(abs(diff) <= thr_ang+comp_2thetas*thr_ang_scale)
- theta_OK=true ;
- [mini,thetatype]=min(abs(diff));
- angle_diff=diff(thetatype);
- else
- theta_OK=false ;
- angle_diff=[];
- thetatype=[] ;
- end
-
-end
diff --git a/3_pairmatching/gtMATCHEvaluate.m b/3_pairmatching/gtMATCHEvaluate.m
deleted file mode 100644
index aa5fcfdb37ac2fac80c20aefc0db32bc48a87d0d..0000000000000000000000000000000000000000
--- a/3_pairmatching/gtMATCHEvaluate.m
+++ /dev/null
@@ -1,479 +0,0 @@
-%
-% 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
-
-
-
diff --git a/3_pairmatching/gtMatchDefaultParameters.m b/3_pairmatching/gtMatchDefaultParameters.m
deleted file mode 100644
index 828dc6b5f274ad0a18a8582c9b17d987e5be9fcd..0000000000000000000000000000000000000000
--- a/3_pairmatching/gtMatchDefaultParameters.m
+++ /dev/null
@@ -1,31 +0,0 @@
-%
-% 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
diff --git a/3_pairmatching/gtMatchDefaultParametersCalib.m b/3_pairmatching/gtMatchDefaultParametersCalib.m
deleted file mode 100644
index 380944932016769fabf207f77bfbdbde10c85b61..0000000000000000000000000000000000000000
--- a/3_pairmatching/gtMatchDefaultParametersCalib.m
+++ /dev/null
@@ -1,33 +0,0 @@
-%
-% 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
diff --git a/3_pairmatching/gtMatchDiffVecInLab.m b/3_pairmatching/gtMatchDiffVecInLab.m
deleted file mode 100644
index 500cf1d53654b07da93820be6eb4b57dea2a16a2..0000000000000000000000000000000000000000
--- a/3_pairmatching/gtMatchDiffVecInLab.m
+++ /dev/null
@@ -1,38 +0,0 @@
-%
-% 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
diff --git a/3_pairmatching/gtMatchEtaOfPoint.m b/3_pairmatching/gtMatchEtaOfPoint.m
deleted file mode 100644
index 2db28c8e616a950b3511f64a95337eb540bb5eda..0000000000000000000000000000000000000000
--- a/3_pairmatching/gtMatchEtaOfPoint.m
+++ /dev/null
@@ -1,56 +0,0 @@
-%
-% 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
diff --git a/3_pairmatching/gtMatchLab2Sam.m b/3_pairmatching/gtMatchLab2Sam.m
deleted file mode 100644
index 34efa6ae26809744058db6d188bfedd13c9deb0b..0000000000000000000000000000000000000000
--- a/3_pairmatching/gtMatchLab2Sam.m
+++ /dev/null
@@ -1,32 +0,0 @@
-%
-% 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
diff --git a/3_pairmatching/gtMatchPlaneNormInLab.m b/3_pairmatching/gtMatchPlaneNormInLab.m
deleted file mode 100644
index 3f9906809a6c06a6f6a2654b06a9116eab9babfa..0000000000000000000000000000000000000000
--- a/3_pairmatching/gtMatchPlaneNormInLab.m
+++ /dev/null
@@ -1,37 +0,0 @@
-%
-% 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
diff --git a/3_pairmatching/gtMatchProjectedSampleVol.m b/3_pairmatching/gtMatchProjectedSampleVol.m
deleted file mode 100644
index f33d5e0a18884d0d847f96bba56795f5147deff4..0000000000000000000000000000000000000000
--- a/3_pairmatching/gtMatchProjectedSampleVol.m
+++ /dev/null
@@ -1,112 +0,0 @@
-%
-% 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;
-
-
diff --git a/3_pairmatching/gtMatchSc2Lab.m b/3_pairmatching/gtMatchSc2Lab.m
deleted file mode 100644
index b4701f37537ec8294c8f4f31e56d22b368101961..0000000000000000000000000000000000000000
--- a/3_pairmatching/gtMatchSc2Lab.m
+++ /dev/null
@@ -1,48 +0,0 @@
-%
-% 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
diff --git a/3_pairmatching/gtMatchShowPairFigure.m b/3_pairmatching/gtMatchShowPairFigure.m
deleted file mode 100644
index b5ba8ffab79ddf6afd1c62958b7412a059f638d3..0000000000000000000000000000000000000000
--- a/3_pairmatching/gtMatchShowPairFigure.m
+++ /dev/null
@@ -1,244 +0,0 @@
-%
-% 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
diff --git a/3_pairmatching/gtMatchThetaOfPair.m b/3_pairmatching/gtMatchThetaOfPair.m
deleted file mode 100644
index c2ae0e351ece15fce30969fc5902861e4c4c46f7..0000000000000000000000000000000000000000
--- a/3_pairmatching/gtMatchThetaOfPair.m
+++ /dev/null
@@ -1,16 +0,0 @@
-%
-% 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
diff --git a/3_pairmatching/gtSetupPairMatching.m b/3_pairmatching/gtSetupPairMatching.m
deleted file mode 100644
index e4b56a1ed5a8569b461e4076a328200d0de393a9..0000000000000000000000000000000000000000
--- a/3_pairmatching/gtSetupPairMatching.m
+++ /dev/null
@@ -1,247 +0,0 @@
-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
-