Commit 5ddb41c5 authored by Yaroslav Gevorkov's avatar Yaroslav Gevorkov
Browse files

added maxResolutionForIndexing_1_per_A

parent bb4c144a
......@@ -42,14 +42,14 @@ class PinkIndexer
};
PinkIndexer(const ExperimentSettings& experimentSettings, ConsideredPeaksCount consideredPeaksCount, AngleResolution angleResolution,
RefinementType refinementType);
RefinementType refinementType, float maxResolutionForIndexing_1_per_A);
void indexPattern(Lattice& indexedLattice, const Eigen::Matrix2Xf& detectorPeaks_m, const Eigen::RowVectorXf& intensities, int threadCount);
void indexPattern(Lattice& indexedLattice, const Eigen::Matrix3Xf& meanReciprocalPeaks_1_per_A, const Eigen::RowVectorXf& intensities, int threadCount);
void indexPattern(Lattice& indexedLattice, Eigen::RowVectorXf& intensities, const Eigen::Matrix2Xf& detectorPeaks_m, int threadCount);
void indexPattern(Lattice& indexedLattice, Eigen::RowVectorXf& intensities, const Eigen::Matrix3Xf& meanReciprocalPeaks_1_per_A, int threadCount);
private:
void raducePeakCount(Eigen::Matrix3Xf& ucsDirections, Eigen::Array2Xf& ucsBorderNorms, const Eigen::Matrix2Xf& detectorPeaks_m,
const Eigen::RowVectorXf& intensities);
void raducePeakCount(Eigen::Matrix3Xf& ucsDirections, Eigen::Array2Xf& ucsBorderNorms, Eigen::RowVectorXf& intensities,
const Eigen::Matrix2Xf& detectorPeaks_m);
float getAngleResolution();
int getConsideredPeaksCount();
......@@ -63,4 +63,5 @@ class PinkIndexer
ConsideredPeaksCount consideredPeaksCount;
AngleResolution angleResolution;
RefinementType refinementType;
float maxResolutionForIndexing_1_per_A;
};
\ No newline at end of file
......@@ -41,11 +41,11 @@ typedef enum {
typedef struct PinkIndexer PinkIndexer;
PinkIndexer* PinkIndexer_new(ExperimentSettings* experimentSettings, consideredPeaksCount_t consideredPeaksCount, angleResolution_t angleResolution,
refinementType_t refinementType);
refinementType_t refinementType, float maxResolutionForIndexing_1_per_A);
void PinkIndexer_delete(PinkIndexer* pinkIndexer);
void PinkIndexer_indexPattern(PinkIndexer* pinkIndexer, Lattice_t* indexedLattice, const reciprocalPeaks_1_per_A_t meanReciprocalPeaks_1_per_A, const float* intensities,
int peakCount, int threadCount);
void PinkIndexer_indexPattern(PinkIndexer* pinkIndexer, Lattice_t* indexedLattice, const reciprocalPeaks_1_per_A_t meanReciprocalPeaks_1_per_A,
const float* intensities, int peakCount, int threadCount);
#ifdef __cplusplus
......
......@@ -9,13 +9,14 @@ using namespace Eigen;
PinkIndexer::PinkIndexer(const ExperimentSettings& experimentSettings, ConsideredPeaksCount consideredPeaksCount, AngleResolution angleResolution,
RefinementType refinementType)
RefinementType refinementType, float maxResolutionForIndexing_1_per_A)
: reciprocalToRealProjection(experimentSettings)
, backprojection(experimentSettings)
, sinogram(experimentSettings.getSampleReciprocalLattice_1A())
, consideredPeaksCount(consideredPeaksCount)
, angleResolution(angleResolution)
, refinementType(refinementType)
, maxResolutionForIndexing_1_per_A(maxResolutionForIndexing_1_per_A)
{
sampleLattice = experimentSettings.getSampleReciprocalLattice_1A();
......@@ -23,14 +24,14 @@ PinkIndexer::PinkIndexer(const ExperimentSettings& experimentSettings, Considere
sinogram.setSinogramAngleResolution(angleResolution_deg);
}
void PinkIndexer::indexPattern(Lattice& indexedLattice, const Matrix3Xf& meanReciprocalPeaks_1_per_A, const Eigen::RowVectorXf& intensities, int threadCount)
void PinkIndexer::indexPattern(Lattice& indexedLattice, Eigen::RowVectorXf& intensities, const Matrix3Xf& meanReciprocalPeaks_1_per_A, int threadCount)
{
Matrix2Xf detectorPeaks_m;
reciprocalToRealProjection.project(detectorPeaks_m, meanReciprocalPeaks_1_per_A);
indexPattern(indexedLattice, detectorPeaks_m, intensities, threadCount);
indexPattern(indexedLattice, intensities, detectorPeaks_m, threadCount);
}
void PinkIndexer::indexPattern(Lattice& indexedLattice, const Matrix2Xf& detectorPeaks_m, const Eigen::RowVectorXf& intensities, int threadCount)
void PinkIndexer::indexPattern(Lattice& indexedLattice, Eigen::RowVectorXf& intensities, const Matrix2Xf& detectorPeaks_m, int threadCount)
{
Matrix3Xf ucsDirections;
Array2Xf ucsBorderNorms;
......@@ -38,7 +39,7 @@ void PinkIndexer::indexPattern(Lattice& indexedLattice, const Matrix2Xf& detecto
Matrix3Xf ucsDirections_reduced = ucsDirections;
Array2Xf ucsBorderNorms_reduced = ucsBorderNorms;
raducePeakCount(ucsDirections_reduced, ucsBorderNorms_reduced, detectorPeaks_m, intensities);
raducePeakCount(ucsDirections_reduced, ucsBorderNorms_reduced, intensities, detectorPeaks_m);
if (threadCount == 1)
{
......@@ -48,7 +49,7 @@ void PinkIndexer::indexPattern(Lattice& indexedLattice, const Matrix2Xf& detecto
{
sinogram.computeSinogramParallel(ucsDirections_reduced, ucsBorderNorms_reduced, threadCount);
}
else //probably not run in parellel to other instances
else // probably not run in parellel to other instances
{
int slaveThreadCount = threadCount - 1;
sinogram.computeSinogramParallel2(ucsDirections_reduced, ucsBorderNorms_reduced, slaveThreadCount);
......@@ -82,16 +83,30 @@ void PinkIndexer::indexPattern(Lattice& indexedLattice, const Matrix2Xf& detecto
// sinogram.saveToFile("C:\\DesyFiles\\workspaces\\VisualStudio_workspace\\pinkIndexer\\workfolder\\sinogram");
}
void PinkIndexer::raducePeakCount(Matrix3Xf& ucsDirections, Array2Xf& ucsBorderNorms, const Eigen::Matrix2Xf& detectorPeaks_m, const RowVectorXf& intensities)
void PinkIndexer::raducePeakCount(Matrix3Xf& ucsDirections, Array2Xf& ucsBorderNorms, RowVectorXf& intensities, const Eigen::Matrix2Xf& detectorPeaks_m)
{
//first clear all above resolution limit
int peaksKept_res = 0;
for (int i = 0; i < intensities.size(); i++)
{
if (ucsBorderNorms(1, i) >= maxResolutionForIndexing_1_per_A)
{
ucsDirections.col(peaksKept_res) = ucsDirections.col(i);
ucsBorderNorms.col(peaksKept_res) = ucsBorderNorms.col(i);
intensities[peaksKept_res] = intensities[i];
peaksKept_res++;
}
}
ucsDirections.conservativeResize(3, peaksKept_res);
ucsBorderNorms.conservativeResize(2, peaksKept_res);
intensities.conservativeResize(peaksKept_res);
//then clear some at high and some at low resolution
int consideredPeaksCount = getConsideredPeaksCount();
if (consideredPeaksCount >= intensities.size())
return;
RowVectorXf intensities_reduced = intensities;
Matrix2Xf detectorPeaks_debug = detectorPeaks_m;
int peaksToRemoveCount = intensities.size() - consideredPeaksCount;
int peaksToRemoveByNormMin = min(0.1 * peaksToRemoveCount, 20.0);
int peaksToRemoveByNormMax = 0.75 * peaksToRemoveCount;
......@@ -110,13 +125,12 @@ void PinkIndexer::raducePeakCount(Matrix3Xf& ucsDirections, Array2Xf& ucsBorderN
{
ucsDirections.col(peaksKept_norms) = ucsDirections.col(i);
ucsBorderNorms.col(peaksKept_norms) = ucsBorderNorms.col(i);
detectorPeaks_debug.col(peaksKept_norms) = detectorPeaks_debug.col(i);
intensities_reduced[peaksKept_norms] = intensities_reduced[i];
intensities[peaksKept_norms] = intensities[i];
peaksKept_norms++;
}
}
RowVectorXf intensities_sorted = intensities_reduced;
RowVectorXf intensities_sorted = intensities;
nth_element((float*)intensities_sorted.data(), (float*)intensities_sorted.data() + peaksToRemoveByIntensity,
(float*)intensities_sorted.data() + peaksKept_norms);
float minIntensity = intensities_sorted[peaksToRemoveByIntensity];
......@@ -124,23 +138,18 @@ void PinkIndexer::raducePeakCount(Matrix3Xf& ucsDirections, Array2Xf& ucsBorderN
int peaksKept_normsIntensities = 0;
for (int i = 0; i < peaksKept_norms; i++)
{
if (intensities_reduced[i] >= minIntensity)
if (intensities[i] >= minIntensity)
{
ucsDirections.col(peaksKept_normsIntensities) = ucsDirections.col(i);
ucsBorderNorms.col(peaksKept_normsIntensities) = ucsBorderNorms.col(i);
detectorPeaks_debug.col(peaksKept_normsIntensities) = detectorPeaks_debug.col(i);
// intensities_reduced[peaksKept_normsIntensities] = intensities_reduced[i];
intensities[peaksKept_normsIntensities] = intensities[i];
peaksKept_normsIntensities++;
}
}
ucsDirections.conservativeResize(3, consideredPeaksCount);
ucsBorderNorms.conservativeResize(2, consideredPeaksCount);
detectorPeaks_debug.conservativeResize(2, consideredPeaksCount);
ofstream myfile("C:\\DesyFiles\\workspaces\\VisualStudio_workspace\\pinkIndexer\\workfolder\\detectorPeaks_m_reduced");
myfile << detectorPeaks_debug;
myfile.close();
intensities.conservativeResize(consideredPeaksCount);
}
float PinkIndexer::getAngleResolution()
......
......@@ -3,7 +3,7 @@
extern "C" PinkIndexer* PinkIndexer_new(ExperimentSettings* experimentSettings, consideredPeaksCount_t consideredPeaksCount, angleResolution_t angleResolution,
refinementType_t refinementType)
refinementType_t refinementType, float maxResolutionForIndexing_1_per_A)
{
PinkIndexer::ConsideredPeaksCount consideredPeaksCount_enumClass;
switch (consideredPeaksCount)
......@@ -74,7 +74,8 @@ extern "C" PinkIndexer* PinkIndexer_new(ExperimentSettings* experimentSettings,
break;
}
return new PinkIndexer(*experimentSettings, consideredPeaksCount_enumClass, angleResolution_enumClass, refinementType_enumClass);
return new PinkIndexer(*experimentSettings, consideredPeaksCount_enumClass, angleResolution_enumClass, refinementType_enumClass,
maxResolutionForIndexing_1_per_A);
}
extern "C" void PinkIndexer_delete(PinkIndexer* pinkIndexer)
......@@ -83,7 +84,7 @@ extern "C" void PinkIndexer_delete(PinkIndexer* pinkIndexer)
}
extern "C" void PinkIndexer_indexPattern(PinkIndexer* pinkIndexer, Lattice_t* indexedLattice, const reciprocalPeaks_1_per_A_t meanReciprocalPeaks_1_per_A,
const float* intensities, int peakCount, int threadCount)
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++)
......@@ -94,7 +95,7 @@ extern "C" void PinkIndexer_indexPattern(PinkIndexer* pinkIndexer, Lattice_t* in
Lattice indexedLattice_class;
Eigen::RowVectorXf intensities_class = Eigen::Map<Eigen::RowVectorXf>((float*)intensities, 1, peakCount);
pinkIndexer->indexPattern(indexedLattice_class, meanReciprocalPeaks_1_per_A_matrix, intensities_class, threadCount);
pinkIndexer->indexPattern(indexedLattice_class, intensities_class, meanReciprocalPeaks_1_per_A_matrix, threadCount);
Eigen::Matrix3f basis = indexedLattice_class.getBasis();
indexedLattice->ax = basis(0, 0);
......
......@@ -49,13 +49,14 @@ 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 = numeric_limits<float>::max();
Lattice indexedLattice;
PinkIndexer pinkIndexer(experimentSettings, consideredPeaksCount, angleResolution, refinementType);
PinkIndexer pinkIndexer(experimentSettings, consideredPeaksCount, angleResolution, refinementType, maxResolutionForIndexing_1_per_A);
{
auto tmp = Chronometer("sinogram");
int threadCount = 4;
pinkIndexer.indexPattern(indexedLattice, peaksOnDetector_m, intensities, threadCount);
pinkIndexer.indexPattern(indexedLattice, intensities, peaksOnDetector_m, threadCount);
}
SimpleDiffractionPatternPrediction simpleDiffractionPatternPrediction(experimentSettings);
......
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