Merge branch '2404-indexing' into 'main'

Indexing improvements

See merge request jungfraujoch/nextgendcu!47
This commit is contained in:
2024-04-20 13:41:41 +02:00
56 changed files with 742 additions and 187 deletions

View File

@@ -21,6 +21,7 @@ inline SpotFindingSettings Convert(const org::openapitools::server::model::Spot_
ret.low_resolution_limit = input.getLowResolutionLimit();
ret.enable = input.isEnable();
ret.indexing = input.isIndexing();
ret.indexing_tolerance = input.getIndexingTolerance();
return ret;
}
@@ -34,6 +35,7 @@ inline org::openapitools::server::model::Spot_finding_settings Convert(const Spo
ret.setLowResolutionLimit(input.low_resolution_limit);
ret.setEnable(input.enable);
ret.setIndexing(input.indexing);
ret.setIndexingTolerance(input.indexing_tolerance);
return ret;
}
@@ -47,9 +49,14 @@ inline org::openapitools::server::model::Measurement_statistics Convert(const Me
ret.setImagesCollected(input.images_collected);
ret.setImagesSent(input.images_sent);
ret.setMaxImageNumberSent(input.max_image_number_sent);
if (input.collection_efficiency >= 0.0)
ret.setCollectionEfficiency(input.collection_efficiency);
ret.setCompressionRatio(input.compression_ratio);
if (input.collection_efficiency)
ret.setCollectionEfficiency(input.collection_efficiency.value());
if (input.compression_ratio)
ret.setCompressionRatio(input.compression_ratio.value());
if (input.beam_center_drift_pxl) {
ret.setBeamCenterDriftXPxl(input.beam_center_drift_pxl->first);
ret.setBeamCenterDriftYPxl(input.beam_center_drift_pxl->second);
}
ret.setCancelled(input.cancelled);
ret.setMaxReceiverDelay(input.max_receive_delay);
@@ -58,11 +65,11 @@ inline org::openapitools::server::model::Measurement_statistics Convert(const Me
ret.setDetectorHeight(input.detector_height);
ret.setDetectorPixelDepth(input.detector_pixel_depth);
if (input.indexing_rate >= 0.0)
ret.setIndexingRate(input.indexing_rate);
if (input.indexing_rate)
ret.setIndexingRate(input.indexing_rate.value());
if (input.bkg_estimate >= 0.0)
ret.setBkgEstimate(input.bkg_estimate);
if (input.bkg_estimate)
ret.setBkgEstimate(input.bkg_estimate.value());
return ret;
}
@@ -404,6 +411,10 @@ inline DatasetSettings Convert(const org::openapitools::server::model::Dataset_s
ret.ImageAppendix(input.getImageAppendix());
ret.SaveCalibration(input.isSaveCalibration());
ret.ImagesPerFile(input.getImagesPerFile());
if (input.dataReductionFactorSerialmxIsSet())
ret.DataReductionFactorSerialMX(input.getDataReductionFactorSerialmx());
return ret;
}
@@ -814,3 +825,12 @@ void JFJochBrokerHttp::plot_roi_valid_pixels_post(const org::openapitools::serve
auto plot = state_machine.GetPlots(req);
ProcessOutput(Convert(plot), response);
}
void JFJochBrokerHttp::plot_beam_center_drift_post(const org::openapitools::server::model::Plot_request &plotRequest,
Pistache::Http::ResponseWriter &response) {
PlotRequest req{.type = PlotType::IndexingDrift, .binning = 0};
if (plotRequest.binningIsSet())
req.binning = plotRequest.getBinning();
auto plot = state_machine.GetPlots(req);
ProcessOutput(Convert(plot), response);
}

View File

@@ -51,6 +51,9 @@ class JFJochBrokerHttp : public org::openapitools::server::api::DefaultApi {
void plot_receiver_free_send_buffers_post(const org::openapitools::server::model::Plot_request &plotRequest,
Pistache::Http::ResponseWriter &response) override;
void plot_beam_center_drift_post(const org::openapitools::server::model::Plot_request &plotRequest,
Pistache::Http::ResponseWriter &response) override;
void plot_receiver_delay_post(const org::openapitools::server::model::Plot_request &plotRequest,
Pistache::Http::ResponseWriter &response) override;

View File

@@ -68,10 +68,16 @@ JFJochServicesOutput JFJochServices::Stop() {
logger.Info("Wait for receiver done");
ret.receiver_output = receiver->Stop();
if (ret.receiver_output.status.compressed_ratio)
logger.Info(" ... Receiver efficiency: {} % Max delay: {} Compression ratio {}x",
static_cast<int>(ret.receiver_output.efficiency * 100.0),
ret.receiver_output.status.max_receive_delay,
static_cast<int>(std::round(ret.receiver_output.status.compressed_ratio)));
static_cast<int>(std::round(ret.receiver_output.status.compressed_ratio.value())));
else
logger.Info(" ... Receiver efficiency: {} % Max delay: {}",
static_cast<int>(ret.receiver_output.efficiency * 100.0),
ret.receiver_output.status.max_receive_delay);
if (ret.receiver_output.efficiency < 1.0) {
for (int i = 0; i < ret.receiver_output.received_packets.size(); i++) {
if (ret.receiver_output.received_packets[i] != ret.receiver_output.expected_packets[i])

View File

@@ -411,6 +411,7 @@ void JFJochStateMachine::SetFullMeasurementOutput(const JFJochServicesOutput &ou
tmp.max_receive_delay = output.receiver_output.status.max_receive_delay;
tmp.indexing_rate = output.receiver_output.status.indexing_rate;
tmp.bkg_estimate = output.receiver_output.status.bkg_estimate;
tmp.beam_center_drift_pxl = output.receiver_output.status.beam_center_drift_pxl;
measurement_statistics = tmp;
}
@@ -445,7 +446,7 @@ std::optional<MeasurementStatistics> JFJochStateMachine::GetMeasurementStatistic
tmp.images_expected = experiment.GetImageNum();
tmp.compression_ratio = rcv_status->compressed_ratio;
tmp.collection_efficiency = -1.0;
tmp.beam_center_drift_pxl = rcv_status->beam_center_drift_pxl;
tmp.images_collected = rcv_status->images_collected;
tmp.images_sent = rcv_status->images_sent;
tmp.cancelled = rcv_status->cancelled;

View File

@@ -41,19 +41,20 @@ struct MeasurementStatistics {
int64_t images_collected;
int64_t images_sent;
int64_t max_image_number_sent;
float collection_efficiency;
float compression_ratio;
std::optional<float> collection_efficiency;
std::optional<float> compression_ratio;
bool cancelled;
int64_t max_receive_delay;
float indexing_rate;
std::optional<float> indexing_rate;
int64_t detector_width;
int64_t detector_height;
int64_t detector_pixel_depth;
float bkg_estimate;
std::optional<float> bkg_estimate;
std::optional<std::pair<float, float>> beam_center_drift_pxl;
};
struct DetectorSettings {

View File

@@ -47,6 +47,7 @@ void DefaultApi::setupRoutes() {
Routes::Get(*router, base + "/detector/status", Routes::bind(&DefaultApi::detector_status_get_handler, this));
Routes::Post(*router, base + "/initialize", Routes::bind(&DefaultApi::initialize_post_handler, this));
Routes::Post(*router, base + "/pedestal", Routes::bind(&DefaultApi::pedestal_post_handler, this));
Routes::Post(*router, base + "/plot/beam_center_drift", Routes::bind(&DefaultApi::plot_beam_center_drift_post_handler, this));
Routes::Post(*router, base + "/plot/bkg_estimate", Routes::bind(&DefaultApi::plot_bkg_estimate_post_handler, this));
Routes::Post(*router, base + "/plot/error_pixel", Routes::bind(&DefaultApi::plot_error_pixel_post_handler, this));
Routes::Post(*router, base + "/plot/image_collection_efficiency", Routes::bind(&DefaultApi::plot_image_collection_efficiency_post_handler, this));
@@ -431,6 +432,39 @@ void DefaultApi::pedestal_post_handler(const Pistache::Rest::Request &, Pistache
response.send(Pistache::Http::Code::Internal_Server_Error, e.what());
}
}
void DefaultApi::plot_beam_center_drift_post_handler(const Pistache::Rest::Request &request, Pistache::Http::ResponseWriter response) {
try {
// Getting the body param
Plot_request plotRequest;
try {
nlohmann::json::parse(request.body()).get_to(plotRequest);
plotRequest.validate();
} catch (std::exception &e) {
const std::pair<Pistache::Http::Code, std::string> errorInfo = this->handleParsingException(e);
response.send(errorInfo.first, errorInfo.second);
return;
}
try {
this->plot_beam_center_drift_post(plotRequest, response);
} catch (Pistache::Http::HttpError &e) {
response.send(static_cast<Pistache::Http::Code>(e.code()), e.what());
return;
} catch (std::exception &e) {
const std::pair<Pistache::Http::Code, std::string> errorInfo = this->handleOperationException(e);
response.send(errorInfo.first, errorInfo.second);
return;
}
} catch (std::exception &e) {
response.send(Pistache::Http::Code::Internal_Server_Error, e.what());
}
}
void DefaultApi::plot_bkg_estimate_post_handler(const Pistache::Rest::Request &request, Pistache::Http::ResponseWriter response) {
try {

View File

@@ -72,6 +72,7 @@ private:
void detector_status_get_handler(const Pistache::Rest::Request &request, Pistache::Http::ResponseWriter response);
void initialize_post_handler(const Pistache::Rest::Request &request, Pistache::Http::ResponseWriter response);
void pedestal_post_handler(const Pistache::Rest::Request &request, Pistache::Http::ResponseWriter response);
void plot_beam_center_drift_post_handler(const Pistache::Rest::Request &request, Pistache::Http::ResponseWriter response);
void plot_bkg_estimate_post_handler(const Pistache::Rest::Request &request, Pistache::Http::ResponseWriter response);
void plot_error_pixel_post_handler(const Pistache::Rest::Request &request, Pistache::Http::ResponseWriter response);
void plot_image_collection_efficiency_post_handler(const Pistache::Rest::Request &request, Pistache::Http::ResponseWriter response);
@@ -224,6 +225,14 @@ private:
/// </remarks>
virtual void pedestal_post(Pistache::Http::ResponseWriter &response) = 0;
/// <summary>
/// Generate plot of beam center drift (from indexing)
/// </summary>
/// <remarks>
/// Max count of beam center drift per image; binning is configurable
/// </remarks>
/// <param name="plotRequest"> (optional)</param>
virtual void plot_beam_center_drift_post(const org::openapitools::server::model::Plot_request &plotRequest, Pistache::Http::ResponseWriter &response) = 0;
/// <summary>
/// Generate background estimate plot
/// </summary>
/// <remarks>

View File

@@ -55,6 +55,8 @@ Dataset_settings::Dataset_settings()
m_Image_appendixIsSet = false;
m_Photon_energy_multiplier = 1.0f;
m_Photon_energy_multiplierIsSet = false;
m_Data_reduction_factor_serialmx = 1.0f;
m_Data_reduction_factor_serialmxIsSet = false;
m_Unit_cellIsSet = false;
}
@@ -224,6 +226,25 @@ bool Dataset_settings::validate(std::stringstream& msg, const std::string& pathP
}
}
if (dataReductionFactorSerialmxIsSet())
{
const float& value = m_Data_reduction_factor_serialmx;
const std::string currentValuePath = _pathPrefix + ".dataReductionFactorSerialmx";
if (value < static_cast<float>(0.0))
{
success = false;
msg << currentValuePath << ": must be greater than or equal to 0.0;";
}
if (value > static_cast<float>(1.0))
{
success = false;
msg << currentValuePath << ": must be less than or equal to 1.0;";
}
}
return success;
}
@@ -294,6 +315,9 @@ bool Dataset_settings::operator==(const Dataset_settings& rhs) const
((!photonEnergyMultiplierIsSet() && !rhs.photonEnergyMultiplierIsSet()) || (photonEnergyMultiplierIsSet() && rhs.photonEnergyMultiplierIsSet() && getPhotonEnergyMultiplier() == rhs.getPhotonEnergyMultiplier())) &&
((!dataReductionFactorSerialmxIsSet() && !rhs.dataReductionFactorSerialmxIsSet()) || (dataReductionFactorSerialmxIsSet() && rhs.dataReductionFactorSerialmxIsSet() && getDataReductionFactorSerialmx() == rhs.getDataReductionFactorSerialmx())) &&
((!unitCellIsSet() && !rhs.unitCellIsSet()) || (unitCellIsSet() && rhs.unitCellIsSet() && getUnitCell() == rhs.getUnitCell()))
;
@@ -342,6 +366,8 @@ void to_json(nlohmann::json& j, const Dataset_settings& o)
j["image_appendix"] = o.m_Image_appendix;
if(o.photonEnergyMultiplierIsSet())
j["photon_energy_multiplier"] = o.m_Photon_energy_multiplier;
if(o.dataReductionFactorSerialmxIsSet())
j["data_reduction_factor_serialmx"] = o.m_Data_reduction_factor_serialmx;
if(o.unitCellIsSet())
j["unit_cell"] = o.m_Unit_cell;
@@ -429,6 +455,11 @@ void from_json(const nlohmann::json& j, Dataset_settings& o)
j.at("photon_energy_multiplier").get_to(o.m_Photon_energy_multiplier);
o.m_Photon_energy_multiplierIsSet = true;
}
if(j.find("data_reduction_factor_serialmx") != j.end())
{
j.at("data_reduction_factor_serialmx").get_to(o.m_Data_reduction_factor_serialmx);
o.m_Data_reduction_factor_serialmxIsSet = true;
}
if(j.find("unit_cell") != j.end())
{
j.at("unit_cell").get_to(o.m_Unit_cell);
@@ -732,6 +763,23 @@ void Dataset_settings::unsetPhoton_energy_multiplier()
{
m_Photon_energy_multiplierIsSet = false;
}
float Dataset_settings::getDataReductionFactorSerialmx() const
{
return m_Data_reduction_factor_serialmx;
}
void Dataset_settings::setDataReductionFactorSerialmx(float const value)
{
m_Data_reduction_factor_serialmx = value;
m_Data_reduction_factor_serialmxIsSet = true;
}
bool Dataset_settings::dataReductionFactorSerialmxIsSet() const
{
return m_Data_reduction_factor_serialmxIsSet;
}
void Dataset_settings::unsetData_reduction_factor_serialmx()
{
m_Data_reduction_factor_serialmxIsSet = false;
}
org::openapitools::server::model::Dataset_settings_unit_cell Dataset_settings::getUnitCell() const
{
return m_Unit_cell;

View File

@@ -191,6 +191,13 @@ public:
bool photonEnergyMultiplierIsSet() const;
void unsetPhoton_energy_multiplier();
/// <summary>
/// Rate at which non-indexed images are accepted to be forwarded to writer. Value of 1.0 (default) means that all images are written. Values below zero mean that non-indexed images will be accepted with a given probability.
/// </summary>
float getDataReductionFactorSerialmx() const;
void setDataReductionFactorSerialmx(float const value);
bool dataReductionFactorSerialmxIsSet() const;
void unsetData_reduction_factor_serialmx();
/// <summary>
///
/// </summary>
org::openapitools::server::model::Dataset_settings_unit_cell getUnitCell() const;
@@ -241,6 +248,8 @@ protected:
bool m_Image_appendixIsSet;
float m_Photon_energy_multiplier;
bool m_Photon_energy_multiplierIsSet;
float m_Data_reduction_factor_serialmx;
bool m_Data_reduction_factor_serialmxIsSet;
org::openapitools::server::model::Dataset_settings_unit_cell m_Unit_cell;
bool m_Unit_cellIsSet;

View File

@@ -49,6 +49,10 @@ Measurement_statistics::Measurement_statistics()
m_Detector_pixel_depthIsSet = false;
m_Bkg_estimate = 0.0f;
m_Bkg_estimateIsSet = false;
m_Beam_center_drift_x_pxl = 0.0f;
m_Beam_center_drift_x_pxlIsSet = false;
m_Beam_center_drift_y_pxl = 0.0f;
m_Beam_center_drift_y_pxlIsSet = false;
}
@@ -104,7 +108,7 @@ bool Measurement_statistics::validate(std::stringstream& msg, const std::string&
}
}
return success;
}
@@ -153,7 +157,13 @@ bool Measurement_statistics::operator==(const Measurement_statistics& rhs) const
((!detectorPixelDepthIsSet() && !rhs.detectorPixelDepthIsSet()) || (detectorPixelDepthIsSet() && rhs.detectorPixelDepthIsSet() && getDetectorPixelDepth() == rhs.getDetectorPixelDepth())) &&
((!bkgEstimateIsSet() && !rhs.bkgEstimateIsSet()) || (bkgEstimateIsSet() && rhs.bkgEstimateIsSet() && getBkgEstimate() == rhs.getBkgEstimate()))
((!bkgEstimateIsSet() && !rhs.bkgEstimateIsSet()) || (bkgEstimateIsSet() && rhs.bkgEstimateIsSet() && getBkgEstimate() == rhs.getBkgEstimate())) &&
((!beamCenterDriftXPxlIsSet() && !rhs.beamCenterDriftXPxlIsSet()) || (beamCenterDriftXPxlIsSet() && rhs.beamCenterDriftXPxlIsSet() && getBeamCenterDriftXPxl() == rhs.getBeamCenterDriftXPxl())) &&
((!beamCenterDriftYPxlIsSet() && !rhs.beamCenterDriftYPxlIsSet()) || (beamCenterDriftYPxlIsSet() && rhs.beamCenterDriftYPxlIsSet() && getBeamCenterDriftYPxl() == rhs.getBeamCenterDriftYPxl()))
;
}
@@ -194,6 +204,10 @@ void to_json(nlohmann::json& j, const Measurement_statistics& o)
j["detector_pixel_depth"] = o.m_Detector_pixel_depth;
if(o.bkgEstimateIsSet())
j["bkg_estimate"] = o.m_Bkg_estimate;
if(o.beamCenterDriftXPxlIsSet())
j["beam_center_drift_x_pxl"] = o.m_Beam_center_drift_x_pxl;
if(o.beamCenterDriftYPxlIsSet())
j["beam_center_drift_y_pxl"] = o.m_Beam_center_drift_y_pxl;
}
@@ -269,6 +283,16 @@ void from_json(const nlohmann::json& j, Measurement_statistics& o)
j.at("bkg_estimate").get_to(o.m_Bkg_estimate);
o.m_Bkg_estimateIsSet = true;
}
if(j.find("beam_center_drift_x_pxl") != j.end())
{
j.at("beam_center_drift_x_pxl").get_to(o.m_Beam_center_drift_x_pxl);
o.m_Beam_center_drift_x_pxlIsSet = true;
}
if(j.find("beam_center_drift_y_pxl") != j.end())
{
j.at("beam_center_drift_y_pxl").get_to(o.m_Beam_center_drift_y_pxl);
o.m_Beam_center_drift_y_pxlIsSet = true;
}
}
@@ -510,6 +534,40 @@ void Measurement_statistics::unsetBkg_estimate()
{
m_Bkg_estimateIsSet = false;
}
float Measurement_statistics::getBeamCenterDriftXPxl() const
{
return m_Beam_center_drift_x_pxl;
}
void Measurement_statistics::setBeamCenterDriftXPxl(float const value)
{
m_Beam_center_drift_x_pxl = value;
m_Beam_center_drift_x_pxlIsSet = true;
}
bool Measurement_statistics::beamCenterDriftXPxlIsSet() const
{
return m_Beam_center_drift_x_pxlIsSet;
}
void Measurement_statistics::unsetBeam_center_drift_x_pxl()
{
m_Beam_center_drift_x_pxlIsSet = false;
}
float Measurement_statistics::getBeamCenterDriftYPxl() const
{
return m_Beam_center_drift_y_pxl;
}
void Measurement_statistics::setBeamCenterDriftYPxl(float const value)
{
m_Beam_center_drift_y_pxl = value;
m_Beam_center_drift_y_pxlIsSet = true;
}
bool Measurement_statistics::beamCenterDriftYPxlIsSet() const
{
return m_Beam_center_drift_y_pxlIsSet;
}
void Measurement_statistics::unsetBeam_center_drift_y_pxl()
{
m_Beam_center_drift_y_pxlIsSet = false;
}
} // namespace org::openapitools::server::model

