diff --git a/image_analysis/scale_merge/ScaleAndMerge.cpp b/image_analysis/scale_merge/ScaleAndMerge.cpp index 2d813178..f6ddef69 100644 --- a/image_analysis/scale_merge/ScaleAndMerge.cpp +++ b/image_analysis/scale_merge/ScaleAndMerge.cpp @@ -235,6 +235,7 @@ namespace { int img_id = 0; int hkl_slot = -1; double sigma = 0.0; + mutable bool selected = true; }; struct CorrectedObs { @@ -243,6 +244,104 @@ namespace { double sigma_corr; }; + void select_reflections_by_quasi_random(const std::vector &obs, + const ScaleMergeOptions &opt, + std::vector &hkl_selected, + double Isigma_cutoff = 1.0, + int min_per_bin = 10000, + int max_per_bin = 80000, + int n_resolution_bins = 20) { + float stat_d_min = std::numeric_limits::max(); + float stat_d_max = 0.0f; + + struct HKLStats { + int n = 0; + float d = std::numeric_limits::max(); + int shell_id = 0; + }; + const int nhkl = static_cast(hkl_selected.size()); + std::vector per_hkl(nhkl); + int reflection_above_cutoff = 0; + + for (const auto &o: obs) { + if (o.r->I / o.r->sigma < Isigma_cutoff) { + o.selected = false; + continue; + } + + const auto d = o.r->d; + reflection_above_cutoff += 1; + if (std::isfinite(d) && d > 0.0f) { + if (opt.d_min_limit_A > 0.0 && d < static_cast(opt.d_min_limit_A)) + continue; + stat_d_min = std::min(stat_d_min, d); + stat_d_max = std::max(stat_d_max, d); + auto &hs = per_hkl[o.hkl_slot]; + hs.n += 1; + hs.d = d; + } + } + + std::cout << "Reflections of I/sigma > " << Isigma_cutoff << " : " << reflection_above_cutoff << std::endl; + if (reflection_above_cutoff < min_per_bin * n_resolution_bins) { + std::cout << "No additional selection applied before scaling" << std::endl; + return; + } + + if (stat_d_min < stat_d_max && stat_d_min > 0.0f) { + const float d_min_pad = stat_d_min * 0.999f; + const float d_max_pad = stat_d_max * 1.001f; + ResolutionShells scaling_shells(d_min_pad, d_max_pad, n_resolution_bins); + + for (int h = 0; h < nhkl; ++h) { + const auto d = per_hkl[h].d; + if (std::isfinite(d) && d > 0.0f) { + if (opt.d_min_limit_A > 0.0 && d < static_cast(opt.d_min_limit_A)) + continue; + auto s = scaling_shells.GetShell(d); + if (s.has_value()) + per_hkl[h].shell_id = s.value(); + } + } + + // Accumulators per shell + struct ShellAccum { + int obs_unique = 0; + int obs_total = 0; + bool selected = true; + }; + std::vector shell_acc(n_resolution_bins); + + for (int h = 0; h < nhkl; ++h) { + auto &sa = shell_acc[per_hkl[h].shell_id]; + if (sa.obs_unique > min_per_bin * 1.2 || sa.obs_total > max_per_bin) + hkl_selected[h] = false; + else + sa.obs_unique += 1; + sa.obs_total += per_hkl[h].n; + } + + const auto shell_min_res = scaling_shells.GetShellMinRes(); + + std::cout << "| d-mean | n_refl_tot | n_refl_uni | selection |" << std::endl; + for (int n=0; n < n_resolution_bins; ++n) { + if (shell_acc[n].obs_unique < min_per_bin) + shell_acc[n].selected = false; + if (shell_min_res[n] <= 0.0f) continue; + std::cout << std::setw(8) << std::fixed << std::setprecision(3) << shell_min_res[n] + << std::setw(12) << std::fixed << std::setprecision(0) << shell_acc[n].obs_unique + << std::setw(12) << shell_acc[n].obs_total + << " " << shell_acc[n].selected << std::endl; + } + + for (int h = 0; h < nhkl; ++h) { + auto &sa = shell_acc[per_hkl[h].shell_id]; + if (!sa.selected) + hkl_selected[h] = false; + } + } + } + void scale(const ScaleMergeOptions &opt, std::vector &g, std::vector &mosaicity, @@ -250,10 +349,12 @@ namespace { const std::vector &image_slot_used, bool rotation_crystallography, size_t nhkl, - const std::vector &obs) { + const std::vector &obs, + bool selection) { ceres::Problem problem; std::vector Itrue(nhkl, 0.0); + std::vector hkl_selected(nhkl, true); // Initialize Itrue from per-HKL median of observed intensities { @@ -277,8 +378,12 @@ namespace { double wedge = opt.wedge_deg.value_or(0.0); + select_reflections_by_quasi_random(obs, opt, hkl_selected); + std::vector is_valid_hkl_slot(nhkl, false); for (const auto &o: obs) { + if (!o.selected) continue; + if (!hkl_selected[o.hkl_slot]) continue; switch (opt.partiality_model) { case ScaleMergeOptions::PartialityModel::Rotation: { auto *cost = new ceres::AutoDiffCostFunction( @@ -398,6 +503,7 @@ namespace { options.function_tolerance = 1e-4; ceres::Solver::Summary summary; + std::cout << "Now start the ceres-solver with residual blocks: " << problem.NumResidualBlocks() << std::endl; ceres::Solve(options, &problem, &summary); std::cout << summary.FullReport() << std::endl; } @@ -694,6 +800,7 @@ namespace { o.img_id = img_id; o.hkl_slot = hkl_slot; o.sigma = sigma; + o.selected = true; obs.push_back(o); } } @@ -707,6 +814,7 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector Fix reference unit cell: -C\"a,b,c,alpha,beta,gamma\" (comma-separated, no spaces; quotes optional)"); logger.Info(" -c Max spot count (default: 250)"); - logger.Info(" -W HDF5 file with analysis results is written"); + logger.Info(" -W HDF5 file with analysis results is written. 'l' or 'light' deactivates image-output"); logger.Info(" -T Noise sigma level for spot finding (default: 3.0)"); + logger.Info(" -a Use all reflection for scaling without filtering (default: false)"); } void trim_in_place(std::string& t) { @@ -163,6 +164,8 @@ int main(int argc, char **argv) { std::optional space_group_number; std::optional fixed_reference_unit_cell; bool write_output = false; + bool write_output_noimage = false; + bool filtering = true; std::optional max_spot_count_override; float sigma_spot_finding = 3.0; std::optional merging_threshold; @@ -179,7 +182,7 @@ int main(int argc, char **argv) { } int opt; - while ((opt = getopt(argc, argv, "o:N:s:e:vc:R::FX:xd:S:M::P:AD:C:T:W")) != -1) { + while ((opt = getopt(argc, argv, "o:N:s:e:vc:R::FX:xd:S:M::P:AD:C:T:W:a")) != -1) { switch (opt) { case 'o': output_prefix = optarg; @@ -195,6 +198,10 @@ int main(int argc, char **argv) { break; case 'W': write_output = true; + if (strcmp(optarg, "light") == 0 || strcmp(optarg, "l") == 0) { + write_output_noimage = true; + logger.Warning("Image data will not be saved."); + } break; case 'v': verbose = true; @@ -237,6 +244,9 @@ int main(int argc, char **argv) { case 'x': refine_beam_center = false; break; + case 'a': + filtering = false; + break; case 'D': d_min_scale_merge = atof(optarg); logger.Info("High resolution limit for scaling/merging set to {:.2f} A", d_min_spot_finding); @@ -334,8 +344,11 @@ int main(int argc, char **argv) { experiment.OverwriteExistingFiles(true); experiment.PolarizationFactor(0.99); - if (fixed_reference_unit_cell.has_value()) + if (fixed_reference_unit_cell.has_value()) { experiment.SetUnitCell(*fixed_reference_unit_cell); + } else { + experiment.SetUnitCell({}); + } if (max_spot_count_override.has_value()) { experiment.MaxSpotCount(max_spot_count_override.value()); @@ -346,6 +359,8 @@ int main(int argc, char **argv) { IndexingSettings indexing_settings; indexing_settings.Algorithm(indexing_algorithm); indexing_settings.RotationIndexing(rotation_indexing); + if (rotation_indexing) + logger.Info("Rotation indexing is activated."); if (rotation_indexing_range.has_value()) indexing_settings.RotationIndexingMinAngularRange_deg(rotation_indexing_range.value()); @@ -355,6 +370,22 @@ int main(int argc, char **argv) { indexing_settings.GeomRefinementAlgorithm(GeomRefinementAlgorithmEnum::None); experiment.ImportIndexingSettings(indexing_settings); + switch (experiment.GetIndexingAlgorithm()) { + case IndexingAlgorithmEnum::FFBIDX: + logger.Info("Indexer used: FFBIDX"); + break; + case IndexingAlgorithmEnum::FFTW: + logger.Info("Indexer used: FFTW"); + break; + case IndexingAlgorithmEnum::FFT: + logger.Info("Indexer used: FFT (CUDA)"); + break; + case IndexingAlgorithmEnum::None: + logger.Warning("Indexer not defined!"); + return 0; + default: ; + } + SpotFindingSettings spot_settings; spot_settings.enable = true; spot_settings.indexing = true; @@ -435,6 +466,10 @@ int main(int argc, char **argv) { compressed_buffer.resize(MaxCompressedSize(experiment.GetCompressionAlgorithm(), experiment.GetPixelsNum(), experiment.GetByteDepthImage())); + auto size = compressor.Compress(compressed_buffer.data(), + compressed_buffer.data(), + experiment.GetPixelsNum(), + sizeof(uint8_t)); // Thread-local analysis resources MXAnalysisWithoutFPGA analysis(experiment, mapping, pixel_mask, indexer); @@ -485,10 +520,11 @@ int main(int argc, char **argv) { auto image_end_time = std::chrono::high_resolution_clock::now(); std::chrono::duration image_duration = image_end_time - image_start_time; - auto size = compressor.Compress(compressed_buffer.data(), - img->Image().data(), - experiment.GetPixelsNum(), - sizeof(int32_t)); + if (!write_output_noimage) + size = compressor.Compress(compressed_buffer.data(), + img->Image().data(), + experiment.GetPixelsNum(), + sizeof(int32_t)); msg.image = CompressedImage(compressed_buffer.data(), size, experiment.GetXPixelsNum(), @@ -521,7 +557,7 @@ int main(int argc, char **argv) { } // Progress log - if (current_idx_offset > 0 && current_idx_offset % 100 == 0) { + if ((current_idx_offset > 0 && (current_idx_offset+1) % 100 == 0) || image_idx == end_image - 1) { std::optional indexing_rate; { std::lock_guard lock(plots_mutex); @@ -530,11 +566,11 @@ int main(int argc, char **argv) { if (indexing_rate.has_value()) { logger.Info("Processed {} / {} images (indexing rate {:.1f}%)", - current_idx_offset, images_to_process, + current_idx_offset+1, images_to_process, indexing_rate.value() * 100.0f); } else { logger.Info("Processed {} / {} images (indexing rate N/A)", - current_idx_offset, images_to_process); + current_idx_offset+1, images_to_process); } if (image_idx == end_image - 1) { if (!merging_threshold.has_value()) merging_threshold = 0.5f; @@ -596,7 +632,6 @@ int main(int argc, char **argv) { logger.Info("Rotation Indexing found lattice"); } - // --- Optional: run scaling (mosaicity refinement) on accumulated reflections --- // --- Optional: run scaling (mosaicity refinement) on accumulated reflections --- if (run_scaling) { logger.Info("Running scaling (mosaicity refinement) ..."); @@ -607,6 +642,7 @@ int main(int argc, char **argv) { scale_opts.max_solver_time_s = 240.0; // generous cutoff for now scale_opts.merge_friedel = !anomalous_mode; scale_opts.d_min_limit_A = d_min_scale_merge.value_or(0.0); + scale_opts.selection = filtering; const bool fixed_space_group = space_group || experiment.GetGemmiSpaceGroup().has_value(); @@ -728,6 +764,8 @@ int main(int argc, char **argv) { cif_meta.unit_cell = rotation_indexer_ret->lattice.GetUnitCell(); } else if (experiment.GetUnitCell().has_value()) { cif_meta.unit_cell = experiment.GetUnitCell().value(); + } else { + logger.Warning("No UnitCell output"); } if (scale_opts.space_group.has_value()) {