VERSION: 1.0.0-rc.144

image_analysis/scale_merge/ScaleAndMerge.cpp

@@ -122,13 +122,11 @@ namespace {
               delta_phi_(r.delta_phi_deg),
               lp_(SafeInv(r.rlp, 1.0)),
               c1_(r.zeta / std::sqrt(2.0)),
-              inv_2d2_(0.5 / (static_cast<double>(r.d) * static_cast<double>(r.d))),
               partiality_(r.partiality) {
         }
 
         template<typename T>
         bool operator()(const T *const G,
-                        const T *const B,
                         const T *const mosaicity,
                         const T *const Itrue,
                         const T *const wedge,
@@ -142,8 +140,7 @@ namespace {
             } else
                 partiality = T(1.0);
 
-            const T bscale = ceres::exp(-B[0] * T(inv_2d2_));
-            const T Ipred = G[0] * bscale * partiality * T(lp_) * Itrue[0];
+            const T Ipred = G[0] * partiality * T(lp_) * Itrue[0];
             residual[0] = (Ipred - T(Iobs_)) * T(weight_);
             return true;
         }
@@ -153,7 +150,6 @@ namespace {
         double delta_phi_;
         double lp_;
         double c1_;
-        double inv_2d2_;
         double partiality_;
     };
 
@@ -162,19 +158,16 @@ namespace {
             : Iobs_(static_cast<double>(r.I)),
               weight_(SafeInv(sigma_obs, 1.0)),
               lp_(SafeInv(r.rlp, 1.0)),
-              dist_ewald_sq_(r.dist_ewald * r.dist_ewald),
-              inv_2d2_(0.5 / (static_cast<double>(r.d) * static_cast<double>(r.d))) {
+              dist_ewald_sq_(r.dist_ewald * r.dist_ewald) {
         }
 
         template<typename T>
         bool operator()(const T *const G,
-                        const T *const B,
                         const T *const R,
                         const T *const Itrue,
                         T *residual) const {
-            const T partiality = ceres::exp(-T(dist_ewald_sq_) / R[0]);
-            const T bscale = ceres::exp(-B[0] * T(inv_2d2_));
-            const T Ipred = G[0] * bscale * partiality * T(lp_) * Itrue[0];
+            const T partiality = ceres::exp(-T(dist_ewald_sq_)/R[0]);
+            const T Ipred = G[0] * partiality * T(lp_) * Itrue[0];
             residual[0] = (Ipred - T(Iobs_)) * T(weight_);
             return true;
         }
@@ -183,24 +176,20 @@ namespace {
         double weight_;
         double lp_;
         double dist_ewald_sq_;
-        double inv_2d2_;
     };
 
     struct IntensityFixedResidual {
         IntensityFixedResidual(const Reflection &r, double sigma_obs, double partiality)
             : Iobs_(static_cast<double>(r.I)),
               weight_(SafeInv(sigma_obs, 1.0)),
-              corr_(partiality * SafeInv(r.rlp, 1.0)),
-              inv_2d2_(0.5 / (static_cast<double>(r.d) * static_cast<double>(r.d))) {
-        }
+              corr_(partiality * SafeInv(r.rlp, 1.0))
+        {}
 
         template<typename T>
         bool operator()(const T *const G,
-                        const T *const B,
                         const T *const Itrue,
                         T *residual) const {
-            const T bscale = ceres::exp(-B[0] * T(inv_2d2_));
-            const T Ipred = T(corr_) * G[0] * bscale * Itrue[0];
+            const T Ipred = T(corr_) * G[0] * Itrue[0];
             residual[0] = (Ipred - T(Iobs_)) * T(weight_);
             return true;
         }
@@ -208,7 +197,6 @@ namespace {
         double Iobs_;
         double weight_;
         double corr_;
-        double inv_2d2_;
     };
 
     struct ScaleRegularizationResidual {
@@ -242,21 +230,6 @@ namespace {
         double inv_sigma_;
     };
 
-    struct ValueRegularizationResidual {
-        ValueRegularizationResidual(double target, double sigma)
-            : target_(target), inv_sigma_(SafeInv(sigma, 1.0)) {
-        }
-
-        template<typename T>
-        bool operator()(const T *const x, T *residual) const {
-            residual[0] = (x[0] - T(target_)) * T(inv_sigma_);
-            return true;
-        }
-
-        double target_;
-        double inv_sigma_;
-    };
-
     struct ObsRef {
         const Reflection *r = nullptr;
         int img_id = 0;
@@ -275,7 +248,6 @@ namespace {
                std::vector<double> &g,
                std::vector<double> &mosaicity,
                std::vector<double> &R_sq,
-               std::vector<double> &image_bfactor,
                const std::vector<uint8_t> &image_slot_used,
                bool rotation_crystallography,
                size_t nhkl,
@@ -310,45 +282,41 @@ namespace {
         for (const auto &o: obs) {
             switch (opt.partiality_model) {
                 case ScaleMergeOptions::PartialityModel::Rotation: {
-                    auto *cost = new ceres::AutoDiffCostFunction<IntensityRotResidual, 1, 1, 1, 1, 1, 1>(
+                    auto *cost = new ceres::AutoDiffCostFunction<IntensityRotResidual, 1, 1, 1, 1, 1>(
                         new IntensityRotResidual(*o.r, o.sigma, opt.wedge_deg.value_or(0.0)));
                     problem.AddResidualBlock(cost,
                                              nullptr,
                                              &g[o.img_id],
-                                             &image_bfactor[o.img_id],
                                              &mosaicity[o.img_id],
                                              &Itrue[o.hkl_slot],
                                              &wedge);
                 }
                 break;
                 case ScaleMergeOptions::PartialityModel::Still: {
-                    auto *cost = new ceres::AutoDiffCostFunction<IntensityStillResidual, 1, 1, 1, 1, 1>(
+                    auto *cost = new ceres::AutoDiffCostFunction<IntensityStillResidual, 1, 1, 1, 1>(
                         new IntensityStillResidual(*o.r, o.sigma));
                     problem.AddResidualBlock(cost,
                                              nullptr,
                                              &g[o.img_id],
-                                             &image_bfactor[o.img_id],
                                              &R_sq[o.img_id],
                                              &Itrue[o.hkl_slot]);
                 }
                 break;
                 case ScaleMergeOptions::PartialityModel::Unity: {
-                    auto *cost = new ceres::AutoDiffCostFunction<IntensityFixedResidual, 1, 1, 1, 1>(
+                    auto *cost = new ceres::AutoDiffCostFunction<IntensityFixedResidual, 1, 1, 1>(
                         new IntensityFixedResidual(*o.r, o.sigma, 1.0));
                     problem.AddResidualBlock(cost,
                                              nullptr,
                                              &g[o.img_id],
-                                             &image_bfactor[o.img_id],
                                              &Itrue[o.hkl_slot]);
                 }
                 break;
                 case ScaleMergeOptions::PartialityModel::Fixed: {
-                    auto *cost = new ceres::AutoDiffCostFunction<IntensityFixedResidual, 1, 1, 1, 1>(
+                    auto *cost = new ceres::AutoDiffCostFunction<IntensityFixedResidual, 1, 1, 1>(
                         new IntensityFixedResidual(*o.r, o.sigma, o.r->partiality));
                     problem.AddResidualBlock(cost,
                                              nullptr,
                                              &g[o.img_id],
-                                             &image_bfactor[o.img_id],
                                              &Itrue[o.hkl_slot]);
                 }
                 break;
@@ -359,12 +327,8 @@ namespace {
         for (int i = 0; i < g.size(); ++i) {
             if (image_slot_used[i]) {
                 auto *cost = new ceres::AutoDiffCostFunction<ScaleRegularizationResidual, 1, 1>(
-                    new ScaleRegularizationResidual(opt.scale_regularization_sigma));
+                    new ScaleRegularizationResidual(0.05));
                 problem.AddResidualBlock(cost, nullptr, &g[i]);
-
-                auto *bcost = new ceres::AutoDiffCostFunction<ValueRegularizationResidual, 1, 1>(
-                    new ValueRegularizationResidual(0.0, opt.image_bfactor_regularization_sigma_A2));
-                problem.AddResidualBlock(bcost, nullptr, &image_bfactor[i]);
             }
         }
 
@@ -400,20 +364,13 @@ namespace {
             }
         }
 
+        // Scaling factors must be always positive
         for (int i = 0; i < g.size(); i++) {
             if (image_slot_used[i])
                 problem.SetParameterLowerBound(&g[i], 0, 1e-12);
         }
 
-        for (int i = 0; i < image_bfactor.size(); ++i) {
-            if (!image_slot_used[i])
-                continue;
-            problem.SetParameterLowerBound(&image_bfactor[i], 0, opt.image_bfactor_min_A2);
-            problem.SetParameterUpperBound(&image_bfactor[i], 0, opt.image_bfactor_max_A2);
-            if (!opt.refine_image_bfactor)
-                problem.SetParameterBlockConstant(&image_bfactor[i]);
-        }
-
+        // Mosaicity refinement + bounds
         if (opt.partiality_model == ScaleMergeOptions::PartialityModel::Rotation) {
             for (int i = 0; i < mosaicity.size(); ++i) {
                 if (image_slot_used[i]) {
@@ -427,9 +384,10 @@ namespace {
                 problem.SetParameterLowerBound(&wedge, 0, 0.0);
         }
 
+        // use all available threads
         unsigned int hw = std::thread::hardware_concurrency();
         if (hw == 0)
-            hw = 1;
+            hw = 1; // fallback
 
         ceres::Solver::Options options;
 
@@ -630,14 +588,14 @@ namespace {
             }
         }
     }
-
     void calc_obs(const ScaleMergeOptions &opt,
-                  std::vector<double> &g,
-                  std::vector<double> &mosaicity,
-                  std::vector<double> &R_sq,
-                  std::vector<double> &image_bfactor,
-                  const std::vector<ObsRef> &obs,
-                  std::vector<CorrectedObs> &corr_obs) {
+               std::vector<double> &g,
+               std::vector<double> &mosaicity,
+               std::vector<double> &R_sq,
+               const std::vector<ObsRef> &obs,
+               std::vector<CorrectedObs> &corr_obs) {
+
+        // ---- Compute corrected observations once (used for both merging and statistics) ----
         const double half_wedge = opt.wedge_deg.value_or(0.0) / 2.0;
 
         for (const auto &o: obs) {
@@ -645,6 +603,7 @@ namespace {
             const double lp = SafeInv(static_cast<double>(r.rlp), 1.0);
             const double G_i = g[o.img_id];
 
+            // Compute partiality with refined mosaicity
             double partiality = 1.0;
 
             switch (opt.partiality_model) {
@@ -657,7 +616,7 @@ namespace {
                     const double arg_minus = (r.delta_phi_deg - half_wedge) * c1 / mosaicity[o.img_id];
                     partiality = (std::erf(arg_plus) - std::erf(arg_minus)) / 2.0;
                 }
-                break;
+                    break;
                 case ScaleMergeOptions::PartialityModel::Still:
                     partiality = std::exp(-r.dist_ewald * r.dist_ewald / R_sq[o.img_id]);
                     break;
@@ -668,8 +627,7 @@ namespace {
             if (partiality <= opt.min_partiality_for_merge)
                 continue;
 
-            const double bscale = std::exp(-image_bfactor[o.img_id] / (2.0 * r.d * r.d));
-            const double correction = G_i * bscale * partiality * lp;
+            const double correction = G_i * partiality * lp;
             if (correction <= 0.0)
                 continue;
 
@@ -767,32 +725,28 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector<std::vector<Ref
     std::vector<double> g(n_image_slots, 1.0);
     std::vector<double> mosaicity(n_image_slots, opt.mosaicity_init_deg);
     std::vector<double> R_sq(n_image_slots, 0.001 * 0.001);
-    std::vector<double> image_bfactor(n_image_slots, opt.image_bfactor_init_A2);
 
     for (int i = 0; i < n_image_slots; i++) {
         if (!image_slot_used[i]) {
             mosaicity[i] = NAN;
             g[i] = NAN;
             R_sq[i] = NAN;
-            image_bfactor[i] = NAN;
         } else if (opt.mosaicity_init_deg_vec.size() > i && std::isfinite(opt.mosaicity_init_deg_vec[i])) {
             mosaicity[i] = opt.mosaicity_init_deg_vec[i];
         }
     }
 
-    scale(opt, g, mosaicity, R_sq, image_bfactor, image_slot_used, rotation_crystallography, nhkl, obs);
+    scale(opt, g, mosaicity, R_sq, image_slot_used, rotation_crystallography, nhkl, obs);
 
     ScaleMergeResult out;
 
     out.image_scale_g.resize(observations.size(), NAN);
     out.mosaicity_deg.resize(observations.size(), NAN);
-    out.image_bfactor.resize(observations.size(), NAN);
     for (int i = 0; i < observations.size(); i++) {
         size_t img_slot = i / opt.image_cluster;
         if (image_slot_used[img_slot]) {
             out.image_scale_g[i] = g[img_slot];
             out.mosaicity_deg[i] = mosaicity[img_slot];
-            out.image_bfactor[i] = image_bfactor[img_slot];
         }
     }
 
@@ -810,6 +764,7 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector<std::vector<Ref
         out.merged[h].d = 0.0;
     }
 
+    // Populate d from median of observations per HKL
     {
         std::vector<std::vector<double> > per_hkl_d(nhkl);
         for (const auto &o: obs) {
@@ -826,7 +781,7 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector<std::vector<Ref
         }
     }
 
-    calc_obs(opt, g, mosaicity, R_sq, image_bfactor, obs, corr_obs);
+    calc_obs(opt, g, mosaicity, R_sq, obs, corr_obs);
     merge(nhkl, out, corr_obs);
     stats(opt, nhkl, out, corr_obs);
 

image_analysis/scale_merge/ScaleAndMerge.h

@@ -50,12 +50,6 @@ struct ScaleMergeOptions {
 
     bool refine_wedge = false;
 
-    bool refine_image_bfactor = false;
-    double image_bfactor_init_A2 = 0.0;
-    double image_bfactor_min_A2 = -50.0;
-    double image_bfactor_max_A2 = 100.0;
-    double image_bfactor_regularization_sigma_A2 = 20.0;
-
     enum class PartialityModel {Fixed, Rotation, Unity, Still} partiality_model = PartialityModel::Fixed;
 };
 
@@ -91,7 +85,7 @@ struct ScaleMergeResult {
     std::vector<MergedReflection> merged;
     std::vector<float> image_scale_g;
     std::vector<float> mosaicity_deg;
-    std::vector<float> image_bfactor;
+
     /// Per-shell and overall merging statistics (populated after merging)
     MergeStatistics statistics;
 };

tools/jfjoch_process.cpp

@@ -37,22 +37,21 @@ void print_usage(Logger &logger) {
     logger.Info("   -s<num>    Start image number (default: 0)");
     logger.Info("   -e<num>    End image number (default: all)");
     logger.Info("   -v         Verbose output");
-    logger.Info("   -d<num>    High resolution limit for spot finding (default: 1.5)");
-    logger.Info("   -T<num>    Noise sigma level for spot finding (default: 3.0)");
-    logger.Info("   -t<num>    Photon count threshold for spot finding (default: 10)");
-    logger.Info("   -c<num>    Max spot count (default: 250)");
     logger.Info("   -R[num]    Rotation indexing (optional: min angular range deg)");
     logger.Info("   -F         Use FFT indexing algorithm (shortcut for -XFFT)");
     logger.Info("   -X<txt>    Indexing algorithm (FFBIDX|FFT|FFTW|Auto|None)");
     logger.Info("   -x         No least-square beam center refinement");
-    logger.Info("   -C<cell>   Fix reference unit cell: -C\"a,b,c,alpha,beta,gamma\" (comma-separated, no spaces; quotes optional)");
-    logger.Info("   -M         Scale and merge (refine mosaicity) and write scaled.hkl + image.dat");
-    logger.Info("   -B         Refine per image B-factor in scaling");
-    logger.Info("   -S<num>    Space group number");
-    logger.Info("   -P<txt>    Partiality refinement fixed|rot|unity (default: fixed)");
+    logger.Info("   -d<num>    High resolution limit for spot finding (default: 1.5)");
     logger.Info("   -D<num>    High resolution limit for scaling/merging (default: 0.0; no limit)");
+    logger.Info("   -S<num>    Space group number");
+    logger.Info("   -M         Scale and merge (refine mosaicity) and write scaled.hkl + image.dat");
+    logger.Info("   -P<txt>    Partiality refinement fixed|rot|unity (default: fixed)");
     logger.Info("   -A         Anomalous mode (don't merge Friedel pairs)");
+    logger.Info("   -C<cell>   Fix reference unit cell: -C\"a,b,c,alpha,beta,gamma\" (comma-separated, no spaces; quotes optional)");
+    logger.Info("   -c<num>    Max spot count (default: 250)");
     logger.Info("   -W         HDF5 file with analysis results is written");
+    logger.Info("   -T<num>    Noise sigma level for spot finding (default: 3.0)");
+    logger.Info("   -t<num>    Photon count threshold for spot finding (default: 10)");
 }
 
 void trim_in_place(std::string& t) {
@@ -138,7 +137,6 @@ int main(int argc, char **argv) {
     std::optional<int64_t> max_spot_count_override;
     float sigma_spot_finding = 3.0;
     int64_t photon_count_threshold_spot_finding = 10;
-    bool refine_b_factor = false;
 
     IndexingAlgorithmEnum indexing_algorithm = IndexingAlgorithmEnum::Auto;
 
@@ -153,7 +151,7 @@ int main(int argc, char **argv) {
     }
 
     int opt;
-    while ((opt = getopt(argc, argv, "o:N:s:e:vc:R::FX:xd:S:MP:AD:C:T:t:WB")) != -1) {
+    while ((opt = getopt(argc, argv, "o:N:s:e:vc:R::FX:xd:S:MP:AD:C:T:t:W")) != -1) {
         switch (opt) {
             case 'o':
                 output_prefix = optarg;
@@ -179,9 +177,6 @@ int main(int argc, char **argv) {
             case 'c':
                 max_spot_count_override = atoll(optarg);
                 break;
-            case 'B':
-                refine_b_factor = true;
-                break;
             case 'R':
                 rotation_indexing = true;
                 if (optarg) rotation_indexing_range = atof(optarg);
@@ -552,7 +547,6 @@ 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.refine_image_bfactor = refine_b_factor;
 
         const bool fixed_space_group = space_group || experiment.GetGemmiSpaceGroup().has_value();