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

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

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