Branch data Line data Source code
1 : : /*=========================================================================
2 : :
3 : : Module: $RCSfile: V_QuadMetric.cpp,v $
4 : :
5 : : Copyright (c) 2006 Sandia Corporation.
6 : : All rights reserved.
7 : : See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
8 : :
9 : : This software is distributed WITHOUT ANY WARRANTY; without even
10 : : the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
11 : : PURPOSE. See the above copyright notice for more information.
12 : :
13 : : =========================================================================*/
14 : :
15 : : /*
16 : : *
17 : : * QuadMetric.cpp contains quality calculations for Quads
18 : : *
19 : : * This file is part of VERDICT
20 : : *
21 : : */
22 : :
23 : : #define VERDICT_EXPORTS
24 : :
25 : : #include "moab/verdict.h"
26 : : #include "VerdictVector.hpp"
27 : : #include "verdict_defines.hpp"
28 : : #include "V_GaussIntegration.hpp"
29 : : #include <memory.h>
30 : :
31 : : //! the average area of a quad
32 : : double verdict_quad_size = 0;
33 : :
34 : : /*!
35 : : weights based on the average size of a quad
36 : : */
37 : 35 : int get_weight( double& m11, double& m21, double& m12, double& m22 )
38 : : {
39 : :
40 : 35 : m11 = 1;
41 : 35 : m21 = 0;
42 : 35 : m12 = 0;
43 : 35 : m22 = 1;
44 : :
45 : 35 : double scale = sqrt( verdict_quad_size / ( m11 * m22 - m21 * m12 ) );
46 : :
47 : 35 : m11 *= scale;
48 : 35 : m21 *= scale;
49 : 35 : m12 *= scale;
50 : 35 : m22 *= scale;
51 : :
52 : 35 : return 1;
53 : : }
54 : :
55 : : //! return the average area of a quad
56 : 36 : C_FUNC_DEF void v_set_quad_size( double size )
57 : : {
58 : 36 : verdict_quad_size = size;
59 : 36 : }
60 : :
61 : : //! returns whether the quad is collapsed or not
62 : 80 : VerdictBoolean is_collapsed_quad( double coordinates[][3] )
63 : : {
64 [ + - ][ + + ]: 80 : if( coordinates[3][0] == coordinates[2][0] && coordinates[3][1] == coordinates[2][1] &&
[ - + ]
65 : 36 : coordinates[3][2] == coordinates[2][2] )
66 : 0 : return VERDICT_TRUE;
67 : :
68 : : else
69 : 80 : return VERDICT_FALSE;
70 : : }
71 : :
72 : 276 : void make_quad_edges( VerdictVector edges[4], double coordinates[][3] )
73 : : {
74 : :
75 : 552 : edges[0].set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
76 : 276 : coordinates[1][2] - coordinates[0][2] );
77 : 552 : edges[1].set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
78 : 276 : coordinates[2][2] - coordinates[1][2] );
79 : 552 : edges[2].set( coordinates[3][0] - coordinates[2][0], coordinates[3][1] - coordinates[2][1],
80 : 276 : coordinates[3][2] - coordinates[2][2] );
81 : 552 : edges[3].set( coordinates[0][0] - coordinates[3][0], coordinates[0][1] - coordinates[3][1],
82 : 276 : coordinates[0][2] - coordinates[3][2] );
83 : 276 : }
84 : :
85 : 124 : void signed_corner_areas( double areas[4], double coordinates[][3] )
86 : : {
87 [ + - ][ + + ]: 620 : VerdictVector edges[4];
88 [ + - ]: 124 : make_quad_edges( edges, coordinates );
89 : :
90 [ + - ][ + + ]: 620 : VerdictVector corner_normals[4];
91 [ + - ][ + - ]: 124 : corner_normals[0] = edges[3] * edges[0];
92 [ + - ][ + - ]: 124 : corner_normals[1] = edges[0] * edges[1];
93 [ + - ][ + - ]: 124 : corner_normals[2] = edges[1] * edges[2];
94 [ + - ][ + - ]: 124 : corner_normals[3] = edges[2] * edges[3];
95 : :
96 : : // principal axes
97 [ + - ][ + + ]: 372 : VerdictVector principal_axes[2];
98 [ + - ][ + - ]: 124 : principal_axes[0] = edges[0] - edges[2];
99 [ + - ][ + - ]: 124 : principal_axes[1] = edges[1] - edges[3];
100 : :
101 : : // quad center unit normal
102 [ + - ]: 124 : VerdictVector unit_center_normal;
103 [ + - ][ + - ]: 124 : unit_center_normal = principal_axes[0] * principal_axes[1];
104 [ + - ]: 124 : unit_center_normal.normalize();
105 : :
106 [ + - ]: 124 : areas[0] = unit_center_normal % corner_normals[0];
107 [ + - ]: 124 : areas[1] = unit_center_normal % corner_normals[1];
108 [ + - ]: 124 : areas[2] = unit_center_normal % corner_normals[2];
109 [ + - ]: 124 : areas[3] = unit_center_normal % corner_normals[3];
110 : 124 : }
111 : :
112 : : /*!
113 : : localize the coordinates of a quad
114 : :
115 : : localizing puts the centriod of the quad
116 : : at the orgin and also rotates the quad
117 : : such that edge (0,1) is aligned with the x axis
118 : : and the quad normal lines up with the y axis.
119 : :
120 : : */
121 : 0 : void localize_quad_coordinates( VerdictVector nodes[4] )
122 : : {
123 : : int i;
124 [ # # ][ # # ]: 0 : VerdictVector global[4] = { nodes[0], nodes[1], nodes[2], nodes[3] };
[ # # ][ # # ]
125 : :
126 [ # # ][ # # ]: 0 : VerdictVector center = ( global[0] + global[1] + global[2] + global[3] ) / 4.0;
[ # # ][ # # ]
127 [ # # ]: 0 : for( i = 0; i < 4; i++ )
128 [ # # ]: 0 : global[i] -= center;
129 : :
130 [ # # ]: 0 : VerdictVector vector_diff;
131 [ # # ]: 0 : VerdictVector vector_sum;
132 [ # # ]: 0 : VerdictVector ref_point( 0.0, 0.0, 0.0 );
133 [ # # ][ # # ]: 0 : VerdictVector tmp_vector, normal( 0.0, 0.0, 0.0 );
134 [ # # ][ # # ]: 0 : VerdictVector vector1, vector2;
135 : :
136 [ # # ]: 0 : for( i = 0; i < 4; i++ )
137 : : {
138 [ # # ]: 0 : vector1 = global[i];
139 [ # # ]: 0 : vector2 = global[( i + 1 ) % 4];
140 [ # # ][ # # ]: 0 : vector_diff = vector2 - vector1;
141 [ # # ]: 0 : ref_point += vector1;
142 [ # # ][ # # ]: 0 : vector_sum = vector1 + vector2;
143 : :
144 [ # # ][ # # ]: 0 : tmp_vector.set( vector_diff.y() * vector_sum.z(), vector_diff.z() * vector_sum.x(),
[ # # ][ # # ]
145 [ # # ][ # # ]: 0 : vector_diff.x() * vector_sum.y() );
[ # # ]
146 [ # # ]: 0 : normal += tmp_vector;
147 : : }
148 : :
149 [ # # ]: 0 : normal.normalize();
150 [ # # ]: 0 : normal *= -1.0;
151 : :
152 [ # # ]: 0 : VerdictVector local_x_axis = global[1] - global[0];
153 [ # # ]: 0 : local_x_axis.normalize();
154 : :
155 [ # # ]: 0 : VerdictVector local_y_axis = normal * local_x_axis;
156 [ # # ]: 0 : local_y_axis.normalize();
157 : :
158 [ # # ]: 0 : for( i = 0; i < 4; i++ )
159 : : {
160 [ # # ][ # # ]: 0 : nodes[i].x( global[i] % local_x_axis );
161 [ # # ][ # # ]: 0 : nodes[i].y( global[i] % local_y_axis );
162 [ # # ][ # # ]: 0 : nodes[i].z( global[i] % normal );
163 : : }
164 : 0 : }
165 : :
166 : : /*!
167 : : moves and rotates the quad such that it enables us to
168 : : use components of ef's
169 : : */
170 : 0 : void localize_quad_for_ef( VerdictVector node_pos[4] )
171 : : {
172 : :
173 [ # # ]: 0 : VerdictVector centroid( node_pos[0] );
174 [ # # ]: 0 : centroid += node_pos[1];
175 [ # # ]: 0 : centroid += node_pos[2];
176 [ # # ]: 0 : centroid += node_pos[3];
177 : :
178 [ # # ]: 0 : centroid /= 4.0;
179 : :
180 [ # # ]: 0 : node_pos[0] -= centroid;
181 [ # # ]: 0 : node_pos[1] -= centroid;
182 [ # # ]: 0 : node_pos[2] -= centroid;
183 [ # # ]: 0 : node_pos[3] -= centroid;
184 : :
185 [ # # ][ # # ]: 0 : VerdictVector rotate = node_pos[1] + node_pos[2] - node_pos[3] - node_pos[0];
[ # # ]
186 [ # # ]: 0 : rotate.normalize();
187 : :
188 [ # # ]: 0 : double cosine = rotate.x();
189 [ # # ]: 0 : double sine = rotate.y();
190 : :
191 : : double xnew;
192 : :
193 [ # # ]: 0 : for( int i = 0; i < 4; i++ )
194 : : {
195 [ # # ][ # # ]: 0 : xnew = cosine * node_pos[i].x() + sine * node_pos[i].y();
196 [ # # ][ # # ]: 0 : node_pos[i].y( -sine * node_pos[i].x() + cosine * node_pos[i].y() );
[ # # ]
197 [ # # ]: 0 : node_pos[i].x( xnew );
198 : : }
199 : 0 : }
200 : :
201 : : /*!
202 : : returns the normal vector of a quad
203 : : */
204 : 17 : VerdictVector quad_normal( double coordinates[][3] )
205 : : {
206 : : // get normal at node 0
207 [ + - ][ + - ]: 17 : VerdictVector edge0, edge1;
208 : :
209 : 34 : edge0.set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
210 [ + - ]: 17 : coordinates[1][2] - coordinates[0][2] );
211 : :
212 : 34 : edge1.set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
213 [ + - ]: 17 : coordinates[3][2] - coordinates[0][2] );
214 : :
215 [ + - ]: 17 : VerdictVector norm0 = edge0 * edge1;
216 [ + - ]: 17 : norm0.normalize();
217 : :
218 : : // because some faces may have obtuse angles, check another normal at
219 : : // node 2 for consistent sense
220 : :
221 : 34 : edge0.set( coordinates[2][0] - coordinates[3][0], coordinates[2][1] - coordinates[3][1],
222 [ + - ]: 17 : coordinates[2][2] - coordinates[3][2] );
223 : :
224 : 34 : edge1.set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
225 [ + - ]: 17 : coordinates[2][2] - coordinates[1][2] );
226 : :
227 [ + - ]: 17 : VerdictVector norm2 = edge0 * edge1;
228 [ + - ]: 17 : norm2.normalize();
229 : :
230 : : // if these two agree, we are done, else test a third to decide
231 : :
232 [ + - ][ + - ]: 17 : if( ( norm0 % norm2 ) > 0.0 )
233 : : {
234 [ + - ]: 17 : norm0 += norm2;
235 [ + - ]: 17 : norm0 *= 0.5;
236 [ + - ]: 17 : return norm0;
237 : : }
238 : :
239 : : // test normal at node1
240 : :
241 : 0 : edge0.set( coordinates[1][0] - coordinates[2][0], coordinates[1][1] - coordinates[2][1],
242 [ # # ]: 0 : coordinates[1][2] - coordinates[2][2] );
243 : :
244 : 0 : edge1.set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
245 [ # # ]: 0 : coordinates[1][2] - coordinates[0][2] );
246 : :
247 [ # # ]: 0 : VerdictVector norm1 = edge0 * edge1;
248 [ # # ]: 0 : norm1.normalize();
249 : :
250 [ # # ][ # # ]: 0 : if( ( norm0 % norm1 ) > 0.0 )
251 : : {
252 [ # # ]: 0 : norm0 += norm1;
253 [ # # ]: 0 : norm0 *= 0.5;
254 [ # # ]: 0 : return norm0;
255 : : }
256 : : else
257 : : {
258 [ # # ]: 0 : norm2 += norm1;
259 [ # # ]: 0 : norm2 *= 0.5;
260 [ # # ]: 17 : return norm2;
261 : : }
262 : : }
263 : :
264 : : /*!
265 : : the edge ratio of a quad
266 : :
267 : : NB (P. Pebay 01/19/07):
268 : : Hmax / Hmin where Hmax and Hmin are respectively the maximum and the
269 : : minimum edge lengths
270 : : */
271 : 16 : C_FUNC_DEF double v_quad_edge_ratio( int /*num_nodes*/, double coordinates[][3] )
272 : : {
273 [ + - ][ + + ]: 80 : VerdictVector edges[4];
274 [ + - ]: 16 : make_quad_edges( edges, coordinates );
275 : :
276 [ + - ]: 16 : double a2 = edges[0].length_squared();
277 [ + - ]: 16 : double b2 = edges[1].length_squared();
278 [ + - ]: 16 : double c2 = edges[2].length_squared();
279 [ + - ]: 16 : double d2 = edges[3].length_squared();
280 : :
281 : : double mab, Mab, mcd, Mcd, m2, M2;
282 [ - + ]: 16 : if( a2 < b2 )
283 : : {
284 : 0 : mab = a2;
285 : 0 : Mab = b2;
286 : : }
287 : : else // b2 <= a2
288 : : {
289 : 16 : mab = b2;
290 : 16 : Mab = a2;
291 : : }
292 [ - + ]: 16 : if( c2 < d2 )
293 : : {
294 : 0 : mcd = c2;
295 : 0 : Mcd = d2;
296 : : }
297 : : else // d2 <= c2
298 : : {
299 : 16 : mcd = d2;
300 : 16 : Mcd = c2;
301 : : }
302 [ - + ]: 16 : m2 = mab < mcd ? mab : mcd;
303 [ - + ]: 16 : M2 = Mab > Mcd ? Mab : Mcd;
304 : :
305 [ - + ]: 16 : if( m2 < VERDICT_DBL_MIN )
306 : 0 : return (double)VERDICT_DBL_MAX;
307 : : else
308 : : {
309 : 16 : double edge_ratio = sqrt( M2 / m2 );
310 : :
311 [ + - ][ + - ]: 16 : if( edge_ratio > 0 ) return (double)VERDICT_MIN( edge_ratio, VERDICT_DBL_MAX );
312 [ # # ]: 16 : return (double)VERDICT_MAX( edge_ratio, -VERDICT_DBL_MAX );
313 : : }
314 : : }
315 : :
316 : : /*!
317 : : maximum of edge ratio of a quad
318 : :
319 : : maximum edge length ratio at quad center
320 : : */
321 : 8 : C_FUNC_DEF double v_quad_max_edge_ratio( int /*num_nodes*/, double coordinates[][3] )
322 : : {
323 [ + - ][ + + ]: 40 : VerdictVector quad_nodes[4];
324 [ + - ]: 8 : quad_nodes[0].set( coordinates[0][0], coordinates[0][1], coordinates[0][2] );
325 [ + - ]: 8 : quad_nodes[1].set( coordinates[1][0], coordinates[1][1], coordinates[1][2] );
326 [ + - ]: 8 : quad_nodes[2].set( coordinates[2][0], coordinates[2][1], coordinates[2][2] );
327 [ + - ]: 8 : quad_nodes[3].set( coordinates[3][0], coordinates[3][1], coordinates[3][2] );
328 : :
329 [ + - ][ + + ]: 24 : VerdictVector principal_axes[2];
330 [ + - ][ + - ]: 8 : principal_axes[0] = quad_nodes[1] + quad_nodes[2] - quad_nodes[0] - quad_nodes[3];
[ + - ][ + - ]
331 [ + - ][ + - ]: 8 : principal_axes[1] = quad_nodes[2] + quad_nodes[3] - quad_nodes[0] - quad_nodes[1];
[ + - ][ + - ]
332 : :
333 [ + - ]: 8 : double len1 = principal_axes[0].length();
334 [ + - ]: 8 : double len2 = principal_axes[1].length();
335 : :
336 [ + - ][ - + ]: 8 : if( len1 < VERDICT_DBL_MIN || len2 < VERDICT_DBL_MIN ) return (double)VERDICT_DBL_MAX;
337 : :
338 [ - + ]: 8 : double max_edge_ratio = VERDICT_MAX( len1 / len2, len2 / len1 );
339 : :
340 [ + - ][ + - ]: 8 : if( max_edge_ratio > 0 ) return (double)VERDICT_MIN( max_edge_ratio, VERDICT_DBL_MAX );
341 [ # # ]: 8 : return (double)VERDICT_MAX( max_edge_ratio, -VERDICT_DBL_MAX );
342 : : }
343 : :
344 : : /*!
345 : : the aspect ratio of a quad
346 : :
347 : : NB (P. Pebay 01/20/07):
348 : : this is a generalization of the triangle aspect ratio
349 : : using Heron's formula.
350 : : */
351 : 17 : C_FUNC_DEF double v_quad_aspect_ratio( int /*num_nodes*/, double coordinates[][3] )
352 : : {
353 : :
354 [ + - ][ + + ]: 85 : VerdictVector edges[4];
355 [ + - ]: 17 : make_quad_edges( edges, coordinates );
356 : :
357 [ + - ]: 17 : double a1 = edges[0].length();
358 [ + - ]: 17 : double b1 = edges[1].length();
359 [ + - ]: 17 : double c1 = edges[2].length();
360 [ + - ]: 17 : double d1 = edges[3].length();
361 : :
362 [ + + ]: 17 : double ma = a1 > b1 ? a1 : b1;
363 [ - + ]: 17 : double mb = c1 > d1 ? c1 : d1;
364 [ + + ]: 17 : double hm = ma > mb ? ma : mb;
365 : :
366 [ + - ]: 17 : VerdictVector ab = edges[0] * edges[1];
367 [ + - ]: 17 : VerdictVector cd = edges[2] * edges[3];
368 [ + - ][ + - ]: 17 : double denominator = ab.length() + cd.length();
369 : :
370 [ - + ]: 17 : if( denominator < VERDICT_DBL_MIN ) return (double)VERDICT_DBL_MAX;
371 : :
372 : 17 : double aspect_ratio = .5 * hm * ( a1 + b1 + c1 + d1 ) / denominator;
373 : :
374 [ + - ][ + - ]: 17 : if( aspect_ratio > 0 ) return (double)VERDICT_MIN( aspect_ratio, VERDICT_DBL_MAX );
375 [ # # ]: 17 : return (double)VERDICT_MAX( aspect_ratio, -VERDICT_DBL_MAX );
376 : : }
377 : :
378 : : /*!
379 : : the radius ratio of a quad
380 : :
381 : : NB (P. Pebay 01/19/07):
382 : : this metric is called "radius ratio" by extension of a concept that does
383 : : not exist in general with quads -- although a different name should probably
384 : : be used in the future.
385 : : */
386 : 16 : C_FUNC_DEF double v_quad_radius_ratio( int /*num_nodes*/, double coordinates[][3] )
387 : : {
388 : : static const double normal_coeff = 1. / ( 2. * sqrt( 2. ) );
389 : :
390 [ + - ][ + + ]: 80 : VerdictVector edges[4];
391 [ + - ]: 16 : make_quad_edges( edges, coordinates );
392 : :
393 [ + - ]: 16 : double a2 = edges[0].length_squared();
394 [ + - ]: 16 : double b2 = edges[1].length_squared();
395 [ + - ]: 16 : double c2 = edges[2].length_squared();
396 [ + - ]: 16 : double d2 = edges[3].length_squared();
397 : :
398 [ + - ]: 16 : VerdictVector diag;
399 : 32 : diag.set( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
400 [ + - ]: 16 : coordinates[2][2] - coordinates[0][2] );
401 [ + - ]: 16 : double m2 = diag.length_squared();
402 : :
403 : 32 : diag.set( coordinates[3][0] - coordinates[1][0], coordinates[3][1] - coordinates[1][1],
404 [ + - ]: 16 : coordinates[3][2] - coordinates[1][2] );
405 [ + - ]: 16 : double n2 = diag.length_squared();
406 : :
407 [ - + ]: 16 : double t0 = a2 > b2 ? a2 : b2;
408 [ - + ]: 16 : double t1 = c2 > d2 ? c2 : d2;
409 [ - + ]: 16 : double t2 = m2 > n2 ? m2 : n2;
410 [ - + ]: 16 : double h2 = t0 > t1 ? t0 : t1;
411 [ - + ]: 16 : h2 = h2 > t2 ? h2 : t2;
412 : :
413 [ + - ]: 16 : VerdictVector ab = edges[0] * edges[1];
414 [ + - ]: 16 : VerdictVector bc = edges[1] * edges[2];
415 [ + - ]: 16 : VerdictVector cd = edges[2] * edges[3];
416 [ + - ]: 16 : VerdictVector da = edges[3] * edges[0];
417 : :
418 [ + - ]: 16 : t0 = da.length();
419 [ + - ]: 16 : t1 = ab.length();
420 [ + - ]: 16 : t2 = bc.length();
421 [ + - ]: 16 : double t3 = cd.length();
422 : :
423 [ - + ]: 16 : t0 = t0 < t1 ? t0 : t1;
424 [ - + ]: 16 : t2 = t2 < t3 ? t2 : t3;
425 [ - + ]: 16 : t0 = t0 < t2 ? t0 : t2;
426 : :
427 [ - + ]: 16 : if( t0 < VERDICT_DBL_MIN ) return (double)VERDICT_DBL_MAX;
428 : :
429 : 16 : double radius_ratio = normal_coeff * sqrt( ( a2 + b2 + c2 + d2 ) * h2 ) / t0;
430 : :
431 [ + - ][ + - ]: 16 : if( radius_ratio > 0 ) return (double)VERDICT_MIN( radius_ratio, VERDICT_DBL_MAX );
432 [ # # ]: 16 : return (double)VERDICT_MAX( radius_ratio, -VERDICT_DBL_MAX );
433 : : }
434 : :
435 : : /*!
436 : : the average Frobenius aspect of a quad
437 : :
438 : : NB (P. Pebay 01/20/07):
439 : : this metric is calculated by averaging the 4 Frobenius aspects at
440 : : each corner of the quad, when the reference triangle is right isosceles.
441 : : */
442 : 16 : C_FUNC_DEF double v_quad_med_aspect_frobenius( int /*num_nodes*/, double coordinates[][3] )
443 : : {
444 : :
445 [ + - ][ + + ]: 80 : VerdictVector edges[4];
446 [ + - ]: 16 : make_quad_edges( edges, coordinates );
447 : :
448 [ + - ]: 16 : double a2 = edges[0].length_squared();
449 [ + - ]: 16 : double b2 = edges[1].length_squared();
450 [ + - ]: 16 : double c2 = edges[2].length_squared();
451 [ + - ]: 16 : double d2 = edges[3].length_squared();
452 : :
453 [ + - ]: 16 : VerdictVector ab = edges[0] * edges[1];
454 [ + - ]: 16 : VerdictVector bc = edges[1] * edges[2];
455 [ + - ]: 16 : VerdictVector cd = edges[2] * edges[3];
456 [ + - ]: 16 : VerdictVector da = edges[3] * edges[0];
457 : :
458 [ + - ]: 16 : double ab1 = ab.length();
459 [ + - ]: 16 : double bc1 = bc.length();
460 [ + - ]: 16 : double cd1 = cd.length();
461 [ + - ]: 16 : double da1 = da.length();
462 : :
463 [ + - ][ + - ]: 16 : if( ab1 < VERDICT_DBL_MIN || bc1 < VERDICT_DBL_MIN || cd1 < VERDICT_DBL_MIN || da1 < VERDICT_DBL_MIN )
[ + - ][ - + ]
464 : 0 : return (double)VERDICT_DBL_MAX;
465 : :
466 : 16 : double qsum = ( a2 + b2 ) / ab1;
467 : 16 : qsum += ( b2 + c2 ) / bc1;
468 : 16 : qsum += ( c2 + d2 ) / cd1;
469 : 16 : qsum += ( d2 + a2 ) / da1;
470 : :
471 : 16 : double med_aspect_frobenius = .125 * qsum;
472 : :
473 [ + - ][ + - ]: 16 : if( med_aspect_frobenius > 0 ) return (double)VERDICT_MIN( med_aspect_frobenius, VERDICT_DBL_MAX );
474 [ # # ]: 16 : return (double)VERDICT_MAX( med_aspect_frobenius, -VERDICT_DBL_MAX );
475 : : }
476 : :
477 : : /*!
478 : : the maximum Frobenius aspect of a quad
479 : :
480 : : NB (P. Pebay 01/20/07):
481 : : this metric is calculated by taking the maximum of the 4 Frobenius aspects at
482 : : each corner of the quad, when the reference triangle is right isosceles.
483 : : */
484 : 16 : C_FUNC_DEF double v_quad_max_aspect_frobenius( int /*num_nodes*/, double coordinates[][3] )
485 : : {
486 : :
487 [ + - ][ + + ]: 80 : VerdictVector edges[4];
488 [ + - ]: 16 : make_quad_edges( edges, coordinates );
489 : :
490 [ + - ]: 16 : double a2 = edges[0].length_squared();
491 [ + - ]: 16 : double b2 = edges[1].length_squared();
492 [ + - ]: 16 : double c2 = edges[2].length_squared();
493 [ + - ]: 16 : double d2 = edges[3].length_squared();
494 : :
495 [ + - ]: 16 : VerdictVector ab = edges[0] * edges[1];
496 [ + - ]: 16 : VerdictVector bc = edges[1] * edges[2];
497 [ + - ]: 16 : VerdictVector cd = edges[2] * edges[3];
498 [ + - ]: 16 : VerdictVector da = edges[3] * edges[0];
499 : :
500 [ + - ]: 16 : double ab1 = ab.length();
501 [ + - ]: 16 : double bc1 = bc.length();
502 [ + - ]: 16 : double cd1 = cd.length();
503 [ + - ]: 16 : double da1 = da.length();
504 : :
505 [ + - ][ + - ]: 16 : if( ab1 < VERDICT_DBL_MIN || bc1 < VERDICT_DBL_MIN || cd1 < VERDICT_DBL_MIN || da1 < VERDICT_DBL_MIN )
[ + - ][ - + ]
506 : 0 : return (double)VERDICT_DBL_MAX;
507 : :
508 : 16 : double qmax = ( a2 + b2 ) / ab1;
509 : :
510 : 16 : double qcur = ( b2 + c2 ) / bc1;
511 [ - + ]: 16 : qmax = qmax > qcur ? qmax : qcur;
512 : :
513 : 16 : qcur = ( c2 + d2 ) / cd1;
514 [ - + ]: 16 : qmax = qmax > qcur ? qmax : qcur;
515 : :
516 : 16 : qcur = ( d2 + a2 ) / da1;
517 [ - + ]: 16 : qmax = qmax > qcur ? qmax : qcur;
518 : :
519 : 16 : double max_aspect_frobenius = .5 * qmax;
520 : :
521 [ + - ][ + - ]: 16 : if( max_aspect_frobenius > 0 ) return (double)VERDICT_MIN( max_aspect_frobenius, VERDICT_DBL_MAX );
522 [ # # ]: 16 : return (double)VERDICT_MAX( max_aspect_frobenius, -VERDICT_DBL_MAX );
523 : : }
524 : :
525 : : /*!
526 : : skew of a quad
527 : :
528 : : maximum ||cos A|| where A is the angle between edges at quad center
529 : : */
530 : 9 : C_FUNC_DEF double v_quad_skew( int /*num_nodes*/, double coordinates[][3] )
531 : : {
532 [ + - ][ + + ]: 45 : VerdictVector node_pos[4];
533 [ + + ]: 45 : for( int i = 0; i < 4; i++ )
534 [ + - ]: 36 : node_pos[i].set( coordinates[i][0], coordinates[i][1], coordinates[i][2] );
535 : :
536 [ + - ][ + + ]: 27 : VerdictVector principle_axes[2];
537 [ + - ][ + - ]: 9 : principle_axes[0] = node_pos[1] + node_pos[2] - node_pos[3] - node_pos[0];
[ + - ][ + - ]
538 [ + - ][ + - ]: 9 : principle_axes[1] = node_pos[2] + node_pos[3] - node_pos[0] - node_pos[1];
[ + - ][ + - ]
539 : :
540 [ + - ][ - + ]: 9 : if( principle_axes[0].normalize() < VERDICT_DBL_MIN ) return 0.0;
541 [ + - ][ - + ]: 9 : if( principle_axes[1].normalize() < VERDICT_DBL_MIN ) return 0.0;
542 : :
543 [ + - ]: 9 : double skew = fabs( principle_axes[0] % principle_axes[1] );
544 : :
545 [ + - ]: 9 : return (double)VERDICT_MIN( skew, VERDICT_DBL_MAX );
546 : : }
547 : :
548 : : /*!
549 : : taper of a quad
550 : :
551 : : maximum ratio of lengths derived from opposite edges
552 : : */
553 : 9 : C_FUNC_DEF double v_quad_taper( int /*num_nodes*/, double coordinates[][3] )
554 : : {
555 [ + - ][ + + ]: 45 : VerdictVector node_pos[4];
556 [ + + ]: 45 : for( int i = 0; i < 4; i++ )
557 [ + - ]: 36 : node_pos[i].set( coordinates[i][0], coordinates[i][1], coordinates[i][2] );
558 : :
559 [ + - ][ + + ]: 27 : VerdictVector principle_axes[2];
560 [ + - ][ + - ]: 9 : principle_axes[0] = node_pos[1] + node_pos[2] - node_pos[3] - node_pos[0];
[ + - ][ + - ]
561 [ + - ][ + - ]: 9 : principle_axes[1] = node_pos[2] + node_pos[3] - node_pos[0] - node_pos[1];
[ + - ][ + - ]
562 : :
563 [ + - ][ + - ]: 9 : VerdictVector cross_derivative = node_pos[0] + node_pos[2] - node_pos[1] - node_pos[3];
[ + - ]
564 : :
565 : : double lengths[2];
566 [ + - ]: 9 : lengths[0] = principle_axes[0].length();
567 [ + - ]: 9 : lengths[1] = principle_axes[1].length();
568 : :
569 : : // get min length
570 [ + + ]: 9 : lengths[0] = VERDICT_MIN( lengths[0], lengths[1] );
571 : :
572 [ - + ]: 9 : if( lengths[0] < VERDICT_DBL_MIN ) return VERDICT_DBL_MAX;
573 : :
574 [ + - ]: 9 : double taper = cross_derivative.length() / lengths[0];
575 [ + - ]: 9 : return (double)VERDICT_MIN( taper, VERDICT_DBL_MAX );
576 : : }
577 : :
578 : : /*!
579 : : warpage of a quad
580 : :
581 : : deviation of element from planarity
582 : : */
583 : 9 : C_FUNC_DEF double v_quad_warpage( int /*num_nodes*/, double coordinates[][3] )
584 : : {
585 : :
586 [ + - ][ + + ]: 45 : VerdictVector edges[4];
587 [ + - ]: 9 : make_quad_edges( edges, coordinates );
588 : :
589 [ + - ][ + + ]: 45 : VerdictVector corner_normals[4];
590 [ + - ][ + - ]: 9 : corner_normals[0] = edges[3] * edges[0];
591 [ + - ][ + - ]: 9 : corner_normals[1] = edges[0] * edges[1];
592 [ + - ][ + - ]: 9 : corner_normals[2] = edges[1] * edges[2];
593 [ + - ][ + - ]: 9 : corner_normals[3] = edges[2] * edges[3];
594 : :
595 [ + - ][ + - ]: 36 : if( corner_normals[0].normalize() < VERDICT_DBL_MIN || corner_normals[1].normalize() < VERDICT_DBL_MIN ||
[ + - ][ + - ]
[ - + ]
596 [ + - ][ + - ]: 27 : corner_normals[2].normalize() < VERDICT_DBL_MIN || corner_normals[3].normalize() < VERDICT_DBL_MIN )
[ + - ][ - + ]
597 : 0 : return (double)VERDICT_DBL_MIN;
598 : :
599 : : double warpage =
600 [ + - ][ + - ]: 9 : pow( VERDICT_MIN( corner_normals[0] % corner_normals[2], corner_normals[1] % corner_normals[3] ), 3 );
[ + + ][ + - ]
[ + - ]
601 : :
602 [ + - ][ + - ]: 9 : if( warpage > 0 ) return (double)VERDICT_MIN( warpage, VERDICT_DBL_MAX );
603 [ # # ]: 9 : return (double)VERDICT_MAX( warpage, -VERDICT_DBL_MAX );
604 : : }
605 : :
606 : : /*!
607 : : the area of a quad
608 : :
609 : : jacobian at quad center
610 : : */
611 : 44 : C_FUNC_DEF double v_quad_area( int /*num_nodes*/, double coordinates[][3] )
612 : : {
613 : :
614 : : double corner_areas[4];
615 [ + - ]: 44 : signed_corner_areas( corner_areas, coordinates );
616 : :
617 : 44 : double area = 0.25 * ( corner_areas[0] + corner_areas[1] + corner_areas[2] + corner_areas[3] );
618 : :
619 [ + - ][ + - ]: 44 : if( area > 0 ) return (double)VERDICT_MIN( area, VERDICT_DBL_MAX );
620 [ # # ]: 44 : return (double)VERDICT_MAX( area, -VERDICT_DBL_MAX );
621 : : }
622 : :
623 : : /*!
624 : : the stretch of a quad
625 : :
626 : : sqrt(2) * minimum edge length / maximum diagonal length
627 : : */
628 : 9 : C_FUNC_DEF double v_quad_stretch( int /*num_nodes*/, double coordinates[][3] )
629 : : {
630 [ + - ][ + + ]: 45 : VerdictVector edges[4], temp;
[ + - ]
631 [ + - ]: 9 : make_quad_edges( edges, coordinates );
632 : :
633 : : double lengths_squared[4];
634 [ + - ]: 9 : lengths_squared[0] = edges[0].length_squared();
635 [ + - ]: 9 : lengths_squared[1] = edges[1].length_squared();
636 [ + - ]: 9 : lengths_squared[2] = edges[2].length_squared();
637 [ + - ]: 9 : lengths_squared[3] = edges[3].length_squared();
638 : :
639 : 18 : temp.set( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
640 [ + - ]: 9 : coordinates[2][2] - coordinates[0][2] );
641 [ + - ]: 9 : double diag02 = temp.length_squared();
642 : :
643 : 18 : temp.set( coordinates[3][0] - coordinates[1][0], coordinates[3][1] - coordinates[1][1],
644 [ + - ]: 9 : coordinates[3][2] - coordinates[1][2] );
645 [ + - ]: 9 : double diag13 = temp.length_squared();
646 : :
647 : : static const double QUAD_STRETCH_FACTOR = sqrt( 2.0 );
648 : :
649 : : // 'diag02' is now the max diagonal of the quad
650 [ - + ]: 9 : diag02 = VERDICT_MAX( diag02, diag13 );
651 : :
652 [ - + ]: 9 : if( diag02 < VERDICT_DBL_MIN )
653 : 0 : return (double)VERDICT_DBL_MAX;
654 : : else
655 : : {
656 : : double stretch =
657 [ - + ][ + + ]: 9 : (double)( QUAD_STRETCH_FACTOR * sqrt( VERDICT_MIN( VERDICT_MIN( lengths_squared[0], lengths_squared[1] ),
658 : : VERDICT_MIN( lengths_squared[2], lengths_squared[3] ) ) /
659 [ - + ][ # # ]: 9 : diag02 ) );
[ + + ]
660 : :
661 [ + - ]: 9 : return (double)VERDICT_MIN( stretch, VERDICT_DBL_MAX );
662 : : }
663 : : }
664 : :
665 : : /*!
666 : : the largest angle of a quad
667 : :
668 : : largest included quad area (degrees)
669 : : */
670 : 9 : C_FUNC_DEF double v_quad_maximum_angle( int /*num_nodes*/, double coordinates[][3] )
671 : : {
672 : :
673 : : // if this is a collapsed quad, just pass it on to
674 : : // the tri_largest_angle routine
675 [ + - ][ - + ]: 9 : if( is_collapsed_quad( coordinates ) == VERDICT_TRUE ) return v_tri_maximum_angle( 3, coordinates );
[ # # ]
676 : :
677 : : double angle;
678 : 9 : double max_angle = 0.0;
679 : :
680 [ + - ][ + + ]: 45 : VerdictVector edges[4];
681 : 18 : edges[0].set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
682 [ + - ]: 9 : coordinates[1][2] - coordinates[0][2] );
683 : 18 : edges[1].set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
684 [ + - ]: 9 : coordinates[2][2] - coordinates[1][2] );
685 : 18 : edges[2].set( coordinates[3][0] - coordinates[2][0], coordinates[3][1] - coordinates[2][1],
686 [ + - ]: 9 : coordinates[3][2] - coordinates[2][2] );
687 : 18 : edges[3].set( coordinates[0][0] - coordinates[3][0], coordinates[0][1] - coordinates[3][1],
688 [ + - ]: 9 : coordinates[0][2] - coordinates[3][2] );
689 : :
690 : : // go around each node and calculate the angle
691 : : // at each node
692 : : double length[4];
693 [ + - ]: 9 : length[0] = edges[0].length();
694 [ + - ]: 9 : length[1] = edges[1].length();
695 [ + - ]: 9 : length[2] = edges[2].length();
696 [ + - ]: 9 : length[3] = edges[3].length();
697 : :
698 [ + - ][ + - ]: 9 : if( length[0] <= VERDICT_DBL_MIN || length[1] <= VERDICT_DBL_MIN || length[2] <= VERDICT_DBL_MIN ||
[ + - ][ - + ]
699 : 9 : length[3] <= VERDICT_DBL_MIN )
700 : 0 : return 0.0;
701 : :
702 [ + - ]: 9 : angle = acos( -( edges[0] % edges[1] ) / ( length[0] * length[1] ) );
703 [ + - ]: 9 : max_angle = VERDICT_MAX( angle, max_angle );
704 : :
705 [ + - ]: 9 : angle = acos( -( edges[1] % edges[2] ) / ( length[1] * length[2] ) );
706 [ + + ]: 9 : max_angle = VERDICT_MAX( angle, max_angle );
707 : :
708 [ + - ]: 9 : angle = acos( -( edges[2] % edges[3] ) / ( length[2] * length[3] ) );
709 [ - + ]: 9 : max_angle = VERDICT_MAX( angle, max_angle );
710 : :
711 [ + - ]: 9 : angle = acos( -( edges[3] % edges[0] ) / ( length[3] * length[0] ) );
712 [ - + ]: 9 : max_angle = VERDICT_MAX( angle, max_angle );
713 : :
714 : 9 : max_angle = max_angle * 180.0 / VERDICT_PI;
715 : :
716 : : // if any signed areas are < 0, then you are getting the wrong angle
717 : : double areas[4];
718 [ + - ]: 9 : signed_corner_areas( areas, coordinates );
719 : :
720 [ + - ][ + - ]: 9 : if( areas[0] < 0 || areas[1] < 0 || areas[2] < 0 || areas[3] < 0 ) { max_angle = 360 - max_angle; }
[ + - ][ - + ]
721 : :
722 [ + - ][ + - ]: 9 : if( max_angle > 0 ) return (double)VERDICT_MIN( max_angle, VERDICT_DBL_MAX );
723 [ # # ]: 9 : return (double)VERDICT_MAX( max_angle, -VERDICT_DBL_MAX );
724 : : }
725 : :
726 : : /*!
727 : : the smallest angle of a quad
728 : :
729 : : smallest included quad angle (degrees)
730 : : */
731 : 1 : C_FUNC_DEF double v_quad_minimum_angle( int /*num_nodes*/, double coordinates[][3] )
732 : : {
733 : : // if this quad is a collapsed quad, then just
734 : : // send it to the tri_smallest_angle routine
735 [ + - ][ - + ]: 1 : if( is_collapsed_quad( coordinates ) == VERDICT_TRUE ) return v_tri_minimum_angle( 3, coordinates );
[ # # ]
736 : :
737 : : double angle;
738 : 1 : double min_angle = 360.0;
739 : :
740 [ + - ][ + + ]: 5 : VerdictVector edges[4];
741 : 2 : edges[0].set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
742 [ + - ]: 1 : coordinates[1][2] - coordinates[0][2] );
743 : 2 : edges[1].set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
744 [ + - ]: 1 : coordinates[2][2] - coordinates[1][2] );
745 : 2 : edges[2].set( coordinates[3][0] - coordinates[2][0], coordinates[3][1] - coordinates[2][1],
746 [ + - ]: 1 : coordinates[3][2] - coordinates[2][2] );
747 : 2 : edges[3].set( coordinates[0][0] - coordinates[3][0], coordinates[0][1] - coordinates[3][1],
748 [ + - ]: 1 : coordinates[0][2] - coordinates[3][2] );
749 : :
750 : : // go around each node and calculate the angle
751 : : // at each node
752 : : double length[4];
753 [ + - ]: 1 : length[0] = edges[0].length();
754 [ + - ]: 1 : length[1] = edges[1].length();
755 [ + - ]: 1 : length[2] = edges[2].length();
756 [ + - ]: 1 : length[3] = edges[3].length();
757 : :
758 [ + - ][ + - ]: 1 : if( length[0] <= VERDICT_DBL_MIN || length[1] <= VERDICT_DBL_MIN || length[2] <= VERDICT_DBL_MIN ||
[ + - ][ - + ]
759 : 1 : length[3] <= VERDICT_DBL_MIN )
760 : 0 : return 360.0;
761 : :
762 [ + - ]: 1 : angle = acos( -( edges[0] % edges[1] ) / ( length[0] * length[1] ) );
763 [ + - ]: 1 : min_angle = VERDICT_MIN( angle, min_angle );
764 : :
765 [ + - ]: 1 : angle = acos( -( edges[1] % edges[2] ) / ( length[1] * length[2] ) );
766 [ - + ]: 1 : min_angle = VERDICT_MIN( angle, min_angle );
767 : :
768 [ + - ]: 1 : angle = acos( -( edges[2] % edges[3] ) / ( length[2] * length[3] ) );
769 [ - + ]: 1 : min_angle = VERDICT_MIN( angle, min_angle );
770 : :
771 [ + - ]: 1 : angle = acos( -( edges[3] % edges[0] ) / ( length[3] * length[0] ) );
772 [ - + ]: 1 : min_angle = VERDICT_MIN( angle, min_angle );
773 : :
774 : 1 : min_angle = min_angle * 180.0 / VERDICT_PI;
775 : :
776 [ + - ][ + - ]: 1 : if( min_angle > 0 ) return (double)VERDICT_MIN( min_angle, VERDICT_DBL_MAX );
777 [ # # ]: 1 : return (double)VERDICT_MAX( min_angle, -VERDICT_DBL_MAX );
778 : : }
779 : :
780 : : /*!
781 : : the oddy of a quad
782 : :
783 : : general distortion measure based on left Cauchy-Green Tensor
784 : : */
785 : 8 : C_FUNC_DEF double v_quad_oddy( int /*num_nodes*/, double coordinates[][3] )
786 : : {
787 : :
788 : 8 : double max_oddy = 0.;
789 : :
790 [ + - ][ + - ]: 40 : VerdictVector first, second, node_pos[4];
[ + - ][ + + ]
791 : :
792 : : double g, g11, g12, g22, cur_oddy;
793 : : int i;
794 : :
795 [ + + ]: 40 : for( i = 0; i < 4; i++ )
796 [ + - ]: 32 : node_pos[i].set( coordinates[i][0], coordinates[i][1], coordinates[i][2] );
797 : :
798 [ + + ]: 40 : for( i = 0; i < 4; i++ )
799 : : {
800 [ + - ][ + - ]: 32 : first = node_pos[i] - node_pos[( i + 1 ) % 4];
801 [ + - ][ + - ]: 32 : second = node_pos[i] - node_pos[( i + 3 ) % 4];
802 : :
803 [ + - ]: 32 : g11 = first % first;
804 [ + - ]: 32 : g12 = first % second;
805 [ + - ]: 32 : g22 = second % second;
806 : 32 : g = g11 * g22 - g12 * g12;
807 : :
808 [ - + ]: 32 : if( g < VERDICT_DBL_MIN )
809 : 0 : cur_oddy = VERDICT_DBL_MAX;
810 : : else
811 : 32 : cur_oddy = ( ( g11 - g22 ) * ( g11 - g22 ) + 4. * g12 * g12 ) / 2. / g;
812 : :
813 [ - + ]: 32 : max_oddy = VERDICT_MAX( max_oddy, cur_oddy );
814 : : }
815 : :
816 [ - + ][ # # ]: 8 : if( max_oddy > 0 ) return (double)VERDICT_MIN( max_oddy, VERDICT_DBL_MAX );
817 [ + - ]: 8 : return (double)VERDICT_MAX( max_oddy, -VERDICT_DBL_MAX );
818 : : }
819 : :
820 : : /*!
821 : : the condition of a quad
822 : :
823 : : maximum condition number of the Jacobian matrix at 4 corners
824 : : */
825 : 9 : C_FUNC_DEF double v_quad_condition( int /*num_nodes*/, double coordinates[][3] )
826 : : {
827 : :
828 [ + - ][ - + ]: 9 : if( is_collapsed_quad( coordinates ) == VERDICT_TRUE ) return v_tri_condition( 3, coordinates );
[ # # ]
829 : :
830 : : double areas[4];
831 [ + - ]: 9 : signed_corner_areas( areas, coordinates );
832 : :
833 : 9 : double max_condition = 0.;
834 : :
835 [ + - ][ + - ]: 9 : VerdictVector xxi, xet;
836 : :
837 : : double condition;
838 : :
839 [ + + ]: 45 : for( int i = 0; i < 4; i++ )
840 : : {
841 : :
842 : 72 : xxi.set( coordinates[i][0] - coordinates[( i + 1 ) % 4][0], coordinates[i][1] - coordinates[( i + 1 ) % 4][1],
843 [ + - ]: 36 : coordinates[i][2] - coordinates[( i + 1 ) % 4][2] );
844 : :
845 : 72 : xet.set( coordinates[i][0] - coordinates[( i + 3 ) % 4][0], coordinates[i][1] - coordinates[( i + 3 ) % 4][1],
846 [ + - ]: 36 : coordinates[i][2] - coordinates[( i + 3 ) % 4][2] );
847 : :
848 [ - + ]: 36 : if( areas[i] < VERDICT_DBL_MIN )
849 : 0 : condition = VERDICT_DBL_MAX;
850 : : else
851 [ + - ][ + - ]: 36 : condition = ( xxi % xxi + xet % xet ) / areas[i];
852 : :
853 [ + + ]: 36 : max_condition = VERDICT_MAX( max_condition, condition );
854 : : }
855 : :
856 : 9 : max_condition /= 2;
857 : :
858 [ + - ][ + - ]: 9 : if( max_condition > 0 ) return (double)VERDICT_MIN( max_condition, VERDICT_DBL_MAX );
859 [ # # ]: 9 : return (double)VERDICT_MAX( max_condition, -VERDICT_DBL_MAX );
860 : : }
861 : :
862 : : /*!
863 : : the jacobian of a quad
864 : :
865 : : minimum pointwise volume of local map at 4 corners and center of quad
866 : : */
867 : 9 : C_FUNC_DEF double v_quad_jacobian( int /*num_nodes*/, double coordinates[][3] )
868 : : {
869 : :
870 [ + - ][ - + ]: 9 : if( is_collapsed_quad( coordinates ) == VERDICT_TRUE ) return (double)( v_tri_area( 3, coordinates ) * 2.0 );
[ # # ]
871 : :
872 : : double areas[4];
873 [ + - ]: 9 : signed_corner_areas( areas, coordinates );
874 : :
875 [ + + ][ - + ]: 9 : double jacobian = VERDICT_MIN( VERDICT_MIN( areas[0], areas[1] ), VERDICT_MIN( areas[2], areas[3] ) );
[ - + ][ # # ]
[ - + ]
876 [ + - ][ + - ]: 9 : if( jacobian > 0 ) return (double)VERDICT_MIN( jacobian, VERDICT_DBL_MAX );
877 [ # # ]: 9 : return (double)VERDICT_MAX( jacobian, -VERDICT_DBL_MAX );
878 : : }
879 : :
880 : : /*!
881 : : scaled jacobian of a quad
882 : :
883 : : Minimum Jacobian divided by the lengths of the 2 edge vector
884 : : */
885 : 27 : C_FUNC_DEF double v_quad_scaled_jacobian( int /*num_nodes*/, double coordinates[][3] )
886 : : {
887 [ + - ][ - + ]: 27 : if( is_collapsed_quad( coordinates ) == VERDICT_TRUE ) return v_tri_scaled_jacobian( 3, coordinates );
[ # # ]
888 : :
889 : 27 : double corner_areas[4], min_scaled_jac = VERDICT_DBL_MAX, scaled_jac;
890 [ + - ]: 27 : signed_corner_areas( corner_areas, coordinates );
891 : :
892 [ + - ][ + + ]: 135 : VerdictVector edges[4];
893 [ + - ]: 27 : make_quad_edges( edges, coordinates );
894 : :
895 : : double length[4];
896 [ + - ]: 27 : length[0] = edges[0].length();
897 [ + - ]: 27 : length[1] = edges[1].length();
898 [ + - ]: 27 : length[2] = edges[2].length();
899 [ + - ]: 27 : length[3] = edges[3].length();
900 : :
901 [ + - ][ + - ]: 27 : if( length[0] < VERDICT_DBL_MIN || length[1] < VERDICT_DBL_MIN || length[2] < VERDICT_DBL_MIN ||
[ + - ][ - + ]
902 : 27 : length[3] < VERDICT_DBL_MIN )
903 : 0 : return 0.0;
904 : :
905 : 27 : scaled_jac = corner_areas[0] / ( length[0] * length[3] );
906 [ + - ]: 27 : min_scaled_jac = VERDICT_MIN( scaled_jac, min_scaled_jac );
907 : :
908 : 27 : scaled_jac = corner_areas[1] / ( length[1] * length[0] );
909 [ + + ]: 27 : min_scaled_jac = VERDICT_MIN( scaled_jac, min_scaled_jac );
910 : :
911 : 27 : scaled_jac = corner_areas[2] / ( length[2] * length[1] );
912 [ + + ]: 27 : min_scaled_jac = VERDICT_MIN( scaled_jac, min_scaled_jac );
913 : :
914 : 27 : scaled_jac = corner_areas[3] / ( length[3] * length[2] );
915 [ + + ]: 27 : min_scaled_jac = VERDICT_MIN( scaled_jac, min_scaled_jac );
916 : :
917 [ + - ][ + - ]: 27 : if( min_scaled_jac > 0 ) return (double)VERDICT_MIN( min_scaled_jac, VERDICT_DBL_MAX );
918 [ # # ]: 27 : return (double)VERDICT_MAX( min_scaled_jac, -VERDICT_DBL_MAX );
919 : : }
920 : :
921 : : /*!
922 : : the shear of a quad
923 : :
924 : : 2/Condition number of Jacobian Skew matrix
925 : : */
926 : 18 : C_FUNC_DEF double v_quad_shear( int /*num_nodes*/, double coordinates[][3] )
927 : : {
928 : 18 : double scaled_jacobian = v_quad_scaled_jacobian( 4, coordinates );
929 : :
930 [ - + ]: 18 : if( scaled_jacobian <= VERDICT_DBL_MIN )
931 : 0 : return 0.0;
932 : : else
933 [ + - ]: 18 : return (double)VERDICT_MIN( scaled_jacobian, VERDICT_DBL_MAX );
934 : : }
935 : :
936 : : /*!
937 : : the shape of a quad
938 : :
939 : : 2/Condition number of weighted Jacobian matrix
940 : : */
941 : 18 : C_FUNC_DEF double v_quad_shape( int /*num_nodes*/, double coordinates[][3] )
942 : : {
943 : :
944 : 18 : double corner_areas[4], min_shape = VERDICT_DBL_MAX, shape;
945 [ + - ]: 18 : signed_corner_areas( corner_areas, coordinates );
946 : :
947 [ + - ][ + + ]: 90 : VerdictVector edges[4];
948 [ + - ]: 18 : make_quad_edges( edges, coordinates );
949 : :
950 : : double length_squared[4];
951 [ + - ]: 18 : length_squared[0] = edges[0].length_squared();
952 [ + - ]: 18 : length_squared[1] = edges[1].length_squared();
953 [ + - ]: 18 : length_squared[2] = edges[2].length_squared();
954 [ + - ]: 18 : length_squared[3] = edges[3].length_squared();
955 : :
956 [ + - ][ + - ]: 18 : if( length_squared[0] <= VERDICT_DBL_MIN || length_squared[1] <= VERDICT_DBL_MIN ||
[ + - ]
957 [ - + ]: 18 : length_squared[2] <= VERDICT_DBL_MIN || length_squared[3] <= VERDICT_DBL_MIN )
958 : 0 : return 0.0;
959 : :
960 : 18 : shape = corner_areas[0] / ( length_squared[0] + length_squared[3] );
961 [ + - ]: 18 : min_shape = VERDICT_MIN( shape, min_shape );
962 : :
963 : 18 : shape = corner_areas[1] / ( length_squared[1] + length_squared[0] );
964 [ + + ]: 18 : min_shape = VERDICT_MIN( shape, min_shape );
965 : :
966 : 18 : shape = corner_areas[2] / ( length_squared[2] + length_squared[1] );
967 [ + + ]: 18 : min_shape = VERDICT_MIN( shape, min_shape );
968 : :
969 : 18 : shape = corner_areas[3] / ( length_squared[3] + length_squared[2] );
970 [ - + ]: 18 : min_shape = VERDICT_MIN( shape, min_shape );
971 : :
972 : 18 : min_shape *= 2;
973 : :
974 [ - + ]: 18 : if( min_shape < VERDICT_DBL_MIN ) min_shape = 0;
975 : :
976 [ + - ][ + - ]: 18 : if( min_shape > 0 ) return (double)VERDICT_MIN( min_shape, VERDICT_DBL_MAX );
977 [ # # ]: 18 : return (double)VERDICT_MAX( min_shape, -VERDICT_DBL_MAX );
978 : : }
979 : :
980 : : /*!
981 : : the relative size of a quad
982 : :
983 : : Min( J, 1/J ), where J is determinant of weighted Jacobian matrix
984 : : */
985 : 27 : C_FUNC_DEF double v_quad_relative_size_squared( int /*num_nodes*/, double coordinates[][3] )
986 : : {
987 : :
988 [ + - ]: 27 : double quad_area = v_quad_area( 4, coordinates );
989 : 27 : double rel_size = 0;
990 : :
991 [ + - ]: 27 : v_set_quad_size( quad_area );
992 : : double w11, w21, w12, w22;
993 [ + - ]: 27 : get_weight( w11, w21, w12, w22 );
994 [ + - ]: 27 : double avg_area = determinant( w11, w21, w12, w22 );
995 : :
996 [ + - ]: 27 : if( avg_area > VERDICT_DBL_MIN )
997 : : {
998 : :
999 : 27 : w11 = quad_area / avg_area;
1000 : :
1001 [ + - ]: 27 : if( w11 > VERDICT_DBL_MIN )
1002 : : {
1003 [ - + ]: 27 : rel_size = VERDICT_MIN( w11, 1 / w11 );
1004 : 27 : rel_size *= rel_size;
1005 : : }
1006 : : }
1007 : :
1008 [ + - ][ + - ]: 27 : if( rel_size > 0 ) return (double)VERDICT_MIN( rel_size, VERDICT_DBL_MAX );
1009 [ # # ]: 27 : return (double)VERDICT_MAX( rel_size, -VERDICT_DBL_MAX );
1010 : : }
1011 : :
1012 : : /*!
1013 : : the relative shape and size of a quad
1014 : :
1015 : : Product of Shape and Relative Size
1016 : : */
1017 : 9 : C_FUNC_DEF double v_quad_shape_and_size( int num_nodes, double coordinates[][3] )
1018 : : {
1019 : : double shape, size;
1020 : 9 : size = v_quad_relative_size_squared( num_nodes, coordinates );
1021 : 9 : shape = v_quad_shape( num_nodes, coordinates );
1022 : :
1023 : 9 : double shape_and_size = shape * size;
1024 : :
1025 [ + - ][ + - ]: 9 : if( shape_and_size > 0 ) return (double)VERDICT_MIN( shape_and_size, VERDICT_DBL_MAX );
1026 [ # # ]: 0 : return (double)VERDICT_MAX( shape_and_size, -VERDICT_DBL_MAX );
1027 : : }
1028 : :
1029 : : /*!
1030 : : the shear and size of a quad
1031 : :
1032 : : product of shear and relative size
1033 : : */
1034 : 9 : C_FUNC_DEF double v_quad_shear_and_size( int num_nodes, double coordinates[][3] )
1035 : : {
1036 : : double shear, size;
1037 : 9 : shear = v_quad_shear( num_nodes, coordinates );
1038 : 9 : size = v_quad_relative_size_squared( num_nodes, coordinates );
1039 : :
1040 : 9 : double shear_and_size = shear * size;
1041 : :
1042 [ + - ][ + - ]: 9 : if( shear_and_size > 0 ) return (double)VERDICT_MIN( shear_and_size, VERDICT_DBL_MAX );
1043 [ # # ]: 0 : return (double)VERDICT_MAX( shear_and_size, -VERDICT_DBL_MAX );
1044 : : }
1045 : :
1046 : : /*!
1047 : : the distortion of a quad
1048 : : */
1049 : 17 : C_FUNC_DEF double v_quad_distortion( int num_nodes, double coordinates[][3] )
1050 : : {
1051 : : // To calculate distortion for linear and 2nd order quads
1052 : : // distortion = {min(|J|)/actual area}*{parent area}
1053 : : // parent area = 4 for a quad.
1054 : : // min |J| is the minimum over nodes and gaussian integration points
1055 : : // created by Ling Pan, CAT on 4/30/01
1056 : :
1057 : 17 : double element_area = 0.0, distrt, thickness_gauss;
1058 : 17 : double cur_jacobian = 0., sign_jacobian, jacobian;
1059 [ + - ][ + - ]: 17 : VerdictVector aa, bb, cc, normal_at_point, xin;
[ + - ][ + - ]
[ + - ]
1060 : :
1061 : : // use 2x2 gauss points for linear quads and 3x3 for 2nd order quads
1062 : 17 : int number_of_gauss_points = 0;
1063 [ + - ]: 17 : if( num_nodes == 4 )
1064 : : { // 2x2 quadrature rule
1065 : 17 : number_of_gauss_points = 2;
1066 : : }
1067 [ # # ]: 0 : else if( num_nodes == 8 )
1068 : : { // 3x3 quadrature rule
1069 : 0 : number_of_gauss_points = 3;
1070 : : }
1071 : :
1072 : 17 : int total_number_of_gauss_points = number_of_gauss_points * number_of_gauss_points;
1073 : :
1074 [ + - ]: 17 : VerdictVector face_normal = quad_normal( coordinates );
1075 : :
1076 : 17 : double distortion = VERDICT_DBL_MAX;
1077 : :
1078 [ + - ][ + - ]: 17 : VerdictVector first, second;
1079 : :
1080 : : int i;
1081 : : // Will work out the case for collapsed quad later
1082 [ + - ][ - + ]: 17 : if( is_collapsed_quad( coordinates ) == VERDICT_TRUE )
1083 : : {
1084 [ # # ]: 0 : for( i = 0; i < 3; i++ )
1085 : : {
1086 : :
1087 : 0 : first.set( coordinates[i][0] - coordinates[( i + 1 ) % 3][0],
1088 : 0 : coordinates[i][1] - coordinates[( i + 1 ) % 3][1],
1089 [ # # ]: 0 : coordinates[i][2] - coordinates[( i + 1 ) % 3][2] );
1090 : :
1091 : 0 : second.set( coordinates[i][0] - coordinates[( i + 2 ) % 3][0],
1092 : 0 : coordinates[i][1] - coordinates[( i + 2 ) % 3][1],
1093 [ # # ]: 0 : coordinates[i][2] - coordinates[( i + 2 ) % 3][2] );
1094 : :
1095 [ # # ][ # # ]: 0 : sign_jacobian = ( face_normal % ( first * second ) ) > 0 ? 1. : -1.;
[ # # ]
1096 [ # # ][ # # ]: 0 : cur_jacobian = sign_jacobian * ( first * second ).length();
1097 [ # # ]: 0 : distortion = VERDICT_MIN( distortion, cur_jacobian );
1098 : : }
1099 [ # # ][ # # ]: 0 : element_area = ( first * second ).length() / 2.0;
1100 : 0 : distortion /= element_area;
1101 : : }
1102 : : else
1103 : : {
1104 : : double shape_function[maxTotalNumberGaussPoints][maxNumberNodes];
1105 : : double dndy1[maxTotalNumberGaussPoints][maxNumberNodes];
1106 : : double dndy2[maxTotalNumberGaussPoints][maxNumberNodes];
1107 : : double weight[maxTotalNumberGaussPoints];
1108 : :
1109 : : // create an object of GaussIntegration
1110 [ + - ]: 17 : GaussIntegration::initialize( number_of_gauss_points, num_nodes );
1111 [ + - ]: 17 : GaussIntegration::calculate_shape_function_2d_quad();
1112 [ + - ]: 17 : GaussIntegration::get_shape_func( shape_function[0], dndy1[0], dndy2[0], weight );
1113 : :
1114 : : // calculate element area
1115 : : int ife, ja;
1116 [ + + ]: 85 : for( ife = 0; ife < total_number_of_gauss_points; ife++ )
1117 : : {
1118 [ + - ]: 68 : aa.set( 0.0, 0.0, 0.0 );
1119 [ + - ]: 68 : bb.set( 0.0, 0.0, 0.0 );
1120 : :
1121 [ + + ]: 340 : for( ja = 0; ja < num_nodes; ja++ )
1122 : : {
1123 [ + - ]: 272 : xin.set( coordinates[ja][0], coordinates[ja][1], coordinates[ja][2] );
1124 [ + - ][ + - ]: 272 : aa += dndy1[ife][ja] * xin;
1125 [ + - ][ + - ]: 272 : bb += dndy2[ife][ja] * xin;
1126 : : }
1127 [ + - ][ + - ]: 68 : normal_at_point = aa * bb;
1128 [ + - ]: 68 : jacobian = normal_at_point.length();
1129 : 68 : element_area += weight[ife] * jacobian;
1130 : : }
1131 : :
1132 : : double dndy1_at_node[maxNumberNodes][maxNumberNodes];
1133 : : double dndy2_at_node[maxNumberNodes][maxNumberNodes];
1134 : :
1135 [ + - ]: 17 : GaussIntegration::calculate_derivative_at_nodes( dndy1_at_node, dndy2_at_node );
1136 : :
1137 [ + - ][ + + ]: 170 : VerdictVector normal_at_nodes[9];
1138 : :
1139 : : // evaluate normal at nodes and distortion values at nodes
1140 : : int jai;
1141 [ + + ]: 85 : for( ja = 0; ja < num_nodes; ja++ )
1142 : : {
1143 [ + - ]: 68 : aa.set( 0.0, 0.0, 0.0 );
1144 [ + - ]: 68 : bb.set( 0.0, 0.0, 0.0 );
1145 [ + + ]: 340 : for( jai = 0; jai < num_nodes; jai++ )
1146 : : {
1147 [ + - ]: 272 : xin.set( coordinates[jai][0], coordinates[jai][1], coordinates[jai][2] );
1148 [ + - ][ + - ]: 272 : aa += dndy1_at_node[ja][jai] * xin;
1149 [ + - ][ + - ]: 272 : bb += dndy2_at_node[ja][jai] * xin;
1150 : : }
1151 [ + - ][ + - ]: 68 : normal_at_nodes[ja] = aa * bb;
1152 [ + - ]: 68 : normal_at_nodes[ja].normalize();
1153 : : }
1154 : :
1155 : : // determine if element is flat
1156 : 17 : bool flat_element = true;
1157 : : double dot_product;
1158 : :
1159 [ + + ]: 82 : for( ja = 0; ja < num_nodes; ja++ )
1160 : : {
1161 [ + - ]: 66 : dot_product = normal_at_nodes[0] % normal_at_nodes[ja];
1162 [ + + ]: 66 : if( fabs( dot_product ) < 0.99 )
1163 : : {
1164 : 1 : flat_element = false;
1165 : 1 : break;
1166 : : }
1167 : : }
1168 : :
1169 : : // take into consideration of the thickness of the element
1170 : : double thickness;
1171 : : // get_quad_thickness(face, element_area, thickness );
1172 : 17 : thickness = 0.001 * sqrt( element_area );
1173 : :
1174 : : // set thickness gauss point location
1175 : 17 : double zl = 0.5773502691896;
1176 [ + + ]: 17 : if( flat_element ) zl = 0.0;
1177 : :
1178 [ + + ]: 17 : int no_gauss_pts_z = ( flat_element ) ? 1 : 2;
1179 : : double thickness_z;
1180 : : int igz;
1181 : : // loop on Gauss points
1182 [ + + ]: 85 : for( ife = 0; ife < total_number_of_gauss_points; ife++ )
1183 : : {
1184 : : // loop on the thickness direction gauss points
1185 [ + + ]: 140 : for( igz = 0; igz < no_gauss_pts_z; igz++ )
1186 : : {
1187 : 72 : zl = -zl;
1188 : 72 : thickness_z = zl * thickness / 2.0;
1189 : :
1190 [ + - ]: 72 : aa.set( 0.0, 0.0, 0.0 );
1191 [ + - ]: 72 : bb.set( 0.0, 0.0, 0.0 );
1192 [ + - ]: 72 : cc.set( 0.0, 0.0, 0.0 );
1193 : :
1194 [ + + ]: 360 : for( ja = 0; ja < num_nodes; ja++ )
1195 : : {
1196 [ + - ]: 288 : xin.set( coordinates[ja][0], coordinates[ja][1], coordinates[ja][2] );
1197 [ + - ][ + - ]: 288 : xin += thickness_z * normal_at_nodes[ja];
1198 [ + - ][ + - ]: 288 : aa += dndy1[ife][ja] * xin;
1199 [ + - ][ + - ]: 288 : bb += dndy2[ife][ja] * xin;
1200 : 288 : thickness_gauss = shape_function[ife][ja] * thickness / 2.0;
1201 [ + - ][ + - ]: 288 : cc += thickness_gauss * normal_at_nodes[ja];
1202 : : }
1203 : :
1204 [ + - ][ + - ]: 72 : normal_at_point = aa * bb;
1205 : : // jacobian = normal_at_point.length();
1206 [ + - ]: 72 : distrt = cc % normal_at_point;
1207 [ + + ]: 72 : if( distrt < distortion ) distortion = distrt;
1208 : : }
1209 : : }
1210 : :
1211 : : // loop through nodal points
1212 [ + + ]: 85 : for( ja = 0; ja < num_nodes; ja++ )
1213 : : {
1214 [ + + ]: 140 : for( igz = 0; igz < no_gauss_pts_z; igz++ )
1215 : : {
1216 : 72 : zl = -zl;
1217 : 72 : thickness_z = zl * thickness / 2.0;
1218 : :
1219 [ + - ]: 72 : aa.set( 0.0, 0.0, 0.0 );
1220 [ + - ]: 72 : bb.set( 0.0, 0.0, 0.0 );
1221 [ + - ]: 72 : cc.set( 0.0, 0.0, 0.0 );
1222 : :
1223 [ + + ]: 360 : for( jai = 0; jai < num_nodes; jai++ )
1224 : : {
1225 [ + - ]: 288 : xin.set( coordinates[jai][0], coordinates[jai][1], coordinates[jai][2] );
1226 [ + - ][ + - ]: 288 : xin += thickness_z * normal_at_nodes[ja];
1227 [ + - ][ + - ]: 288 : aa += dndy1_at_node[ja][jai] * xin;
1228 [ + - ][ + - ]: 288 : bb += dndy2_at_node[ja][jai] * xin;
1229 [ + + ]: 288 : if( jai == ja )
1230 : 72 : thickness_gauss = thickness / 2.0;
1231 : : else
1232 : 216 : thickness_gauss = 0.;
1233 [ + - ][ + - ]: 288 : cc += thickness_gauss * normal_at_nodes[jai];
1234 : : }
1235 : : }
1236 [ + - ][ + - ]: 68 : normal_at_point = aa * bb;
1237 [ + - ][ + - ]: 68 : sign_jacobian = ( face_normal % normal_at_point ) > 0 ? 1. : -1.;
1238 [ + - ]: 68 : distrt = sign_jacobian * ( cc % normal_at_point );
1239 : :
1240 [ - + ]: 68 : if( distrt < distortion ) distortion = distrt;
1241 : : }
1242 : :
1243 [ + - ]: 17 : if( element_area * thickness != 0 )
1244 : 17 : distortion *= 8. / ( element_area * thickness );
1245 : : else
1246 : 17 : distortion *= 8.;
1247 : : }
1248 : :
1249 : 17 : return (double)distortion;
1250 : : }
1251 : :
1252 : : /*!
1253 : : multiple quality measures of a quad
1254 : : */
1255 : 8 : C_FUNC_DEF void v_quad_quality( int num_nodes, double coordinates[][3], unsigned int metrics_request_flag,
1256 : : QuadMetricVals* metric_vals )
1257 : : {
1258 : :
1259 : 8 : memset( metric_vals, 0, sizeof( QuadMetricVals ) );
1260 : :
1261 : : // for starts, lets set up some basic and common information
1262 : :
1263 : : /* node numbers and side numbers used below
1264 : :
1265 : : 2
1266 : : 3 +--------- 2
1267 : : / +
1268 : : / |
1269 : : 3 / | 1
1270 : : / |
1271 : : + |
1272 : : 0 -------------+ 1
1273 : : 0
1274 : : */
1275 : :
1276 : : // vectors for each side
1277 [ + - ][ + + ]: 40 : VerdictVector edges[4];
1278 [ + - ]: 8 : make_quad_edges( edges, coordinates );
1279 : :
1280 : : double areas[4];
1281 [ + - ]: 8 : signed_corner_areas( areas, coordinates );
1282 : :
1283 : : double lengths[4];
1284 [ + - ]: 8 : lengths[0] = edges[0].length();
1285 [ + - ]: 8 : lengths[1] = edges[1].length();
1286 [ + - ]: 8 : lengths[2] = edges[2].length();
1287 [ + - ]: 8 : lengths[3] = edges[3].length();
1288 : :
1289 [ + - ]: 8 : VerdictBoolean is_collapsed = is_collapsed_quad( coordinates );
1290 : :
1291 : : // handle collapsed quads metrics here
1292 [ - + ][ # # ]: 8 : if( is_collapsed == VERDICT_TRUE && metrics_request_flag & ( V_QUAD_MINIMUM_ANGLE | V_QUAD_MAXIMUM_ANGLE |
1293 : : V_QUAD_JACOBIAN | V_QUAD_SCALED_JACOBIAN ) )
1294 : : {
1295 [ # # ]: 0 : if( metrics_request_flag & V_QUAD_MINIMUM_ANGLE )
1296 [ # # ]: 0 : metric_vals->minimum_angle = v_tri_minimum_angle( 3, coordinates );
1297 [ # # ]: 0 : if( metrics_request_flag & V_QUAD_MAXIMUM_ANGLE )
1298 [ # # ]: 0 : metric_vals->maximum_angle = v_tri_maximum_angle( 3, coordinates );
1299 [ # # ]: 0 : if( metrics_request_flag & V_QUAD_JACOBIAN )
1300 [ # # ]: 0 : metric_vals->jacobian = (double)( v_tri_area( 3, coordinates ) * 2.0 );
1301 [ # # ]: 0 : if( metrics_request_flag & V_QUAD_SCALED_JACOBIAN )
1302 [ # # ]: 0 : metric_vals->jacobian = (double)( v_tri_scaled_jacobian( 3, coordinates ) * 2.0 );
1303 : : }
1304 : :
1305 : : // calculate both largest and smallest angles
1306 [ + - ][ + - ]: 8 : if( metrics_request_flag & ( V_QUAD_MINIMUM_ANGLE | V_QUAD_MAXIMUM_ANGLE ) && is_collapsed == VERDICT_FALSE )
1307 : : {
1308 : : // gather the angles
1309 : : double angles[4];
1310 [ + - ]: 8 : angles[0] = acos( -( edges[0] % edges[1] ) / ( lengths[0] * lengths[1] ) );
1311 [ + - ]: 8 : angles[1] = acos( -( edges[1] % edges[2] ) / ( lengths[1] * lengths[2] ) );
1312 [ + - ]: 8 : angles[2] = acos( -( edges[2] % edges[3] ) / ( lengths[2] * lengths[3] ) );
1313 [ + - ]: 8 : angles[3] = acos( -( edges[3] % edges[0] ) / ( lengths[3] * lengths[0] ) );
1314 : :
1315 [ + - ][ + - ]: 8 : if( lengths[0] <= VERDICT_DBL_MIN || lengths[1] <= VERDICT_DBL_MIN || lengths[2] <= VERDICT_DBL_MIN ||
[ + - ][ - + ]
1316 : 8 : lengths[3] <= VERDICT_DBL_MIN )
1317 : : {
1318 : 0 : metric_vals->minimum_angle = 360.0;
1319 : 0 : metric_vals->maximum_angle = 0.0;
1320 : : }
1321 : : else
1322 : : {
1323 : : // if smallest angle, find the smallest angle
1324 [ + - ]: 8 : if( metrics_request_flag & V_QUAD_MINIMUM_ANGLE )
1325 : : {
1326 : 8 : metric_vals->minimum_angle = VERDICT_DBL_MAX;
1327 [ + + ]: 40 : for( int i = 0; i < 4; i++ )
1328 [ + + ]: 32 : metric_vals->minimum_angle = VERDICT_MIN( angles[i], metric_vals->minimum_angle );
1329 : 8 : metric_vals->minimum_angle *= 180.0 / VERDICT_PI;
1330 : : }
1331 : : // if largest angle, find the largest angle
1332 [ + - ]: 8 : if( metrics_request_flag & V_QUAD_MAXIMUM_ANGLE )
1333 : : {
1334 : 8 : metric_vals->maximum_angle = 0.0;
1335 [ + + ]: 40 : for( int i = 0; i < 4; i++ )
1336 [ + + ]: 32 : metric_vals->maximum_angle = VERDICT_MAX( angles[i], metric_vals->maximum_angle );
1337 : 8 : metric_vals->maximum_angle *= 180.0 / VERDICT_PI;
1338 : :
1339 [ + - ][ + - ]: 8 : if( areas[0] < 0 || areas[1] < 0 || areas[2] < 0 || areas[3] < 0 )
[ + - ][ - + ]
1340 : 8 : metric_vals->maximum_angle = 360 - metric_vals->maximum_angle;
1341 : : }
1342 : : }
1343 : : }
1344 : :
1345 : : // handle max_edge_ratio, skew, taper, and area together
1346 [ + - ]: 8 : if( metrics_request_flag & ( V_QUAD_MAX_EDGE_RATIO | V_QUAD_SKEW | V_QUAD_TAPER ) )
1347 : : {
1348 : : // get principle axes
1349 [ + - ][ + + ]: 24 : VerdictVector principal_axes[2];
1350 [ + - ][ + - ]: 8 : principal_axes[0] = edges[0] - edges[2];
1351 [ + - ][ + - ]: 8 : principal_axes[1] = edges[1] - edges[3];
1352 : :
1353 [ + - ]: 8 : if( metrics_request_flag & ( V_QUAD_MAX_EDGE_RATIO | V_QUAD_SKEW | V_QUAD_TAPER ) )
1354 : : {
1355 [ + - ]: 8 : double len1 = principal_axes[0].length();
1356 [ + - ]: 8 : double len2 = principal_axes[1].length();
1357 : :
1358 : : // calculate the max_edge_ratio ratio
1359 [ + - ]: 8 : if( metrics_request_flag & V_QUAD_MAX_EDGE_RATIO )
1360 : : {
1361 [ + - ][ - + ]: 8 : if( len1 < VERDICT_DBL_MIN || len2 < VERDICT_DBL_MIN )
1362 : 0 : metric_vals->max_edge_ratio = VERDICT_DBL_MAX;
1363 : : else
1364 [ - + ]: 8 : metric_vals->max_edge_ratio = VERDICT_MAX( len1 / len2, len2 / len1 );
1365 : : }
1366 : :
1367 : : // calculate the taper
1368 [ + - ]: 8 : if( metrics_request_flag & V_QUAD_TAPER )
1369 : : {
1370 [ - + ]: 8 : double min_length = VERDICT_MIN( len1, len2 );
1371 : :
1372 [ + - ]: 8 : VerdictVector cross_derivative = edges[1] + edges[3];
1373 : :
1374 [ - + ]: 8 : if( min_length < VERDICT_DBL_MIN )
1375 : 0 : metric_vals->taper = VERDICT_DBL_MAX;
1376 : : else
1377 [ + - ]: 8 : metric_vals->taper = cross_derivative.length() / min_length;
1378 : : }
1379 : :
1380 : : // calculate the skew
1381 [ + - ]: 8 : if( metrics_request_flag & V_QUAD_SKEW )
1382 : : {
1383 [ + - ][ + - ]: 8 : if( principal_axes[0].normalize() < VERDICT_DBL_MIN || principal_axes[1].normalize() < VERDICT_DBL_MIN )
[ + - ][ - + ]
[ - + ]
1384 : 0 : metric_vals->skew = 0.0;
1385 : : else
1386 [ + - ]: 8 : metric_vals->skew = fabs( principal_axes[0] % principal_axes[1] );
1387 : : }
1388 : : }
1389 : : }
1390 : :
1391 : : // calculate the area
1392 [ + - ]: 8 : if( metrics_request_flag & ( V_QUAD_AREA | V_QUAD_RELATIVE_SIZE_SQUARED ) )
1393 : 8 : { metric_vals->area = 0.25 * ( areas[0] + areas[1] + areas[2] + areas[3] ); }
1394 : :
1395 : : // calculate the relative size
1396 [ + - ]: 8 : if( metrics_request_flag & ( V_QUAD_RELATIVE_SIZE_SQUARED | V_QUAD_SHAPE_AND_SIZE | V_QUAD_SHEAR_AND_SIZE ) )
1397 : : {
1398 [ + - ]: 8 : double quad_area = v_quad_area( 4, coordinates );
1399 [ + - ]: 8 : v_set_quad_size( quad_area );
1400 : : double w11, w21, w12, w22;
1401 [ + - ]: 8 : get_weight( w11, w21, w12, w22 );
1402 [ + - ]: 8 : double avg_area = determinant( w11, w21, w12, w22 );
1403 : :
1404 [ - + ]: 8 : if( avg_area < VERDICT_DBL_MIN )
1405 : 0 : metric_vals->relative_size_squared = 0.0;
1406 : : else
1407 : : metric_vals->relative_size_squared =
1408 [ - + ]: 8 : pow( VERDICT_MIN( metric_vals->area / avg_area, avg_area / metric_vals->area ), 2 );
1409 : : }
1410 : :
1411 : : // calculate the jacobian
1412 [ + - ]: 8 : if( metrics_request_flag & V_QUAD_JACOBIAN )
1413 [ - + ][ - + ]: 8 : { metric_vals->jacobian = VERDICT_MIN( VERDICT_MIN( areas[0], areas[1] ), VERDICT_MIN( areas[2], areas[3] ) ); }
[ - + ][ # # ]
[ - + ]
1414 : :
1415 [ + - ]: 8 : if( metrics_request_flag & ( V_QUAD_SCALED_JACOBIAN | V_QUAD_SHEAR | V_QUAD_SHEAR_AND_SIZE ) )
1416 : : {
1417 : 8 : double scaled_jac, min_scaled_jac = VERDICT_DBL_MAX;
1418 : :
1419 [ + - ][ + - ]: 8 : if( lengths[0] < VERDICT_DBL_MIN || lengths[1] < VERDICT_DBL_MIN || lengths[2] < VERDICT_DBL_MIN ||
[ + - ][ - + ]
1420 : 8 : lengths[3] < VERDICT_DBL_MIN )
1421 : : {
1422 : 0 : metric_vals->scaled_jacobian = 0.0;
1423 : 0 : metric_vals->shear = 0.0;
1424 : : }
1425 : : else
1426 : : {
1427 : 8 : scaled_jac = areas[0] / ( lengths[0] * lengths[3] );
1428 [ + - ]: 8 : min_scaled_jac = VERDICT_MIN( scaled_jac, min_scaled_jac );
1429 : :
1430 : 8 : scaled_jac = areas[1] / ( lengths[1] * lengths[0] );
1431 [ - + ]: 8 : min_scaled_jac = VERDICT_MIN( scaled_jac, min_scaled_jac );
1432 : :
1433 : 8 : scaled_jac = areas[2] / ( lengths[2] * lengths[1] );
1434 [ - + ]: 8 : min_scaled_jac = VERDICT_MIN( scaled_jac, min_scaled_jac );
1435 : :
1436 : 8 : scaled_jac = areas[3] / ( lengths[3] * lengths[2] );
1437 [ - + ]: 8 : min_scaled_jac = VERDICT_MIN( scaled_jac, min_scaled_jac );
1438 : :
1439 : 8 : metric_vals->scaled_jacobian = min_scaled_jac;
1440 : :
1441 : : // what the heck...set shear as well
1442 [ - + ]: 8 : if( min_scaled_jac <= VERDICT_DBL_MIN )
1443 : 0 : metric_vals->shear = 0.0;
1444 : : else
1445 : 8 : metric_vals->shear = min_scaled_jac;
1446 : : }
1447 : : }
1448 : :
1449 [ + - ]: 8 : if( metrics_request_flag & ( V_QUAD_WARPAGE | V_QUAD_ODDY ) )
1450 : : {
1451 [ + - ][ + + ]: 40 : VerdictVector corner_normals[4];
1452 [ + - ][ + - ]: 8 : corner_normals[0] = edges[3] * edges[0];
1453 [ + - ][ + - ]: 8 : corner_normals[1] = edges[0] * edges[1];
1454 [ + - ][ + - ]: 8 : corner_normals[2] = edges[1] * edges[2];
1455 [ + - ][ + - ]: 8 : corner_normals[3] = edges[2] * edges[3];
1456 : :
1457 [ + - ]: 8 : if( metrics_request_flag & V_QUAD_ODDY )
1458 : : {
1459 : 8 : double oddy, max_oddy = 0.0;
1460 : :
1461 : : double diff, dot_prod;
1462 : :
1463 : : double length_squared[4];
1464 [ + - ]: 8 : length_squared[0] = corner_normals[0].length_squared();
1465 [ + - ]: 8 : length_squared[1] = corner_normals[1].length_squared();
1466 [ + - ]: 8 : length_squared[2] = corner_normals[2].length_squared();
1467 [ + - ]: 8 : length_squared[3] = corner_normals[3].length_squared();
1468 : :
1469 [ + - ][ + - ]: 8 : if( length_squared[0] < VERDICT_DBL_MIN || length_squared[1] < VERDICT_DBL_MIN ||
[ + - ]
1470 [ - + ]: 8 : length_squared[2] < VERDICT_DBL_MIN || length_squared[3] < VERDICT_DBL_MIN )
1471 : 0 : metric_vals->oddy = VERDICT_DBL_MAX;
1472 : : else
1473 : : {
1474 : 8 : diff = ( lengths[0] * lengths[0] ) - ( lengths[1] * lengths[1] );
1475 [ + - ]: 8 : dot_prod = edges[0] % edges[1];
1476 : 8 : oddy = ( ( diff * diff ) + 4 * dot_prod * dot_prod ) / ( 2 * length_squared[1] );
1477 [ - + ]: 8 : max_oddy = VERDICT_MAX( oddy, max_oddy );
1478 : :
1479 : 8 : diff = ( lengths[1] * lengths[1] ) - ( lengths[2] * lengths[2] );
1480 [ + - ]: 8 : dot_prod = edges[1] % edges[2];
1481 : 8 : oddy = ( ( diff * diff ) + 4 * dot_prod * dot_prod ) / ( 2 * length_squared[2] );
1482 [ - + ]: 8 : max_oddy = VERDICT_MAX( oddy, max_oddy );
1483 : :
1484 : 8 : diff = ( lengths[2] * lengths[2] ) - ( lengths[3] * lengths[3] );
1485 [ + - ]: 8 : dot_prod = edges[2] % edges[3];
1486 : 8 : oddy = ( ( diff * diff ) + 4 * dot_prod * dot_prod ) / ( 2 * length_squared[3] );
1487 [ - + ]: 8 : max_oddy = VERDICT_MAX( oddy, max_oddy );
1488 : :
1489 : 8 : diff = ( lengths[3] * lengths[3] ) - ( lengths[0] * lengths[0] );
1490 [ + - ]: 8 : dot_prod = edges[3] % edges[0];
1491 : 8 : oddy = ( ( diff * diff ) + 4 * dot_prod * dot_prod ) / ( 2 * length_squared[0] );
1492 [ - + ]: 8 : max_oddy = VERDICT_MAX( oddy, max_oddy );
1493 : :
1494 : 8 : metric_vals->oddy = max_oddy;
1495 : : }
1496 : : }
1497 : :
1498 [ + - ]: 8 : if( metrics_request_flag & V_QUAD_WARPAGE )
1499 : : {
1500 [ + - ][ + - ]: 32 : if( corner_normals[0].normalize() < VERDICT_DBL_MIN || corner_normals[1].normalize() < VERDICT_DBL_MIN ||
[ + - ][ + - ]
[ - + ]
1501 [ + - ][ + - ]: 24 : corner_normals[2].normalize() < VERDICT_DBL_MIN || corner_normals[3].normalize() < VERDICT_DBL_MIN )
[ + - ][ - + ]
1502 : 0 : metric_vals->warpage = VERDICT_DBL_MAX;
1503 : : else
1504 : : {
1505 : : metric_vals->warpage =
1506 [ + - ][ + - ]: 8 : pow( VERDICT_MIN( corner_normals[0] % corner_normals[2], corner_normals[1] % corner_normals[3] ),
1507 [ - + ][ # # ]: 8 : 3 );
[ + - ]
1508 : : }
1509 : : }
1510 : : }
1511 : :
1512 [ + - ]: 8 : if( metrics_request_flag & V_QUAD_STRETCH )
1513 : : {
1514 [ + - ]: 8 : VerdictVector temp;
1515 : :
1516 : 16 : temp.set( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
1517 [ + - ]: 8 : coordinates[2][2] - coordinates[0][2] );
1518 [ + - ]: 8 : double diag02 = temp.length_squared();
1519 : :
1520 : 16 : temp.set( coordinates[3][0] - coordinates[1][0], coordinates[3][1] - coordinates[1][1],
1521 [ + - ]: 8 : coordinates[3][2] - coordinates[1][2] );
1522 [ + - ]: 8 : double diag13 = temp.length_squared();
1523 : :
1524 : : static const double QUAD_STRETCH_FACTOR = sqrt( 2.0 );
1525 : :
1526 : : // 'diag02' is now the max diagonal of the quad
1527 [ - + ]: 8 : diag02 = VERDICT_MAX( diag02, diag13 );
1528 : :
1529 [ - + ]: 8 : if( diag02 < VERDICT_DBL_MIN )
1530 : 0 : metric_vals->stretch = VERDICT_DBL_MAX;
1531 : : else
1532 : : metric_vals->stretch =
1533 [ - + ][ # # ]: 8 : QUAD_STRETCH_FACTOR *
[ - + ]
1534 [ - + ][ - + ]: 8 : VERDICT_MIN( VERDICT_MIN( lengths[0], lengths[1] ), VERDICT_MIN( lengths[2], lengths[3] ) ) /
1535 : 8 : sqrt( diag02 );
1536 : : }
1537 : :
1538 [ + - ]: 8 : if( metrics_request_flag & ( V_QUAD_CONDITION | V_QUAD_SHAPE | V_QUAD_SHAPE_AND_SIZE ) )
1539 : : {
1540 : : double lengths_squared[4];
1541 [ + - ]: 8 : lengths_squared[0] = edges[0].length_squared();
1542 [ + - ]: 8 : lengths_squared[1] = edges[1].length_squared();
1543 [ + - ]: 8 : lengths_squared[2] = edges[2].length_squared();
1544 [ + - ]: 8 : lengths_squared[3] = edges[3].length_squared();
1545 : :
1546 [ + - ][ + - ]: 8 : if( areas[0] < VERDICT_DBL_MIN || areas[1] < VERDICT_DBL_MIN || areas[2] < VERDICT_DBL_MIN ||
[ + - ][ - + ]
1547 : 8 : areas[3] < VERDICT_DBL_MIN )
1548 : : {
1549 : 0 : metric_vals->condition = VERDICT_DBL_MAX;
1550 : 0 : metric_vals->shape = 0.0;
1551 : : }
1552 : : else
1553 : : {
1554 : 8 : double max_condition = 0.0, condition;
1555 : :
1556 : 8 : condition = ( lengths_squared[0] + lengths_squared[3] ) / areas[0];
1557 [ - + ]: 8 : max_condition = VERDICT_MAX( max_condition, condition );
1558 : :
1559 : 8 : condition = ( lengths_squared[1] + lengths_squared[0] ) / areas[1];
1560 [ - + ]: 8 : max_condition = VERDICT_MAX( max_condition, condition );
1561 : :
1562 : 8 : condition = ( lengths_squared[2] + lengths_squared[1] ) / areas[2];
1563 [ - + ]: 8 : max_condition = VERDICT_MAX( max_condition, condition );
1564 : :
1565 : 8 : condition = ( lengths_squared[3] + lengths_squared[2] ) / areas[3];
1566 [ - + ]: 8 : max_condition = VERDICT_MAX( max_condition, condition );
1567 : :
1568 : 8 : metric_vals->condition = 0.5 * max_condition;
1569 : 8 : metric_vals->shape = 2 / max_condition;
1570 : : }
1571 : : }
1572 : :
1573 [ + - ]: 8 : if( metrics_request_flag & V_QUAD_AREA )
1574 : : {
1575 [ + - ][ + - ]: 8 : if( metric_vals->area > 0 ) metric_vals->area = (double)VERDICT_MIN( metric_vals->area, VERDICT_DBL_MAX );
1576 [ + - ]: 8 : metric_vals->area = (double)VERDICT_MAX( metric_vals->area, -VERDICT_DBL_MAX );
1577 : : }
1578 : :
1579 [ + - ]: 8 : if( metrics_request_flag & V_QUAD_MAX_EDGE_RATIO )
1580 : : {
1581 [ + - ]: 8 : if( metric_vals->max_edge_ratio > 0 )
1582 [ + - ]: 8 : metric_vals->max_edge_ratio = (double)VERDICT_MIN( metric_vals->max_edge_ratio, VERDICT_DBL_MAX );
1583 [ + - ]: 8 : metric_vals->max_edge_ratio = (double)VERDICT_MAX( metric_vals->max_edge_ratio, -VERDICT_DBL_MAX );
1584 : : }
1585 : :
1586 [ + - ]: 8 : if( metrics_request_flag & V_QUAD_CONDITION )
1587 : : {
1588 [ + - ]: 8 : if( metric_vals->condition > 0 )
1589 [ + - ]: 8 : metric_vals->condition = (double)VERDICT_MIN( metric_vals->condition, VERDICT_DBL_MAX );
1590 [ + - ]: 8 : metric_vals->condition = (double)VERDICT_MAX( metric_vals->condition, -VERDICT_DBL_MAX );
1591 : : }
1592 : :
1593 : : // calculate distortion
1594 [ + - ]: 8 : if( metrics_request_flag & V_QUAD_DISTORTION )
1595 : : {
1596 [ + - ]: 8 : metric_vals->distortion = v_quad_distortion( num_nodes, coordinates );
1597 : :
1598 [ + - ]: 8 : if( metric_vals->distortion > 0 )
1599 [ + - ]: 8 : metric_vals->distortion = (double)VERDICT_MIN( metric_vals->distortion, VERDICT_DBL_MAX );
1600 [ + - ]: 8 : metric_vals->distortion = (double)VERDICT_MAX( metric_vals->distortion, -VERDICT_DBL_MAX );
1601 : : }
1602 : :
1603 [ + - ]: 8 : if( metrics_request_flag & V_QUAD_JACOBIAN )
1604 : : {
1605 [ + - ]: 8 : if( metric_vals->jacobian > 0 )
1606 [ + - ]: 8 : metric_vals->jacobian = (double)VERDICT_MIN( metric_vals->jacobian, VERDICT_DBL_MAX );
1607 [ + - ]: 8 : metric_vals->jacobian = (double)VERDICT_MAX( metric_vals->jacobian, -VERDICT_DBL_MAX );
1608 : : }
1609 : :
1610 [ + - ]: 8 : if( metrics_request_flag & V_QUAD_MAXIMUM_ANGLE )
1611 : : {
1612 [ + - ]: 8 : if( metric_vals->maximum_angle > 0 )
1613 [ + - ]: 8 : metric_vals->maximum_angle = (double)VERDICT_MIN( metric_vals->maximum_angle, VERDICT_DBL_MAX );
1614 [ + - ]: 8 : metric_vals->maximum_angle = (double)VERDICT_MAX( metric_vals->maximum_angle, -VERDICT_DBL_MAX );
1615 : : }
1616 : :
1617 [ + - ]: 8 : if( metrics_request_flag & V_QUAD_MINIMUM_ANGLE )
1618 : : {
1619 [ + - ]: 8 : if( metric_vals->minimum_angle > 0 )
1620 [ + - ]: 8 : metric_vals->minimum_angle = (double)VERDICT_MIN( metric_vals->minimum_angle, VERDICT_DBL_MAX );
1621 [ + - ]: 8 : metric_vals->minimum_angle = (double)VERDICT_MAX( metric_vals->minimum_angle, -VERDICT_DBL_MAX );
1622 : : }
1623 : :
1624 [ + - ]: 8 : if( metrics_request_flag & V_QUAD_ODDY )
1625 : : {
1626 [ - + ][ # # ]: 8 : if( metric_vals->oddy > 0 ) metric_vals->oddy = (double)VERDICT_MIN( metric_vals->oddy, VERDICT_DBL_MAX );
1627 [ + - ]: 8 : metric_vals->oddy = (double)VERDICT_MAX( metric_vals->oddy, -VERDICT_DBL_MAX );
1628 : : }
1629 : :
1630 [ + - ]: 8 : if( metrics_request_flag & V_QUAD_RELATIVE_SIZE_SQUARED )
1631 : : {
1632 [ + - ]: 8 : if( metric_vals->relative_size_squared > 0 )
1633 : : metric_vals->relative_size_squared =
1634 [ + - ]: 8 : (double)VERDICT_MIN( metric_vals->relative_size_squared, VERDICT_DBL_MAX );
1635 : : metric_vals->relative_size_squared =
1636 [ + - ]: 8 : (double)VERDICT_MAX( metric_vals->relative_size_squared, -VERDICT_DBL_MAX );
1637 : : }
1638 : :
1639 [ + - ]: 8 : if( metrics_request_flag & V_QUAD_SCALED_JACOBIAN )
1640 : : {
1641 [ + - ]: 8 : if( metric_vals->scaled_jacobian > 0 )
1642 [ + - ]: 8 : metric_vals->scaled_jacobian = (double)VERDICT_MIN( metric_vals->scaled_jacobian, VERDICT_DBL_MAX );
1643 [ + - ]: 8 : metric_vals->scaled_jacobian = (double)VERDICT_MAX( metric_vals->scaled_jacobian, -VERDICT_DBL_MAX );
1644 : : }
1645 : :
1646 [ + - ]: 8 : if( metrics_request_flag & V_QUAD_SHEAR )
1647 : : {
1648 [ + - ][ + - ]: 8 : if( metric_vals->shear > 0 ) metric_vals->shear = (double)VERDICT_MIN( metric_vals->shear, VERDICT_DBL_MAX );
1649 [ + - ]: 8 : metric_vals->shear = (double)VERDICT_MAX( metric_vals->shear, -VERDICT_DBL_MAX );
1650 : : }
1651 : :
1652 : : // calculate shear and size
1653 : : // reuse values from above
1654 [ + - ]: 8 : if( metrics_request_flag & V_QUAD_SHEAR_AND_SIZE )
1655 : : {
1656 : 8 : metric_vals->shear_and_size = metric_vals->shear * metric_vals->relative_size_squared;
1657 : :
1658 [ + - ]: 8 : if( metric_vals->shear_and_size > 0 )
1659 [ + - ]: 8 : metric_vals->shear_and_size = (double)VERDICT_MIN( metric_vals->shear_and_size, VERDICT_DBL_MAX );
1660 [ + - ]: 8 : metric_vals->shear_and_size = (double)VERDICT_MAX( metric_vals->shear_and_size, -VERDICT_DBL_MAX );
1661 : : }
1662 : :
1663 [ + - ]: 8 : if( metrics_request_flag & V_QUAD_SHAPE )
1664 : : {
1665 [ + - ][ + - ]: 8 : if( metric_vals->shape > 0 ) metric_vals->shape = (double)VERDICT_MIN( metric_vals->shape, VERDICT_DBL_MAX );
1666 [ + - ]: 8 : metric_vals->shape = (double)VERDICT_MAX( metric_vals->shape, -VERDICT_DBL_MAX );
1667 : : }
1668 : :
1669 : : // calculate shape and size
1670 : : // reuse values from above
1671 [ + - ]: 8 : if( metrics_request_flag & V_QUAD_SHAPE_AND_SIZE )
1672 : : {
1673 : 8 : metric_vals->shape_and_size = metric_vals->shape * metric_vals->relative_size_squared;
1674 : :
1675 [ + - ]: 8 : if( metric_vals->shape_and_size > 0 )
1676 [ + - ]: 8 : metric_vals->shape_and_size = (double)VERDICT_MIN( metric_vals->shape_and_size, VERDICT_DBL_MAX );
1677 [ + - ]: 8 : metric_vals->shape_and_size = (double)VERDICT_MAX( metric_vals->shape_and_size, -VERDICT_DBL_MAX );
1678 : : }
1679 : :
1680 [ + - ]: 8 : if( metrics_request_flag & V_QUAD_SKEW )
1681 : : {
1682 [ - + ][ # # ]: 8 : if( metric_vals->skew > 0 ) metric_vals->skew = (double)VERDICT_MIN( metric_vals->skew, VERDICT_DBL_MAX );
1683 [ + - ]: 8 : metric_vals->skew = (double)VERDICT_MAX( metric_vals->skew, -VERDICT_DBL_MAX );
1684 : : }
1685 : :
1686 [ + - ]: 8 : if( metrics_request_flag & V_QUAD_STRETCH )
1687 : : {
1688 [ + - ]: 8 : if( metric_vals->stretch > 0 )
1689 [ + - ]: 8 : metric_vals->stretch = (double)VERDICT_MIN( metric_vals->stretch, VERDICT_DBL_MAX );
1690 [ + - ]: 8 : metric_vals->stretch = (double)VERDICT_MAX( metric_vals->stretch, -VERDICT_DBL_MAX );
1691 : : }
1692 : :
1693 [ + - ]: 8 : if( metrics_request_flag & V_QUAD_TAPER )
1694 : : {
1695 [ - + ][ # # ]: 8 : if( metric_vals->taper > 0 ) metric_vals->taper = (double)VERDICT_MIN( metric_vals->taper, VERDICT_DBL_MAX );
1696 [ + - ]: 8 : metric_vals->taper = (double)VERDICT_MAX( metric_vals->taper, -VERDICT_DBL_MAX );
1697 : : }
1698 : :
1699 [ + - ]: 8 : if( metrics_request_flag & V_QUAD_WARPAGE )
1700 : : {
1701 [ + - ]: 8 : if( metric_vals->warpage > 0 )
1702 [ + - ]: 8 : metric_vals->warpage = (double)VERDICT_MIN( metric_vals->warpage, VERDICT_DBL_MAX );
1703 [ + - ]: 8 : metric_vals->warpage = (double)VERDICT_MAX( metric_vals->warpage, -VERDICT_DBL_MAX );
1704 : : }
1705 : :
1706 [ + - ]: 8 : if( metrics_request_flag & V_QUAD_MED_ASPECT_FROBENIUS )
1707 [ + - ]: 8 : metric_vals->med_aspect_frobenius = v_quad_med_aspect_frobenius( 4, coordinates );
1708 [ + - ]: 8 : if( metrics_request_flag & V_QUAD_MAX_ASPECT_FROBENIUS )
1709 [ + - ]: 8 : metric_vals->max_aspect_frobenius = v_quad_max_aspect_frobenius( 4, coordinates );
1710 [ + - ][ + - ]: 8 : if( metrics_request_flag & V_QUAD_EDGE_RATIO ) metric_vals->edge_ratio = v_quad_edge_ratio( 4, coordinates );
1711 [ + - ][ + - ]: 8 : if( metrics_request_flag & V_QUAD_ASPECT_RATIO ) metric_vals->aspect_ratio = v_quad_aspect_ratio( 4, coordinates );
1712 [ + - ][ + - ]: 8 : if( metrics_request_flag & V_QUAD_RADIUS_RATIO ) metric_vals->radius_ratio = v_quad_radius_ratio( 4, coordinates );
1713 : 8 : }
|
