~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

Linux Cross Reference
Tina4/src/snake/kw_snake.c

Version: ~
Architecture: ~ [ i386 ] ~ [ alpha ] ~ [ m68k ] ~ [ mips ] ~ [ ppc ] ~ [ sparc ] ~ [ sparc64 ] ~

  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 

~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

This page was automatically generated by the LXR engine.
Visit the LXR main site for more information.