blob: 1018a66e4f758557887a8f17b51a8491e39d1e1f (
plain)
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.c
*
* Parameterise features in an image for reconstruction
*
* (c) 2007 Thomas White <taw27@cam.ac.uk>
* Gordon Ball <gfb21@cam.ac.uk>
*
* dtr - Diffraction Tomography Reconstruction
*
*/
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
#include <math.h>
#include "image.h"
#include "itrans-threshold.h"
#include "itrans-zaefferer.h"
#include "itrans-stat.h"
/* Return a feature list for the given image */
ImageFeatureList *itrans_process_image(ImageRecord *image, PeakSearchMode psmode) {
ImageFeatureList *list;
switch ( psmode ) {
case PEAKSEARCH_THRESHOLD : { list = itrans_peaksearch_threshold(image); break; }
case PEAKSEARCH_ADAPTIVE_THRESHOLD : { list = itrans_peaksearch_adaptive_threshold(image); break; }
case PEAKSEARCH_ZAEFFERER : { list = itrans_peaksearch_zaefferer(image); break; }
case PEAKSEARCH_STAT : { list = itrans_peaksearch_stat(image); break; }
default: list = NULL;
}
return list;
}
/* Quantification radius */
#define R 20
/* Quantify each feature in the image's feature list */
void itrans_quantify_features(ImageRecord *image) {
int i;
for ( i=0; i<image->features->n_features; i++ ) {
double theta;
double background, total;
int n;
int xi, yi;
background = 0.0;
n = 0;
for ( theta=0; theta<2*M_PI; theta+=2*M_PI/10 ) {
double x, y;
long int xd, yd;
x = image->features->features[i].x + R*sin(theta);
y = image->features->features[i].y + R*cos(theta);
xd = rint(x);
yd = rint(y);
if ( (xd>=0) && (yd>=0) && (xd<image->width) && (yd<image->height) ) {
background += image->image[xd + yd*image->width];
n++;
}
}
background /= n;
total = 0.0;
for ( xi=-20; xi<=20; xi++ ) {
for ( yi=-20; yi<=20; yi++ ) {
double x, y;
long int xd, yd;
if ( xi*xi + yi*yi > 20 ) continue;
x = image->features->features[i].x + xi;
y = image->features->features[i].y + yi;
xd = rint(x);
yd = rint(y);
if ( (xd>=0) && (yd>=0) && (xd<image->width) && (yd<image->height) ) {
total += image->image[xd + yd*image->width] - background;
}
}
}
image->features->features[i].intensity = total;
}
}
|
