Commit 6f6298b8 authored by Yaroslav Gevorkov's avatar Yaroslav Gevorkov
Browse files

added and tested PinkIndexer class, refactoring, added crystfel integration, not tested yet

parent 4bf2073d
......@@ -14,9 +14,7 @@ endif()
find_package(Eigen3 3.3 NO_MODULE)
include_directories(include
include/iniReader
)
include_directories(include)
set(SOURCES src/Chronometer.cpp
src/ReciprocalToRealProjection.cpp
......@@ -28,6 +26,13 @@ set(SOURCES src/Chronometer.cpp
src/Sinogram.cpp
src/Refinement.cpp
src/SimpleDiffractionPatternPrediction.cpp
src/PinkIndexer.cpp
src/adaptions/crystfel/ExperimentSettings.cpp
src/adaptions/crystfel/PinkIndexer.cpp
src/adaptions/crystfel/SimpleDiffractionPatternPrediction.cpp
src/adaptions/crystfel/projectionData.cpp
src/adaptions/crystfel/SimpleProjection.cpp
)
set(SOURCES_test src/main.cpp
......
#pragma once
#include "Backprojection.h"
#include "ExperimentSettings.h"
#include "Lattice.h"
#include "Refinement.h"
#include "SimpleDiffractionPatternPrediction.h"
#include "SimpleProjection.h"
#include "Sinogram.h"
#include <Eigen/Dense>
#include <stdint.h>
class PinkIndexer
{
public:
enum class ConsideredPeaksCount
{
veryFew,
few,
standard,
many,
manyMany
};
enum class AngleResolution
{
extremelyLoose,
loose,
standard,
dense,
extremelyDense
};
enum class RefinementType
{
none,
fixedLatticeParameters,
variableLatticeParameters,
firstFixedThenVariableLatticeParameters
};
PinkIndexer(const ExperimentSettings& experimentSettings, ConsideredPeaksCount consideredPeaksCount, AngleResolution angleResolution,
RefinementType refinementType);
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);
private:
void raducePeakCount(Eigen::Matrix3Xf& ucsDirections, Eigen::Array2Xf& ucsBorderNorms, const Eigen::Matrix2Xf& detectorPeaks_m,
const Eigen::RowVectorXf& intensities);
float getAngleResolution();
int getConsideredPeaksCount();
SimpleProjection reciprocalToRealProjection;
Backprojection backprojection;
Sinogram sinogram;
Refinement refinement;
Lattice sampleLattice;
ConsideredPeaksCount consideredPeaksCount;
AngleResolution angleResolution;
RefinementType refinementType;
};
\ No newline at end of file
#pragma once
#include <Eigen/Dense>
#include "ExperimentSettings.h"
#include <Eigen/Dense>
class ReciprocalToRealProjection
{
public:
ReciprocalToRealProjection(const ExperimentSettings& experimentSettings);
virtual ~ReciprocalToRealProjection() = default;
public:
ReciprocalToRealProjection(const ExperimentSettings& experimentSettings);
virtual ~ReciprocalToRealProjection() = default;
virtual void project(const Eigen::Matrix3Xf& reciprocalPoints, Eigen::Matrix2Xf& projectedPoints) = 0;
virtual void project(Eigen::Matrix2Xf& projectedPoints, const Eigen::Matrix3Xf& reciprocalPoints) = 0;
protected:
ExperimentSettings experimentSettings;
protected:
ExperimentSettings experimentSettings;
};
......@@ -8,5 +8,5 @@ class SimpleProjection : public ReciprocalToRealProjection
public:
SimpleProjection(const ExperimentSettings& experimentSettings);
void project(const Eigen::Matrix3Xf& reciprocalPeaks, Eigen::Matrix2Xf& projectedPeaks);
void project(Eigen::Matrix2Xf& projectedPeaks, const Eigen::Matrix3Xf& reciprocalPeaks);
};
......@@ -16,9 +16,9 @@ class Sinogram
void setSinogramAngleResolution(float angleResolution_deg);
void computeSinogram(const Eigen::Matrix3Xf& ucsDirections, EigenSTL::vector_Matrix3Xf& candidateReflectionsDirections);
void computeSinogramParallel(const Eigen::Matrix3Xf& ucsDirections, EigenSTL::vector_Matrix3Xf& candidateReflectionsDirections, int slaveThreadCount);
void computeSinogramParallel2(const Eigen::Matrix3Xf& ucsDirections, EigenSTL::vector_Matrix3Xf& candidateReflectionsDirections, int slaveThreadCount);
void computeSinogram(const Eigen::Matrix3Xf& ucsDirections, const Eigen::Matrix2Xf ucsBorderNorms);
void computeSinogramParallel(const Eigen::Matrix3Xf& ucsDirections, const Eigen::Matrix2Xf ucsBorderNorms, int slaveThreadCount);
void computeSinogramParallel2(const Eigen::Matrix3Xf& ucsDirections, const Eigen::Matrix2Xf ucsBorderNorms, int slaveThreadCount);
void getBestRotation(Eigen::AngleAxisf& bestRotation);
......
#ifndef ADAPTIONS_CRYSTFEL_EXPERIMENT_SETTINGS_H
#define ADAPTIONS_CRYSTFEL_EXPERIMENT_SETTINGS_H
#include "Lattice.h"
#ifdef __cplusplus
extern "C" {
#endif
typedef struct ExperimentSettings ExperimentSettings;
ExperimentSettings* ExperimentSettings_new_nolatt(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* ExperimentSettings_new(float beamEenergy_eV, float detectorDistance_m, float detectorRadius_m, float divergenceAngle_deg,
float nonMonochromaticity, const Lattice_t sampleReciprocalLattice_1A, float tolerance,
float reflectionRadius_1_per_A);
void ExperimentSettings_delete(ExperimentSettings* experimentSettings);
#ifdef __cplusplus
}
#endif
#endif
\ No newline at end of file
#ifndef ADAPTIONS_CRYSTFEL_LATTICE_H
#define ADAPTIONS_CRYSTFEL_LATTICE_H
typedef struct
{
float ax;
float ay;
float az;
float bx;
float by;
float bz;
float cx;
float cy;
float cz;
} Lattice_t;
#endif
\ No newline at end of file
#ifndef ADAPTIONS_CRYSTFEL_PINK_INDEXER_H
#define ADAPTIONS_CRYSTFEL_PINK_INDEXER_H
#include "ExperimentSettings.h"
#include "indexerData.h"
#ifdef __cplusplus
extern "C" {
#endif
typedef enum {
CONSIDERED_PEAKS_COUNT_veryFew = 0,
CONSIDERED_PEAKS_COUNT_few = 1,
CONSIDERED_PEAKS_COUNT_standard = 2,
CONSIDERED_PEAKS_COUNT_many = 3,
CONSIDERED_PEAKS_COUNT_manyMany = 4,
CONSIDERED_PEAKS_COUNT_lastEnum
} consideredPeaksCount_t;
typedef enum {
ANGLE_RESOLUTION_extremelyLoose = 0,
ANGLE_RESOLUTION_loose = 1,
ANGLE_RESOLUTION_standard = 2,
ANGLE_RESOLUTION_dense = 3,
ANGLE_RESOLUTION_extremelyDense = 4,
ANGLE_RESOLUTION_lastEnum
} angleResolution_t;
typedef enum {
REFINEMENT_TYPE_none = 0,
REFINEMENT_TYPE_fixedLatticeParameters = 1,
REFINEMENT_TYPE_variableLatticeParameters = 2,
REFINEMENT_TYPE_firstFixedThenVariableLatticeParameters = 3,
REFINEMENT_TYPE_lastEnum
} refinementType_t;
typedef struct PinkIndexer PinkIndexer;
PinkIndexer* PinkIndexer_new(ExperimentSettings* experimentSettings, consideredPeaksCount_t consideredPeaksCount, angleResolution_t angleResolution,
refinementType_t refinementType);
void PinkIndexer_delete(PinkIndexer* pinkIndexer);
void 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
}
#endif
#endif
\ No newline at end of file
#ifndef ADAPTIONS_CRYSTFEL_SIMPLE_DIFFRACTION_PATTERN_PREDICTION_H
#define ADAPTIONS_CRYSTFEL_SIMPLE_DIFFRACTION_PATTERN_PREDICTION_H
#include "ExperimentSettings.h"
#include "indexerData.h"
#include "projectionData.h"
#ifdef __cplusplus
extern "C" {
#endif
typedef struct SimpleDiffractionPatternPrediction SimpleDiffractionPatternPrediction;
SimpleDiffractionPatternPrediction* SimpleDiffractionPrediction_new(ExperimentSettings* experimentSettings);
void SimpleDiffractionPatternPrediction_delete(SimpleDiffractionPatternPrediction* simpleDiffractionPatternPrediction);
void getPeaksOnEwaldSphere(SimpleDiffractionPatternPrediction* simpleDiffractionPatternPrediction, reciprocalPeaks_1_per_A_t reciprocalPeaks_1_per_A,
Lattice_t lattice);
void predictPattern(SimpleDiffractionPatternPrediction* simpleDiffractionPatternPrediction, detectorPeaks_m_t predictedPeaks, Lattice_t lattice);
#ifdef __cplusplus
}
#endif
#endif
#ifndef ADAPTIONS_CRYSTFEL_SIMPLE_PROJECTION_H
#define ADAPTIONS_CRYSTFEL_SIMPLE_PROJECTION_H
#include "ExperimentSettings.h"
#include "indexerData.h"
#include "projectionData.h"
#ifdef __cplusplus
extern "C" {
#endif
typedef struct SimpleProjection SimpleProjection;
SimpleProjection* SimpleProjection_new(ExperimentSettings* experimentSettings);
void SimpleProjection_delete(SimpleProjection* simpleProjection);
void project(SimpleProjection* simpleProjection, detectorPeaks_m_t projectedPeaks_m, const reciprocalPeaks_1_per_A_t reciprocalPeaks_1_per_A);
#ifdef __cplusplus
}
#endif
#endif
#ifndef ADAPTIONS_CRYSTFEL_INDEXER_DATA_H
#define ADAPTIONS_CRYSTFEL_INDEXER_DATA_H
#define MAX_PEAK_COUNT_FOR_INDEXER 20000
typedef struct
{
float* coordinates_x;
float* coordinates_y;
float* coordinates_z;
int peakCount;
} reciprocalPeaks_1_per_A_t;
#ifdef __cplusplus
extern "C" {
#endif
void allocReciprocalPeaks(reciprocalPeaks_1_per_A_t* reciprocalPeaks_1_per_A);
void freeReciprocalPeaks(reciprocalPeaks_1_per_A_t reciprocalPeaks_1_per_A);
#ifdef __cplusplus
}
#endif
#endif
\ No newline at end of file
#ifndef ADAPTIONS_CRYSTFEL_PROJECTION_DATA_H
#define ADAPTIONS_CRYSTFEL_PROJECTION_DATA_H
#define MAX_PEAK_COUNT_FOR_PROJECTION 40000
typedef struct
{
float* coordinates_x;
float* coordinates_y;
int peakCount;
} detectorPeaks_m_t;
#ifdef __cplusplus
extern "C" {
#endif
void allocDetectorPeaks(detectorPeaks_m_t* detectorPeaks_m);
void freeDetectorPeaks(detectorPeaks_m_t detectorPeaks_m);
#ifdef __cplusplus
}
#endif
#endif
// Read an INI file into easy-to-access name/value pairs.
// inih and INIReader are released under the New BSD license (see LICENSE.txt).
// Go to the project home page for more info:
//
// https://github.com/benhoyt/inih
/* inih -- simple .INI file parser
inih is released under the New BSD license (see LICENSE.txt). Go to the project
home page for more info:
https://github.com/benhoyt/inih
*/
#ifndef __INI_H__
#define __INI_H__
/* Make this header file easier to include in C++ code */
#ifdef __cplusplus
extern "C" {
#endif
#include <stdio.h>
/* Typedef for prototype of handler function. */
typedef int (*ini_handler)(void* user, const char* section,
const char* name, const char* value);
/* Typedef for prototype of fgets-style reader function. */
typedef char* (*ini_reader)(char* str, int num, void* stream);
/* Parse given INI-style file. May have [section]s, name=value pairs
(whitespace stripped), and comments starting with ';' (semicolon). Section
is "" if name=value pair parsed before any section heading. name:value
pairs are also supported as a concession to Python's configparser.
For each name=value pair parsed, call handler function with given user
pointer as well as section, name, and value (data only valid for duration
of handler call). Handler should return nonzero on success, zero on error.
Returns 0 on success, line number of first error on parse error (doesn't
stop on first error), -1 on file open error, or -2 on memory allocation
error (only when INI_USE_STACK is zero).
*/
int ini_parse(const char* filename, ini_handler handler, void* user);
/* Same as ini_parse(), but takes a FILE* instead of filename. This doesn't
close the file when it's finished -- the caller must do that. */
int ini_parse_file(FILE* file, ini_handler handler, void* user);
/* Same as ini_parse(), but takes an ini_reader function pointer instead of
filename. Used for implementing custom or string-based I/O. */
int ini_parse_stream(ini_reader reader, void* stream, ini_handler handler,
void* user);
/* Nonzero to allow multi-line value parsing, in the style of Python's
configparser. If allowed, ini_parse() will call the handler with the same
name for each subsequent line parsed. */
#ifndef INI_ALLOW_MULTILINE
#define INI_ALLOW_MULTILINE 1
#endif
/* Nonzero to allow a UTF-8 BOM sequence (0xEF 0xBB 0xBF) at the start of
the file. See http://code.google.com/p/inih/issues/detail?id=21 */
#ifndef INI_ALLOW_BOM
#define INI_ALLOW_BOM 1
#endif
/* Nonzero to allow inline comments (with valid inline comment characters
specified by INI_INLINE_COMMENT_PREFIXES). Set to 0 to turn off and match
Python 3.2+ configparser behaviour. */
#ifndef INI_ALLOW_INLINE_COMMENTS
#define INI_ALLOW_INLINE_COMMENTS 1
#endif
#ifndef INI_INLINE_COMMENT_PREFIXES
#define INI_INLINE_COMMENT_PREFIXES ";"
#endif
/* Nonzero to use stack, zero to use heap (malloc/free). */
#ifndef INI_USE_STACK
#define INI_USE_STACK 1
#endif
/* Stop parsing on first error (default is to keep parsing). */
#ifndef INI_STOP_ON_FIRST_ERROR
#define INI_STOP_ON_FIRST_ERROR 0
#endif
/* Maximum line length for any line in INI file. */
#ifndef INI_MAX_LINE
#define INI_MAX_LINE 200
#endif
#ifdef __cplusplus
}
#endif
/* inih -- simple .INI file parser
inih is released under the New BSD license (see LICENSE.txt). Go to the project
home page for more info:
https://github.com/benhoyt/inih
*/
#if defined(_MSC_VER) && !defined(_CRT_SECURE_NO_WARNINGS)
#define _CRT_SECURE_NO_WARNINGS
#endif
#include <stdio.h>
#include <ctype.h>
#include <string.h>
#if !INI_USE_STACK
#include <stdlib.h>
#endif
#define MAX_SECTION 50
#define MAX_NAME 50
/* Strip whitespace chars off end of given string, in place. Return s. */
inline static char* rstrip(char* s)
{
char* p = s + strlen(s);
while (p > s && isspace((unsigned char)(*--p)))
*p = '\0';
return s;
}
/* Return pointer to first non-whitespace char in given string. */
inline static char* lskip(const char* s)
{
while (*s && isspace((unsigned char)(*s)))
s++;
return (char*)s;
}
/* Return pointer to first char (of chars) or inline comment in given string,
or pointer to null at end of string if neither found. Inline comment must
be prefixed by a whitespace character to register as a comment. */
inline static char* find_chars_or_comment(const char* s, const char* chars)
{
#if INI_ALLOW_INLINE_COMMENTS
int was_space = 0;
while (*s && (!chars || !strchr(chars, *s)) &&
!(was_space && strchr(INI_INLINE_COMMENT_PREFIXES, *s))) {
was_space = isspace((unsigned char)(*s));
s++;
}
#else
while (*s && (!chars || !strchr(chars, *s))) {
s++;
}
#endif
return (char*)s;
}
/* Version of strncpy that ensures dest (size bytes) is null-terminated. */
inline static char* strncpy0(char* dest, const char* src, size_t size)
{
strncpy(dest, src, size);
dest[size - 1] = '\0';
return dest;
}
/* See documentation in header file. */
inline int ini_parse_stream(ini_reader reader, void* stream, ini_handler handler,
void* user)
{
/* Uses a fair bit of stack (use heap instead if you need to) */
#if INI_USE_STACK
char line[INI_MAX_LINE];
#else
char* line;
#endif
char section[MAX_SECTION] = "";
char prev_name[MAX_NAME] = "";
char* start;
char* end;
char* name;
char* value;
int lineno = 0;
int error = 0;
#if !INI_USE_STACK
line = (char*)malloc(INI_MAX_LINE);
if (!line) {
return -2;
}
#endif
/* Scan through stream line by line */
while (reader(line, INI_MAX_LINE, stream) != NULL) {
lineno++;
start = line;
#if INI_ALLOW_BOM
if (lineno == 1 && (unsigned char)start[0] == 0xEF &&
(unsigned char)start[1] == 0xBB &&
(unsigned char)start[2] == 0xBF) {
start += 3;
}
#endif
start = lskip(rstrip(start));
if (*start == ';' || *start == '#') {
/* Per Python configparser, allow both ; and # comments at the
start of a line */
}
#if INI_ALLOW_MULTILINE
else if (*prev_name && *start && start > line) {
#if INI_ALLOW_INLINE_COMMENTS
end = find_chars_or_comment(start, NULL);
if (*end)
*end = '\0';
rstrip(start);
#endif
/* Non-blank line with leading whitespace, treat as continuation
of previous name's value (as per Python configparser). */
if (!handler(user, section, prev_name, start) && !error)
error = lineno;
}
#endif
else if (*start == '[') {
/* A "[section]" line */
end = find_chars_or_comment(start + 1, "]");
if (*end == ']') {
*end = '\0';
strncpy0(section, start + 1, sizeof(section));
*prev_name = '\0';
}
else if (!error) {
/* No ']' found on section line */
error = lineno;
}
}
else if (*start) {
/* Not a comment, must be a name[=:]value pair */
end = find_chars_or_comment(start, "=:");
if (*end == '=' || *end == ':') {
*end = '\0';
name = rstrip(start);
value = lskip(end + 1);
#if INI_ALLOW_INLINE_COMMENTS
end = find_chars_or_comment(value, NULL);
if (*end)
*end = '\0';
#endif
rstrip(value);
/* Valid name[=:]value pair found, call handler */
strncpy0(prev_name, name, sizeof(prev_name));