1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
|
/*
* itrans-zaefferer.c
*
* Zaefferer peak search
*
* (c) 2007 Thomas White <taw27@cam.ac.uk>
*
* dtr - Diffraction Tomography Reconstruction
*
*/
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
#include <stdint.h>
#include "control.h"
#include "imagedisplay.h"
#include "reflections.h"
#include "utils.h"
#define PEAK_WINDOW_SIZE 20
unsigned int itrans_peaksearch_zaefferer(int16_t *image, ControlContext *ctx, double tilt_degrees, ImageDisplay *imagedisplay) {
int x, y;
unsigned int n_reflections;
n_reflections = 0;
for ( x=1; x<ctx->width-1; x++ ) {
for ( y=1; y<ctx->height-1; y++ ) {
double dx1, dx2, dy1, dy2;
double dxs, dys;
double grad;
if ( ctx->max_d ) {
if ( (x-ctx->x_centre)*(x-ctx->x_centre) + (y-ctx->y_centre)*(y-ctx->y_centre) > ctx->max_d*ctx->max_d ) {
continue;
}
}
/* Get gradients */
dx1 = image[x+ctx->width*y] - image[(x+1)+ctx->width*y];
dx2 = image[(x-1)+ctx->width*y] - image[x+ctx->width*y];
dy1 = image[x+ctx->width*y] - image[(x+1)+ctx->width*(y+1)];
dy2 = image[x+ctx->width*(y-1)] - image[x+ctx->width*y];
/* Average gradient measurements from both sides */
dxs = ((dx1*dx1) + (dx2*dx2)) / 2;
dys = ((dy1*dy1) + (dy2*dy2)) / 2;
/* Calculate overall gradient */
grad = dxs + dys;
if ( grad > 400 ) {
unsigned int mask_x, mask_y;
unsigned int sx, sy;
double max;
unsigned int did_something = 1;
mask_x = x;
mask_y = y;
while ( (did_something) && (distance(mask_x, mask_y, x, y)<50) ) {
max = image[mask_x+ctx->width*mask_y];
did_something = 0;
for ( sy=biggest(mask_y-PEAK_WINDOW_SIZE/2, 0); sy<smallest(mask_y+PEAK_WINDOW_SIZE/2, ctx->height); sy++ ) {
for ( sx=biggest(mask_x-PEAK_WINDOW_SIZE/2, 0); sx<smallest(mask_x+PEAK_WINDOW_SIZE/2, ctx->width); sx++ ) {
if ( image[sx+ctx->width*sy] > max ) {
max = image[sx+ctx->width*sy];
mask_x = sx;
mask_y = sy;
did_something = 1;
}
}
}
}
if ( !did_something ) {
if ( ctx->fmode == FORMULATION_PIXELSIZE ) {
reflection_add_from_reciprocal(ctx, (signed)(mask_x-ctx->x_centre)*ctx->pixel_size,
(signed)(mask_y-ctx->y_centre)*ctx->pixel_size,
tilt_degrees, image[mask_x + ctx->width*mask_y]);
} else {
reflection_add_from_dp(ctx, (signed)(mask_x-ctx->x_centre), (signed)(mask_y-ctx->y_centre),
tilt_degrees, image[mask_x + ctx->width*mask_y]);
}
if ( ctx->first_image ) {
imagedisplay_mark_point(imagedisplay, mask_x, mask_y);
}
n_reflections++;
}
}
}
}
return n_reflections;
}
|