View File

@@ -156,6 +156,20 @@ public:
void setBkgEstimate(float const value);
bool bkgEstimateIsSet() const;
void unsetBkg_estimate();
/// <summary>
///
/// </summary>
float getBeamCenterDriftXPxl() const;
void setBeamCenterDriftXPxl(float const value);
bool beamCenterDriftXPxlIsSet() const;
void unsetBeam_center_drift_x_pxl();
/// <summary>
///
/// </summary>
float getBeamCenterDriftYPxl() const;
void setBeamCenterDriftYPxl(float const value);
bool beamCenterDriftYPxlIsSet() const;
void unsetBeam_center_drift_y_pxl();
friend void to_json(nlohmann::json& j, const Measurement_statistics& o);
friend void from_json(const nlohmann::json& j, Measurement_statistics& o);
@@ -188,6 +202,10 @@ protected:
bool m_Detector_pixel_depthIsSet;
float m_Bkg_estimate;
bool m_Bkg_estimateIsSet;
float m_Beam_center_drift_x_pxl;
bool m_Beam_center_drift_x_pxlIsSet;
float m_Beam_center_drift_y_pxl;
bool m_Beam_center_drift_y_pxlIsSet;
};

View File

@@ -29,6 +29,7 @@ Spot_finding_settings::Spot_finding_settings()
m_Max_pix_per_spot = 0L;
m_High_resolution_limit = 0.0f;
m_Low_resolution_limit = 0.0f;
m_Indexing_tolerance = 0.0f;
}
@@ -107,7 +108,26 @@ bool Spot_finding_settings::validate(std::stringstream& msg, const std::string&
}
}
/* Indexing_tolerance */ {
const float& value = m_Indexing_tolerance;
const std::string currentValuePath = _pathPrefix + ".indexingTolerance";
if (value < static_cast<float>(0.0))
{
success = false;
msg << currentValuePath << ": must be greater than or equal to 0.0;";
}
if (value > static_cast<float>(1.0))
{
success = false;
msg << currentValuePath << ": must be less than or equal to 1.0;";
}
}
return success;
}
@@ -138,6 +158,9 @@ bool Spot_finding_settings::operator==(const Spot_finding_settings& rhs) const
&&
(getLowResolutionLimit() == rhs.getLowResolutionLimit())
&&
(getIndexingTolerance() == rhs.getIndexingTolerance())
;
@@ -159,6 +182,7 @@ void to_json(nlohmann::json& j, const Spot_finding_settings& o)
j["max_pix_per_spot"] = o.m_Max_pix_per_spot;
j["high_resolution_limit"] = o.m_High_resolution_limit;
j["low_resolution_limit"] = o.m_Low_resolution_limit;
j["indexing_tolerance"] = o.m_Indexing_tolerance;
}
@@ -172,6 +196,7 @@ void from_json(const nlohmann::json& j, Spot_finding_settings& o)
j.at("max_pix_per_spot").get_to(o.m_Max_pix_per_spot);
j.at("high_resolution_limit").get_to(o.m_High_resolution_limit);
j.at("low_resolution_limit").get_to(o.m_Low_resolution_limit);
j.at("indexing_tolerance").get_to(o.m_Indexing_tolerance);
}
@@ -239,6 +264,14 @@ void Spot_finding_settings::setLowResolutionLimit(float const value)
{
m_Low_resolution_limit = value;
}
float Spot_finding_settings::getIndexingTolerance() const
{
return m_Indexing_tolerance;
}
void Spot_finding_settings::setIndexingTolerance(float const value)
{
m_Indexing_tolerance = value;
}
} // namespace org::openapitools::server::model

