/***********************************************************************/ /* Open Visualization Data Explorer */ /* (C) Copyright IBM Corp. 1989,1999 */ /* ALL RIGHTS RESERVED */ /* This code licensed under the */ /* "IBM PUBLIC LICENSE - Open Visualization Data Explorer" */ /***********************************************************************/ /* * $Header: /home/gda/dxcvs/dx/src/exec/dxmods/_maptoplane.c,v 1.3 1999/05/10 15:45:20 gda Exp $ */ #include #include #include #include #include "vectors.h" #include "cases.h" #define FUZZ 0.001 #define POSITIVE 0x001 #define EQ_ZERO 0x010 #define NEGATIVE 0x100 #define SIDE0 0 #define ON_PLANE 1 #define SIDE1 2 #define I_COUNT 3 #define P 0x1000 #define Q 0x0100 #define R 0x0010 #define S 0x0001 #define PRIME1 18913 #define PRIME2 10607 #define MULVECTOR(a,b,c) { \ (c).x = (a).x * (b); \ (c).y = (a).y * (b); \ (c).z = (a).z * (b); \ } #define ADDVECTOR(a,b,c) { \ (c).x = (a).x + (b).x; \ (c).y = (a).y + (b).y; \ (c).z = (a).z + (b).z; \ } #define ADDMULVECTOR(a, b, c, d) { \ (d).x = ((a).x * b) + (c).x; \ (d).y = ((a).y * b) + (c).y; \ (d).z = ((a).z * b) + (c).z; \ } #define DOTVECTOR(a,b) \ ((a).x*(b).x + (a).y*(b).y + (a).z*(b).z) #define PLANEDISTANCE(pt, pl) \ (DOTVECTOR((pt), (pl).normal) + (pl).offset) typedef struct { Vector normal; float offset; } Plane; static Error MapToPlaneObject(Object, Plane, Vector *); static Error MTP_Field(Pointer); static Field MTP_tetras(Field, Plane); static Field MTP_cubes(Field, Plane, Matrix); static Field MTP_cubes_reg(Field, Plane, Matrix); static Field MTP_cubes_reg_AxisAligned(Field, Plane, Matrix); static Field MTP_cubes_irreg(Field, Plane); static Object MTP_Cleanup(Object); static Error AddNormalsComponent(Field, Vector); static Error AddColorsComponent(Field); static int DoConnectionsMap(Field); static int Cmp(Key, Key); static PseudoKey Hash(Key); #define MIS_ALIGNED 2 #define ALIGNED 1 #define INDETERMINATE 0 static int Orientation(float *, float *, float *, Vector *); static Error OrientTriangles(Field, Vector *); static Error MakePlane(Vector *, Vector *, Plane *); typedef struct { int i0, i1; int index; float d; } Intercept; typedef struct { long hash; Intercept *ptr; } HashItem; typedef struct { Intercept x; int axis; int i, j; float intercept; } RegIntercept; static int RegAddIntercept(SegList *, int, int, int); static int IrregAddIntercept(HashTable, SegList *, int, int); static Error GetRegRegInfo(Field, Plane, Plane *, Matrix *); static Error MapArrays(Field, SegList *, SegList *, Matrix *); static Array MapArrays_PDep(Array, SegList *, Matrix *); static Array MapArrays_Constant(Array, int); static Array MapArrays_CDep(Array, SegList *); static Array MapArrays_CDep_Irreg(Array, SegList *); static Array MapArrays_PDep_Grid_Reg(Array, SegList *, Matrix *); static Array MapArrays_PDep_Grid_Irreg(Array, SegList *); static Array MapArrays_PDep_Irreg(Array, SegList *); static int IsRegReg(Field); static Matrix PermuteMatrix(Matrix, int *); Object _dxfMapToPlane(Object object, Vector *point, Vector *normal) { Object copy = NULL; Object result = NULL; Plane plane; if (! MakePlane(normal, point, &plane)) goto error; copy = DXCopy(object, COPY_STRUCTURE); if (! copy) goto error; if (! DXCreateTaskGroup()) goto error; if (! MapToPlaneObject(copy, plane, normal)) { DXAbortTaskGroup(); goto error; } if (! DXExecuteTaskGroup() || DXGetError() != ERROR_NONE) goto error; if (DXGetObjectClass(copy) != CLASS_FIELD) { result = MTP_Cleanup(copy); DXDelete(copy); if (DXGetError() != ERROR_NONE) goto error; } else result = copy; return result; error: DXDelete((Object)result); if (copy != object) DXDelete(copy); return NULL; } static Object MTP_Cleanup(Object object) { Class class; class = DXGetObjectClass(object); if (class == CLASS_GROUP) class = DXGetGroupClass((Group)object); switch(class) { case CLASS_FIELD: if (DXEmptyField((Field)object)) return NULL; else return DXCopy(object, COPY_STRUCTURE); case CLASS_GROUP: case CLASS_MULTIGRID: case CLASS_COMPOSITEFIELD: { Group old, new; Object child; int i, empty = 1; char *name; old = (Group)object; new = (Group)DXCopy(object, COPY_ATTRIBUTES); if (! new) return NULL; i = 0; while (NULL != (child = DXGetEnumeratedMember(old, i++, &name))) { child = MTP_Cleanup(child); if (child) { empty = 0; if (! DXSetMember(new, name, child)) { DXDelete(child); return NULL; } } else if (DXGetError() != ERROR_NONE) { DXDelete((Object)new); return NULL; } } if (empty) { DXDelete((Object)new); return NULL; } return (Object)new; } case CLASS_SERIES: { Series old, new; Object child; int i; float pos; int empty = 1; old = (Series)object; new = (Series)DXCopy(object, COPY_ATTRIBUTES); if (! new) return NULL; i = -1; while (NULL != (child = DXGetSeriesMember(old, ++i, &pos))) { child = MTP_Cleanup(child); if (! child && DXGetError() != ERROR_NONE) { DXDelete((Object)new); return NULL; } if (child) empty = 0; else child = (Object)DXEndField(DXNewField()); if (! DXSetSeriesMember(new, i, (double)pos, child)) { DXDelete(child); DXDelete((Object)new); return NULL; } } if (empty) { DXDelete((Object)new); return NULL; } return (Object)new; } case CLASS_XFORM: { Object child; Matrix matrix; Object new; DXGetXformInfo((Xform)object, &child, &matrix); child = MTP_Cleanup(child); if (child) { new = (Object)DXNewXform(child, matrix); if (! new) DXDelete((Object)child); return new; } else return NULL; } case CLASS_CLIPPED: { Object clipped; Object child; Object new; DXGetClippedInfo((Clipped)object, &child, &clipped); child = MTP_Cleanup(child); if (child) { new = (Object)DXNewClipped(child, clipped); if (! new) DXDelete((Object)child); return new; } else return NULL; } default: DXSetError(ERROR_BAD_CLASS, "#11381"); return ERROR; } } typedef struct { Plane plane; Vector normal; Field field; Matrix matrix; } MTPArgs; static Error MapToPlaneObject(Object object, Plane plane, Vector *normal) { MTPArgs args, iargs; Class class; class = DXGetObjectClass(object); if (class == CLASS_GROUP) class = DXGetGroupClass((Group)object); args.normal = iargs.normal = *normal; switch(class) { case CLASS_FIELD: { if (IsRegReg((Field)object)) GetRegRegInfo((Field)object, plane, &args.plane, &args.matrix); else args.plane = plane; args.field = (Field)object; return DXAddTask(MTP_Field, (Pointer)&args, sizeof(args), 1.0); } case CLASS_COMPOSITEFIELD: { Group group = (Group)object; Object child; int i, firstRegReg = 1; args.plane = plane; i = 0; while (NULL != (child = DXGetEnumeratedMember(group, i++, NULL))) { if (DXGetObjectClass(child) != CLASS_FIELD) { DXSetError(ERROR_BAD_CLASS, "#11381"); return ERROR; } if (! DXEmptyField((Field)child)) { if (IsRegReg((Field)child)) { if (firstRegReg) { if (! GetRegRegInfo((Field)child, plane, &iargs.plane, &iargs.matrix)) return ERROR; firstRegReg = 0; } iargs.field = (Field)child; if (!DXAddTask(MTP_Field, (Pointer)&iargs, sizeof(iargs), 1.0)) return ERROR; } else { args.field = (Field)child; if (! DXAddTask(MTP_Field, (Pointer)&args, sizeof(args), 1.0)) return ERROR; } } } return OK; } case CLASS_MULTIGRID: case CLASS_SERIES: case CLASS_GROUP: { Group group = (Group)object; Object child; int i; i = 0; while (NULL != (child = DXGetEnumeratedMember(group, i++, NULL))) if (! MapToPlaneObject(child, plane, normal)) return ERROR; return OK; } case CLASS_XFORM: { Matrix matrix; Plane iPlane; Object child; /* * (PtR + O)Pl + d * = DXPt(RPl) + (OPl + d) */ DXGetXformInfo((Xform)object, &child, &matrix); iPlane.normal.x = matrix.A[0][0]*plane.normal.x + matrix.A[0][1]*plane.normal.y + matrix.A[0][2]*plane.normal.z; iPlane.normal.y = matrix.A[1][0]*plane.normal.x + matrix.A[1][1]*plane.normal.y + matrix.A[1][2]*plane.normal.z; iPlane.normal.z = matrix.A[2][0]*plane.normal.x + matrix.A[2][1]*plane.normal.y + matrix.A[2][2]*plane.normal.z; iPlane.offset = matrix.b[0]*plane.normal.x + matrix.b[1]*plane.normal.y + matrix.b[2]*plane.normal.z + plane.offset; return MapToPlaneObject(child, iPlane, normal); } case CLASS_CLIPPED: { Object child; DXGetClippedInfo((Clipped)object, &child, NULL); return MapToPlaneObject(child, plane, normal); } default: DXSetError(ERROR_BAD_CLASS, "#11381"); return ERROR; } } static Error MTP_Field(Pointer ptr) { Plane plane; Field field; Matrix matrix; Vector normal; Object attr; Type type; Category cat; int rank, shape[32]; Array inData; int i, cDep; char *name; plane = ((MTPArgs *)ptr)->plane; field = ((MTPArgs *)ptr)->field; matrix = ((MTPArgs *)ptr)->matrix; normal = ((MTPArgs *)ptr)->normal; if (! DXInvalidateConnections((Object)field)) goto error; DXDeleteComponent(field, "invalid positions"); /* * Check first whether the input field is empty. */ if (DXEmptyField(field)) return OK; if (!DXGetComponentValue(field,"connections")) { DXSetError(ERROR_MISSING_DATA,"#10240", "connections"); goto error; } attr = DXGetComponentAttribute(field, "connections", "element type"); if (! attr) { DXSetError(ERROR_MISSING_DATA, "#10255", "connections", "element type"); goto error; } if (DXGetObjectClass(attr) != CLASS_STRING) { DXSetError(ERROR_INVALID_DATA, "#10200", "element type attribute"); goto error; } if (DXGetComponentValue(field, "neighbors")) DXDeleteComponent(field, "neighbors"); if (DXGetComponentValue(field, "data statistics")) DXDeleteComponent(field, "data statistics"); if (! strcmp(DXGetString((String)attr), "tetrahedra")) { if (! MTP_tetras(field, plane)) goto error; } else if (! strcmp(DXGetString((String)attr), "cubes")) { if (! MTP_cubes(field, plane, matrix)) goto error; } else { DXSetError(ERROR_INVALID_DATA, "#11380", DXGetString((String)attr)); goto error; } /* * Verify the orientation of the triangles */ if (! OrientTriangles(field, &normal)) goto error; if (! AddColorsComponent(field)) goto error; if (! AddNormalsComponent(field, normal)) goto error; if (! DXEndField(field)) goto error; return OK; error: return ERROR; } static Field MTP_cubes(Field field, Plane plane, Matrix matrix) { Array pA, cA; cA = (Array)DXGetComponentValue(field, "connections"); pA = (Array)DXGetComponentValue(field, "positions"); if (DXQueryGridConnections(cA, NULL, NULL) && DXQueryGridPositions(pA, NULL, NULL, NULL, NULL)) { if (! MTP_cubes_reg(field, plane, matrix)) return NULL; } else { if (! MTP_cubes_irreg(field, plane)) return NULL; } return field; } #define NEXTEDGE \ int next = start[(v == (nVerts-1)) ? 0 : v+1]; #define SKIP(i,j) \ { \ if (i == next || j == next) continue; \ } static Field MTP_tetras(Field field, Plane plane) { Array cArray, pArray; int nTetras, nPts; int *connections, *tet; float *hilo; HashTable ht; int i; int nTotTris = 0; SegList *tBuf = NULL; SegList *iBuf = NULL; SegList *mBuf = NULL; int count[3]; Vector origin, delta[3]; int nAlloc; ArrayHandle pHandle = NULL; Pointer pScratch = NULL; int size; float *pPtr; InvalidComponentHandle icHandle = NULL; icHandle = DXCreateInvalidComponentHandle((Object)field, NULL, "connections"); if (! icHandle) goto error; hilo = NULL; ht = NULL; cArray = (Array)DXGetComponentValue(field, "connections"); if (! cArray) { DXSetError(ERROR_MISSING_DATA, "#10240", "connections"); goto error; } pArray = (Array)DXGetComponentValue(field, "positions"); if (! pArray) { DXSetError(ERROR_MISSING_DATA, "#10240", "positions"); goto error; } DXGetArrayInfo (cArray, &nTetras, NULL, NULL, NULL, NULL); DXGetArrayInfo (pArray, &nPts, NULL, NULL, NULL, NULL); hilo = (float *)DXAllocateLocal(nPts*sizeof(float)); if (! hilo) { DXResetError(); hilo = (float *)DXAllocate(nPts*sizeof(float)); } if (! hilo) goto error; pHandle = DXCreateArrayHandle(pArray); if (! pHandle) goto error; size = DXGetItemSize(pArray); pScratch = DXAllocate(size); if (! pScratch) goto error; for (i = 0; i < nPts; i++) { float d; pPtr = (float *)DXIterateArray(pHandle, i, ((Pointer)pPtr), pScratch); d = pPtr[0]*plane.normal.x + pPtr[1]*plane.normal.y + pPtr[2]*plane.normal.z + plane.offset; hilo[i] = d; } DXFreeArrayHandle(pHandle); pHandle = NULL; DXFree(pScratch); pScratch = NULL; ht = DXCreateHash(sizeof(HashItem), NULL, Cmp); if (! ht) goto error; nAlloc = 0.1*nTetras; if (nAlloc < 100) nAlloc = 100; tBuf = DXNewSegList(3*sizeof(int), nAlloc, nAlloc); if (! tBuf) goto error; iBuf = DXNewSegList(sizeof(Intercept), nAlloc, nAlloc); if (! tBuf) goto error; if (DoConnectionsMap(field)) { mBuf = DXNewSegList(sizeof(int), nAlloc, nAlloc); if (! mBuf) goto error; } /* * First, generate the triangles. * * For each tetrahedron... */ tet = (int *)DXGetArrayData(cArray); for (i = 0; i < nTetras; i++) { int ip, iq, ir, is; float fp, fq, fr, fs; int v, flags, nVerts, loop, *start; int ptIndices[4], nRealVerts, vertex; int pos = 0, neg = 0; if (DXIsElementInvalid(icHandle, i)) { tet += 4; continue; } ip = *tet++; iq = *tet++; ir = *tet++; is = *tet++; flags = 0; if ((fp = hilo[ip]) > 0.0) { flags |= 0x1; pos ++; } else if (fp < 0.0) neg ++; if ((fq = hilo[iq]) > 0.0) { flags |= 0x2; pos ++; } else if (fq < 0.0) neg ++; if ((fr = hilo[ir]) > 0.0) { flags |= 0x4; pos ++; } else if (fr < 0.0) neg ++; if ((fs = hilo[is]) > 0.0) { flags |= 0x8; pos ++; } else if (fs < 0.0) neg ++; if (! pos && ! neg) continue; loop = tetras_case[flags]; if (tetras_case[flags+1] == loop) continue; nVerts = tetras_loops[loop].knt; start = tetras_edges + tetras_loops[loop].start; nRealVerts = 0; for (v = 0; v < nVerts; v++) { int p0, p1; NEXTEDGE; switch(start[v]) { case 0: if (fp == 0.0) { SKIP(1,2); p0 = ip; p1 = -1; } else if (fq == 0.0) { SKIP(3,4); p0 = -1; p1 = iq; } else { p0 = ip; p1 = iq; } break; case 1: if (fp == 0.0) { SKIP(0,2) p0 = ip; p1 = -1; } else if (fr == 0.0) { SKIP(3,5); p0 = -1; p1 = ir; } else { p0 = ip; p1 = ir; } break; case 2: if (fp == 0.0) { SKIP(0,1) p0 = ip; p1 = -1; } else if (fs == 0.0) { SKIP(4,5); p0 = -1; p1 = is; } else { p0 = ip; p1 = is; } break; case 3: if (fq == 0.0) { SKIP(0,4) p0 = iq; p1 = -1; } else if (fr == 0.0) { SKIP(1,5); p0 = -1; p1 = ir; } else { p0 = iq; p1 = ir; } break; case 4: if (fq == 0.0) { SKIP(0,3) p0 = iq; p1 = -1; } else if (fs == 0.0) { SKIP(2,5); p0 = -1; p1 = is; } else { p0 = iq; p1 = is; } break; case 5: if (fr == 0.0) { SKIP(1,3) p0 = ir; p1 = -1; } else if (fs == 0.0) { SKIP(2,4); p0 = -1; p1 = is; } else { p0 = ir; p1 = is; } break; } vertex = IrregAddIntercept(ht, iBuf, p0, p1); if (vertex == -1) goto error; ptIndices[nRealVerts++] = vertex; } for (v = 2; v < nRealVerts; v++) { int p = ptIndices[0]; int q = ptIndices[v-1]; int r = ptIndices[v]; if (p != q && p != r && q != r) { int *tri = (int *)DXNewSegListItem(tBuf); if (! tri) goto error; *tri++ = p; *tri++ = q; *tri++ = r; if (mBuf) { int *m = (int *)DXNewSegListItem(mBuf); if (! m) goto error; *m = i; } nTotTris ++; } } } /* * Now generate the points. */ if (nTotTris) { float a, b, c, d; Intercept *i; Array tArray; SegListSegment *slist; int tot; tArray = DXNewArray(TYPE_INT, CATEGORY_REAL, 1, 3); if (! tArray) goto error; if (! DXAllocateArray(tArray, nTotTris)) { DXDelete((Object)tArray); goto error; } DXInitGetNextSegListSegment(tBuf); tot = 0; while (NULL != (slist = DXGetNextSegListSegment(tBuf))) { int n = DXGetSegListSegmentItemCount(slist); if (! DXAddArrayData(tArray, tot, n, DXGetSegListSegmentPointer(slist))) { DXDelete((Object)tArray); goto error; } tot += n; } DXSetComponentValue(field, "connections", (Object)tArray); if (! DXSetComponentAttribute(field, "connections", "element type", (Object)DXNewString("triangles"))) goto error; if (! DXInitGetNextSegListItem(iBuf)) goto error; while (NULL != (i = (Intercept *)DXGetNextSegListItem(iBuf))) { float d0, d1; if (i->i1 != -1) { d0 = hilo[i->i0]; d1 = hilo[i->i1]; i->d = d0 / (d0 - d1); } } if (! MapArrays(field, iBuf, mBuf, NULL)) goto error; if (! DXEndField(field)) goto error; } else { char *name; while (DXGetEnumeratedComponentValue(field, 0, &name)) DXDeleteComponent(field, name); if (!DXEndField(field)) goto error; } ok: if (ht) DXDestroyHash(ht); DXFree((Pointer)hilo); if (iBuf) DXDeleteSegList(iBuf); if (tBuf) DXDeleteSegList(tBuf); if (mBuf) DXDeleteSegList(mBuf); if (icHandle) DXFreeInvalidComponentHandle(icHandle); return field; error: if (icHandle) DXFreeInvalidComponentHandle(icHandle); if (ht) DXDestroyHash(ht); DXFree((Pointer)hilo); if (iBuf) DXDeleteSegList(iBuf); if (tBuf) DXDeleteSegList(tBuf); if (mBuf) DXDeleteSegList(mBuf); return NULL; } static Field MTP_cubes_irreg(Field field, Plane plane) { Array cArray, pArray; int nCubes, nPts; int *cube; float *hilo; HashTable ht; int i; int nTotTris = 0; SegList *tBuf = NULL; SegList *iBuf = NULL; SegList *mBuf = NULL; int count[3]; Vector origin, delta[3]; int nAlloc; ArrayHandle cHandle = NULL, pHandle = NULL; Pointer cScratch = NULL, pScratch = NULL; int size; float *pPtr; InvalidComponentHandle icHandle = NULL; icHandle = DXCreateInvalidComponentHandle((Object)field, NULL, "connections"); if (! icHandle) goto error; hilo = NULL; ht = NULL; cArray = (Array)DXGetComponentValue(field, "connections"); if (! cArray) { DXSetError(ERROR_MISSING_DATA, "#10240", "connectionscomponent"); goto error; } pArray = (Array)DXGetComponentValue(field, "positions"); if (! pArray) { DXSetError(ERROR_MISSING_DATA, "#10240", "positions"); goto error; } DXGetArrayInfo (cArray, &nCubes, NULL, NULL, NULL, NULL); DXGetArrayInfo (pArray, &nPts, NULL, NULL, NULL, NULL); if (! hilo) { DXResetError(); hilo = (float *)DXAllocate(nPts*sizeof(float)); } if (! hilo) goto error; pHandle = DXCreateArrayHandle(pArray); if (! pHandle) goto error; size = DXGetItemSize(pArray); pScratch = DXAllocate(size); if (! pScratch) goto error; for (i = 0; i < nPts; i++) { float d; pPtr = (float *)DXIterateArray(pHandle, i, ((Pointer)pPtr), pScratch); d = pPtr[0]*plane.normal.x + pPtr[1]*plane.normal.y + pPtr[2]*plane.normal.z + plane.offset; hilo[i] = d; } DXFreeArrayHandle(pHandle); pHandle = NULL; DXFree(pScratch); pScratch = NULL; ht = DXCreateHash(sizeof(HashItem), NULL, Cmp); if (! ht) goto error; nAlloc = 0.1*nCubes; if (nAlloc < 100) nAlloc = 100; tBuf = DXNewSegList(3*sizeof(int), nAlloc, nAlloc); if (! tBuf) goto error; iBuf = DXNewSegList(sizeof(Intercept), nAlloc, nAlloc); if (! tBuf) goto error; if (DoConnectionsMap(field)) { mBuf = DXNewSegList(sizeof(int), nAlloc, nAlloc); if (! mBuf) goto error; } /* * First, generate the triangles. * * For each cube... */ cHandle = DXCreateArrayHandle(cArray); if (! cHandle) goto error; size = DXGetItemSize(cArray); cScratch = DXAllocate(size); if (! cScratch) goto error; for (i = 0; i < nCubes; i++) { int ip, iq, ir, is, it, iu, iv, iw; float fp, fq, fr, fs, ft, fu, fv, fw; int flags, nVerts, loop, *start, nLoops; int ptIndices[12], nRealVerts, vertex; int v, l; cube = (int *)DXIterateArray(cHandle, i, (Pointer)cube, cScratch); if (DXIsElementInvalid(icHandle, i)) continue; ip = cube[0]; iq = cube[1]; ir = cube[2]; is = cube[3]; it = cube[4]; iu = cube[5]; iv = cube[6]; iw = cube[7]; flags = 0; if ((fp = hilo[ip]) > 0.0) flags |= 0x01; if ((fq = hilo[iq]) > 0.0) flags |= 0x02; if ((fr = hilo[ir]) > 0.0) flags |= 0x04; if ((fs = hilo[is]) > 0.0) flags |= 0x08; if ((ft = hilo[it]) > 0.0) flags |= 0x10; if ((fu = hilo[iu]) > 0.0) flags |= 0x20; if ((fv = hilo[iv]) > 0.0) flags |= 0x40; if ((fw = hilo[iw]) > 0.0) flags |= 0x80; loop = cubes_case[flags]; nLoops = cubes_case[flags+1] - loop; for (l = 0; l < nLoops; l++) { nVerts = cubes_loops[loop].knt; start = cubes_edges + cubes_loops[loop].start; nRealVerts = 0; for (v = 0; v < nVerts; v++) { int p0, p1; NEXTEDGE; switch(start[v]) { case 0: if (fp == 0.0) { SKIP(4,8) p0 = ip; p1 = -1; } else if (fq == 0.0) { SKIP(5,9); p0 = -1; p1 = iq; } else { p0 = ip; p1 = iq; } break; case 1: if (fr == 0.0) { SKIP(4,10) p0 = ir; p1 = -1; } else if (fs == 0.0) { SKIP(5,11); p0 = -1; p1 = is; } else { p0 = ir; p1 = is; } break; case 2: if (ft == 0.0) { SKIP(6,8) p0 = it; p1 = -1; } else if (fu == 0.0) { SKIP(7,9); p0 = -1; p1 = iu; } else { p0 = it; p1 = iu; } break; case 3: if (fv == 0.0) { SKIP(6,10) p0 = iv; p1 = -1; } else if (fw == 0.0) { SKIP(7,11); p0 = -1; p1 = iw; } else { p0 = iv; p1 = iw; } break; case 4: if (fp == 0.0) { SKIP(0,8) p0 = ip; p1 = -1; } else if (fr == 0.0) { SKIP(1,10); p0 = -1; p1 = ir; } else { p0 = ip; p1 = ir; } break; case 5: if (fq == 0.0) { SKIP(0,9) p0 = iq; p1 = -1; } else if (fs == 0.0) { SKIP(1,11); p0 = -1; p1 = is; } else { p0 = iq; p1 = is; } break; case 6: if (ft == 0.0) { SKIP(2,8) p0 = it; p1 = -1; } else if (fv == 0.0) { SKIP(3,10); p0 = -1; p1 = iv; } else { p0 = it; p1 = iv; } break; case 7: if (fu == 0.0) { SKIP(2,9) p0 = iu; p1 = -1; } else if (fw == 0.0) { SKIP(3,11); p0 = -1; p1 = iw; } else { p0 = iu; p1 = iw; } break; case 8: if (fp == 0.0) { SKIP(0,4) p0 = ip; p1 = -1; } else if (ft == 0.0) { SKIP(2,6); p0 = -1; p1 = it; } else { p0 = ip; p1 = it; } break; case 9: if (fu == 0.0) { SKIP(2,7) p0 = iu; p1 = -1; } else if (fq == 0.0) { SKIP(0,5); p0 = -1; p1 = iq; } else { p0 = iu; p1 = iq; } break; case 10: if (fr == 0.0) { SKIP(1,4) p0 = ir; p1 = -1; } else if (fv == 0.0) { SKIP(3,6); p0 = -1; p1 = iv; } else { p0 = ir; p1 = iv; } break; case 11: if (fs == 0.0) { SKIP(1,5) p0 = is; p1 = -1; } else if (fw == 0.0) { SKIP(3,7); p0 = -1; p1 = iw; } else { p0 = is; p1 = iw; } break; } vertex = IrregAddIntercept(ht, iBuf, p0, p1); if (vertex == -1) goto error; ptIndices[nRealVerts++] = vertex; } for (v = 2; v < nRealVerts; v++) { int p, q, r; p = ptIndices[0]; q = ptIndices[v-1]; r = ptIndices[v]; if (p != q && p != r && q != r) { int *tri = (int *)DXNewSegListItem(tBuf); if (! tri) goto error; *tri++ = p; *tri++ = q; *tri++ = r; if (mBuf) { int *m = (int *)DXNewSegListItem(mBuf); if (! m) goto error; *m = i; } nTotTris ++; } } } } DXFreeArrayHandle(cHandle); cHandle = NULL; DXFree(cScratch); cScratch = NULL; /* * Now generate the points. */ if (nTotTris) { float a, b, c, d; Intercept *i; Array tArray; SegListSegment *slist; int tot; tArray = DXNewArray(TYPE_INT, CATEGORY_REAL, 1, 3); if (! tArray) goto error; if (! DXAllocateArray(tArray, nTotTris)) { DXDelete((Object)tArray); goto error; } DXInitGetNextSegListSegment(tBuf); tot = 0; while (NULL != (slist = DXGetNextSegListSegment(tBuf))) { int n = DXGetSegListSegmentItemCount(slist); if (! DXAddArrayData(tArray, tot, n, DXGetSegListSegmentPointer(slist))) { DXDelete((Object)tArray); goto error; } tot += n; } DXSetComponentValue(field, "connections", (Object)tArray); if (! DXSetComponentAttribute(field, "connections", "element type", (Object)DXNewString("triangles"))) goto error; if (! DXInitGetNextSegListItem(iBuf)) goto error; while (NULL != (i = (Intercept *)DXGetNextSegListItem(iBuf))) { float d0, d1; if (i->i1 != -1) { d0 = hilo[i->i0]; d1 = hilo[i->i1]; i->d = d0 / (d0 - d1); } } if (! MapArrays(field, iBuf, mBuf, NULL)) goto error; if (! DXEndField(field)) goto error; } else { char *name; while (DXGetEnumeratedComponentValue(field, 0, &name)) DXDeleteComponent(field, name); if (!DXEndField(field)) goto error; } ok: if (ht) DXDestroyHash(ht); DXFree((Pointer)hilo); if (tBuf) DXDeleteSegList(tBuf); if (iBuf) DXDeleteSegList(iBuf); if (mBuf) DXDeleteSegList(mBuf); if (cHandle) DXFreeArrayHandle(cHandle); if (pHandle) DXFreeArrayHandle(pHandle); if (icHandle) DXFreeInvalidComponentHandle(icHandle); return field; error: if (icHandle) DXFreeInvalidComponentHandle(icHandle); if (cHandle) DXFreeArrayHandle(cHandle); if (pHandle) DXFreeArrayHandle(pHandle); DXFree(pScratch); DXFree(cScratch); if (ht) DXDestroyHash(ht); DXFree((Pointer)hilo); if (iBuf) DXDeleteSegList(iBuf); if (tBuf) DXDeleteSegList(tBuf); if (mBuf) DXDeleteSegList(mBuf); return NULL; } static int IrregAddIntercept(HashTable ht, SegList *iBuf, int r, int s) { HashItem hi, *hi_ptr; Intercept i; if (r > s) { i.i0 = r; i.i1 = s; } else { i.i0 = s; i.i1 = r; } hi.hash = i.i0*PRIME1 + i.i1*PRIME2; hi.ptr = &i; if ((hi_ptr = (HashItem *)DXQueryHashElement(ht, (Key)&hi)) != NULL) return hi_ptr->ptr->index; hi.ptr = (Intercept *)DXNewSegListItem(iBuf); if (! hi.ptr) return -1; hi.ptr->i0 = i.i0; hi.ptr->i1 = i.i1; hi.ptr->index = DXGetSegListItemCount(iBuf) - 1; if (! DXInsertHashElement(ht, (Element)&hi)) return -1; return hi.ptr->index; } #define XAXIS 0 #define YAXIS 1 #define ZAXIS 2 #define LAST 0 #define NEXT 1 /* * The following cases define the action taken when the corresponding * edge is determined to intersect the plane in the regular algorithm below. */ #define REG_CASE0 \ { \ if (i == ixstart && j == jystart) \ { \ vertex = RegAddIntercept(iBuf,ZAXIS,i,j); \ if (vertex == -1) goto error; \ state[ZAXIS][LAST][j] = vertex; \ } \ else \ vertex = state[ZAXIS][LAST][j]; \ } #define REG_CASE1 \ { \ if (i == ixstart) \ { \ vertex = RegAddIntercept(iBuf,ZAXIS,i,j+1); \ if (vertex == -1) goto error; \ state[ZAXIS][LAST][j+1] = vertex; \ } \ else \ vertex = state[ZAXIS][LAST][j+1]; \ } #define REG_CASE2 \ { \ if (j == jystart) \ { \ vertex = RegAddIntercept(iBuf,ZAXIS,i+1,j); \ if (vertex == -1) goto error; \ state[ZAXIS][NEXT][j] = vertex; \ } \ else \ vertex = state[ZAXIS][NEXT][j]; \ } #define REG_CASE3 \ { \ vertex = RegAddIntercept(iBuf,ZAXIS,i+1,j+1); \ if (vertex == -1) goto error; \ state[ZAXIS][NEXT][j+1] = vertex; \ } #define REG_CASE4 \ { \ if (i == ixstart && k == kzstart) \ { \ vertex = RegAddIntercept(iBuf,YAXIS,i,k); \ if (vertex == -1) goto error; \ state[YAXIS][LAST][k] = vertex; \ } \ else \ vertex = state[YAXIS][LAST][k]; \ } #define REG_CASE5 \ { \ if (i == ixstart) \ { \ vertex = RegAddIntercept(iBuf,YAXIS,i,k+1); \ if (vertex == -1) goto error; \ state[YAXIS][LAST][k+1] = vertex; \ } \ else \ vertex = state[YAXIS][LAST][k+1]; \ } #define REG_CASE6 \ { \ if (k == kzstart) \ { \ vertex = RegAddIntercept(iBuf,YAXIS,i+1,k); \ if (vertex == -1) goto error; \ state[YAXIS][NEXT][k] = vertex; \ } \ else \ vertex = state[YAXIS][NEXT][k]; \ } #define REG_CASE7 \ { \ vertex = RegAddIntercept(iBuf,YAXIS,i+1,k+1); \ if (vertex == -1) goto error; \ state[YAXIS][NEXT][k+1] = vertex; \ } #define REG_CASE8 \ { \ if (j == jystart && k == kzstart) \ { \ vertex = RegAddIntercept(iBuf,XAXIS,j,k); \ if (vertex == -1) goto error; \ state[XAXIS][LAST][k] = vertex; \ } \ else \ vertex = state[XAXIS][LAST][k]; \ } #define REG_CASE9 \ { \ if (j == jystart) \ { \ vertex = RegAddIntercept(iBuf,XAXIS,j,k+1); \ if (vertex == -1) goto error; \ state[XAXIS][LAST][k+1] = vertex; \ } \ else \ vertex = state[XAXIS][LAST][k+1]; \ } #define REG_CASE10 \ { \ if (k == kzstart) \ { \ vertex = RegAddIntercept(iBuf,XAXIS,j+1,k); \ if (vertex == -1) goto error; \ state[XAXIS][NEXT][k] = vertex; \ } \ else \ vertex = state[XAXIS][NEXT][k]; \ } #define REG_CASE11 \ { \ vertex = RegAddIntercept(iBuf,XAXIS,j+1,k+1); \ if (vertex == -1) goto error; \ state[XAXIS][NEXT][k+1] = vertex; \ } /* * The following are used when the plane exactly hits * grid vertices. */ #define EXACT_HIT_INDEX(J,K) ((J)*count[2]+(K)) /* * The logic in the following is as follows: * * Since we are not in the axis-aligned case, we know that * the plane is not perpendicular to any one axis and therefore * is single-valued when projected onto any of the planes * perpendicular to an axis. We can therefore project the exact * match points all along the same axis. * * Two buffers are used and are swapped when the outer loop * bumps. The buffers correspond to the projection of the * last and next X-faces onto the (X,Y) plane. Since we have * the single-value property, the X-faces project onto a vertical * line in the (X,Y) plane. * * Exact matching vertices are created when the vertex is first * encountered and retrieved from the buffers when re-encountered. * The following flag = evaluates the condition deciding whether * this is the first time the vertex is encountered. The c(ijk)flgs * are set if this is the first time through the (ijk)-loops. * For example, vertex 0 is encountered for the first time only * if this is the first cube in all three directions. Vertex 1 * is encountered for the first time if this is the first * time through the i and j loops (since it is on a face shared * by a cube that precedes it in the looping) but regardless * of whether this is the first time through the k-loop. Vertex * 7 is always added, since none of the cubes it is shared by * will have been processed yet. * * * The last/next buffer is chosen by the vertex number: vertices * 0 - 3 are in the x face (eg. last) while 4 - 7 are on the x+1 face * (eg next). The actual intercept is projected along the Z axis at a * point (I,J), determined by the (i,j) location of the cube and the * vertex number. */ #define REG_EXACT(i, j, k, v) \ { \ int *buf, index; \ int ciflg = (i) == ixstart; \ int cjflg = (j) == jystart; \ int ckflg = (k) == kzstart; \ int flag = (v); \ int I, J; \ \ flag = ((v) == 0 && ciflg && cjflg && ckflg) || \ ((v) == 1 && ciflg && cjflg ) || \ ((v) == 2 && ciflg && ckflg) || \ ((v) == 3 && ciflg ) || \ ((v) == 4 && cjflg && ckflg) || \ ((v) == 5 && cjflg ) || \ ((v) == 6 && ckflg) || \ ((v) == 7 ); \ \ if ((v) > 3) buf = nextExact; \ else buf = lastExact; \ \ I = (i) + (((v) & 0x4) ? 1 : 0); \ J = (j) + (((v) & 0x2) ? 1 : 0); \ \ if (flag) \ { \ vertex = RegAddIntercept(iBuf, ZAXIS, I, J); \ buf[J] = vertex; \ } \ else vertex = buf[J]; \ } #define REG_EXACT_0 REG_EXACT(i, j, k, 0) #define REG_EXACT_1 REG_EXACT(i, j, k, 1) #define REG_EXACT_2 REG_EXACT(i, j, k, 2) #define REG_EXACT_3 REG_EXACT(i, j, k, 3) #define REG_EXACT_4 REG_EXACT(i, j, k, 4) #define REG_EXACT_5 REG_EXACT(i, j, k, 5) #define REG_EXACT_6 REG_EXACT(i, j, k, 6) #define REG_EXACT_7 REG_EXACT(i, j, k, 7) #define SNAP(t, fuzz) \ { \ int flag = t < 0; \ float d; \ if (flag) t = -t; \ d = t - (int)t; \ if (d < fuzz) \ t = (int)t; \ else if (d > (1-fuzz)) \ t = (int)t + 1; \ if (flag) t = -t; \ } #define SWAP(t, a, b) \ { \ t temp = a; \ a = b; \ b = temp; \ } static Field MTP_cubes_reg(Field field, Plane plane, Matrix matrix) { int i, j, k, l, m, n, v; Array cArray; float c00, c01, c10, c11; int slabPoints[3]; int slabCubes[3]; int count[3]; int ixstart, jystart, kzstart; int ixend, jyend, kzend; float xstart, ystart, zstart; float xend, yend, zend; int xcube, ycube, zcube; int flags; int nTotPts = 0, nTotTris = 0; int *state[3][2]; int *stateBuf = NULL; int *triBuf; int *map; int ptIndices[16]; int meshOffsets[3]; int ix, jy, kz; int ixMin, jyMin, kzMin; int ixMax, jyMax, kzMax; int ixSiz, jySiz, kzSiz; float aInv, bInv, cInv; Field new = NULL; int doit; int xproj, yproj, zproj; int planeSize; SegList *iBuf = NULL; SegList *tBuf = NULL; SegList *mBuf = NULL; float a, b, c, d; int permute[3]; float xfuzz, yfuzz, zfuzz; int cubesIgnored = 0; InvalidComponentHandle icHandle = NULL; RegIntercept *ri; int *exactHits = NULL, *lastExact, *nextExact; /* * normalize the planar equation */ d = sqrt(plane.normal.x*plane.normal.x + plane.normal.y*plane.normal.y + plane.normal.z*plane.normal.z); if (d == 0.0) { DXSetError(ERROR_BAD_PARAMETER, "degenerate planar equation"); goto error; } d = 1.0 / d; plane.normal.x *= d; plane.normal.y *= d; plane.normal.z *= d; plane.offset *= d; if (fabs(plane.normal.x) < 0.00001) plane.normal.x = 0.0; if (fabs(plane.normal.y) < 0.00001) plane.normal.y = 0.0; if (fabs(plane.normal.z) < 0.00001) plane.normal.z = 0.0; if (plane.normal.x==0.0 || plane.normal.y==0.0 || plane.normal.z==0.0) return MTP_cubes_reg_AxisAligned(field, plane, matrix); icHandle = DXCreateInvalidComponentHandle((Object)field, NULL, "connections"); if (! icHandle) goto error; triBuf = NULL; map = NULL; stateBuf = NULL; cArray = (Array)DXGetComponentValue(field, "connections"); DXQueryGridConnections(cArray, NULL, count); DXGetMeshOffsets((MeshArray)cArray, meshOffsets); slabPoints[0] = count[1]*count[2]; slabPoints[1] = count[2]; slabPoints[2] = 1; slabCubes[0] = (count[1] - 1)*(count[2] - 1); slabCubes[1] = (count[2] - 1); slabCubes[2] = 1; /* * normalize the planar equation */ d = 1.0 / sqrt(plane.normal.x*plane.normal.x + plane.normal.y*plane.normal.y + plane.normal.z*plane.normal.z); plane.normal.x *= d; plane.normal.y *= d; plane.normal.z *= d; plane.offset *= d; /* * Choose the loop ordering so that the innermost loop tends to * be as long as possible (remember that we handle each cube that * intersects the plane exactly once - this ordering minimizes the * number of loops that have to be initialized). */ permute[0] = 0; permute[1] = 1; permute[2] = 2; a = count[0]*plane.normal.x; if (a < 0) a = -a; b = count[1]*plane.normal.y; if (b < 0) b = -b; c = count[2]*plane.normal.z; if (c < 0) c = -c; if (b > a && b > c) { SWAP(float, plane.normal.x, plane.normal.y); SWAP(int, count[0], count[1]); SWAP(int, meshOffsets[0], meshOffsets[1]); SWAP(int, slabPoints[0], slabPoints[1]); SWAP(int, slabCubes[0], slabCubes[1]); SWAP(int, permute[0], permute[1]); SWAP(float, a, b); } else if (c > a && c > b) { SWAP(float, plane.normal.x, plane.normal.z); SWAP(int, count[0], count[2]); SWAP(int, meshOffsets[0], meshOffsets[2]); SWAP(int, slabPoints[0], slabPoints[2]); SWAP(int, slabCubes[0], slabCubes[2]); SWAP(int, permute[0], permute[2]); SWAP(float, a, c); } if (c > b) { SWAP(float, plane.normal.y, plane.normal.z); SWAP(int, count[1], count[2]); SWAP(int, meshOffsets[1], meshOffsets[2]); SWAP(int, slabPoints[1], slabPoints[2]); SWAP(int, slabCubes[1], slabCubes[2]); SWAP(int, permute[1], permute[2]); } matrix = PermuteMatrix(matrix, permute); if (plane.normal.x != 0) aInv = -1.0 / plane.normal.x; else aInv = 0.0; if (plane.normal.y != 0) bInv = -1.0 / plane.normal.y; else bInv = 0.0; if (plane.normal.z != 0) cInv = -1.0 / plane.normal.z; else cInv = 0.0; xfuzz = (plane.normal.x != 0.0) ? FUZZ / plane.normal.x : FUZZ; yfuzz = (plane.normal.y != 0.0) ? FUZZ / plane.normal.y : FUZZ; zfuzz = (plane.normal.z != 0.0) ? FUZZ / plane.normal.z : FUZZ; xfuzz = fabs(xfuzz); yfuzz = fabs(yfuzz); zfuzz = fabs(zfuzz); /* * An axis either intersects the plane once or is parallel to the plane. * The following flags are one if the plane intersects the axis once. */ xproj = (plane.normal.x != 0.0); yproj = (plane.normal.y != 0.0); zproj = (plane.normal.z != 0.0); /* * Get the extents of this part of the nesh. */ ixSiz = count[0] - 1; ixMin = meshOffsets[0]; ixMax = ixMin + ixSiz; jySiz = count[1] - 1; jyMin = meshOffsets[1]; jyMax = jyMin + jySiz; kzSiz = count[2] - 1; kzMin = meshOffsets[2]; kzMax = kzMin + kzSiz; stateBuf = (int *)DXAllocate((4*count[2]+2*count[1])*sizeof(int)); if (! stateBuf) goto error; state[XAXIS][LAST] = stateBuf; state[XAXIS][NEXT] = state[XAXIS][LAST] + count[2]; state[YAXIS][LAST] = state[XAXIS][NEXT] + count[2]; state[YAXIS][NEXT] = state[YAXIS][LAST] + count[2]; state[ZAXIS][LAST] = state[YAXIS][NEXT] + count[2]; state[ZAXIS][NEXT] = state[ZAXIS][LAST] + count[1]; exactHits = (int *)DXAllocate(2*count[1]*sizeof(int)); if (! exactHits) goto error; lastExact = exactHits; nextExact = lastExact + count[1]; if (count[0] >= count[1] && count[0] >= count[2]) { if (count[1] >= count[2]) planeSize = count[0]*count[1]; else planeSize = count[0]*count[2]; } else if (count[1] >= count[0] && count[1] >= count[2]) { if (count[0] >= count[2]) planeSize = count[1]*count[0]; else planeSize = count[1]*count[2]; } else { if (count[0] >= count[1]) planeSize = count[2]*count[0]; else planeSize = count[2]*count[1]; } iBuf = DXNewSegList(sizeof(RegIntercept), 0.25*planeSize, 0.25*planeSize); if (! iBuf) goto error; tBuf = DXNewSegList(3*sizeof(int), 0.5*planeSize, 0.25*planeSize); if (! tBuf) goto error; if (DoConnectionsMap(field)) { mBuf = DXNewSegList(sizeof(int), planeSize, 0.5*planeSize); if (! mBuf) goto error; } /* * Determine the slab along the X mesh axis that intersects the * plane. In the following loops, ix will index the current grid * position in the overall (composite field) mesh, while i will index * the current position in the local grid. */ if (plane.normal.x == 0.0) { xstart = ixstart = ixMin; xend = ixend = ixMax - 1; doit = 1; } else { float t, denom; if (plane.normal.y == 0.0 && plane.normal.z == 0.0) { xstart = xend = plane.offset * aInv; if (xstart == (int)xstart && plane.normal.x < 0.0) xstart = xend = xstart - 1; } else { xstart = xend = (jyMin*plane.normal.y + kzMin*plane.normal.z + plane.offset) * aInv; t = (jyMax*plane.normal.y + kzMin*plane.normal.z + plane.offset) * aInv; if (t > xend) xend = t; if (t < xstart) xstart = t; t = (jyMin*plane.normal.y + kzMax*plane.normal.z + plane.offset) * aInv; if (t > xend) xend = t; if (t < xstart) xstart = t; t = (jyMax*plane.normal.y + kzMax*plane.normal.z + plane.offset) * aInv; if (t > xend) xend = t; if (t < xstart) xstart = t; } SNAP(xend, xfuzz); SNAP(xstart, xfuzz); if (xstart == xend) { xend = ixend = xstart = ixstart = (int)xstart; } else { if (xstart < ixMin) xstart = ixMin; if (xend > ixMax) xend = ixMax; ixstart = (int)xstart; if (floor(xend) == xend) ixend = (int)(xend - 1); else ixend = (int)xend; } if (xstart >= ixMax || xend < ixMin) doit = 0; else doit = 1; } /* * Here's the X loop. For each YZ slab of the intersection area in * X, process the slab. */ if (doit) { ix = ixstart; ixstart -= ixMin; ixend -= ixMin; xcube = ixstart * slabCubes[0]; for (i = ixstart; i <= ixend; i++, ix++) { if (plane.normal.y == 0.0) { ystart = jystart = jyMin; yend = jyend = jyMax - 1; doit = 1; } else { float t; if (plane.normal.x == 0.0 && plane.normal.z == 0.0) { ystart = yend = plane.offset * bInv; if (ystart == (int)ystart && plane.normal.y < 0.0) ystart = yend = ystart - 1; } else { t = (ix*plane.normal.x + kzMin*plane.normal.z + plane.offset) * bInv; ystart = yend = t; t = (ix*plane.normal.x + kzMax*plane.normal.z + plane.offset) * bInv; if (t > yend) yend = t; if (t < ystart) ystart = t; t = ((ix+1)*plane.normal.x + kzMin*plane.normal.z + plane.offset) * bInv; if (t > yend) yend = t; if (t < ystart) ystart = t; t = ((ix+1)*plane.normal.x + kzMax*plane.normal.z + plane.offset) * bInv; if (t > yend) yend = t; if (t < ystart) ystart = t; } SNAP(yend, yfuzz); SNAP(ystart, yfuzz); if (ystart == yend) { yend = jyend = ystart = jystart = (int)ystart; } else { if (ystart < jyMin) ystart = jyMin; if (yend > jyMax) yend = jyMax; jystart = (int)ystart; if (floor(yend) == yend) jyend = (int)(yend - 1); else jyend = (int)yend; } if (ystart >= jyMax || yend < jyMin) doit = 0; else doit = 1; } /* * Here's the Y loop. For each row of a YZ slab that * intersects the plane, process the row. * In the following loops, jy will index the current grid * position in the overall (composite field) mesh, while j * will index the current position in the local grid. */ if (doit) { jy = jystart; jystart -= jyMin; jyend -= jyMin; /* * The logic depends on knowing whether this is the first * slab in the intersection volume. Because the intersection * along the slab is calculated differently from the * determination of the intersection volume, we may find no * intersections in this slab. In this case, this won't be * the first slab anymore. */ if (jystart > jyend) /* then no intersections */ ixstart ++; else { ycube = xcube + jystart*slabCubes[1]; for (j = jystart; j <= jyend; j++, jy++) { if (plane.normal.z == 0.0) { float t00, t01, t10, t11; int nVerts; int *start; int non, npos, nneg; int start1[12], nVerts1; int loops, nLoops; kzstart = 0; kzend = kzSiz - 1; t00 = ix*plane.normal.x + jy*plane.normal.y + plane.offset; SNAP(t00, zfuzz); t10 = (ix+1)*plane.normal.x + jy*plane.normal.y + plane.offset; SNAP(t10, zfuzz); t01 = ix*plane.normal.x + (jy+1)*plane.normal.y + plane.offset; SNAP(t01, zfuzz); t11 = (ix+1)*plane.normal.x + (jy+1)*plane.normal.y + plane.offset; SNAP(t11, zfuzz); non = npos = nneg = 0; if (t00 == 0.0) non ++; else if (t00 > 0) npos ++; else nneg ++; if (t10 == 0.0) non ++; else if (t10 > 0) npos ++; else nneg ++; if (t01 == 0.0) non ++; else if (t01 > 0) npos ++; else nneg ++; if (t11 == 0.0) non ++; else if (t11 > 0) npos ++; else nneg ++; if (non == 1 && (npos == 0 || nneg == 0)) { jystart ++; continue; } flags = 0; if (t00 > 0) flags |= 0x03; if (t10 > 0) flags |= 0x30; if (t01 > 0) flags |= 0x0C; if (t11 > 0) flags |= 0xC0; zcube = ycube + kzstart; loops = cubes_case[flags]; nLoops = cubes_case[flags+1] - loops; if (nLoops != 1) { DXSetError(ERROR_INTERNAL, "cube case %d", flags); goto error; } nVerts = cubes_loops[loops].knt; start = cubes_edges + cubes_loops[loops].start; for (v = 0, nVerts1 = 0; v < nVerts; v++) { NEXTEDGE; switch(start[v]) { case 4: if (t00 == 0.0) SKIP(0,8) else if (t01 == 0.0) SKIP(1,10); break; case 5: if (t00 == 0.0) SKIP(0,9) else if (t01 == 0.0) SKIP(1,11); break; case 6: if (t10 == 0.0) SKIP(2,8) else if (t11 == 0.0) SKIP(3,10); break; case 7: if (t10 == 0.0) SKIP(2,9) else if (t11 == 0.0) SKIP(3,11); break; case 8: if (t00 == 0.0) SKIP(0,4) else if (t10 == 0.0) SKIP(2,6); break; case 9: if (t00 == 0.0) SKIP(0,5) else if (t10 == 0.0) SKIP(2,7); break; case 10: if (t01 == 0.0) SKIP(1,4) else if (t11 == 0.0) SKIP(3,6); break; case 11: if (t01 == 0.0) SKIP(1,5) else if (t11 == 0.0) SKIP(3,7); break; } start1[nVerts1++] = start[v]; } for (k=kzstart; k <= kzend; k++, zcube += slabCubes[2]) { for (v = 0; v < nVerts1; v++) { int vertex; /* * Call an appropriate action for the edge. */ switch(start1[v]) { case 4: if (t00 == 0.0 && xproj) REG_CASE8 else if (t01 == 0.0 && xproj) REG_CASE10 else REG_CASE4; break; case 5: if (t00 == 0.0 && xproj) REG_CASE9 else if (t01 == 0.0 && xproj) REG_CASE11 else REG_CASE5; break; case 6: if (t11 == 0.0 && xproj) REG_CASE10 else if (t10 == 0.0 && xproj) REG_CASE8 else REG_CASE6; break; case 7: if (t11 == 0.0 && xproj) REG_CASE11 else if (t10 == 0.0 && xproj) REG_CASE9 else REG_CASE7; break; case 8: if (t00 == 0.0 && !xproj) REG_CASE4 else if (t10 == 0.0 && !xproj) REG_CASE6 else REG_CASE8; break; case 9: if (t00 == 0.0 && !xproj) REG_CASE5 else if (t01 == 0.0 && !xproj) REG_CASE7 else REG_CASE9; break; case 10: if (t01 == 0.0 && !xproj) REG_CASE4 else if (t11 == 0.0 && !xproj) REG_CASE6 else REG_CASE10; break; case 11: if (t10 == 0.0 && !xproj) REG_CASE5 else if (t11 == 0.0 && !xproj) REG_CASE7 else REG_CASE11; break; } if (vertex < 0) goto error; ptIndices[v] = vertex; } if (DXIsElementInvalid(icHandle, zcube)) { cubesIgnored ++; continue; } for (v = 2; v < nVerts; v++) { int p = ptIndices[0]; int q = ptIndices[v-1]; int r = ptIndices[v]; if (p != q && p != r && q != r) { int *tri = (int *)DXNewSegListItem(tBuf); if (! tri) goto error; *tri++ = ptIndices[0]; *tri++ = ptIndices[v-1]; *tri++ = ptIndices[v]; if (mBuf) { int *m = (int *)DXNewSegListItem(mBuf); if (! m) goto error; *m = zcube; } nTotTris ++; } } } } else { float t; if (j == jystart) { c00 = cInv * (ix*plane.normal.x + jy*plane.normal.y + plane.offset); #define ZSNAP(c) \ { \ if ((c) < kzMin) \ c = kzMin; \ else if ((c) > kzMax) \ c = kzMax; \ else \ SNAP(c, zfuzz); \ } SNAP(c00, zfuzz); c10 = cInv * ((ix+1)*plane.normal.x + jy*plane.normal.y + plane.offset); SNAP(c10, zfuzz); } else { c00 = c01; c10 = c11; } c01 = cInv * (ix*plane.normal.x + (jy+1)*plane.normal.y + plane.offset); SNAP(c01, zfuzz); c11 = cInv * ((ix+1)*plane.normal.x + (jy+1)*plane.normal.y + plane.offset); SNAP(c11, zfuzz); zstart = zend = c00; if (c10 > zend) zend = c10; if (c10 < zstart) zstart = c10; if (c01 > zend) zend = c01; if (c01 < zstart) zstart = c01; if (c11 > zend) zend = c11; if (c11 < zstart) zstart = c11; if (zstart > kzMax || zend < kzMin) continue; if (zstart == zend) { zend = kzend = zstart = kzstart = (int)zstart; } else { kzstart = (int)zstart; if (zstart < kzMin) kzstart = kzMin; if (floor(zend) == zend) kzend = (int)(zend - 1); else kzend = (int)zend; if (kzend > kzMax-1) kzend = kzMax-1; } kz = kzstart; kzstart -= kzMin; kzend -= kzMin; zcube = ycube + kzstart; flags = 0; if (kz < c00) flags |= 0x03; if (kz < c01) flags |= 0x0C; if (kz < c10) flags |= 0x30; if (kz < c11) flags |= 0xC0; /* * The logic depends on knowing whether this * is the first row in the intersection slab. * Because the intersection along the row is * calculated differently from the determination * of the intersection slab, we may find no * intersections in this row. In this case, this * won't be the first row anymore. */ if (kzstart > kzend) jystart ++; else { for (k = kzstart; k <= kzend; k++, kz++, zcube += slabCubes[2]) { int loops, nLoops; int nVerts, nRealVerts; int *start; flags = (flags >> 1) & ~0xaa; if ((kz + 1) < c00) flags |= 0x02; if ((kz + 1) < c01) flags |= 0x08; if ((kz + 1) < c10) flags |= 0x20; if ((kz + 1) < c11) flags |= 0x80; nRealVerts = 0; loops = cubes_case[flags]; nLoops = cubes_case[flags+1] - loops; if (nLoops != 1) { DXSetError(ERROR_INTERNAL, "cube case %d", flags); goto error; } nVerts = cubes_loops[loops].knt; start = cubes_edges + cubes_loops[loops].start; for (v = 0; v < nVerts; v++) { int vertex; NEXTEDGE; /* * Call an appropriate action for the * edge. If an edge endpoint lies * precisely in the plane, find an * axis that intersects the plane at * most once and project the vertex * to that plane. */ switch(start[v]) { case 0: if (kz == c00) { SKIP(4,8); REG_EXACT_0; } else if ((kz+1) == c00) { SKIP(5,9); REG_EXACT_1; } else REG_CASE0; break; case 1: if (kz == c01) { SKIP(4,10); REG_EXACT_2; } else if ((kz+1) == c01) { SKIP(5,11); REG_EXACT_3; } else REG_CASE1; break; case 2: if (kz == c10) { SKIP(8,6); REG_EXACT_4; } else if ((kz+1) == c10) { SKIP(9,7); REG_EXACT_5; } else REG_CASE2; break; case 3: if (kz == c11) { SKIP(10,6); REG_EXACT_6; } else if ((kz+1) == c11) { SKIP(11,7); REG_EXACT_7; } else REG_CASE3; break; case 4: if (kz == c00) { SKIP(0,8); REG_EXACT_0; } else if (kz == c01) { SKIP(1,10); REG_EXACT_2; } else REG_CASE4; break; case 5: if ((kz+1) == c00) { SKIP(0,9); REG_EXACT_1; } else if ((kz+1) == c01) { SKIP(1,11); REG_EXACT_3; } else REG_CASE5; break; case 6: if (kz == c10) { SKIP(2,8); REG_EXACT_4; } else if (kz == c11) { SKIP(3,10); REG_EXACT_6; } else REG_CASE6; break; case 7: if ((kz+1) == c10) { SKIP(2,9); REG_EXACT_5; } else if ((kz+1) == c11) { SKIP(3,11); REG_EXACT_7; } else REG_CASE7; break; case 8: if (kz == c00) { SKIP(0,4); REG_EXACT_0; } else if (kz == c10) { SKIP(2,6); REG_EXACT_4; } else REG_CASE8; break; case 9: if ((kz+1) == c00) { SKIP(0,5); REG_EXACT_1; } else if ((kz+1) == c10) { SKIP(2,7); REG_EXACT_5; } else REG_CASE9; break; case 10: if (kz == c01) { SKIP(1,4); REG_EXACT_2; } else if (kz == c11) { SKIP(3,6); REG_EXACT_6; } else REG_CASE10; break; case 11: if ((kz+1) == c01) { SKIP(1,5); REG_EXACT_3; } else if ((kz+1) == c11) { SKIP(3,7); REG_EXACT_7; } else REG_CASE11; break; } if (vertex < 0) goto error; ptIndices[nRealVerts++] = vertex; } if (DXIsElementInvalid(icHandle, zcube)) { cubesIgnored ++; continue; } for (v = 2; v < nRealVerts; v++) { int p = ptIndices[0]; int q = ptIndices[v-1]; int r = ptIndices[v]; if (p != q && p != r && q != r) { int *tri = (int *)DXNewSegListItem(tBuf); if (! tri) goto error; *tri++ = p; *tri++ = q; *tri++ = r; if (mBuf) { int *m = (int *) DXNewSegListItem(mBuf); if (! m) goto error; *m = zcube; } nTotTris ++; } } } } } { int *t = state[XAXIS][LAST]; state[XAXIS][LAST] = state[XAXIS][NEXT]; state[XAXIS][NEXT] = t; } ycube += slabCubes[1]; } } } { int *t = state[YAXIS][LAST]; state[YAXIS][LAST] = state[YAXIS][NEXT]; state[YAXIS][NEXT] = t; t = state[ZAXIS][LAST]; state[ZAXIS][LAST] = state[ZAXIS][NEXT]; state[ZAXIS][NEXT] = t; t = lastExact; lastExact = nextExact; nextExact = t; } xcube += slabCubes[0]; } } if (nTotTris) { float a, b, c, d; Array tArray; SegListSegment *slist; int tot; tArray = DXNewArray(TYPE_INT, CATEGORY_REAL, 1, 3); if (! tArray) goto error; if (! DXAllocateArray(tArray, nTotTris)) { DXDelete((Object)tArray); goto error; } DXInitGetNextSegListSegment(tBuf); tot = 0; while (NULL != (slist = DXGetNextSegListSegment(tBuf))) { int n = DXGetSegListSegmentItemCount(slist); if (! DXAddArrayData(tArray, tot, n, DXGetSegListSegmentPointer(slist))) { DXDelete((Object)tArray); goto error; } tot += n; } DXSetComponentValue(field, "connections", (Object)tArray); if (! DXSetComponentAttribute(field, "connections", "element type", (Object)DXNewString("triangles"))) goto error; if (! DXInitGetNextSegListItem(iBuf)) goto error; a = plane.normal.x; b = plane.normal.y; c = plane.normal.z; d = plane.offset; while (NULL != (ri = (RegIntercept *)DXGetNextSegListItem(iBuf))) { float t; int it; int I, J; I = ri->i; J = ri->j; switch(ri->axis) { case XAXIS: ri->i += meshOffsets[1]; ri->j += meshOffsets[2]; t = (b*ri->i + c*ri->j + d) * aInv; if (t < meshOffsets[0]) t = meshOffsets[0]; else if (t > (meshOffsets[0]+count[0]-1)) t = meshOffsets[0]+count[0]-1; ri->intercept = t; if ((it = floor(t)) == t && it != meshOffsets[0]) { ri->x.d = 1.0; it -= 1; } else ri->x.d = t - it; ri->x.i0 = (it-meshOffsets[0])*slabPoints[0] + I*slabPoints[1] + J*slabPoints[2]; ri->x.i1 = ri->x.i0 + slabPoints[0]; break; case YAXIS: ri->i += meshOffsets[0]; ri->j += meshOffsets[2]; t = (a*ri->i + c*ri->j + d) * bInv; if (t < meshOffsets[1]) t = meshOffsets[1]; else if (t > (meshOffsets[1]+count[1]-1)) t = meshOffsets[1]+count[1]-1; ri->intercept = t; if ((it = floor(t)) == t && it != meshOffsets[1]) { ri->x.d = 1.0; it -= 1; } else ri->x.d = t - it; ri->x.i0 = I*slabPoints[0] + (it-meshOffsets[1])*slabPoints[1] + J*slabPoints[2]; ri->x.i1 = ri->x.i0 + slabPoints[1]; break; case ZAXIS: ri->i += meshOffsets[0]; ri->j += meshOffsets[1]; t = (a*ri->i + b*ri->j + d) * cInv; if (t < meshOffsets[2]) t = meshOffsets[2]; else if (t > (meshOffsets[2]+count[2]-1)) t = meshOffsets[2]+count[2]-1; ri->intercept = t; if ((it = floor(t)) == t && it != meshOffsets[2]) { ri->x.d = 1.0; it -= 1; } else ri->x.d = t - it; ri->x.i0 = I*slabPoints[0] + J*slabPoints[1] + (it-meshOffsets[2])*slabPoints[2]; ri->x.i1 = ri->x.i0 + slabPoints[2]; break; } } if (! MapArrays(field, iBuf, mBuf, &matrix)) goto error; if (DXGetComponentValue(field, "invalid positions")) DXDeleteComponent(field, "invalid posititions"); if (DXGetComponentValue(field, "invalid connections")) DXDeleteComponent(field, "invalid connections"); if (cubesIgnored) { if (! DXInvalidateUnreferencedPositions((Object)field)) goto error; if (! DXCull((Object)field)) goto error; if (! DXEndField(field)) goto error; } } else { char *name; while (DXGetEnumeratedComponentValue(field, 0, &name)) DXDeleteComponent(field, name); if (!DXEndField(field)) goto error; } DXFree((Pointer)stateBuf); DXDeleteSegList(iBuf); DXDeleteSegList(mBuf); DXDeleteSegList(tBuf); DXFree((Pointer)exactHits); if (icHandle) DXFreeInvalidComponentHandle(icHandle); return field; error: if (icHandle) DXFreeInvalidComponentHandle(icHandle); DXFree((Pointer)stateBuf); DXFree((Pointer)exactHits); DXDeleteSegList(iBuf); DXDeleteSegList(mBuf); DXDeleteSegList(tBuf); return NULL; } static Error MapArrays(Field field, SegList *iBuf, SegList *mBuf, Matrix *mat) { int i, nIntercepts; Array old, new = NULL; Object attr; char *name; if (DXGetComponentValue(field, "data statistics")) DXDeleteComponent(field, "data statistics"); if (DXGetComponentValue(field, "box")) DXDeleteComponent(field, "box"); nIntercepts = DXGetSegListItemCount(iBuf); i = 0; while ((old = (Array)DXGetEnumeratedComponentValue(field, i++, &name)) != NULL) { if (! strcmp(name, "connections")) continue; if (NULL == (attr = DXGetComponentAttribute(field, name, "dep"))) continue; if (! strcmp(DXGetString((String)attr), "positions")) new = MapArrays_PDep(old, iBuf, mat); else new = MapArrays_CDep(old, mBuf); if (! new) goto error; if (! DXSetComponentValue(field, name, (Object)new)) goto error; new = NULL; } return OK; error: DXDelete((Object)new); return ERROR; } static Array MapArrays_PDep(Array old, SegList *iBuf, Matrix *mat) { Array new; if (DXQueryGridPositions(old, NULL, NULL, NULL, NULL)) { if (mat) { new = MapArrays_PDep_Grid_Reg(old, iBuf, mat); } else { new = MapArrays_PDep_Grid_Irreg(old, iBuf); } } else if (DXQueryConstantArray(old, NULL, NULL)) { new = MapArrays_Constant(old, DXGetSegListItemCount(iBuf)); } else { new = MapArrays_PDep_Irreg(old, iBuf); } return new; } static Array MapArrays_PDep_Grid_Reg(Array old, SegList *iBuf, Matrix *mat) { Array new; Vector *base; RegIntercept *x; float A00, A01, A02, A10, A11, A12, A20, A21, A22; float B0, B1, B2; float X, Y, Z; new = DXNewArray(TYPE_FLOAT, CATEGORY_REAL, 1, 3); if (! new) goto error; if (! DXAddArrayData(new, 0, DXGetSegListItemCount(iBuf), NULL)) goto error; base = (Vector *)DXGetArrayData(new); if (! DXInitGetNextSegListItem(iBuf)) goto error; A00 = mat->A[0][0]; A01 = mat->A[0][1]; A02 = mat->A[0][2]; A10 = mat->A[1][0]; A11 = mat->A[1][1]; A12 = mat->A[1][2]; A20 = mat->A[2][0]; A21 = mat->A[2][1]; A22 = mat->A[2][2]; B0 = mat->b[0]; B1 = mat->b[1]; B2 = mat->b[2]; while (NULL != (x = (RegIntercept *)DXGetNextSegListItem(iBuf))) { Vector v, *dst = base + x->x.index; switch(x->axis) { case XAXIS: X = x->intercept; Y = x->i; Z = x->j; break; case YAXIS: X = x->i; Y = x->intercept; Z = x->j; break; case ZAXIS: X = x->i; Y = x->j; Z = x->intercept; break; } dst->x = A00*X + A10*Y + A20*Z + B0; dst->y = A01*X + A11*Y + A21*Z + B1; dst->z = A02*X + A12*Y + A22*Z + B2; } return new; error: DXDelete((Object)new); return NULL; } static Array MapArrays_PDep_Grid_Irreg(Array old, SegList *iBuf) { Array new; Vector *base; Intercept *ix; int c[3]; Vector o, d[3]; int zknt, yknt; if (! DXQueryGridPositions(old, NULL, c, (float *)&o, (float *)d)) goto error; zknt = c[2]; yknt = c[1]; new = DXNewArray(TYPE_FLOAT, CATEGORY_REAL, 1, 3); if (! new) goto error; if (! DXAddArrayData(new, 0, DXGetSegListItemCount(iBuf), NULL)) goto error; base = (Vector *)DXGetArrayData(new); if (! DXInitGetNextSegListItem(iBuf)) goto error; while (NULL != (ix = (Intercept *)DXGetNextSegListItem(iBuf))) { Vector p0, p1, *dst = base + ix->index; int jj; float x, y, z; z = ix->i0 % zknt; jj = ix->i0 / zknt; y = jj % yknt; x = jj / yknt; ADDMULVECTOR(d[0], x, o, p0); ADDMULVECTOR(d[1], y, p0, p0); ADDMULVECTOR(d[2], z, p0, p0); if (ix->i1 != -1) { z = ix->i1 % zknt; jj = ix->i1 / zknt; y = jj % yknt; x = jj / yknt; ADDMULVECTOR(d[0], x, o, p1); ADDMULVECTOR(d[1], y, p1, p1); ADDMULVECTOR(d[2], z, p1, p1); dst->x = p0.x + ix->d*(p1.x - p0.x); dst->y = p0.y + ix->d*(p1.y - p0.y); dst->z = p0.z + ix->d*(p1.z - p0.z); } else *dst = p0; } return new; error: DXDelete((Object)new); return NULL; } static Array MapArrays_Constant(Array old, int n) { int sz = DXGetItemSize(old); Pointer a; Array new = NULL; Type t; Category c; int r; int s[32]; a = DXAllocate(sz); if (! a) goto error; DXGetArrayInfo(old, NULL, &t, &c, &r, s); DXQueryConstantArray(old, NULL, a); new = (Array)DXNewConstantArrayV(n, a, t, c, r, s); error: DXFree(a); return new; } #define INTERPOLATE_INTERCEPT(type, rnd) \ { \ type *nbase = (type *)DXGetArrayData(new); \ type *p0, *p1, *p; \ int i, itemSize; \ float d; \ \ if (!nbase) \ goto error; \ \ itemSize = DXGetItemSize(old) / DXTypeSize(t); \ \ if (! DXInitGetNextSegListItem(iBuf)) \ goto error; \ \ while (NULL != (x = (Intercept *)DXGetNextSegListItem(iBuf))) \ { \ if (x->i1 != -1) \ { \ p0 = (type *)DXGetArrayEntry(oHandle, x->i0, oScratch0); \ p1 = (type *)DXGetArrayEntry(oHandle, x->i1, oScratch1); \ p = nbase + itemSize*x->index; \ d = x->d; \ \ for (i = 0; i < itemSize; i++) \ { \ *p = (float)*p0 + d*((float)*p1 - (float)*p0) + rnd; \ p++, p0++, p1++; \ } \ } \ else \ { \ p0 = (type *)DXGetArrayEntry(oHandle, x->i0, oScratch0); \ p = nbase + itemSize*x->index; \ \ for (i = 0; i < itemSize; i++) \ { \ *p = *p0; \ p++, p0++; \ } \ } \ } \ } static Array MapArrays_PDep_Irreg(Array old, SegList *iBuf) { Array new = NULL; Intercept *x; Type t; Category c; int r, s[32], itemSize; ArrayHandle oHandle = NULL; Pointer oScratch0 = NULL, oScratch1; DXGetArrayInfo(old, NULL, &t, &c, &r, s); new = DXNewArrayV(t, c, r, s); if (! new) goto error; if (! DXAddArrayData(new, 0, DXGetSegListItemCount(iBuf), NULL)) goto error; oHandle = DXCreateArrayHandle(old); if (! oHandle) goto error; oScratch0 = DXAllocate(2*DXGetItemSize(old)); if (! oScratch0) goto error; oScratch1 = (Pointer)((char *)oScratch0 + DXGetItemSize(old)); switch(t) { case TYPE_DOUBLE: INTERPOLATE_INTERCEPT(double, 0.0); break; case TYPE_FLOAT: INTERPOLATE_INTERCEPT(float, 0.0); break; case TYPE_INT: INTERPOLATE_INTERCEPT(int, 0.5); break; case TYPE_UINT: INTERPOLATE_INTERCEPT(uint, 0.5); break; case TYPE_SHORT: INTERPOLATE_INTERCEPT(short, 0.5); break; case TYPE_USHORT: INTERPOLATE_INTERCEPT(ushort, 0.5); break; case TYPE_BYTE: INTERPOLATE_INTERCEPT(byte, 0.5); break; case TYPE_UBYTE: INTERPOLATE_INTERCEPT(ubyte, 0.5); break; default: DXSetError(ERROR_INVALID_DATA, "#10320", "data"); goto error; } DXFreeArrayHandle(oHandle); DXFree(oScratch0); return new; error: DXFreeArrayHandle(oHandle); DXFree(oScratch0); DXDelete((Object)new); return NULL; } static Array MapArrays_CDep(Array old, SegList *mBuf) { Array new; if (DXGetArrayClass(old) == CLASS_REGULARARRAY) { new = MapArrays_Constant(old, DXGetSegListItemCount(mBuf)); } else { new = MapArrays_CDep_Irreg(old, mBuf); } return new; } static Array MapArrays_CDep_Irreg(Array old, SegList *mBuf) { Array new = NULL; Type t; Category c; int r, s[32], itemSize; int i, size; char *src, *dst; SegListSegment *slist; Pointer oScratch = NULL; ArrayHandle oHandle = NULL; DXGetArrayInfo(old, NULL, &t, &c, &r, s); new = DXNewArrayV(t, c, r, s); if (! new) goto error; if (! DXAddArrayData(new, 0, DXGetSegListItemCount(mBuf), NULL)) goto error; oHandle = DXCreateArrayHandle(old); if (! oHandle) goto error; size = DXGetItemSize(old); oScratch = DXAllocate(size); if (! oScratch) goto error; dst = (char *)DXGetArrayData(new); DXInitGetNextSegListSegment(mBuf); while (NULL != (slist = DXGetNextSegListSegment(mBuf))) { int *m = (int *)DXGetSegListSegmentPointer(slist); int n = DXGetSegListSegmentItemCount(slist); for (i = 0; i < n; i++, m++) { src = (char *)DXGetArrayEntry(oHandle, (*m), oScratch); memcpy(dst, src, size); dst += size; } } DXFreeArrayHandle(oHandle); DXFree(oScratch); return new; error: DXDelete((Object)new); return NULL; } static Error AddNormalsComponent(Field field, Vector normal) { Array array; int n; Object attr; char *name; if (DXEmptyField(field)) return OK; vector_normalize(&normal, &normal); name = "colors"; array = (Array)DXGetComponentValue(field, name); if (! array) { name = "data"; array = (Array)DXGetComponentValue(field, name); } if (! array) { array = (Array)DXGetComponentValue(field, "positions"); if (! array) return OK; DXGetArrayInfo(array, &n, NULL, NULL, NULL, NULL); attr = (Object)DXNewString("positions"); } else { attr = DXGetAttribute((Object)array, "dep"); if (!attr || DXGetObjectClass(attr) != CLASS_STRING) { DXSetError(ERROR_INTERNAL, "#10241", name); return ERROR; } if (! strcmp(DXGetString((String)attr), "positions")) { array = (Array)DXGetComponentValue(field, "positions"); if (! array) return OK; DXGetArrayInfo(array, &n, NULL, NULL, NULL, NULL); } else if (! strcmp(DXGetString((String)attr), "connections")) { array = (Array)DXGetComponentValue(field, "connections"); if (! array) return OK; DXGetArrayInfo(array, &n, NULL, NULL, NULL, NULL); } else { DXSetError(ERROR_INTERNAL, "#11250", "data"); return ERROR; } } array = (Array)DXNewConstantArray(n, (Pointer)&normal, TYPE_FLOAT, CATEGORY_REAL, 1, 3); if (! array) return ERROR; if (! DXSetComponentValue(field, "normals", (Object)array)) { DXDelete((Object)array); return ERROR; } if (! DXSetComponentAttribute(field, "normals", "dep", attr)) return ERROR; return OK; } static Error AddColorsComponent(Field field) { if (DXEmptyField(field)) return OK; if (! DXGetComponentValue(field, "colors") && ! DXGetComponentValue(field, "front colors")) { Array array; float color[3]; int n; Object attr; if (!DXGetComponentValue(field, "data")) { array = (Array)DXGetComponentValue(field, "positions"); if (! array) return OK; DXGetArrayInfo(array, &n, NULL, NULL, NULL, NULL); attr = (Object)DXNewString("positions"); } else { attr = DXGetComponentAttribute(field, "data", "dep"); if (!attr || DXGetObjectClass(attr) != CLASS_STRING) { DXSetError(ERROR_INTERNAL, "#10241", "data"); return ERROR; } if (! strcmp(DXGetString((String)attr), "positions")) { array = (Array)DXGetComponentValue(field, "positions"); if (! array) return OK; DXGetArrayInfo(array, &n, NULL, NULL, NULL, NULL); } else if (! strcmp(DXGetString((String)attr), "connections")) { array = (Array)DXGetComponentValue(field, "connections"); if (! array) return OK; DXGetArrayInfo(array, &n, NULL, NULL, NULL, NULL); } else { DXSetError(ERROR_INTERNAL, "#11250", "data"); return ERROR; } } color[0] = 0.5; color[1] = 0.7; color[2] = 1.0; array = (Array)DXNewConstantArray(n, (Pointer)color, TYPE_FLOAT, CATEGORY_REAL, 1, 3); if (! array) return ERROR; if (! DXSetComponentValue(field, "colors", (Object)array)) { DXDelete((Object)array); return ERROR; } if (! DXSetComponentAttribute(field, "colors", "dep", attr)) return ERROR; } return OK; } static PseudoKey Hash(Key key) { Intercept *i; i = *(Intercept **)key; return i->i0*PRIME1 + i->i1*PRIME2; } static int Cmp(Key k0, Key k1) { Intercept *i0 = ((HashItem *)k0)->ptr; Intercept *i1 = ((HashItem *)k1)->ptr; return ! (i0->i0 == i1->i0 && i0->i1 == i1->i1); } static int Orientation(float *p, float *q, float *r, Vector *normal) { Vector qp, rp, sp, cross; float d; qp.x = q[0] - p[0]; qp.y = q[1] - p[1]; qp.z = q[2] - p[2]; rp.x = r[0] - p[0]; rp.y = r[1] - p[1]; rp.z = r[2] - p[2]; _dxfvector_cross_3D(&qp, &rp, &cross); d = _dxfvector_dot_3D(&cross, normal); if (d == 0.0) return INDETERMINATE; else if (d > 0) return ALIGNED; else return MIS_ALIGNED; } static Error OrientTriangles(Field field, Vector *n) { Array pA, cA; float *points; int *tri, i; int nPoints, nTriangles; if (DXEmptyField(field)) return OK; pA = (Array)DXGetComponentValue(field, "positions"); cA = (Array)DXGetComponentValue(field, "connections"); if (! DXGetArrayInfo(pA, &nPoints, NULL, NULL, NULL, NULL)) goto error; if (! DXGetArrayInfo(cA, &nTriangles, NULL, NULL, NULL, NULL)) goto error; points = (float *)DXGetArrayData(pA); tri = (int *)DXGetArrayData(cA); for (i = 0; i < nTriangles; i++, tri += 3) { float *p = points+3*tri[0]; float *q = points+3*tri[1]; float *r = points+3*tri[2]; if (Orientation(p, q, r, n) == MIS_ALIGNED) { int t = tri[0]; tri[0] = tri[1]; tri[1] = t; } } return OK; error: return ERROR; } static int DoConnectionsMap(Field field) { int i, j; char *name; Object attr; for (i = 0; DXGetEnumeratedComponentValue(field, i, &name); i++) { if (!strcmp(name, "connections")) continue; attr = DXGetComponentAttribute(field, name, "dep"); if (attr && !strcmp(DXGetString((String)attr), "connections")) return 1; } return 0; } static Error MakePlane(Vector *normal, Vector *point, Plane *plane) { if (vector_zero(normal)) { DXSetError(ERROR_BAD_PARAMETER, "#11822", "normal"); return ERROR; } if (! vector_normalize(normal, &plane->normal)) return ERROR; plane->offset = -vector_dot(&plane->normal, point); return OK; } static Error GetRegRegInfo(Field field, Plane plane, Plane *iPlane, Matrix *matrix) { int i; Array pA; Array cA; int meshOffsets[3]; Vector d[3]; Vector o; double A, B, C, D, l; /* * Get the regular positions information */ pA = (Array)DXGetComponentValue(field, "positions"); if (!DXQueryGridPositions(pA, NULL, NULL, (float *)&o, (float *)d)) return ERROR; /* * Get the planar equation in mesh coordinates relative to the * overall mesh origin. */ cA = (Array)DXGetComponentValue(field, "connections"); if (! DXQueryGridConnections(cA, NULL, NULL)) return ERROR; DXGetMeshOffsets((MeshArray)cA, meshOffsets); ADDMULVECTOR(d[0], -meshOffsets[0], o, o); ADDMULVECTOR(d[1], -meshOffsets[1], o, o); ADDMULVECTOR(d[2], -meshOffsets[2], o, o); /* * If I is a point in index space, M = (Mr, Mt) the transformation * from index space to positions space defined by a 3x3 matrix (the * deltas) and a 3-vector translation (the offset), and (Pabc, Pd) * a planar equation expressed in positions space, then the result * of applying the planar equation to the point I is: * P(IM) = (IM)Pabc + Pd * = (IMr + Mt)Pabc + Pd * = IMrPabc + (MtPabc + Pd) * = I(MrPabc) + (MtPabc + Pd) * In other words, we can define a new planar equation in index space * by the 3-vector MrPabc and the offset MtPabc + Pd. */ A = d[0].x*plane.normal.x + d[0].y*plane.normal.y + d[0].z*plane.normal.z; B = d[1].x*plane.normal.x + d[1].y*plane.normal.y + d[1].z*plane.normal.z; C = d[2].x*plane.normal.x + d[2].y*plane.normal.y + d[2].z*plane.normal.z; D = o.x*plane.normal.x + o.y*plane.normal.y + o.z*plane.normal.z + plane.offset; l = A*A + B*B + C*C; if (l == 0.0) { DXSetError(ERROR_INTERNAL, "#11831", "positions grid"); return ERROR; } iPlane->normal.x = A / l; iPlane->normal.y = B / l; iPlane->normal.z = C / l; iPlane->offset = D / l; matrix->A[0][0] = d[0].x; matrix->A[0][1] = d[0].y; matrix->A[0][2] = d[0].z; matrix->A[1][0] = d[1].x; matrix->A[1][1] = d[1].y; matrix->A[1][2] = d[1].z; matrix->A[2][0] = d[2].x; matrix->A[2][1] = d[2].y; matrix->A[2][2] = d[2].z; matrix->b[0] = o.x; matrix->b[1] = o.y; matrix->b[2] = o.z; return OK; } static int IsRegReg(Field field) { Array a; int n; a = (Array)DXGetComponentValue(field, "connections"); if (! a) return 0; if (! DXQueryGridConnections(a, &n, NULL)) return 0; if (n != 3) return 0; a = (Array)DXGetComponentValue(field, "positions"); if (! a) return 0; if (! DXQueryGridPositions(a, NULL, NULL, NULL, NULL)) return 0; return 1; } static int RegAddIntercept(SegList *intercepts, int axis, int i, int j) { RegIntercept *intercept; intercept = (RegIntercept *)DXNewSegListItem((SegList *)intercepts); if (! intercept) goto error; intercept->axis = axis; intercept->i = i; intercept->j = j; intercept->x.index = DXGetSegListItemCount((SegList *)intercepts) - 1; return intercept->x.index; error: return -1; } static Matrix PermuteMatrix(Matrix m, int *permute) { Matrix r; int i, j, k; i = permute[0]; j = permute[1]; k = permute[2]; r.A[0][0] = m.A[i][0]; r.A[0][1] = m.A[i][1]; r.A[0][2] = m.A[i][2]; r.A[1][0] = m.A[j][0]; r.A[1][1] = m.A[j][1]; r.A[1][2] = m.A[j][2]; r.A[2][0] = m.A[k][0]; r.A[2][1] = m.A[k][1]; r.A[2][2] = m.A[k][2]; r.b[0] = m.b[0]; r.b[1] = m.b[1]; r.b[2] = m.b[2]; return r; } #define ADDINTERCEPT(_axis, _i, _j, _knt) \ { \ RegIntercept *xcept = (RegIntercept *)DXNewSegListItem(iBuf); \ if (! xcept) \ goto error; \ xcept->axis = _axis; \ xcept->i = _i; \ xcept->j = _j; \ xcept->x.index = _knt; \ } #define ADDTRIANGLE(_i, _j, _k) \ { \ int *tri = (int *)DXNewSegListItem(tBuf); \ if (! tri) \ goto error; \ *tri++ = _i; \ *tri++ = _j; \ *tri++ = _k; \ } #define ADDMBUF(k) \ { \ int *m; \ m = (int *)DXNewSegListItem(mBuf); \ if (! m) \ goto error; \ *m = k; \ } #define ADDINTERCEPTROW(_axis, _index, _base) \ { \ int _k; \ for (_k = 0; _k < kknt; _k++) \ ADDINTERCEPT(_axis, _index, _k, (_base)+_k); \ } #define ADDTRIANGLEROW(_pbase, _cbase) \ { \ int _k, _pb = _pbase, _cb = _cbase; \ for (_k = 0; _k < (kknt-1); _k++, _pb++, _cb += slabCubes[2]) \ { \ if (DXIsElementInvalid(icHandle, _cb)) \ { \ cubesIgnored = 1; \ continue; \ } \ ADDTRIANGLE(_pb, _pb+kknt, _pb+kknt+1); \ ADDTRIANGLE(_pb, _pb+kknt+1, _pb+1); \ if (mBuf) { \ ADDMBUF(_cb); \ ADDMBUF(_cb); \ } \ } \ } static Field MTP_cubes_reg_AxisAligned(Field field, Plane plane, Matrix matrix) { int i, j, k, l, m, n, v; Array array; int slabPoints[3]; int slabCubes[3]; int count[3]; int meshOffsets[3]; int ix, jy; int ixMin, jyMin; int ixMax, jyMax; double aInv, bInv; SegList *iBuf = NULL; SegList *tBuf = NULL; SegList *mBuf = NULL; int permute[3]; double xfuzz, yfuzz, zfuzz; int cubesIgnored = 0; double dydx; int kknt; float x0, x1, y0, y1; double a, b, c, d; RegIntercept *xcept; SegListSegment *slist; int tot; InvalidComponentHandle icHandle = NULL; icHandle = DXCreateInvalidComponentHandle((Object)field, NULL, "connections"); if (! icHandle) goto error; array = (Array)DXGetComponentValue(field, "connections"); DXQueryGridConnections(array, NULL, count); DXGetMeshOffsets((MeshArray)array, meshOffsets); slabPoints[0] = count[1]*count[2]; slabPoints[1] = count[2]; slabPoints[2] = 1; slabCubes[0] = (count[1] - 1)*(count[2] - 1); slabCubes[1] = (count[2] - 1); slabCubes[2] = 1; /* * Choose the loop ordering so that the innermost loop tends to * be as long as possible (remember that we handle each cube that * intersects the plane exactly once - this ordering minimizes the * number of loops that have to be initialized). */ permute[0] = 0; permute[1] = 1; permute[2] = 2; a = plane.normal.x; if (a < 0) a = -a; b = plane.normal.y; if (b < 0) b = -b; c = plane.normal.z; if (c < 0) c = -c; if (b > a && b >= c) { SWAP(double, plane.normal.x, plane.normal.y); SWAP(int, count[0], count[1]); SWAP(int, meshOffsets[0], meshOffsets[1]); SWAP(int, slabPoints[0], slabPoints[1]); SWAP(int, slabCubes[0], slabCubes[1]); SWAP(int, permute[0], permute[1]); SWAP(double, a, b); } else if (c > a && c > b) { SWAP(double, plane.normal.x, plane.normal.z); SWAP(int, count[0], count[2]); SWAP(int, meshOffsets[0], meshOffsets[2]); SWAP(int, slabPoints[0], slabPoints[2]); SWAP(int, slabCubes[0], slabCubes[2]); SWAP(int, permute[0], permute[2]); SWAP(double, a, c); } if (c > b) { SWAP(double, plane.normal.y, plane.normal.z); SWAP(int, count[1], count[2]); SWAP(int, meshOffsets[1], meshOffsets[2]); SWAP(int, slabPoints[1], slabPoints[2]); SWAP(int, slabCubes[1], slabCubes[2]); SWAP(int, permute[1], permute[2]); } matrix = PermuteMatrix(matrix, permute); if (plane.normal.x != 0) aInv = -1.0 / plane.normal.x; else aInv = 0.0; if (plane.normal.y != 0) bInv = -1.0 / plane.normal.y; else bInv = 0.0; xfuzz = (plane.normal.x != 0.0) ? FUZZ / plane.normal.x : FUZZ; yfuzz = (plane.normal.y != 0.0) ? FUZZ / plane.normal.y : FUZZ; xfuzz = fabs(xfuzz); yfuzz = fabs(yfuzz); /* * Get the extents of this part of the nesh. */ ixMin = meshOffsets[0]; ixMax = ixMin + (count[0] - 1); jyMin = meshOffsets[1]; jyMax = jyMin + (count[1] - 1); kknt = count[2]; tBuf = DXNewSegList(3*sizeof(int), 2*(kknt-1), 2*(kknt-1)); if (! tBuf) goto error; iBuf = DXNewSegList(sizeof(RegIntercept), kknt, kknt); if (! tBuf) goto error; if (DoConnectionsMap(field)) { mBuf = DXNewSegList(sizeof(int), 2*(kknt-1), 2*(kknt-1)); if (! mBuf) goto error; } /* * Get endpoints of line in XY box */ if (plane.normal.x == 0.0) { x0 = ixMin; x1 = ixMax; y0 = plane.offset * bInv; SNAP(y0, yfuzz); y1 = y0; if (y0 < jyMin || y0 >= jyMax) goto return_empty; } else if (plane.normal.y == 0.0) { y0 = jyMin; y1 = jyMax; x0 = plane.offset * aInv; SNAP(x0, xfuzz); x1 = x0; if (x0 < ixMin || x0 >= ixMax) goto return_empty; } else { /* * clip line against 4 sides in XY plane */ x0 = ixMin; y0 = (plane.normal.x * x0 + plane.offset) * bInv; SNAP(y0, yfuzz); if (y0 < jyMin) { y0 = jyMin; x0 = (plane.normal.y * y0 + plane.offset) * aInv; SNAP(x0, xfuzz); if (x0 < ixMin || x0 > ixMax) goto return_empty; } else if (y0 > jyMax) { y0 = jyMax; x0 = (plane.normal.y * y0 + plane.offset) * aInv; SNAP(x0, xfuzz); if (x0 < ixMin || x0 > ixMax) goto return_empty; } x1 = ixMax; y1 = (plane.normal.x * x1 + plane.offset) * bInv; SNAP(y1, yfuzz); if (y1 < jyMin) { y1 = jyMin; x1 = (plane.normal.y * y1 + plane.offset) * aInv; SNAP(x1, xfuzz); if (x1 < ixMin || x1 > ixMax) goto return_empty; } else if (y1 > jyMax) { y1 = jyMax; x1 = (plane.normal.y * y1 + plane.offset) * aInv; SNAP(x1, xfuzz); if (x1 < ixMin || x1 > ixMax) goto return_empty; } } if ((x0 == x1 && ((plane.normal.x > 0 && x0 >= ixMax) || (plane.normal.x < 0 && x0 <= ixMin))) || (y0 == y1 && ((plane.normal.y > 0 && y0 >= jyMax) || (plane.normal.y < 0 && y0 <= jyMin))) || (x0 == x1 && y0 == y1)) goto return_empty; if (x0 == x1) { int pbase = 0; int cbase = ((int)x0 - meshOffsets[0])*slabCubes[0]; /* * if x0 == x1, its a slice parallel to the YZ plane. */ ADDINTERCEPTROW(XAXIS, jyMin, pbase); for (i = jyMin+1; i <= jyMax; i++, pbase += kknt, cbase += slabCubes[1]) { ADDINTERCEPTROW(XAXIS, i, pbase+kknt); ADDTRIANGLEROW(pbase, cbase); } } else { int pbase = 0; int cbase; int ix, iy; /* * swap endpoints so we increment up the x axis. */ if (x0 > x1) { float t = y0; y0 = y1; y1 = t; t = x0; x0 = x1; x1 = t; } dydx = (y1 - y0) / (x1 - x0); if (y0 == (int)y0 && dydx < 0) { ix = x0 - meshOffsets[0]; iy = (y0 - meshOffsets[1]) - 1; } else { ix = x0 - meshOffsets[0]; iy = y0 - meshOffsets[1]; } cbase = ix*slabCubes[0] + iy*slabCubes[1]; /* * Add the initial row of intercepts. The initial point (x0,y0) * is on the boundary of the XY box, so either x0 is ixMin or * y0 is jyMin or jyMax */ if (x0 == ixMin) ADDINTERCEPTROW(YAXIS, ixMin, pbase) else ADDINTERCEPTROW(XAXIS, (int)y0, pbase); /* * We will do an outer loop for each x-axis interval. * From the current point (x0, y0), we determine the end of * the current interval: either the next integral x value or, * if this exceeds x1, then x1. Then we determine the * y value at the interval end. */ while (x0 != x1 && y0 != y1) { float xn, yn; xn = floor(x0 + 1); if (xn >= x1) { xn = x1; yn = y1; } else { yn = y0 + (xn - x0)*dydx; SNAP(yn, yfuzz); } /* * Now for each integral y value between y0 and yn, set a row of * intercepts and triangles, bumping x0 and y0 appropriately. */ if (dydx > 0.0) { for (y0 = floor(y0+1); y0 < yn; y0 += 1, pbase += kknt) { ADDINTERCEPTROW(XAXIS, y0, pbase+kknt); ADDTRIANGLEROW(pbase, cbase); cbase += slabCubes[1]; } } else { for (y0 = ceil(y0-1); y0 > yn; y0 -= 1, pbase += kknt) { ADDINTERCEPTROW(XAXIS, y0, pbase+kknt); ADDTRIANGLEROW(pbase, cbase); cbase -= slabCubes[1]; } } /* * Add the intercepts at the X critical point. */ if (xn == (int)xn) ADDINTERCEPTROW(YAXIS, (int)xn, pbase+kknt) else ADDINTERCEPTROW(XAXIS, (int)yn, pbase+kknt); ADDTRIANGLEROW(pbase, cbase); pbase += kknt; cbase += slabCubes[0]; x0 = xn; y0 = yn; } } array = DXNewArray(TYPE_INT, CATEGORY_REAL, 1, 3); if (! array) goto error; if (! DXAllocateArray(array, DXGetSegListItemCount(tBuf))) { DXDelete((Object)array); goto error; } DXInitGetNextSegListSegment(tBuf); tot = 0; while (NULL != (slist = DXGetNextSegListSegment(tBuf))) { int n = DXGetSegListSegmentItemCount(slist); if (! DXAddArrayData(array, tot, n, DXGetSegListSegmentPointer(slist))) { DXDelete((Object)array); goto error; } tot += n; } DXSetComponentValue(field, "connections", (Object)array); array = NULL; if (! DXSetComponentAttribute(field, "connections", "element type", (Object)DXNewString("triangles"))) goto error; if (! DXInitGetNextSegListItem(iBuf)) goto error; a = plane.normal.x; b = plane.normal.y; c = plane.normal.z; d = plane.offset; while (NULL != (xcept = (RegIntercept *)DXGetNextSegListItem(iBuf))) { double t; int it; int I, J; switch(xcept->axis) { case XAXIS: I = (xcept->i - meshOffsets[1]); J = xcept->j; xcept->j += meshOffsets[2]; t = (b*xcept->i + d) * aInv; SNAP(t, xfuzz); xcept->intercept = t; it = floor(t); if ((it = floor(t)) == t) { if (it == ixMin) { xcept->x.d = 0.0; } else { xcept->x.d = 1.0; it -= 1; } } else xcept->x.d = t - it; xcept->x.i0 = (it - meshOffsets[0])*slabPoints[0] + I*slabPoints[1] + J*slabPoints[2]; SNAP(xcept->x.i0, xfuzz); xcept->x.i1 = xcept->x.i0 + slabPoints[0]; break; case YAXIS: I = (xcept->i - meshOffsets[0]); J = xcept->j; xcept->j += meshOffsets[2]; t = (a*xcept->i + d) * bInv; SNAP(t, yfuzz); xcept->intercept = t; if ((it = floor(t)) == t) { if (it == jyMin) { xcept->x.d = 0.0; } else { xcept->x.d = 1.0; it -= 1; } } else xcept->x.d = t - it; xcept->x.i0 = I*slabPoints[0] + (it - meshOffsets[1])*slabPoints[1] + J*slabPoints[2]; xcept->x.i1 = xcept->x.i0 + slabPoints[1]; break; } } if (! MapArrays(field, iBuf, mBuf, &matrix)) goto error; if (DXGetComponentValue(field, "invalid positions")) DXDeleteComponent(field, "invalid posititions"); if (DXGetComponentValue(field, "invalid connections")) DXDeleteComponent(field, "invalid connections"); if (cubesIgnored) { if (! DXInvalidateUnreferencedPositions((Object)field)) goto error; if (! DXCull((Object)field)) goto error; } if (! DXEndField(field)) goto error; if (icHandle) DXFreeInvalidComponentHandle(icHandle); if (iBuf) DXDeleteSegList(iBuf); if (mBuf) DXDeleteSegList(mBuf); if (tBuf) DXDeleteSegList(tBuf); return field; return_empty: if (icHandle) DXFreeInvalidComponentHandle(icHandle); if (iBuf) DXDeleteSegList(iBuf); if (mBuf) DXDeleteSegList(mBuf); if (tBuf) DXDeleteSegList(tBuf); { char *name; while (DXGetEnumeratedComponentValue(field, 0, &name)) DXDeleteComponent(field, name); if (!DXEndField(field)) goto error; } return field; error: if (icHandle) DXFreeInvalidComponentHandle(icHandle); DXDeleteSegList(iBuf); DXDeleteSegList(mBuf); DXDeleteSegList(tBuf); return NULL; }