aboutsummaryrefslogtreecommitdiff
path: root/libcrystfel/src/integration.c
diff options
context:
space:
mode:
Diffstat (limited to 'libcrystfel/src/integration.c')
-rw-r--r--libcrystfel/src/integration.c133
1 files changed, 25 insertions, 108 deletions
diff --git a/libcrystfel/src/integration.c b/libcrystfel/src/integration.c
index 4a743c63..16c27bbe 100644
--- a/libcrystfel/src/integration.c
+++ b/libcrystfel/src/integration.c
@@ -54,94 +54,6 @@
#include "integration.h"
-static void check_eigen(gsl_vector *e_val)
-{
- int i;
- double vmax, vmin;
- const int n = e_val->size;
- const double max_condition = 1e6;
- const int verbose = 0;
-
- if ( verbose ) STATUS("Eigenvalues:\n");
- vmin = +INFINITY;
- vmax = 0.0;
- for ( i=0; i<n; i++ ) {
- double val = gsl_vector_get(e_val, i);
- if ( verbose ) STATUS("%i: %e\n", i, val);
- if ( val > vmax ) vmax = val;
- if ( val < vmin ) vmin = val;
- }
-
- for ( i=0; i<n; i++ ) {
- double val = gsl_vector_get(e_val, i);
- if ( val < vmax/max_condition ) {
- gsl_vector_set(e_val, i, 0.0);
- }
- }
-
- vmin = +INFINITY;
- vmax = 0.0;
- for ( i=0; i<n; i++ ) {
- double val = gsl_vector_get(e_val, i);
- if ( val == 0.0 ) continue;
- if ( val > vmax ) vmax = val;
- if ( val < vmin ) vmin = val;
- }
- if ( verbose ) {
- STATUS("Condition number: %e / %e = %5.2f\n",
- vmax, vmin, vmax/vmin);
- }
-}
-
-
-static gsl_vector *solve_svd(gsl_vector *v, gsl_matrix *Mp)
-{
- gsl_matrix *s_vec;
- gsl_vector *s_val;
- int err, n;
- gsl_vector *shifts;
- gsl_matrix *M;
-
- n = v->size;
- if ( v->size != Mp->size1 ) return NULL;
- if ( v->size != Mp->size2 ) return NULL;
-
- M = gsl_matrix_alloc(n, n);
- if ( M == NULL ) return NULL;
- gsl_matrix_memcpy(M, Mp);
-
- s_val = gsl_vector_calloc(n);
- s_vec = gsl_matrix_calloc(n, n);
-
- err = gsl_linalg_SV_decomp_jacobi(M, s_vec, s_val);
- if ( err ) {
- ERROR("SVD failed: %s\n", gsl_strerror(err));
- gsl_matrix_free(s_vec);
- gsl_vector_free(s_val);
- return NULL;
- }
- /* "M" is now "U" */
-
- check_eigen(s_val);
-
- shifts = gsl_vector_calloc(n);
- err = gsl_linalg_SV_solve(M, s_vec, s_val, v, shifts);
- if ( err ) {
- ERROR("Matrix solution failed: %s\n", gsl_strerror(err));
- gsl_matrix_free(s_vec);
- gsl_vector_free(s_val);
- gsl_vector_free(shifts);
- return NULL;
- }
-
- gsl_matrix_free(s_vec);
- gsl_vector_free(s_val);
- gsl_matrix_free(M);
-
- return shifts;
-}
-
-
enum boxmask_val
{
BM_IG, /* "Soft" ignore */
@@ -340,14 +252,14 @@ static void show_reference_profile(struct intcontext *ic, int i)
static void show_peak_box(struct intcontext *ic, struct peak_box *bx,
- int results_pipe)
+ pthread_mutex_t *term_lock)
{
#ifdef HAVE_CURSES_COLOR
int q;
signed int h, k, l;
double fs, ss;
- if ( results_pipe != 0 ) write(results_pipe, "SUSPEND\n", 8);
+ if ( term_lock != NULL ) pthread_mutex_lock(term_lock);
initscr();
clear();
@@ -411,7 +323,10 @@ static void show_peak_box(struct intcontext *ic, struct peak_box *bx,
getch();
endwin();
- if ( results_pipe != 0 ) write(results_pipe, "RELEASE\n", 8);
+ if ( term_lock != NULL ) pthread_mutex_unlock(term_lock);
+#else
+ STATUS("Not showing peak box because CrystFEL was compiled without "
+ "ncurses.\n");
#endif
}
@@ -442,8 +357,8 @@ static void fit_bg(struct intcontext *ic, struct peak_box *bx)
}
}
- /* SVD is massive overkill here */
- ans = solve_svd(v, bx->bgm);
+ /* FIXME: SVD is massive overkill here */
+ ans = solve_svd(v, bx->bgm, NULL, 0);
gsl_vector_free(v);
bx->a = gsl_vector_get(ans, 0);
@@ -500,7 +415,7 @@ static int init_intcontext(struct intcontext *ic)
}
/* How many reference profiles? */
- ic->n_reference_profiles = ic->image->det->n_panels;
+ ic->n_reference_profiles = 1;
ic->reference_profiles = calloc(ic->n_reference_profiles,
sizeof(double *));
if ( ic->reference_profiles == NULL ) return 1;
@@ -1255,7 +1170,7 @@ static int get_int_diag(struct intcontext *ic, Reflection *refl)
static void integrate_prof2d_once(struct intcontext *ic, struct peak_box *bx,
- int results_pipe)
+ pthread_mutex_t *term_lock)
{
bx->intensity = fit_intensity(ic, bx);
bx->sigma = calc_sigma(ic, bx);
@@ -1278,15 +1193,16 @@ static void integrate_prof2d_once(struct intcontext *ic, struct peak_box *bx,
get_detector_pos(bx->refl, &pfs, &pss);
pfs += bx->offs_fs;
pss += bx->offs_ss;
- set_detector_pos(bx->refl, 0.0, pfs, pss);
+ set_detector_pos(bx->refl, pfs, pss);
if ( bx->intensity < -5.0*bx->sigma ) {
ic->n_implausible++;
set_redundancy(bx->refl, 0);
}
- if ( get_int_diag(ic, bx->refl) ) show_peak_box(ic, bx,
- results_pipe);
+ if ( get_int_diag(ic, bx->refl) ) {
+ show_peak_box(ic, bx, term_lock);
+ }
} else {
@@ -1379,7 +1295,7 @@ static void integrate_prof2d(IntegrationMethod meth,
Crystal *cr, struct image *image, IntDiag int_diag,
signed int idh, signed int idk, signed int idl,
double ir_inn, double ir_mid, double ir_out,
- int results_pipe, int **masks)
+ pthread_mutex_t *term_lock, int **masks)
{
RefList *list;
UnitCell *cell;
@@ -1423,7 +1339,7 @@ static void integrate_prof2d(IntegrationMethod meth,
for ( i=0; i<ic.n_boxes; i++ ) {
struct peak_box *bx;
bx = &ic.boxes[i];
- integrate_prof2d_once(&ic, bx, results_pipe);
+ integrate_prof2d_once(&ic, bx, term_lock);
}
//refine_rigid_groups(&ic);
@@ -1436,7 +1352,7 @@ static void integrate_prof2d(IntegrationMethod meth,
static void integrate_rings_once(Reflection *refl, struct image *image,
struct intcontext *ic, UnitCell *cell,
- int results_pipe)
+ pthread_mutex_t *term_lock)
{
double pfs, pss;
struct peak_box *bx;
@@ -1482,6 +1398,7 @@ static void integrate_rings_once(Reflection *refl, struct image *image,
r = check_box(ic, bx, &saturated);
if ( !r ) {
fit_bg(ic, bx);
+ if ( !bg_ok(bx) ) r = 1;
}
bx->offs_fs = 0.0;
bx->offs_ss = 0.0;
@@ -1532,9 +1449,9 @@ static void integrate_rings_once(Reflection *refl, struct image *image,
/* Update position */
pfs += bx->offs_fs;
pss += bx->offs_ss;
- set_detector_pos(refl, 0.0, pfs, pss);
+ set_detector_pos(refl, pfs, pss);
- if ( get_int_diag(ic, refl) ) show_peak_box(ic, bx, results_pipe);
+ if ( get_int_diag(ic, refl) ) show_peak_box(ic, bx, term_lock);
if ( intensity < -5.0*sigma ) {
ic->n_implausible++;
@@ -1634,7 +1551,7 @@ static void integrate_rings(IntegrationMethod meth,
Crystal *cr, struct image *image, IntDiag int_diag,
signed int idh, signed int idk, signed int idl,
double ir_inn, double ir_mid, double ir_out,
- int results_pipe, int **masks)
+ pthread_mutex_t *term_lock, int **masks)
{
RefList *list;
Reflection *refl;
@@ -1670,7 +1587,7 @@ static void integrate_rings(IntegrationMethod meth,
refl != NULL;
refl = next_refl(refl, iter) )
{
- integrate_rings_once(refl, image, &ic, cell, results_pipe);
+ integrate_rings_once(refl, image, &ic, cell, term_lock);
}
//refine_rigid_groups(&ic);
@@ -1687,7 +1604,7 @@ void integrate_all_4(struct image *image, IntegrationMethod meth,
double ir_inn, double ir_mid, double ir_out,
IntDiag int_diag,
signed int idh, signed int idk, signed int idl,
- int results_pipe)
+ pthread_mutex_t *term_lock)
{
int i;
int *masks[image->det->n_panels];
@@ -1727,14 +1644,14 @@ void integrate_all_4(struct image *image, IntegrationMethod meth,
integrate_rings(meth, cr, image,
int_diag, idh, idk, idl,
ir_inn, ir_mid, ir_out,
- results_pipe, masks);
+ term_lock, masks);
break;
case INTEGRATION_PROF2D :
integrate_prof2d(meth, cr, image,
int_diag, idh, idk, idl,
ir_inn, ir_mid, ir_out,
- results_pipe, masks);
+ term_lock, masks);
break;
default :
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-17 00:47:03 +0000