View File

@@ -97,6 +97,11 @@ public:
/// </summary>
float getLowResolutionLimit() const;
void setLowResolutionLimit(float const value);
/// <summary>
/// Acceptance tolerance for spots after the indexing run - the larger the number, the more spots will be accepted
/// </summary>
float getIndexingTolerance() const;
void setIndexingTolerance(float const value);
friend void to_json(nlohmann::json& j, const Spot_finding_settings& o);
friend void from_json(const nlohmann::json& j, Spot_finding_settings& o);
@@ -117,6 +122,8 @@ protected:
float m_Low_resolution_limit;
float m_Indexing_tolerance;
};

View File

@@ -167,6 +167,16 @@ components:
minimum: 0.015625
maximum: 4.0
description: For JUNGFRAU conversion it is possible to multiply energy by a given factor to get fractional/multiplied photon counts
data_reduction_factor_serialmx:
type: number
format: float
default: 1.0
minimum: 0.0
maximum: 1.0
description: |
Rate at which non-indexed images are accepted to be forwarded to writer.
Value of 1.0 (default) means that all images are written.
Values below zero mean that non-indexed images will be accepted with a given probability.
unit_cell:
type: object
description: Units of angstrom and degree
@@ -322,6 +332,7 @@ components:
- min_pix_per_spot
- high_resolution_limit
- low_resolution_limit
- indexing_tolerance
properties:
enable:
type: boolean
@@ -353,6 +364,12 @@ components:
low_resolution_limit:
type: number
format: float
indexing_tolerance:
type: number
format: float
minimum: 0.0
maximum: 1.0
description: Acceptance tolerance for spots after the indexing run - the larger the number, the more spots will be accepted
rad_int_settings:
type: object
required:
@@ -478,6 +495,12 @@ components:
bkg_estimate:
type: number
format: float
beam_center_drift_x_pxl:
type: number
format: float
beam_center_drift_y_pxl:
type: number
format: float
broker_status:
type: object
required:
@@ -1276,6 +1299,30 @@ paths:
schema:
type: string
description: Exception error
/plot/beam_center_drift:
post:
summary: Generate plot of beam center drift (from indexing)
description: Max count of beam center drift per image; binning is configurable
requestBody:
content:
application/json:
schema:
$ref: '#/components/schemas/plot_request'
responses:
"200":
description: Everything OK
content:
application/json:
schema:
$ref: '#/components/schemas/plots'
"400":
description: Input parsing or validation error
content:
text/plain:
schema:
type: string
description:
Exception error
/plot/roi_max_count:
post:
summary: Generate plot of ROI max count

File diff suppressed because one or more lines are too long

View File

