Commit 17f4230b authored by Yaroslav Gevorkov's avatar Yaroslav Gevorkov
Browse files

better handling of too high-q peaks

parent 83fd095f
...@@ -8,9 +8,9 @@ using namespace Eigen; ...@@ -8,9 +8,9 @@ using namespace Eigen;
ReflectionsInRangeFinder::ReflectionsInRangeFinder(const Lattice& lattice) ReflectionsInRangeFinder::ReflectionsInRangeFinder(const Lattice& lattice)
{ {
Matrix3f basis = lattice.getBasis(); const Matrix3f basis = lattice.getBasis();
int maxMillerIndex = 100; int maxMillerIndex = 175;
maxRadius = maxMillerIndex * lattice.getBasisVectorNorms().minCoeff(); maxRadius = maxMillerIndex * lattice.getBasisVectorNorms().minCoeff();
int millerIndicesPerDirection = 2 * maxMillerIndex + 1; int millerIndicesPerDirection = 2 * maxMillerIndex + 1;
Matrix3Xf reflections(3, (int)pow(millerIndicesPerDirection, 3)); Matrix3Xf reflections(3, (int)pow(millerIndicesPerDirection, 3));
...@@ -33,13 +33,13 @@ ReflectionsInRangeFinder::ReflectionsInRangeFinder(const Lattice& lattice) ...@@ -33,13 +33,13 @@ ReflectionsInRangeFinder::ReflectionsInRangeFinder(const Lattice& lattice)
iota(sortIndices.begin(), sortIndices.end(), 0); iota(sortIndices.begin(), sortIndices.end(), 0);
sort(sortIndices.begin(), sortIndices.end(), [&norms](uint32_t i, uint32_t j) { return norms[i] < norms[j]; }); sort(sortIndices.begin(), sortIndices.end(), [&norms](uint32_t i, uint32_t j) { return norms[i] < norms[j]; });
reflectionsDirections_sorted.resize(3, reflections.cols() - 1); //leave out reflection (0,0,0) reflectionsDirections_sorted.resize(3, reflections.cols() - 1); // leave out reflection (0,0,0)
norms_sorted.resize(reflections.cols()); norms_sorted.resize(reflections.cols());
for (int i = 1; i < reflections.cols(); ++i) for (int i = 1; i < reflections.cols(); ++i)
{ {
uint32_t sortIndex = sortIndices[i]; uint32_t sortIndex = sortIndices[i];
norms_sorted[i] = norms(sortIndex); norms_sorted[i] = norms(sortIndex);
reflectionsDirections_sorted.col(i-1) = reflections.col(sortIndex) / norms_sorted[i]; reflectionsDirections_sorted.col(i - 1) = reflections.col(sortIndex) / norms_sorted[i];
} }
} }
...@@ -48,17 +48,25 @@ void ReflectionsInRangeFinder::getReflectionsInRanges(EigenSTL::vector_Matrix3Xf ...@@ -48,17 +48,25 @@ void ReflectionsInRangeFinder::getReflectionsInRanges(EigenSTL::vector_Matrix3Xf
{ {
if (ranges.row(1).maxCoeff() > maxRadius) if (ranges.row(1).maxCoeff() > maxRadius)
{ {
throw BadInputException("getReflectionsInRanges is called with too large ranges! Limit the maximum resolution of the Bragg peaks!"); cerr << "the maximum resolution of Bragg spots exceeds the hardcoded limit of " << maxRadius
<< "A^-1. These Bragg spots will not be used for indexing, but they will be used for refinement!" << endl;
} }
candidateReflectionsDirections.resize(ranges.cols()); candidateReflectionsDirections.resize(ranges.cols());
for (int i = 0; i < ranges.cols(); i++) for (int i = 0; i < ranges.cols(); i++)
{ {
vector<float>::iterator low, up; if (ranges(1, i) > maxRadius)
low = lower_bound(norms_sorted.begin(), norms_sorted.end(), ranges(0, i)); {
up = upper_bound(norms_sorted.begin(), norms_sorted.end(), ranges(1, i)); candidateReflectionsDirections[i] = Matrix3Xf(3, 0);
}
else
{
vector<float>::iterator low, up;
low = lower_bound(norms_sorted.begin(), norms_sorted.end(), ranges(0, i));
up = upper_bound(norms_sorted.begin(), norms_sorted.end(), ranges(1, i));
candidateReflectionsDirections[i] = reflectionsDirections_sorted.block(0, low - norms_sorted.begin(), 3, up - low); candidateReflectionsDirections[i] = reflectionsDirections_sorted.block(0, low - norms_sorted.begin(), 3, up - low);
}
} }
} }
\ No newline at end of file
...@@ -86,7 +86,7 @@ void Sinogram::computeSinogram(const Eigen::Matrix3Xf& ucsDirections, const Eige ...@@ -86,7 +86,7 @@ void Sinogram::computeSinogram(const Eigen::Matrix3Xf& ucsDirections, const Eige
Matrix<uint32_t, 28, 1> validSinogramPoints_matrix_lin_dilated; Matrix<uint32_t, 28, 1> validSinogramPoints_matrix_lin_dilated;
int measuredPeaksCount = ucsDirections.cols(); int measuredPeaksCount = ucsDirections.cols();
cout << "peak count used for indexing: " << measuredPeaksCount << endl; cout << "peak count used for indexing: " << measuredPeaksCount << endl;
for (int measuredPeakNumber = 0; measuredPeakNumber < measuredPeaksCount; measuredPeakNumber++) for (int measuredPeakNumber = 0; measuredPeakNumber < measuredPeaksCount; measuredPeakNumber++)
{ {
if (measuredPeakNumber % 16 == 0) if (measuredPeakNumber % 16 == 0)
...@@ -97,6 +97,8 @@ void Sinogram::computeSinogram(const Eigen::Matrix3Xf& ucsDirections, const Eige ...@@ -97,6 +97,8 @@ void Sinogram::computeSinogram(const Eigen::Matrix3Xf& ucsDirections, const Eige
Vector3f l = ucsDirections.col(measuredPeakNumber); // rotation axis Vector3f l = ucsDirections.col(measuredPeakNumber); // rotation axis
Matrix3Xf& candidateReflectionDirections = candidateReflectionsDirections[measuredPeakNumber]; Matrix3Xf& candidateReflectionDirections = candidateReflectionsDirections[measuredPeakNumber];
if (candidateReflectionDirections.cols() == 0)
continue;
fill(sinogram_oneMeasuredPeak.begin(), sinogram_oneMeasuredPeak.end(), 0); fill(sinogram_oneMeasuredPeak.begin(), sinogram_oneMeasuredPeak.end(), 0);
int candidatePeaksCount = candidateReflectionDirections.cols(); int candidatePeaksCount = candidateReflectionDirections.cols();
...@@ -193,17 +195,19 @@ void Sinogram::computeSinogramParallel(const Eigen::Matrix3Xf& ucsDirections, co ...@@ -193,17 +195,19 @@ void Sinogram::computeSinogramParallel(const Eigen::Matrix3Xf& ucsDirections, co
sinogram_oneMeasuredPeak_matrixMapped(sinogram_oneMeasuredPeak.data(), sinogram_oneMeasuredPeak.size()); sinogram_oneMeasuredPeak_matrixMapped(sinogram_oneMeasuredPeak.data(), sinogram_oneMeasuredPeak.size());
int measuredPeaksCount = ucsDirections.cols(); int measuredPeaksCount = ucsDirections.cols();
cout << "peak count used for indexing: " << measuredPeaksCount << endl; cout << "peak count used for indexing: " << measuredPeaksCount << endl;
for (int measuredPeakNumber = 0; measuredPeakNumber < measuredPeaksCount; measuredPeakNumber++) for (int measuredPeakNumber = 0; measuredPeakNumber < measuredPeaksCount; measuredPeakNumber++)
{ {
if (measuredPeakNumber % 16 == 0) if (measuredPeakNumber % 16 == 0)
{ {
cout << "computed " << 100 * (float)measuredPeakNumber / measuredPeaksCount << "%" << endl; cout << "computed " << 100 * (float)measuredPeakNumber / measuredPeaksCount << "%" << endl;
} }
Vector3f l = ucsDirections.col(measuredPeakNumber); // rotation axis Vector3f l = ucsDirections.col(measuredPeakNumber); // rotation axis
Matrix3Xf& candidateReflectionDirections = candidateReflectionsDirections[measuredPeakNumber]; Matrix3Xf& candidateReflectionDirections = candidateReflectionsDirections[measuredPeakNumber];
if (candidateReflectionDirections.cols() == 0)
continue;
fill(sinogram_oneMeasuredPeak.begin(), sinogram_oneMeasuredPeak.end(), 0); fill(sinogram_oneMeasuredPeak.begin(), sinogram_oneMeasuredPeak.end(), 0);
...@@ -241,15 +245,17 @@ void Sinogram::computeSinogramParallel2(const Eigen::Matrix3Xf& ucsDirections, c ...@@ -241,15 +245,17 @@ void Sinogram::computeSinogramParallel2(const Eigen::Matrix3Xf& ucsDirections, c
cout << "peak count used for indexing: " << measuredPeaksCount << endl; cout << "peak count used for indexing: " << measuredPeaksCount << endl;
for (int measuredPeakNumber = 0; measuredPeakNumber < measuredPeaksCount; measuredPeakNumber++) for (int measuredPeakNumber = 0; measuredPeakNumber < measuredPeaksCount; measuredPeakNumber++)
{ {
if (measuredPeakNumber % 16 == 0) if (measuredPeakNumber % 16 == 0)
{ {
cout << "computed " << 100 * (float)measuredPeakNumber / measuredPeaksCount << "%" << endl; cout << "computed " << 100 * (float)measuredPeakNumber / measuredPeaksCount << "%" << endl;
} }
Vector3f l = ucsDirections.col(measuredPeakNumber); // rotation axis Vector3f l = ucsDirections.col(measuredPeakNumber); // rotation axis
Matrix3Xf& candidateReflectionDirections = candidateReflectionsDirections[measuredPeakNumber]; Matrix3Xf& candidateReflectionDirections = candidateReflectionsDirections[measuredPeakNumber];
if (candidateReflectionDirections.cols() == 0)
continue;
fill(sinogram_oneMeasuredPeak.begin(), sinogram_oneMeasuredPeak.end(), 0); fill(sinogram_oneMeasuredPeak.begin(), sinogram_oneMeasuredPeak.end(), 0);
......
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