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

Linux Cross Reference
Tina6/tina-libs/tina/geometry/geomCurve_curvature.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/geomCurve_curvature.c,v $
 27  * Date    :  $Date: 2005/02/07 23:27:41 $
 28  * Version :  $Revision: 1.3 $
 29  * CVS Id  :  $Id: geomCurve_curvature.c,v 1.3 2005/02/07 23:27:41 paul Exp $
 30  *
 31  * Author  : Legacy TINA
 32  *
 33  * Notes :
 34  *
 35  *********
 36 */
 37 
 38 #include "geomCurve_curvature.h"
 39 
 40 #if HAVE_CONFIG_H
 41   #include <config.h>
 42 #endif
 43 
 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/geometry/geomDef.h>
 50 
 51 
 52 static double qspline_curvature(Vec2 p0, Vec2 p1, Vec2 p2)
 53 {
 54     double  x0 = vec2_x(p0), x1 = vec2_x(p1), x2 = vec2_x(p2);
 55     double  y0 = vec2_y(p0), y1 = vec2_y(p1), y2 = vec2_y(p2);
 56     double  dx, dy, ddx, ddy, ds;
 57 
 58     dx = 0.5 * (x2 - x0);
 59     ddx = x2 - 2.0 * x1 + x0;
 60     dy = 0.5 * (y2 - y0);
 61     ddy = y2 - 2.0 * y1 + y0;
 62 
 63     ds = sqrt(dx * dx + dy * dy);
 64     return ((dx * ddy - dy * ddx) / (ds * ds * ds));
 65 }
 66 
 67 double  dds_curvature(List * s, int n)
 68 {
 69     int     i;
 70     List *ptr0;
 71     List *ptr2;
 72     Vec2    p0 = {Vec2_id};
 73     Vec2    p1 = {Vec2_id};
 74     Vec2    p2 = {Vec2_id};
 75 
 76     /** try to go n points forward && back **/
 77 
 78     for (ptr0 = s, i = 0; ptr0->last != NULL && i < n; ptr0 = ptr0->last, i++)
 79         ;
 80     for (ptr2 = s, i = 0; ptr2->next != NULL && i < n; ptr2 = ptr2->next, i++)
 81         ;
 82 
 83     DD_GET_POS2(ptr0, p0);
 84     DD_GET_POS2(s, p1);
 85     DD_GET_POS2(ptr2, p2);
 86 
 87     return (qspline_curvature(p0, p1, p2));
 88 }
 89 
 90 double  es_curvature(Tstring * es, List * p, double arcratio, double arcmin, double arcmax)
 91 /* edge type string */
 92 /* point on the string at which to measure curvature */
 93 /* try to achieve this ratio of arc length to radius */
 94 /* min and max arc length to use */
 95 {
 96     List *start = es->start;
 97     List *end = es->end;
 98     List *ptr1;
 99     List *ptr2;
100     List *ptrm;
101     Vec2    p1 = {Vec2_id};
102     Vec2    p2 = {Vec2_id};
103     Vec2    pm = {Vec2_id};
104     double  arc, curvature=0.0;
105     int     extend1, extend2;
106 
107     arc = arcmin;
108     ptr1 = ptr2 = p;
109     DD_GET_POS2(ptr1, p1);
110     DD_GET_POS2(ptr2, p2);
111     extend1 = (ptr1 == start) ? 0 : 1;
112     extend2 = (ptr2 == end) ? 0 : 1;
113     while (extend1 || extend2)
114     {
115         if (extend1)
116         {
117             DD_GET_POS2(ptr1->last, p1);
118             if (vec2_dist(p1, p2) > arc)
119                 extend1 = 0;
120             else
121                 ptr1 = ptr1->last;
122             if (ptr1 == start)
123                 extend1 = 0;
124         }
125         if (extend2)
126         {
127             DD_GET_POS2(ptr2->next, p2);
128             if (vec2_dist(p1, p2) > arc)
129                 extend2 = 0;
130             else
131                 ptr2 = ptr2->next;
132             if (ptr2 == end)
133                 extend2 = 0;
134         }
135         if (!extend1 && !extend2)
136         {
137             double  arcnew;
138 
139             ptrm = ddstr_mid_point(ptr1, ptr2);
140             if (ptr1 == ptrm || ptr2 == ptrm)
141                 return (0.0);
142             DD_GET_POS2(ptrm, pm);
143             curvature = qspline_curvature(p1, pm, p2);
144             arcnew = ((arcnew = arcratio / curvature) > arcmax) ? arcmax : arcnew;
145             if (arcnew > arc)
146             {
147                 arc = arcnew;
148                 if (ptr1 != start)
149                     extend1 = 1;
150                 if (ptr2 != end)
151                     extend2 = 1;
152             }
153         }
154     }
155     return (curvature);
156 }
157 
158 static double qspline_diffgeom(Vec2 p0, Vec2 p1, Vec2 p2, Vec2 * p, Vec2 * t, double *k)
159 {
160     double  x0 = vec2_x(p0), x1 = vec2_x(p1), x2 = vec2_x(p2);
161     double  y0 = vec2_y(p0), y1 = vec2_y(p1), y2 = vec2_y(p2);
162     double  dx, dy, ddx, ddy, ds;
163 
164     dx = 0.5 * (x2 - x0);
165     ddx = x2 - 2.0 * x1 + x0;
166     dy = 0.5 * (y2 - y0);
167     ddy = y2 - 2.0 * y1 + y0;
168 
169     ds = sqrt(dx * dx + dy * dy);
170     *p = p1;
171     *t = vec2(dx / ds, dy / ds);
172     *k = (dx * ddy - dy * ddx) / (ds * ds * ds);
173     return(ds);
174 }
175 
176 void  dds_diffgeom(List * s, Vec2 * p, Vec2 * t, double *k, int n)
177 {
178     int     i;
179     List *ptr0;
180     List *ptr2;
181     Vec2    p0 = {Vec2_id};
182     Vec2    p1 = {Vec2_id};
183     Vec2    p2 = {Vec2_id};
184 
185     /** try to go n points forward && back **/
186 
187     for (ptr0 = s, i = 0; ptr0->last != NULL && i < n; ptr0 = ptr0->last, i++)
188         ;
189     for (ptr2 = s, i = 0; ptr2->next != NULL && i < n; ptr2 = ptr2->next, i++)
190         ;
191     DD_GET_POS2(ptr0, p0);
192     DD_GET_POS2(s, p1);
193     DD_GET_POS2(ptr2, p2);
194     qspline_diffgeom(p0, p1, p2, p, t, k);
195 }
196 
197 void  es_diffgeom(Tstring * es, List * p, double arcratio, double arcmin, double arcmax, Vec2 * pos, Vec2 * t, double *k)
198 /* edge type string */
199 /* point on the string at which to measure curvature */
200 /* try to achieve this ratio of arc length to radius */
201 /* min and max arc length to use */
202 {
203     List *start = es->start;
204     List *end = es->end;
205     List *ptr1;
206     List *ptr2;
207     List *ptrm;
208     Vec2    p1 = {Vec2_id};
209     Vec2    p2 = {Vec2_id};
210     Vec2    pm = {Vec2_id};
211     double  arc, curvature;
212     int     extend1, extend2;
213 
214     arc = arcmin;
215     ptr1 = ptr2 = p;
216     DD_GET_POS2(ptr1, p1);
217     DD_GET_POS2(ptr2, p2);
218     extend1 = (ptr1 == start) ? 0 : 1;
219     extend2 = (ptr2 == end) ? 0 : 1;
220     while (extend1 || extend2)
221     {
222         if (extend1)
223         {
224             DD_GET_POS2(ptr1->last, p1);
225             if (vec2_dist(p1, p2) > arc)
226                 extend1 = 0;
227             else
228                 ptr1 = ptr1->last;
229             if (ptr1 == start)
230                 extend1 = 0;
231         }
232         if (extend2)
233         {
234             DD_GET_POS2(ptr2->next, p2);
235             if (vec2_dist(p1, p2) > arc)
236                 extend2 = 0;
237             else
238                 ptr2 = ptr2->next;
239             if (ptr2 == end)
240                 extend2 = 0;
241         }
242         if (!extend1 && !extend2)
243         {
244             double  arcnew;
245 
246             ptrm = ddstr_mid_point(ptr1, ptr2);
247             DD_GET_POS2(ptrm, pm);
248             qspline_diffgeom(p1, pm, p2, pos, t, k);
249             curvature = *k;
250             arcnew = ((arcnew = arcratio / curvature) > arcmax) ? arcmax : arcnew;
251             if (arcnew > arc)
252             {
253                 arc = arcnew;
254                 if (ptr1 != start)
255                     extend1 = 1;
256                 if (ptr2 != end)
257                     extend2 = 1;
258             }
259         }
260     }
261 }
262 
263 static Vec2 *str_pos2_to_array(Tstring * s, int *n)
264 {
265     int     i, m;
266     Vec2   *va;
267     List *dd;
268 
269     if (s == NULL)
270     {
271         /* BUGFIX JB 22/11/93 was *n == -1; */
272         *n = -1;
273         return (NULL);
274     }
275     m = str_length(s);
276 
277     va = ts_nralloc(m, Vec2);
278     for (dd = s->start, i = 0; i < m; dd = dd->next, i++)
279         DD_GET_POS2(dd, va[i]);
280 
281     *n = m;
282     return (va);
283 }
284 
285 static void str_pos2_from_array(Tstring * s, Vec2 * va)
286 {
287     List *dd;
288     int     i;
289 
290     if (s == NULL || va == NULL)
291         return;
292 
293     for (dd = s->start, i = 0;; dd = dd->next, i++)
294     {
295         DD_SET_POS2(dd, va[i]);
296         if (dd == s->end)
297             break;
298     }
299 }
300 
301 Vec2    vec2_smooth(Vec2 v0, Vec2 v1, Vec2 v2)
302 {
303     double  x0 = vec2_x(v0), y0 = vec2_y(v0);
304     double  x1 = vec2_x(v1), y1 = vec2_y(v1);
305     double  x2 = vec2_x(v2), y2 = vec2_y(v2);
306 
307     x1 = 0.25 * (x0 + 2.0 * x1 + x2);
308     y1 = 0.25 * (y0 + 2.0 * y1 + y2);
309 
310     return (vec2(x1, y1));
311 }
312 
313 void    loop_smooth(int l, Vec2 * array, Vec2 * store)
314 {
315     int     i;
316 
317     store[0] = vec2_smooth(array[l - 1], array[0], array[1]);
318     for (i = 1; i < l - 1; i++)
319         store[i] = vec2_smooth(array[i - 1], array[i], array[i + 1]);
320     store[l - 1] = vec2_smooth(array[l - 2], array[l - 1], array[0]);
321 
322     for (i = 0; i < l; i++)
323         array[i] = store[i];
324 }
325 
326 void    endfix_smooth(int l, Vec2 * array, Vec2 * store)
327 {
328     int     i;
329 
330     for (i = 1; i < l - 1; i++)
331         store[i] = vec2_smooth(array[i - 1], array[i], array[i + 1]);
332 
333     for (i = 1; i < l - 1; i++)
334         array[i] = store[i];
335 }
336 
337 void    str_smooth(Tstring * s, int n)
338 {
339     Vec2   *array;
340     Vec2   *store;
341     int     i, l;
342 
343     array = str_pos2_to_array(s, &l);
344     store = ts_nralloc(l, Vec2);
345 
346     for (i = 0; i < n; i++)
347     {
348         if (s->type == LOOP)
349             loop_smooth(l, array, store);
350         else
351             endfix_smooth(l, array, store);
352     }
353     str_pos2_from_array(s, array);
354 
355     rfree((void *) array);
356     rfree((void *) store);
357 }
358 

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