@@ -27,6 +27,7 @@ DatasetSettings::DatasetSettings() {
photon_energy_multiplier = 1.0f;
omega_start = 0.0f;
images_per_file = 1000;
data_reduction_factor_serialmx = 1.0;
}
DatasetSettings &DatasetSettings::ImagesPerTrigger(int64_t input) {
@@ -329,4 +330,15 @@ DatasetSettings &DatasetSettings::ImagesPerFile(int64_t input) {
int64_t DatasetSettings::GetImagesPerFile() const {
return images_per_file;
}
}
DatasetSettings &DatasetSettings::DataReductionFactorSerialMX(float input) {
check_min("Data reduction factor for serial MX", input, 0.0);
check_max("Data reduction factor for serial MX", input, 1.0);
data_reduction_factor_serialmx = input;
return *this;
}
float DatasetSettings::GetDataReductionFactorSerialMX() const {
return data_reduction_factor_serialmx;
}

View File

@@ -47,6 +47,8 @@ class DatasetSettings {
std::string header_appendix;
std::optional<Coord> omega_rotation_axis;
float data_reduction_factor_serialmx;
public:
DatasetSettings();
@@ -74,6 +76,7 @@ public:
DatasetSettings& Summation(int64_t input);
DatasetSettings& FPGAOutputMode(FPGAPixelOutput input);
DatasetSettings& ImagesPerFile(int64_t input);
DatasetSettings& DataReductionFactorSerialMX(float input);
std::optional<float> GetAttenuatorTransmission() const;
std::optional<float> GetTotalFlux() const;
@@ -103,6 +106,7 @@ public:
int64_t GetNumTriggers() const;
int64_t GetImageNumPerTrigger() const;
int64_t GetImagesPerFile() const;
float GetDataReductionFactorSerialMX() const;
};
#endif //JUNGFRAUJOCH_DATASETSETTINGS_H

View File

@@ -29,6 +29,6 @@
#define DEFAULT_G2_FACTOR (-0.1145)
#define DEFAULT_HG0_FACTOR (100.0)
#define MAX_SPOT_COUNT (100)
#define MAX_SPOT_COUNT (250)
#endif //DEFINITIONS_H

View File

@@ -66,7 +66,7 @@ DiffractionExperiment::DiffractionExperiment(const DetectorSetup& det_setup)
mode = DetectorMode::Conversion;
max_spot_count = MAX_SPOT_COUNT;
max_spot_count = 100;
}
// setter functions
@@ -640,7 +640,7 @@ float DiffractionExperiment::GetQSpacingForAzimInt_recipA() const {
DiffractionExperiment &DiffractionExperiment::MaxSpotCount(int64_t input) {
check_min("Max spot count", input, 10);
check_max("Max spot count", input, 500);
check_max("Max spot count", input, MAX_SPOT_COUNT);
max_spot_count = input;
return *this;
}
@@ -669,6 +669,8 @@ void DiffractionExperiment::CheckDataProcessingSettings(const SpotFindingSetting
check_min("Spot finding low resolution limit", settings.low_resolution_limit, 1.0);
check_max("Spot finding low resolution limit", settings.low_resolution_limit, 50.0);
}
check_min("indexing tolerance", settings.indexing_tolerance, 0.0);
check_max("indexing tolerance", settings.indexing_tolerance, 1.0);
}
SpotFindingSettings DiffractionExperiment::DefaultDataProcessingSettings() {
@@ -1174,3 +1176,6 @@ int64_t DiffractionExperiment::GetSendBufferLocationSize() const {
return GetMaxCompressedSize() + 1024 * 1024;
}
float DiffractionExperiment::GetDataReductionFactorSerialMX() const {
return dataset.GetDataReductionFactorSerialMX();
}

View File

@@ -324,6 +324,8 @@ public:
const ROIMask& ROI() const;
void ExportROIMask(uint16_t *v, size_t module_number) const;
int64_t GetImagesPerFile() const;
float GetDataReductionFactorSerialMX() const;
};
#endif //DIFFRACTIONEXPERIMENT_H

View File

