From c6115928bb875b0ed408655ff12b9ca00dae017c Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 7 Jun 2023 14:02:26 +0200 Subject: Move Mille stuff to separate file --- libcrystfel/src/crystfel-mille.c | 168 +++++++++++++++++++++++++++++++++++++++ libcrystfel/src/crystfel-mille.h | 61 ++++++++++++++ libcrystfel/src/predict-refine.c | 123 +++------------------------- libcrystfel/src/predict-refine.h | 23 ++++-- 4 files changed, 257 insertions(+), 118 deletions(-) create mode 100644 libcrystfel/src/crystfel-mille.c create mode 100644 libcrystfel/src/crystfel-mille.h (limited to 'libcrystfel/src') diff --git a/libcrystfel/src/crystfel-mille.c b/libcrystfel/src/crystfel-mille.c new file mode 100644 index 00000000..c06dff01 --- /dev/null +++ b/libcrystfel/src/crystfel-mille.c @@ -0,0 +1,168 @@ +/* + * crystfel-mille.c + * + * Interface to Millepede geometry refinement + * + * Copyright © 2023 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. + * + * Authors: + * 2023 Thomas White + * + * This file is part of CrystFEL. + * + * CrystFEL is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * CrystFEL is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with CrystFEL. If not, see . + * + */ + +#include + +#include +#include + +#include "image.h" +#include "geometry.h" +#include "cell-utils.h" +#include "predict-refine.h" +#include "profile.h" + +#include + + +static const enum gparam rv[] = +{ + GPARAM_ASX, + GPARAM_ASY, + GPARAM_ASZ, + GPARAM_BSX, + GPARAM_BSY, + GPARAM_BSZ, + GPARAM_CSX, + GPARAM_CSY, + GPARAM_CSZ, +}; + + +int mille_label(int hierarchy_level, int member_index, char param) +{ + int label; + + assert(member_index < 1000); + + label = 100000*hierarchy_level + 100*member_index; + switch ( param ) { + case 'x' : return label+1; + case 'y' : return label+2; + case 'z' : return label+3; + default : abort(); + } +} + + +void write_mille(Mille *mille, int n, UnitCell *cell, + struct reflpeak *rps, struct image *image, + double dx, double dy, + const struct detgeom_panel_group *group, + int hierarchy_level, int member_index) +{ + int i; + float local_gradients[9]; + float global_gradients[6]; + int labels[6]; + + /* Spot x-position terms */ + for ( i=0; idetgeom, dx, dy), + 0.65*rps[i].panel->pixel_pitch); + } + + /* Spot y-position terms */ + for ( i=0; idetgeom, dx, dy), + 0.65*rps[i].panel->pixel_pitch); + } + + /* Next level of hierarchy */ + for ( i=0; in_children; i++ ) { + write_mille(mille, n, cell, rps, image, dx, dy, + group->children[i], hierarchy_level+1, i); + } + +} + + +Mille *crystfel_mille_new(const char *outFileName, + int asBinary, + int writeZero) +{ + return mille_new(outFileName, asBinary, writeZero); +} + + +void crystfel_mille_free(Mille *m) +{ + mille_free(m); +} + + +void crystfel_mille_delete_last_record(Mille *m) +{ + mille_delete_last_record(m); +} + + +void crystfel_mille_write_record(Mille *m) +{ + mille_write_record(m); +} diff --git a/libcrystfel/src/crystfel-mille.h b/libcrystfel/src/crystfel-mille.h new file mode 100644 index 00000000..7ba8de61 --- /dev/null +++ b/libcrystfel/src/crystfel-mille.h @@ -0,0 +1,61 @@ +/* + * crystfel-mille.h + * + * Interface to Millepede geometry refinement + * + * Copyright © 2023 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. + * + * Authors: + * 2023 Thomas White + * + * This file is part of CrystFEL. + * + * CrystFEL is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * CrystFEL is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with CrystFEL. If not, see . + * + */ + +#ifndef CRYSTFEL_MILLE_H +#define CRYSTFEL_MILLE_H + +typedef void *Mille; + +#include "cell.h" +#include "image.h" +#include "predict-refine.h" + +/** + * \file crystfel-mille.h + * Detector geometry refinement using Millepede + */ + +extern Mille *crystfel_mille_new(const char *outFileName, + int asBinary, + int writeZero); + +extern void crystfel_mille_free(Mille *m); + +extern int mille_label(int hierarchy_level, int member_index, char param); + +extern void write_mille(Mille *mille, int n, UnitCell *cell, + struct reflpeak *rps, struct image *image, + double dx, double dy, + const struct detgeom_panel_group *group, + int hierarchy_level, int member_index); + +extern void crystfel_mille_delete_last_record(Mille *m); + +extern void crystfel_mille_write_record(Mille *m); + +#endif /* CRYSTFEL_MILLE_H */ diff --git a/libcrystfel/src/predict-refine.c b/libcrystfel/src/predict-refine.c index 86a37014..585df4f8 100644 --- a/libcrystfel/src/predict-refine.c +++ b/libcrystfel/src/predict-refine.c @@ -39,10 +39,7 @@ #include "cell-utils.h" #include "predict-refine.h" #include "profile.h" - -#ifdef HAVE_MILLEPEDE -#include -#endif +#include "crystfel-mille.h" /** \file predict-refine.h */ @@ -70,15 +67,6 @@ static const enum gparam rv[] = GPARAM_DETY, }; -struct reflpeak { - Reflection *refl; - struct imagefeature *peak; - double Ih; /* normalised */ - struct detgeom_panel *panel; /* panel the reflection appears on - * (we assume this never changes) */ -}; - - static void twod_mapping(double fs, double ss, double *px, double *py, struct detgeom_panel *p, double dx, double dy) { @@ -99,8 +87,8 @@ static double r_dev(struct reflpeak *rp) } -static double x_dev(struct reflpeak *rp, struct detgeom *det, - double dx, double dy) +double x_dev(struct reflpeak *rp, struct detgeom *det, + double dx, double dy) { /* Peak position term */ double xpk, ypk, xh, yh; @@ -112,8 +100,8 @@ static double x_dev(struct reflpeak *rp, struct detgeom *det, } -static double y_dev(struct reflpeak *rp, struct detgeom *det, - double dx, double dy) +double y_dev(struct reflpeak *rp, struct detgeom *det, + double dx, double dy) { /* Peak position term */ double xpk, ypk, xh, yh; @@ -573,76 +561,6 @@ static void free_rps_noreflist(struct reflpeak *rps, int n) } -static void write_mille(Mille *mille, int n, UnitCell *cell, - struct reflpeak *rps, struct image *image, - double dx, double dy) -{ -#ifdef HAVE_MILLEPEDE - int i; - float local_gradients[9]; - float global_gradients[6]; - int labels[6]; - - profile_start("mille-calc"); - - /* Spot x-position terms */ - for ( i=0; idetgeom, dx, dy), - 0.65*rps[i].panel->pixel_pitch); - } - - /* Spot y-position terms */ - for ( i=0; idetgeom, dx, dy), - 0.65*rps[i].panel->pixel_pitch); - } - - profile_end("mille-calc"); - -#endif /* HAVE_MILLEPEDE */ -} - - int refine_prediction(struct image *image, Crystal *cr, Mille *mille) { int n; @@ -725,10 +643,15 @@ int refine_prediction(struct image *image, Crystal *cr, Mille *mille) pred_residual(rps, n, image->detgeom, total_x, total_y)); crystal_add_notes(cr, tmp); + #ifdef HAVE_MILLEPEDE if ( mille != NULL ) { + profile_start("mille-calc"); write_mille(mille, n, crystal_get_cell(cr), rps, image, - total_x, total_y); + total_x, total_y, + image->detgeom->top_group, 0, 0); + profile_end("mille-calc"); } + #endif crystal_set_det_shift(cr, total_x, total_y); @@ -741,7 +664,7 @@ int refine_prediction(struct image *image, Crystal *cr, Mille *mille) crystal_set_det_shift(cr, orig_shift_x, orig_shift_y); #ifdef HAVE_MILLEPEDE if ( mille != NULL ) { - mille_delete_last_record(mille); + crystfel_mille_delete_last_record(mille); } #endif return 1; @@ -750,30 +673,10 @@ int refine_prediction(struct image *image, Crystal *cr, Mille *mille) #ifdef HAVE_MILLEPEDE if ( mille != NULL ) { profile_start("mille-write"); - mille_write_record(mille); + crystfel_mille_write_record(mille); profile_end("mille-write"); } #endif return 0; } - - -Mille *crystfel_mille_new(const char *outFileName, - int asBinary, - int writeZero) -{ - #ifdef HAVE_MILLEPEDE - return mille_new(outFileName, asBinary, writeZero); - #else - return NULL; - #endif -} - - -void crystfel_mille_free(Mille *m) -{ - #ifdef HAVE_MILLEPEDE - mille_free(m); - #endif -} diff --git a/libcrystfel/src/predict-refine.h b/libcrystfel/src/predict-refine.h index 91b44e07..5e2552e6 100644 --- a/libcrystfel/src/predict-refine.h +++ b/libcrystfel/src/predict-refine.h @@ -29,24 +29,31 @@ #ifndef PREDICT_REFINE_H #define PREDICT_REFINE_H +struct reflpeak; #include "crystal.h" +#include "crystfel-mille.h" -struct image; - -typedef void *Mille; +struct reflpeak { + Reflection *refl; + struct imagefeature *peak; + double Ih; /* normalised */ + struct detgeom_panel *panel; /* panel the reflection appears on + * (we assume this never changes) */ +}; /** * \file predict-refine.h * Prediction refinement: refinement of indexing solutions before integration. */ -extern Mille *crystfel_mille_new(const char *outFileName, - int asBinary, - int writeZero); -extern void crystfel_mille_free(Mille *m); - extern int refine_prediction(struct image *image, Crystal *cr, Mille *mille); + extern int refine_radius(Crystal *cr, struct image *image); +extern double x_dev(struct reflpeak *rp, struct detgeom *det, + double dx, double dy); + +extern double y_dev(struct reflpeak *rp, struct detgeom *det, + double dx, double dy); #endif /* PREDICT_REFINE_H */ -- cgit v1.2.3