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

Linux Cross Reference
Tina6/tina-libs/tina/geometry/geomSpline_kws_run.c

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

  1 /**********
  2  * 
  3  * This file is part of the TINA Open Source Image Analysis Environment
  4  * henceforth known as TINA
  5  *
  6  * TINA is free software; you can redistribute it and/or modify
  7  * it under the terms of the GNU General Public License as 
  8  * published by the Free Software Foundation.
  9  *
 10  * TINA is distributed in the hope that it will be useful,
 11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
 12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 13  * GNU General Public License for more details.
 14  *
 15  * You should have received a copy of the GNU General Public License
 16  * along with TINA; if not, write to the Free Software Foundation, Inc., 
 17  * 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 18  *
 19  * ANY users of TINA who require exemption from the existing licence must
 20  * negotiate a new licence with Dr. Neil.A.Thacker, the sole agent for
 21  * the University of Manchester.
 22  *
 23  **********
 24  * 
 25  * Program :    TINA
 26  * File    :  $Source: /home/tina/cvs/tina-libs/tina/geometry/geomSpline_kws_run.c,v $
 27  * Date    :  $Date: 2005/01/23 19:10:21 $
 28  * Version :  $Revision: 1.3 $
 29  * CVS Id  :  $Id: geomSpline_kws_run.c,v 1.3 2005/01/23 19:10:21 paul Exp $
 30  *
 31  * Notes :
 32  *
 33  *********
 34 */
 35 
 36 
 37 #include "geomSpline_kws_run.h"
 38 
 39 #if HAVE_CONFIG_H
 40   #include <config.h>
 41 #endif
 42 
 43 #include <stdio.h>
 44 #include <math.h>
 45 #include <tina/sys/sysDef.h>
 46 #include <tina/sys/sysPro.h>
 47 #include <tina/math/mathDef.h>
 48 #include <tina/math/mathPro.h>
 49 #include <tina/image/imgDef.h>
 50 #include <tina/image/imgPro.h>
 51 #include <tina/geometry/geom_SplineDef.h>
 52 #include <tina/geometry/geomSpline_kws.h>
 53 
 54 void    kwsnake_correct(Kwsnake * kwsnake, double *fx, double *fy, double dt)
 55 {
 56     int     i;
 57 
 58     for (i = 0; i < kwsnake->n; i++)
 59     {
 60         kwsnake->x[i] -= dt * fx[i];
 61         kwsnake->y[i] -= dt * fy[i];
 62     }
 63 }
 64 
 65 void    kwsnake_correct_maxstep(Kwsnake * kwsnake, double *dx, double *dy, double maxstep)
 66 {
 67     int     i;
 68     double  step = 0.0, k;
 69 
 70     for (i = 0; i < kwsnake->n; i++)
 71     {
 72         step = MAX(step, fabs(dx[i]));
 73         step = MAX(step, fabs(dy[i]));
 74     }
 75     if (step == 0.0)
 76         k = 0.0;
 77     else
 78         k = maxstep / step;
 79 
 80     kwsnake_correct(kwsnake, dx, dy, k);
 81 }
 82 
 83 void    kwsnake_correct_fullstep(Kwsnake * kwsnake, double *dx, double *dy, double fullstep)
 84 {
 85     int     i;
 86 
 87     for (i = 0; i < kwsnake->n; i++)
 88     {
 89         dx[i] = (dx[i] == 0) ? 0 : SGN(dx[i]);
 90         dy[i] = (dy[i] == 0) ? 0 : SGN(dy[i]);
 91     }
 92     kwsnake_correct(kwsnake, dx, dy, fullstep);
 93 }
 94 
 95 static Vec2 grad_orth(float *rast_orth, Vec2 p, Vec2 v, Vec2 q, int n1, int n2)
 96 {
 97     double  d, g;
 98     int     i;
 99 
100     d = vec2_dot(v, vec2_diff(q, p));
101     i = (int)floor(d);
102 
103     if (i < n1 || i > n2)
104         return (vec2(0.0, 0.0));
105 
106     if (i > n1 && i < n2)
107         g = (rast_orth[i + 1] - rast_orth[i - 1]) * 0.5;
108     else if (i == n1)
109         g = rast_orth[i + 1] - rast_orth[i];
110     else
111         g = rast_orth[i] - rast_orth[i - 1];
112     return (vec2_times(g, v));
113 }
114 
115 Vec2    kwsnake_orth(Kwsnake * kwsnake, int i)
116 {
117     int     n = kwsnake->n;
118     double *x = kwsnake->x, *y = kwsnake->y;
119     Vec2    v = {Vec2_id};
120 
121     if (i != 0 && i != n - 1)
122         v = vec2_unit(vec2_diff(vec2(x[i + 1], y[i + 1]), vec2(x[i - 1], y[i - 1
123 ])));
124     else if (i == 0)
125     {
126         if (kwsnake->type == LOOP)
127             v = vec2_unit(vec2_diff(vec2(x[1], y[1]), vec2(x[n - 1], y[n - 1])))
128 ;
129         else
130             v = vec2_unit(vec2_diff(vec2(x[1], y[1]), vec2(x[0], y[0])));
131     } else
132     {
133         if (kwsnake->type == LOOP)
134             v = vec2_unit(vec2_diff(vec2(x[0], y[0]), vec2(x[n - 2], y[n - 2])))
135 ;
136         else
137             v = vec2_unit(vec2_diff(vec2(x[n - 1], y[n - 1]), vec2(x[n - 2], y[n
138  - 2])));
139     }
140 
141     return (vec2(-vec2_y(v), vec2_x(v)));
142 }
143 
144 
145 void    kwsnake_external_step_orth(Kwsnake * kwsnake, float **im_orth, Vec2 * p, Vec2 * v, int n1, int n2, double step, int type)
146 {
147     static int n = 0;                                                   /* static data! */
148     double *x = kwsnake->x;
149     double *y = kwsnake->y;
150     static double *fx = NULL, *fy = NULL;               /* static data! */
151     int     i;
152 
153     if (kwsnake == NULL)
154         return;
155 
156     if (n != kwsnake->n)
157     {
158         n = kwsnake->n;
159         /* BUG  dvector_free(fx, 0); dvector_free(fy, 0); */
160         fx = dvector_alloc(0, n);
161         fy = dvector_alloc(0, n);
162     }
163     for (i = 0; i < n; i++)
164     {
165         Vec2    grad = {Vec2_id};
166 
167         grad = grad_orth(im_orth[i], p[i], v[i], vec2(x[i], y[i]), n1, n2);
168         fx[i] = vec2_x(grad);
169         fy[i] = vec2_y(grad);
170     }
171 
172     switch (type)
173     {
174     case KW_STEP:
175         kwsnake_correct(kwsnake, fx, fy, step);
176         break;
177     case KW_MAXSTEP:
178         kwsnake_correct_maxstep(kwsnake, fx, fy, step);
179         break;
180     case KW_FULLSTEP:
181         kwsnake_correct_fullstep(kwsnake, fx, fy, step);
182         break;
183     }
184 }
185 
186 static void im_gradsubpixf(Imrect * im, double y, double x, double *gy, double *gx)
187 {
188     double  g1, g2;
189 
190     g1 = im_sub_pixf(im, y - 1, x);
191     g2 = im_sub_pixf(im, y + 1, x);
192     *gy = 0.5 * (g2 - g1);
193     g1 = im_sub_pixf(im, y, x - 1);
194     g2 = im_sub_pixf(im, y, x + 1);
195     *gx = 0.5 * (g2 - g1);
196 }
197 
198 static void im_d2subpixf(Imrect * im, double y, double x, double *ixx, double *ixy, double *iyy)
199 {
200     *ixx = 0.5 * (im_sub_pixf(im, y, x + 1)
201                   - 2.0 * im_sub_pixf(im, y, x)
202                   + im_sub_pixf(im, y, x - 1));
203     *iyy = 0.5 * (im_sub_pixf(im, y + 1, x)
204                   - 2.0 * im_sub_pixf(im, y, x)
205                   + im_sub_pixf(im, y - 1, x));
206     *ixy = 0.25 * (im_sub_pixf(im, y + 1, x + 1)
207                    - im_sub_pixf(im, y + 1, x - 1)
208                    - im_sub_pixf(im, y - 1, x + 1)
209                    + im_sub_pixf(im, y - 1, x - 1));
210 }
211 
212 void    kwsnake_external_step_from_pot(Kwsnake * kwsnake, Imrect * pot, double step, int type)
213 {
214     static int n = 0;                                           /* static data! */
215     double *x = kwsnake->x;
216     double *y = kwsnake->y;
217     static double *fx = NULL, *fy = NULL;       /* static data! */
218     int     i;
219 
220     if (kwsnake == NULL)
221         return;
222 
223     if (n != kwsnake->n)
224     {
225         n = kwsnake->n;
226         dvector_free(fx, 0);
227         dvector_free(fy, 0);
228         fx = dvector_alloc(0, n);
229         fy = dvector_alloc(0, n);
230     }
231     for (i = 0; i < n; i++)
232     {
233         double  gx, gy;
234 
235         im_gradsubpixf(pot, y[i], x[i], &gy, &gx);
236         fx[i] = gx;
237         fy[i] = gy;
238     }
239     switch (type)
240     {
241     case KW_STEP:
242         kwsnake_correct(kwsnake, fx, fy, step);
243         break;
244     case KW_MAXSTEP:
245         kwsnake_correct_maxstep(kwsnake, fx, fy, step);
246         break;
247     case KW_FULLSTEP:
248         kwsnake_correct_fullstep(kwsnake, fx, fy, step);
249         break;
250     }
251 }
252 
253 static void pentadiag0(int n_call, double a0, double b0, double c0, double d0, double e0, double *x, int cyclic)
254 {
255     static double *a, *b, *c, *d, *e;           /* static data! */
256     static int n = 0;                                           /* static data! */
257     int     i;
258 
259     if (n != n_call)
260     {
261         n = n_call;
262         dvector_free(a, 0);
263         dvector_free(b, 0);
264         dvector_free(c, 0);
265         dvector_free(d, 0);
266         dvector_free(e, 0);
267         a = dvector_alloc(0, n);
268         b = dvector_alloc(0, n);
269         c = dvector_alloc(0, n);
270         d = dvector_alloc(0, n);
271         e = dvector_alloc(0, n);
272     }
273     for (i = 0; i < n; i++)
274     {
275         a[i] = a0;
276         b[i] = b0;
277         c[i] = c0;
278         d[i] = d0;
279         e[i] = e0;
280     }
281     if (!cyclic)
282     {
283         a[0] = a[1] = e[n - 1] = e[n - 2] = b[0] = d[n - 1] = 0;
284         c[0] += a0 + b0;
285         c[1] += a0;
286         c[n - 1] += e0 + d0;
287         c[n - 2] += e0;
288     }
289     pentadiag(n, a, b, c, d, e, x);
290 }
291 
292 void    kwsnake_internal_step_orth(Kwsnake * kwsnake, Vec2 * p, Vec2 * v, double alpha, double beta)
293 {
294     double *x = kwsnake->x;
295     double *y = kwsnake->y;
296     static double *o;                   /* static data! */
297     double  a, b, c, d, e;
298     static int n = 0;                   /* static data! */
299     int     i;
300 
301     if (n != kwsnake->n)
302     {
303         n = kwsnake->n;
304         dvector_free(o, 0);
305         o = dvector_alloc(0, n);
306     }
307     for (i = 0; i < n; ++i)
308         o[i] = vec2_dot(vec2_diff(vec2(x[i], y[i]), p[i]), v[i]);
309 
310     e = a = beta;
311     d = b = -alpha - 4.0 * beta;
312     c = 1.0 + 2.0 * alpha + 6.0 * beta;
313 
314     pentadiag0(n, a, b, c, d, e, o, kwsnake->type == LOOP);
315 
316     for (i = 0; i < n; ++i)
317     {
318         Vec2    pi = {Vec2_id};
319 
320         pi = vec2_sum(p[i], vec2_times(o[i], v[i]));
321         x[i] = vec2_x(pi);
322         y[i] = vec2_y(pi);
323     }
324 }
325 
326 void    kwsnake_dt_step(Kwsnake * kwsnake, Imrect * pot, double alpha, double beta, double trans_dt, double shape_dt)
327 {
328     static int n = -1;                                                                          /* static data! */
329     static double *fx = NULL, *fy = NULL, fxmean, fymean;       /* static data! */
330     double *x = kwsnake->x;
331     double *y = kwsnake->y;
332     double  a, b, c, d, e, fmax;
333     int     i;
334 
335     if (kwsnake == NULL)
336         return;
337 
338     if (n != kwsnake->n)
339     {
340         n = kwsnake->n;
341         /* BUG  dvector_alloc(fx, 0); dvector_alloc(fy, 0); */
342         fx = dvector_alloc(0, n);
343         fy = dvector_alloc(0, n);
344     }
345     fmax = 0.0;
346     for (i = 0; i < n; i++)
347     {
348         double  gx, gy, k;
349         Vec2    n = {Vec2_id};
350         Vec2    kwsnake_orth();
351 
352         im_gradsubpixf(pot, y[i], x[i], &gy, &gx);
353 
354         n = kwsnake_orth(kwsnake, i);
355         k = fabs(vec2_dot(n, vec2_unit(vec2(gx, gy))));
356         gx *= k;
357         gy *= k;
358 
359         fx[i] = -gx;
360         fxmean += fx[i];
361         fmax = MAX(fmax, fabs(gx));
362         fy[i] = -gy;
363         fymean += fy[i];
364         fmax = MAX(fmax, fabs(gy));
365     }
366     fxmean /= n;
367     fymean /= n;
368 
369     if (fmax == 0.0)
370         fmax = 1.0;
371     for (i = 0; i < n; i++)
372     {
373         fx[i] -= fxmean;
374         fx[i] /= fmax;
375         fy[i] -= fymean;
376         fy[i] /= fmax;
377     }
378 
379     fmax = MAX(fabs(fxmean), fabs(fymean));
380     if (fmax == 0.0)
381         fmax = 1.0;
382     fxmean /= fmax;
383     fymean /= fmax;
384 
385     for (i = 0; i < n; i++)
386     {
387         x[i] += trans_dt * fxmean + shape_dt * fx[i];
388         y[i] += trans_dt * fymean + shape_dt * fy[i];
389     }
390 
391     e = a = shape_dt * beta;
392     d = b = -shape_dt * (alpha + 4.0 * beta);
393     c = 1.0 + shape_dt * (2.0 * alpha + 6.0 * beta);
394 
395     pentadiag0(n, a, b, c, d, e, x, kwsnake->type == LOOP);
396     pentadiag0(n, a, b, c, d, e, y, kwsnake->type == LOOP);
397 }
398 
399 void    kwsnake_dt_step2(Kwsnake * kwsnake, Imrect * pot, double alpha, double beta, double trans_dt, double shape_dt)
400 {
401     static int n = -1;                                                                                  /* static data! */
402     static double *fx = NULL, *fy = NULL, fxmean, fymean;               /* static data! */
403     double *x = kwsnake->x;
404     double *y = kwsnake->y;
405     double  a, b, c, d, e, fmax;
406     int     i;
407 
408     if (kwsnake == NULL)
409         return;
410 
411     if (n != kwsnake->n)
412     {
413         n = kwsnake->n;
414         /* BUG  dvector_alloc(fx, 0); dvector_alloc(fy, 0); */
415         fx = dvector_alloc(0, n);
416         fy = dvector_alloc(0, n);
417     }
418     fmax = 0.0;
419     for (i = 0; i < n; i++)
420     {
421         double  ixx, ixy, iyy, gx, gy;
422         Vec2    normal = {Vec2_id};
423         Vec2    kwsnake_orth();
424 
425         im_d2subpixf((Imrect *) pot, *y, *x, &ixx, &ixy, &iyy);
426         normal = kwsnake_orth(kwsnake, i);
427         gx = ixx * vec2_x(normal) + ixy * vec2_y(normal);
428         gy = ixy * vec2_x(normal) + iyy * vec2_y(normal);
429 
430         fx[i] = -gx;
431         fxmean += fx[i];
432         fmax = MAX(fmax, fabs(gx));
433         fy[i] = -gy;
434         fymean += fy[i];
435         fmax = MAX(fmax, fabs(gy));
436     }
437     fxmean /= n;
438     fymean /= n;
439 
440     if (fmax == 0.0)
441         fmax = 1.0;
442     for (i = 0; i < n; i++)
443     {
444         fx[i] -= fxmean;
445         fx[i] /= fmax;
446         fy[i] -= fymean;
447         fy[i] /= fmax;
448     }
449 
450     fmax = MAX(fabs(fxmean), fabs(fymean));
451     if (fmax == 0.0)
452         fmax = 1.0;
453     fxmean /= fmax;
454     fymean /= fmax;
455 
456     for (i = 0; i < n; i++)
457     {
458         x[i] += trans_dt * fxmean + shape_dt * fx[i];
459         y[i] += trans_dt * fymean + shape_dt * fy[i];
460     }
461 
462     e = a = shape_dt * beta;
463     d = b = -shape_dt * (alpha + 4.0 * beta);
464     c = 1.0 + shape_dt * (2.0 * alpha + 6.0 * beta);
465 
466     pentadiag0(n, a, b, c, d, e, x, kwsnake->type == LOOP);
467     pentadiag0(n, a, b, c, d, e, y, kwsnake->type == LOOP);
468 }
469 
470 double  kwsnake_dist(Kwsnake * kws1, Kwsnake * kws2)
471 {
472     int     n = kws1->n;
473     double *x1 = kws1->x, *x2 = kws2->x;
474     double *y1 = kws1->y, *y2 = kws2->y;
475     double  d, dmax = 0.0;
476     int     i;
477 
478     for (i = 0; i < n; i++)
479     {
480         d = fabs(x1[i] - x2[i]) + fabs(y1[i] - y2[i]);
481         if (dmax < d)
482             dmax = d;
483     }
484     return (dmax);
485 }
486 
487 double  kwsnake_dt_run(Kwsnake * kwsnake, Imrect * pot, double alpha, double beta, double trans_dt, double shape_dt, int steps)
488 {
489     int     i;
490     double  d;
491     Kwsnake *copy = kwsnake_copy(kwsnake);
492 
493     for (i = 0; i < steps; i++)
494         kwsnake_dt_step(kwsnake, pot, alpha, beta, trans_dt, shape_dt);
495     d = kwsnake_dist(kwsnake, copy);
496     kwsnake_free(copy);
497     return (d);
498 }
499 
500 double  kwsnake_dt_run2(Kwsnake * kwsnake, Imrect * pot, double alpha, double beta, double trans_dt, double shape_dt, int steps)
501 {
502     int     i;
503     double  d;
504     Kwsnake *copy = kwsnake_copy(kwsnake);
505 
506     for (i = 0; i < steps; i++)
507         kwsnake_dt_step2(kwsnake, pot, alpha, beta, trans_dt, shape_dt);
508     d = kwsnake_dist(kwsnake, copy);
509     kwsnake_free(copy);
510     return (d);
511 }
512 
513 void    kwsnake_internal_step(Kwsnake * kwsnake, double alpha, double beta)
514 {
515     double *x = kwsnake->x;
516     double *y = kwsnake->y;
517     double  a, b, c, d, e;
518     int     n = kwsnake->n;
519 
520     e = a = beta;
521     d = b = -alpha - 4.0 * beta;
522     c = 1.0 + 2.0 * alpha + 6.0 * beta;
523     pentadiag0(n, a, b, c, d, e, x, kwsnake->type == LOOP);
524     pentadiag0(n, a, b, c, d, e, y, kwsnake->type == LOOP);
525 }
526 
527 void    kwsnake_step_from_pot(Kwsnake * kwsnake, Imrect * pot, double alpha, double beta, double step, int type)
528 {
529     kwsnake_external_step_from_pot(kwsnake, pot, step, type);
530     kwsnake_internal_step(kwsnake, alpha, beta);
531 }
532 
533 void    kwsnake_reorth(Kwsnake * kwsnake, Vec2 * p, Vec2 * v)
534 {
535     double *x, *y;
536     int     i, n;
537 
538     if (kwsnake == NULL)
539         return;
540 
541     x = kwsnake->x;
542     y = kwsnake->y;
543     n = kwsnake->n;
544 
545     for (i = 0; i < n; ++i)
546     {
547         double  dp = vec2_dot(v[i], vec2_diff(vec2(x[i], y[i]), p[i]));
548         Vec2    new = {Vec2_id};
549 
550         new = vec2_sum(p[i], vec2_times(dp, v[i]));
551         x[i] = vec2_x(new);
552         y[i] = vec2_y(new);
553     }
554 }
555 
556 void    kwsnake_orth_step(Kwsnake * kwsnake, float **im_orth, Vec2 * p, Vec2 * v, int n1, int n2, double alpha, double beta, double step, int type)
557 {
558     kwsnake_external_step_orth(kwsnake, im_orth, p, v, n1, n2, step, type);
559     kwsnake_internal_step(kwsnake, alpha, beta);
560 }
561 
562 void    kwsnake_step_orth(Kwsnake * kwsnake, float **im_orth, Vec2 * p, Vec2 * v, int n1, int n2, double alpha, double beta, double step, int type)
563 {
564     kwsnake_external_step_orth(kwsnake, im_orth, p, v, n1, n2, step, type);
565     kwsnake_internal_step(kwsnake, alpha, beta);
566     kwsnake_reorth(kwsnake, p, v);
567 }
568 
569 void    kwsnake_inflate(Kwsnake * kwsnake, Vec2 * p, Vec2 * v, double step)
570 {
571     double *x, *y;
572     int     i, n;
573 
574     if (kwsnake == NULL)
575         return;
576 
577     x = kwsnake->x;
578     y = kwsnake->y;
579     n = kwsnake->n;
580 
581     for (i = 0; i < n; ++i)
582     {
583         Vec2    pi = {Vec2_id};
584 
585         pi = vec2_sum(p[i], vec2_times(step, v[i]));
586         x[i] = vec2_x(pi);
587         y[i] = vec2_y(pi);
588     }
589 }
590 
591 void    kwsnake_grow(Kwsnake * kwsnake, Vec2 * p, Vec2 * v, Imrect * im, double mean, double sd, double step)
592 {
593     double *x, *y;
594     double  gl, gl_low, gl_up;
595     int     i, n;
596 
597     if (kwsnake == NULL)
598         return;
599 
600     x = kwsnake->x;
601     y = kwsnake->y;
602     n = kwsnake->n;
603 
604     gl_low = mean - 2 * sd;
605     gl_up = mean + 2 * sd;
606 
607     for (i = 0; i < n; ++i)
608     {
609         Vec2    pi = {Vec2_id};
610 
611         gl = im_sub_pixf(im, y[i], x[i]);
612         if (gl < gl_low || gl > gl_up)
613         {
614             pi = vec2_sum(p[i], vec2_times(-step, v[i]));
615             x[i] = vec2_x(pi);
616             y[i] = vec2_y(pi);
617             continue;
618         }
619         pi = vec2_sum(p[i], vec2_times(step, v[i]));
620         gl = im_sub_pixf(im, vec2_y(pi), vec2_x(pi));
621         if (gl > gl_low && gl < gl_up)
622         {
623             x[i] = vec2_x(pi);
624             y[i] = vec2_y(pi);
625         }
626     }
627 }
628 
629 /* ARGSUSED quieten lint */
630 void    kwsnake_region(Kwsnake * kwsnake, Imrect * im, double mean, double sd, Vec2 * p, Vec2 * v, double alpha, double beta, double step, int type)
631 
632 
633 
634 
635 
636 /* unused */
637 {
638     kwsnake_grow(kwsnake, p, v, im, mean, sd, step);
639     kwsnake_internal_step(kwsnake, alpha, beta);
640 }
641 

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