@@ -9,7 +9,7 @@
enum class PlotType {BkgEstimate, RadInt, RadIntPerFile, SpotCount, IndexingRate, IndexingRatePerFile,
ErrorPixels, ImageCollectionEfficiency, ReceiverDelay, ReceiverFreeSendBuf, StrongPixels,
ROISum, ROIMaxCount, ROIPixels,
ROISum, ROIMaxCount, ROIPixels, IndexingDrift,
ResEstimation};
struct PlotRequest {

View File

@@ -5,6 +5,11 @@
ZMQContext::ZMQContext() {
context = zmq_ctx_new();
// Default is to have 2 I/O threads per ZMQ context
if (zmq_ctx_set(context, ZMQ_IO_THREADS, 2) != 0)
throw JFJochException(JFJochExceptionCategory::ZeroMQ,
"Cannot set number of I/O threads");
}
ZMQContext &ZMQContext::NumThreads(int32_t threads) {
@@ -99,19 +104,6 @@ void ZMQSocket::SendZeroCopy(const void *buf, size_t buf_size, ZeroCopyReturnVal
}
}
int64_t ZMQSocket::Receive(bool blocking) {
std::vector<uint8_t> msg;
return Receive(msg, blocking, true);
}
int64_t ZMQSocket::Receive(std::string &s, bool blocking) {
std::vector<uint8_t> v;
int64_t rc = Receive(v, blocking, true);
if (rc > 0)
s = std::string(v.begin(), v.end());
return rc;
}
void ZMQSocket::SetSocketOption(int32_t option_name, int32_t value) {
if (zmq_setsockopt(socket, option_name, &value, sizeof(value)) != 0)
throw JFJochException(JFJochExceptionCategory::ZeroMQ, "Cannot set socket option");
@@ -200,3 +192,30 @@ std::string ZMQSocket::GetEndpointName() {
zmq_getsockopt(socket, ZMQ_LAST_ENDPOINT, tmp, &len);
return tmp;
}
bool ZMQSocket::Receive(ZMQMessage &msg, bool blocking) {
return (zmq_msg_recv(msg.GetMsg(), socket, blocking ? 0 : ZMQ_DONTWAIT) >= 0);
}
ZMQMessage::ZMQMessage() {
int rc = zmq_msg_init(&msg);
if (rc) throw JFJochException(JFJochExceptionCategory::ZeroMQ, "Cannot init zmq_msg");
}
ZMQMessage::~ZMQMessage() {
zmq_msg_close(&msg);
}
zmq_msg_t *ZMQMessage::GetMsg() {
return &msg;
}
size_t ZMQMessage::size() {
return zmq_msg_size(&msg);
}
const uint8_t *ZMQMessage::data() {
return (uint8_t *) zmq_msg_data(&msg);
}

View File

@@ -28,6 +28,18 @@ public:
enum class ZMQSocketType : int {Push = ZMQ_PUSH, Pull = ZMQ_PULL, Req = ZMQ_REQ, Rep = ZMQ_REP,
Pub = ZMQ_PUB, Sub = ZMQ_SUB};
class ZMQMessage {
zmq_msg_t msg;
public:
ZMQMessage();
~ZMQMessage();
ZMQMessage(ZMQMessage &other) = delete;
ZMQMessage& operator=(ZMQMessage &other) = delete;
zmq_msg_t *GetMsg();
size_t size();
const uint8_t *data();
};
class ZMQSocket {
std::mutex m;
ZMQSocketType socket_type;
@@ -54,33 +66,7 @@ public:
ZMQSocket &SendBufferSize(int32_t bytes);
ZMQSocket &ReceiverBufferSize(int32_t bytes);
int64_t Receive(bool blocking = true);
int64_t Receive(std::string &j, bool blocking = true);
template <class T> int64_t Receive(std::vector<T> &vector, bool blocking = true, bool resize = true) {
std::unique_lock<std::mutex> ul(m);
zmq_msg_t zmq_msg;
zmq_msg_init(&zmq_msg);
int64_t msg_size = zmq_msg_recv(&zmq_msg, socket, blocking ? 0 : ZMQ_DONTWAIT);
if (msg_size < 0) {
if ((errno == EAGAIN) || (errno == EINTR))
return -1;
else
throw JFJochException(JFJochExceptionCategory::ZeroMQ, "zmq_msg_recv failed "
+ std::string(strerror(errno)));
} else if (msg_size == 0) {
zmq_msg_close (&zmq_msg);
return 0;
} else if (resize) {
vector.resize(msg_size / sizeof(T) + ((msg_size % sizeof(T) != 0) ? 1 : 0));
} else {
zmq_msg_close (&zmq_msg);
throw JFJochException(JFJochExceptionCategory::ZeroMQ, "Receive buffer too small");
}
memcpy(vector.data(), zmq_msg_data(&zmq_msg), msg_size);
zmq_msg_close (&zmq_msg);
return msg_size;
}
bool Receive(ZMQMessage& msg, bool blocking = true);
void Send();
bool Send(const void *buf, size_t buf_size, bool blocking = true, bool multipart = false);

View File

@@ -11,7 +11,7 @@
std::string WriteJPEGToMem(const std::vector<uint8_t> &input, size_t width, size_t height, int quality) {
unsigned char *buf;
unsigned char *buf = nullptr;
unsigned long buf_size = 0;
if (width * height * 3 != input.size())
@@ -45,8 +45,10 @@ std::string WriteJPEGToMem(const std::vector<uint8_t> &input, size_t width, size
jpeg_finish_compress(&cinfo);
jpeg_destroy_compress(&cinfo);
std::string ret;
if (buf_size > 0)
return {(char *)buf, buf_size};
else
return "";
ret = std::string((char *)buf, buf_size);
free(buf);
return ret;
}

View File

@@ -48,6 +48,9 @@ struct DataMessage {
float bkg_estimate;
uint64_t indexing_result; // 0 - not tried, 1 - tried and failed, 2 - tried and success
std::vector<float> indexing_lattice;
std::optional<float> indexing_drift_x_pxl;
std::optional<float> indexing_drift_y_pxl;
UnitCell indexing_unit_cell;
std::vector<uint64_t> adu_histogram;

View File

@@ -9,7 +9,7 @@ type MyProps = {};
class BkgEstimatePlot extends Component<MyProps> {
render() {
return <Paper style={{textAlign: 'center'}} sx={{height: 400, width: "100%"}}>
return <Paper style={{textAlign: 'center'}} sx={{height: 450, width: "100%"}}>
<Box sx={{width: "100%", height: 50}}>
<br/>
<center><strong>Background estimate</strong></center>

View File

@@ -17,7 +17,8 @@ export enum PlotType {
ROI_SUM,
ROI_MAX_COUNT,
RES_ESTIMATION,
RECEIVER_FREE_SEND_BUFS
RECEIVER_FREE_SEND_BUFS,
BEAM_CENTER_DRIFT
}
type MyProps = {
@@ -157,8 +158,17 @@ class DataProcessingPlot extends Component<MyProps, MyState> {
.catch(error => {
this.setState({connection_error: true});
});
break;
case PlotType.BEAM_CENTER_DRIFT:
DefaultService.postPlotBeamCenterDrift({binning: this.props.binning})
.then(data => this.setState({plots: data, connection_error: false}))
.catch(error => {
this.setState({connection_error: true});
});
break;
}
}
componentDidMount() {
this.getValues();
this.interval = setInterval(() => this.getValues(), 500);

View File

@@ -75,6 +75,9 @@ class DataProcessingPlots extends Component<MyProps, MyState> {
case "13":
this.setState({type: PlotType.RECEIVER_FREE_SEND_BUFS, xlabel: "Image number", ylabel: "Number of buffers"});
break;
case "14":
this.setState({type: PlotType.BEAM_CENTER_DRIFT, xlabel: "Image number", ylabel: "Drift [pxl]"});
break;
}
};
@@ -102,6 +105,7 @@ class DataProcessingPlots extends Component<MyProps, MyState> {
<MenuItem value={6}>ROI area max count</MenuItem>
<MenuItem value={12}>Crystal resolution histogram</MenuItem>
<MenuItem value={4}>Indexing rate (per file)</MenuItem>
<MenuItem value={14}>Beam center drift (based on indexing)</MenuItem>
<MenuItem value={5}>Azimuthal integration profile (per file)</MenuItem>
<MenuItem value={8}>Error and saturated pixels</MenuItem>
<MenuItem value={10}>Strong pixels (for spot finding)</MenuItem>

View File

@@ -22,7 +22,8 @@ class DataProcessingSettings extends Component<MyProps, MyState> {
min_pix_per_spot: 2,
max_pix_per_spot: 50,
high_resolution_limit: 2.5,
low_resolution_limit: 20.0
low_resolution_limit: 20.0,
indexing_tolerance: 0.05
},
connection_error: true
}
@@ -91,6 +92,19 @@ class DataProcessingSettings extends Component<MyProps, MyState> {
this.getValues();
}
setIndexingTolerance = (event: Event, newValue: number | number[]) => {
this.setState(prevState => (
{
s : {
...prevState.s,
indexing_tolerance: newValue as number
}
}
));
this.putValues();
this.getValues();
}
setHighResolutionLimit = (event: Event, newValue: number | number[]) => {
this.setState(prevState => (
{
@@ -154,6 +168,11 @@ class DataProcessingSettings extends Component<MyProps, MyState> {
value={Number(this.state.s.high_resolution_limit)}
onChange={this.setHighResolutionLimit}
min={1} max={5} step={0.2} valueLabelDisplay="auto" />
<Typography> Indexing spot acceptance tolerance </Typography>
<Slider disabled={this.state.connection_error || !this.state.s.enable}
value={Number(this.state.s.indexing_tolerance)}
onChange={this.setIndexingTolerance}
min={0.0} max={0.3} step={0.01} valueLabelDisplay="auto" />
<br/><br/>
</Grid>
<Grid item xs={1}/>

View File

@@ -47,7 +47,7 @@ class MeasurementStatistics extends Component<MyProps, MyState> {
}
render() {
return <Paper style={{textAlign: 'center'}} sx={{ height: 400, width: '100%' }}>
return <Paper style={{textAlign: 'center'}} sx={{ height: 450, width: '100%' }}>
<Grid container spacing={0}>
<Grid item xs={1}/>
<Grid item xs={10}>
@@ -75,23 +75,30 @@ class MeasurementStatistics extends Component<MyProps, MyState> {
<TableRow>
<TableCell component="th" scope="row"> Compression ratio: </TableCell>
<TableCell align="right">{(this.state.s.compression_ratio !== undefined)
? (this.state.s.compression_ratio.toFixed(1)) + "x" : "" } </TableCell>
? (this.state.s.compression_ratio.toFixed(1)) + "x" : "-" } </TableCell>
</TableRow>
<TableRow>
<TableCell component="th" scope="row"> Data acquisition efficiency: </TableCell>
<TableCell align="right">{(this.state.s.collection_efficiency !== undefined)
? (Math.floor(this.state.s.collection_efficiency * 1000.0) / 1000.0).toFixed(3) : ""}</TableCell>
? (Math.floor(this.state.s.collection_efficiency * 1000.0) / 1000.0).toFixed(3) : "-"}</TableCell>
</TableRow>
<TableRow>
<TableCell component="th" scope="row"> Indexing rate: </TableCell>
<TableCell align="right">{(this.state.s.indexing_rate !== undefined)
? this.state.s.indexing_rate.toFixed(2) : ""}</TableCell>
? this.state.s.indexing_rate.toFixed(2) : "-"}</TableCell>
</TableRow>
<TableRow>
<TableCell component="th" scope="row"> Beam center drift (x/y) [pxl]: </TableCell>
<TableCell align="right">{(this.state.s.beam_center_drift_x_pxl !== undefined)
? this.state.s.beam_center_drift_x_pxl.toFixed(2) : "-"}&nbsp;/&nbsp;
{(this.state.s.beam_center_drift_y_pxl !== undefined)
? this.state.s.beam_center_drift_y_pxl.toFixed(2) : "-"}</TableCell>
</TableRow>
<TableRow>
<TableCell component="th" scope="row"> Background estimate: </TableCell>
<TableCell align="right">{(this.state.s.bkg_estimate !== undefined)
? this.state.s.bkg_estimate.toPrecision(7) : ""}</TableCell>
? this.state.s.bkg_estimate.toPrecision(7) : "-"}</TableCell>
</TableRow>
</TableBody>
</Table>

View File

@@ -97,6 +97,13 @@ export type dataset_settings = {
* For JUNGFRAU conversion it is possible to multiply energy by a given factor to get fractional/multiplied photon counts
*/
photon_energy_multiplier?: number;
/**
* Rate at which non-indexed images are accepted to be forwarded to writer.
* Value of 1.0 (default) means that all images are written.
* Values below zero mean that non-indexed images will be accepted with a given probability.
*
*/
data_reduction_factor_serialmx?: number;
/**
* Units of angstrom and degree
*/

View File

@@ -18,6 +18,8 @@ export type measurement_statistics = {
detector_height?: number;
detector_pixel_depth?: measurement_statistics.detector_pixel_depth;
bkg_estimate?: number;
beam_center_drift_x_pxl?: number;
beam_center_drift_y_pxl?: number;
};
export namespace measurement_statistics {

View File

@@ -18,5 +18,9 @@ export type spot_finding_settings = {
max_pix_per_spot: number;
high_resolution_limit: number;
low_resolution_limit: number;
/**
* Acceptance tolerance for spots after the indexing run - the larger the number, the more spots will be accepted
*/
indexing_tolerance: number;
};

View File

@@ -596,6 +596,27 @@ export class DefaultService {
});
}
/**
* Generate plot of beam center drift (from indexing)
* Max count of beam center drift per image; binning is configurable
* @param requestBody
* @returns plots Everything OK
* @throws ApiError
*/
public static postPlotBeamCenterDrift(
requestBody?: plot_request,
): CancelablePromise<plots> {
return __request(OpenAPI, {
method: 'POST',
url: '/plot/beam_center_drift',
body: requestBody,
mediaType: 'application/json',
errors: {
400: `Input parsing or validation error`,
},
});
}
/**
* Generate plot of ROI max count
* Max count of ROI per image; binning is configurable

View File

@@ -28,6 +28,13 @@ IF (CMAKE_CUDA_COMPILER)
fast-feedback-indexer/indexer/src/ffbidx/log.h
fast-feedback-indexer/indexer/src/ffbidx/exception.h)
ADD_CUSTOM_TARGET(version_txt
${CMAKE_COMMAND} -D SRC=${CMAKE_CURRENT_SOURCE_DIR}/fast-feedback-indexer/indexer/src/version.h.in
-D DST=${CMAKE_CURRENT_SOURCE_DIR}/fast-feedback-indexer/indexer/src/ffbidx/version.h
-P ${CMAKE_CURRENT_SOURCE_DIR}/fast-feedback-indexer/indexer/src/GenerateVersionH.cmake
)
ADD_DEPENDENCIES(JFJochImageAnalysis version_txt)
TARGET_COMPILE_DEFINITIONS(JFJochImageAnalysis PUBLIC -DVECTOR_CANDIDATE_REFINEMENT=1)
TARGET_INCLUDE_DIRECTORIES(JFJochImageAnalysis PUBLIC fast-feedback-indexer/indexer/src/)
TARGET_LINK_LIBRARIES(JFJochImageAnalysis ${CUDART_LIBRARY} ${CMAKE_DL_LIBS} rt)

View File

@@ -1,6 +1,7 @@
// Copyright (2019-2023) Paul Scherrer Institute
// Copyright (2019-2024) Paul Scherrer Institute
#include "IndexerWrapper.h"
#define MAX_SPOT_COUNT_INDEXING 50
void IndexerWrapper::Setup(const UnitCell &cell) {
#ifdef JFJOCH_USE_CUDA
@@ -8,16 +9,18 @@ void IndexerWrapper::Setup(const UnitCell &cell) {
#endif
}
std::vector<IndexingResult> IndexerWrapper::Run(const std::vector<Coord> &coord) {
std::vector<IndexingResult> IndexerWrapper::Run(const std::vector<Coord> &coord, float indexing_threshold) {
#ifdef JFJOCH_USE_CUDA
std::vector<IndexingResult> ret;
if (coord.size() <= min_spots)
return ret;
size_t nspots = std::min<size_t>(MAX_SPOT_COUNT, coord.size());
assert(coord.size() < MAX_SPOT_COUNT);
for (int i = 0; i < nspots; i++) {
size_t nspots = std::min<size_t>(MAX_SPOT_COUNT_INDEXING, coord.size());
for (int i = 0; i < coord.size(); i++) {
indexer.spotX(i) = coord[i].x;
indexer.spotY(i) = coord[i].y;
indexer.spotZ(i) = coord[i].z;
@@ -26,28 +29,40 @@ std::vector<IndexingResult> IndexerWrapper::Run(const std::vector<Coord> &coord)
// Index
indexer.index(1, nspots);
fast_feedback::refine::indexer_ifssr<float>::refine(indexer.spotM().topRows(coord.size()),
indexer.oCellM(),
indexer.oScoreV(),
fast_feedback::refine::config_ifssr<float>{});
// Find cell most similar to the current one
for (unsigned i=0u; i<cpers.max_output_cells; i++)
indexer.oScore(i) += fast_feedback::refine::cell_similarity(indexer.oCell(i), indexer.iCell(0), 0.02f);
// Get best cell
auto id = fast_feedback::refine::best_cell(indexer.oScoreV());
// get indexed spots
using M3x = Eigen::MatrixX3<float>;
M3x resid = indexer.Spots() * indexer.oCell(id).transpose();
const M3x miller = round(resid.array());
resid -= miller;
auto arr = resid.rowwise().norm().array() < threshold;
auto indexed_spot_count = arr.count();
bool indexed = fast_feedback::refine::is_viable_cell(indexer.oCell(id), indexer.Spots(), indexing_threshold, 9);
// Check if result is viable
if (indexed_spot_count > min_spots) {
if (indexed) {
auto cell = indexer.oCell(id).colwise().reverse();
fast_feedback::refine::make_right_handed(cell);
// get indexed spots
using M3x = Eigen::MatrixX3<float>;
M3x resid = indexer.spotM().topRows(coord.size()) * cell.transpose();
const M3x miller = round(resid.array());
const M3x predicted = miller * cell.transpose().inverse();
IndexingResult result;
result.l = CrystalLattice(indexer.oCell(id));
result.l = CrystalLattice(cell);
result.indexed_spots.resize(coord.size());
for (int i = 0; i < arr.size(); i++)
result.indexed_spots[i] = arr[i];
result.predicted_spots.resize(coord.size());
result.indexed_spots_count = indexed_spot_count;
for (int i = 0; i < coord.size(); i++)
result.predicted_spots[i] = Coord(predicted.coeff(i, 0),
predicted.coeff(i, 1),
predicted.coeff(i, 2));
ret.emplace_back(result);
}

View File

@@ -1,4 +1,4 @@
// Copyright (2019-2023) Paul Scherrer Institute
// Copyright (2019-2024) Paul Scherrer Institute
#ifndef JUNGFRAUJOCH_INDEXERWRAPPER_H
#define JUNGFRAUJOCH_INDEXERWRAPPER_H
@@ -15,26 +15,24 @@
struct IndexingResult {
CrystalLattice l;
std::vector<uint8_t> indexed_spots;
uint64_t indexed_spots_count;
std::vector<Coord> predicted_spots;
};
class IndexerWrapper {
#ifdef JFJOCH_USE_CUDA
fast_feedback::config_runtime<float> crt{};
fast_feedback::config_persistent<float> cpers{
.max_output_cells = 4,
.max_output_cells = 16,
.max_input_cells = 1,
.max_spots = MAX_SPOT_COUNT
}; // default persistent config
fast_feedback::refine::config_ifss<float> conf_ifss{};
fast_feedback::refine::indexer_ifss<float> indexer{cpers, crt, conf_ifss};
fast_feedback::refine::indexer<float> indexer{cpers, crt};
#endif
constexpr const static uint64_t min_spots = 9;
constexpr const static float threshold = 0.05f;
public:
void Setup(const UnitCell &cell);
std::vector<IndexingResult> Run(const std::vector<Coord> &coord);
std::vector<IndexingResult> Run(const std::vector<Coord> &coord, float indexing_threshold);
};

View File

@@ -4,6 +4,24 @@
#include "CPUSpotFinder.h"
#include "../common/DiffractionGeometry.h"
double stddev(const std::vector<float> &v) {
if (v.size() <= 1)
return 0.0;
double mean = 0.0f;
for (const auto &i: v)
mean += i;
mean /= v.size();
double stddev = 0.0f;
for (const auto &i: v)
stddev += (i - mean) * (i - mean);
return sqrt(stddev / (v.size() - 1));
}
MXAnalyzer::MXAnalyzer(const DiffractionExperiment &in_experiment)
: experiment(in_experiment) {
auto uc = experiment.GetUnitCell();
@@ -39,8 +57,7 @@ void MXAnalyzer::ReadFromCPU(const int16_t *image, const SpotFindingSettings &se
ReadFromFPGA(&output, settings, module_number);
}
bool MXAnalyzer::Process(DataMessage &message,
const SpotFindingSettings& settings) {
bool MXAnalyzer::Process(DataMessage &message, const SpotFindingSettings& settings) {
message.indexing_result = 0;
if (!find_spots)
return false;
@@ -55,18 +72,33 @@ bool MXAnalyzer::Process(DataMessage &message,
if (do_indexing && settings.indexing) {
std::vector<Coord> recip;
recip.reserve(spots_out.size());
for (const auto &i: spots_out)
recip.push_back(i.ReciprocalCoord(experiment));
auto indexer_result = indexer.Run(recip);
auto indexer_result = indexer.Run(recip, settings.indexing_tolerance);
float x_drift = 0.0f, y_drift = 0.0f;
if (!indexer_result.empty()) {
message.indexing_result = 1;
for (int i = 0; i < recip.size(); i++)
message.spots[i].indexed = indexer_result[0].indexed_spots[i];
for (int i = 0; i < recip.size(); i++) {
auto predicted_pos = RecipToDector(experiment, indexer_result[0].predicted_spots[i]);
float x_diff = predicted_pos.first - spots_out[i].RawCoord().x;
float y_diff = predicted_pos.second - spots_out[i].RawCoord().y;
x_drift += x_diff;
y_drift += y_diff;
message.spots[i].indexed = (x_diff * x_diff + y_diff * y_diff < spot_distance_threshold);
}
indexer_result[0].l.Save(message.indexing_lattice);
message.indexing_unit_cell = indexer_result[0].l.GetUnitCell();
indexed = true;
if (!recip.empty()) {
message.indexing_drift_x_pxl = x_drift / recip.size();
message.indexing_drift_y_pxl = y_drift / recip.size();
}
}
}
spots.clear();

View File

@@ -13,6 +13,7 @@ class MXAnalyzer {
bool do_indexing = false;
bool find_spots = false;
std::vector<DiffractionSpot> spots;
constexpr static const float spot_distance_threshold = 2.0f * 2.0f;
public:
explicit MXAnalyzer(const DiffractionExperiment& experiment);

View File

@@ -14,6 +14,8 @@ struct SpotFindingSettings {
int64_t max_pix_per_spot = 50; // Maximum pixels per spot
float high_resolution_limit = 2.5;
float low_resolution_limit = 50.0;
float indexing_tolerance = 0.05;
};
#endif //JUNGFRAUJOCH_SPOTFINDINGSETTINGS_H

View File

@@ -21,13 +21,10 @@ struct res_id {
void FilterSpotsByCount(const DiffractionExperiment& experiment,
const std::vector<DiffractionSpot> &input,
std::vector<DiffractionSpot> &output) {
if (input.size() < experiment.GetMaxSpotCount())
output = input;
size_t output_size = std::min<size_t>(input.size(), experiment.GetMaxSpotCount());
std::vector<intensity_id> intensity_id_vector(input.size());
for (int i = 0; i < input.size(); i++) {
intensity_id_vector[i].I = input[i].Count();
intensity_id_vector[i].id = i;

View File

@@ -36,7 +36,8 @@ JFJochReceiver::JFJochReceiver(const DiffractionExperiment& in_experiment,
preview_counter_indexed(experiment.GetPreviewPeriod()),
spot_finding_settings(in_spot_finding_settings),
preview_image(experiment),
preview_image_indexed(experiment)
preview_image_indexed(experiment),
serialmx_filter(experiment)
{
if (experiment.GetDetectorSetup().GetDetectorType() == DetectorType::JUNGFRAU)
calibration = in_calibration;
@@ -361,7 +362,7 @@ void JFJochReceiver::FrameTransformationThread() {
if (indexed && preview_counter_indexed.GeneratePreview())
preview_image_indexed.UpdateImage(transformation.GetImage(), message.spots);
if (push_images_to_writer) {
if (push_images_to_writer && serialmx_filter.ApplyFilter(message)) {
auto loc = send_buf_ctrl.GetBufLocation();
if (loc != nullptr) {
auto writer_buffer = (uint8_t *)loc->get_ptr();
@@ -555,35 +556,29 @@ JFJochReceiverOutput JFJochReceiver::GetStatistics() const {
}
JFJochReceiverStatus JFJochReceiver::GetStatus() const {
float indexing_rate = -1, bkg_estimate = -1, compressed_ratio = 0.0;
auto tmp = plots.GetIndexingRate();
if (!std::isnan(tmp))
indexing_rate = tmp;
JFJochReceiverStatus ret;
tmp = plots.GetBkgEstimate();
if (!std::isnan(tmp))
bkg_estimate = tmp;
ret.indexing_rate = plots.GetIndexingRate();
ret.bkg_estimate = plots.GetBkgEstimate();
ret.beam_center_drift_pxl = plots.GetBeamCenterDriftPxl();
if ((experiment.GetImageNum() > 0) && (compressed_size > 0)) {
compressed_ratio = static_cast<double> (images_sent
ret.compressed_ratio = static_cast<double> (images_sent
* experiment.GetPixelDepth()
* experiment.GetModulesNum()
* RAW_MODULE_SIZE)
/ static_cast<double> (compressed_size);
}
return {
.progress = GetProgress(),
.indexing_rate = indexing_rate,
.bkg_estimate = bkg_estimate,
.compressed_ratio = compressed_ratio,
.compressed_size = compressed_size,
.max_receive_delay = max_delay,
.max_image_number_sent = max_image_number_sent,
.images_collected = images_collected,
.images_sent = images_sent,
.cancelled = cancelled
};
ret.progress = GetProgress();
ret.compressed_size = compressed_size;
ret.max_receive_delay = max_delay;
ret.max_image_number_sent = max_image_number_sent;
ret.images_collected = images_collected;
ret.images_sent = images_sent;
ret.cancelled = cancelled;
return ret;
}
std::string JFJochReceiver::GetTIFF(bool calibration_run) const {

View File

@@ -31,12 +31,14 @@
#include "JFJochReceiverPlots.h"
#include "PreviewCounter.h"
#include "../export_images/PreviewImage.h"
#include "LossyFilter.h"
struct JFJochReceiverStatus {
float progress;
float indexing_rate;
float bkg_estimate;
float compressed_ratio;
std::optional<float> indexing_rate;
std::optional<float> bkg_estimate;
std::optional<float> compressed_ratio;
std::optional<std::pair<float, float>> beam_center_drift_pxl;
uint64_t compressed_size;
uint64_t max_receive_delay;
uint64_t max_image_number_sent;
@@ -113,6 +115,8 @@ class JFJochReceiver {
PreviewCounter preview_counter_indexed;
PreviewImage preview_image_indexed;
LossyFilter serialmx_filter;
void AcquireThread(uint16_t data_stream);
void FrameTransformationThread();
void Cancel(const JFJochException &e);

View File

@@ -36,6 +36,11 @@ void JFJochReceiverPlots::Add(const DataMessage &msg, uint64_t file_number, cons
indexing_solution.AddElement(msg.number, msg.indexing_result);
indexing_solution_per_file.Add(file_number, msg.indexing_result);
if (msg.indexing_drift_x_pxl)
indexing_drift_x_pxl.AddElement(msg.number, msg.indexing_drift_x_pxl.value());
if (msg.indexing_drift_y_pxl)
indexing_drift_y_pxl.AddElement(msg.number, msg.indexing_drift_y_pxl.value());
if (msg.resolution_estimation)
resolution_estimation.Add(msg.resolution_estimation.value());
@@ -79,6 +84,14 @@ MultiLinePlot JFJochReceiverPlots::GetPlots(const PlotRequest &request) {
return bkg_estimate.GetMeanPlot(nbins);
case PlotType::IndexingRatePerFile:
return indexing_solution_per_file.GetPlot();
case PlotType::IndexingDrift:
tmp = indexing_drift_x_pxl.GetMeanPlot(nbins);
tmp.at(0).title = "X beam drift (pxl)";
ret.emplace_back(tmp.at(0));
tmp = indexing_drift_y_pxl.GetMeanPlot(nbins);
tmp.at(0).title = "Y beam drift (pxl)";
ret.emplace_back(tmp.at(0));
return ret;
case PlotType::ErrorPixels:
tmp = saturated_pixels.GetMeanPlot(nbins);
tmp.at(0).title = "Saturated pixels (mean)";
@@ -138,12 +151,29 @@ MultiLinePlot JFJochReceiverPlots::GetPlots(const PlotRequest &request) {
}
}
float JFJochReceiverPlots::GetIndexingRate() const {
return indexing_solution.Mean();
std::optional<float> JFJochReceiverPlots::GetIndexingRate() const {
auto tmp = indexing_solution.Mean();
if (std::isfinite(tmp))
return tmp;
else
return {};
}
float JFJochReceiverPlots::GetBkgEstimate() const {
return bkg_estimate.Mean();
std::optional<float> JFJochReceiverPlots::GetBkgEstimate() const {
auto tmp = bkg_estimate.Mean();
if (std::isfinite(tmp))
return tmp;
else
return {};
}
std::optional<std::pair<float, float>> JFJochReceiverPlots::GetBeamCenterDriftPxl() const {
auto diffx = indexing_drift_x_pxl.Mean();
auto diffy = indexing_drift_y_pxl.Mean();
if (std::isfinite(diffx) && std::isfinite(diffy))
return std::make_pair(diffx, diffy);
else
return {};
}
void JFJochReceiverPlots::GetXFELPulseID(std::vector<uint64_t> &v) const {
@@ -164,3 +194,5 @@ std::vector<float> JFJochReceiverPlots::GetAzIntProfilePerFile(uint64_t file_num
else
throw JFJochException(JFJochExceptionCategory::ArrayOutOfBounds, "JFJochReceiverPlots::GetAzIntProfilePerFile out of bounds");
}

View File

@@ -27,6 +27,9 @@ class JFJochReceiverPlots {
StatusVector<float> image_collection_efficiency;
StatusVector<float> unit_cell[6];
SetAverage<uint64_t> indexing_solution_per_file;
StatusVector<float> indexing_drift_x_pxl;
StatusVector<float> indexing_drift_y_pxl;
std::map<std::string, std::unique_ptr<StatusVector<int64_t>>> roi_sum;
std::map<std::string, std::unique_ptr<StatusVector<int64_t>>> roi_max_count;
std::map<std::string, std::unique_ptr<StatusVector<uint64_t>>> roi_pixels;
@@ -43,10 +46,13 @@ public:
void GetXFELPulseID(std::vector<uint64_t>& v) const;
void GetXFELEventCode(std::vector<uint64_t>& v) const;
float GetIndexingRate() const;
float GetBkgEstimate() const;
std::optional<float> GetIndexingRate() const;
std::optional<float> GetBkgEstimate() const;
std::optional<std::pair<float, float>> GetBeamCenterDriftPxl() const;
std::vector<float> GetAzIntProfile() const;
std::vector<float> GetAzIntProfilePerFile(uint64_t file_number);
};

View File

@@ -2,20 +2,26 @@
#include "LossyFilter.h"
LossyFilter::LossyFilter(bool in_enabled, float in_p)
: enabled(in_enabled), p(in_p) {}
LossyFilter::LossyFilter(float in_p)
: p(in_p) {}
bool LossyFilter::RollDice() {
std::unique_lock<std::mutex> ul(random_m);
if ((p > 0.0) && (distr(mt) < p))
if (distr(mt) < p)
return true;
else
return false;
}
bool LossyFilter::ApplyFilter(DataMessage &message) {
if (!enabled)
if ((p == 1.0) || message.indexing_result)
return true;
return (message.indexing_result || RollDice());
else if (p == 0.0)
return message.indexing_result;
return RollDice();
}
LossyFilter::LossyFilter(const DiffractionExperiment &x)
: p (x.GetDataReductionFactorSerialMX()){}

View File

@@ -4,22 +4,24 @@
#define JUNGFRAUJOCH_LOSSYFILTER_H
#include <random>
#include <atomic>
#include <mutex>
#include <cstdint>
#include "../frame_serialize/JFJochMessages.h"
#include "../common/DiffractionExperiment.h"
class LossyFilter {
std::mutex random_m;
std::random_device rdev;
std::mt19937 mt{rdev()};
std::uniform_real_distribution<float> distr{0.0, 1.0};
mutable std::mutex random_m;
bool enabled;
float p;
bool RollDice();
public:
LossyFilter(bool enabled, float p);
explicit LossyFilter(float p);
explicit LossyFilter(const DiffractionExperiment& x);
bool ApplyFilter(DataMessage& message);
};

View File

@@ -6,6 +6,10 @@ PreviewCounter::PreviewCounter(std::chrono::microseconds in_period) : period(in_
bool PreviewCounter::GeneratePreview() {
std::unique_lock<std::mutex> ul(m);
if (period.count() == 0)
return false;
auto now = std::chrono::system_clock::now();
if (now > last_preview + period) {

View File

@@ -194,7 +194,8 @@ int main(int argc, char **argv) {
logger.Info("Efficiency: {:.2f}%", output.efficiency * 100.f);
logger.Info("Max delay: {}",output.status.max_receive_delay);
logger.Info("Compression factor: {}x", output.status.compressed_ratio);
if (output.status.compressed_ratio)
logger.Info("Compression factor: {}x", output.status.compressed_ratio.value());
logger.Info("Receiving time: {} s", receiving_time);
logger.Info("Frame rate: {} Hz", static_cast<double>(nimages)/receiving_time);
logger.Info("Total throughput: {:.2f} GB/s",

View File

@@ -86,14 +86,13 @@ TEST_CASE("FastFeedbackIndexer","[Indexing]") {
IndexerWrapper wrapper;
wrapper.Setup(c);
auto ret = wrapper.Run(recip);
auto ret = wrapper.Run(recip, 0.05f);
REQUIRE(!ret.empty());
//auto uc = ret[0].GetUnitCell();
//REQUIRE(c.a == Approx(uc.a));
//REQUIRE(c.b == Approx(uc.b));
//REQUIRE(c.c == Approx(uc.c));
REQUIRE(ret[0].indexed_spots_count == recip.size());
double err[3] = {0.0, 0.0, 0.0};
for (const auto &iter: recip) {

View File

@@ -4,7 +4,7 @@
#include "../receiver/LossyFilter.h"
TEST_CASE("LossyFilterDisabled","[LossyFilter]") {
LossyFilter filter(false, 0.0);
LossyFilter filter(1.0);
DataMessage message{
.number = 123,
@@ -17,7 +17,7 @@ TEST_CASE("LossyFilterDisabled","[LossyFilter]") {
}
TEST_CASE("LossyFilterEnable","[LossyFilter]") {
LossyFilter filter(true, 0.0);
LossyFilter filter(0.0);
DataMessage message_1{
.number = 123,
@@ -36,7 +36,7 @@ TEST_CASE("LossyFilterEnable","[LossyFilter]") {
}
TEST_CASE("LossyFilterEnable_Prob","[LossyFilter]") {
LossyFilter filter(true, 0.5);
LossyFilter filter(0.5);
int64_t accepted = 0;
int64_t count = 0;

View File

@@ -18,3 +18,14 @@ TEST_CASE("PreviewCount", "") {
REQUIRE(counter.GeneratePreview());
REQUIRE(!counter.GeneratePreview()); // immediately after not
}
TEST_CASE("PreviewCount_Off", "") {
PreviewCounter counter(std::chrono::seconds(0));
REQUIRE(!counter.GeneratePreview()); // at first run it should generate preview
std::this_thread::sleep_for(std::chrono::seconds(1));
REQUIRE(!counter.GeneratePreview());
REQUIRE(!counter.GeneratePreview()); // immediately after not
REQUIRE(!counter.GeneratePreview()); // immediately after not
REQUIRE(!counter.GeneratePreview()); // immediately after not
}

View File

@@ -29,6 +29,9 @@ void print_usage() {
logger.Info("-p<int> min pixels per spot");
logger.Info("-M<int> max spots saved per image");
logger.Info("-o<string> output prefix");
logger.Info("-i<int> image number");
logger.Info("-I enable indexing");
logger.Info("-t<float> indexing tolerance (0.0-1.0)");
logger.Info("");
}
@@ -41,6 +44,15 @@ void GetGeometry(DiffractionExperiment &x, HDF5Object &master_file) {
HDF5DataSet(master_file, "/entry/instrument/detector/detector_distance").ReadScalar<double>() *1e3);
x.PhotonEnergy_keV(WVL_1A_IN_KEV /
HDF5DataSet(master_file, "/entry/instrument/beam/incident_wavelength").ReadScalar<double>());
try {
std::vector<float> v;
HDF5DataSet(master_file, "/entry/sample/unit_cell").ReadVector(v);
if (v.size() == 6) {
logger.Info("Unit cell {} {} {} {} {} {}", v[0], v[1], v[2], v[3], v[4], v[5]);
x.SetUnitCell(UnitCell(v[0], v[1], v[2], v[3], v[4], v[5]));
}
} catch (...) {}
}
int parse_options(DiffractionExperiment& experiment, SpotFindingSettings& settings,
@@ -57,12 +69,15 @@ int parse_options(DiffractionExperiment& experiment, SpotFindingSettings& settin
{"jpeg", no_argument, 0, 'j'},
{"help", no_argument, 0, 'h'},
{"max-spots", required_argument, 0, 'S'},
{"indexing", no_argument, 0, 'I'},
{"indexing_tolerance", required_argument, 0, 't'},
{"nimages", required_argument, 0, 'n'},
{0, 0, 0, 0}
};
int option_index = 0;
int opt;
while ((opt = getopt_long(argc, argv, "c:s:o:p:d:D:j?hS:",long_options, &option_index)) != -1 ) {
while ((opt = getopt_long(argc, argv, "c:s:o:p:d:D:j?hS:Ii:t:",long_options, &option_index)) != -1 ) {
switch (opt) {
case 'j':
write_jpeg = true;
@@ -88,6 +103,15 @@ int parse_options(DiffractionExperiment& experiment, SpotFindingSettings& settin
case 'S':
experiment.MaxSpotCount(atol(optarg));
break;
case 'I':
settings.indexing = true;
break;
case 'i':
experiment.ImagesPerTrigger(atol(optarg));
break;
case 't':
settings.indexing_tolerance = atof(optarg);
break;
case '?':
case 'h':
print_usage();
@@ -105,9 +129,10 @@ int main(int argc, char **argv) {
RegisterHDF5Filter();
DiffractionExperiment x;
x.FilePrefix("output");
x.FilePrefix("output").ImagesPerTrigger(100000).Summation(1).NumTriggers(1);
SpotFindingSettings settings = DiffractionExperiment::DefaultDataProcessingSettings();
settings.indexing = false;
int first_argc = parse_options(x, settings, argc, argv);
try {
@@ -145,7 +170,7 @@ int main(int argc, char **argv) {
GetGeometry(x, master_file);
uint64_t nimages = file_space.GetDimensions()[0];
uint64_t nimages = std::min<uint64_t>(file_space.GetDimensions()[0], x.GetImageNum());
logger.Info("Number of images in the original dataset: " + std::to_string(nimages));
x.Mode(DetectorMode::Conversion).ImagesPerTrigger(nimages);
@@ -173,6 +198,7 @@ int main(int argc, char **argv) {
MXAnalyzer analyzer(x);
PreviewImage preview(x);
uint64_t indexed_images = 0;
for (int i = 0; i < nimages; i++) {
std::vector<hsize_t> start = {(hsize_t) i,0,0};
@@ -198,6 +224,8 @@ int main(int argc, char **argv) {
message.number = i;
analyzer.Process(message, settings);
if (message.indexing_result)
indexed_images++;
if (write_jpeg) {
preview.UpdateImage(transformation.GetImage(), message.spots);
@@ -215,5 +243,7 @@ int main(int argc, char **argv) {
fileset->Write(message);
}
fileset.reset();
if (settings.indexing)
logger.Info("Indexing rate: {:0.1f}%", static_cast<float>(indexed_images) / static_cast<float>(nimages) * 100.0f);
}

View File

@@ -7,7 +7,6 @@ ZMQImagePuller::ZMQImagePuller(ZMQContext &context, const std::string &repub_add
socket (context, ZMQSocketType::Pull) {
socket.ReceiveWaterMark(ReceiverWaterMark);
socket.ReceiveTimeout(ReceiveTimeout);
zmq_recv_buffer.reserve(2*1024*1024); // Reasonable size 2 MiB
if (!repub_address.empty()) {
repub_socket = std::make_unique<ZMQSocket>(context, ZMQSocketType::Push);
@@ -36,16 +35,15 @@ void ZMQImagePuller::Abort() {
}
bool ZMQImagePuller::WaitForImage() {
int64_t msg_size = -1;
bool received;
while ((msg_size < 0) && (!abort))
msg_size = socket.Receive(zmq_recv_buffer, true, true);
do
received = socket.Receive(msg, true);
while (!received && !abort);
if (msg_size > 0) {
deserializer.Process(zmq_recv_buffer);
if (received) {
deserializer.Process(msg.data(), msg.size());
switch (deserializer.GetType()) {
case CBORStream2Deserializer::Type::START:
start_message = std::make_unique<StartMessage>(deserializer.GetStartMessage());
end_message.reset();
@@ -66,8 +64,7 @@ bool ZMQImagePuller::WaitForImage() {
if (repub_socket) {
// Republishing is non-blocking for images
// and blocking (with 100ms timeout) for START/END
repub_socket->Send(zmq_recv_buffer.data(), zmq_recv_buffer.size(),
deserializer.GetType() != CBORStream2Deserializer::Type::IMAGE);
repub_socket->Send(msg.data(), msg.size(), deserializer.GetType() != CBORStream2Deserializer::Type::IMAGE);
}
return true;

View File

@@ -11,7 +11,6 @@
#include "../frame_serialize/CBORStream2Deserializer.h"
class ZMQImagePuller {
std::vector<uint8_t> zmq_recv_buffer;
CBORStream2Deserializer deserializer;
@@ -29,6 +28,8 @@ class ZMQImagePuller {
std::unique_ptr<CompressedImage> calibration_message;
std::unique_ptr<ZMQSocket> repub_socket;
ZMQMessage msg;
public:
ZMQImagePuller(ZMQContext &context, const std::string &repub_address = "");
void Connect(const std::string &in_address);