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

Linux Cross Reference
Tina4/src/covira/statsnake.c

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

  1 /*
  2  @(#) Statistical Snakes: Active Region Models.  By Jim Ivins, AIVRU.
  3  */
  4 #include <tina/all_tina.h>
  5 #include <tina/brain.h>
  6 #include <tina/brainfuncs.h>
  7 
  8 /*
  9 Stat snake implemented as Tstring
 10 */
 11 
 12 static Tstring *ss_make(int n)
 13 {
 14     Tstring *ss;
 15     List *start, *end, *ptr;
 16 
 17     if (n < 0)
 18         return (NULL);
 19 
 20     start = dd_list_make(n, VEC2);
 21     end = dd_get_end(start);
 22     ss = str_make(LOOP, start, end);
 23     for (ptr = start;; ptr = ptr->next)
 24     {
 25         ptr->to = vec2_alloc();
 26         if (ptr == end)
 27             break;
 28     }
 29     return (ss);
 30 }
 31 
 32 void ss_free(Tstring * ss)
 33 {
 34     str2_free(ss);
 35 }
 36 
 37 /*
 38 Make JIM snake from JP snake
 39 */
 40 
 41 Tstring *ss_from_snake(Snake * snake)
 42 {
 43     int i;
 44     List *ptr;
 45     Tstring *ss;
 46 
 47     if (snake == NULL)
 48         return (NULL);
 49 
 50     ss = ss_make(snake->n);
 51     for (i = 0, ptr = ss->start;; i++, ptr = ptr->next)
 52     {
 53         DD_VEC2(ptr) = snake->p[i];
 54         if (ptr == ss->end)
 55             break;
 56     }
 57     return (ss);
 58 }
 59 
 60 /*
 61 Make JP snake from JIM snake
 62 */
 63 
 64 Snake *ss_to_snake(Tstring * ss)
 65 {
 66     int i, n;
 67     Snake *snake;
 68     List *ptr;
 69 
 70     if (ss == NULL)
 71         return (NULL);
 72 
 73     /* superstition */
 74     n = ddstr_length(ss->start, ss->end);
 75     if(ss->count != n)
 76     {
 77         fprintf(stderr, "Problem\n");
 78         ss->count = n;
 79     }
 80 
 81     snake = snake_make(ss->count);
 82     for (i = 0, ptr = ss->start;; i++, ptr = ptr->next)
 83     {
 84         snake->p[i] = DD_VEC2(ptr);
 85         if (ptr == ss->end)
 86             break;
 87     }
 88     return (snake);
 89 }
 90 
 91 Tstring *ss_copy(Tstring * ss)
 92 {
 93     return (str2_copy(ss));
 94 }
 95 
 96 double ss_dist(Tstring * ss1, Tstring * ss2)
 97 {
 98     List *ptr1, *ptr2;
 99     double dmax = -MAXFLOAT;
100 
101     if (ss1 == NULL || ss2 == NULL || ss1->count != ss2->count)
102         return (dmax);
103     for (ptr1 = ss1->start, ptr2 = ss2->start;;
104          ptr1 = ptr1->next, ptr2 = ptr2->next)
105     {
106         Vec2 p1, p2;
107         double d;
108 
109         DD_GET_POS2(ptr1, p1);
110         DD_GET_POS2(ptr2, p2);
111         d = vec2_dist(p1, p2);
112         dmax = MAX(dmax, d);
113         if (ptr1 == ss1->end)
114             break;
115     }
116     return (dmax);
117 }
118 
119 /*
120 move statsnake ddlink by step, stopping at roi boundaries
121 */
122 static void ss_el_step(List * ptr, Vec2 step, int lx, int ly, int ux, int uy)
123 {
124 
125     Vec2 p;
126 
127     p = DD_VEC2(ptr);
128     p = vec2_sum(p, step);
129     vec2_x(p) = MAX(lx, vec2_x(p));
130     vec2_x(p) = MIN(ux - 1.0, vec2_x(p));
131     vec2_y(p) = MAX(ly, vec2_y(p));
132     vec2_y(p) = MIN(uy - 1.0, vec2_y(p));
133     DD_VEC2(ptr) = p;
134 }
135 
136 /*
137 increment accumulator position by plot_val (called during line_plot)
138 */
139 
140 static int plot_val;
141 static void cell_plot(int x, int y, Imrect * im)
142 {
143     IM_CHAR(im, y, x) += plot_val;
144 }
145 
146 /*
147 increment line in accumulator by val
148 */
149 
150 static void line_plot(Imrect * accum, Vec2 p1, Vec2 p2, int val)
151 {
152     plot_val = val;
153     apply_bresenham_line(p1, p2, LINE8, cell_plot, accum);
154 }
155 
156 /*
157 increment string segments in accumulator by val
158 */
159 
160 static void dd_string_plot(Imrect * accum, List * ptr1, List * ptr2, int val)
161 {
162     List *ptr;
163 
164     for (ptr = ptr1; ptr != ptr2; ptr = ptr->next)
165     {
166         Vec2 p1, p2;
167 
168         p1 = DD_VEC2(ptr);
169         p2 = DD_VEC2(ptr->next);
170         line_plot(accum, p1, p2, val);
171     }
172 }
173 
174 /*
175 sum accumulator position into total (called during line_total)
176 */
177 
178 static int total;
179 static void cell_test(int x, int y, Imrect * im)
180 {
181     total += IM_CHAR(im, y, x);
182 }
183 
184 /*
185 sum accumulator on line into total
186 */
187 
188 static int line_total(Imrect * accum, Vec2 p1, Vec2 p2)
189 {
190     total = 0;
191     apply_bresenham_line(p1, p2, LINE4, cell_test, accum);
192     return (total);
193 }
194 
195 /*
196 sum accumulator on string segments into total
197 */
198 
199 static int dd_string_total(Imrect * accum, List * ptr1, List * ptr2)
200 {
201     double tot;
202     List *ptr;
203 
204     for (tot = 0.0, ptr = ptr1; ptr != ptr2; ptr = ptr->next)
205     {
206         Vec2 p1, p2;
207 
208         p1 = DD_VEC2(ptr);
209         p2 = DD_VEC2(ptr->next);
210         tot += line_total(accum, p1, p2);
211     }
212     return (tot);
213 }
214 
215 /*
216 Force due to tension: second derivative
217 */
218 
219 static Vec2 grad_ten(List * ptr)
220 {
221     Vec2 p0, p1, p2;
222 
223     p0 = DD_VEC2(ptr->last);
224     p1 = DD_VEC2(ptr);
225     p2 = DD_VEC2(ptr->next);
226     return vec2_sum3(p0, vec2_times(-2.0, p1), p2);
227 }
228 
229 /*
230 Force due to bending: fourth derivative
231 */
232 
233 static Vec2 grad_sti(List * ptr)
234 {
235     Vec2 p0, p1, p2, p3, p4, temp;
236 
237     p0 = DD_VEC2(ptr->last->last);
238     p1 = DD_VEC2(ptr->last);
239     p2 = DD_VEC2(ptr);
240     p3 = DD_VEC2(ptr->next);
241     p4 = DD_VEC2(ptr->next->next);
242     temp = vec2_sum4(p0, vec2_times(-4.0, p1), vec2_times(-4.0, p3), p4);
243     return (vec2_sum(temp, vec2_times(6.0, p2)));
244 }
245 
246 /*
247 rotate vector through 90 degrees
248 */
249 
250 static Vec2 vec2_star(Vec2 v)
251 {
252     return (vec2(-vec2_y(v), vec2_x(v)));
253 }
254 
255 /*
256 smoothed subpixel value
257 */
258 
259 static double smooth_pix(Imrect * im, double y, double x)
260 {
261     double g[3][3], sum = 0;
262     int i, j;
263 
264     for (i = 0; i < 3; i++)
265         for (j = 0; j < 3; j++)
266             g[i][j] = im_sub_pixf(im, y + j - 1, x + i - 1);
267     sum += g[0][0] + 2 * g[0][1] + g[0][2];
268     sum += 2 * g[1][0] + 4 * g[1][1] + 2 * g[1][2];
269     sum += g[2][0] + 2 * g[2][1] + g[2][2];
270     return (sum / 16);
271 }
272 
273 /*
274 binary pressure reverses at pix = mean
275 */
276 
277 static Vec2 grad_pre_bin(Imrect * im, List * ptr, Stats * stats, double k)
278 {
279     double pix;
280     Vec2 p0, p1, p2, n;
281 
282     if (stats == NULL)
283         return (vec2_zero());
284 
285     p0 = DD_VEC2(ptr->last);
286     p1 = DD_VEC2(ptr);
287     p2 = DD_VEC2(ptr->next);
288     n = vec2_times(-0.5, vec2_star(vec2_diff(p2, p0)));
289 
290     pix = smooth_pix(im, vec2_y(p1), vec2_x(p1));
291     if (fabs(stats->mean - pix) > k * stats->sd)
292         return (n);
293     else
294         return (vec2_minus(n));
295 }
296 
297 /*
298 Linear pressure is linear across pix = mean
299 */
300 
301 static Vec2 grad_pre_lin(Imrect * im, List * ptr, Stats * stats, double k)
302 {
303     double pix, range, scale;
304     Vec2 p0, p1, p2, n;
305 
306     if (stats == NULL)
307         return (vec2_zero());
308 
309     range = k * stats->sd;
310     if (range <= 0)
311         return grad_pre_bin(im, ptr, stats, k);
312 
313     p0 = DD_VEC2(ptr->last);
314     p1 = DD_VEC2(ptr);
315     p2 = DD_VEC2(ptr->next);
316     n = vec2_times(0.5, vec2_star(vec2_diff(p2, p0)));
317 
318     pix = smooth_pix(im, vec2_y(p1), vec2_x(p1));
319     scale = range - fabs(stats->mean - pix);
320     return (vec2_times(scale / range, n));
321 }
322 
323 static void stats_pixel(int x, int y, void **data)
324 {
325     Stats *stats = data[0];
326     Imrect *im = data[1];
327     int pix = im_get_pix(im, y, x);
328 
329     ++stats->n;
330     stats->tot1 += pix;
331     stats->tot2 += pix * pix;
332 }
333 
334 static Vec2 grad_pre_ave(Imrect * im, List * ptr, Stats * stats, double k)
335 {
336     int n;
337     double range, scale;
338     Vec2 diff, p0, p1, p2;
339     Stats bound_stats;
340     void *data[2];
341 
342     if (!stats)
343         return vec2_zero();
344 
345     p0 = DD_VEC2(ptr->last);
346     p1 = DD_VEC2(ptr);
347     p2 = DD_VEC2(ptr->next);
348 
349     diff = vec2_diff(p0, p2);
350 
351     bound_stats.n = 0;
352     bound_stats.tot1 = 0.0;
353     bound_stats.tot2 = 0.0;
354     data[0] = &bound_stats;
355     data[1] = im;
356     apply_symmetric_line(p0, p1, stats_pixel, data);
357     apply_symmetric_line(p1, p2, stats_pixel, data);
358     if ((n = bound_stats.n) > 1)
359     {
360         bound_stats.mean = bound_stats.tot1 / n;
361         bound_stats.var = (bound_stats.tot2 - sqr(bound_stats.tot1) / n) / (n - 1);
362         bound_stats.sd = sqrt(bound_stats.var);
363     }
364 
365     if ((range = k * stats->sd) <= 0)
366     {
367         if (fabs(stats->mean - bound_stats.mean) > range)
368             return vec2_times(0.5, vec2(-vec2_y(diff), +vec2_x(diff)));
369         else
370             return vec2_times(0.5, vec2(+vec2_y(diff), -vec2_x(diff)));
371     }
372     else
373     {
374         scale = range - fabs(stats->mean - bound_stats.mean);
375         return vec2_times(scale / range, vec2(+vec2_y(diff), -vec2_x(diff)));
376     }
377 }
378 
379 /*
380 make and initialise snake into accumulator array
381 */
382 
383 Imrect *accum_make(Imrect * im, Tstring * ss)
384 {
385     Imrect *accum = im_alloc(im->height, im->width, NULL, char_v);
386     List *ptr;
387 
388     for (ptr = ss->start;; ptr = ptr->next)
389     {
390         dd_string_plot(accum, ptr, ptr->next, 1);
391         if (ptr == ss->end)
392             break;
393     }
394     return (accum);
395 }
396 
397 /*
398 translational snake step
399 */
400 
401 void stat_snake_trans(Tstring * ss, Imrect * im, Stats * stats,
402   double alpha, double beta, double pressure, double k_s_d, double dt)
403 {
404     int width, height, n;
405     Vec2 g_ten, g_sti, g_pre, gradient, trans;
406     List *ptr;
407 
408     width = im->width;
409     height = im->height;
410     trans = vec2_zero();
411     n = ss->count;
412     alpha *= n * n / (32.0 * 32.0);
413     beta *= n * n * n * n / (32.0 * 32.0 * 32.0 * 32.0);
414     for (ptr = ss->start;; ptr = ptr->next)
415     {
416         g_ten = vec2_times(-alpha, grad_ten(ptr));
417         g_sti = vec2_times(beta, grad_sti(ptr));
418         g_pre = vec2_times(pressure, grad_pre_ave(im, ptr, stats, k_s_d));
419         gradient = vec2_unit(vec2_sum3(g_ten, g_sti, g_pre));
420         trans = vec2_sum(trans, gradient);
421         if (ptr == ss->end)
422             break;
423     }
424     trans = vec2_times(-dt, vec2_unit(trans));
425     for (ptr = ss->start;; ptr = ptr->next)
426     {
427         ss_el_step(ptr, trans, 0, 0, width, height);
428         if (ptr == ss->end)
429             break;
430     }
431 }
432 
433 /*
434 all modes snake step
435 */
436 
437 void stat_snake_step(Tstring * ss, Imrect * accum, Imrect * im, Stats * stats,
438   double alpha, double beta, double pressure, double k_s_d, double dt)
439 {
440     int n, check, width, height;
441     Vec2 g_ten, g_sti, g_pre, gradient;
442     List *ptr;
443 
444     width = im->width;
445     height = im->height;
446     n = ss->count;
447     alpha *= n * n / (32.0 * 32.0);
448     beta *= n * n * n * n / (32.0 * 32.0 * 32.0 * 32.0);
449     dd_string_plot(accum, ss->start->last, ss->start, -1);
450     for (ptr = ss->start;; ptr = ptr->next)
451     {
452         Vec2 p = DD_VEC2(ptr);
453 
454         dd_string_plot(accum, ptr, ptr->next, -1);
455         g_ten = vec2_times(-alpha, grad_ten(ptr));
456         g_sti = vec2_times(beta, grad_sti(ptr));
457         g_pre = vec2_times(pressure, grad_pre_ave(im, ptr, stats, k_s_d));
458         gradient = vec2_unit(vec2_sum3(g_ten, g_sti, g_pre));
459 
460         ss_el_step(ptr, vec2_times(-dt, gradient), 0, 0, width, height);
461         check = dd_string_total(accum, ptr->last, ptr->next);
462         if (check != MOVE_OK)
463         {
464             DD_VEC2(ptr) = p;
465             ss_el_step(ptr, vec2_times(dt, gradient), 0, 0, width, height);
466             check = dd_string_total(accum, ptr->last, ptr->next);
467             if (check != MOVE_OK)
468                 DD_VEC2(ptr) = p;
469         }
470         dd_string_plot(accum, ptr->last, ptr, 1);
471 
472         if (ptr == ss->end)
473             break;
474     }
475     dd_string_plot(accum, ss->start->last, ss->start, 1);
476 }
477 
478 /*
479 free a ddlist link (-> libraries)
480 */
481 
482 static void dd_link_free(List * ptr, void (*freefunc) ())
483 {
484     freefunc(ptr->to, ptr->type);
485     rfree(ptr);
486 }
487 
488 /*
489 delete snake elements to stay above minimum gap
490 */
491 
492 void stat_snake_mingap(Tstring * ss, Imrect * accum, int min_gap)
493 {
494     List *ptr;
495 
496     if (ss->count < 4)
497         return;
498 
499     for (ptr = ss->start;; ptr = ptr->next)
500     {
501         Vec2 p1, p2;
502         double dist;
503         int check;
504 
505         p1 = DD_VEC2(ptr);
506         p2 = DD_VEC2(ptr->next);
507         dist = vec2_dist(p1, p2);
508 
509         if ((dist < min_gap) && (ptr != ss->start) && ss->count > 2)
510         {
511             dd_string_plot(accum, ptr->last, ptr->next, -1);
512 
513             /* try removing node */
514             ptr->last->next = ptr->next;
515             ptr->next->last = ptr->last;
516             check = dd_string_total(accum, ptr->last, ptr->last->next);
517 
518             /* check if ok */
519             if (check == DEL_OK)
520             {
521                 /* if so really remove */
522                 if (ptr == ss->end)
523                     ss->end = ptr->last;
524                 --(ss->count);
525                 dd_link_free(ptr, vec2_free);
526                 ptr = ptr->last;
527                 dd_string_plot(accum, ptr, ptr->next, 1);
528             }
529             else
530             {
531                 /* else replace */
532                 ptr->last->next = ptr;
533                 ptr->next->last = ptr;
534                 dd_string_plot(accum, ptr->last, ptr->next, 1);
535             }
536         }
537         if (ptr == ss->end)
538             break;
539     }
540 }
541 
542 /*
543 add snake elements to stay below maximum gap
544 */
545 
546 void stat_snake_maxgap(Tstring * ss, Imrect * accum, int min_gap, int max_gap)
547 {
548     List *ptr;
549 
550     for (ptr = ss->start;; ptr = ptr->next)
551     {
552         Vec2 p1, p2;
553         double dist;
554 
555         p1 = DD_VEC2(ptr);
556         p2 = DD_VEC2(ptr->next);
557         dist = vec2_dist(p1, p2);
558         if (dist > max_gap && dist > 2 * min_gap)
559         {
560             Vec2 m = vec2_midpoint(p1, p2);
561             List *new;
562 
563             dd_string_plot(accum, ptr, ptr->next, -1);
564             new = dd_link_alloc(vec2_make(m), VEC2);
565             dd_link_addafter(ptr, new);
566             if (ss->end == ptr)
567                 ss->end = new;
568             ++(ss->count);
569             dd_string_plot(accum, ptr, ptr->next->next, 1);
570         }
571         if (ptr == ss->end)
572             break;
573     }
574 }
575 
576 /*
577 combination of above (slightly faster)
578 */
579 
580 void stat_snake_edit(Tstring * ss, Imrect * accum, int min_gap, int max_gap)
581 {
582     List *ptr;
583 
584     for (ptr = ss->start;; ptr = ptr->next)
585     {
586         Vec2 p1, p2;
587         double dist;
588         int check;
589 
590         p1 = DD_VEC2(ptr);
591         p2 = DD_VEC2(ptr->next);
592         dist = vec2_dist(p1, p2);
593 
594         if ((dist < min_gap) && (ptr != ss->start) && ss->count > 2)
595         {
596             dd_string_plot(accum, ptr->last, ptr->next, -1);
597 
598             /* try removing node */
599             ptr->last->next = ptr->next;
600             ptr->next->last = ptr->last;
601             check = dd_string_total(accum, ptr->last, ptr->last->next);
602 
603             /* check if ok */
604             if (check == DEL_OK)
605             {
606                 /* if so really remove */
607                 if (ptr == ss->end)
608                     ss->end = ptr->last;
609                 --(ss->count);
610                 dd_link_free(ptr, vec2_free);
611                 ptr = ptr->last;
612                 dd_string_plot(accum, ptr, ptr->next, 1);
613             }
614             else
615             {
616                 /* else replace */
617                 ptr->last->next = ptr;
618                 ptr->next->last = ptr;
619                 dd_string_plot(accum, ptr->last, ptr->next, 1);
620             }
621         }
622         else if (dist > max_gap)
623         {
624             Vec2 m = vec2_midpoint(p1, p2);
625 
626             dd_string_plot(accum, ptr, ptr->next, -1);
627             check = line_total(accum, p1, m) + line_total(accum, m, p2);
628             if (check == INS_OK)
629             {
630                 List *new = dd_link_alloc(vec2_make(m), VEC2);
631 
632                 dd_link_addafter(ptr, new);
633                 if (ss->end == ptr)
634                     ss->end = new;
635                 ++(ss->count);
636                 dd_string_plot(accum, ptr, ptr->next->next, 1);
637             }
638             else
639                 dd_string_plot(accum, ptr, ptr->next, 1);
640         }
641         if (ptr == ss->end)
642             break;
643     }
644 }
645 
646 Stats *stats_make()
647 {
648     Stats *stats = talloc(Stats);
649 
650     stats->n = 0;
651     stats->tot1 = 0.0;
652     stats->tot2 = 0.0;
653     return (stats);
654 }
655 

~ [ 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.