Commit a4dfe325 authored by Yaroslav Gevorkov's avatar Yaroslav Gevorkov
Browse files

added a bit lazy evaluation in refinement

parent 95c21f32
......@@ -27,8 +27,10 @@ class Refinement
}
private:
void getDefects(Eigen::ArrayXf& defects, const Eigen::Matrix3f& basis, const Eigen::Matrix3Xf& ucsDirections, const Eigen::Matrix2Xf& ucsBorderNorms);
double getMeanDefect(const Eigen::Matrix3f& basis, const Eigen::Matrix3Xf& ucsDirections, const Eigen::Matrix2Xf& ucsBorderNorms);
void getDefects(Eigen::ArrayXf& defects, const Eigen::Matrix3f& basis, const Eigen::Matrix3Xf& ucsDirections, const Eigen::Matrix2Xf& ucsBorderNorms,
bool significantChangesToPreviousCall = true);
double getMeanDefect(const Eigen::Matrix3f& basis, const Eigen::Matrix3Xf& ucsDirections, const Eigen::Matrix2Xf& ucsBorderNorms,
bool significantChangesToPreviousCall = true);
float tolerance;
......@@ -45,4 +47,5 @@ class Refinement
Eigen::Matrix3Xf defectVectors_relative;
Eigen::Array<bool, 1, Eigen::Dynamic> notPredictablePeaks;
Eigen::ArrayXf meanDefects;
int validCandidatePeaksCount;
};
......@@ -45,7 +45,7 @@ PinkIndexer* PinkIndexer_new(ExperimentSettings* experimentSettings, consideredP
void PinkIndexer_delete(PinkIndexer* pinkIndexer);
int PinkIndexer_indexPattern(PinkIndexer* pinkIndexer, Lattice_t* indexedLattice, reciprocalPeaks_1_per_A_t* meanReciprocalPeaks_1_per_A,
const float* intensities, int threadCount);
const float* intensities, float maxRefinementDisbalance, int threadCount);
#ifdef __cplusplus
......
......@@ -36,22 +36,22 @@ void Refinement::refineVariableLattice(Lattice& lattice, const Matrix3Xf& ucsDir
Array33f gradient;
Matrix3f offsetBasis = basis;
offsetBasis(0, 0) += delta;
gradient(0, 0) = getMeanDefect(offsetBasis, ucsDirections, ucsBorderNorms) - meanDefects[i];
gradient(0, 0) = getMeanDefect(offsetBasis, ucsDirections, ucsBorderNorms, false) - meanDefects[i];
offsetBasis(0, 0) = basis(0, 0);
offsetBasis(1, 0) += delta;
gradient(1, 0) = getMeanDefect(offsetBasis, ucsDirections, ucsBorderNorms) - meanDefects[i];
gradient(1, 0) = getMeanDefect(offsetBasis, ucsDirections, ucsBorderNorms, false) - meanDefects[i];
offsetBasis(1, 0) = basis(1, 0);
offsetBasis(2, 0) += delta;
gradient(2, 0) = getMeanDefect(offsetBasis, ucsDirections, ucsBorderNorms) - meanDefects[i];
gradient(2, 0) = getMeanDefect(offsetBasis, ucsDirections, ucsBorderNorms, false) - meanDefects[i];
offsetBasis(2, 0) = basis(2, 0);
offsetBasis(0, 1) += delta;
gradient(0, 1) = getMeanDefect(offsetBasis, ucsDirections, ucsBorderNorms) - meanDefects[i];
gradient(0, 1) = getMeanDefect(offsetBasis, ucsDirections, ucsBorderNorms, false) - meanDefects[i];
offsetBasis(0, 1) = basis(0, 1);
offsetBasis(1, 1) += delta;
gradient(1, 1) = getMeanDefect(offsetBasis, ucsDirections, ucsBorderNorms) - meanDefects[i];
gradient(1, 1) = getMeanDefect(offsetBasis, ucsDirections, ucsBorderNorms, false) - meanDefects[i];
offsetBasis(1, 1) = basis(1, 1);
offsetBasis(2, 1) += delta;
gradient(2, 1) = getMeanDefect(offsetBasis, ucsDirections, ucsBorderNorms) - meanDefects[i];
gradient(2, 1) = getMeanDefect(offsetBasis, ucsDirections, ucsBorderNorms, false) - meanDefects[i];
offsetBasis(2, 1) = basis(2, 1);
offsetBasis(0, 2) += delta;
gradient(0, 2) = getMeanDefect(offsetBasis, ucsDirections, ucsBorderNorms) - meanDefects[i];
......@@ -150,15 +150,20 @@ int Refinement::getFittedPeaks(Lattice& lattice, Eigen::Array<bool, Eigen::Dynam
}
void Refinement::getDefects(ArrayXf& defects, const Matrix3f& basis, const Matrix3Xf& ucsDirections, const Matrix2Xf& ucsBorderNorms)
void Refinement::getDefects(ArrayXf& defects, const Matrix3f& basis, const Matrix3Xf& ucsDirections, const Matrix2Xf& ucsBorderNorms,
bool significantChangesToPreviousCall)
{
ucsBorderNormsSquared = ucsBorderNorms.array().square().matrix();
Matrix3f basis_inverse = basis.inverse();
millerIndices_close = basis_inverse * (ucsDirections.array().rowwise() * ucsBorderNorms.array().row(0)).matrix();
roundTowardsZero(millerIndices_close);
millerIndices_far = basis_inverse * (ucsDirections.array().rowwise() * ucsBorderNorms.array().row(1)).matrix();
roundAwayFromZero(millerIndices_far);
if (significantChangesToPreviousCall)
{
ucsBorderNormsSquared = ucsBorderNorms.array().square().matrix();
millerIndices_close = basis_inverse * (ucsDirections.array().rowwise() * ucsBorderNorms.array().row(0)).matrix();
roundTowardsZero(millerIndices_close);
millerIndices_far = basis_inverse * (ucsDirections.array().rowwise() * ucsBorderNorms.array().row(1)).matrix();
roundAwayFromZero(millerIndices_far);
}
int peakCount = ucsDirections.cols();
defects.resize(peakCount);
......@@ -182,10 +187,10 @@ void Refinement::getDefects(ArrayXf& defects, const Matrix3f& basis, const Matri
candidatePeaks.noalias() = basis * Map<Matrix3Xf>((float*)millerIndices.data(), 3, millerIndices.size());
candidatePeaksNormsSquared = candidatePeaks.colwise().squaredNorm();
// clear peaks that exceed the borders
int validCandidatePeaksCount = 0;
validCandidatePeaksCount = 0;
for (int j = 0, end = candidatePeaksNormsSquared.size(); j < end; j++)
{
if (candidatePeaksNormsSquared(j) > ucsBorderNormsSquared(0, i) && candidatePeaksNormsSquared(j) < ucsBorderNormsSquared(1, i))
if (candidatePeaksNormsSquared(j) > ucsBorderNormsSquared(0, i) & candidatePeaksNormsSquared(j) < ucsBorderNormsSquared(1, i))
{
candidatePeaks.col(validCandidatePeaksCount) = candidatePeaks.col(j);
validCandidatePeaksCount++;
......@@ -205,9 +210,9 @@ void Refinement::getDefects(ArrayXf& defects, const Matrix3f& basis, const Matri
}
}
double Refinement::getMeanDefect(const Matrix3f& basis, const Matrix3Xf& ucsDirections, const Matrix2Xf& ucsBorderNorms)
double Refinement::getMeanDefect(const Matrix3f& basis, const Matrix3Xf& ucsDirections, const Matrix2Xf& ucsBorderNorms, bool significantChangesToPreviousCall)
{
getDefects(defects, basis, ucsDirections, ucsBorderNorms);
getDefects(defects, basis, ucsDirections, ucsBorderNorms, significantChangesToPreviousCall);
notPredictablePeaks = defects > tolerance;
int16_t notPredictablePeaksCount = notPredictablePeaks.count();
......
#include "PinkIndexer.h"
#include "adaptions/crystfel/PinkIndexer.h"
static float getRefinementDisbalance(Eigen::Array<bool, Eigen::Dynamic, 1>& fittedPeaks, Eigen::Matrix3Xf reciprocalPeaks);
static float getRefinementDisbalance2(Eigen::Array<bool, Eigen::Dynamic, 1>& fittedPeaks, Eigen::Matrix3Xf reciprocalPeaks);
extern "C" PinkIndexer* PinkIndexer_new(ExperimentSettings* experimentSettings, consideredPeaksCount_t consideredPeaksCount, angleResolution_t angleResolution,
refinementType_t refinementType, float maxResolutionForIndexing_1_per_A)
......@@ -84,7 +86,7 @@ extern "C" void PinkIndexer_delete(PinkIndexer* pinkIndexer)
}
extern "C" int PinkIndexer_indexPattern(PinkIndexer* pinkIndexer, Lattice_t* indexedLattice, reciprocalPeaks_1_per_A_t* meanReciprocalPeaks_1_per_A,
const float* intensities, int threadCount)
const float* intensities, float maxRefinementDisbalance, 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++)
......@@ -109,6 +111,13 @@ extern "C" int PinkIndexer_indexPattern(PinkIndexer* pinkIndexer, Lattice_t* ind
indexedLattice->cy = basis(1, 2);
indexedLattice->cz = basis(2, 2);
float refinementDisbalance = getRefinementDisbalance(fittedPeaks, meanReciprocalPeaks_1_per_A_matrix);
float refinementDisbalance2 = getRefinementDisbalance2(fittedPeaks, meanReciprocalPeaks_1_per_A_matrix);
if (refinementDisbalance > maxRefinementDisbalance || refinementDisbalance2 > maxRefinementDisbalance * 0.5)
{
return 0;
}
int unfittedPeaksCount = 0;
for (int i = 0; i < fittedPeaks.size(); i++)
{
......@@ -124,4 +133,69 @@ extern "C" int PinkIndexer_indexPattern(PinkIndexer* pinkIndexer, Lattice_t* ind
meanReciprocalPeaks_1_per_A->peakCount = unfittedPeaksCount;
return matchedPeaksCount;
}
static float getRefinementDisbalance(Eigen::Array<bool, Eigen::Dynamic, 1>& fittedPeaks, Eigen::Matrix3Xf reciprocalPeaks)
{
int fittedCount = fittedPeaks.count();
int totalCount = fittedPeaks.size();
if (fittedCount == 0)
{
return 2;
}
Eigen::Vector2f fittedCenterOfGravity, totalCenterOfGravity;
totalCenterOfGravity.setZero();
fittedCenterOfGravity.setZero();
for (int i = 0; i < totalCount; ++i)
{
Eigen::Vector2f normalizedVector = reciprocalPeaks.col(i).head(2).normalized();
totalCenterOfGravity += normalizedVector;
if (fittedPeaks[i])
{
fittedCenterOfGravity += normalizedVector;
}
}
totalCenterOfGravity /= totalCount;
fittedCenterOfGravity /= fittedCount;
float disbalance = (totalCenterOfGravity - fittedCenterOfGravity).norm();
return disbalance;
}
static float getRefinementDisbalance2(Eigen::Array<bool, Eigen::Dynamic, 1>& fittedPeaks, Eigen::Matrix3Xf reciprocalPeaks)
{
int fittedCount = fittedPeaks.count();
int totalCount = fittedPeaks.size();
if (fittedCount == 0)
{
return 2;
}
Eigen::Matrix3Xf fittedReciprocalPeaks;
fittedReciprocalPeaks.resize(3, fittedCount);
fittedCount = 0;
for (int i = 0; i < totalCount; ++i)
{
if (fittedPeaks[i])
{
fittedReciprocalPeaks.col(fittedCount) = reciprocalPeaks.col(i);
fittedCount++;
}
}
Eigen::Vector3f totalCenterOfGravity = reciprocalPeaks.rowwise().mean();
float totalMeanDist = (reciprocalPeaks.colwise() - totalCenterOfGravity).cwiseAbs().rowwise().mean().norm();
Eigen::Vector3f fittedCenterOfGravity = fittedReciprocalPeaks.rowwise().mean();
float fittedMeanDist = (fittedReciprocalPeaks.colwise() - fittedCenterOfGravity).cwiseAbs().rowwise().mean().norm();
float disbalance = 1.0f - std::abs(fittedMeanDist / totalMeanDist);
return disbalance;
}
\ No newline at end of file
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