Commit 8da3bb67 authored by Yaroslav Gevorkov's avatar Yaroslav Gevorkov
Browse files

pinkIndexer now returns indexed status

parent e9f48ed2
......@@ -12,7 +12,7 @@ if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES)
"MinSizeRel" "RelWithDebInfo")
endif()
find_package(Eigen3 3.3 NO_MODULE)
find_package(Eigen3 3.4 NO_MODULE)
include_directories(include)
......
Take into account:
- different detector counts for different lambda (numbers of electrons released by photon differs)
- differen crossection dependent on lambda
- possibility of SAD/MAD phasing with pink
\ No newline at end of file
......@@ -19,13 +19,6 @@ class ExperimentSettings
ExperimentSettings(float beamEenergy_eV, float detectorDistance_m, float detectorRadius_m, float divergenceAngle_deg, float nonMonochromaticity,
float minRealLatticeVectorLength_A, float maxRealLatticeVectorLength_A, float reflectionRadius_1_per_A);
// ExperimentSettings(float coffset_m, float clen_mm, float beamEenergy_eV, float divergenceAngle_deg, float nonMonochromaticity, float pixelLength_m,
// float detectorRadius_pixel, float minRealLatticeVectorLength_A, float maxRealLatticeVectorLength_A, float minRealLatticeDeterminant_A3,
// float maxRealLatticeDeterminant_A3, float reflectionRadius_1_per_A);
// ExperimentSettings(float beamEenergy_eV, float detectorDistance_m, float detectorRadius_m, float divergenceAngle_deg, float nonMonochromaticity,
// float minRealLatticeVectorLength_A, float maxRealLatticeVectorLength_A, float minRealLatticeDeterminant_A3,
// float maxRealLatticeDeterminant_A3, float reflectionRadius_1_per_A);
ExperimentSettings(float coffset_m, float clen_mm, float beamEenergy_eV, float divergenceAngle_deg, float nonMonochromaticity, float pixelLength_m,
float detectorRadius_pixel, const Lattice& sampleReciprocalLattice_1A, float tolerance, float reflectionRadius_1_per_A);
ExperimentSettings(float beamEenergy_eV, float detectorDistance_m, float detectorRadius_m, float divergenceAngle_deg, float nonMonochromaticity,
......
......@@ -44,8 +44,8 @@ class PinkIndexer
PinkIndexer(const ExperimentSettings& experimentSettings, ConsideredPeaksCount consideredPeaksCount, AngleResolution angleResolution,
RefinementType refinementType, float maxResolutionForIndexing_1_per_A);
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);
bool indexPattern(Lattice& indexedLattice, Eigen::RowVectorXf& intensities, const Eigen::Matrix2Xf& detectorPeaks_m, int threadCount);
bool indexPattern(Lattice& indexedLattice, Eigen::RowVectorXf& intensities, const Eigen::Matrix3Xf& meanReciprocalPeaks_1_per_A, int threadCount);
private:
void reducePeakCount(Eigen::Matrix3Xf& ucsDirections, Eigen::Array2Xf& ucsBorderNorms, Eigen::RowVectorXf& intensities,
......
......@@ -3,15 +3,35 @@
#include "Lattice.h"
#include <Eigen/Dense>
#include <vector>
class Refinement
{
public:
Refinement();
Refinement(float tolerance);
void refineFixedLattice(Lattice& lattice, const Eigen::Matrix3Xf& ucsDirections, const Eigen::Matrix2Xf& ucsBorderNorms);
void refineVariableLattice(Lattice& lattice, const Eigen::Matrix3Xf& ucsDirections, const Eigen::Matrix2Xf& ucsBorderNorms);
int getFittedPeaksCount(Lattice& lattice, const Eigen::Matrix3Xf& ucsDirections, const Eigen::Matrix2Xf& ucsBorderNorms);
private:
double getDefect(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);
double getMeanDefect(const Eigen::Matrix3f& basis, const Eigen::Matrix3Xf& ucsDirections, const Eigen::Matrix2Xf& ucsBorderNorms);
float tolerance;
// preallocation
Eigen::ArrayXf defects;
Eigen::Matrix2Xf ucsBorderNormsSquared;
Eigen::Matrix3Xf millerIndices_close;
Eigen::Matrix3Xf millerIndices_far;
Eigen::Matrix3Xf candidatePeaks;
Eigen::RowVectorXf candidatePeaksNormsSquared;
std::vector<Eigen::Vector3f> millerIndices;
Eigen::RowVectorXf projectedVectorNorms;
Eigen::Matrix3Xf defectVectors_absolute;
Eigen::Matrix3Xf defectVectors_relative;
Eigen::Array<bool, 1, Eigen::Dynamic> notPredictablePeaks;
Eigen::ArrayXf meanDefects;
};
......@@ -44,8 +44,8 @@ PinkIndexer* PinkIndexer_new(ExperimentSettings* experimentSettings, consideredP
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);
int 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
......
/*
* eigenSTLContainers.h
*
* Created on: 21.12.2016
* Author: Yaro
*/
#ifndef EIGENSTLCONTAINERS_H_
#define EIGENSTLCONTAINERS_H_
......
......@@ -62,48 +62,6 @@ ExperimentSettings::ExperimentSettings(float beamEenergy_eV, float detectorDista
differentRealLatticeVectorLengths_A[1] = maxRealLatticeVectorLength_A;
}
// ExperimentSettings::ExperimentSettings(float coffset_m, float clen_mm, float beamEenergy_eV, float divergenceAngle_deg, float nonMonochromaticity,
// float pixelLength_m, float detectorRadius_pixel, float minRealLatticeVectorLength_A, float
// maxRealLatticeVectorLength_A, float minRealLatticeDeterminant_A3, float maxRealLatticeDeterminant_A3)
// : latticeParametersKnown(false)
// , minRealLatticeVectorLength_A(minRealLatticeVectorLength_A)
// , maxRealLatticeVectorLength_A(maxRealLatticeVectorLength_A)
// , minRealLatticeDeterminant_A3(minRealLatticeDeterminant_A3)
// , maxRealLatticeDeterminant_A3(maxRealLatticeDeterminant_A3)
//{
// constructFromGeometryFileValues(coffset_m, clen_mm, beamEenergy_eV, divergenceAngle_deg, nonMonochromaticity, pixelLength_m, detectorRadius_pixel);
//
// minReciprocalLatticeVectorLength_1A = 1 / maxRealLatticeVectorLength_A;
// maxReciprocalLatticeVectorLength_1A = 1 / minRealLatticeVectorLength_A;
// minReciprocalLatticeDeterminant_1A3 = 1 / maxRealLatticeDeterminant_A3;
// maxReciprocalLatticeDeterminant_1A3 = 1 / minRealLatticeDeterminant_A3;
//
// differentRealLatticeVectorLengths_A.resize(2);
// differentRealLatticeVectorLengths_A[0] = minRealLatticeVectorLength_A;
// differentRealLatticeVectorLengths_A[1] = maxRealLatticeVectorLength_A;
//}
//
// ExperimentSettings::ExperimentSettings(float beamEenergy_eV, float detectorDistance_m, float detectorRadius_m, float divergenceAngle_deg,
// float nonMonochromaticity, float minRealLatticeVectorLength_A, float maxRealLatticeVectorLength_A,
// float minRealLatticeDeterminant_A3, float maxRealLatticeDeterminant_A3)
// : latticeParametersKnown(false)
// , minRealLatticeVectorLength_A(minRealLatticeVectorLength_A)
// , maxRealLatticeVectorLength_A(maxRealLatticeVectorLength_A)
// , minRealLatticeDeterminant_A3(minRealLatticeDeterminant_A3)
// , maxRealLatticeDeterminant_A3(maxRealLatticeDeterminant_A3)
//{
// constructFromPrecomputedValues(beamEenergy_eV, detectorDistance_m, detectorRadius_m, divergenceAngle_deg, nonMonochromaticity);
//
// minReciprocalLatticeVectorLength_1A = 1 / maxRealLatticeVectorLength_A;
// maxReciprocalLatticeVectorLength_1A = 1 / minRealLatticeVectorLength_A;
// minReciprocalLatticeDeterminant_1A3 = 1 / maxRealLatticeDeterminant_A3;
// maxReciprocalLatticeDeterminant_1A3 = 1 / minRealLatticeDeterminant_A3;
//
// differentRealLatticeVectorLengths_A.resize(2);
// differentRealLatticeVectorLengths_A[0] = minRealLatticeVectorLength_A;
// differentRealLatticeVectorLengths_A[1] = maxRealLatticeVectorLength_A;
//}
ExperimentSettings::ExperimentSettings(float coffset_m, float clen_mm, float beamEenergy_eV, float divergenceAngle_deg, float nonMonochromaticity,
float pixelLength_m, float detectorRadius_pixel, const Lattice& sampleReciprocalLattice_1A, float tolerance,
float reflectionRadius_1_per_A)
......
......@@ -13,6 +13,7 @@ PinkIndexer::PinkIndexer(const ExperimentSettings& experimentSettings, Considere
: reciprocalToRealProjection(experimentSettings)
, backprojection(experimentSettings)
, sinogram(experimentSettings.getSampleReciprocalLattice_1A())
, refinement(experimentSettings.getTolerance())
, consideredPeaksCount(consideredPeaksCount)
, angleResolution(angleResolution)
, refinementType(refinementType)
......@@ -24,14 +25,14 @@ PinkIndexer::PinkIndexer(const ExperimentSettings& experimentSettings, Considere
sinogram.setSinogramAngleResolution(angleResolution_deg);
}
void PinkIndexer::indexPattern(Lattice& indexedLattice, Eigen::RowVectorXf& intensities, const Matrix3Xf& meanReciprocalPeaks_1_per_A, int threadCount)
bool 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, intensities, detectorPeaks_m, threadCount);
return indexPattern(indexedLattice, intensities, detectorPeaks_m, threadCount);
}
void PinkIndexer::indexPattern(Lattice& indexedLattice, Eigen::RowVectorXf& intensities, const Matrix2Xf& detectorPeaks_m, int threadCount)
bool PinkIndexer::indexPattern(Lattice& indexedLattice, Eigen::RowVectorXf& intensities, const Matrix2Xf& detectorPeaks_m, int threadCount)
{
Matrix3Xf ucsDirections;
Array2Xf ucsBorderNorms;
......@@ -80,14 +81,24 @@ void PinkIndexer::indexPattern(Lattice& indexedLattice, Eigen::RowVectorXf& inte
break;
}
int fittedPeaksCount = refinement.getFittedPeaksCount(indexedLattice, ucsDirections, ucsBorderNorms);
if (fittedPeaksCount >= 20 || fittedPeaksCount >= detectorPeaks_m.cols() * 0.5)
{
return true;
}
else
{
return false;
}
// sinogram.saveToFile("C:\\DesyFiles\\workspaces\\VisualStudio_workspace\\pinkIndexer\\workfolder\\sinogram");
}
void PinkIndexer::reducePeakCount(Matrix3Xf& ucsDirections, Array2Xf& ucsBorderNorms, RowVectorXf& intensities, const Eigen::Matrix2Xf& detectorPeaks_m)
{
Matrix2Xf detectorPeaks_m_copy = detectorPeaks_m;
Matrix2Xf detectorPeaks_m_copy = detectorPeaks_m;
//first clear all above resolution limit
// first clear all above resolution limit
int peaksKept_res = 0;
for (int i = 0; i < intensities.size(); i++)
{
......@@ -96,16 +107,16 @@ void PinkIndexer::reducePeakCount(Matrix3Xf& ucsDirections, Array2Xf& ucsBorderN
ucsDirections.col(peaksKept_res) = ucsDirections.col(i);
ucsBorderNorms.col(peaksKept_res) = ucsBorderNorms.col(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++;
}
}
ucsDirections.conservativeResize(NoChange, peaksKept_res);
ucsBorderNorms.conservativeResize(NoChange, peaksKept_res);
intensities.conservativeResize(peaksKept_res);
detectorPeaks_m_copy.conservativeResize(NoChange, peaksKept_res);
detectorPeaks_m_copy.conservativeResize(NoChange, peaksKept_res);
//then clear some at high and some at low resolution
// then clear some at high and some at low resolution
int consideredPeaksCount = getConsideredPeaksCount();
if (consideredPeaksCount >= intensities.size())
......
......@@ -13,7 +13,11 @@ using namespace Eigen;
static void roundTowardsZero(Matrix3Xf& x);
static void roundAwayFromZero(Matrix3Xf& x);
Refinement::Refinement() {}
Refinement::Refinement(float tolerance)
: tolerance(tolerance)
{
millerIndices.reserve(500);
}
void Refinement::refineVariableLattice(Lattice& lattice, const Matrix3Xf& ucsDirections, const Matrix2Xf& ucsBorderNorms)
{
......@@ -23,40 +27,40 @@ void Refinement::refineVariableLattice(Lattice& lattice, const Matrix3Xf& ucsDir
float stepSize = lattice.getBasisVectorNorms().maxCoeff() * 0.002;
float minStepSize = lattice.getBasisVectorNorms().minCoeff() * 0.00001;
int maxStepsCount = 200;
ArrayXf defects(maxStepsCount);
defects[0] = getDefect(basis, ucsDirections, ucsBorderNorms);
meanDefects.resize(maxStepsCount);
meanDefects[0] = getMeanDefect(basis, ucsDirections, ucsBorderNorms);
for (int i = 0; i < maxStepsCount; i++)
{
// cout << defects[i] << endl;
cout << meanDefects[i] << endl;
Array33f gradient;
Matrix3f offsetBasis = basis;
offsetBasis(0, 0) += delta;
gradient(0, 0) = getDefect(offsetBasis, ucsDirections, ucsBorderNorms) - defects[i];
gradient(0, 0) = getMeanDefect(offsetBasis, ucsDirections, ucsBorderNorms) - meanDefects[i];
offsetBasis(0, 0) = basis(0, 0);
offsetBasis(1, 0) += delta;
gradient(1, 0) = getDefect(offsetBasis, ucsDirections, ucsBorderNorms) - defects[i];
gradient(1, 0) = getMeanDefect(offsetBasis, ucsDirections, ucsBorderNorms) - meanDefects[i];
offsetBasis(1, 0) = basis(1, 0);
offsetBasis(2, 0) += delta;
gradient(2, 0) = getDefect(offsetBasis, ucsDirections, ucsBorderNorms) - defects[i];
gradient(2, 0) = getMeanDefect(offsetBasis, ucsDirections, ucsBorderNorms) - meanDefects[i];
offsetBasis(2, 0) = basis(2, 0);
offsetBasis(0, 1) += delta;
gradient(0, 1) = getDefect(offsetBasis, ucsDirections, ucsBorderNorms) - defects[i];
gradient(0, 1) = getMeanDefect(offsetBasis, ucsDirections, ucsBorderNorms) - meanDefects[i];
offsetBasis(0, 1) = basis(0, 1);
offsetBasis(1, 1) += delta;
gradient(1, 1) = getDefect(offsetBasis, ucsDirections, ucsBorderNorms) - defects[i];
gradient(1, 1) = getMeanDefect(offsetBasis, ucsDirections, ucsBorderNorms) - meanDefects[i];
offsetBasis(1, 1) = basis(1, 1);
offsetBasis(2, 1) += delta;
gradient(2, 1) = getDefect(offsetBasis, ucsDirections, ucsBorderNorms) - defects[i];
gradient(2, 1) = getMeanDefect(offsetBasis, ucsDirections, ucsBorderNorms) - meanDefects[i];
offsetBasis(2, 1) = basis(2, 1);
offsetBasis(0, 2) += delta;
gradient(0, 2) = getDefect(offsetBasis, ucsDirections, ucsBorderNorms) - defects[i];
gradient(0, 2) = getMeanDefect(offsetBasis, ucsDirections, ucsBorderNorms) - meanDefects[i];
offsetBasis(0, 2) = basis(0, 2);
offsetBasis(1, 2) += delta;
gradient(1, 2) = getDefect(offsetBasis, ucsDirections, ucsBorderNorms) - defects[i];
gradient(1, 2) = getMeanDefect(offsetBasis, ucsDirections, ucsBorderNorms) - meanDefects[i];
offsetBasis(1, 2) = basis(1, 2);
offsetBasis(2, 2) += delta;
gradient(2, 2) = getDefect(offsetBasis, ucsDirections, ucsBorderNorms) - defects[i];
gradient(2, 2) = getMeanDefect(offsetBasis, ucsDirections, ucsBorderNorms) - meanDefects[i];
float norm = gradient.matrix().norm();
gradient = gradient / norm * stepSize;
......@@ -66,13 +70,13 @@ void Refinement::refineVariableLattice(Lattice& lattice, const Matrix3Xf& ucsDir
}
basis = basis - gradient.matrix();
defects[i + 1] = getDefect(basis, ucsDirections, ucsBorderNorms);
meanDefects[i + 1] = getMeanDefect(basis, ucsDirections, ucsBorderNorms);
if (defects[i + 1] > defects[i])
if (meanDefects[i + 1] > meanDefects[i])
{
stepSize = stepSize * 0.9;
if (i > 10 && (defects.segment(i - 4, 4).maxCoeff() - defects.segment(i - 4, 4).minCoeff()) / defects[i] < 0.01) // settled down
if (i > 10 && (meanDefects.segment(i - 4, 4).maxCoeff() - meanDefects.segment(i - 4, 4).minCoeff()) / meanDefects[i] < 0.01) // settled down
stepSize = stepSize * 0.2;
if (stepSize < minStepSize)
......@@ -97,26 +101,26 @@ void Refinement::refineFixedLattice(Lattice& lattice, const Matrix3Xf& ucsDirect
float stepSize = 0.1 / 180 * M_PI;
float minStepSize = 0.001 / 180 * M_PI;
int maxStepsCount = 200;
ArrayXf defects(maxStepsCount);
defects[0] = getDefect(basis, ucsDirections, ucsBorderNorms);
meanDefects.resize(maxStepsCount);
meanDefects[0] = getMeanDefect(basis, ucsDirections, ucsBorderNorms);
for (int i = 0; i < maxStepsCount; i++)
{
// cout << defects[i] << endl;
cout << meanDefects[i] << endl;
Vector3f gradient;
gradient(0) = getDefect(rotX * basis, ucsDirections, ucsBorderNorms) - defects[i];
gradient(1) = getDefect(rotY * basis, ucsDirections, ucsBorderNorms) - defects[i];
gradient(2) = getDefect(rotZ * basis, ucsDirections, ucsBorderNorms) - defects[i];
gradient(0) = getMeanDefect(rotX * basis, ucsDirections, ucsBorderNorms) - meanDefects[i];
gradient(1) = getMeanDefect(rotY * basis, ucsDirections, ucsBorderNorms) - meanDefects[i];
gradient(2) = getMeanDefect(rotZ * basis, ucsDirections, ucsBorderNorms) - meanDefects[i];
gradient = -gradient.normalized() * stepSize;
basis = AngleAxisf(gradient(0), Vector3f::UnitX()) * AngleAxisf(gradient(1), Vector3f::UnitY()) * AngleAxisf(gradient(2), Vector3f::UnitZ()) * basis;
defects[i + 1] = getDefect(basis, ucsDirections, ucsBorderNorms);
meanDefects[i + 1] = getMeanDefect(basis, ucsDirections, ucsBorderNorms);
if (defects[i + 1] > defects[i])
if (meanDefects[i + 1] > meanDefects[i])
{
stepSize = stepSize * 0.9;
if (i > 10 && (defects.segment(i - 4, 4).maxCoeff() - defects.segment(i - 4, 4).minCoeff()) / defects[i] < 0.001) // settled down
if (i > 10 && (meanDefects.segment(i - 4, 4).maxCoeff() - meanDefects.segment(i - 4, 4).minCoeff()) / meanDefects[i] < 0.001) // settled down
break;
if (stepSize < minStepSize)
......@@ -127,25 +131,29 @@ void Refinement::refineFixedLattice(Lattice& lattice, const Matrix3Xf& ucsDirect
lattice = Lattice(basis);
}
double Refinement::getDefect(const Matrix3f& basis, const Matrix3Xf& ucsDirections, const Matrix2Xf& ucsBorderNorms)
int Refinement::getFittedPeaksCount(Lattice& lattice, const Eigen::Matrix3Xf& ucsDirections, const Eigen::Matrix2Xf& ucsBorderNorms)
{
getDefects(defects, lattice.getBasis(), ucsDirections, ucsBorderNorms);
getDefects(defects, lattice.getBasis(), ucsDirections, ucsBorderNorms);
int fittedPeaksCount = (defects < tolerance).count();
return fittedPeaksCount;
}
void Refinement::getDefects(ArrayXf& defects, const Matrix3f& basis, const Matrix3Xf& ucsDirections, const Matrix2Xf& ucsBorderNorms)
{
Matrix2Xf ucsBorderNormsSquared = ucsBorderNorms.array().square().matrix();
ucsBorderNormsSquared = ucsBorderNorms.array().square().matrix();
Matrix3f basis_inverse = basis.inverse();
Matrix3Xf millerIndices_close = basis_inverse * (ucsDirections.array().rowwise() * ucsBorderNorms.array().row(0)).matrix();
millerIndices_close = basis_inverse * (ucsDirections.array().rowwise() * ucsBorderNorms.array().row(0)).matrix();
roundTowardsZero(millerIndices_close);
Matrix3Xf millerIndices_far = basis_inverse * (ucsDirections.array().rowwise() * ucsBorderNorms.array().row(1)).matrix();
millerIndices_far = basis_inverse * (ucsDirections.array().rowwise() * ucsBorderNorms.array().row(1)).matrix();
roundAwayFromZero(millerIndices_far);
Matrix3Xf candidatePeaks;
RowVectorXf candidatePeaksNormsSquared;
vector<Vector3f> millerIndices;
millerIndices.reserve(500);
int peakCount = ucsDirections.cols();
RowVectorXf projectedVectorNorms;
Matrix3Xf defectVectors;
ArrayXf defects(peakCount);
int notPredictablePeaksCount = 0;
defects.resize(peakCount);
for (int i = 0; i < peakCount; i++)
{
// create miller indices close to UCS
......@@ -177,20 +185,29 @@ double Refinement::getDefect(const Matrix3f& basis, const Matrix3Xf& ucsDirectio
}
if (validCandidatePeaksCount == 0)
{
notPredictablePeaksCount++;
defects(i) = numeric_limits<float>::max();
defects(i) = 1;
continue;
}
projectedVectorNorms.noalias() = candidatePeaks.leftCols(validCandidatePeaksCount).transpose() * ucsDirections.col(i);
defectVectors = candidatePeaks.leftCols(validCandidatePeaksCount);
defectVectors.noalias() -= ucsDirections.col(i) * projectedVectorNorms;
defects(i) = defectVectors.colwise().norm().minCoeff();
defectVectors_absolute = candidatePeaks.leftCols(validCandidatePeaksCount);
defectVectors_absolute.noalias() -= ucsDirections.col(i) * projectedVectorNorms;
defectVectors_relative.noalias() = basis_inverse * defectVectors_absolute;
defects(i) = defectVectors_relative.cwiseAbs().colwise().maxCoeff().minCoeff();
}
}
double Refinement::getMeanDefect(const Matrix3f& basis, const Matrix3Xf& ucsDirections, const Matrix2Xf& ucsBorderNorms)
{
getDefects(defects, basis, ucsDirections, ucsBorderNorms);
notPredictablePeaks = defects > tolerance;
int16_t notPredictablePeaksCount = notPredictablePeaks.count();
// cout << "np " << notPredictablePeaksCount << endl;
if (defects.size() == notPredictablePeaksCount)
if (notPredictablePeaksCount == defects.size())
{
return (basis.colwise().norm()).maxCoeff();
return 1;
}
sort((float*)defects.data(), (float*)defects.data() + defects.size());
......
......@@ -153,7 +153,7 @@ void Sinogram::computePartOfSinogramOnePeak(Matrix3Xf* candidateReflectionDirect
temp2 = sinah * (-l->dot(m));
gammah = acos(temp2);
temp1 = m * cosah.matrix();
temp1.noalias() = m * cosah.matrix();
temp1.noalias() += l->cross(m) * sinah.matrix();
n = temp1.array().rowwise() * rsqrt(1 - temp2.square());
// gammah = acos(sinah * (-l->dot(m)));
......
......@@ -83,8 +83,8 @@ extern "C" void PinkIndexer_delete(PinkIndexer* pinkIndexer)
delete 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)
extern "C" int PinkIndexer_indexPattern(PinkIndexer* pinkIndexer, Lattice_t* indexedLattice, const reciprocalPeaks_1_per_A_t meanReciprocalPeaks_1_per_A,
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++)
......@@ -95,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, intensities_class, meanReciprocalPeaks_1_per_A_matrix, threadCount);
bool indexed = pinkIndexer->indexPattern(indexedLattice_class, intensities_class, meanReciprocalPeaks_1_per_A_matrix, threadCount);
Eigen::Matrix3f basis = indexedLattice_class.getBasis();
indexedLattice->ax = basis(0, 0);
......@@ -107,4 +107,6 @@ extern "C" void PinkIndexer_indexPattern(PinkIndexer* pinkIndexer, Lattice_t* in
indexedLattice->cx = basis(0, 2);
indexedLattice->cy = basis(1, 2);
indexedLattice->cz = basis(2, 2);
return (int)indexed;
}
\ No newline at end of file
......@@ -21,11 +21,11 @@ int main()
// testReflectionsInRangeFinder();
// testSinogram();
// testSinogram2();
//testSinogramComplete();
// testSinogramComplete();
// testRefinementGetDefect();
// testRefinement();
testRefinement();
// testPatternPrediction();
testPinkIndexer();
// testPinkIndexer();
}
catch (exception& e)
{
......
......@@ -55,7 +55,7 @@ void testPinkIndexer()
{
auto tmp = Chronometer("sinogram");
int threadCount = 4;
int threadCount = 6;
pinkIndexer.indexPattern(indexedLattice, intensities, peaksOnDetector_m, threadCount);
}
......@@ -101,8 +101,9 @@ void testRefinement()
loadEigenMatrixFromDisk(ucsDirections, "C:\\DesyFiles\\workspaces\\VisualStudio_workspace\\pinkIndexer\\workfolder\\ucsDirections");
loadEigenMatrixFromDisk(basis, "C:\\DesyFiles\\workspaces\\VisualStudio_workspace\\pinkIndexer\\workfolder\\bestFitBasis");
float tolerance = 0.15;
Lattice lattice(basis);
Refinement refinement;
Refinement refinement(tolerance);
cout << lattice << endl;
refinement.refineFixedLattice(lattice, ucsDirections, ucsBorderNorms);
......
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