blob: fb5f6497404a59a98222233dcb4ece5a15952a35 (
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
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
|
/*
* itrans.c
*
* Parameterise features in an image for reconstruction
*
* (c) 2007-2009 Thomas White <taw27@cam.ac.uk>
* (c) 2007 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 an empty feature list */
ImageFeatureList *itrans_peaksearch_none(ImageRecord *imagerecord)
{
return image_feature_list_new();
}
/* Return a feature list for the given image */
ImageFeatureList *itrans_process_image(ImageRecord *image,
PeakSearchMode psmode)
{
ImageFeatureList *list;
switch ( psmode ) {
case PEAKSEARCH_NONE : {
list = itrans_peaksearch_none(image);
break;
}
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;
}
}
|
