Commit 47c99f8a authored by Yaroslav Gevorkov's avatar Yaroslav Gevorkov
Browse files

documentation (with Oleksandr)

parent 917dd159
...@@ -49,12 +49,15 @@ namespace pinkIndexer ...@@ -49,12 +49,15 @@ namespace pinkIndexer
PinkIndexer(const ExperimentSettings& experimentSettings, ConsideredPeaksCount consideredPeaksCount, AngleResolution angleResolution, PinkIndexer(const ExperimentSettings& experimentSettings, ConsideredPeaksCount consideredPeaksCount, AngleResolution angleResolution,
RefinementType refinementType, float maxResolutionForIndexing_1_per_A); RefinementType refinementType, float maxResolutionForIndexing_1_per_A);
// main function
int indexPattern(Lattice& indexedLattice, Eigen::Vector2f& centerShift, Eigen::Array<bool, Eigen::Dynamic, 1>& fittedPeaks, int indexPattern(Lattice& indexedLattice, Eigen::Vector2f& centerShift, Eigen::Array<bool, Eigen::Dynamic, 1>& fittedPeaks,
Eigen::RowVectorXf& intensities, const Eigen::Matrix2Xf& detectorPeaks_m, int threadCount); Eigen::RowVectorXf& intensities, const Eigen::Matrix2Xf& detectorPeaks_m, int threadCount);
int indexPattern(Lattice& indexedLattice, Eigen::Vector2f& centerShift, Eigen::Array<bool, Eigen::Dynamic, 1>& fittedPeaks, // for crystfel
int indexPattern(Lattice& indexedLattice, Eigen::Vector2f& centerShift, Eigen::Array<bool, Eigen::Dynamic, 1>& fittedPeaks,
Eigen::RowVectorXf& intensities, const Eigen::Matrix3Xf& meanReciprocalPeaks_1_per_A, int threadCount); Eigen::RowVectorXf& intensities, const Eigen::Matrix3Xf& meanReciprocalPeaks_1_per_A, int threadCount);
private: private:
// reduces number of peaks according to parameter "ConsideredPeaksCount"
void reducePeakCount(Eigen::Matrix3Xf& ucsDirections, Eigen::Array2Xf& ucsBorderNorms, Eigen::RowVectorXf& intensities, void reducePeakCount(Eigen::Matrix3Xf& ucsDirections, Eigen::Array2Xf& ucsBorderNorms, 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& ucsDirections, const Eigen::Array2Xf& ucsBorderNorms,
...@@ -73,6 +76,6 @@ namespace pinkIndexer ...@@ -73,6 +76,6 @@ namespace pinkIndexer
AngleResolution angleResolution; AngleResolution angleResolution;
RefinementType refinementType; RefinementType refinementType;
float maxResolutionForIndexing_1_per_A; float maxResolutionForIndexing_1_per_A;
float finalRefinementTolerance; float finalRefinementTolerance; //from command line in crystfel
}; };
} // namespace pinkIndexer } // namespace pinkIndexer
\ No newline at end of file
...@@ -48,15 +48,15 @@ namespace pinkIndexer ...@@ -48,15 +48,15 @@ namespace pinkIndexer
if (detectorPeaks_m.cols() < 2) if (detectorPeaks_m.cols() < 2)
return 0; return 0;
Matrix3Xf ucsDirections; Matrix3Xf ucsDirections; // Unit vector of ULS
Array2Xf ucsBorderNorms; Array2Xf ucsBorderNorms; // 2 borders between intersection with the 2 Ewald spheres
backprojection.backProject(detectorPeaks_m, ucsDirections, ucsBorderNorms); backprojection.backProject(detectorPeaks_m, ucsDirections, ucsBorderNorms);
Matrix3Xf ucsDirections_reduced = ucsDirections; Matrix3Xf ucsDirections_reduced = ucsDirections;
Array2Xf ucsBorderNorms_reduced = ucsBorderNorms; Array2Xf ucsBorderNorms_reduced = ucsBorderNorms;
reducePeakCount(ucsDirections_reduced, ucsBorderNorms_reduced, intensities, detectorPeaks_m); reducePeakCount(ucsDirections_reduced, ucsBorderNorms_reduced, intensities, detectorPeaks_m);
if (threadCount == 1) if (threadCount == 1) //threads from crystfel (--pinkIndexer_thread_count)
{ {
sinogram.computeSinogram(ucsDirections_reduced, ucsBorderNorms_reduced); sinogram.computeSinogram(ucsDirections_reduced, ucsBorderNorms_reduced);
} }
...@@ -79,8 +79,9 @@ namespace pinkIndexer ...@@ -79,8 +79,9 @@ namespace pinkIndexer
refine(indexedLattice, centerShift, ucsDirections, ucsBorderNorms, detectorPeaks_m, threadCount); refine(indexedLattice, centerShift, ucsDirections, ucsBorderNorms, detectorPeaks_m, threadCount);
indexedLattice.minimize(); indexedLattice.minimize();
indexedLattice.reorder(sampleLattice); indexedLattice.reorder(sampleLattice); //reorder lattice vectors to match sample lattice (from crystfel)
// return how many peaks have been fitted correctly (according to tolerance)
return refinement.getFittedPeaks(indexedLattice, fittedPeaks, ucsDirections, ucsBorderNorms); return refinement.getFittedPeaks(indexedLattice, fittedPeaks, ucsDirections, ucsBorderNorms);
// sinogram.saveToFile("C:\\DesyFiles\\workspaces\\VisualStudio_workspace\\pinkIndexer\\workfolder\\sinogram"); // sinogram.saveToFile("C:\\DesyFiles\\workspaces\\VisualStudio_workspace\\pinkIndexer\\workfolder\\sinogram");
...@@ -327,6 +328,7 @@ namespace pinkIndexer ...@@ -327,6 +328,7 @@ namespace pinkIndexer
float PinkIndexer::getAngleResolution() float PinkIndexer::getAngleResolution()
{ {
// angluar step in sinogram (in degree)
switch (angleResolution) switch (angleResolution)
{ {
case AngleResolution::extremelyLoose: case AngleResolution::extremelyLoose:
......
...@@ -74,15 +74,15 @@ namespace pinkIndexer ...@@ -74,15 +74,15 @@ namespace pinkIndexer
{ {
sinogram.setZero(); sinogram.setZero();
EigenSTL::vector_Matrix3Xf candidateReflectionsDirections; EigenSTL::vector_Matrix3Xf candidateReflectionsDirections; //unit vectors to all candidate reflections
reflectionsInRangeFinder.getReflectionsInRanges(candidateReflectionsDirections, ucsBorderNorms); reflectionsInRangeFinder.getReflectionsInRanges(candidateReflectionsDirections, ucsBorderNorms); //get candidate reflections
Matrix3Xf n(3, anglesCount); Matrix3Xf n(3, anglesCount); //resulting rotation axes
Array<float, 1, Eigen::Dynamic> gammah(anglesCount); // gamma half Array<float, 1, Eigen::Dynamic> gammah(anglesCount); // gamma half - resulting rotation angle (around n)
Matrix<uint32_t, 3, Eigen::Dynamic> validSinogramPoints_matrix(3, anglesCount); Matrix<uint32_t, 3, Eigen::Dynamic> validSinogramPoints_matrix(3, anglesCount); // coordinates (3D) of sampled curve in sinogram
Matrix<uint32_t, 1, Eigen::Dynamic> validSinogramPoints_matrix_lin(1, anglesCount); Matrix<uint32_t, 1, Eigen::Dynamic> validSinogramPoints_matrix_lin(1, anglesCount); //linear indices (1D) of validSinogramPoints_matrix
Matrix<uint32_t, 28, 1> validSinogramPoints_matrix_lin_dilated; Matrix<uint32_t, 28, 1> validSinogramPoints_matrix_lin_dilated; //all voxels neigboring one single voxel (in sinogram)
int measuredPeaksCount = ucsDirections.cols(); int measuredPeaksCount = ucsDirections.cols();
cout << "peak count used for indexing: " << measuredPeaksCount << endl; cout << "peak count used for indexing: " << measuredPeaksCount << endl;
...@@ -99,9 +99,9 @@ namespace pinkIndexer ...@@ -99,9 +99,9 @@ namespace pinkIndexer
if (candidateReflectionDirections.cols() == 0) if (candidateReflectionDirections.cols() == 0)
continue; continue;
sinogram_oneMeasuredPeak.setZero(); sinogram_oneMeasuredPeak.setZero(); //sinogram for one Bragg spot (from detector)
int candidatePeaksCount = candidateReflectionDirections.cols(); int candidatePeaksCount = candidateReflectionDirections.cols();
for (int candidatePeakNumber = 0; candidatePeakNumber < candidatePeaksCount; candidatePeakNumber++) for (int candidatePeakNumber = 0; candidatePeakNumber < candidatePeaksCount; candidatePeakNumber++) //calculation of curves for each candidate peak
{ {
const auto& candidateReflectionDirection = candidateReflectionDirections.col(candidatePeakNumber); const auto& candidateReflectionDirection = candidateReflectionDirections.col(candidatePeakNumber);
...@@ -116,6 +116,7 @@ namespace pinkIndexer ...@@ -116,6 +116,7 @@ namespace pinkIndexer
validSinogramPoints_matrix = ((validSinogramPoints * realToMatrixScaling).array() + realToMatrixOffset).matrix().cast<uint32_t>(); validSinogramPoints_matrix = ((validSinogramPoints * realToMatrixScaling).array() + realToMatrixOffset).matrix().cast<uint32_t>();
validSinogramPoints_matrix_lin = validSinogramPoints_matrix.row(0) + strides * validSinogramPoints_matrix.bottomRows(2); validSinogramPoints_matrix_lin = validSinogramPoints_matrix.row(0) + strides * validSinogramPoints_matrix.bottomRows(2);
// for each point from the curve, dilate it and insert into sinogram
for (uint32_t *linIdx = validSinogramPoints_matrix_lin.data(), for (uint32_t *linIdx = validSinogramPoints_matrix_lin.data(),
*end = validSinogramPoints_matrix_lin.data() + validSinogramPoints_matrix_lin.size(); *end = validSinogramPoints_matrix_lin.data() + validSinogramPoints_matrix_lin.size();
linIdx < end; linIdx++) linIdx < end; linIdx++)
...@@ -287,6 +288,7 @@ namespace pinkIndexer ...@@ -287,6 +288,7 @@ namespace pinkIndexer
maxElementSub(1) = maxElementIdx / strides(0); maxElementSub(1) = maxElementIdx / strides(0);
maxElementSub(0) = maxElementIdx - maxElementSub(1) * strides(0); maxElementSub(0) = maxElementIdx - maxElementSub(1) * strides(0);
// TODO: might be improved by proper search of centor of mass
Vector3f centerOfMassSub_f; Vector3f centerOfMassSub_f;
Matrix<uint32_t, 3, 1> centerOfMassSub, oldCenterOfMassSub = maxElementSub; Matrix<uint32_t, 3, 1> centerOfMassSub, oldCenterOfMassSub = maxElementSub;
int maxIterationCount = 5; int maxIterationCount = 5;
...@@ -299,6 +301,7 @@ namespace pinkIndexer ...@@ -299,6 +301,7 @@ namespace pinkIndexer
oldCenterOfMassSub = centerOfMassSub; oldCenterOfMassSub = centerOfMassSub;
} }
// calculation from siogram coordinates to vector coordinates
Vector3f bestRotationVector = centerOfMassSub_f - Vector3f(sinogramCenter, sinogramCenter, sinogramCenter); Vector3f bestRotationVector = centerOfMassSub_f - Vector3f(sinogramCenter, sinogramCenter, sinogramCenter);
Vector3f bestAxis = bestRotationVector.normalized(); Vector3f bestAxis = bestRotationVector.normalized();
float bestAngle = tan(bestRotationVector.norm() * sinogramScale * atan(M_PI / 4)) * 4; float bestAngle = tan(bestRotationVector.norm() * sinogramScale * atan(M_PI / 4)) * 4;
......
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