Skip to content

GitLab

Commit 54dad96a authored by Yaroslav Gevorkov's avatar Yaroslav Gevorkov
Browse files

Merge branch 'master' of https://stash.desy.de/scm/~gevorkov/pinkindexer

parents 0e3742ce 8281cc6d
Hide whitespace changes
Inline Side-by-side
include/Backprojection.h
View file @ 54dad96a
...@@ -11,7 +11,7 @@ namespace pinkIndexer ...@@ -11,7 +11,7 @@ namespace pinkIndexer
public: public:
Backprojection(const ExperimentSettings& experimentSettings); Backprojection(const ExperimentSettings& experimentSettings);
void backProject(const Eigen::Matrix2Xf& detectorPeaks_m, Eigen::Matrix3Xf& ucsDirections, Eigen::Array2Xf& ucsBorderNorms) const; void backProject(const Eigen::Matrix2Xf& detectorPeaks_m, Eigen::Matrix3Xf& ulsDirections, Eigen::Array2Xf& ulsBorderNorms) const;
private: private:
ExperimentSettings experimentSettings; ExperimentSettings experimentSettings;
......
include/PinkIndexer.h
View file @ 54dad96a
...@@ -58,9 +58,9 @@ namespace pinkIndexer ...@@ -58,9 +58,9 @@ namespace pinkIndexer
private: private:
// reduces number of peaks according to parameter "ConsideredPeaksCount" // reduces number of peaks according to parameter "ConsideredPeaksCount"
void reducePeakCount(Eigen::Matrix3Xf& ucsDirections, Eigen::Array2Xf& ucsBorderNorms, Eigen::RowVectorXf& intensities, void reducePeakCount(Eigen::Matrix3Xf& ulsDirections, Eigen::Array2Xf& ulsBorderNorms, Eigen::RowVectorXf& intensities,
const Eigen::Matrix2Xf& detectorPeaks_m); const Eigen::Matrix2Xf& detectorPeaks_m);
void refine(Lattice& indexedLattice, Eigen::Vector2f& centerShift, const Eigen::Matrix3Xf& ucsDirections, const Eigen::Array2Xf& ucsBorderNorms, void refine(Lattice& indexedLattice, Eigen::Vector2f& centerShift, const Eigen::Matrix3Xf& ulsDirections, const Eigen::Array2Xf& ulsBorderNorms,
const Eigen::Matrix2Xf& detectorPeaks_m, int threadCount); const Eigen::Matrix2Xf& detectorPeaks_m, int threadCount);
float getAngleResolution(); float getAngleResolution();
......
include/Refinement.h
View file @ 54dad96a
...@@ -15,16 +15,16 @@ namespace pinkIndexer ...@@ -15,16 +15,16 @@ namespace pinkIndexer
Refinement(float tolerance); Refinement(float tolerance);
Refinement(float tolerance, const Backprojection& backprojection); Refinement(float tolerance, const Backprojection& backprojection);
void refineFixedLattice(Lattice& lattice, const Eigen::Matrix3Xf& ucsDirections, const Eigen::Array2Xf& ucsBorderNorms); void refineFixedLattice(Lattice& lattice, const Eigen::Matrix3Xf& ulsDirections, const Eigen::Array2Xf& ulsBorderNorms);
void refineVariableLattice(Lattice& lattice, const Eigen::Matrix3Xf& ucsDirections, const Eigen::Array2Xf& ucsBorderNorms); void refineVariableLattice(Lattice& lattice, const Eigen::Matrix3Xf& ulsDirections, const Eigen::Array2Xf& ulsBorderNorms);
void refineVariableLatticeWithCenter(Lattice& lattice, Eigen::Vector2f& centerShift, const Eigen::Matrix2Xf& detectorPeaks_m); void refineVariableLatticeWithCenter(Lattice& lattice, Eigen::Vector2f& centerShift, const Eigen::Matrix2Xf& detectorPeaks_m);
void refineCenter(Eigen::Vector2f& centerShift, const Eigen::Matrix3f& basis, const Eigen::Matrix2Xf& detectorPeaks_m, float startStepSize = 20e-6); void refineCenter(Eigen::Vector2f& centerShift, const Eigen::Matrix3f& basis, const Eigen::Matrix2Xf& detectorPeaks_m, float startStepSize = 20e-6);
int getFittedPeaksCount(Lattice& lattice, const Eigen::Matrix3Xf& ucsDirections, const Eigen::Array2Xf& ucsBorderNorms); int getFittedPeaksCount(Lattice& lattice, const Eigen::Matrix3Xf& ulsDirections, const Eigen::Array2Xf& ulsBorderNorms);
int getFittedPeaks(Lattice& lattice, Eigen::Array<bool, Eigen::Dynamic, 1>& fittedPeaks, const Eigen::Matrix3Xf& ucsDirections, int getFittedPeaks(Lattice& lattice, Eigen::Array<bool, Eigen::Dynamic, 1>& fittedPeaks, const Eigen::Matrix3Xf& ulsDirections,
const Eigen::Array2Xf& ucsBorderNorms); const Eigen::Array2Xf& ulsBorderNorms);
double getMeanDefect(const Eigen::Matrix3f& basis, const Eigen::Matrix3Xf& ucsDirections, const Eigen::Array2Xf& ucsBorderNorms, double getMeanDefect(const Eigen::Matrix3f& basis, const Eigen::Matrix3Xf& ulsDirections, const Eigen::Array2Xf& ulsBorderNorms,
bool significantChangesToPreviousCall = true); bool significantChangesToPreviousCall = true);
void setTolerance(float tolerance) void setTolerance(float tolerance)
...@@ -37,10 +37,10 @@ namespace pinkIndexer ...@@ -37,10 +37,10 @@ namespace pinkIndexer
} }
private: private:
void getDefects(Eigen::ArrayXf& defects, const Eigen::Matrix3f& basis, const Eigen::Matrix3Xf& ucsDirections, const Eigen::Array2Xf& ucsBorderNorms, void getDefects(Eigen::ArrayXf& defects, const Eigen::Matrix3f& basis, const Eigen::Matrix3Xf& ulsDirections, const Eigen::Array2Xf& ulsBorderNorms,
bool significantChangesToPreviousCall = true); bool significantChangesToPreviousCall = true);
void getCenterShiftedBackprojection(Eigen::Matrix3Xf& ucsDirections, Eigen::Array2Xf& ucsBorderNorms, const Eigen::Matrix2Xf& detectorPeaks_m, void getCenterShiftedBackprojection(Eigen::Matrix3Xf& ulsDirections, Eigen::Array2Xf& ulsBorderNorms, const Eigen::Matrix2Xf& detectorPeaks_m,
const Eigen::Vector2f& centerShift); const Eigen::Vector2f& centerShift);
float tolerance; float tolerance;
...@@ -49,7 +49,7 @@ namespace pinkIndexer ...@@ -49,7 +49,7 @@ namespace pinkIndexer
// preallocation // preallocation
Eigen::ArrayXf defects; Eigen::ArrayXf defects;
Eigen::Array2Xf ucsBorderNormsSquared; Eigen::Array2Xf ulsBorderNormsSquared;
Eigen::Matrix3Xf millerIndices_close; Eigen::Matrix3Xf millerIndices_close;
Eigen::Matrix3Xf millerIndices_far; Eigen::Matrix3Xf millerIndices_far;
Eigen::Matrix3Xf candidatePeaks; Eigen::Matrix3Xf candidatePeaks;
...@@ -63,7 +63,7 @@ namespace pinkIndexer ...@@ -63,7 +63,7 @@ namespace pinkIndexer
Eigen::ArrayXf meanDefects_centerAdjustment; Eigen::ArrayXf meanDefects_centerAdjustment;
int validCandidatePeaksCount; int validCandidatePeaksCount;
Eigen::Matrix2Xf detectorPeaks_m_shifted; Eigen::Matrix2Xf detectorPeaks_m_shifted;
Eigen::Matrix3Xf ucsDirections; Eigen::Matrix3Xf ulsDirections;
Eigen::Array2Xf ucsBorderNorms; Eigen::Array2Xf ulsBorderNorms;
}; };
} // namespace pinkIndexer } // namespace pinkIndexer
\ No newline at end of file
include/Sinogram.h
View file @ 54dad96a
...@@ -18,9 +18,9 @@ namespace pinkIndexer ...@@ -18,9 +18,9 @@ namespace pinkIndexer
void setSinogramAngleResolution(float angleResolution_deg); void setSinogramAngleResolution(float angleResolution_deg);
void computeSinogram(const Eigen::Matrix3Xf& ucsDirections, const Eigen::Matrix2Xf ucsBorderNorms); void computeSinogram(const Eigen::Matrix3Xf& ulsDirections, const Eigen::Matrix2Xf ulsBorderNorms);
void computeSinogramParallel(const Eigen::Matrix3Xf& ucsDirections, const Eigen::Matrix2Xf ucsBorderNorms, int slaveThreadCount); void computeSinogramParallel(const Eigen::Matrix3Xf& ulsDirections, const Eigen::Matrix2Xf ulsBorderNorms, int slaveThreadCount);
void computeSinogramParallel2(const Eigen::Matrix3Xf& ucsDirections, const Eigen::Matrix2Xf ucsBorderNorms, int slaveThreadCount); void computeSinogramParallel2(const Eigen::Matrix3Xf& ulsDirections, const Eigen::Matrix2Xf ulsBorderNorms, int slaveThreadCount);
void getBestRotation(Eigen::AngleAxisf& bestRotation); void getBestRotation(Eigen::AngleAxisf& bestRotation);
......
src/Backprojection.cpp
View file @ 54dad96a
...@@ -14,7 +14,7 @@ namespace pinkIndexer ...@@ -14,7 +14,7 @@ namespace pinkIndexer
{ {
} }
void Backprojection::backProject(const Matrix2Xf& detectorPeaks_m, Matrix3Xf& ucsDirections, Array2Xf& ucsBorderNorms) const void Backprojection::backProject(const Matrix2Xf& detectorPeaks_m, Matrix3Xf& ulsDirections, Array2Xf& ulsBorderNorms) const
{ {
Matrix3Xf projectionDirections(3, detectorPeaks_m.cols()); Matrix3Xf projectionDirections(3, detectorPeaks_m.cols());
projectionDirections << RowVectorXf::Constant(1, detectorPeaks_m.cols(), experimentSettings.getDetectorDistance_m()), detectorPeaks_m.row(0), projectionDirections << RowVectorXf::Constant(1, detectorPeaks_m.cols(), experimentSettings.getDetectorDistance_m()), detectorPeaks_m.row(0),
...@@ -26,11 +26,11 @@ namespace pinkIndexer ...@@ -26,11 +26,11 @@ namespace pinkIndexer
Matrix3Xf backprojectedPointsFar = projectionDirections * experimentSettings.getReciprocalLambdaShort_1A(); Matrix3Xf backprojectedPointsFar = projectionDirections * experimentSettings.getReciprocalLambdaShort_1A();
backprojectedPointsFar.row(0) -= RowVectorXf::Constant(backprojectedPointsFar.cols(), experimentSettings.getReciprocalLambdaShort_1A()); backprojectedPointsFar.row(0) -= RowVectorXf::Constant(backprojectedPointsFar.cols(), experimentSettings.getReciprocalLambdaShort_1A());
ucsBorderNorms.resize(2, detectorPeaks_m.cols()); ulsBorderNorms.resize(2, detectorPeaks_m.cols());
ucsBorderNorms.row(0) = backprojectedPointsClose.colwise().norm().array() - experimentSettings.getReflectionRadius(); ulsBorderNorms.row(0) = backprojectedPointsClose.colwise().norm().array() - experimentSettings.getReflectionRadius();
ucsBorderNorms.row(1) = backprojectedPointsFar.colwise().norm().array() + experimentSettings.getReflectionRadius(); ulsBorderNorms.row(1) = backprojectedPointsFar.colwise().norm().array() + experimentSettings.getReflectionRadius();
ucsDirections = backprojectedPointsFar.colwise().normalized(); ulsDirections = backprojectedPointsFar.colwise().normalized();
} }
} // namespace pinkIndexer } // namespace pinkIndexer
\ No newline at end of file
src/PinkIndexer.cpp
View file @ 54dad96a
...@@ -48,26 +48,26 @@ namespace pinkIndexer ...@@ -48,26 +48,26 @@ namespace pinkIndexer
if (detectorPeaks_m.cols() < 2) if (detectorPeaks_m.cols() < 2)
return 0; return 0;
Matrix3Xf ucsDirections; // Unit vector of ULS Matrix3Xf ulsDirections; // Unit vector of ULS
Array2Xf ucsBorderNorms; // 2 borders between intersection with the 2 Ewald spheres Array2Xf ulsBorderNorms; // 2 borders between intersection with the 2 Ewald spheres
backprojection.backProject(detectorPeaks_m, ucsDirections, ucsBorderNorms); backprojection.backProject(detectorPeaks_m, ulsDirections, ulsBorderNorms);
Matrix3Xf ucsDirections_reduced = ucsDirections; Matrix3Xf ulsDirections_reduced = ulsDirections;
Array2Xf ucsBorderNorms_reduced = ucsBorderNorms; Array2Xf ulsBorderNorms_reduced = ulsBorderNorms;
reducePeakCount(ucsDirections_reduced, ucsBorderNorms_reduced, intensities, detectorPeaks_m); reducePeakCount(ulsDirections_reduced, ulsBorderNorms_reduced, intensities, detectorPeaks_m);
if (threadCount == 1) //threads from crystfel (--pinkIndexer_thread_count) if (threadCount == 1) //threads from crystfel (--pinkIndexer_thread_count)
{ {
sinogram.computeSinogram(ucsDirections_reduced, ucsBorderNorms_reduced); sinogram.computeSinogram(ulsDirections_reduced, ulsBorderNorms_reduced);
} }
else if (threadCount < 32) else if (threadCount < 32)
{ {
sinogram.computeSinogramParallel(ucsDirections_reduced, ucsBorderNorms_reduced, threadCount); sinogram.computeSinogramParallel(ulsDirections_reduced, ulsBorderNorms_reduced, threadCount);
} }
else // probably not ran in parellel to other instances else // probably not ran in parellel to other instances
{ {
int slaveThreadCount = threadCount - 1; int slaveThreadCount = threadCount - 1;
sinogram.computeSinogramParallel2(ucsDirections_reduced, ucsBorderNorms_reduced, slaveThreadCount); sinogram.computeSinogramParallel2(ulsDirections_reduced, ulsBorderNorms_reduced, slaveThreadCount);
} }
AngleAxisf bestRotation; AngleAxisf bestRotation;
...@@ -76,18 +76,18 @@ namespace pinkIndexer ...@@ -76,18 +76,18 @@ namespace pinkIndexer
Matrix3f bestBasis = bestRotation * sampleLattice.getBasis(); Matrix3f bestBasis = bestRotation * sampleLattice.getBasis();
indexedLattice = Lattice(bestBasis); indexedLattice = Lattice(bestBasis);
refine(indexedLattice, centerShift, ucsDirections, ucsBorderNorms, detectorPeaks_m, threadCount); refine(indexedLattice, centerShift, ulsDirections, ulsBorderNorms, detectorPeaks_m, threadCount);
indexedLattice.minimize(); indexedLattice.minimize();
indexedLattice.reorder(sampleLattice); //reorder lattice vectors to match sample lattice (from crystfel) indexedLattice.reorder(sampleLattice); //reorder lattice vectors to match sample lattice (from crystfel)
// return how many peaks have been fitted correctly (according to tolerance) // return how many peaks have been fitted correctly (according to tolerance)
return refinement.getFittedPeaks(indexedLattice, fittedPeaks, ucsDirections, ucsBorderNorms); return refinement.getFittedPeaks(indexedLattice, fittedPeaks, ulsDirections, ulsBorderNorms);
// sinogram.saveToFile("C:\\DesyFiles\\workspaces\\VisualStudio_workspace\\pinkIndexer\\workfolder\\sinogram"); // sinogram.saveToFile("C:\\DesyFiles\\workspaces\\VisualStudio_workspace\\pinkIndexer\\workfolder\\sinogram");
} }
void PinkIndexer::refine(Lattice& indexedLattice, Vector2f& centerShift, const Matrix3Xf& ucsDirections, const Array2Xf& ucsBorderNorms, void PinkIndexer::refine(Lattice& indexedLattice, Vector2f& centerShift, const Matrix3Xf& ulsDirections, const Array2Xf& ulsBorderNorms,
const Matrix2Xf& detectorPeaks_m, int threadCount) const Matrix2Xf& detectorPeaks_m, int threadCount)
{ {
centerShift.setZero(); centerShift.setZero();
...@@ -98,32 +98,32 @@ namespace pinkIndexer ...@@ -98,32 +98,32 @@ namespace pinkIndexer
break; break;
case RefinementType::fixedLatticeParameters: 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.refineFixedLattice(indexedLattice, ulsDirections, ulsBorderNorms);
refinement.setTolerance(finalRefinementTolerance); refinement.setTolerance(finalRefinementTolerance);
refinement.refineFixedLattice(indexedLattice, ucsDirections, ucsBorderNorms); refinement.refineFixedLattice(indexedLattice, ulsDirections, ulsBorderNorms);
break; break;
case RefinementType::variableLatticeParameters: case RefinementType::variableLatticeParameters:
refinement.setTolerance(min(finalRefinementTolerance * 2.0, 0.12)); refinement.setTolerance(min(finalRefinementTolerance * 2.0, 0.12));
refinement.refineVariableLattice(indexedLattice, ucsDirections, ucsBorderNorms); refinement.refineVariableLattice(indexedLattice, ulsDirections, ulsBorderNorms);
refinement.setTolerance(finalRefinementTolerance); refinement.setTolerance(finalRefinementTolerance);
refinement.refineVariableLattice(indexedLattice, ucsDirections, ucsBorderNorms); refinement.refineVariableLattice(indexedLattice, ulsDirections, ulsBorderNorms);
break; break;
case RefinementType::firstFixedThenVariableLatticeParameters: case RefinementType::firstFixedThenVariableLatticeParameters:
refinement.setTolerance(min(finalRefinementTolerance * 2.5, 0.12)); refinement.setTolerance(min(finalRefinementTolerance * 2.5, 0.12));
refinement.refineFixedLattice(indexedLattice, ucsDirections, ucsBorderNorms); refinement.refineFixedLattice(indexedLattice, ulsDirections, ulsBorderNorms);
refinement.setTolerance(min(finalRefinementTolerance * 1.8, 0.10)); refinement.setTolerance(min(finalRefinementTolerance * 1.8, 0.10));
refinement.refineVariableLattice(indexedLattice, ucsDirections, ucsBorderNorms); refinement.refineVariableLattice(indexedLattice, ulsDirections, ulsBorderNorms);
refinement.setTolerance(finalRefinementTolerance); refinement.setTolerance(finalRefinementTolerance);
refinement.refineVariableLattice(indexedLattice, ucsDirections, ucsBorderNorms); refinement.refineVariableLattice(indexedLattice, ulsDirections, ulsBorderNorms);
break; break;
case RefinementType::firstFixedThenVariableLatticeParametersMultiSeed: case RefinementType::firstFixedThenVariableLatticeParametersMultiSeed:
{ {
constexpr int refinementTries = 1000; constexpr int refinementTries = 1000; //TODO: increase if better refinement needed
int fittedNodesCount[refinementTries]; int fittedNodesCount[refinementTries];
double fittedNodesMeanDefects[refinementTries]; double fittedNodesMeanDefects[refinementTries];
Lattice fittedLattices[refinementTries]; Lattice fittedLattices[refinementTries];
float maxRelativeDeviation = 0.0125; float maxRelativeDeviation = 0.0125; //TODO: Parameter??
Array<float, 1, 3> columnDeviationNorms = indexedLattice.getBasis().colwise().norm() * maxRelativeDeviation; Array<float, 1, 3> columnDeviationNorms = indexedLattice.getBasis().colwise().norm() * maxRelativeDeviation;
#ifdef OPENMP_AVAILABLE #ifdef OPENMP_AVAILABLE
...@@ -146,14 +146,14 @@ namespace pinkIndexer ...@@ -146,14 +146,14 @@ namespace pinkIndexer
fittedLattices[i] = Lattice(currentBasis); fittedLattices[i] = Lattice(currentBasis);
refinement.setTolerance(min(finalRefinementTolerance * 2.5, 0.12)); refinement.setTolerance(min(finalRefinementTolerance * 2.5, 0.12));
refinement.refineFixedLattice(fittedLattices[i], ucsDirections, ucsBorderNorms); refinement.refineFixedLattice(fittedLattices[i], ulsDirections, ulsBorderNorms);
refinement.setTolerance(min(finalRefinementTolerance * 1.8, 0.10)); refinement.setTolerance(min(finalRefinementTolerance * 1.8, 0.10));
refinement.refineVariableLattice(fittedLattices[i], ucsDirections, ucsBorderNorms); refinement.refineVariableLattice(fittedLattices[i], ulsDirections, ulsBorderNorms);
refinement.setTolerance(finalRefinementTolerance); refinement.setTolerance(finalRefinementTolerance);
refinement.refineVariableLattice(fittedLattices[i], ucsDirections, ucsBorderNorms); refinement.refineVariableLattice(fittedLattices[i], ulsDirections, ulsBorderNorms);
fittedNodesCount[i] = refinement.getFittedPeaksCount(fittedLattices[i], ucsDirections, ucsBorderNorms); fittedNodesCount[i] = refinement.getFittedPeaksCount(fittedLattices[i], ulsDirections, ulsBorderNorms);
fittedNodesMeanDefects[i] = refinement.getMeanDefect(fittedLattices[i].getBasis(), ucsDirections, ucsBorderNorms); fittedNodesMeanDefects[i] = refinement.getMeanDefect(fittedLattices[i].getBasis(), ulsDirections, ulsBorderNorms);
} }
int maxFittedNodesCount = 1; int maxFittedNodesCount = 1;
...@@ -164,7 +164,7 @@ namespace pinkIndexer ...@@ -164,7 +164,7 @@ namespace pinkIndexer
// (fittedNodesCount[i] == maxFittedNodesCount && fittedNodesMeanDefects[i] < minFittedNodesMeanDefect)) // (fittedNodesCount[i] == maxFittedNodesCount && fittedNodesMeanDefects[i] < minFittedNodesMeanDefect))
if (fittedNodesCount[i] / maxFittedNodesCount > 1.2 || if (fittedNodesCount[i] / maxFittedNodesCount > 1.2 ||
((fittedNodesCount[i] > maxFittedNodesCount || fittedNodesMeanDefects[i] < minFittedNodesMeanDefect) && ((fittedNodesCount[i] > maxFittedNodesCount || fittedNodesMeanDefects[i] < minFittedNodesMeanDefect) &&
(fittedNodesCount[i] / maxFittedNodesCount - 1) * 3 > (minFittedNodesMeanDefect / fittedNodesMeanDefects[i] - 1))) (fittedNodesCount[i] / maxFittedNodesCount - 1) * 3 > (minFittedNodesMeanDefect / fittedNodesMeanDefects[i] - 1))) //TODO: can be optimized
{ {
maxFittedNodesCount = fittedNodesCount[i]; maxFittedNodesCount = fittedNodesCount[i];
minFittedNodesMeanDefect = fittedNodesMeanDefects[i]; minFittedNodesMeanDefect = fittedNodesMeanDefects[i];
...@@ -175,7 +175,7 @@ namespace pinkIndexer ...@@ -175,7 +175,7 @@ namespace pinkIndexer
break; break;
case RefinementType::firstFixedThenVariableLatticeParametersCenterAdjustmentMultiSeed: case RefinementType::firstFixedThenVariableLatticeParametersCenterAdjustmentMultiSeed:
{ {
constexpr int refinementTries = 2000; constexpr int refinementTries = 2000; //TODO: increase if better refinement needed
int fittedNodesCount[refinementTries]; int fittedNodesCount[refinementTries];
double fittedNodesMeanDefects[refinementTries]; double fittedNodesMeanDefects[refinementTries];
Lattice fittedLattices[refinementTries]; Lattice fittedLattices[refinementTries];
...@@ -203,17 +203,17 @@ namespace pinkIndexer ...@@ -203,17 +203,17 @@ namespace pinkIndexer
currentBasis = indexedLattice.getBasis(); currentBasis = indexedLattice.getBasis();
} }
fittedLattices[i] = Lattice(currentBasis); fittedLattices[i] = Lattice(currentBasis);
centerShifts[i] = Vector2f::Random() * 80e-6; centerShifts[i] = Vector2f::Random() * 80e-6; //TODO: hardcoded JUNGFRAU pixel size. Can be improved (maybe parameter??)
refinement.setTolerance(min(finalRefinementTolerance * 2.5, 0.12)); refinement.setTolerance(min(finalRefinementTolerance * 2.5, 0.12));
refinement.refineFixedLattice(fittedLattices[i], ucsDirections, ucsBorderNorms); refinement.refineFixedLattice(fittedLattices[i], ulsDirections, ulsBorderNorms);
refinement.setTolerance(min(finalRefinementTolerance * 1.8, 0.10)); refinement.setTolerance(min(finalRefinementTolerance * 1.8, 0.10));
refinement.refineVariableLatticeWithCenter(fittedLattices[i], centerShifts[i], detectorPeaks_m); refinement.refineVariableLatticeWithCenter(fittedLattices[i], centerShifts[i], detectorPeaks_m);
refinement.setTolerance(finalRefinementTolerance); refinement.setTolerance(finalRefinementTolerance);
refinement.refineVariableLatticeWithCenter(fittedLattices[i], centerShifts[i], detectorPeaks_m); refinement.refineVariableLatticeWithCenter(fittedLattices[i], centerShifts[i], detectorPeaks_m);
fittedNodesCount[i] = refinement.getFittedPeaksCount(fittedLattices[i], ucsDirections, ucsBorderNorms); fittedNodesCount[i] = refinement.getFittedPeaksCount(fittedLattices[i], ulsDirections, ulsBorderNorms);
fittedNodesMeanDefects[i] = refinement.getMeanDefect(fittedLattices[i].getBasis(), ucsDirections, ucsBorderNorms); fittedNodesMeanDefects[i] = refinement.getMeanDefect(fittedLattices[i].getBasis(), ulsDirections, ulsBorderNorms);
} }
int maxFittedNodesCount = 1; int maxFittedNodesCount = 1;
...@@ -239,7 +239,7 @@ namespace pinkIndexer ...@@ -239,7 +239,7 @@ namespace pinkIndexer
} }
} }
void PinkIndexer::reducePeakCount(Matrix3Xf& ucsDirections, Array2Xf& ucsBorderNorms, RowVectorXf& intensities, const Eigen::Matrix2Xf& detectorPeaks_m) void PinkIndexer::reducePeakCount(Matrix3Xf& ulsDirections, Array2Xf& ulsBorderNorms, RowVectorXf& intensities, const Eigen::Matrix2Xf& detectorPeaks_m)
{ {
Matrix2Xf detectorPeaks_m_copy = detectorPeaks_m; Matrix2Xf detectorPeaks_m_copy = detectorPeaks_m;
...@@ -247,17 +247,17 @@ namespace pinkIndexer ...@@ -247,17 +247,17 @@ namespace pinkIndexer
int peaksKept_res = 0; int peaksKept_res = 0;
for (int i = 0; i < intensities.size(); i++) for (int i = 0; i < intensities.size(); i++)
{ {
if (ucsBorderNorms(1, i) <= maxResolutionForIndexing_1_per_A) if (ulsBorderNorms(1, i) <= maxResolutionForIndexing_1_per_A)
{ {
ucsDirections.col(peaksKept_res) = ucsDirections.col(i); ulsDirections.col(peaksKept_res) = ulsDirections.col(i);
ucsBorderNorms.col(peaksKept_res) = ucsBorderNorms.col(i); ulsBorderNorms.col(peaksKept_res) = ulsBorderNorms.col(i);
intensities[peaksKept_res] = intensities[i]; intensities[peaksKept_res] = intensities[i];
detectorPeaks_m_copy.col(peaksKept_res) = detectorPeaks_m_copy.col(i); detectorPeaks_m_copy.col(peaksKept_res) = detectorPeaks_m_copy.col(i);
peaksKept_res++; peaksKept_res++;
} }
} }
ucsDirections.conservativeResize(NoChange, peaksKept_res); ulsDirections.conservativeResize(NoChange, peaksKept_res);
ucsBorderNorms.conservativeResize(NoChange, peaksKept_res); ulsBorderNorms.conservativeResize(NoChange, peaksKept_res);
intensities.conservativeResize(peaksKept_res); intensities.conservativeResize(peaksKept_res);
detectorPeaks_m_copy.conservativeResize(NoChange, peaksKept_res); detectorPeaks_m_copy.conservativeResize(NoChange, peaksKept_res);
...@@ -283,8 +283,8 @@ namespace pinkIndexer ...@@ -283,8 +283,8 @@ namespace pinkIndexer
{ {
if (norms[i] >= minNorm && norms[i] <= maxNorm) if (norms[i] >= minNorm && norms[i] <= maxNorm)
{ {
ucsDirections.col(peaksKept_norms) = ucsDirections.col(i); ulsDirections.col(peaksKept_norms) = ulsDirections.col(i);
ucsBorderNorms.col(peaksKept_norms) = ucsBorderNorms.col(i); ulsBorderNorms.col(peaksKept_norms) = ulsBorderNorms.col(i);
intensities[peaksKept_norms] = intensities[i]; intensities[peaksKept_norms] = intensities[i];
peaksKept_norms++; peaksKept_norms++;
} }
...@@ -300,8 +300,8 @@ namespace pinkIndexer ...@@ -300,8 +300,8 @@ namespace pinkIndexer
{ {
if (intensities[i] > minIntensity) if (intensities[i] > minIntensity)
{ {
ucsDirections.col(peaksKept_normsIntensities) = ucsDirections.col(i); ulsDirections.col(peaksKept_normsIntensities) = ulsDirections.col(i);
ucsBorderNorms.col(peaksKept_normsIntensities) = ucsBorderNorms.col(i); ulsBorderNorms.col(peaksKept_normsIntensities) = ulsBorderNorms.col(i);
intensities[peaksKept_normsIntensities] = intensities[i]; intensities[peaksKept_normsIntensities] = intensities[i];
peaksKept_normsIntensities++; peaksKept_normsIntensities++;
} }
...@@ -313,16 +313,16 @@ namespace pinkIndexer ...@@ -313,16 +313,16 @@ namespace pinkIndexer
{ {
if (intensities[i] == minIntensity) if (intensities[i] == minIntensity)
{ {
ucsDirections.col(peaksKept_normsIntensities) = ucsDirections.col(i); ulsDirections.col(peaksKept_normsIntensities) = ulsDirections.col(i);
ucsBorderNorms.col(peaksKept_normsIntensities) = ucsBorderNorms.col(i); ulsBorderNorms.col(peaksKept_normsIntensities) = ulsBorderNorms.col(i);
intensities[peaksKept_normsIntensities] = intensities[i]; intensities[peaksKept_normsIntensities] = intensities[i];
peaksKept_normsIntensities++; peaksKept_normsIntensities++;
} }
} }
} }
ucsDirections.conservativeResize(NoChange, peaksKept_normsIntensities); ulsDirections.conservativeResize(NoChange, peaksKept_normsIntensities);
ucsBorderNorms.conservativeResize(NoChange, peaksKept_normsIntensities); ulsBorderNorms.conservativeResize(NoChange, peaksKept_normsIntensities);
intensities.conservativeResize(peaksKept_normsIntensities); intensities.conservativeResize(peaksKept_normsIntensities);
} }
......
src/Refinement.cpp
View file @ 54dad96a
...@@ -29,48 +29,48 @@ namespace pinkIndexer ...@@ -29,48 +29,48 @@ namespace pinkIndexer
millerIndices.reserve(500); millerIndices.reserve(500);
} }
void Refinement::refineVariableLattice(Lattice& lattice, const Matrix3Xf& ucsDirections, const Array2Xf& ucsBorderNorms) void Refinement::refineVariableLattice(Lattice& lattice, const Matrix3Xf& ulsDirections, const Array2Xf& ulsBorderNorms)
{ {
Matrix3f basis = lattice.getBasis(); Matrix3f basis = lattice.getBasis();
float delta = 1e-8; float delta = 1e-8; //for numerical differentiation, in A^-1
float stepSize = lattice.getBasisVectorNorms().maxCoeff() * 0.002; float stepSize = lattice.getBasisVectorNorms().maxCoeff() * 0.002;
float minStepSize = lattice.getBasisVectorNorms().minCoeff() * 0.00001; float minStepSize = lattice.getBasisVectorNorms().minCoeff() * 0.00001;
int maxStepsCount = 200; int maxStepsCount = 200;
meanDefects.resize(maxStepsCount); meanDefects.resize(maxStepsCount);
meanDefects[0] = getMeanDefect(basis, ucsDirections, ucsBorderNorms); meanDefects[0] = getMeanDefect(basis, ulsDirections, ulsBorderNorms);
for (int i = 0; i < maxStepsCount; i++) for (int i = 0; i < maxStepsCount; i++)
{ {
// cout << meanDefects[i] << endl; // cout << meanDefects[i] << endl;
Array33f gradient; Array33f gradient; //gradient for change of each basis matrix element
Matrix3f offsetBasis = basis; Matrix3f offsetBasis = basis;
offsetBasis(0, 0) += delta; offsetBasis(0, 0) += delta;
gradient(0, 0) = getMeanDefect(offsetBasis, ucsDirections, ucsBorderNorms, false) - meanDefects[i]; gradient(0, 0) = getMeanDefect(offsetBasis, ulsDirections, ulsBorderNorms, false) - meanDefects[i];
offsetBasis(0, 0) = basis(0, 0); offsetBasis(0, 0) = basis(0, 0);
offsetBasis(1, 0) += delta; offsetBasis(1, 0) += delta;
gradient(1, 0) = getMeanDefect(offsetBasis, ucsDirections, ucsBorderNorms, false) - meanDefects[i]; gradient(1, 0) = getMeanDefect(offsetBasis, ulsDirections, ulsBorderNorms, false) - meanDefects[i];
offsetBasis(1, 0) = basis(1, 0); offsetBasis(1, 0) = basis(1, 0);
offsetBasis(2, 0) += delta; offsetBasis(2, 0) += delta;
gradient(2, 0) = getMeanDefect(offsetBasis, ucsDirections, ucsBorderNorms, false) - meanDefects[i]; gradient(2, 0) = getMeanDefect(offsetBasis, ulsDirections, ulsBorderNorms, false) - meanDefects[i];
offsetBasis(2, 0) = basis(2, 0); offsetBasis(2, 0) = basis(2, 0);
offsetBasis(0, 1) += delta; offsetBasis(0, 1) += delta;
gradient(0, 1) = getMeanDefect(offsetBasis, ucsDirections, ucsBorderNorms, false) - meanDefects[i]; gradient(0, 1) = getMeanDefect(offsetBasis, ulsDirections, ulsBorderNorms, false) - meanDefects[i];
offsetBasis(0, 1) = basis(0, 1); offsetBasis(0, 1) = basis(0, 1);
offsetBasis(1, 1) += delta; offsetBasis(1, 1) += delta;
gradient(1, 1) = getMeanDefect(offsetBasis, ucsDirections, ucsBorderNorms, false) - meanDefects[i]; gradient(1, 1) = getMeanDefect(offsetBasis, ulsDirections, ulsBorderNorms, false) - meanDefects[i];
offsetBasis(1, 1) = basis(1, 1); offsetBasis(1, 1) = basis(1, 1);
offsetBasis(2, 1) += delta; offsetBasis(2, 1) += delta;
gradient(2, 1) = getMeanDefect(offsetBasis, ucsDirections, ucsBorderNorms, false) - meanDefects[i]; gradient(2, 1) = getMeanDefect(offsetBasis, ulsDirections, ulsBorderNorms, false) - meanDefects[i];
offsetBasis(2, 1) = basis(2, 1); offsetBasis(2, 1) = basis(2, 1);
offsetBasis(0, 2) += delta; offsetBasis(0, 2) += delta;
gradient(0, 2) = getMeanDefect(offsetBasis, ucsDirections, ucsBorderNorms) - meanDefects[i]; gradient(0, 2) = getMeanDefect(offsetBasis, ulsDirections, ulsBorderNorms) - meanDefects[i];
offsetBasis(0, 2) = basis(0, 2); offsetBasis(0, 2) = basis(0, 2);
offsetBasis(1, 2) += delta; offsetBasis(1, 2) += delta;
gradient(1, 2) = getMeanDefect(offsetBasis, ucsDirections, ucsBorderNorms) - meanDefects[i]; gradient(1, 2) = getMeanDefect(offsetBasis, ulsDirections, ulsBorderNorms) - meanDefects[i];
offsetBasis(1, 2) = basis(1, 2); offsetBasis(1, 2) = basis(1, 2);
offsetBasis(2, 2) += delta; offsetBasis(2, 2) += delta;
gradient(2, 2) = getMeanDefect(offsetBasis, ucsDirections, ucsBorderNorms) - meanDefects[i]; gradient(2, 2) = getMeanDefect(offsetBasis, ulsDirections, ulsBorderNorms) - meanDefects[i];
float norm = gradient.matrix().norm(); float norm = gradient.matrix().norm();
gradient = gradient / norm * stepSize; gradient = gradient / norm * stepSize;
...@@ -81,7 +81,7 @@ namespace pinkIndexer ...@@ -81,7 +81,7 @@ namespace pinkIndexer
}