1 /** @(#)Snake handling (create, attract, repel, internal parameters set )
2 */
3 #include <stdio.h>
4 #include <math.h>
5 #include <tina/sys.h>
6 #include <tina/sysfuncs.h>
7 #include <tina/math.h>
8 #include <tina/mathfuncs.h>
9 #include <tina/vision.h>
10 #include <tina/visionfuncs.h>
11 #include <tina/cvr.h>
12 #include <tina/snakefuncs.h>
13
14 /* STATICS */
15 static Kwsnake *kwsnake; /* Current snake */
16 static double alpha = 0.08; /* Snake tension */
17 static double beta = 0.03; /* Snake stiffness */
18 static double sigma = 5.0; /* Snake scale */
19 static Vec2 *p; /* Fixed snake knot points */
20 static Vec2 *v; /* Fixed knot point normal vectors */
21 static float **im_orth; /* Image in normal direction */
22 static float ***im_temp; /* Image template 2*corr_size+1 */
23 static int corr_size = 3; /* Correlation size */
24 static int n1; /* Low im_orth index value */
25 static int n2; /* High im_orth index value */
26 static double step = 0.5;
27 static int step_count = 10;
28 static int kspace = 5;
29 static double maxerror = 20.0; /* Region growing parameters */
30 static double zspace = 8.0; /* Nominal thickness of slice in pixels */
31
32
33
34 /* EXTERNS */
35 extern Tstring *cvr_string_get();
36 extern void contour_changed_notify();
37
38 /* FORWARD REFS */
39 double kw_alpha_get(void);
40 double kw_beta_get(void);
41 double kw_sigma_get(void);
42
43 /* Correlate saved image template against an image about position p in
44 * 1D along direction v over the range n1 to n2 and store the result in
45 * rast */
46 void corr_orth(Imrect * im, float **temp, int corr_size, float *rast, Vec2 p, Vec2 v, int n1, int n2)
47 {
48 int i, j, k;
49 int n = SQR(2 * corr_size + 1);
50
51 for (i = n1; i < n2; ++i)
52 {
53 Vec2 pi = {Vec2_id};
54 double x, y;
55 double m1 = 0, m2 = 0, m12 = 0, s1 = 0, s2 = 0;
56
57 pi = vec2_sum(p, vec2_times((double) i, v));
58 x = vec2_x(pi);
59 y = vec2_y(pi);
60
61 for (j = -corr_size; j <= corr_size; ++j)
62 for (k = -corr_size; k <= corr_size; ++k)
63 {
64 float gl = im_sub_pixf(im, y + j, x + k);
65
66 m1 += temp[j][k];
67 m2 += gl;
68 m12 += gl * temp[j][k];
69 s1 += SQR(temp[j][k]);
70 s2 += SQR(gl);
71 }
72 m1 /= n;
73 m2 /= n;
74 m12 /= n;
75 s1 /= n;
76 s2 /= n;
77 s1 = sqrt(s1 - SQR(m1));
78 s2 = sqrt(s2 - SQR(m2));
79 if (s1 == 0 || s2 == 0)
80 rast[i] = 0;
81 else
82 rast[i] = (m12 - m1 * m2) / s1 / s2;
83 }
84 }
85
86 /* Gross attraction of snake */
87 Vec2 kw_attract(Vec2 mouse_pos)
88 {
89 if (kwsnake)
90 {
91 double *x = kwsnake->x;
92 double *y = kwsnake->y;
93 int i;
94 int n = kwsnake->n;
95
96 for (i = 0; i < n; ++i)
97 {
98 Vec2 v = {Vec2_id};
99 Vec2 pi = {Vec2_id};
100
101 pi = vec2(x[i], y[i]);
102 v = vec2_diff(pi, mouse_pos);
103 pi = vec2_sum(pi, vec2_times(-10.0 / vec2_mod(v), vec2_unit(v)));
104 x[i] = vec2_x(pi);
105 y[i] = vec2_y(pi);
106 }
107 kwsnake_internal_step(kwsnake, kw_alpha_get(), kw_beta_get());
108 contour_changed_notify();
109 } else
110 {
111 format("kw_attract: kwsnake is NULL\n");
112 }
113 return mouse_pos;
114 }
115
116 /* Get snake tension */
117 double kw_alpha_get(void)
118 {
119 return (alpha);
120 }
121
122 /* Set snake tension */
123 void kw_alpha_set(double x)
124 {
125 alpha = x;
126 }
127
128 /* Get snake stiffness */
129 double kw_beta_get(void)
130 {
131 return (beta);
132 }
133
134 /* Set snake stiffness */
135 void kw_beta_set(double x)
136 {
137 beta = x;
138 }
139
140 /* Knot spacing along strings */
141 int kw_kspace_get(void)
142 {
143 return (kspace);
144 }
145
146 void kw_kspace_set(int k)
147 {
148 kspace = k;
149 }
150
151 /* Region finding snake */
152 void kw_region_fn(void)
153 {
154 Imrect *im = imrect_get();
155 Kwsnake *kwsnake = kwsnake_get();
156
157 if (kwsnake && im)
158 {
159 Imrect *rr;
160 Tstring *es;
161 double alpha = kw_alpha_get();
162 double beta = kw_beta_get();
163 double mean, sd = kw_maxerror_get();
164 double step = kw_step_get();
165 int i, count = kw_step_count_get();
166
167 es = es_of_kwsnake(kwsnake);
168 rr = es_region(es);
169 str_free(es, rfree);
170 region_mean_sd(im, rr, &mean, &sd);
171 im_free(rr);
172
173 for (i = 0; i < count; ++i)
174 {
175 kwsnake_pos_and_orth(kwsnake, p, v);
176 kwsnake_region(kwsnake, im, mean, sd, p, v, alpha, beta, step, KW_FULLSTEP);
177 }
178 contour_changed_notify();
179 }
180 }
181
182
183 /* Gross repel of snake */
184 Vec2 kw_repel(Vec2 mouse_pos)
185 {
186 if (kwsnake)
187 {
188 double *x = kwsnake->x;
189 double *y = kwsnake->y;
190 int i;
191 int n = kwsnake->n;
192
193 for (i = 0; i < n; ++i)
194 {
195 Vec2 v = {Vec2_id};
196 Vec2 pi = {Vec2_id};
197
198 pi = vec2(x[i], y[i]);
199 v = vec2_diff(pi, mouse_pos);
200 pi = vec2_sum(pi, vec2_times(10.0 / vec2_mod(v), vec2_unit(v)));
201 x[i] = vec2_x(pi);
202 y[i] = vec2_y(pi);
203 }
204 kwsnake_internal_step(kwsnake, kw_alpha_get(), kw_beta_get());
205 contour_changed_notify();
206 } else
207 {
208 format("kw_repel: kwsnake is NULL\n");
209 }
210 return mouse_pos;
211 }
212
213 /* Get snake scale */
214 double kw_sigma_get(void)
215 {
216 return sigma;
217 }
218
219 /* Set snake scale */
220 void kw_sigma_set(double s)
221 {
222 sigma = s;
223 }
224
225 /* Region growing parameters */
226 double kw_maxerror_get(void)
227 {
228 return (maxerror);
229 }
230
231 void kw_maxerror_set(double newerror)
232 {
233 maxerror = newerror;
234 }
235
236 /* Nominal thickness of slice in pixels */
237 double kw__zspace_get(void)
238 {
239 return (zspace);
240 }
241
242 void kw__zspace_set(double z)
243 {
244 zspace = z;
245 }
246
247 /* Perform 1D correlation between saved snake templates and the
248 * underlying image */
249 void kwsnake_corr_fn(void)
250 {
251 Imrect *im = imrect_get();
252 Prof1 *prof;
253 int i, n;
254
255 if (kwsnake && im_temp)
256 {
257 n = kwsnake->n;
258 prof = prof_gauss_simple(kw_sigma_get(), 0.05);
259
260 kwsnake_pos_and_orth(kwsnake, p, v);
261 for (i = 0; i < n; ++i)
262 {
263 corr_orth(im, im_temp[i], corr_size, im_orth[i], p[i], v[i], n1, n2);
264 raster_envelope_blurr(im_orth[i], prof, n1, n2);
265 }
266 prof1_free(prof);
267 step_orth_fn();
268 }
269 }
270
271 /* Return current snake */
272 Kwsnake *kwsnake_get(void)
273 {
274 return kwsnake;
275 }
276
277 /* Find aprox othogonal direction to ith node in snake */
278 /*
279 Vec2 kwsnake_orth(Kwsnake * kwsnake, int i)
280 {
281 int n = kwsnake->n;
282 double *x = kwsnake->x, *y = kwsnake->y;
283 Vec2 v = {Vec2_id};
284
285 if (i != 0 && i != n - 1)
286 v = vec2_unit(vec2_diff(vec2(x[i + 1], y[i + 1]), vec2(x[i - 1], y[i - 1])));
287 else if (i == 0)
288 {
289 if (kwsnake->type == LOOP)
290 v = vec2_unit(vec2_diff(vec2(x[1], y[1]), vec2(x[n - 1], y[n - 1])));
291 else
292 v = vec2_unit(vec2_diff(vec2(x[1], y[1]), vec2(x[0], y[0])));
293 } else
294 {
295 if (kwsnake->type == LOOP)
296 v = vec2_unit(vec2_diff(vec2(x[0], y[0]), vec2(x[n - 2], y[n - 2])));
297 else
298 v = vec2_unit(vec2_diff(vec2(x[n - 1], y[n - 1]), vec2(x[n - 2], y[n - 2])));
299 }
300
301 return (vec2(-vec2_y(v), vec2_x(v)));
302 }
303 */
304
305 /* Save position and approx orthogonal direction of snake */
306 void kwsnake_pos_and_orth(Kwsnake * kwsnake, Vec2 * p, Vec2 * v)
307 {
308 double *x, *y;
309 int i, n;
310
311 if (kwsnake == NULL)
312 return;
313
314 n = kwsnake->n;
315 x = kwsnake->x;
316 y = kwsnake->y;
317
318 for (i = 0; i < n; ++i)
319 {
320 p[i] = vec2(x[i], y[i]);
321 v[i] = kwsnake_orth(kwsnake, i);
322 }
323 }
324
325 /* Number os steps to take */
326 int kw_step_count_get(void)
327 {
328 return (step_count);
329 }
330
331 void kw_step_count_set(int c)
332 {
333 step_count = c;
334 }
335
336 double kw_step_get(void)
337 {
338 return (step);
339 }
340
341 void kw_step_set(double x)
342 {
343 step = x;
344 }
345
346
347
348 /* Construct a set of square image templates around the nodes of the
349 * snake */
350 void kwsnake_temp_fn(void)
351 {
352 Imrect *im = imrect_get();
353 int i, j, k;
354 int n;
355
356 if (kwsnake == NULL)
357 return;
358
359 n = kwsnake->n;
360
361 im_temp = (float ***) ralloc((unsigned) n * sizeof(void *));
362
363 for (i = 0; i < n; ++i)
364 {
365 float **temp;
366 float x = kwsnake->x[i], y = kwsnake->y[i];
367
368 temp = farray_alloc(-corr_size, -corr_size, corr_size + 1, corr_size + 1);
369 im_temp[i] = temp;
370
371 for (j = -corr_size; j <= corr_size; ++j)
372 for (k = -corr_size; k <= corr_size; ++k)
373 temp[j][k] = im_sub_pixf(im, y + j, x + k);
374 }
375 }
376
377 /* Set up a local edge attracting potential orthogonal to the snake */
378 void local_pot_proc(void)
379 {
380 Imrect *im = imrect_get();
381 Prof1 *prof_blurr;
382 Prof1 *prof_1;
383 int n, i;
384
385 if (kwsnake)
386 {
387 n = kwsnake->n;
388
389 prof_1 = prof_gauss(1.0, 0.05);
390 prof_blurr = prof_gauss_simple(kw_sigma_get(), 0.001); /* wide blurr */
391 kwsnake_pos_and_orth(kwsnake, p, v);
392 orthog_image_data(im, im_orth, p, v, n, n1, n2, 1.0);
393 for (i = 0; i < n; ++i)
394 {
395 raster_pot_proc(im_orth[i], prof_1, n1, n2);
396 raster_envelope_blurr(im_orth[i], prof_blurr, n1, n2);
397 }
398 prof1_free(prof_1);
399 prof1_free(prof_blurr);
400 }
401 }
402
403
404 /* Save image data along orthogonal directions */
405 void orthog_image_data(Imrect * im, float **im_orth, Vec2 * p, Vec2 * v, int n, int n1, int n2, double ds)
406 /* original image */
407 /* saved image storage */
408 /* position and direction vectors */
409 /* count and orthogonal range data */
410 /* image step size */
411 {
412 int i;
413
414 if (im_orth == NULL || im == NULL)
415 return;
416
417 for (i = 0; i < n; ++i)
418 im_get_interp_rast(im_orth[i], im, p[i], v[i], n1, n2, ds);
419 }
420
421 /* Apply an envelope bluring function to a raster */
422 void raster_envelope_blurr(float *r, Prof1 * prof, int n1, int n2)
423 {
424 int i, j, k;
425 int m1 = prof->n1, m2 = prof->n2;
426 float *wm;
427
428 wm = fvector_alloc(n1, n2);
429
430 for (i = n1; i < n2; ++i)
431 {
432 double c, maxc = r[i];
433
434 for (j = m1; j < m2; ++j)
435 {
436 k = i + j;
437 if (k < n1 || k >= n2)
438 continue;
439
440 c = r[k] * prof->el.float_v[j];
441 if (c > maxc)
442 maxc = c;
443 }
444 wm[i] = maxc;
445 }
446
447 for (i = n1; i < n2; ++i)
448 r[i] = -wm[i];
449 }
450
451 /* Compute an edge attracting potential along a raster */
452 void raster_pot_proc(float *r, Prof1 * prof, int n1, int n2)
453 /* the raster */
454 /* a gaussian smoothing profile */
455 /* a range */
456 {
457 float *wm;
458 double pot;
459 int i;
460
461 if (r == NULL || n1 == n2)
462 return;
463
464 smooth_1d_sym(r, n1, n2, prof);
465
466 wm = fvector_alloc(n1, n2);
467
468 pot = (r[n1 + 1] - r[n1]) * 2;
469 wm[n1] = fabs(pot);
470 for (i = n1 + 1; i < n2 - 1; ++i)
471 {
472 pot = r[i + 1] - r[i - 1];
473 wm[i] = fabs(pot);
474 }
475 pot = (r[n2 - 1] - r[n2 - 2]) * 2;
476 wm[n2 - 1] = fabs(pot);
477
478 r[n1] = r[n2 - 1] = 0;
479 for (i = n1 + 1; i < n2 - 1; ++i)
480 {
481 if (wm[i] > 4 && wm[i] > wm[i - 1] && wm[i] >= wm[i + 1]) /* above thres edge */
482 r[i] = wm[i];
483 else
484 r[i] = 0;
485 }
486
487 fvector_free((void *) wm, n1);
488 }
489
490 /* Compute mean and standard deviation of image over a given region */
491 void region_mean_sd(Imrect * im, Imrect * rr, double *mean, double *sd)
492 {
493 int i, j;
494 int fi, li, fj, lj;
495 double gl, x, x2;
496 int n;
497 Imregion *roi;
498
499 *mean = *sd = 0;
500
501 if (im == NULL || rr == NULL)
502 return;
503
504 roi = roi_inter(im->region, rr->region);
505 if (roi == NULL)
506 return;
507
508 fi = roi->ly;
509 fj = roi->lx;
510 li = roi->uy;
511 lj = roi->ux;
512
513 x = x2 = n = 0;
514
515 for (i = fi; i < li; ++i)
516 for (j = fj; j < lj; ++j)
517 {
518 if (IM_UCHAR(rr, i, j))
519 continue;
520
521 gl = IM_UCHAR(im, i, j);
522
523 x += gl;
524 x2 += gl * gl;
525 n++;
526 }
527
528 if (n == 0)
529 return;
530
531 *mean = x / n;
532 *sd = sqrt(x2 / n - SQR(*mean));
533 }
534
535 /* Create snake from edge string */
536 void snake_create(Tstring * edge_string)
537 {
538 int k = kw_kspace_get();
539 int n;
540
541 if (edge_string)
542 {
543 if (kwsnake)
544 {
545 nvector_free((void *) p, 0, sizeof(Vec2));
546 nvector_free((void *) v, 0, sizeof(Vec2));
547 farray_free((char **) im_orth, 0, n1, kwsnake->n, n2);
548 kwsnake_free(kwsnake);
549 }
550 kwsnake = kwsnake_of_es(edge_string, k);
551
552 if (kwsnake)
553 {
554 n = kwsnake->n;
555 n1 = -k;
556 n2 = k + 1;
557 im_orth = farray_alloc(0, n1, n, n2);
558 p = ts_nvector_alloc(0, n, Vec2);
559 v = ts_nvector_alloc(0, n, Vec2);
560 }
561 contour_changed_notify();
562 }
563 }
564
565 /* Perform a series of steps of the snake */
566 void kw_orth_fn(void)
567 {
568 local_pot_proc();
569
570 step_orth_fn();
571 }
572
573 void step_orth_fn(void)
574 {
575 double alpha = kw_alpha_get();
576 double beta = kw_beta_get();
577 double step = kw_step_get();
578 int count = kw_step_count_get();
579 int i;
580
581 if (kwsnake)
582 {
583 local_pot_proc();
584 for (i = 0; i < count; ++i)
585 {
586 kwsnake_step_orth(kwsnake, im_orth, p, v, n1, n2, alpha, beta, step,
587 KW_FULLSTEP);
588 }
589 contour_changed_notify();
590 }
591 }
592
This page was automatically generated by the
LXR engine.
Visit the LXR main site for more
information.