diff --git a/common/ScalingSettings.h b/common/ScalingSettings.h index a20bcd21..f7c724f8 100644 --- a/common/ScalingSettings.h +++ b/common/ScalingSettings.h @@ -6,7 +6,7 @@ #include #include "JFJochException.h" -enum class PartialityModel { Fixed, Rotation, Unity, Postrefinement }; +enum class PartialityModel { Fixed, Rotation, Unity }; enum class IntensityFormat { Text, mmCIF, MTZ}; class ScalingSettings { diff --git a/image_analysis/scale_merge/ScaleOnTheFly.cpp b/image_analysis/scale_merge/ScaleOnTheFly.cpp index f5765374..46fc6cea 100644 --- a/image_analysis/scale_merge/ScaleOnTheFly.cpp +++ b/image_analysis/scale_merge/ScaleOnTheFly.cpp @@ -105,145 +105,6 @@ namespace { double partiality; }; -struct ScalingPostRefResidual : public ScalingResidual { - // Postrefinement for still images - // - // In this algorithm at the point of post-refinement we don't anymore care for where maximum of - // the reflection was located and if it fits the observed position. - // This reflections was already integrated and we cannot integrate it better at this point. - // But, we could adjust partiality to indicate that this reflection was wrongly predicted. - // I.e., integrated position was far away from true reflection, so partiality must be low. - // This is an empiric model and need to see if this will work in practice at all. - // I hope it will allow the model to find that reflections were misindexed. - // We assume at this point that initial indexing was done properly and integration was generally OK => most low resolution reflections fit correctly - // Yet we know, that small errors in indexing are inducing misalignment at high resolution - sometimes it is visible that high-resolution reflections - // are away from shoe-boxes observed in the image, if we can catch this at post-refinement/scaling step, this would be great. - // Next logical step is to do this pixel-wise - for each pixel refine partiality and merge pixels - // This should work for per-image scaling, or even, maybe, for full rotation datasets (3600 images) - // Then we could properly take into account misalignment of shoe-box center vs. partiality and also remove most pixels - // in the shoe-box that don't really contribute to the reflection. - // But...it could also drift to downweighting partiality for all high resolution reflections to make loss function "fake happy". - - // We assume rot3 == 0. Rot3 is not really helping much in crystallography (other than fixing polarization correction) - ScalingPostRefResidual(const Reflection &r, double Itrue, double sigma, - const DiffractionGeometry &geometry, - const CrystalLattice &lattice) - : ScalingResidual(r, Itrue, sigma), - integration_center_x(r.predicted_x), - integration_center_y(r.predicted_y), - inv_lambda(SafeInv(geometry.GetWavelength_A(), 0.0)), - pixel_size(geometry.GetPixelSize_mm()), - det_dist_mm(geometry.GetDetectorDistance_mm()), - beam_x(geometry.GetBeamX_pxl()), - beam_y(geometry.GetBeamY_pxl()), - exp_h(r.h), - exp_k(r.k), - exp_l(r.l), - Astar(lattice.Astar()), - Bstar(lattice.Bstar()), - Cstar(lattice.Cstar()), - c1(std::cos(geometry.GetPoniRot1_rad())), - s1(std::sin(geometry.GetPoniRot1_rad())), - c2(std::cos(geometry.GetPoniRot2_rad())), - s2(std::sin(geometry.GetPoniRot2_rad())) { - } - - template - T CalcPartiality(const T *const R, - const T *const beam_corr, - const T *const p0) const { - // Detector coordinates in mm - const T det_x = (T(integration_center_x) - beam_x - beam_corr[0]) * T(pixel_size); - const T det_y = (T(integration_center_y) - beam_y - beam_corr[1]) * T(pixel_size); - const T det_z = T(det_dist_mm); - - // Apply Ry(rot1) first: rotate around Y - const T t1_x = T(c1) * det_x + T(s1) * det_z; - const T t1_y = det_y; - const T t1_z = T(-s1) * det_x + T(c1) * det_z; - - // Then apply Rx(-rot2): rotate around X - const T x = t1_x; - const T y = T(c2) * t1_y + T(s2) * t1_z; - const T z = - T(s2) * t1_y + T(c2) * t1_z; - - // convert to recip space - const T lab_norm = ceres::sqrt(x * x + y * y + z * z); - const T inv_norm = T(1) / lab_norm; - - T recip_obs[3]; - recip_obs[0] = x * inv_norm * inv_lambda; - recip_obs[1] = y * inv_norm * inv_lambda; - recip_obs[2] = (z * inv_norm - T(1.0)) * inv_lambda; - - const Eigen::Matrix e_obs_recip(recip_obs[0], recip_obs[1], recip_obs[2]); - - const T astar_unrot[3] = {T(Astar.x), T(Astar.y), T(Astar.z)}; - const T bstar_unrot[3] = {T(Bstar.x), T(Bstar.y), T(Bstar.z)}; - const T cstar_unrot[3] = {T(Cstar.x), T(Cstar.y), T(Cstar.z)}; - - T astar_rot[3], bstar_rot[3], cstar_rot[3]; - - ceres::AngleAxisRotatePoint(p0, astar_unrot, astar_rot); - ceres::AngleAxisRotatePoint(p0, bstar_unrot, bstar_rot); - ceres::AngleAxisRotatePoint(p0, cstar_unrot, cstar_rot); - - const Eigen::Matrix e_pred_recip(T(exp_h) * astar_rot[0] + T(exp_k) * bstar_rot[0] + T(exp_l) * cstar_rot[0], - T(exp_h) * astar_rot[1] + T(exp_k) * bstar_rot[1] + T(exp_l) * cstar_rot[1], - T(exp_h) * astar_rot[2] + T(exp_k) * bstar_rot[2] + T(exp_l) * cstar_rot[2] - ); - - // Ewald sphere centre is at -k_i = (0, 0, -inv_lambda) - // Radial direction: outward normal at g_hkl - const Eigen::Matrix S_pred( - e_pred_recip[0], - e_pred_recip[1], - e_pred_recip[2] + T(inv_lambda) // g_hkl + k_i - ); - const T S_pred_norm = S_pred.norm(); - if (S_pred_norm < T(1e-10)) - return T(0); - - const Eigen::Matrix n_radial = S_pred / S_pred_norm; - - const Eigen::Matrix delta_q = e_obs_recip - e_pred_recip; - const T eps_radial = delta_q.dot(n_radial); - const Eigen::Matrix dq_tang = delta_q - eps_radial * n_radial; - const T eps_tangential_sq = dq_tang.squaredNorm(); // guaranteed ≥ 0 - // ───────────────────────────────────────────────────────────── - - return ceres::exp(- eps_radial * eps_radial / (R[0] * R[0]) - eps_tangential_sq / (R[1] * R[1])); - } - - template - bool operator()(const T *const G, - const T *const B, - const T *const R, - const T *const beam_corr, - const T *const p0, - T *residual) const { - if (R[0] < T(1e-10) || R[1] < T(1e-10)) - return false; - - const T B_term = ceres::exp(B[0] * T(b_resolution_coeff)); - - const T partiality = CalcPartiality(R, beam_corr, p0); - residual[0] = (G[0] * partiality * B_term * T(lp) * T(Itrue) - - T(Iobs)) * T(weight); - return true; - } - - const double integration_center_x, integration_center_y; - const double inv_lambda; - const double pixel_size; - const double det_dist_mm; - const double beam_x, beam_y; - const double exp_h; - const double exp_k; - const double exp_l; - const Coord Astar, Bstar, Cstar; - const double c1,s1,c2,s2; -}; struct RotationNormRegularizer { explicit RotationNormRegularizer(double weight) : weight(weight) {} @@ -282,7 +143,6 @@ bool ScaleOnTheFly::Accept(const Reflection &r) const { return std::isfinite(r.zeta) && r.zeta > 0.0f; case PartialityModel::Fixed: case PartialityModel::Unity: - case PartialityModel::Postrefinement: return true; } @@ -395,12 +255,6 @@ void ScaleOnTheFly::Scale(IntegrationOutcome &integration_outcome) const { problem.AddResidualBlock(cost, nullptr, &result.G, &result.B); } break; - case PartialityModel::Postrefinement: { - auto *cost = new ceres::AutoDiffCostFunction( - new ScalingPostRefResidual(r, Itrue, sigma, integration_outcome.geom, integration_outcome.latt)); - problem.AddResidualBlock(cost, nullptr, &result.G, &result.B, result.R, result.beam_corr, result.p0); - } - break; case PartialityModel::Rotation: { auto *cost = new ceres::AutoDiffCostFunction( new ScalingRotationResidual(r, Itrue, sigma)); @@ -460,21 +314,6 @@ void ScaleOnTheFly::Scale(IntegrationOutcome &integration_outcome) const { problem.SetParameterUpperBound(&result.mos, 0, s.GetMaxMosaicity()); } - if (model == PartialityModel::Postrefinement) { - problem.SetParameterLowerBound(result.R, 0, 1e-6); - problem.SetParameterLowerBound(result.R, 1, 1e-6); - - // Beam center can be off by max 1 pixel - problem.SetParameterLowerBound(result.beam_corr, 0, -1); - problem.SetParameterUpperBound(result.beam_corr, 0, 1); - problem.SetParameterLowerBound(result.beam_corr, 1, -1); - problem.SetParameterUpperBound(result.beam_corr, 1, 1); - - auto *angle_reg_cost = new ceres::AutoDiffCostFunction( - new RotationNormRegularizer(0.05)); - problem.AddResidualBlock(angle_reg_cost, nullptr, result.p0); - } - ceres::Solver::Options options; options.linear_solver_type = ceres::DENSE_QR; options.minimizer_progress_to_stdout = false; @@ -495,19 +334,12 @@ void ScaleOnTheFly::Scale(IntegrationOutcome &integration_outcome) const { r.partiality = RotationPartiality(r.delta_phi_deg, r.zeta, result.mos, result.wedge); break; } - case PartialityModel::Postrefinement: { - ScalingPostRefResidual residual(r, 0, 0, integration_outcome.geom, integration_outcome.latt); - r.partiality = static_cast(residual.CalcPartiality(result.R, result.beam_corr, result.p0)); - } - break; default: // For fixed partiality there is no need to change anything break; } const double denom = B_term * r.partiality * result.G; - if (std::isfinite(r.rlp) && - std::isfinite(denom) && - denom > 0.0) { + if (std::isfinite(r.rlp) && std::isfinite(denom) && denom > 0.0) { r.image_scale_corr = static_cast(r.rlp / denom); } else { r.image_scale_corr = NAN; diff --git a/tools/jfjoch_process.cpp b/tools/jfjoch_process.cpp index 47cd0736..5876418b 100644 --- a/tools/jfjoch_process.cpp +++ b/tools/jfjoch_process.cpp @@ -469,8 +469,6 @@ int main(int argc, char **argv) { partiality_model = PartialityModel::Fixed; else if (strcmp(optarg, "rot") == 0) partiality_model = PartialityModel::Rotation; - else if (strcmp(optarg, "postref") == 0) - partiality_model = PartialityModel::Postrefinement; else { logger.Error("Invalid partiality mode: {}", optarg); print_usage();