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
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
|
#include <cstdlib>
#include <cmath>
#include <cstring>
#include <stack>
#include <vector>
#include <iostream>
const float MYINFINITY = 1.0e7f;
inline float SQR(float x) { return x * x; }
inline int SQR(int x) { return x * x; }
inline float INTERP(float a, float b, float t) { return a * (1 - t) + b * t; }
struct Coord {
int i, j;
Coord(int i_, int j_)
: i(i_)
, j(j_){};
Coord() {}
int operator==(const Coord &c) const { return i == c.i && j == c.j; }
};
using namespace std;
const float ALIVE = -10, NARROW_BAND = ALIVE + 1,
FAR_AWAY = ALIVE + 2; // Classification values for pixels in an
// image. The values are arbitrary.
int skelft2DSize(int nx, int ny) {
int dx = (int) floor(
pow(2.0f, int(log(float(nx)) / log(2.0f) +
1))); // Find minimal pow of 2 which fits the input image
int dy = (int) floor(pow(2.0f, int(log(float(ny)) / log(2.0f) + 1)));
int fboSize = max(dx, dy); // DT/FT/skeleton image size (pow of 2, should be
// able to contain the input image)
return fboSize;
}
// Classify 'f' into inside, outside, and boundary pixels according to isovalue
// 'iso'
// Inside pixels are those with value > iso (if thr_upper==true), else those
// with value < iso.
template <class T>
void classify(vector<Coord> &s, const T *f, int xm, int ym, int xM, int yM,
int size, T iso, bool thr_upper, float *siteParam,
bool scan_bot_to_top) {
for (int i = 0; i < size * size; ++i)
siteParam[i] = ALIVE; // First: all points are ALIVE, i.e., sites
for (int j = ym; j < yM; ++j)
for (int i = xm; i < xM; ++i) {
T fptr = f[i + j * size]; // Mark subset of points which are
// FAR_AWAY, i.e. NOT sites
if (((thr_upper && fptr >= iso) || (!thr_upper && fptr <= iso)))
siteParam[i + size * j] = FAR_AWAY;
}
int js, je, jd;
if (scan_bot_to_top) {
js = ym;
je = yM;
jd = 1;
} else {
js = yM - 1;
je = ym - 1;
jd = -1;
}
for (int j = js; j != je; j += jd)
for (int i = xm; i < xM; ++i) {
float &cv = siteParam[i + size * j];
if (cv == FAR_AWAY)
if (((i - 1 + size * j) > 0 &&
siteParam[i - 1 + size * j] == ALIVE) ||
siteParam[i + 1 + size * j] == ALIVE ||
((i + size * (j - 1)) > 0 &&
siteParam[i + size * (j - 1)] == ALIVE) ||
siteParam[i + size * (j + 1)] == ALIVE) {
cv = NARROW_BAND;
s.push_back(Coord(i, j));
}
}
}
// Encode input image (classified as inside/outside/boundary into a
// boundary-length parameterization image, needed for skeleton computations.
//
float encodeBoundary(vector<Coord> &s, int xm, int ym, int xM, int yM,
float *siteParam, int size) {
const float LARGE_BOUNDARY_VAL =
1000; // Value that the boundary-param will jump with between disjoint,
// self-connected, boundary pixel chains
float length = -LARGE_BOUNDARY_VAL + 1; // Compute the boundary length
while (s.size()) // Process all narrowband points, add their
// arc-length-parameterization in siteParam
{
Coord pt = s[s.size() - 1];
s.pop_back(); // Pick narrowband-point, test if already numbered
if (siteParam[pt.i + size * pt.j] > 0)
continue;
int ci, cj;
float cc = 0; // Point pt not numbered
for (int dj = -1; dj < 2; ++dj) // Find min-numbered-neighbour of pt
for (int di = -1; di < 2; ++di) {
int ii = pt.i + di, jj = pt.j + dj;
if (ii < xm || ii >= xM || jj < ym || jj >= yM)
continue;
float val = siteParam[ii + size * jj];
if (val < 0)
continue;
if (cc && val > cc)
continue;
ci = ii;
cj = jj;
cc = val;
}
float c = (!cc) ? length + LARGE_BOUNDARY_VAL
: cc + sqrt(float((pt.i - ci) * (pt.i - ci) +
(pt.j - cj) * (pt.j - cj)));
siteParam[pt.i + size * pt.j] =
c; // Also store the sites' parameterization in siteParam[] for CUDA
if (c > length)
length = c; // Determine length of boundary, used for
// wraparound-distance-computations
for (int dj = -1; dj < 2; ++dj)
for (int di = -1; di < 2; ++di) {
if (di == 0 && dj == 0)
continue;
int ii = pt.i + di, jj = pt.j + dj;
if (ii >= xm && ii < xM && jj >= ym && jj < yM &&
siteParam[ii + size * jj] == NARROW_BAND)
s.push_back(Coord(ii, jj));
}
}
for (int i = 0, iend = size * size; i < iend;
++i) // All ALIVE points are marked now as non-sites:
if (siteParam[i] < 0)
siteParam[i] = 0;
return length;
}
float skelft2DMakeBoundary(const float *f, int xm, int ym, int xM, int yM,
float *siteParam, int size, float iso,
bool thr_upper) {
bool scan_bot_to_top = false;
vector<Coord> s;
s.reserve((xM - xm) * (yM - ym));
classify(s, f, xm, ym, xM, yM, size, iso, thr_upper, siteParam,
scan_bot_to_top);
return encodeBoundary(s, xm, ym, xM, yM, siteParam, size);
}
float skelft2DMakeBoundary(const unsigned char *f, int xm, int ym, int xM,
int yM, float *siteParam, int size, short iso,
bool thr_upper) {
bool scan_bot_to_top = false;
vector<Coord> s;
s.reserve((xM - xm) * (yM - ym));
classify(s, f, xm, ym, xM, yM, size, (unsigned char) iso, thr_upper,
siteParam, scan_bot_to_top);
return encodeBoundary(s, xm, ym, xM, yM, siteParam, size);
}
void skelft2DSave(short *outputFT, int dx, int dy, const char *f) {
FILE *fp = fopen(f, "wb");
fprintf(fp, "P5 %d %d 255\n", dx, dy);
const int SIZE = 3000;
unsigned char buf[SIZE];
int size = dx * dy;
int bb = 0;
float range = max(dx, dy) / 255.0f;
for (short *v = outputFT, *vend = outputFT + size; v < vend; ++v) {
short val = *v;
buf[bb++] = (unsigned char) (val / range);
if (bb == SIZE) {
fwrite(buf, sizeof(unsigned char), SIZE, fp);
bb = 0;
}
}
if (bb)
fwrite(buf, sizeof(unsigned char), bb, fp);
fclose(fp);
}
|
