1 /***********************************************************************/
2 /* Open Visualization Data Explorer */
3 /* (C) Copyright IBM Corp. 1989,1999 */
4 /* ALL RIGHTS RESERVED */
5 /* This code licensed under the */
6 /* "IBM PUBLIC LICENSE - Open Visualization Data Explorer" */
7 /***********************************************************************/
8 /*
9 * $Header: /src/master/dx/src/exec/dxmods/_normals.c,v 1.7 2006/01/04 22:00:50 davidt Exp $
10 */
11
12 #include <dxconfig.h>
13
14
15 #include <stdio.h>
16 #include <math.h>
17 #include <string.h>
18 #include <dx/dx.h>
19 #include "vectors.h"
20 #include "_normals.h"
21
22 typedef struct
23 {
24 Object obj;
25 float rad;
26 } AreaTask;
27
28 typedef struct
29 {
30 Field field;
31 float radius;
32 } AreaArgs;
33
34 #define NO_DEPENDENCY 1
35 #define DEP_ON_POSITIONS 2
36 #define DEP_ON_CONNECTIONS 3
37 #define DEP_ON_FACES 4
38
39 static Error DetermineFunction(Object, char *, int *, float *);
40 static Error _dxfNormalsTraverse(Object, char *);
41 static Error _dxfPositionNormals(Pointer);
42 static Error _dxfConnectionNormals(Pointer);
43 static Error _dxfFaceNormals(Pointer);
44 static Error _dxfAreaNormals(Pointer);
45 static int CompareDependency(Object, char *);
46 static CompositeField GridNormalsCompositeField(CompositeField, int);
47 static Field GridNormalsField(Field, int);
48 static Error _DetermineFunction(Object, char *, int *);
49 static Error _dxfConnectionNormals_Field(Field);
50 static Array _dxf_CDep_triangles(Field);
51 static Array _dxf_CDep_quads(Field);
52 static Error _dxfPositionNormals_Field(Field);
53 static Array _dxf_PDep_triangles(Field);
54 static Array _dxf_PDep_quads(Field);
55 static Error _dxfAreaNormals_Field(Field, float);
56 static float _dxfTriangleArea(float *, float *, float *);
57 static int _dxfInside(float *, float, float *);
58 static float _dxfDistance(float *, float *);
59 static Error _dxfCheckTriangle(int, float, int, int *,
60 int *, int *, float *, float *, float *);
61 static float _dxfIntersect(float *, float, float *, float *);
62 static float _dxfTriangleArea(float *, float *, float *);
63
64
65 Object
_dxfNormals(Object object,char * method)66 _dxfNormals(Object object, char *method)
67 {
68 Object copy = NULL;
69
70 if (method == NULL)
71 {
72 DXSetError(ERROR_BAD_PARAMETER, "#10000", "target dependency");
73 goto error;
74 }
75
76 copy = DXCopy(object, COPY_STRUCTURE);
77 if (! copy)
78 goto error;
79
80 if (! _dxfNormalsObject(copy, method))
81 goto error;
82
83 return copy;
84
85 error:
86 DXDelete(copy);
87 return NULL;
88 }
89
90 Object
_dxfNormalsObject(Object o,char * m)91 _dxfNormalsObject(Object o, char *m)
92 {
93 int n;
94 char *gm = NULL;
95
96 /*
97 * make sure method string is in global memory
98 */
99 n = strlen(m)+1;
100 gm = (char *)DXAllocate(n);
101 if (! gm)
102 goto error;
103 strcpy(gm, m);
104 gm[n-1] = '\0';
105
106 if (! DXCreateTaskGroup())
107 goto error;
108
109 if (! _dxfNormalsTraverse(o, gm))
110 {
111 DXAbortTaskGroup();
112 goto error;
113 }
114
115 if (! DXExecuteTaskGroup() || DXGetError() != ERROR_NONE)
116 goto error;
117
118 DXFree((Pointer)gm);
119
120 return o;
121
122 error:
123 DXFree((Pointer)gm);
124 return NULL;
125 }
126
127 static Error
_dxfNormalsTraverse(Object object,char * method)128 _dxfNormalsTraverse(Object object, char *method)
129 {
130 Class class;
131 int i;
132 Object child;
133 int dep;
134 float rad;
135 Object comp;
136
137 class = DXGetObjectClass(object);
138 if (class == CLASS_GROUP)
139 class = DXGetGroupClass((Group)object);
140
141 switch(class)
142 {
143 case CLASS_FIELD:
144 if (DXEmptyField((Field)object))
145 return OK;
146
147 if (NULL != (comp = DXGetComponentValue((Field)object, "connections")))
148 {
149 char *str;
150 Object eType = DXGetAttribute(comp, "element type");
151
152 if (! eType || DXGetObjectClass(eType) != CLASS_STRING)
153 {
154 DXSetError(ERROR_DATA_INVALID,
155 "invalid or missing element type attribute");
156 return ERROR;
157 }
158
159 str = DXGetString((String)eType);
160
161 if (strcmp(str, "quads") && strcmp(str, "triangles"))
162 return OK;
163 }
164 else if (! DXGetComponentValue((Field)object, "faces"))
165 return OK;
166
167 i = CompareDependency(object, method);
168 if (i == -1)
169 goto error;
170 else if (i == 1)
171 return OK;
172
173 if (! DetermineFunction(object, method, &dep, &rad))
174 goto error;
175
176 if (dep == NO_DEPENDENCY)
177 {
178 DXSetError(ERROR_DATA_INVALID, "#10450", method);
179 goto error;
180 }
181
182 if (! GridNormalsField((Field)object, dep))
183 {
184 if (dep == DEP_ON_POSITIONS && rad <= 0.0)
185 {
186 return DXAddTask(_dxfPositionNormals,
187 (Pointer)object, 0, 1.0);
188 }
189 else if (dep == DEP_ON_POSITIONS && rad > 0.0)
190 {
191 AreaTask task;
192
193 task.obj = object;
194 task.rad = rad;
195
196 return DXAddTask(_dxfAreaNormals,
197 (Pointer)&task, sizeof(task), 1.0);
198 }
199 else if (dep == DEP_ON_CONNECTIONS)
200 {
201 return DXAddTask(_dxfConnectionNormals,(Pointer)object,0,1.0);
202 }
203 else
204 {
205 return DXAddTask(_dxfFaceNormals,(Pointer)object,0,1.0);
206 }
207 }
208
209 break;
210
211 case CLASS_COMPOSITEFIELD:
212
213 i = CompareDependency(object, method);
214 if (i == -1)
215 goto error;
216 else if (i == 1)
217 return OK;
218
219 if (! DetermineFunction(object, method, &dep, &rad))
220 goto error;
221
222 if (GridNormalsCompositeField((CompositeField)object, dep))
223 return OK;
224
225 if (dep == NO_DEPENDENCY)
226 {
227 DXSetError(ERROR_DATA_INVALID, "#10450", "method");
228 goto error;
229 }
230
231 /*
232 * positions and area require special stuff (eg. Grow) at
233 * the CompositeField level
234 */
235
236 if (dep == DEP_ON_POSITIONS && rad <= 0.0)
237 {
238 return DXAddTask(_dxfPositionNormals,(Pointer)object,0,1.0);
239 }
240 else if (dep == DEP_ON_POSITIONS && rad > 0.0)
241 {
242 AreaTask task;
243
244 task.obj = object;
245 task.rad = rad;
246
247 return DXAddTask(_dxfAreaNormals,
248 (Pointer)&task, sizeof(task), 1.0);
249 }
250
251 /*
252 * else connections or faces, and since we don't do anything
253 * special at the composite field level, so fall through to
254 * the generic group case.
255 */
256
257 case CLASS_MULTIGRID:
258 case CLASS_SERIES:
259 case CLASS_GROUP:
260
261 i = 0;
262 for (;;)
263 {
264 child = DXGetEnumeratedMember((Group)object, i++, NULL);
265 if (! child)
266 break;
267
268 if (! _dxfNormalsTraverse(child, method))
269 goto error;
270 }
271
272 break;
273
274 case CLASS_XFORM:
275
276 if (! DXGetXformInfo((Xform)object, &child, 0))
277 goto error;
278
279 if (! _dxfNormalsObject(child, method))
280 goto error;
281
282 break;
283
284 case CLASS_CLIPPED:
285
286 if (! DXGetClippedInfo((Clipped)object, &child, 0))
287 goto error;
288
289 if (! _dxfNormalsObject(child, method))
290 goto error;
291
292 break;
293
294 default:
295 break;
296 }
297
298 return OK;
299
300 error:
301 return ERROR;
302 }
303
304 static Error
DetermineFunction(Object o,char * method,int * dep,float * radius)305 DetermineFunction(Object o, char *method, int *dep, float *radius)
306 {
307 int i;
308
309 if (method == NULL || !strcmp(method, "manhattan"))
310 {
311 *radius = -1.0;
312 *dep = DEP_ON_POSITIONS;
313 return OK;
314 }
315
316 for (i=0; method[i]; i++)
317 if ((method[i] < '0' || method[i] > '9') && method[i] != '.')
318 break;
319
320 if (method[i])
321 {
322 *radius = 0;
323 return _DetermineFunction(o, method, dep);
324 }
325 else
326 {
327 if (sscanf(method, "%f", radius) != 1)
328 return ERROR;
329 *dep = DEP_ON_POSITIONS;
330 return OK;
331 }
332 }
333
334 static Error
_DetermineFunction(Object o,char * name,int * dep)335 _DetermineFunction(Object o, char *name, int *dep)
336 {
337 Class class = DXGetObjectClass(o);
338 if (class == CLASS_GROUP)
339 class = DXGetGroupClass((Group)o);
340
341 switch(class)
342 {
343 case CLASS_FIELD:
344 {
345 char *str; Array a;
346
347 if (DXEmptyField((Field)o))
348 {
349 *dep = NO_DEPENDENCY;
350 return OK;
351 }
352
353 /*
354 * make sure the component is there
355 */
356 a = (Array)DXGetComponentValue((Field)o, name);
357 if (a == NULL)
358 {
359 DXSetError(ERROR_MISSING_DATA, "#10240", name);
360 *dep = NO_DEPENDENCY;
361 return ERROR;
362 }
363
364 if (! strcmp(name, "connections"))
365 {
366 *dep = DEP_ON_CONNECTIONS;
367 return OK;
368 }
369
370 if (! strcmp(name, "faces"))
371 {
372 *dep = DEP_ON_FACES;
373 return OK;
374 }
375
376 if (! DXGetStringAttribute((Object)a, "dep", &str))
377 {
378 DXSetError(ERROR_DATA_INVALID, "#10255", name, "dep");
379 *dep = NO_DEPENDENCY;
380 return ERROR;
381 }
382
383 if (! strcmp(str, "positions"))
384 {
385 *dep = DEP_ON_POSITIONS;
386 return OK;
387 }
388 else if (! strcmp(str, "connections"))
389 {
390 *dep = DEP_ON_CONNECTIONS;
391 return OK;
392 }
393 else if (! strcmp(str, "faces"))
394 {
395 *dep = DEP_ON_CONNECTIONS;
396 return OK;
397 }
398 else
399 {
400 DXSetError(ERROR_DATA_INVALID, "#11250", name);
401 *dep = NO_DEPENDENCY;
402 return ERROR;
403 }
404 }
405
406 case CLASS_COMPOSITEFIELD:
407 {
408 Object c;
409 int i;
410
411 i = 0;
412 for (*dep = NO_DEPENDENCY; *dep == NO_DEPENDENCY; )
413 {
414 c = DXGetEnumeratedMember((Group)o, i++, NULL);
415 if (c == NULL)
416 return OK;
417
418 if (! _DetermineFunction(c, name, dep))
419 return ERROR;
420 }
421
422 return OK;
423 }
424
425 default:
426 {
427 DXSetError(ERROR_DATA_INVALID, "#11381");
428 return ERROR;
429 }
430 }
431 }
432
433 static Error
_dxfConnectionNormals(Pointer ptr)434 _dxfConnectionNormals(Pointer ptr)
435 {
436 return _dxfConnectionNormals_Field((Field)ptr);
437 }
438
439
440 #define NORMAL(nv, nd, norm) \
441 { \
442 float length; \
443 int i; \
444 \
445 norm->x = 0.0; \
446 norm->y = 0.0; \
447 norm->z = 0.0; \
448 \
449 for (i = 0; i < nv; i++) \
450 { \
451 int j = (i == nv-1) ? 0 : i + 1; \
452 norm->x += (ptrs[i][1] - ptrs[j][1]) * \
453 (ptrs[i][2] + ptrs[j][2]); \
454 norm->y += (ptrs[i][2] - ptrs[j][2]) * \
455 (ptrs[i][0] + ptrs[j][0]); \
456 norm->z += (ptrs[i][0] - ptrs[j][0]) * \
457 (ptrs[i][1] + ptrs[j][1]); \
458 } \
459 \
460 length = vector_length(norm); \
461 if (length == 0.0) \
462 { \
463 norm->x = 0.0; \
464 norm->y = 0.0; \
465 norm->z = 1.0; \
466 goto degenerate; \
467 } \
468 else \
469 { \
470 float d = 1.0 / length; \
471 vector_scale(norm, d, norm); \
472 } \
473 }
474
475 static Error
_dxfConnectionNormals_Field(Field field)476 _dxfConnectionNormals_Field(Field field)
477 {
478 Array array;
479 int rank, shape[10];
480 Type type;
481 Category cat;
482 int n;
483 char *primitive;
484 Object attr;
485
486 /*
487 * Check first whether the input field is empty.
488 */
489
490 if (DXEmptyField(field))
491 return OK;
492
493 array = (Array)DXGetComponentValue(field, "connections");
494 if (! array)
495 return OK;
496
497 DXGetArrayInfo(array, &n, NULL, NULL, NULL,NULL);
498
499 DXDeleteComponent(field, "normals");
500
501 array = (Array)DXGetComponentValue(field, "positions");
502 DXGetArrayInfo(array, NULL, &type, &cat, &rank, shape);
503
504 if (rank != 1 || shape[0] > 3 || shape[0] < 2)
505 {
506 DXSetError(ERROR_DATA_INVALID, "positions must be 2 or 3-D");
507 return ERROR;
508 }
509
510 if (shape[0] == 2)
511 {
512 float normal[3] = {0.0, 0.0, 1.0};
513 array = (Array)DXNewConstantArray(n, (Pointer)normal,
514 TYPE_FLOAT, CATEGORY_REAL, 1, 3);
515 }
516 else
517 {
518 attr = DXGetComponentAttribute(field, "connections", "element type");
519 if (!attr)
520 {
521 DXSetError(ERROR_DATA_INVALID, "#10255",
522 "connections", "element type");
523 return ERROR;
524 }
525
526 if (DXGetObjectClass(attr) != CLASS_STRING)
527 {
528 DXSetError(ERROR_DATA_INVALID, "#10200", "element type attribute");
529 return ERROR;
530 }
531
532 /*
533 * Check the topological primitive.
534 */
535 primitive = DXGetString((String)attr);
536
537 if (! strcmp(primitive, "triangles"))
538 array = _dxf_CDep_triangles(field);
539 else if (! strcmp(primitive, "quads"))
540 array = _dxf_CDep_quads(field);
541 else
542 array = NULL;
543 }
544
545 if (array)
546 {
547 /*
548 * Set the components of the output fields.
549 */
550 if (! DXSetComponentValue(field, "normals",(Object) array))
551 return ERROR;
552
553 if (! DXSetComponentAttribute(field, "normals", "dep",
554 (Object) DXNewString("connections")))
555 return ERROR;
556
557 if (! DXEndField(field))
558 return ERROR;
559 }
560
561 return OK;
562 }
563
564 static Array
_dxf_CDep_triangles(Field field)565 _dxf_CDep_triangles(Field field)
566 {
567 Array a_connections, a_positions;
568 Array a_normals = NULL;
569 int n_connections;
570 Vector *normals;
571 int i, j;
572 int degenerate_count;
573 ArrayHandle cHandle = NULL, pHandle = NULL;
574 int nd;
575
576 /*
577 * Do quick type checks to ensure that the components are correct.
578 */
579 a_connections = (Array)DXGetComponentValue(field, "connections");
580 if (!a_connections)
581 goto error;
582
583 a_positions = (Array)DXGetComponentValue(field, "positions");
584 if (!a_positions)
585 {
586 DXSetError(ERROR_MISSING_DATA, "#10240", "positions");
587 goto error;
588 }
589
590 if (! DXGetArrayInfo(a_positions, &n_connections, NULL, NULL, NULL, &nd))
591 goto error;
592
593 if (! DXGetArrayInfo(a_connections, &n_connections, NULL, NULL, NULL, NULL))
594 goto error;
595
596 a_normals = DXNewArray(TYPE_FLOAT, CATEGORY_REAL, 1, 3);
597 if (!a_normals)
598 goto error;
599
600 /*
601 * We now assume that there is at least one connection.
602 */
603
604 if (! DXTypeCheck(a_connections, TYPE_INT, CATEGORY_REAL, 1, 3))
605 goto error;
606
607 /*
608 * Setup the final array information.
609 */
610
611 if (! DXAddArrayData(a_normals, 0, n_connections,(Pointer) NULL))
612 goto error;
613
614 /*
615 * Setup and initialize the components.
616 */
617
618 if (NULL ==(normals =(Vector *)DXGetArrayData(a_normals)))
619 goto error;
620
621 if (NULL ==(cHandle = DXCreateArrayHandle(a_connections)))
622 goto error;
623
624 if (NULL ==(pHandle = DXCreateArrayHandle(a_positions)))
625 goto error;
626
627
628 /*
629 * For each connection...
630 */
631 degenerate_count = 0;
632 for(i = 0; i < n_connections; i++, normals++)
633 {
634 int *element, ebuf[3];
635 float *ptrs[3], buf[9];
636
637 element =(int *)DXCalculateArrayEntry(cHandle, i,(Pointer)ebuf);
638
639 for(j = 0; j < 3; j++)
640 ptrs[j] = (float *)DXGetArrayEntry(pHandle, element[j], buf+j*nd);
641
642 NORMAL(3, nd, normals);
643
644 continue;
645
646 degenerate:
647 degenerate_count ++;
648 }
649
650 /*
651 * Set the components of the output fields.
652 */
653 if (! DXSetComponentValue(field, "normals",(Object) a_normals))
654 goto error;
655
656 if (! DXSetComponentAttribute(field, "normals", "dep",
657 (Object) DXNewString("connections")))
658 goto error;
659
660 if (! DXEndField(field))
661 goto error;
662
663 if (degenerate_count > 0)
664 DXMessage("Number of degenerate triangles: %d", degenerate_count);
665
666 DXFreeArrayHandle(cHandle);
667 DXFreeArrayHandle(pHandle);
668
669 return a_normals;
670
671 error:
672 if (cHandle)
673 DXFreeArrayHandle(cHandle);
674
675 if (pHandle)
676 DXFreeArrayHandle(pHandle);
677
678 if (a_normals)
679 DXDelete((Object)a_normals);
680
681 return NULL;
682
683 }
684
685 static Array
_dxf_CDep_quads(Field field)686 _dxf_CDep_quads(Field field)
687 {
688 Array a_connections, a_positions;
689 Array a_normals = NULL;
690 int n_connections;
691 Vector *normals;
692 int i, nd;
693 int degenerate_count;
694 ArrayHandle cHandle = NULL, pHandle = NULL;
695
696 /*
697 * Do quick type checks to ensure that the components are correct.
698 */
699 a_connections = (Array)DXGetComponentValue(field, "connections");
700 if (!a_connections)
701 goto error;
702
703 a_positions = (Array)DXGetComponentValue(field, "positions");
704 if (!a_positions)
705 {
706 DXSetError(ERROR_MISSING_DATA, "#10240", "positions");
707 goto error;
708 }
709
710 if (! DXGetArrayInfo(a_positions, NULL, NULL, NULL, NULL, &nd))
711 goto error;
712
713 if (! DXGetArrayInfo(a_connections, &n_connections, NULL, NULL, NULL, NULL))
714 goto error;
715
716 a_normals = DXNewArray(TYPE_FLOAT, CATEGORY_REAL, 1, 3);
717 if (!a_normals)
718 goto error;
719
720 /*
721 * We now assume that there is at least one connection.
722 */
723
724 if (! DXTypeCheck(a_connections, TYPE_INT, CATEGORY_REAL, 1, 4))
725 goto error;
726
727 /*
728 * Setup the final array information.
729 */
730
731 if (! DXAddArrayData(a_normals, 0, n_connections,(Pointer) NULL))
732 goto error;
733
734 /*
735 * Setup and initialize the components.
736 */
737
738 if (NULL ==(normals =(Vector *)DXGetArrayData(a_normals)))
739 goto error;
740
741 if (NULL ==(cHandle = DXCreateArrayHandle(a_connections)))
742 goto error;
743
744 if (NULL ==(pHandle = DXCreateArrayHandle(a_positions)))
745 goto error;
746
747 /*
748 * For each connection...
749 */
750 degenerate_count = 0;
751 for(i = 0; i < n_connections; i++, normals++)
752 {
753 int *element, ebuf[4];
754 float *ptrs[4], buf[12];
755
756 element =(int *)DXCalculateArrayEntry(cHandle, i,(Pointer)ebuf);
757
758 ptrs[0] = (float *)DXGetArrayEntry(pHandle, element[0], buf+0*nd);
759 ptrs[1] = (float *)DXGetArrayEntry(pHandle, element[2], buf+1*nd);
760 ptrs[2] = (float *)DXGetArrayEntry(pHandle, element[3], buf+2*nd);
761 ptrs[3] = (float *)DXGetArrayEntry(pHandle, element[1], buf+3*nd);
762
763 NORMAL(4, nd, normals);
764
765 continue;
766
767 degenerate:
768 degenerate_count ++;
769 }
770
771 if (degenerate_count > 0)
772 DXMessage("Number of degenerate triangles: %d", degenerate_count);
773
774 DXFreeArrayHandle(cHandle);
775 DXFreeArrayHandle(pHandle);
776
777 return a_normals;
778
779 error:
780 if (cHandle)
781 DXFreeArrayHandle(cHandle);
782
783 if (pHandle)
784 DXFreeArrayHandle(pHandle);
785
786 if (a_normals)
787 DXDelete((Object)a_normals);
788
789 return NULL;
790
791 }
792
793 static Error
positionsNormalsTask(Pointer ptr)794 positionsNormalsTask(Pointer ptr)
795 {
796 Field field = (Field)ptr;
797
798 return _dxfPositionNormals_Field(field);
799 }
800
801 static Error
_dxfPositionNormals(Pointer ptr)802 _dxfPositionNormals(Pointer ptr)
803 {
804 Object object = (Object)ptr;
805 Object child;
806 int i;
807 Class class = DXGetObjectClass(object);
808
809 if (class == CLASS_GROUP)
810 {
811 if (! DXGrow(object, 1, NULL, NULL))
812 return ERROR;
813
814 if (! DXCreateTaskGroup())
815 return ERROR;
816
817 for (i = 0; ;i++)
818 {
819 child = DXGetEnumeratedMember((Group)object, i, NULL);
820 if (! child)
821 break;
822
823 if (! DXAddTask(positionsNormalsTask, (Pointer)child, 0, 1.0))
824 {
825 DXAbortTaskGroup();
826 return ERROR;
827 }
828 }
829
830 if (! DXExecuteTaskGroup() || DXGetError() != ERROR_NONE)
831 return ERROR;
832
833 if (! DXShrink(object))
834 return ERROR;
835
836 return OK;
837 }
838 else
839 return _dxfPositionNormals_Field((Field)object);
840 }
841
842 static Error
_dxfPositionNormals_Field(Field field)843 _dxfPositionNormals_Field(Field field)
844 {
845 Array array;
846 int rank, shape[10];
847 Type type;
848 Category cat;
849 int n;
850 char *primitive;
851 Object attr;
852
853 /*
854 * Check first whether the input field is empty.
855 */
856
857 if (DXEmptyField(field))
858 return OK;
859
860 array = (Array)DXGetComponentValue(field, "connections");
861 if (! array)
862 return OK;
863
864 DXGetArrayInfo(array, &n, NULL, NULL, NULL, NULL);
865 if (n <= 0)
866 return OK;
867
868 DXDeleteComponent(field, "normals");
869
870 array = (Array)DXGetComponentValue(field, "positions");
871 if (! array)
872 {
873 DXSetError(ERROR_DATA_INVALID, "#10200", "positions component");
874 return ERROR;
875 }
876
877 DXGetArrayInfo(array, &n, &type, &cat, &rank, shape);
878
879 if (rank != 1 || shape[0] < 2 || shape[0] > 3)
880 {
881 DXSetError(ERROR_DATA_INVALID, "positions must be 2 or 3-D");
882 return ERROR;
883 }
884
885 if (shape[0] == 2)
886 {
887 float normal[3] = {0.0, 0.0, 1.0};
888 array = (Array)DXNewConstantArray(n, (Pointer)normal,
889 TYPE_FLOAT, CATEGORY_REAL, 1, 3);
890 }
891 else
892 {
893 attr = DXGetComponentAttribute(field, "connections", "element type");
894 if (!attr)
895 {
896 DXSetError(ERROR_DATA_INVALID, "#10255",
897 "connections", "element type");
898 return ERROR;
899 }
900
901 if (DXGetObjectClass(attr) != CLASS_STRING)
902 {
903 DXSetError(ERROR_DATA_INVALID, "#10200", "element type attribute");
904 return ERROR;
905 }
906
907 /*
908 * Check the topological primitive.
909 */
910
911 primitive = DXGetString((String)attr);
912
913 if (! strcmp(primitive, "triangles"))
914 array =_dxf_PDep_triangles(field);
915 else if (! strcmp(primitive, "quads"))
916 array =_dxf_PDep_quads(field);
917 else
918 array = NULL;
919 }
920
921 if (array)
922 {
923 /*
924 * Set the components of the output fields.
925 */
926
927 if (DXGetComponentValue(field, "original normals"))
928 DXDeleteComponent(field, "original normals");
929
930 if (! (DXSetComponentValue(field, "normals",(Object)array) &&
931 DXChangedComponentValues(field, "normals") &&
932 DXEndField(field)))
933 return ERROR;
934 }
935
936 return OK;
937 }
938
939 static Array
_dxf_PDep_triangles(Field field)940 _dxf_PDep_triangles(Field field)
941 {
942 Array nArray;
943 Array a_connections, a_positions;
944 int n_connections, n_positions;
945 float centroid[3];
946 float dist;
947 float invdist;
948 Vector *normals, *n;
949 int i, nd;
950 int degenerate_count;
951 ArrayHandle pHandle = NULL;
952 ArrayHandle cHandle = NULL;
953
954 degenerate_count = 0;
955
956 /*
957 * Initialize outputs.
958 */
959
960 nArray = NULL;
961 if (NULL == (nArray = DXNewArray(TYPE_FLOAT, CATEGORY_REAL, 1, 3)))
962 goto error;
963
964 a_connections = (Array)DXGetComponentValue(field, "connections");
965 a_positions = (Array)DXGetComponentValue(field, "positions");
966
967 DXGetArrayInfo(a_connections, &n_connections, NULL, NULL, NULL, NULL);
968 DXGetArrayInfo(a_positions, &n_positions, NULL, NULL, NULL, &nd);
969
970 if (! DXTypeCheck(a_connections, TYPE_INT, CATEGORY_REAL, 1, 3))
971 goto error;
972
973 /*
974 * Setup the final array information.
975 */
976 if (! DXAddArrayData(nArray, 0, n_positions, (Pointer)NULL))
977 goto error;
978
979 /*
980 * Setup and initialize the components.
981 */
982
983 normals = (Vector *)DXGetArrayData(nArray);
984
985 for(i = 0, n = normals; i < n_positions; i++, n++)
986 n->x = n->y = n->z = 0.0;
987
988 /*
989 * For each connection...
990 */
991
992 pHandle = DXCreateArrayHandle(a_positions);
993 cHandle = DXCreateArrayHandle(a_connections);
994 if (! pHandle || ! cHandle)
995 goto error;
996
997 for(i = 0; i < n_connections; i++)
998 {
999 Vector enorm;
1000 int j, *element, ebuf[3];
1001 float *ptrs[3], buf[9];
1002
1003 element =(int *)DXCalculateArrayEntry(cHandle, i,(Pointer)ebuf);
1004
1005 for(j = 0; j < 3; j++)
1006 ptrs[j] = (float *)DXGetArrayEntry(pHandle, element[j], buf+j*nd);
1007
1008 NORMAL(3, nd, (&enorm));
1009
1010 /*
1011 * Determine the coordinates of the centroid of the triangle.
1012 */
1013
1014 if (nd == 3)
1015 {
1016 centroid[0] = (ptrs[0][0]+ptrs[1][0]+ptrs[2][0])*ONE_THIRD;
1017 centroid[1] = (ptrs[0][1]+ptrs[1][1]+ptrs[2][1])*ONE_THIRD;
1018 centroid[2] = (ptrs[0][2]+ptrs[1][2]+ptrs[2][2])*ONE_THIRD;
1019
1020 for (j = 0; j < 3; j++)
1021 {
1022 float x, y, z;
1023 x = ptrs[j][0] - centroid[0];
1024 y = ptrs[j][1] - centroid[1];
1025 z = ptrs[j][2] - centroid[2];
1026 dist = ABS(x)+ ABS(y)+ ABS(z);
1027 if (dist > 0.0)
1028 {
1029 invdist = 1.0 / dist;
1030 n = normals + element[j];
1031 n->x += enorm.x * invdist;
1032 n->y += enorm.y * invdist;
1033 n->z += enorm.z * invdist;
1034 }
1035 }
1036 }
1037 else
1038 {
1039 centroid[0] = (ptrs[0][0]+ptrs[1][0]+ptrs[2][0])*ONE_THIRD;
1040 centroid[1] = (ptrs[0][1]+ptrs[1][1]+ptrs[2][1])*ONE_THIRD;
1041
1042 for (j = 0; j < 3; j++)
1043 {
1044 float x, y;
1045 x = ptrs[j][0] - centroid[0];
1046 y = ptrs[j][1] - centroid[1];
1047 dist = ABS(x)+ ABS(y);
1048 if (dist > 0.0)
1049 {
1050 invdist = 1.0 / dist;
1051 n = normals + element[j];
1052 n->z += enorm.z * invdist;
1053 }
1054 }
1055 }
1056
1057 continue;
1058
1059 degenerate:
1060 degenerate_count ++;
1061 }
1062
1063 DXFreeArrayHandle(pHandle);
1064 DXFreeArrayHandle(cHandle);
1065
1066 /*
1067 * Now normalize all the normals
1068 */
1069
1070 for(i = 0; i < n_positions; i++)
1071 {
1072 if (vector_zero((Vector *)normals))
1073 {
1074 normals->x = (float)0.0;
1075 normals->y = (float)0.0;
1076 normals->z = (float)1.0;
1077 }
1078
1079 if (! vector_normalize((Vector *)normals, (Vector *)normals))
1080 goto error;
1081
1082 normals ++;
1083 }
1084
1085 if (degenerate_count > 0)
1086 DXMessage("Number of degenerate triangles: %d", degenerate_count);
1087
1088 return nArray;
1089
1090 error:
1091 DXDelete((Object)nArray);
1092 return NULL;
1093 }
1094
1095 static Array
_dxf_PDep_quads(Field field)1096 _dxf_PDep_quads(Field field)
1097 {
1098 Array nArray;
1099 Array a_connections, a_positions;
1100 int n_connections, n_positions;
1101 float centroid[3];
1102 float dist;
1103 float invdist;
1104 Vector *normals, *n;
1105 int i, nd;
1106 int degenerate_count;
1107 ArrayHandle pHandle = NULL;
1108 ArrayHandle cHandle = NULL;
1109
1110 degenerate_count = 0;
1111
1112 /*
1113 * Initialize outputs.
1114 */
1115
1116 nArray = NULL;
1117
1118 if (NULL == (nArray = DXNewArray(TYPE_FLOAT, CATEGORY_REAL, 1, 3)))
1119 goto error;
1120
1121 a_connections = (Array)DXGetComponentValue(field, "connections");
1122 a_positions = (Array)DXGetComponentValue(field, "positions");
1123
1124 DXGetArrayInfo(a_connections, &n_connections, NULL, NULL, NULL, NULL);
1125 DXGetArrayInfo(a_positions, &n_positions, NULL, NULL, NULL, &nd);
1126
1127 if (! DXTypeCheck(a_connections, TYPE_INT, CATEGORY_REAL, 1, 4))
1128 goto error;
1129
1130 /*
1131 * Setup the final array information.
1132 */
1133 if (! DXAddArrayData(nArray, 0, n_positions, (Pointer)NULL))
1134 goto error;
1135
1136 /*
1137 * Setup and initialize the components.
1138 */
1139
1140 normals = (Vector *)DXGetArrayData(nArray);
1141
1142 for(i = 0, n = normals; i < n_positions; i++, n++)
1143 n->x = n->y = n->z = 0.0;
1144
1145 pHandle = DXCreateArrayHandle(a_positions);
1146 cHandle = DXCreateArrayHandle(a_connections);
1147 if (! pHandle || ! cHandle)
1148 goto error;
1149
1150 for(i = 0; i < n_connections; i++)
1151 {
1152 Vector enorm;
1153 int j, *element, ebuf[4];
1154 float *ptrs[4], buf[12];
1155
1156 element = (int *)DXCalculateArrayEntry(cHandle, i,(Pointer)ebuf);
1157
1158 ptrs[0] = (float *)DXGetArrayEntry(pHandle, element[0], buf+0*nd);
1159 ptrs[1] = (float *)DXGetArrayEntry(pHandle, element[2], buf+1*nd);
1160 ptrs[2] = (float *)DXGetArrayEntry(pHandle, element[3], buf+2*nd);
1161 ptrs[3] = (float *)DXGetArrayEntry(pHandle, element[1], buf+3*nd);
1162
1163 NORMAL(4, nd, (&enorm));
1164
1165 /*
1166 * Determine the addresses of the normals at the
1167 * three points.
1168 */
1169 if (nd == 3)
1170 {
1171 centroid[0] = (ptrs[0][0]+ptrs[1][0]+ptrs[2][0]+ptrs[3][0])*ONE_FOURTH;
1172 centroid[1] = (ptrs[0][1]+ptrs[1][1]+ptrs[2][1]+ptrs[3][1])*ONE_FOURTH;
1173 centroid[2] = (ptrs[0][2]+ptrs[1][2]+ptrs[2][2]+ptrs[3][2])*ONE_FOURTH;
1174
1175 for (j = 0; j < 4; j++)
1176 {
1177 float x, y, z;
1178 x = ptrs[j][0] - centroid[0];
1179 y = ptrs[j][1] - centroid[1];
1180 z = ptrs[j][2] - centroid[2];
1181 dist = ABS(x)+ ABS(y)+ ABS(z);
1182 if (dist > 0.0)
1183 {
1184 invdist = 1.0 / dist;
1185 n = normals + element[j];
1186 n->x += enorm.x * invdist;
1187 n->y += enorm.y * invdist;
1188 n->z += enorm.z * invdist;
1189 }
1190 }
1191 }
1192 else
1193 {
1194 centroid[0] = (ptrs[0][0]+ptrs[1][0]+ptrs[2][0]+ptrs[3][0])*ONE_FOURTH;
1195 centroid[1] = (ptrs[0][1]+ptrs[1][1]+ptrs[2][1]+ptrs[3][1])*ONE_FOURTH;
1196
1197 n->x = 0.0;
1198 n->y = 0.0;
1199
1200 for (j = 0; j < 4; j++)
1201 {
1202 float x, y;
1203 x = ptrs[j][0] - centroid[0];
1204 y = ptrs[j][1] - centroid[1];
1205 dist = ABS(x)+ ABS(y);
1206 if (dist > 0.0)
1207 {
1208 invdist = 1.0 / dist;
1209 n = normals + element[j];
1210 n->z += enorm.z * invdist;
1211 }
1212 }
1213 }
1214
1215 continue;
1216
1217 degenerate:
1218 degenerate_count ++;
1219 }
1220
1221 DXFreeArrayHandle(pHandle);
1222 DXFreeArrayHandle(cHandle);
1223
1224 /*
1225 * Now normalize all the normals
1226 */
1227
1228 for(i = 0; i < n_positions; i++)
1229 {
1230 if (vector_zero((Vector *)normals))
1231 {
1232 normals->x = (float)0.0;
1233 normals->y = (float)0.0;
1234 normals->z = (float)1.0;
1235 }
1236
1237 if (! vector_normalize((Vector *)normals, (Vector *)normals))
1238 goto error;
1239
1240 normals ++;
1241 }
1242
1243 if (degenerate_count > 0)
1244 DXMessage("Number of degenerate triangles: %d", degenerate_count);
1245
1246 return nArray;
1247
1248 error:
1249 DXDelete((Object)nArray);
1250 return NULL;
1251 }
1252
1253 static Error
areaNormalsTask(Pointer ptr)1254 areaNormalsTask(Pointer ptr)
1255 {
1256 AreaArgs *args = (AreaArgs *)ptr;
1257
1258 return _dxfAreaNormals_Field(args->field, args->radius);
1259 }
1260
1261 static Error
_dxfAreaNormals(Pointer ptr)1262 _dxfAreaNormals(Pointer ptr)
1263 {
1264 Object object = ((AreaTask *)ptr)->obj;
1265 float radius = ((AreaTask *)ptr)->rad;
1266 int i;
1267 Object child;
1268
1269 if (DXGetObjectClass(object) == CLASS_COMPOSITEFIELD)
1270 {
1271 if (! DXGrow(object, 1, NULL, NULL))
1272 return ERROR;
1273
1274 if (! DXCreateTaskGroup())
1275 return ERROR;
1276
1277 for (i = 0; ;i++)
1278 {
1279 AreaArgs a;
1280
1281 child = DXGetEnumeratedMember((Group)object, i, NULL);
1282 if (! child)
1283 break;
1284
1285 a.radius = radius;
1286 a.field = (Field)child;
1287
1288 if (! DXAddTask(areaNormalsTask, (Pointer)&a, sizeof(a), 1.0))
1289 {
1290 DXAbortTaskGroup();
1291 return ERROR;
1292 }
1293 }
1294
1295 if (! DXExecuteTaskGroup() || DXGetError() != ERROR_NONE)
1296 return ERROR;
1297
1298 if (! DXShrink(object))
1299 return ERROR;
1300
1301 return OK;
1302 }
1303 else
1304 return _dxfAreaNormals_Field((Field)object, radius);
1305 }
1306
1307 static Error
_dxfAreaNormals_Field(Field field,float radius)1308 _dxfAreaNormals_Field(Field field, float radius)
1309 {
1310 Array triArray, nbrArray, ptArray, nrmArray;
1311 int nTriangles;
1312 int nPoints;
1313 float *pts, *normals, *facets;
1314 int *tris, *nbrs;
1315 int *pointTris=NULL;
1316 int *pointFlags=NULL;
1317 float *f;
1318 int *t;
1319 int p, q, r, i, j;
1320 float v0[3], v1[3];
1321 float totEdgeLen;
1322
1323 facets = NULL;
1324 nrmArray = NULL;
1325
1326 triArray = (Array)DXGetComponentValue(field, "connections");
1327 ptArray = (Array)DXGetComponentValue(field, "positions");
1328
1329 DXMarkTimeLocal("nbrs start");
1330
1331 nbrArray = DXNeighbors(field);
1332 if (! nbrArray)
1333 goto error;
1334
1335 nbrs = (int *)DXGetArrayData(nbrArray);
1336 pts = (float *)DXGetArrayData(ptArray);
1337
1338 if (! DXQueryOriginalSizes(field, &nPoints, NULL))
1339 DXGetArrayInfo(ptArray, &nPoints, NULL, NULL, NULL, NULL);
1340
1341 DXGetArrayInfo(triArray, &nTriangles, NULL, NULL, NULL, NULL);
1342
1343 pointTris = (int *)DXAllocate(nPoints * sizeof(int));
1344 if (! pointTris)
1345 goto error;
1346
1347 memset((int *)pointTris, -1, nPoints * sizeof(int));
1348
1349 tris = (int *)DXGetArrayData(triArray);
1350
1351 facets = (float *)DXAllocate(nTriangles*3*sizeof(float));
1352 if (! facets)
1353 goto error;
1354
1355 memset((int *)facets, 0, nTriangles*3*sizeof(float));
1356
1357 f = facets;
1358 t = tris;
1359 totEdgeLen = 0.0;
1360 for(i = 0; i < nTriangles; i++)
1361 {
1362 p = t[0];
1363 q = t[1];
1364 r = t[2];
1365
1366 totEdgeLen += _dxfDistance(pts+(3*p), pts+(3*q));
1367 totEdgeLen += _dxfDistance(pts+(3*q), pts+(3*r));
1368 totEdgeLen += _dxfDistance(pts+(3*r), pts+(3*p));
1369
1370 v0[0] = pts[3*p + 0] - pts[3*q + 0];
1371 v0[1] = pts[3*p + 1] - pts[3*q + 1];
1372 v0[2] = pts[3*p + 2] - pts[3*q + 2];
1373
1374 v1[0] = pts[3*p + 0] - pts[3*r + 0];
1375 v1[1] = pts[3*p + 1] - pts[3*r + 1];
1376 v1[2] = pts[3*p + 2] - pts[3*r + 2];
1377
1378 f[0] = v0[1] * v1[2] - v1[1] * v0[2];
1379 f[1] = v0[2] * v1[0] - v1[2] * v0[0];
1380 f[2] = v0[0] * v1[1] - v1[0] * v0[1];
1381
1382 vector_normalize((Vector *)f, (Vector *)f);
1383
1384 f += 3;
1385
1386 for (j = 0; j < 3; j++)
1387 {
1388 if (*t < nPoints)
1389 pointTris[*t] = i;
1390 t++;
1391 }
1392 }
1393
1394 radius *= (totEdgeLen /(3*nTriangles));
1395
1396 nrmArray = DXNewArray(TYPE_FLOAT, CATEGORY_REAL, 1, 3);
1397 if (! nrmArray)
1398 goto error;
1399
1400 if (! DXAddArrayData(nrmArray, 0, nPoints, NULL))
1401 goto error;
1402
1403 normals = (float *)DXGetArrayData(nrmArray);
1404
1405 memset((int *)normals, 0, nPoints*3*sizeof(float));
1406
1407 pointFlags = (int *)DXAllocate(nTriangles * sizeof(int));
1408 if (! pointFlags)
1409 goto error;
1410
1411 memset((int *)pointFlags, -1, nTriangles * sizeof(int));
1412
1413 for(i = 0; i < nPoints; i++)
1414 {
1415 if (pointTris[i] >= 0)
1416 {
1417 _dxfCheckTriangle(i, radius, pointTris[i], pointFlags,
1418 tris, nbrs, pts, facets, normals);
1419 if (! vector_normalize((Vector *)(normals+3*i),
1420 (Vector *)(normals+3*i)))
1421 {
1422 normals[3*i+0] = 0.0;
1423 normals[3*i+1] = 1.0;
1424 normals[3*i+2] = 0.0;
1425 }
1426 }
1427 else
1428 {
1429 normals[3*i+0] = 0.0;
1430 normals[3*i+1] = 1.0;
1431 normals[3*i+2] = 0.0;
1432 }
1433 }
1434
1435 DXFree((Pointer)facets);
1436 DXFree((Pointer)pointTris);
1437 DXFree((Pointer)pointFlags);
1438
1439 /*
1440 * Set the components of the output fields.
1441 */
1442
1443 if (DXGetComponentValue(field, "original normals"))
1444 DXDeleteComponent(field, "original normals");
1445
1446 if (! (DXSetComponentValue(field, "normals",(Object)nrmArray) &&
1447 DXSetComponentAttribute(field, "normals", "dep",
1448 (Object)DXNewString("positions")) &&
1449 DXChangedComponentValues(field, "normals") &&
1450 DXEndField(field)))
1451 {
1452 goto error;
1453 }
1454
1455 return OK;
1456
1457 error:
1458 DXFree((Pointer)facets);
1459 DXFree((Pointer)pointTris);
1460 DXFree((Pointer)pointFlags);
1461 DXDelete((Object)nrmArray);
1462
1463 return ERROR;
1464 }
1465
1466 static Error
_dxfCheckTriangle(int centerIndex,float radius,int triIndex,int * flags,int * triangles,int * neighbors,float * positions,float * facets,float * normals)1467 _dxfCheckTriangle(int centerIndex, float radius, int triIndex, int *flags,
1468 int *triangles, int *neighbors, float *positions,
1469 float *facets, float *normals)
1470 {
1471 float *p, *q, *r;
1472 int *triangle;
1473 int *nbrs;
1474 float *normal;
1475 int in[3];
1476 float trap[12];
1477 float *center;
1478 float *facet;
1479 int i, j, next;
1480 int allOut, allIn;
1481 float area, t;
1482
1483 if (flags[triIndex] == centerIndex)
1484 return OK;
1485
1486 flags[triIndex] = centerIndex;
1487
1488 center = positions +(3*centerIndex);
1489 normal = normals +(3*centerIndex);
1490 facet = facets +(3*triIndex);
1491 triangle = triangles +(3*triIndex);
1492 nbrs = neighbors +(3*triIndex);
1493
1494 allIn = 1;
1495 allOut = 1;
1496 for(i = 0; i < 3; i++)
1497 {
1498 p = positions +(3*triangle[i]);
1499 in[i] = _dxfInside(center, radius, p);
1500
1501 if (in[i])
1502 allOut = 0;
1503 else
1504 allIn = 0;
1505 }
1506
1507 if (allOut)
1508 return OK;
1509
1510 if (allIn)
1511 {
1512 p = positions + 3*triangle[0];
1513 q = positions + 3*triangle[1];
1514 r = positions + 3*triangle[2];
1515
1516 area = _dxfTriangleArea(p, q, r);
1517 facet = facets + 3*(triIndex);
1518
1519 normal[0] += area*facet[0];
1520 normal[1] += area*facet[1];
1521 normal[2] += area*facet[2];
1522 }
1523 else
1524 {
1525 j = 0;
1526 for(i = 0; i < 3; i++)
1527 {
1528 if (i == 2)
1529 next = 0;
1530 else
1531 next = i + 1;
1532
1533 p = positions +(3*triangle[i]);
1534 q = positions +(3*triangle[next]);
1535
1536 if (in[i])
1537 {
1538 trap[j++] = p[0];
1539 trap[j++] = p[1];
1540 trap[j++] = p[2];
1541
1542 if (! in[next])
1543 {
1544 t = _dxfIntersect(center, radius, p, q);
1545 if (t == -1.0)
1546 {
1547 DXMessage("Intersection error");
1548 t = 0.0;
1549 }
1550
1551 trap[j++] = p[0] + t*(q[0]-p[0]);
1552 trap[j++] = p[1] + t*(q[1]-p[1]);
1553 trap[j++] = p[2] + t*(q[2]-p[2]);
1554 }
1555 }
1556 else
1557 {
1558 if (in[next])
1559 {
1560 t = _dxfIntersect(center, radius, p, q);
1561 if (-0.001 > t || 1.0001 < t)
1562 {
1563 DXMessage("Intersection error: %f", t);
1564 if (t < 0.5)t = 0.0;
1565 else t = 1.0;
1566 }
1567
1568 trap[j++] = p[0] + t*(q[0]-p[0]);
1569 trap[j++] = p[1] + t*(q[1]-p[1]);
1570 trap[j++] = p[2] + t*(q[2]-p[2]);
1571 }
1572 }
1573 }
1574
1575 if (j > 12)
1576 DXMessage("clip error");
1577
1578 facet = facets + 3*(triIndex);
1579
1580 if (j == 12)
1581 {
1582 area = _dxfTriangleArea(trap+0, trap+3, trap+6);
1583 normal[0] += area*facet[0];
1584 normal[1] += area*facet[1];
1585 normal[2] += area*facet[2];
1586
1587 area = _dxfTriangleArea(trap+0, trap+6, trap+9);
1588 normal[0] += area*facet[0];
1589 normal[1] += area*facet[1];
1590 normal[2] += area*facet[2];
1591 }
1592 else
1593 {
1594 area = _dxfTriangleArea(trap+0, trap+3, trap+6);
1595 normal[0] += area*facet[0];
1596 normal[1] += area*facet[1];
1597 normal[2] += area*facet[2];
1598 }
1599 }
1600
1601 for(i = 0; i < 3; i++)
1602 {
1603 if (i == 2)
1604 next = 0;
1605 else
1606 next = i + 1;
1607
1608 if (nbrs[i] >= 0 &&(in[i] || in[next]))
1609 _dxfCheckTriangle(centerIndex, radius, nbrs[i], flags,
1610 triangles, neighbors, positions, facets, normals);
1611 }
1612
1613 return OK;
1614 }
1615
1616 static int
_dxfInside(float * center,float radius,float * pt)1617 _dxfInside(float *center, float radius, float *pt)
1618 {
1619 return(_dxfDistance(center, pt)< radius);
1620 }
1621
1622 static float
_dxfDistance(float * p,float * q)1623 _dxfDistance(float *p, float *q)
1624 {
1625 float a, b, c;
1626
1627 a = *p++ - *q++;
1628 b = *p++ - *q++;
1629 c = *p++ - *q++;
1630
1631 return sqrt(a*a + b*b + c*c);
1632 }
1633
1634 static float
_dxfTriangleArea(float * p,float * q,float * r)1635 _dxfTriangleArea(float *p, float *q, float *r)
1636 {
1637 float a, b, c, S;
1638
1639 a = _dxfDistance(p, q);
1640 b = _dxfDistance(q, r);
1641 c = _dxfDistance(r, p);
1642
1643 S = (a + b + c)/ 2;
1644
1645 return sqrt(S*(S-a)*(S-b)*(S-c));
1646 }
1647
1648 static float
_dxfIntersect(float * center,float radius,float * p0,float * p1)1649 _dxfIntersect(float *center, float radius, float *p0, float *p1)
1650 {
1651 float x10, y10, z10;
1652 float x0c, y0c, z0c;
1653 float A, B, C, D;
1654 float t;
1655
1656 x10 = p1[0] - p0[0];
1657 y10 = p1[1] - p0[1];
1658 z10 = p1[2] - p0[2];
1659
1660 x0c = p0[0] - center[0];
1661 y0c = p0[1] - center[1];
1662 z0c = p0[2] - center[2];
1663
1664 A = x10*x10 + y10*y10 + z10*z10;
1665 B = 2 *(x10*x0c + y10*y0c + z10*z0c);
1666 C = x0c*x0c + y0c*y0c + z0c*z0c - radius*radius;
1667
1668 D = B*B - 4*A*C;
1669
1670 if (D >= 0.0)
1671 {
1672 D = sqrt(D);
1673
1674 t = (-B + D)/(2 * A);
1675 if (t >= -0.0001 && t <= 1.0001)
1676 return t;
1677
1678 t = (-B - D)/(2 * A);
1679 if (t >= -0.0001 && t <= 1.0001)
1680 return t;
1681 }
1682
1683 return -1.0;
1684 }
1685
1686 static int
CompareDependency(Object o,char * target)1687 CompareDependency(Object o, char *target)
1688 {
1689 Class class = DXGetObjectClass(o);
1690 if (class == CLASS_GROUP)
1691 class = DXGetGroupClass((Group)o);
1692
1693 if (class == CLASS_FIELD)
1694 {
1695 Object attr;
1696
1697 if (! DXGetComponentValue((Field)o, "normals"))
1698 return 0;
1699
1700 attr = DXGetComponentAttribute((Field)o, "normals", "dep");
1701 if (! attr)
1702 {
1703 DXSetError(ERROR_DATA_INVALID,
1704 "normals component is missing the dep attribute");
1705 return -1;
1706 }
1707
1708 if (DXGetObjectClass(attr) != CLASS_STRING)
1709 {
1710 DXSetError(ERROR_BAD_CLASS,
1711 "dep attribute must be class STRING");
1712 return -1;
1713 }
1714
1715 return !strcmp(DXGetString((String)attr), target);
1716 }
1717 else if (class == CLASS_COMPOSITEFIELD)
1718 {
1719 Object child;
1720 int i, j;
1721
1722 i = 0;
1723 while(NULL != (child = DXGetEnumeratedMember((Group)o, i++, NULL)))
1724 {
1725 j = CompareDependency(child, target);
1726 if (j != 1)
1727 return j;
1728 }
1729
1730 return 1;
1731 }
1732 else
1733 {
1734 DXSetError(ERROR_BAD_CLASS, "field or composite field");
1735 return -1;
1736 }
1737 }
1738
1739 static Error
_dxfFaceNormals(Pointer ptr)1740 _dxfFaceNormals(Pointer ptr)
1741 {
1742 Field field = (Field)ptr;
1743 Array fa, la, ea, pa, na;
1744 int *f, *l, *e;
1745 float *p;
1746 int nf, nl, ne, ndim, rank;
1747
1748 na = (Array)DXGetComponentValue(field, "normals");
1749 if (na)
1750 {
1751 char *str;
1752
1753 if (! DXGetStringAttribute((Object)na, "dep", &str))
1754 {
1755 DXSetError(ERROR_DATA_INVALID, "normals dependency attribute");
1756 goto error;
1757 }
1758
1759 if (! strcmp(str, "faces"))
1760 return OK;
1761
1762 na = NULL;
1763 }
1764
1765 fa = (Array)DXGetComponentValue(field, "faces");
1766 la = (Array)DXGetComponentValue(field, "loops");
1767 ea = (Array)DXGetComponentValue(field, "edges");
1768 pa = (Array)DXGetComponentValue(field, "positions");
1769 if (! fa || ! la || ! ea || ! pa)
1770 {
1771 DXSetError(ERROR_MISSING_DATA, "one of faces/loops/edges/positions");
1772 goto error;
1773 }
1774
1775 f = (int *)DXGetArrayData(fa);
1776 DXGetArrayInfo(fa, &nf, NULL, NULL, NULL, NULL);
1777
1778 l = (int *)DXGetArrayData(la);
1779 DXGetArrayInfo(la, &nl, NULL, NULL, NULL, NULL);
1780
1781 e = (int *)DXGetArrayData(ea);
1782 DXGetArrayInfo(ea, &ne, NULL, NULL, NULL, NULL);
1783
1784 p = (float *)DXGetArrayData(pa);
1785 DXGetArrayInfo(pa, NULL, NULL, NULL, &rank, &ndim);
1786
1787 if (rank != 1 || ndim < 2 || ndim > 3)
1788 {
1789 DXSetError(ERROR_DATA_INVALID, "data must be 2- or 3-D");
1790 goto error;
1791 }
1792 else
1793 {
1794 na = _dxf_FLE_Normals(f, nf, l, nl, e, ne, p, ndim);
1795 }
1796
1797 if (! na)
1798 goto error;
1799
1800 if (! DXSetStringAttribute((Object)na, "dep", "faces"))
1801 goto error;
1802
1803 if (! DXSetComponentValue(field, "normals", (Object)na))
1804 goto error;
1805
1806 return OK;
1807
1808 error:
1809 DXDelete((Object)na);
1810 return ERROR;
1811 }
1812
1813 /*
1814 * This is the interface called in triangulation code
1815 */
1816 Array
_dxf_FLE_Normals(int * f,int nf,int * l,int nl,int * e,int ne,float * v,int nd)1817 _dxf_FLE_Normals(int *f, int nf, int *l, int nl, int *e, int ne, float *v, int nd)
1818 {
1819 int i, face, maxloop;
1820 Point *normal;
1821 float **ptrs = NULL;
1822 Array nA = DXNewArray(TYPE_FLOAT, CATEGORY_REAL, 1, 3);
1823 if (! nA)
1824 goto error;
1825
1826 if (! DXAddArrayData(nA, 0, nf, NULL))
1827 goto error;
1828
1829 normal = (Point *)DXGetArrayData(nA);
1830
1831 /*
1832 * For each face, find the longest outer loop
1833 */
1834 maxloop = 0;
1835 for (face = 0; face < nf; face ++)
1836 {
1837 int loop = f[face];
1838 int start = l[loop];
1839 int end = (loop == nl-1) ? ne : l[loop+1];
1840 int knt = end - start;
1841
1842 if (knt > maxloop)
1843 maxloop = knt;
1844 }
1845
1846 ptrs = (float **)DXAllocate(maxloop*sizeof(float *));
1847 if (! ptrs)
1848 goto error;
1849
1850 /*
1851 * For each face...
1852 */
1853 for (face = 0; face < nf; face++, normal++)
1854 {
1855 int loop = f[face];
1856 int start = l[loop];
1857 int end = (loop == nl-1) ? ne : l[loop+1];
1858 int knt = end - start;
1859 float l;
1860
1861 for (i = 0; i < knt; i++)
1862 ptrs[i] = v + e[i+start]*nd;
1863
1864 NORMAL(knt, nd, normal);
1865
1866 l = DXLength(*normal);
1867 if (l != 0.0)
1868 *normal = DXMul(*normal, l);
1869
1870 continue;
1871
1872 degenerate:
1873 normal->x = 0.0;
1874 normal->y = 0.0;
1875 normal->z = 1.0;
1876 continue;
1877 }
1878
1879 DXFree((Pointer)ptrs);
1880
1881 return nA;
1882
1883 error:
1884 DXDelete((Object)nA);
1885 DXFree((Pointer)ptrs);
1886 return NULL;
1887 }
1888
1889 static Field
GridNormalsField(Field field,int dep)1890 GridNormalsField(Field field, int dep)
1891 {
1892 Array pArray, cArray, nArray;
1893 int pDim, cDim, n, counts[3];
1894 float x, y, z, d, *vectors[2], deltas[9], norm[3];
1895 int i, j;
1896 Object attr;
1897
1898 pArray = (Array)DXGetComponentValue(field, "positions");
1899 cArray = (Array)DXGetComponentValue(field, "connections");
1900
1901 if (! pArray || ! cArray)
1902 return NULL;
1903
1904 if (! DXQueryGridConnections(cArray, &cDim, NULL) ||
1905 ! DXQueryGridPositions(pArray, &pDim, counts, NULL, deltas))
1906 return NULL;
1907
1908 if (cDim != 2 || pDim < 2 || pDim > 3)
1909 return NULL;
1910
1911 for (i = j = 0; i < pDim; i++)
1912 if (counts[i] > 1)
1913 vectors[j++] = deltas + i*pDim;
1914
1915 if (j > 2)
1916 return NULL;
1917
1918 if (pDim == 3)
1919 {
1920 x = vectors[0][1]*vectors[1][2] - vectors[0][2]*vectors[1][1];
1921 y = vectors[0][2]*vectors[1][0] - vectors[0][0]*vectors[1][2];
1922 z = vectors[0][0]*vectors[1][1] - vectors[0][1]*vectors[1][0];
1923 }
1924 else
1925 {
1926 x = 0.0;
1927 y = 0.0;
1928 z = vectors[0][0]*vectors[1][1] - vectors[0][1]*vectors[1][0];
1929 }
1930
1931 d = x*x + y*y + z*z;
1932 if (d != 0.0)
1933 {
1934 d = 1.0 / sqrt(d);
1935
1936 norm[0] = x * d;
1937 norm[1] = y * d;
1938 norm[2] = z * d;
1939 }
1940 else norm[0] = norm[1] = norm[2] = 0.0;
1941
1942 if (dep == DEP_ON_POSITIONS)
1943 {
1944 DXGetArrayInfo(pArray, &n, NULL, NULL, NULL, NULL);
1945 attr = (Object)DXNewString("positions");
1946 }
1947 else if (dep == DEP_ON_CONNECTIONS)
1948 {
1949 DXGetArrayInfo(cArray, &n, NULL, NULL, NULL, NULL);
1950 attr = (Object)DXNewString("connections");
1951 }
1952 else
1953 return NULL;
1954
1955 if (! attr)
1956 return NULL;
1957
1958 nArray = (Array)DXNewConstantArray(n, (Pointer)&norm, TYPE_FLOAT,
1959 CATEGORY_REAL, 1, 3);
1960 if (! nArray)
1961 {
1962 DXDelete((Object)attr);
1963 return NULL;
1964 }
1965
1966 if (DXGetComponentValue(field, "normals"))
1967 DXDeleteComponent(field, "normals");
1968
1969 if (! DXSetComponentValue(field, "normals", (Object)nArray))
1970 {
1971 DXDelete((Object)nArray);
1972 DXDelete((Object)attr);
1973 return NULL;
1974 }
1975
1976 if (! DXSetComponentAttribute(field, "normals", "dep", attr))
1977 {
1978 DXDelete((Object)nArray);
1979 DXDelete((Object)attr);
1980 return NULL;
1981 }
1982
1983 return field;
1984 }
1985
1986 static CompositeField
GridNormalsCompositeField(CompositeField o,int dep)1987 GridNormalsCompositeField(CompositeField o, int dep)
1988 {
1989 int i;
1990 Object c;
1991
1992 i = 0;
1993 while(NULL != (c = DXGetEnumeratedMember((Group)o, i++, NULL)))
1994 if (! GridNormalsField((Field)c, dep))
1995 return NULL;
1996
1997 return o;
1998 }
1999
2000