Commit 61319dfb authored by Yaroslav Gevorkov's avatar Yaroslav Gevorkov
Browse files

new feature: delete fitted peaks from provided ones in crystfel adaption

parent a008767d
......@@ -44,8 +44,10 @@ class PinkIndexer
PinkIndexer(const ExperimentSettings& experimentSettings, ConsideredPeaksCount consideredPeaksCount, AngleResolution angleResolution,
RefinementType refinementType, float maxResolutionForIndexing_1_per_A);
int indexPattern(Lattice& indexedLattice, Eigen::RowVectorXf& intensities, const Eigen::Matrix2Xf& detectorPeaks_m, int threadCount);
int indexPattern(Lattice& indexedLattice, Eigen::RowVectorXf& intensities, const Eigen::Matrix3Xf& meanReciprocalPeaks_1_per_A, int threadCount);
int indexPattern(Lattice& indexedLattice, Eigen::Array<bool, Eigen::Dynamic, 1>& fittedPeaks, Eigen::RowVectorXf& intensities,
const Eigen::Matrix2Xf& detectorPeaks_m, int threadCount);
int indexPattern(Lattice& indexedLattice, Eigen::Array<bool, Eigen::Dynamic, 1>& fittedPeaks, Eigen::RowVectorXf& intensities,
const Eigen::Matrix3Xf& meanReciprocalPeaks_1_per_A, int threadCount);
private:
void reducePeakCount(Eigen::Matrix3Xf& ucsDirections, Eigen::Array2Xf& ucsBorderNorms, Eigen::RowVectorXf& intensities,
......
......@@ -14,6 +14,8 @@ class Refinement
void refineVariableLattice(Lattice& lattice, const Eigen::Matrix3Xf& ucsDirections, const Eigen::Matrix2Xf& ucsBorderNorms);
int getFittedPeaksCount(Lattice& lattice, const Eigen::Matrix3Xf& ucsDirections, const Eigen::Matrix2Xf& ucsBorderNorms);
int getFittedPeaks(Lattice& lattice, Eigen::Array<bool, Eigen::Dynamic, 1>& fittedPeaks, const Eigen::Matrix3Xf& ucsDirections,
const Eigen::Matrix2Xf& ucsBorderNorms);
void setTolerance(float tolerance)
{
......
......@@ -44,7 +44,7 @@ PinkIndexer* PinkIndexer_new(ExperimentSettings* experimentSettings, consideredP
refinementType_t refinementType, float maxResolutionForIndexing_1_per_A);
void PinkIndexer_delete(PinkIndexer* pinkIndexer);
int PinkIndexer_indexPattern(PinkIndexer* pinkIndexer, Lattice_t* indexedLattice, const reciprocalPeaks_1_per_A_t meanReciprocalPeaks_1_per_A,
int PinkIndexer_indexPattern(PinkIndexer* pinkIndexer, Lattice_t* indexedLattice, reciprocalPeaks_1_per_A_t* meanReciprocalPeaks_1_per_A,
const float* intensities, int peakCount, int threadCount);
......
......@@ -26,14 +26,16 @@ PinkIndexer::PinkIndexer(const ExperimentSettings& experimentSettings, Considere
sinogram.setSinogramAngleResolution(angleResolution_deg);
}
int PinkIndexer::indexPattern(Lattice& indexedLattice, Eigen::RowVectorXf& intensities, const Matrix3Xf& meanReciprocalPeaks_1_per_A, int threadCount)
int PinkIndexer::indexPattern(Lattice& indexedLattice, Eigen::Array<bool, Eigen::Dynamic, 1>& fittedPeaks, Eigen::RowVectorXf& intensities,
const Matrix3Xf& meanReciprocalPeaks_1_per_A, int threadCount)
{
Matrix2Xf detectorPeaks_m;
reciprocalToRealProjection.project(detectorPeaks_m, meanReciprocalPeaks_1_per_A);
return indexPattern(indexedLattice, intensities, detectorPeaks_m, threadCount);
return indexPattern(indexedLattice, fittedPeaks, intensities, detectorPeaks_m, threadCount);
}
int PinkIndexer::indexPattern(Lattice& indexedLattice, Eigen::RowVectorXf& intensities, const Matrix2Xf& detectorPeaks_m, int threadCount)
int PinkIndexer::indexPattern(Lattice& indexedLattice, Eigen::Array<bool, Eigen::Dynamic, 1>& fittedPeaks, Eigen::RowVectorXf& intensities,
const Matrix2Xf& detectorPeaks_m, int threadCount)
{
Matrix3Xf ucsDirections;
Array2Xf ucsBorderNorms;
......@@ -69,7 +71,7 @@ int PinkIndexer::indexPattern(Lattice& indexedLattice, Eigen::RowVectorXf& inten
case RefinementType::none:
break;
case RefinementType::fixedLatticeParameters:
refinement.setTolerance(min(finalRefinementTolerance*2.0, 0.12));
refinement.setTolerance(min(finalRefinementTolerance * 2.0, 0.12));
refinement.refineFixedLattice(indexedLattice, ucsDirections, ucsBorderNorms);
refinement.setTolerance(finalRefinementTolerance);
refinement.refineFixedLattice(indexedLattice, ucsDirections, ucsBorderNorms);
......@@ -92,7 +94,7 @@ int PinkIndexer::indexPattern(Lattice& indexedLattice, Eigen::RowVectorXf& inten
break;
}
return refinement.getFittedPeaksCount(indexedLattice, ucsDirections, ucsBorderNorms);
return refinement.getFittedPeaks(indexedLattice, fittedPeaks, ucsDirections, ucsBorderNorms);
// int fittedPeaksCount = refinement.getFittedPeaksCount(indexedLattice, ucsDirections, ucsBorderNorms);
// if ((fittedPeaksCount >= 20 && fittedPeaksCount >= detectorPeaks_m.cols() * 0.3) || fittedPeaksCount >= detectorPeaks_m.cols() * 0.5)
......
......@@ -31,7 +31,7 @@ void Refinement::refineVariableLattice(Lattice& lattice, const Matrix3Xf& ucsDir
meanDefects[0] = getMeanDefect(basis, ucsDirections, ucsBorderNorms);
for (int i = 0; i < maxStepsCount; i++)
{
//cout << meanDefects[i] << endl;
// cout << meanDefects[i] << endl;
Array33f gradient;
Matrix3f offsetBasis = basis;
......@@ -66,7 +66,7 @@ void Refinement::refineVariableLattice(Lattice& lattice, const Matrix3Xf& ucsDir
gradient = gradient / norm * stepSize;
if (norm == 0)
{
//throw WrongUsageException("Numerical problems! Delta has been chosen too small for current lattice!\n");
// throw WrongUsageException("Numerical problems! Delta has been chosen too small for current lattice!\n");
break;
}
......@@ -106,7 +106,7 @@ void Refinement::refineFixedLattice(Lattice& lattice, const Matrix3Xf& ucsDirect
meanDefects[0] = getMeanDefect(basis, ucsDirections, ucsBorderNorms);
for (int i = 0; i < maxStepsCount; i++)
{
//cout << meanDefects[i] << endl;
// cout << meanDefects[i] << endl;
Vector3f gradient;
gradient(0) = getMeanDefect(rotX * basis, ucsDirections, ucsBorderNorms) - meanDefects[i];
......@@ -140,6 +140,15 @@ int Refinement::getFittedPeaksCount(Lattice& lattice, const Eigen::Matrix3Xf& uc
return fittedPeaksCount;
}
int Refinement::getFittedPeaks(Lattice& lattice, Eigen::Array<bool, Eigen::Dynamic, 1>& fittedPeaks, const Eigen::Matrix3Xf& ucsDirections,
const Eigen::Matrix2Xf& ucsBorderNorms)
{
getDefects(defects, lattice.getBasis(), ucsDirections, ucsBorderNorms);
fittedPeaks = (defects < tolerance);
return fittedPeaks.count();
}
void Refinement::getDefects(ArrayXf& defects, const Matrix3f& basis, const Matrix3Xf& ucsDirections, const Matrix2Xf& ucsBorderNorms)
{
......
......@@ -83,19 +83,20 @@ extern "C" void PinkIndexer_delete(PinkIndexer* pinkIndexer)
delete pinkIndexer;
}
extern "C" int PinkIndexer_indexPattern(PinkIndexer* pinkIndexer, Lattice_t* indexedLattice, const reciprocalPeaks_1_per_A_t meanReciprocalPeaks_1_per_A,
extern "C" int PinkIndexer_indexPattern(PinkIndexer* pinkIndexer, Lattice_t* indexedLattice, reciprocalPeaks_1_per_A_t* meanReciprocalPeaks_1_per_A,
const float* intensities, int peakCount, int threadCount)
{
Eigen::Matrix3Xf meanReciprocalPeaks_1_per_A_matrix(3, meanReciprocalPeaks_1_per_A.peakCount);
for (int i = 0; i < meanReciprocalPeaks_1_per_A.peakCount; i++)
Eigen::Matrix3Xf meanReciprocalPeaks_1_per_A_matrix(3, meanReciprocalPeaks_1_per_A->peakCount);
for (int i = 0; i < meanReciprocalPeaks_1_per_A->peakCount; i++)
{
meanReciprocalPeaks_1_per_A_matrix.col(i) << meanReciprocalPeaks_1_per_A.coordinates_x[i], meanReciprocalPeaks_1_per_A.coordinates_y[i],
meanReciprocalPeaks_1_per_A.coordinates_z[i];
meanReciprocalPeaks_1_per_A_matrix.col(i) << meanReciprocalPeaks_1_per_A->coordinates_x[i], meanReciprocalPeaks_1_per_A->coordinates_y[i],
meanReciprocalPeaks_1_per_A->coordinates_z[i];
}
Lattice indexedLattice_class;
Eigen::RowVectorXf intensities_class = Eigen::Map<Eigen::RowVectorXf>((float*)intensities, 1, peakCount);
int matchedPeaksCount = pinkIndexer->indexPattern(indexedLattice_class, intensities_class, meanReciprocalPeaks_1_per_A_matrix, threadCount);
Eigen::Array<bool, Eigen::Dynamic, 1> fittedPeaks;
int matchedPeaksCount = pinkIndexer->indexPattern(indexedLattice_class, fittedPeaks, intensities_class, meanReciprocalPeaks_1_per_A_matrix, threadCount);
Eigen::Matrix3f basis = indexedLattice_class.getBasis();
indexedLattice->ax = basis(0, 0);
......@@ -108,5 +109,16 @@ extern "C" int PinkIndexer_indexPattern(PinkIndexer* pinkIndexer, Lattice_t* ind
indexedLattice->cy = basis(1, 2);
indexedLattice->cz = basis(2, 2);
int unfittedPeaksCount = 0;
for (int i = 0; i < fittedPeaks.size(); i++)
{
meanReciprocalPeaks_1_per_A->coordinates_x[unfittedPeaksCount] = meanReciprocalPeaks_1_per_A->coordinates_x[i];
meanReciprocalPeaks_1_per_A->coordinates_y[unfittedPeaksCount] = meanReciprocalPeaks_1_per_A->coordinates_y[i];
meanReciprocalPeaks_1_per_A->coordinates_z[unfittedPeaksCount] = meanReciprocalPeaks_1_per_A->coordinates_z[i];
unfittedPeaksCount++;
}
meanReciprocalPeaks_1_per_A->peakCount = unfittedPeaksCount;
return matchedPeaksCount;
}
\ No newline at end of file
......@@ -49,14 +49,15 @@ void testPinkIndexer()
PinkIndexer::ConsideredPeaksCount consideredPeaksCount = PinkIndexer::ConsideredPeaksCount::standard;
PinkIndexer::AngleResolution angleResolution = PinkIndexer::AngleResolution::standard;
PinkIndexer::RefinementType refinementType = PinkIndexer::RefinementType::firstFixedThenVariableLatticeParameters;
float maxResolutionForIndexing_1_per_A = 0.2;//numeric_limits<float>::max();
float maxResolutionForIndexing_1_per_A = 0.2; // numeric_limits<float>::max();
Lattice indexedLattice;
PinkIndexer pinkIndexer(experimentSettings, consideredPeaksCount, angleResolution, refinementType, maxResolutionForIndexing_1_per_A);
{
auto tmp = Chronometer("sinogram");
int threadCount = 6;
pinkIndexer.indexPattern(indexedLattice, intensities, peaksOnDetector_m, threadCount);
Eigen::Array<bool, Eigen::Dynamic, 1> fittedPeaks;
pinkIndexer.indexPattern(indexedLattice, fittedPeaks, intensities, peaksOnDetector_m, threadCount);
}
SimpleDiffractionPatternPrediction simpleDiffractionPatternPrediction(experimentSettings);
......@@ -101,7 +102,7 @@ void testRefinement()
loadEigenMatrixFromDisk(ucsDirections, "C:\\DesyFiles\\workspaces\\VisualStudio_workspace\\pinkIndexer\\workfolder\\ucsDirections");
loadEigenMatrixFromDisk(basis, "C:\\DesyFiles\\workspaces\\VisualStudio_workspace\\pinkIndexer\\workfolder\\bestFitBasis");
float tolerance = 0.15;
float tolerance = 0.15;
Lattice lattice(basis);
Refinement refinement(tolerance);
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment