/***********************************************************************/ /* 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" */ /***********************************************************************/ #include #include #include #include #include #include #include "vectors.h" #define TRUE 1 #define FALSE 0 #define MAX_V_PER_E 8 static Object _Sample(Object, int); Object Sample(Object o, int n) { Object copy = NULL; copy = DXCopy(o, COPY_STRUCTURE); if (! copy) goto error; if (! DXInvalidateConnections(copy)) goto error; if (! DXCreateTaskGroup()) return NULL; if (! _Sample(copy, n)) { DXAbortTaskGroup(); goto error; } if (! DXExecuteTaskGroup() || DXGetError() != ERROR_NONE) return NULL; return copy; error: if (copy != o) DXDelete((Object)copy); return NULL; } typedef float (*PFF)(); typedef Error (*PFE)(); typedef struct { Field field; int share; Array grid; } SampleArgs; static Object _Sample_CompositeField(CompositeField, int); static Object _Sample_MultiGrid(MultiGrid, int); static Object _Sample_Field(Field, int, Array); static Field _Sample_Field_Regular(Field, int, Array); static Field _Sample_Field_Irregular(Field, int, Array); static Error _Sample_Field_NoConnections(Field, int); static Error SpawnVolumeTasks(Object); static Error SF_Task(Pointer); static int Field_Count(Object, int *); static int FieldVolumeTask(Pointer); static float MemberVolume(Object); static Array _Sample_Grid(Object, int, int); static Array _Sample_Grid_Regular(Object, int, int); static Array _Sample_Grid_Irregular(Object, int, int); static Array _Solve_Grid(int, int, float *, float *, int *, int); static Error Cleanup_Sample_Field(Field); static void _Check_Sample_Count(Object o, int n); static float _dxfFieldVolume(Field); static float _dxfBoxVolume(Object, int *); static PFE _dxfGetSampleMethod(Field); static PFF _dxfGetVolumeMethod(Field); static Field CreateOutputReg(Field, Array); static int invert33(float *, float *); static Error invert(int, float *, float *); static Object _Sample_Field(Field field, int n, Array grid) { Object pos, con; int reg; if (DXEmptyField(field)) return (Object)DXEndField(DXNewField()); con = DXGetComponentValue(field, "connections"); if (! con) con = DXGetComponentValue(field, "faces"); if (! con) con = DXGetComponentValue(field, "polylines"); if (! con) { if (! _Sample_Field_NoConnections(field, n)) return NULL; } else { pos = DXGetComponentValue(field, "positions"); reg = DXQueryGridPositions((Array)pos, NULL, NULL, NULL, NULL) && DXQueryGridConnections((Array)con, NULL, NULL); if (reg && ! DXGetComponentValue(field, "invalid connections")) { if (! _Sample_Field_Regular(field, n, grid)) return NULL; } else { if (! _Sample_Field_Irregular(field, n, grid)) return NULL; } } if (! Cleanup_Sample_Field(field)) goto error; return (Object)field; error: return NULL; } static Field _Sample_Field_Regular(Field field, int n, Array grid) { Array array; int i, j, k, nd, ed; int oCounts[3], pCounts[3]; int pOffsets[3], oOffsets[3], gOffsets[3]; float gDeltas[9], pDeltas[9]; float gOrigin[3], pOrigin[3], oOrigin[3], tOrigin[3]; float bestG[3], bestP[3]; int bestI[3]; float p0[3], p1[3], g0[3], g1[3]; int deleteGrid; char *name; array = NULL; if (DXEmptyField(field)) return field; array = (Array)DXGetComponentValue(field, "positions"); if (! array) goto error; if (! DXQueryGridPositions(array, &nd, pCounts, p0, pDeltas)) goto irreg; ed = 0; for (i = 0; i < nd; i++) if (pCounts[i] != 1) ed ++; if (! grid) { deleteGrid = 1; grid = _Sample_Grid((Object)field, n, ed); } else deleteGrid = 0; if (! grid) goto error; if (! DXQueryGridPositions(grid, &nd, NULL, gOrigin, gDeltas)) goto irreg; /* * Move the non-trivial deltas to the lower slots */ k = 0; for (i = 0; i < nd; i++) if (pCounts[i] > 1) { for (j = 0; j < nd; j++) pDeltas[k*nd+j] = pDeltas[i*nd+j]; pCounts[k] = pCounts[i]; k++; } for (i = k; i < nd; i++) { for (j = 0; j < nd; j++) pDeltas[i*nd+j] = 0.0; pCounts[k] = 1; } /* * Get the corners of the original grid */ for (i = 0; i < nd; i++) { p1[i] = p0[i]; for (j = 0; j < nd; j++) p1[i] += (pCounts[j]-1)*pDeltas[j*nd+i]; } /* * Invert the delta matrix of the grid to get a mapping from (xyz) * space to grid (ijk) space. Subtract the grid origin from the * corners of the original grid and back-translate to get the corners * of the original grid in grid (ijk) coordinates. */ switch(nd) { case 3: { float matrix[9]; if (! invert33(gDeltas, matrix)) goto error; p0[0] -= gOrigin[0]; p0[1] -= gOrigin[1]; p0[2] -= gOrigin[2]; g0[0] = p0[0]*matrix[0]+p0[1]*matrix[3]+p0[2]*matrix[6]; g0[1] = p0[0]*matrix[1]+p0[1]*matrix[4]+p0[2]*matrix[7]; g0[2] = p0[0]*matrix[2]+p0[1]*matrix[5]+p0[2]*matrix[8]; p1[0] -= gOrigin[0]; p1[1] -= gOrigin[1]; p1[2] -= gOrigin[2]; g1[0] = p1[0]*matrix[0]+p1[1]*matrix[3]+p1[2]*matrix[6]; g1[1] = p1[0]*matrix[1]+p1[1]*matrix[4]+p1[2]*matrix[7]; g1[2] = p1[0]*matrix[2]+p1[1]*matrix[5]+p1[2]*matrix[8]; } break; case 2: { float matrix[4], m; m = gDeltas[0]*gDeltas[3] - gDeltas[1]*gDeltas[2]; if (m == 0) { DXSetError(ERROR_INVALID_DATA, "bad delta vectors"); goto error; } m = 1.0 / m; matrix[0] = gDeltas[3]*m; matrix[1] = -gDeltas[1]*m; matrix[2] = -gDeltas[2]*m; matrix[3] = gDeltas[0]*m; p0[0] -= gOrigin[0]; p0[1] -= gOrigin[1]; g0[0] = p0[0]*matrix[0] + p0[1]*matrix[2]; g0[1] = p0[0]*matrix[1] + p0[1]*matrix[3]; p1[0] -= gOrigin[0]; p1[1] -= gOrigin[1]; g1[0] = p1[0]*matrix[0] + p1[1]*matrix[2]; g1[1] = p1[0]*matrix[1] + p1[1]*matrix[3]; } break; case 1: { float m; m = gDeltas[0]; if (m == 0) { DXSetError(ERROR_INVALID_DATA, "bad delta vectors"); goto error; } m = 1.0 / m; g0[0] = (p0[0] - gOrigin[0])*m; g0[1] = (p0[0] - gOrigin[0])*m; g1[0] = (p1[0] - gOrigin[0])*m; g1[1] = (p1[0] - gOrigin[0])*m; } break; } for (i = ed; i < nd; i++) g0[i] = g1[i] = 0.0; for (j = 0; j < nd; j++) { if (ceil(g0[j]) > floor(g1[j])) goto empty; oOffsets[j] = ceil(g0[j]); oCounts[j] = (floor(g1[j]) - oOffsets[j]) + 1; if (oCounts[j] == 0) oCounts[j] = 1; } for (j = 0; j < nd; j++) { oOrigin[j] = gOrigin[j]; for (k = 0; k < nd; k++) oOrigin[j] += oOffsets[k]*gDeltas[k*nd + j]; } array = DXMakeGridPositionsV(nd, oCounts, oOrigin, gDeltas); if (! array) goto error; if (! CreateOutputReg(field, array)) goto error; if (deleteGrid) DXDelete((Object)grid); return field; empty: if (deleteGrid) DXDelete((Object)grid); while (DXGetEnumeratedComponentValue(field, 0, &name)) DXDeleteComponent(field, name); return field; irreg: return _Sample_Field_Irregular(field, n, grid); error: DXDelete((Object)array); return NULL; } static Field _Sample_Field_Irregular(Field field, int nSamples, Array grid) { PFE worker; Array pos; int r, s[16]; Category c; Type t; /* * Check first whether the input field is empty. */ if (DXEmptyField (field)) goto error; pos = (Array) DXGetComponentValue (field, "positions"); if (pos == NULL) { DXSetError (ERROR_MISSING_DATA, "missing positions"); goto error; } /* * Ascertain whether two dimensional or irregular data. * If two dimensional, find the bounding square and generate the * corresponding regular grid. * Otherwise, leave the work to specialized algorithm (worker routines). */ DXGetArrayInfo (pos, NULL, &t, &c, &r, s); if (t != TYPE_FLOAT) { DXSetError (ERROR_INVALID_DATA, "Positions must be of type float"); goto error; } if (c != CATEGORY_REAL) { DXSetError (ERROR_INVALID_DATA, "Positions must be of category real"); goto error; } worker = _dxfGetSampleMethod(field); if (! worker) goto error; if (((* worker) (field, nSamples, grid)) != OK) goto error; return field; error: return NULL; } static float rand1(); static Error _Sample_Field_NoConnections(Field field, int n) { Array old, new = NULL; int size, nDim, nPositions, i, j; InvalidComponentHandle ich = NULL; ArrayHandle handle = NULL; char *name, *oBuf = NULL; if (DXEmptyField(field)) return OK; old = (Array)DXGetComponentValue(field, "positions"); DXGetArrayInfo(old, &nPositions, NULL, NULL, NULL, NULL); ich = DXCreateInvalidComponentHandle((Object)field, NULL, "positions"); if (! ich) goto error; if (DXGetValidCount(ich) == 0) { char *names[256]; for (i = 0; !DXGetEnumeratedComponentValue(field, i, names+i); i++); for (i = i-1; i >= 0; i--) DXDeleteComponent(field, names[i]); goto done; } else { int knt; float f, fskip; float buf[32]; char *toDelete[256]; int nDelete; Type t; Category c; int rank, shape[32], size; n *= (float)nPositions / DXGetValidCount(ich); fskip = (float)nPositions / n; if (fskip <= 1.0) fskip = 1.0; for (knt = 0, f = 0.0; f < nPositions; f += fskip) if (DXIsElementValid(ich, (int)f)) knt ++; i = 0; nDelete = 0; while ((old = (Array)DXGetEnumeratedComponentValue(field, i++, &name)) != NULL) { Object attr; if (DXGetComponentAttribute(field, name, "der") || DXGetComponentAttribute(field, name, "ref") || ! (attr = DXGetComponentAttribute(field, name, "dep"))) { toDelete[nDelete++] = name; continue; } else { char *nPtr, *oPtr; if (DXGetObjectClass(attr) != CLASS_STRING) { DXSetError(ERROR_BAD_CLASS, "%s dependency attribute", name); goto error; } if (strcmp(DXGetString((String)attr), "positions")) { toDelete[nDelete++] = name; continue; } DXGetArrayInfo(old, NULL, &t, &c, &rank, shape); size = DXGetItemSize(old); new = DXNewArrayV(t, c, rank, shape); if (! new) goto error; if (! DXAddArrayData(new, 0, knt, NULL)) goto error; nPtr = DXGetArrayData(new); handle = DXCreateArrayHandle(old); if (! handle) goto error; oBuf = (char *)DXAllocate(size); for (f = 0.0, j = 0; f < nPositions; f += fskip, j = (int)f) if (DXIsElementValid(ich, j)) { oPtr = DXGetArrayEntry(handle, j, oBuf); memcpy(nPtr, oPtr, size); nPtr += size; } DXFreeArrayHandle(handle); handle = NULL; DXFree(oBuf); oBuf = NULL; if (! DXSetComponentValue(field, name, (Object)new)) goto error; } } for (i = 0; i < nDelete; i++) DXDeleteComponent(field, toDelete[i]); } done: DXFreeInvalidComponentHandle(ich); return OK; error: DXFreeArrayHandle(handle); DXFreeInvalidComponentHandle(ich); DXFree(oBuf); DXDelete((Object)new); return ERROR; } static Array _Sample_Grid(Object o, int nSamples, int eltDim) { Field field; Array pA; Class class; if ((class = DXGetObjectClass(o)) == CLASS_GROUP) class = DXGetGroupClass((Group)o); if (class == CLASS_FIELD) { field = (Field)o; } else if (class == CLASS_COMPOSITEFIELD) { int i; for (i = 0; ;i++) { field = (Field)DXGetEnumeratedMember((Group)o, i, NULL); if (! field) break; if (DXGetObjectClass((Object)field) != CLASS_FIELD) { DXSetError(ERROR_INVALID_DATA, "invalid composite field member"); return NULL; } if (! DXEmptyField(field)) break; } if (field == NULL) return NULL; } pA = (Array)DXGetComponentValue(field, "positions"); if (! pA) return NULL; if (DXQueryGridPositions(pA, NULL, NULL, NULL, NULL)) return _Sample_Grid_Regular(o, nSamples, eltDim); else return _Sample_Grid_Irregular(o, nSamples, eltDim); } static Array _Sample_Grid_Irregular(Object object, int nSamples, int eltDim) { int i, j, counts[3]; Point box[8], min, max; float origin[3], deltas[9]; if (! DXBoundingBox(object, box)) { DXSetError(ERROR_INVALID_DATA, "unable to compute bounding box"); return NULL; } min.x = min.y = min.z = DXD_MAX_FLOAT; max.x = max.y = max.z = -DXD_MAX_FLOAT; for (i = 0; i < 8; i++) { if (min.x > box[i].x) min.x = box[i].x; if (max.x < box[i].x) max.x = box[i].x; if (min.y > box[i].y) min.y = box[i].y; if (max.y < box[i].y) max.y = box[i].y; if (min.z > box[i].z) min.z = box[i].z; if (max.z < box[i].z) max.z = box[i].z; } origin[0] = min.x; origin[1] = min.y; origin[2] = min.z; deltas[0] = max.x-min.x; deltas[1] = 0; deltas[2] = 0; deltas[3] = 0; deltas[4] = max.y-min.y; deltas[5] = 0; deltas[6] = 0; deltas[7] = 0; deltas[8] = max.z-min.z; counts[0] = counts[1] = counts[2] = 2; if (deltas[0] == 0.0) { deltas[0] = 1.0; counts[0] = 1; }; if (deltas[4] == 0.0) { deltas[4] = 1.0; counts[1] = 1; }; if (deltas[8] == 0.0) { deltas[8] = 1.0; counts[2] = 1; }; return _Solve_Grid(3, nSamples, origin, deltas, counts, eltDim); } static Array _Solve_Grid(int nd, int nSamples, float *iOrigin, float *iDeltas, int *iCounts, int eltDim) { int tCounts[3], oCounts[3], i, j, k, n, c[3]; float A, B, a, l[3], d, oOrigin[3], oDeltas[9], tDeltas[9], divisor; /* * Compress the n-D grid so the first k are non-trivial */ k = 0; for (i = 0; i < nd; i++) if (iCounts[i] > 1) { for (j = 0; j < nd; j++) tDeltas[k*nd+j] = iDeltas[i*nd+j]; tCounts[k] = iCounts[i]; k++; } if (k == 0) { DXSetError(ERROR_INVALID_DATA, "degenerate grid"); goto error; } /* * Remaining dimensions are trivial. */ for (i = k; i < nd; i++) { for (j = 0; j < nd; j++) tDeltas[i*nd+j] = 0.0; tCounts[i] = 1; } /* * Get edge lengths of grid */ for (i = 0; i < k; i++) { l[i] = 0.0; for (j = 0; j < nd; j++) { d = tDeltas[i*nd + j]; l[i] += d * d; } l[i] = sqrt(l[i]) * (tCounts[i] - 1); } if (eltDim == -1) eltDim = k; switch(k) { case 1: oCounts[0] = nSamples; oCounts[1] = 1; oCounts[2] = 1; break; case 2: switch(eltDim) { case 1: oCounts[0] = nSamples * (l[0] / (l[1] + l[1])); oCounts[1] = nSamples * (l[1] / (l[1] + l[1])); oCounts[2] = 1; break; case 2: A = l[1] / l[0]; a = sqrt(nSamples / A); oCounts[0] = (int)(a + 0.5); if (oCounts[0] == 0) oCounts[0] = 1; oCounts[1] = (int)((nSamples / oCounts[0]) + 0.5); oCounts[2] = 1; break; default: DXSetError(ERROR_INTERNAL, "%d-D elements embedded in 2-D space", eltDim); goto error; } break; case 3: switch(eltDim) { case 1: oCounts[0] = nSamples * (l[0] / (l[1] + l[1] + l[2])); oCounts[1] = nSamples * (l[1] / (l[1] + l[1] + l[2])); oCounts[2] = nSamples * (l[2] / (l[1] + l[1] + l[2])); break; case 2: A = l[1] / l[0]; B = l[2] / l[0]; a = sqrt(3*nSamples / (A + B + A*B)); oCounts[0] = (int)(a + 0.5); oCounts[1] = (int)((a*A) + 0.5); oCounts[2] = (int)((a*B) + 0.5); break; case 3: A = l[0] / l[2]; B = l[1] / l[2]; a = pow(nSamples / (A*B), 0.333333333); oCounts[2] = (int)(a + 0.5); if (oCounts[2] < 1) oCounts[2] = 1; nSamples = ((float)nSamples / (float)oCounts[2]) + 0.5; A = l[1] / l[0]; a = sqrt(nSamples / A); oCounts[0] = (int)(a + 0.5); if (oCounts[0] == 0) oCounts[0] = 1; oCounts[1] = (int)((nSamples / oCounts[0]) + 0.5); break; default: DXSetError(ERROR_INTERNAL, "%d-D elements embedded in 2-D space", eltDim); goto error; } break; } for (i = 0; i < k; i++) if (oCounts[i] < 1) oCounts[i] = 1; for (i = 0; i < nd; i++) { if (oCounts[i] > 0) { d = ((float)tCounts[i] - 1) / oCounts[i]; for (j = 0; j < nd; j++) oDeltas[i*nd + j] = d * tDeltas[i*nd + j]; k++; } else for (j = 0; j < nd; j++) oDeltas[i*nd + j] = 0.0; } for (i = 0; i < nd; i++) { oOrigin[i] = iOrigin[i]; for (j = 0; j < nd; j++) oOrigin[i] += oDeltas[j*nd+i] / 2.0; } return DXMakeGridPositionsV(nd, oCounts, oOrigin, oDeltas); error: return NULL; } #define INTERPOLATE_3D(type, round) \ { \ Vector *p = (Vector *)pList; \ char *sData = (char *)DXGetArrayData(sA); \ type *dPtr = (type *)DXGetArrayData(dA); \ int s; \ \ for (s = 0; s < nSamples; s++, p++) \ { \ float x, y, z, w[8], A, B, C, D; \ int ix, iy, iz, base, j, k; \ type *sPtr; \ \ if (p->x == oCounts[0]-1) {ix = oCounts[0]-2; x = 1.0;} \ else {ix = (int)p->x; x = p->x - ix;} \ if (p->y == oCounts[1]-1) {iy = oCounts[1]-2; y = 1.0;} \ else {iy = (int)p->y; y = p->y - iy;} \ if (p->z == oCounts[2]-1) {iz = oCounts[2]-2; z = 1.0;} \ else {iz = (int)p->z; z = p->z - iz;} \ \ A = x * y; B = x * z; C = y * z; D = A * z; \ \ w[7] = D; \ w[6] = C - D; \ w[5] = B - D; \ w[4] = z - B - C + D; \ w[3] = A - D; \ w[2] = y - C - A + D; \ w[1] = x - B - A + D; \ w[0] = 1.0 - z - y - x + A + B + C - D; \ \ base = ix*skip[0] + iy*skip[1] + iz*skip[2]; \ \ sPtr = (type *)(sData + base*itemSize); \ \ for (j = 0; j < itemValues; j++) \ fBuf[j] = w[0] * sPtr[j]; \ \ for (j = 1; j < 8; j++) \ { \ int v = base + (j & 0x1 ? skip[0] : 0) \ + (j & 0x2 ? skip[1] : 0) \ + (j & 0x4 ? skip[2] : 0); \ type *sPtr = (type *)(sData + v*itemSize); \ for (k = 0; k < itemValues; k++) \ fBuf[k] += w[j] * sPtr[k]; \ } \ \ for (j = 0; j < itemValues; j++) \ dPtr[j] = fBuf[j] + round; \ \ dPtr += itemValues; \ } \ } #define INTERPOLATE_2D(type, round) \ { \ float *p = (float *)pList; \ char *sData = (char *)DXGetArrayData(sA); \ type *dPtr = (type *)DXGetArrayData(dA); \ int s; \ \ for (s = 0; s < nSamples; s++, p += 2) \ { \ float x, y, w[4]; \ int ix, iy, base, j, k; \ type *sPtr; \ \ if (p[0] == oCounts[0]-1) {ix = oCounts[0]-2; x = 1.0;} \ else {ix = (int)p[0]; x = p[0] - ix;} \ if (p[1] == oCounts[1]-1) {iy = oCounts[1]-2; y = 1.0;} \ else {iy = (int)p[1]; y = p[1] - iy;} \ \ w[0] = (1 - x)*(1 - y); \ w[1] = x*(1 - y); \ w[2] = (1 - x)*y; \ w[3] = x*y; \ \ base = ix*skip[0] + iy*skip[1]; \ \ sPtr = (type *)(sData + base*itemSize); \ \ for (j = 0; j < itemValues; j++) \ fBuf[j] = w[0] * sPtr[j]; \ \ for (j = 1; j < 4; j++) \ { \ int v = base + (j & 0x1 ? skip[0] : 0) \ + (j & 0x2 ? skip[1] : 0); \ type *sPtr = (type *)(sData + v*itemSize); \ for (k = 0; k < itemValues; k++) \ fBuf[k] += w[j] * sPtr[k]; \ } \ \ for (j = 0; j < itemValues; j++) \ dPtr[j] = fBuf[j] + round; \ \ dPtr += itemValues; \ } \ } #define INTERPOLATE_1D(type, round) \ { \ float *p = (float *)pList; \ char *sData = (char *)DXGetArrayData(sA); \ type *dPtr = (type *)DXGetArrayData(dA); \ int s; \ \ for (s = 0; s < nSamples; s++, p ++) \ { \ float x; \ int ix, j, k; \ type *src0, *src1; \ \ if (*p == oCounts[0]-1) {ix = oCounts[0]-2; x = 1.0;} \ else {ix = (int)*p; x = p[0] - ix;} \ \ src0 = (type *)sData + (ix * itemSize); \ src1 = src0 + itemSize; \ \ for (j = 0; j < itemValues; j++) \ dPtr[j] = (x*src1[j] + (1 - x)*src0[j]) + round; \ \ dPtr += itemValues; \ } \ } static Error InvertRegular(int, float *, float *, float *, float *); static void MatMult(int, float *, float *, float *, float *, float *, float *); static Field CreateOutputReg(Field f, Array np) { Array op; int size[3], skip[3], oCounts[3], nCounts[3], nd, ndInterp, i, j, k; float oOrigin[3], oDeltas[9]; float nOrigin[3], nDeltas[9]; float mOrigin[3], mDeltas[9]; float oOriginInv[3], oDeltasInv[9]; int *cList, *cc; float *pList, *pp; int indices[3]; int nSamples; Array sA, dA = NULL; char *dptr; char *name, *toDelete[100], nDelete = 0; Object at; float *fBuf = NULL; int *p; InvalidComponentHandle inInvalids = NULL, outInvalids = NULL; indices[0] = 0; indices[1] = 0; indices[2] = 0; op = (Array)DXGetComponentValue(f, "positions"); if (! op) goto error; if (! DXQueryGridPositions(op, &nd, oCounts, oOrigin, oDeltas)) goto error; for (i = j = 0; i < nd; i++) { if (oCounts[i] != 1) { if (i != j) { oCounts[j] = oCounts[i]; for (k = 0; k < nd; k++) oDeltas[(j*nd)+k] = oDeltas[(i*nd)+k]; } j++; } } for (i = j; i < nd; i++) { oCounts[i] = 1; for (k = 0; k < nd; k++) oDeltas[(i*nd)+k] = 0.0; } ndInterp = j; if (! DXQueryGridPositions(np, NULL, nCounts, nOrigin, nDeltas)) goto error; if (! InvertRegular(nd, oOrigin, oDeltas, oOriginInv, oDeltasInv)) goto error; MatMult(nd, nOrigin, nDeltas, oOriginInv, oDeltasInv, mOrigin, mDeltas); DXGetArrayInfo(np, &nSamples, NULL, NULL, NULL, NULL); cList = (int *)DXAllocate(nSamples*sizeof(int)); pList = (float *)DXAllocate(nSamples*ndInterp*sizeof(float)); if (! pList || ! cList) goto error; size[nd-1] = 1; skip[nd-1] = 1; for (i = nd - 2; i >= 0; i --) { size[i] = (oCounts[i+1] == 1 ? 1 : (oCounts[i+1] - 1)) * size[i+1]; skip[i] = (i == nd - 2) ? oCounts[i+1] : skip[i+1] * oCounts[i+1]; } /* * For each output sample, get the index of the connection it * lies in (the c array) and the point in the original position * index space that the sample corresponds to. */ pp = pList; cc = cList; for (i = 0; i < nSamples; i++) { *cc = 0; for (j = 0; j < ndInterp; j++) { pp[j] = mOrigin[j]; for (k = 0; k < nd; k++) pp[j] += indices[k]*mDeltas[(k*nd) + j]; *cc += ((int)pp[j])*size[j]; } cc++; pp += ndInterp; if (i != nSamples-1) { indices[nd-1]++; for (j = nd-1; j >= 0; j--) if (indices[j] >= nCounts[j]) { indices[j] = 0; indices[j-1]++; } else break; } } if (DXGetComponentValue(f, "invalid positions") || DXGetComponentValue(f, "invalid connections")) { inInvalids = DXCreateInvalidComponentHandle((Object)f, NULL, "connections"); outInvalids = DXCreateInvalidComponentHandle((Object)np, NULL, "positions"); if (! inInvalids || ! outInvalids) goto error; for (i = 0; i < nSamples; i++) if (DXIsElementInvalid(inInvalids, cList[i])) DXSetElementInvalid(outInvalids, i); if (! DXSaveInvalidComponent(f, outInvalids)) goto error; DXFreeInvalidComponentHandle(inInvalids); inInvalids = NULL; DXFreeInvalidComponentHandle(outInvalids); outInvalids = NULL; } DXDeleteComponent(f, "connections"); i = 0; while(NULL != (sA = (Array)DXGetEnumeratedComponentValue(f, i++, &name))) { Object at; Type t; Category c; int r, shape[32]; int itemSize, itemValues; if (!strcmp(name, "positions") || !strcmp(name, "invalid positions")) continue; if (DXGetComponentAttribute(f, name, "ref") || DXGetComponentAttribute(f, name, "der") || NULL == (at = DXGetComponentAttribute(f, name, "dep")) || !strcmp(name, "invalid connections")) { toDelete[nDelete++] = name; continue; } if (DXGetObjectClass(at) != CLASS_STRING) { DXSetError(ERROR_BAD_CLASS, "%s dependency attribute", name); goto error; } DXGetArrayInfo(sA, NULL, &t, &c, &r, shape); itemSize = DXGetItemSize(sA); itemValues = itemSize / DXTypeSize(t); if (DXQueryConstantArray(sA, NULL, NULL)) { dA = (Array)DXNewConstantArrayV(nSamples, DXGetConstantArrayData(sA), t, c, r, shape); if (! dA) goto error; } else { dA = DXNewArrayV(t, c, r, shape); if (! dA) goto error; if (! DXAddArrayData(dA, 0, nSamples, NULL)) goto error; if (! strcmp(DXGetString((String)at), "positions")) { int nItems = DXGetItemSize(dA) / DXTypeSize(t); fBuf = (float *)DXAllocate(nItems*sizeof(float)); if (! fBuf) goto error; switch(ndInterp) { case 3: switch(t) { case TYPE_FLOAT: INTERPOLATE_3D(float, 0.0); break; case TYPE_DOUBLE: INTERPOLATE_3D(double, 0.0); break; case TYPE_SHORT: INTERPOLATE_3D(short, 0.5); break; case TYPE_USHORT: INTERPOLATE_3D(ushort, 0.5); break; case TYPE_INT: INTERPOLATE_3D(int, 0.5); break; case TYPE_UINT: INTERPOLATE_3D(uint, 0.5); break; case TYPE_UBYTE: INTERPOLATE_3D(ubyte, 0.5); break; case TYPE_BYTE: INTERPOLATE_3D(byte, 0.5); break; } break; case 2: switch(t) { case TYPE_FLOAT: INTERPOLATE_2D(float, 0.0); break; case TYPE_DOUBLE: INTERPOLATE_2D(double, 0.0); break; case TYPE_SHORT: INTERPOLATE_2D(short, 0.5); break; case TYPE_USHORT: INTERPOLATE_2D(ushort, 0.5); break; case TYPE_INT: INTERPOLATE_2D(int, 0.5); break; case TYPE_UINT: INTERPOLATE_2D(uint, 0.5); break; case TYPE_UBYTE: INTERPOLATE_2D(ubyte, 0.5); break; case TYPE_BYTE: INTERPOLATE_2D(byte, 0.5); break; } break; case 1: switch(t) { case TYPE_FLOAT: INTERPOLATE_1D(float, 0.0); break; case TYPE_DOUBLE: INTERPOLATE_1D(double, 0.0); break; case TYPE_SHORT: INTERPOLATE_1D(short, 0.5); break; case TYPE_USHORT: INTERPOLATE_1D(ushort, 0.5); break; case TYPE_INT: INTERPOLATE_1D(int, 0.5); break; case TYPE_UINT: INTERPOLATE_1D(uint, 0.5); break; case TYPE_UBYTE: INTERPOLATE_1D(ubyte, 0.5); break; case TYPE_BYTE: INTERPOLATE_1D(byte, 0.5); break; } break; } DXFree((Pointer)fBuf); fBuf = 0; } else if (! strcmp(DXGetString((String)at), "connections")) { char *src, *dst; int s; src = (char *)DXGetArrayData(sA); dst = (char *)DXGetArrayData(dA); for (s = 0; s < nSamples; s++) { memcpy(dst, src+(cList[s]*itemSize), itemSize); dst += itemSize; } } else { DXSetError(ERROR_BAD_CLASS, "%s dependency attribute", name); goto error; } } if (! DXSetComponentValue(f, name, (Object)dA)) goto error; dA = NULL; if (! DXSetComponentAttribute(f, name, "dep", (Object)DXNewString("positions"))) goto error; } for (i = 0; i < nDelete; i++) DXDeleteComponent(f, toDelete[i]); if (! DXSetComponentValue(f, "positions", (Object)np)) goto error; if (outInvalids) DXFreeInvalidComponentHandle(outInvalids); if (inInvalids) DXFreeInvalidComponentHandle(inInvalids); DXFree((Pointer)pList); DXFree((Pointer)cList); return f; error: if (outInvalids) DXFreeInvalidComponentHandle(outInvalids); if (inInvalids) DXFreeInvalidComponentHandle(inInvalids); DXFree((Pointer)dA); DXFree((Pointer)fBuf); DXFree((Pointer)pList); DXFree((Pointer)cList); return NULL; } static Error InvertRegular(int nd, float *sO, float *sD, float *dO, float *dD) { float m; switch(nd) { case 1: if (sD[0] == 0.0) goto error; dD[0] = 1.0 / sD[0]; dO[0] = -sO[0]*dD[0]; break; case 2: m = sD[0]*sD[3] - sD[1]*sD[2]; if (m == 0) goto error; m = 1.0 / m; dD[0] = sD[3]*m; dD[1] = -sD[1]*m; dD[2] = -sD[2]*m; dD[3] = sD[0]*m; dO[0] = -(sO[0]*dD[0] + sO[1]*dD[2]); dO[1] = -(sO[0]*dD[1] + sO[1]*dD[3]); break; case 3: if (! invert33(sD, dD)) goto error; dO[0] = -(sO[0]*dD[0] + sO[1]*dD[3] + sO[2]*dD[6]); dO[1] = -(sO[0]*dD[1] + sO[1]*dD[4] + sO[2]*dD[7]); dO[2] = -(sO[0]*dD[2] + sO[1]*dD[5] + sO[2]*dD[8]); break; default: DXSetError(ERROR_NOT_IMPLEMENTED, ">3 dimensions"); return ERROR; } return OK; error: DXSetError(ERROR_INVALID_DATA, "non-invertible data"); return ERROR; } static void MatMult(int nd, float *o0, float *d0, float *o1, float *d1, float *or, float *dr) { int i, j, k; for (i = 0; i < nd; i++) { for (j = 0; j < nd; j++) { dr[i*nd+j] = 0.0; for (k = 0; k < nd; k++) dr[i*nd+j] += d0[i*nd+k]*d1[k*nd+j]; } or[i] = o1[i]; for (j = 0; j < nd; j++) or[i] += o0[j]*d1[j*nd+i]; } } static Object _Sample(Object o, int n) { Class class; Object c, new; int i; SampleArgs args; class = DXGetObjectClass(o); if (class == CLASS_GROUP) class = DXGetGroupClass((Group)o); switch(class) { case CLASS_SERIES: { float position; i = 0; while (NULL != (c = DXGetSeriesMember((Series)o, i, &position))) { new = _Sample(c, n); if (! new) goto error; if (new != c) DXSetSeriesMember((Series)o, i, position, (Object)new); i++; } break; } case CLASS_GROUP: { char *name; i = 0; while (NULL != (c = DXGetEnumeratedMember((Group)o, i, &name))) { new = _Sample(c, n); if (! new) goto error; if (new != c) if (name) DXSetMember((Group)o, name, (Object)new); else DXSetEnumeratedMember((Group)o, i, (Object)new); i++; } break; } case CLASS_COMPOSITEFIELD: _Check_Sample_Count(o, n); if (NULL == (o = _Sample_CompositeField((CompositeField)o, n))) goto error; break; case CLASS_MULTIGRID: _Check_Sample_Count(o, n); if (NULL == (o = _Sample_MultiGrid((MultiGrid)o, n))) goto error; break; case CLASS_FIELD: _Check_Sample_Count(o, n); args.field = (Field)o; args.share = n; args.grid = NULL; if (! DXAddTask(SF_Task, (Pointer)&args, sizeof(args), 1.0)) goto error; break; case CLASS_XFORM: { Object new; DXGetXformInfo((Xform)o, &c, NULL); if (NULL == (new = _Sample(c, n))) goto error; if (new != c) DXSetXformObject((Xform)o, new); break; } case CLASS_CLIPPED: { Object new; DXGetClippedInfo((Clipped)o, &c, NULL); if (NULL == (new = _Sample(c, n))) goto error; if (new != c) DXSetClippedObjects((Clipped)o, new, NULL); break; } default: break; } return o; error: return NULL; } static int Field_Count(Object o, int *n) { int i, tot, t; Object c; switch(DXGetObjectClass(o)) { case CLASS_GROUP: i = 0; tot = 0; while (NULL != (c = DXGetEnumeratedMember((Group)o, i++, NULL))) { t = Field_Count(c, n); if (t < 0) return -1; tot += t; } return tot; case CLASS_FIELD: if (! DXEmptyField((Field)o)) { char *str; Object attr; int m; if (DXGetComponentValue((Field)o, "connections")) { attr = DXGetComponentAttribute((Field)o, "connections", "element type"); if (! attr) { DXSetError(ERROR_MISSING_DATA, "#10255", "connections", "element type"); return -1; } str = DXGetString((String)attr); if (!strcmp(str, "triangles")) m = 2; else if (!strcmp(str, "quads")) m = 2; else if (!strcmp(str, "lines")) m = 1; else if (!strcmp(str, "cubes")) m = 3; else if (!strcmp(str, "tetrahedra")) m = 3; else { DXSetError(ERROR_INVALID_DATA, "#10410", str); return -1; } if (*n == -1) *n = m; else if (*n != m) { DXSetError(ERROR_INVALID_DATA, "#10420"); return -1; } } else { if (*n == -1) *n = 0; else if (*n > 0) { DXSetError(ERROR_INVALID_DATA, "#10420"); return -1; } } } return 1; case CLASS_XFORM: DXGetXformInfo((Xform)o, &c, NULL); return Field_Count(c, n); case CLASS_CLIPPED: DXGetClippedInfo((Clipped)o, &c, NULL); return Field_Count(c, n); default: return 0; } } static Error SpawnVolumeTasks(Object o) { int i, tot; Object c; Pointer ptr; Object attr; switch(DXGetObjectClass(o)) { case CLASS_GROUP: i = 0; while (NULL != (c = DXGetEnumeratedMember((Group)o, i++, NULL))) if (! SpawnVolumeTasks(c)) return ERROR; return OK; case CLASS_FIELD: if (NULL == (attr = DXGetAttribute(o, "volume"))) if (! DXAddTask(FieldVolumeTask, (Pointer)o, 0, 1.0)) return ERROR; return OK; case CLASS_XFORM: DXGetXformInfo((Xform)o, &c, NULL); return SpawnVolumeTasks(c); case CLASS_CLIPPED: DXGetClippedInfo((Clipped)o, &c, NULL); return SpawnVolumeTasks(c); default: return OK; } } static float MemberVolume(Object o) { int i; float tot; Object c; Pointer ptr; Object attr; switch(DXGetObjectClass(o)) { case CLASS_GROUP: i = 0; tot = 0.0; while (NULL != (c = DXGetEnumeratedMember((Group)o, i++, NULL))) tot += MemberVolume(c); return tot; case CLASS_FIELD: if (NULL == (attr = DXGetAttribute(o, "volume"))) { DXSetError(ERROR_INTERNAL, "missing field volume"); return -1; } return *(float *)DXGetConstantArrayData((Array)attr); case CLASS_XFORM: DXGetXformInfo((Xform)o, &c, NULL); return MemberVolume(c); case CLASS_CLIPPED: DXGetClippedInfo((Clipped)o, &c, NULL); return MemberVolume(c); default: return 0.0; } } static Error SF_Task(Pointer ptr) { SampleArgs *args = (SampleArgs *)ptr; if (! _Sample_Field(args->field, args->share, args->grid)) { DXDelete((Object)args->grid); return ERROR; } DXDelete((Object)args->grid); return OK; } static Object _Sample_CompositeField(CompositeField cf, int n) { float tv, bv, mv, v; int nd, i, nElts; SampleArgs args; Array grid; int eltDim = -1; Object c; nElts = Field_Count((Object)cf, &eltDim); if (nElts < 0) return NULL; if (! DXCreateTaskGroup()) goto error; nElts = 0; if (! SpawnVolumeTasks((Object)cf)) { DXAbortTaskGroup(); goto error; } if (! DXExecuteTaskGroup() || DXGetError() != ERROR_NONE) goto error; tv = 0; for (i = 0; NULL != (c = DXGetEnumeratedMember((Group)cf, i, NULL)); i++) { if (DXGetObjectClass(c) != CLASS_FIELD) { DXSetError(ERROR_INVALID_DATA, "composite members must be fields"); goto error; } v = MemberVolume(c); if (v < 0) goto error; tv += v; } if (tv == 0.0) goto empty; if (eltDim != 0) { if (-1 == (bv = _dxfBoxVolume((Object)cf, &nd))) return NULL; if (nd == eltDim) n *= bv / tv; grid = _Sample_Grid((Object)cf, n, eltDim); if (! grid) goto error; } for (i = 0; NULL != (c = DXGetEnumeratedMember((Group)cf, i, NULL)); i++) { args.grid = (Array)DXReference((Object)grid); mv = MemberVolume(c); if (mv > 0.0) { args.share = 1 + n*(mv/tv); args.field = (Field)c; if (! DXAddTask(SF_Task, (Pointer)&args, sizeof(args), 1.0)) goto error; } } return (Object)cf; empty: return (Object)DXEndField(DXNewField()); error: return NULL; } static Object _Sample_MultiGrid(MultiGrid mg, int n) { float tv, bv, v; int nd, i, nElts; SampleArgs args; int eltDim = -1; Object c; if (-1 == (bv = _dxfBoxVolume((Object)mg, &nd))) return NULL; nElts = Field_Count((Object)mg, &eltDim); if (nElts < 0) return NULL; if (! DXCreateTaskGroup()) goto error; nElts = 0; if (! SpawnVolumeTasks((Object)mg)) { DXAbortTaskGroup(); goto error; } if (! DXExecuteTaskGroup() || DXGetError() != ERROR_NONE) goto error; tv = 0; for (i = 0; NULL != (c = DXGetEnumeratedMember((Group)mg, i, NULL)); i++) { v = MemberVolume(c); if (v < 0) goto error; tv += v; } if (tv <= 0.0) { DXSetError(ERROR_INTERNAL, "multigrid volume is <= 0.0"); goto error; } for (i = 0; NULL != (c = DXGetEnumeratedMember((Group)mg, i, NULL)); i++) { if (! _Sample(c, 1+n*MemberVolume(c)/tv)) goto error; } return (Object)mg; error: return NULL; } static Error PolylinesSample(Field, int, Array); static Error FacesSample(Field, int, Array); static float LineLength(float **, int); static Error LineSample(Field, int, Array); static Error LineSample_Irreg(Field, int, Array); static Error LineSample_PathArray(Field, int, Array); static float QuadArea(float **, int); static Error QuadSample(Field, int, Array); static int coords_quad(float *, float **, float *); static float TriangleArea(float **, int); static Error TriangleSample(Field, int, Array); static int coords_triangle(float *, float **, float *); static float CubeVolume(float **, int); static Error CubeSample(Field, int, Array); static int coords_cube(float *, float **, float *); static float TetrahedraVolume(float **, int); static Error TetrahedraSample(Field, int, Array); static int coords_tetrahedra (float *, float **, float *); static int _invert33(float *, float *, float *, float *); static Field CreateOutputIrreg(Field, SegList *, SegList *); static Field CreateOutputPolylines(Field, SegList *, SegList *); static Field CreateOutputFaces(Field, SegList *, SegList *); static Error PlanarSurface(Field, int, int (*inside)(), int, Array); static Error CurvedSurface(Field, int, int (*inside)(), int, Array); static Error Volume(Field, int, int (*inside)(), int, Array); static PFE _dxfGetSampleMethod(Field field) { Object obj; char *primitive; obj = DXGetComponentAttribute (field, "connections", "element type"); if (obj == NULL) { if (DXGetComponentValue(field, "polylines")) primitive = "polylines"; else if (DXGetComponentValue(field, "faces")) primitive = "faces"; else { DXSetError(ERROR_MISSING_DATA, "no element type attribute"); goto error; } } else { primitive = DXGetString ((String) obj); if (primitive == NULL) { DXSetError(ERROR_MISSING_DATA,"unable to access element type attribute"); goto error; } } if (! strcmp (primitive, "tetrahedra")) { return TetrahedraSample; } else if (! strcmp (primitive, "cubes")) { return CubeSample; } else if (! strcmp (primitive, "triangles")) { return TriangleSample; } else if (! strcmp (primitive, "quads")) { return QuadSample; } else if (! strcmp (primitive, "lines")) { return LineSample; } else if (! strcmp (primitive, "polylines")) { return PolylinesSample; } else if (! strcmp (primitive, "faces")) { return FacesSample; } else { DXSetError (ERROR_INVALID_DATA, "unknown primitive: %s", primitive); goto error; } error: return NULL; } static Error TriangleSample(Field field, int nSamples, Array grid) { int i, n, nd; int counts[3]; Array pA; pA = (Array)DXGetComponentValue(field, "positions"); if (! DXQueryGridPositions(pA, &n, counts, NULL, NULL)) { DXGetArrayInfo(pA, NULL, NULL, NULL, NULL, &nd); } else { for (i = nd = 0; i < n; i++) if (counts[i] > 1) nd ++; } if (nd == 2) return PlanarSurface(field, 3, coords_triangle, nSamples, grid); else if (nd == 3) return CurvedSurface(field, 3, coords_triangle, nSamples, grid); else { DXSetError(ERROR_INVALID_DATA, "triangles require 2 or 3D positions"); goto error; } error: return ERROR; } static Error QuadSample(Field field, int nSamples, Array grid) { int i, n, nd; int counts[3]; Array pA; pA = (Array)DXGetComponentValue(field, "positions"); if (! DXQueryGridPositions(pA, &n, counts, NULL, NULL)) { DXGetArrayInfo(pA, NULL, NULL, NULL, NULL, &nd); } else { for (i = nd = 0; i < n; i++) if (counts[i] > 1) nd ++; } if (nd == 2) return PlanarSurface(field, 4, coords_quad, nSamples, grid); else if (nd == 3) return CurvedSurface(field, 4, coords_quad, nSamples, grid); else { DXSetError(ERROR_INVALID_DATA, "triangles require 2 or 3D positions"); goto error; } error: return ERROR; } static Error TetrahedraSample(Field field, int nSamples, Array grid) { return Volume(field, 4, coords_tetrahedra, nSamples, grid); } static Error CubeSample(Field field, int nSamples, Array grid) { return Volume(field, 8, coords_cube, nSamples, grid); } static PFF _dxfGetVolumeMethod(Field field) { Object obj; char *primitive; obj = DXGetComponentAttribute (field, "connections", "element type"); if (obj == NULL) { DXSetError(ERROR_MISSING_DATA, "no element type attribute"); goto error; } primitive = DXGetString ((String) obj); if (primitive == NULL) { DXSetError(ERROR_MISSING_DATA,"unable to access element type attribute"); goto error; } if (! strcmp (primitive, "tetrahedra")) { return TetrahedraVolume; } else if (! strcmp (primitive, "cubes")) { return CubeVolume; } else if (! strcmp (primitive, "triangles")) { return TriangleArea; } else if (! strcmp (primitive, "quads")) { return QuadArea; } else if (! strcmp (primitive, "lines")) { return LineLength; } else { DXSetError (ERROR_INVALID_DATA, "invalid primitive: %s", primitive); goto error; } error: return NULL; } typedef struct { int index; float weights[1]; } Interp; static Error PlanarSurface(Field field, int vPerE, int (*inside)(), int nSamples, Array grid) { Array a_connections, a_positions; int n_connections, n_positions; int nOut; float *eltPoints[MAX_V_PER_E]; int basex, pDim, gDim; float gridx[2], gridy[2]; char *grid_status = NULL; int grid_pt; float mint[2], maxt[2]; float mintrel[2], maxtrel[2]; int mintgrid[2], maxtgrid[2]; int i, j, k, l, n; float gDeltas[9], gOrigin[3]; int gCounts[3]; int deleteGrid; SegList *pList = NULL, *iList = NULL; int size; float pBuf[MAX_V_PER_E][3]; float tBuf[MAX_V_PER_E][2]; ArrayHandle pHandle = NULL; ArrayHandle cHandle = NULL; int cBuf[32], *elt; InvalidComponentHandle iHandle = NULL; float mA[9], mB[3]; float mAinv[9], mBinv[3]; char *name; if (! grid) { float bVol, tVol, rat; float deltas[9]; int nd, counts[3]; Array pA = (Array)DXGetComponentValue(field, "positions"); if (DXQueryGridPositions(pA, &nd, counts, NULL, deltas)) { int i, j; float l, *d; for (bVol = 1.0, d = deltas, i = 0; i < nd; i++) { if (counts[i] <= 1) { d += nd; } else { for (l = 0.0, j = 0; j < nd; j++, d++) l += *d * *d; bVol *= (counts[i]-1)*sqrt(l); } } } else { bVol = _dxfBoxVolume((Object)field, NULL); if (bVol < 0.0) goto error; } tVol = _dxfFieldVolume(field); if (tVol < 0.0) goto error; if (tVol == 0.0) goto empty; rat = bVol / tVol; if (rat > 0.99) rat = 1.0; nSamples *= rat; deleteGrid = 1; grid = _Sample_Grid((Object)field, nSamples, 2); if (! grid) goto error; } else deleteGrid = 0; a_connections = (Array)DXGetComponentValue(field, "connections"); a_positions = (Array)DXGetComponentValue(field, "positions" ); DXGetArrayInfo(a_connections, &n_connections, NULL, NULL, NULL, NULL); DXGetArrayInfo(a_positions, &n_positions, NULL, NULL, NULL, &pDim); size = nSamples < 100 ? nSamples : nSamples / 10; pList = DXNewSegList(pDim*sizeof(float), size, size); iList = DXNewSegList(sizeof(int)+vPerE*sizeof(float), size, size); if (! pList || ! iList) goto error; if (! DXQueryGridPositions(grid, &gDim, gCounts, mB, mA)) goto error; if (! InvertRegular(gDim, mB, mA, mBinv, mAinv)) goto error; /* * DXAllocate and initialize the grid-status array. */ grid_status = (char *)DXAllocateZero(gCounts[0]*gCounts[1]*sizeof(char)); if (grid_status == NULL) goto error; pHandle = DXCreateArrayHandle(a_positions); if (! pHandle) goto error; if (DXGetComponentValue(field, "invalid connections")) { iHandle = DXCreateInvalidComponentHandle((Object)field, NULL, "connections"); if (iHandle == NULL) goto error; } cHandle = DXCreateArrayHandle(a_connections); if (! cHandle) goto error; nOut = 0; for (i = 0; i < n_connections; i++) { elt = (int *)DXGetArrayEntry(cHandle, i, cBuf); if (iHandle && DXIsElementInvalid(iHandle, i)) continue; /* * Determine the extents of the bounding box of the element */ mint[0] = mint[1] = DXD_MAX_FLOAT; maxt[0] = maxt[1] = -DXD_MAX_FLOAT; for (j = 0; j < vPerE; j++) { float *p; int ii, jj; p = (float *)DXGetArrayEntry(pHandle, *elt++, pBuf+j*gDim); for (ii = 0; ii < 2; ii++) { tBuf[j][ii] = mBinv[ii]; for (jj = 0; jj < pDim; jj++) tBuf[j][ii] += mAinv[jj*gDim+ii]*p[jj]; } eltPoints[j] = &tBuf[j][0]; p = eltPoints[j]; if (p[0] < mint[0]) mint[0] = p[0]; if (p[0] > maxt[0]) maxt[0] = p[0]; if (p[1] < mint[1]) mint[1] = p[1]; if (p[1] > maxt[1]) maxt[1] = p[1]; } mintrel[0] = mint[0]; mintrel[1] = mint[1]; maxtrel[0] = maxt[0]; maxtrel[1] = maxt[1]; for (j = 0; j < 2; j++) { mintgrid[j] = (int)(mintrel[j] - 0.1); if (mintgrid[j] < 0) mintgrid[j] = 0; else if (mintgrid[j] >= gCounts[j]) mintgrid[j] = gCounts[j]-1; maxtgrid[j] = ceil(maxtrel[j] + 0.1); if (maxtgrid[j] < 0) maxtgrid[j] = 0; else if (maxtgrid[j] >= gCounts[j]) maxtgrid[j] = gCounts[j]-1; } basex = mintgrid[0] * gCounts[1]; for (j = mintgrid[0]; j <= maxtgrid[0]; j++) { grid_pt = basex + mintgrid[1]; for (k = mintgrid[1]; k <= maxtgrid[1]; k++) { if (! grid_status[grid_pt]) { float weights[10]; gridy[0] = j; gridy[1] = k; if ((*inside)(gridy, (float **)eltPoints, weights)) { float *pptr = (float *)DXNewSegListItem(pList); Interp *iptr = (Interp *)DXNewSegListItem(iList); int ii; if (! pptr || ! iptr) goto error; grid_status[grid_pt] = 1; if (pDim == 3) { pptr[0] = mA[0*gDim+0]*gridy[0] + mA[1*gDim+0]*gridy[1] + mB[0]; pptr[1] = mA[0*gDim+1]*gridy[0] + mA[1*gDim+1]*gridy[1] + mB[1]; pptr[2] = mA[0*gDim+2]*gridy[0] + mA[1*gDim+2]*gridy[1] + mB[2]; } else { pptr[0] = mA[0*gDim+0]*gridy[0] + mA[1*gDim+0]*gridy[1] + mB[0]; pptr[1] = mA[0*gDim+1]*gridy[0] + mA[1*gDim+1]*gridy[1] + mB[1]; } iptr->index = i; for (ii=0; iiweights[ii] = weights[ii]; } } grid_pt ++; } basex += gCounts[1]; } } DXFree((Pointer) grid_status); grid_status = NULL; if (DXGetSegListItemCount(pList) != 0) { if (! CreateOutputIrreg(field, pList, iList)) goto error; goto done; } empty: while (DXGetEnumeratedComponentValue(field, 0, &name)) DXDeleteComponent(field, name); done: if (pHandle) DXFreeArrayHandle(pHandle); if (cHandle) DXFreeArrayHandle(cHandle); if (iHandle) DXFreeInvalidComponentHandle(iHandle); if (iList) DXDeleteSegList(iList); if (pList) DXDeleteSegList(pList); if (deleteGrid) DXDelete((Object)grid); return OK; error: if (pHandle) DXFreeArrayHandle(pHandle); if (cHandle) DXFreeArrayHandle(cHandle); if (iHandle) DXFreeInvalidComponentHandle(iHandle); if (iList) DXDeleteSegList(iList); if (pList) DXDeleteSegList(pList); if (deleteGrid) DXDelete((Object)grid); DXFree ((Pointer)grid_status); return ERROR; } static Error CurvedSurface(Field field, int vPerE, int (*inside)(), int nSamples, Array grid) { Array a_connections, a_positions; int n_connections, n_positions; int nOut; float element[MAX_V_PER_E][3]; float element2D[MAX_V_PER_E][2], *elt2D[MAX_V_PER_E]; int *connections, cBuf[32]; float gridx[2], gridy[2]; float mint[2], maxt[2]; int mintgrid[2], maxtgrid[2]; int i, j, k, l; float gDeltas[9], gOrigin[3]; int gCounts[3], gSizes[3]; int fCounts[2]; float fOrigin[2]; int deleteGrid; float matrix[9], m; float area, maxArea; int p_axis, axis0, axis1; float v0[3], v1[3], normal[3], a, b, c, d, *p; SegList *pList = NULL, *iList = NULL; int size; ArrayHandle pHandle = NULL; ArrayHandle cHandle = NULL; InvalidComponentHandle iHandle = NULL; char *name; size = nSamples < 100 ? nSamples : nSamples / 10; pList = DXNewSegList(3*sizeof(float), size, size); iList = DXNewSegList(sizeof(int)+vPerE*sizeof(float), size, size); if (! pList || ! iList) goto error; for (i = 0; i < vPerE; i++) elt2D[i] = &element2D[i][0]; if (! grid) { deleteGrid = 1; grid = _Sample_Grid((Object)field, nSamples, 2); if (! grid) goto error; } else deleteGrid = 0; a_connections = (Array)DXGetComponentValue(field, "connections"); a_positions = (Array)DXGetComponentValue(field, "positions" ); DXGetArrayInfo(a_connections, &n_connections, NULL, NULL, NULL, NULL); DXGetArrayInfo(a_positions, &n_positions, NULL, NULL, NULL, NULL); if (! DXQueryGridPositions(grid, NULL, gCounts, gOrigin, gDeltas)) goto error; gSizes[2] = 1; gSizes[1] = gCounts[2]; gSizes[0] = gCounts[0]*gSizes[1]; if (! invert33(gDeltas, matrix)) goto error; pHandle = DXCreateArrayHandle(a_positions); if (! pHandle) goto error; if (DXGetComponentValue(field, "invalid connections")) { iHandle = DXCreateInvalidComponentHandle((Object)field, NULL, "connections"); if (iHandle == NULL) goto error; } cHandle = DXCreateArrayHandle(a_connections); if (! cHandle) goto error; nOut = 0; for (i = 0; i < n_connections; i++) { if (iHandle && DXIsElementInvalid(iHandle, i)) { connections += vPerE; continue; } connections = (int *)DXGetArrayEntry(cHandle, i, cBuf); /* * DXTransform element to canonical coordinates: unit steps, * origin at grid origin. */ for (j = 0; j < vPerE; j++) { float *p, p0[3], pBuf[MAX_V_PER_E]; p = (float *)DXGetArrayEntry(pHandle, *connections++, (Pointer)pBuf); p0[0] = p[0] - gOrigin[0]; p0[1] = p[1] - gOrigin[1]; p0[2] = p[2] - gOrigin[2]; element[j][0] = p0[0]*matrix[0]+p0[1]*matrix[3]+p0[2]*matrix[6]; element[j][1] = p0[0]*matrix[1]+p0[1]*matrix[4]+p0[2]*matrix[7]; element[j][2] = p0[0]*matrix[2]+p0[1]*matrix[5]+p0[2]*matrix[8]; } /* * Determine which face element is most perpendicular to. To * do so, compute planar equation. */ v0[0] = element[1][0] - element[0][0]; v0[1] = element[1][1] - element[0][1]; v0[2] = element[1][2] - element[0][2]; v1[0] = element[2][0] - element[0][0]; v1[1] = element[2][1] - element[0][1]; v1[2] = element[2][2] - element[0][2]; if (vector_zero((Vector *)v0) || vector_zero((Vector *)v1)) continue; normal[0] = v0[1]*v1[2] - v1[1]*v0[2]; normal[1] = v0[2]*v1[0] - v1[2]*v0[0]; normal[2] = v0[0]*v1[1] - v1[0]*v0[1]; if (vector_zero((Vector *)normal)) continue; if ((a = normal[0]) < 0.0) a = -a; if ((b = normal[1]) < 0.0) b = -b; if ((c = normal[2]) < 0.0) c = -c; if (a >= b && a >= c) { /* * use the yz plane */ p_axis = 0; axis0 = 1; axis1 = 2; } else if (b >= c) { /* * use the xz plane */ axis0 = 0; p_axis = 1; axis1 = 2; } else { /* * use the xy plane */ axis0 = 0; axis1 = 1; p_axis = 2; } a = normal[0]; b = normal[1]; c = normal[2]; d = -(a*element[0][0] + b*element[0][1] + c*element[0][2]); fCounts[0] = gCounts[axis0]; fCounts[1] = gCounts[axis1]; fOrigin[0] = gOrigin[axis0]; fOrigin[1] = gOrigin[axis1]; mint[0] = mint[1] = DXD_MAX_FLOAT; maxt[0] = maxt[1] = -DXD_MAX_FLOAT; /* * Project onto selected axis and subtract off the face origin * This leaves us in a face-space in which the (0,0) grid point * is at (0.0, 0.0). */ for (j = 0; j < vPerE; j++) { element2D[j][0] = element[j][axis0]; element2D[j][1] = element[j][axis1]; p = &element2D[j][0]; if (p[0] < mint[0]) mint[0] = p[0]; if (p[0] > maxt[0]) maxt[0] = p[0]; if (p[1] < mint[1]) mint[1] = p[1]; if (p[1] > maxt[1]) maxt[1] = p[1]; } for (j = 0; j < 2; j++) { mintgrid[j] = (int)mint[j]; if (mintgrid[j] < 0) mintgrid[j] = 0; else if (mintgrid[j] >= fCounts[j]) mintgrid[j] = fCounts[j] - 1; maxtgrid[j] = (int)maxt[j]; if (maxtgrid[j] < 0) maxtgrid[j] = 0; else if (maxtgrid[j] >= fCounts[j]) maxtgrid[j] = fCounts[j] - 1; } gridx[0] = mintgrid[0]; gridx[1] = 0.0; for (j = mintgrid[0]; j <= maxtgrid[0]; j++) { gridy[0] = gridx[0]; gridy[1] = gridx[1] + mintgrid[1]; for (k = mintgrid[1]; k <= maxtgrid[1]; k++) { float g_pt[MAX_V_PER_E]; float weights[MAX_V_PER_E]; if ((*inside)(gridy, elt2D, weights)) { float *pptr = (float *)DXNewSegListItem(pList); Interp *iptr = (Interp *)DXNewSegListItem(iList); int ii; if (! pptr || ! iptr) goto error; switch(p_axis) { case 0: g_pt[1] = gridy[0]; g_pt[2] = gridy[1]; g_pt[0] = (b*g_pt[1] + c*g_pt[2] + d) / -a; break; case 1: g_pt[0] = gridy[0]; g_pt[2] = gridy[1]; g_pt[1] = (a*g_pt[0] + c*g_pt[2] + d) / -b; break; case 2: g_pt[0] = gridy[0]; g_pt[1] = gridy[1]; g_pt[2] = (a*g_pt[0] + b*g_pt[1] + d) / -c; break; } pptr[0] = g_pt[0]*gDeltas[0] + g_pt[1]*gDeltas[3] + g_pt[2]*gDeltas[6] + gOrigin[0]; pptr[1] = g_pt[0]*gDeltas[1] + g_pt[1]*gDeltas[4] + g_pt[2]*gDeltas[7] + gOrigin[1]; pptr[2] = g_pt[0]*gDeltas[2] + g_pt[1]*gDeltas[5] + g_pt[2]*gDeltas[8] + gOrigin[2]; iptr->index = i; for (ii = 0; ii < vPerE; ii++) iptr->weights[ii] = weights[ii]; } gridy[1] += 1.0; } gridx[0] += 1.0; } } if (DXGetSegListItemCount(pList) != 0) { if (! CreateOutputIrreg(field, pList, iList)) goto error; goto done; } empty: while(DXGetEnumeratedComponentValue(field, 0, &name)) DXDeleteComponent(field, name); done: if (pHandle) DXFreeArrayHandle(pHandle); if (cHandle) DXFreeArrayHandle(cHandle); if (iHandle) DXFreeInvalidComponentHandle(iHandle); if (pList) DXDeleteSegList(pList); if (iList) DXDeleteSegList(iList); if (deleteGrid) DXDelete((Object)grid); return OK; error: if (pHandle) DXFreeArrayHandle(pHandle); if (cHandle) DXFreeArrayHandle(cHandle); if (iHandle) DXFreeInvalidComponentHandle(iHandle); if (pList) DXDeleteSegList(pList); if (iList) DXDeleteSegList(iList); if (deleteGrid) DXDelete((Object)grid); return ERROR; } static Error Volume(Field field, int vPerE, int (*inside)(), int nSamples, Array grid) { Array a_connections, a_positions; int n_connections, n_positions; int nOut; float *element[MAX_V_PER_E]; int *connections; int yz, xyz; int basex, basey; float gridx[3], gridy[3], gridz[3]; char *grid_status = NULL; int grid_pt; float mint[3], maxt[3]; float mintrel[3], maxtrel[3]; int mintgrid[3], maxtgrid[3]; int i, j, k, l; float gDeltas[9], gOrigin[3]; int gCounts[3]; int deleteGrid; float matrix[9]; SegList *pList = NULL, *iList = NULL; int size; ArrayHandle pHandle = NULL, cHandle = NULL; float pBuf[MAX_V_PER_E][3]; int eBuf[MAX_V_PER_E], *elt; InvalidComponentHandle iHandle = NULL; char *name; size = nSamples < 100 ? nSamples : nSamples / 10; pList = DXNewSegList(3*sizeof(float), size, size); iList = DXNewSegList(sizeof(int)+vPerE*sizeof(float), size, size); if (! pList || ! iList) goto error; if (! grid) { float bVol, tVol, rat; deleteGrid = 1; bVol = _dxfBoxVolume((Object)field, NULL); if (bVol < 0.0) goto error; tVol = _dxfFieldVolume(field); if (tVol < 0.0) goto error; if (tVol == 0.0) goto empty; rat = bVol / tVol; if (rat > 0.99) rat = 1.0; nSamples *= rat; grid = _Sample_Grid((Object)field, nSamples, 3); if (! grid) goto error; } else deleteGrid = 0; a_connections = (Array)DXGetComponentValue(field, "connections"); a_positions = (Array)DXGetComponentValue(field, "positions" ); DXGetArrayInfo(a_connections, &n_connections, NULL, NULL, NULL, NULL); DXGetArrayInfo(a_positions, &n_positions, NULL, NULL, NULL, NULL); if (! DXQueryGridPositions(grid, NULL, gCounts, gOrigin, gDeltas)) goto error; yz = gCounts[1] * gCounts[2]; xyz = gCounts[0] * yz; if (! invert33(gDeltas, matrix)) goto error; /* * DXAllocate and initialize the grid-status array. */ grid_status = (char *)DXAllocateZero(xyz * sizeof(char)); if (grid_status == NULL) goto error; pHandle = DXCreateArrayHandle(a_positions); cHandle = DXCreateArrayHandle(a_connections); if (! pHandle || ! cHandle) goto error; if (DXGetComponentValue(field, "invalid connections")) { iHandle = DXCreateInvalidComponentHandle((Object)field, NULL, "connections"); if (! iHandle) goto error; } nOut = 0; for (i = 0; i < n_connections; i++) { if (iHandle && DXIsElementInvalid(iHandle, i)) continue; elt = (int *)DXGetArrayEntry(cHandle, i, (Pointer)eBuf); /* * Determine the extents of the bounding box of the element. */ mint[0] = mint[1] = mint[2] = DXD_MAX_FLOAT; maxt[0] = maxt[1] = maxt[2] = -DXD_MAX_FLOAT; for (j = 0; j < vPerE; j++) { float *p; element[j] = (float *)DXGetArrayEntry(pHandle, elt[j], (Pointer)&pBuf[j][0]); p = element[j]; if (*p < mint[0]) mint[0] = *p; if (*p > maxt[0]) maxt[0] = *p; p++; if (*p < mint[1]) mint[1] = *p; if (*p > maxt[1]) maxt[1] = *p; p++; if (*p < mint[2]) mint[2] = *p; if (*p > maxt[2]) maxt[2] = *p; } for (j = 0; j < 3; j++) { mint[j] = mint[j] - gOrigin[j]; maxt[j] = maxt[j] - gOrigin[j]; } mintrel[0] = mint[0]*matrix[0]+mint[1]*matrix[3]+mint[2]*matrix[6]; mintrel[1] = mint[0]*matrix[1]+mint[1]*matrix[4]+mint[2]*matrix[7]; mintrel[2] = mint[0]*matrix[2]+mint[1]*matrix[5]+mint[2]*matrix[8]; maxtrel[0] = maxt[0]*matrix[0]+maxt[1]*matrix[3]+maxt[2]*matrix[6]; maxtrel[1] = maxt[0]*matrix[1]+maxt[1]*matrix[4]+maxt[2]*matrix[7]; maxtrel[2] = maxt[0]*matrix[2]+maxt[1]*matrix[5]+maxt[2]*matrix[8]; for (j = 0; j < 3; j++) if (mintrel[j] > maxtrel[j]) { float t = maxtrel[j]; maxtrel[j] = mintrel[j]; mintrel[j] = t; } for (j = 0; j < 3; j++) { mintgrid[j] = floor(mintrel[j]-0.1); if (mintgrid[j] < 0) mintgrid[j] = 0; else if (mintgrid[j] > gCounts[j] - 1) mintgrid[j] = gCounts[j] - 1; maxtgrid[j] = ceil(maxtrel[j]+0.1); if (maxtgrid[j] < 0) maxtgrid[j] = 0; else if (maxtgrid[j] > gCounts[j] - 1) maxtgrid[j] = gCounts[j] - 1; } basex = mintgrid[0] * yz; gridx[0] = gOrigin[0] + mintgrid[0]*gDeltas[0]; gridx[1] = gOrigin[1] + mintgrid[0]*gDeltas[1]; gridx[2] = gOrigin[2] + mintgrid[0]*gDeltas[2]; for (j = mintgrid[0]; j <= maxtgrid[0]; j++) { basey = basex + mintgrid[1] * gCounts[2]; gridy[0] = gridx[0] + mintgrid[1]*gDeltas[3]; gridy[1] = gridx[1] + mintgrid[1]*gDeltas[4]; gridy[2] = gridx[2] + mintgrid[1]*gDeltas[5]; for (k = mintgrid[1]; k <= maxtgrid[1]; k++) { grid_pt = basey + mintgrid[2]; gridz[0] = gridy[0] + mintgrid[2]*gDeltas[6]; gridz[1] = gridy[1] + mintgrid[2]*gDeltas[7]; gridz[2] = gridy[2] + mintgrid[2]*gDeltas[8]; for (l = mintgrid[2]; l <= maxtgrid[2]; l++) { /* * If we have already determined the element which * contains the grid point, continue. */ if (! grid_status[grid_pt]) { int ii; float weights[MAX_V_PER_E]; if ((*inside)(gridz, element, weights)) { float *pptr = (float *)DXNewSegListItem(pList); Interp *iptr = (Interp *)DXNewSegListItem(iList); int j; grid_status[grid_pt] = 1; if (! pptr || ! iptr) goto error; pptr[0] = gridz[0]; pptr[1] = gridz[1]; pptr[2] = gridz[2]; iptr->index = i; memcpy(&iptr->weights[0], weights, vPerE*sizeof(float)); } } grid_pt ++; gridz[0] += gDeltas[6]; gridz[1] += gDeltas[7]; gridz[2] += gDeltas[8]; } basey += gCounts[2]; gridy[0] += gDeltas[3]; gridy[1] += gDeltas[4]; gridy[2] += gDeltas[5]; } basex += yz; gridx[0] += gDeltas[0]; gridx[1] += gDeltas[1]; gridx[2] += gDeltas[2]; } } DXFree((Pointer) grid_status); grid_status = NULL; if (DXGetSegListItemCount(pList) != 0) { if (! CreateOutputIrreg(field, pList, iList)) goto error; goto done; } empty: while (DXGetEnumeratedComponentValue(field, 0, &name)) DXDeleteComponent(field, name); done: if (pHandle) DXFreeArrayHandle(pHandle); if (cHandle) DXFreeArrayHandle(cHandle); if (iHandle) DXFreeInvalidComponentHandle(iHandle); if (iList) DXDeleteSegList(iList); if (pList) DXDeleteSegList(pList); if (deleteGrid) DXDelete((Object)grid); return OK; error: if (pHandle) DXFreeArrayHandle(pHandle); if (cHandle) DXFreeArrayHandle(cHandle); if (iHandle) DXFreeInvalidComponentHandle(iHandle); if (iList) DXDeleteSegList(iList); if (pList) DXDeleteSegList(pList); if (deleteGrid) DXDelete((Object)grid); DXFree ((Pointer)grid_status); return ERROR; } typedef struct { int index; int i0, i1; float weights[2]; } PolylineInterp; static Error PolylinesSample(Field field, int nSamples, Array grid) { int i, j, k, nDim, nPts, nPolylines, nEdges; ArrayHandle points = NULL; InvalidComponentHandle ich = NULL; Array pA, plA, eA; float length, totLength; float spacing, next; char *name; int *polylines, *edges; float *lengths; SegList *pList = NULL; SegList *iList = NULL; if (nSamples == 0) goto empty; if (DXEmptyField(field)) goto error; pA = (Array)DXGetComponentValue(field, "positions"); if (! pA) goto error; if (! DXGetArrayInfo(pA, &nPts, NULL, NULL, NULL, &nDim)) goto error; points = DXCreateArrayHandle(pA); if (! points) goto error; plA = (Array)DXGetComponentValue(field, "polylines"); if (! plA) goto error; if (! DXGetArrayInfo(plA, &nPolylines, NULL, NULL, NULL, NULL)) goto error; eA = (Array)DXGetComponentValue(field, "edges"); if (! plA) goto error; if (! DXGetArrayInfo(eA, &nEdges, NULL, NULL, NULL, NULL)) goto error; polylines = (int *)DXGetArrayData(plA); edges = (int *)DXGetArrayData(eA); if (DXGetComponentValue(field, "invalid polylines")) { ich = DXCreateInvalidComponentHandle((Object)field, NULL, "polylines"); if (! ich) goto error; } totLength = 0.0; for (i = 0; i < nPolylines; i++) { float *p0, p0Buf[10]; float *p1, p1Buf[10]; int start, end, knt, *seg; if (ich && DXIsElementInvalid(ich, i)) continue; start = polylines[i]; end = (i == nPolylines-1) ? nEdges-1 : polylines[i+1]-1; knt = end - start; seg = edges + start; p1 = (float *)DXGetArrayEntry(points, *seg, p0Buf); for (j = 0, seg++; j < knt; j++, seg++) { float l; p0 = p1; p1 = (float *)DXGetArrayEntry(points, *seg, ((j & 0x1) ? p0Buf : p1Buf)); for (l = 0.0, k = 0; k < nDim; k++) { float d = p1[k] - p0[k]; l += d*d; } totLength += sqrt(l); } } if (totLength == 0) goto empty; pList = DXNewSegList(nDim*sizeof(float), nSamples, 10); iList = DXNewSegList(sizeof(PolylineInterp), nSamples, 10); if (! pList || ! iList) goto error; spacing = totLength / nSamples; length = 0; next = spacing / 2; totLength = 0.0; for (i = 0; i < nPolylines; i++) { float *p0, p0Buf[10]; float *p1, p1Buf[10]; int start, end, knt, *seg; float l0, l1, l; int p, q; if (ich && DXIsElementInvalid(ich, i)) continue; start = polylines[i]; end = (i == nPolylines-1) ? nEdges-1 : polylines[i+1]-1; knt = end - start; seg = edges + start; q = *seg; p1 = (float *)DXGetArrayEntry(points, q, p0Buf); for (j = 0; j < knt; j++, totLength += l) { p = q; q = *(++seg); p0 = p1; p1 = (float *)DXGetArrayEntry(points, q, ((j & 0x1) ? p0Buf : p1Buf)); for (l = 0.0, k = 0; k < nDim; k++) { float d = p1[k] - p0[k]; l += d*d; } if (l == 0.0) continue; l = sqrt(l); l0 = totLength; l1 = totLength + l; while (next < l1) { float w1 = (next - l0) / (l1 - l0); float w0 = 1.0 - w1; float *pptr=(float *)DXNewSegListItem(pList); PolylineInterp *iptr=(PolylineInterp *)DXNewSegListItem(iList); iptr->index = i; iptr->i0 = p; iptr->i1 = q; iptr->weights[0] = w0; iptr->weights[1] = w1; for (k = 0; k < nDim; k++) *pptr++ = w0*p0[k] + w1*p1[k]; next += spacing; } } } if (DXGetSegListItemCount(pList) != 0) { if (! CreateOutputPolylines(field, pList, iList)) goto error; goto done; } empty: while (DXGetEnumeratedComponentValue(field, 0, &name)) DXDeleteComponent(field, name); done: if (ich) DXFreeInvalidComponentHandle(ich); DXFreeArrayHandle(points); DXDeleteSegList(pList); DXDeleteSegList(iList); return OK; error: if (ich) DXFreeInvalidComponentHandle(ich); DXFreeArrayHandle(points); DXDeleteSegList(pList); DXDeleteSegList(iList); return ERROR; } typedef struct { int index; float weights[2]; } LineInterp; struct segment { struct segment *p0, *p1; byte visited, invalid, reverse; int seg; float length; }; static Error LineSample(Field field, int nSamples, Array grid) { int i, j, nDim, nPts, nSegs; struct segment *segment = NULL; struct segment **hash = NULL; SegList *pList = NULL; SegList *iList = NULL; ArrayHandle points = NULL; ArrayHandle elements = NULL; InvalidComponentHandle ich = NULL; Array pA, cA; float length, totLength, invLength; float spacing, next; char *name; if (DXEmptyField(field)) goto error; pA = (Array)DXGetComponentValue(field, "positions"); if (! pA) goto error; if (! DXGetArrayInfo(pA, &nPts, NULL, NULL, NULL, &nDim)) goto error; points = DXCreateArrayHandle(pA); if (! points) goto error; cA = (Array)DXGetComponentValue(field, "connections"); if (! cA) goto error; if (! DXGetArrayInfo(cA, &nSegs, NULL, NULL, NULL, NULL)) goto error; elements = DXCreateArrayHandle(cA); if (! elements) goto error; segment = (struct segment *)DXAllocate(nSegs * sizeof(struct segment)); if (! segment) goto error; hash = (struct segment **)DXAllocateZero(nPts * sizeof(struct segment *)); if (! hash) goto error; pList = DXNewSegList(nDim*sizeof(float), nSamples, 10); iList = DXNewSegList(sizeof(int)+2*sizeof(float), nSamples, 10); if (! pList || ! iList) goto error; ich = DXCreateInvalidComponentHandle((Object)field, NULL, "connections"); if (! ich) goto error; totLength = invLength = 0.0; for (i = 0; i < nSegs; i++) { struct segment *s; int *elt, eBuf[2]; float *p0, p0Buf[10]; float *p1, p1Buf[10]; float l; segment[i].visited = 0; segment[i].seg = i; elt = (int *)DXGetArrayEntry(elements, i, eBuf); if (elt[0] == elt[1]) { segment[i].p0 = NULL; segment[i].p1 = NULL; segment[i].invalid = 1; continue; } p0 = (float *)DXGetArrayEntry(points, elt[0], p0Buf); p1 = (float *)DXGetArrayEntry(points, elt[1], p1Buf); segment[i].invalid = DXIsElementInvalid(ich, i); for (l = 0.0, j = 0; j < nDim; j++) { float d = p1[j] - p0[j]; l += d*d; } segment[i].length = sqrt(l); if ((s = hash[elt[0]]) == NULL) { hash[elt[0]] = segment + i; segment[i].p0 = NULL; } else { int *nbr, nBuf[2]; nbr = DXGetArrayEntry(elements, s->seg, nBuf); if (nbr[0] == elt[0]) s->p0 = segment + i; else s->p1 = segment + i; segment[i].p0 = s; hash[elt[0]] = NULL; } if ((s = hash[elt[1]]) == NULL) { hash[elt[1]] = segment + i; segment[i].p1 = NULL; } else { int *nbr, nBuf[2]; nbr = DXGetArrayEntry(elements, s->seg, nBuf); if (nbr[0] == elt[1]) s->p0 = segment + i; else s->p1 = segment + i; segment[i].p1 = s; hash[elt[1]] = NULL; } totLength += segment[i].length; if (segment[i].invalid) invLength += segment[i].length; } if (invLength == totLength || nSamples == 0) goto empty; nSamples = nSamples * totLength / (totLength - invLength); spacing = (totLength - invLength) / nSamples; length = 0; next = spacing / 2; for (i = 0; i < nSegs; i++) if (segment[i].visited == 0 && !segment[i].invalid) { int first; struct segment *p, *start, *l; /* * Trace back as far as possible */ l = NULL; start = p = &segment[i]; do { if (l == NULL) /* then go either way */ { l = p; p = p->p0; } else if (l == p->p0) { l = p; p = p->p1; } else { l = p; p = p->p0; } } while (p && (p != start)); /* * Now forward setting reverse flag */ start = p = l; l = NULL; do { if (l == NULL) /* then go either way */ { l = p; if (p->p0) { p->reverse = 1; p = p->p0; } else { p->reverse = 0; p = p->p1; } } else if (l == p->p0) { l = p; p->reverse = 0; p = p->p1; } else { l = p; p->reverse = 1; p = p->p0; } } while (p && (p != start)); p = start; do { float l0 = length, l1 = l0 + p->length; p->visited = 1; length += p->length; first = 1; while (next < l1) { int *elt, eBuf[10]; float *p0, p0Buf[10]; float *p1, p1Buf[10]; if (! p->invalid) { float w1 = (next - l0) / (l1 - l0); float w0 = 1.0 - w1; float *pptr=(float *)DXNewSegListItem(pList); LineInterp *iptr=(LineInterp *)DXNewSegListItem(iList); if (first) { elt=(int *)DXGetArrayEntry(elements, p->seg, eBuf); p0=(float *)DXGetArrayEntry(points, elt[0], p0Buf); p1=(float *)DXGetArrayEntry(points, elt[1], p1Buf); first = 0; } if (! pptr || ! iptr) goto error; iptr->index = p->seg; if (p->reverse) { iptr->weights[1] = w0; iptr->weights[0] = w1; for (j = 0; j < nDim; j++) *pptr++ = w1*p0[j] + w0*p1[j]; } else { iptr->weights[0] = w0; iptr->weights[1] = w1; for (j = 0; j < nDim; j++) *pptr++ = w0*p0[j] + w1*p1[j]; } } next += spacing; } if (p->reverse) p = p->p0; else p = p->p1; } while (p && (p != start)); } if (DXGetSegListItemCount(pList) != 0) { if (! CreateOutputIrreg(field, pList, iList)) goto error; goto done; } empty: while (DXGetEnumeratedComponentValue(field, 0, &name)) DXDeleteComponent(field, name); done: DXFree((Pointer)segment); DXFree((Pointer)hash); DXFreeArrayHandle(points); DXFreeArrayHandle(elements); DXDeleteSegList(pList); DXDeleteSegList(iList); DXFreeInvalidComponentHandle(ich); return OK; error: DXFree((Pointer)segment); DXFree((Pointer)hash); DXFreeArrayHandle(points); DXFreeArrayHandle(elements); DXDeleteSegList(pList); DXDeleteSegList(iList); DXFreeInvalidComponentHandle(ich); return ERROR; } static int _invert33(float *a, float *b, float *c, float *dst) { float m00, m01, m02; float m10, m11, m12; float m20, m21, m22; float c00, c01, c02; float c10, c11, c12; float c20, c21, c22; float det; m00 = *a++; m01 = *a++; m02 = *a++; m10 = *b++; m11 = *b++; m12 = *b++; m20 = *c++; m21 = *c++; m22 = *c++; c00 = m11*m22 - m12*m21; c01 = m10*m22 - m20*m12; c02 = m10*m21 - m20*m11; c10 = m01*m22 - m21*m02; c11 = m00*m22 - m02*m20; c12 = m00*m21 - m20*m01; c20 = m01*m12 - m11*m02; c21 = m00*m12 - m10*m02; c22 = m00*m11 - m10*m01; det = m00*c00 - m10*c10 + m20*c20; if (det == 0.0) return 0; det = 1.0 / det; *dst++ = det * c00; *dst++ = -det * c10; *dst++ = det * c20; *dst++ = -det * c01; *dst++ = det * c11; *dst++ = -det * c21; *dst++ = det * c02; *dst++ = -det * c12; *dst++ = det * c22; return 1; } static int invert33(float *src, float *dst) { int i, j, k; int badVector[3], nBad; if (_invert33(src, src+3, src+6, dst)) return 1; /* * OK. Matrix was singular, so we can replace dependant vectors (or * zero vectors) with a arbitrary independent vectors. First, find bad * vectors. * * First possibility is zero vectors. */ nBad = 0; for (i = 0; i < 3; i++) { int zeroVector = 1; for (j = 0; j < 3 && zeroVector; j++) zeroVector = (src[i*3+j] == 0.0); if ((badVector[i] = zeroVector) == 1) nBad ++; } if (nBad == 3) goto error; /* * Now for each non-zero vector, cross it against the remaining * non-zero vectors and look at the length of the result to see if they * are parallel. */ for (i = 0; i < 2; i++) { if (badVector[i]) continue; else { float *v = src + 3*i; for (j = i+1; j < 3 && !badVector[i]; j++) { if (badVector[j]) continue; else { float *w = src + 3*j; if ((v[2]*w[1] == v[1]*w[2]) && (v[0]*w[2] == v[2]*w[0]) && (v[1]*w[0] == v[0]*w[1])) { badVector[i] = 1; nBad ++; } } } } } /* * There are either 0, 1 or two vectors identified as bad. */ switch(nBad) { /* * If there are zero bad vectors, then any one is a combination * of the other two. So we can replace it by the cross. */ case 0: { float *v = src; float *w = src+3; float c[3]; c[0] = v[2]*w[1] - v[1]*w[2]; c[1] = v[0]*w[2] - v[2]*w[0]; c[2] = v[1]*w[0] - v[0]*w[1]; if (! _invert33(v, w, c, dst)) goto error; } break; /* * If there is one bad vector, again, we replace it by the cross * product of the others. */ #define CROSS(a, b, c) \ { \ c[0] = a[2]*b[1] - a[1]*b[2]; \ c[1] = a[0]*b[2] - a[2]*b[0]; \ c[2] = a[1]*b[0] - a[0]*b[1]; \ } case 1: { float *v, *w, c[3]; if (badVector[0]) { v = src+3; w = src+6; CROSS(v, w, c); if (! _invert33(v, c, w, dst)) goto error; } else if (badVector[1]) { v = src; w = src+6; CROSS(v, w, c); if (! _invert33(v, c, w, dst)) goto error; } else { v = src; w = src+3; CROSS(v, w, c); if (! _invert33(v, w, c, dst)) goto error; } } break; /* * If two are bad, then we find the good one, any non-parallel * to it and the cross of the two. */ #define NON_PARALLEL(a, b) \ { \ if (a[0] == 0.0 && a[1] == 0.0) \ { \ b[0] = 1.0; b[1] = 0.0; b[2] = 0.0; \ } \ else \ { \ b[0] = 0.0; b[1] = 0.0; b[2] = 1.0; \ } \ } case 2: { float *v, c0[3], c1[3]; if (! badVector[0]) { v = src; NON_PARALLEL(v, c0); CROSS(v, c0, c1); if (! _invert33(v, c0, c1, dst)) goto error; } else if (! badVector[1]) { v = src+3; NON_PARALLEL(v, c0); CROSS(v, c0, c1); if (! _invert33(c0, v, c1, dst)) goto error; } else { v = src+6; NON_PARALLEL(v, c0); CROSS(v, c0, c1); if (! _invert33(c0, c1, v, dst)) goto error; } } break; } return 1; error: DXSetError(ERROR_INVALID_DATA, "degenerate regular positions"); return 0; } static float TetrahedraVolume(float **points, int pDim) { float *p, *q, *r, *s, t; float pq[3], pr[3], ps[3]; int i; p = points[0]; q = points[1]; r = points[2]; s = points[3]; for (i = 0; i < 3; i++) { pq[i] = q[i] - p[i]; pr[i] = r[i] - p[i]; ps[i] = s[i] - p[i]; } if (vector_zero((Vector *)pq) || vector_zero((Vector *)pr) || vector_zero((Vector *)ps)) return 0.0; t = 0.16666667 * (pq[0] * (pr[1]*ps[2] - pr[2]*ps[1]) - pr[0] * (pq[1]*ps[2] - pq[2]*ps[1]) + ps[0] * (pq[1]*pr[2] - pq[2]*pr[1])); if (t >= 0.0) return t; else return -t; } static int coords_tetrahedra(float *pt, float **tetra, float *weights) { float vol0, vol1, vol2, vol3; float xt, x0, x1, x2, x3; float yt, y0, y1, y2, y3; float zt, z0, z1, z2, z3; float x01, x02, x03, x0t, xt1, xt2, xt3; float y01, y02, y03, y0t, yt1, yt2, yt3; float z01, z02, z03, z0t, zt1, zt2, zt3; float yz012, yz013, yz023, yzt12, yzt13; float yzt23, yz0t2, yz0t3, yz01t, yz02t; xt = pt[0]; yt = pt[1]; zt = pt[2]; x0 = tetra[0][0]; y0 = tetra[0][1]; z0 = tetra[0][2]; x1 = tetra[1][0]; y1 = tetra[1][1]; z1 = tetra[1][2]; x2 = tetra[2][0]; y2 = tetra[2][1]; z2 = tetra[2][2]; x3 = tetra[3][0]; y3 = tetra[3][1]; z3 = tetra[3][2]; x01 = x1 - x0; y01 = y1 - y0; z01 = z1 - z0; x02 = x2 - x0; y02 = y2 - y0; z02 = z2 - z0; x03 = x3 - x0; y03 = y3 - y0; z03 = z3 - z0; x0t = xt - x0; y0t = yt - y0; z0t = zt - z0; xt1 = x1 - xt; yt1 = y1 - yt; zt1 = z1 - zt; xt2 = x2 - xt; yt2 = y2 - yt; zt2 = z2 - zt; xt3 = x3 - xt; yt3 = y3 - yt; zt3 = z3 - zt; yz012 = y01 * z02 - y02 * z01; yz013 = y01 * z03 - y03 * z01; yz023 = y02 * z03 - y03 * z02; yzt12 = yt1 * zt2 - yt2 * zt1; yzt13 = yt1 * zt3 - yt3 * zt1; yzt23 = yt2 * zt3 - yt3 * zt2; yz0t2 = y0t * z02 - y02 * z0t; yz0t3 = y0t * z03 - y03 * z0t; yz01t = y01 * z0t - y0t * z01; yz02t = y02 * z0t - y0t * z02; vol0 = xt1 * yzt23 - xt2 * yzt13 + xt3 * yzt12; vol1 = x0t * yz023 - x02 * yz0t3 + x03 * yz0t2; vol2 = x01 * yz0t3 - x0t * yz013 + x03 * yz01t; vol3 = x01 * yz02t - x02 * yz01t + x0t * yz012; if (vol0 >= 0 && vol1 >= 0 && vol2 >= 0 && vol3 >= 0 || vol0 <= 0 && vol1 <= 0 && vol2 <= 0 && vol3 <= 0) { if (weights) { float d, v = vol0 + vol1 + vol2 + vol3; if (v != 0.0) d = 1.0 / v; else d = 0.0; weights[0] = vol0 * d; weights[1] = vol1 * d; weights[2] = vol2 * d; weights[3] = vol3 * d; } return 1; } return 0; } static float CubeVolume(float **points, int pDim) { float *tetra[4]; float vol; tetra[0] = points[2]; tetra[1] = points[3]; tetra[2] = points[7]; tetra[3] = points[1]; vol = TetrahedraVolume(tetra, pDim); tetra[0] = points[7]; tetra[1] = points[1]; tetra[2] = points[5]; tetra[3] = points[4]; vol += TetrahedraVolume(tetra, pDim); tetra[0] = points[7]; tetra[1] = points[2]; tetra[2] = points[6]; tetra[3] = points[4]; vol += TetrahedraVolume(tetra, pDim); tetra[0] = points[0]; tetra[1] = points[1]; tetra[2] = points[4]; tetra[3] = points[2]; vol += TetrahedraVolume(tetra, pDim); tetra[0] = points[2]; tetra[1] = points[7]; tetra[2] = points[1]; tetra[3] = points[4]; vol += TetrahedraVolume(tetra, pDim); return vol; } static int coords_cube(float *pt, float **cube, float *weights) { float *tet[4]; float tetweights[4]; if (weights) memset(weights, 0, 8*sizeof(float)); tet[0] = cube[2]; tet[1] = cube[3]; tet[2] = cube[7]; tet[3] = cube[1]; if (coords_tetrahedra(pt, tet, tetweights)) { if (weights) { weights[2] = tetweights[0]; weights[3] = tetweights[1]; weights[7] = tetweights[2]; weights[1] = tetweights[3]; } return 1; } tet[0] = cube[7]; tet[1] = cube[1]; tet[2] = cube[5]; tet[3] = cube[4]; if (coords_tetrahedra(pt, tet, tetweights)) { if (weights) { weights[7] = tetweights[0]; weights[1] = tetweights[1]; weights[5] = tetweights[2]; weights[4] = tetweights[3]; } return 1; } tet[0] = cube[7]; tet[1] = cube[2]; tet[2] = cube[6]; tet[3] = cube[4]; if (coords_tetrahedra(pt, tet, tetweights)) { if (weights) { weights[7] = tetweights[0]; weights[2] = tetweights[1]; weights[6] = tetweights[2]; weights[4] = tetweights[3]; } return 1; } tet[0] = cube[0]; tet[1] = cube[1]; tet[2] = cube[4]; tet[3] = cube[2]; if (coords_tetrahedra(pt, tet, tetweights)) { if (weights) { weights[0] = tetweights[0]; weights[1] = tetweights[1]; weights[4] = tetweights[2]; weights[2] = tetweights[3]; } return 1; } tet[0] = cube[2]; tet[1] = cube[7]; tet[2] = cube[1]; tet[3] = cube[4]; if (coords_tetrahedra(pt, tet, tetweights)) { if (weights) { weights[2] = tetweights[0]; weights[7] = tetweights[1]; weights[1] = tetweights[2]; weights[4] = tetweights[3]; } return 1; } return 0; } static float TriangleArea(float **points, int pDim) { float *p0, *p1, *p2; float v0x, v0y, v0z; float v1x, v1y, v1z; float a; p0 = points[0]; p1 = points[1]; p2 = points[2]; if (pDim == 2) { v0x = p1[0] - p0[0]; v0y = p1[1] - p0[1]; v1x = p2[0] - p0[0]; v1y = p2[1] - p0[1]; a = 0.5 * (v0x*v1y - v0y*v1x); if (a >= 0) return a; else return -a; } else { Vector v0, v1, cross; vector_subtract((Vector *)p0, (Vector *)p1, &v0); vector_subtract((Vector *)p0, (Vector *)p2, &v1); vector_cross(&v0, &v1, &cross); a = 0.5 * vector_length(&cross); } if (a >= 0) return a; else return -a; } static int coords_triangle(float *pt, float **triangle, float *weights) { float *p0, *p1, *p2; double v0x, v0y; double v1x, v1y; double a0, a1, a2; float d, a; p0 = triangle[0]; p1 = triangle[1]; p2 = triangle[2]; v0x = p1[0] - pt[0]; v0y = p1[1] - pt[1]; v1x = p2[0] - pt[0]; v1y = p2[1] - pt[1]; a0 = v0x*v1y - v0y*v1x; v0x = p2[0] - pt[0]; v0y = p2[1] - pt[1]; v1x = p0[0] - pt[0]; v1y = p0[1] - pt[1]; a1 = v0x*v1y - v0y*v1x; v0x = p0[0] - pt[0]; v0y = p0[1] - pt[1]; v1x = p1[0] - pt[0]; v1y = p1[1] - pt[1]; a2 = v0x*v1y - v0y*v1x; a = a0 + a1 + a2; if (a == 0.0) return 0; else d = 1.0 / a; a0 = a0 * d; a1 = a1 * d; a2 = a2 * d; if (a0 < -0.0001 || a1 < -0.0001 || a2 < -0.0001) return 0; if (weights) { weights[0] = a0; weights[1] = a1; weights[2] = a2; } return 1; } static float QuadArea(float **points, int pDim) { float *tri0[3], *tri1[3]; tri0[0] = points[0]; tri0[1] = points[1]; tri0[2] = points[2]; tri1[0] = points[0]; tri1[1] = points[2]; tri1[2] = points[3]; return TriangleArea(tri0, pDim) + TriangleArea(tri1, pDim); } static int coords_quad(float *pt, float **quad, float *weights) { float *tri[3]; float tweights[3]; tri[0] = quad[0]; tri[1] = quad[1]; tri[2] = quad[2]; if (coords_triangle(pt, tri, tweights)) { if (weights) { weights[0] = tweights[0]; weights[1] = tweights[1]; weights[2] = tweights[2]; weights[3] = 0.0; } return 1; } tri[0] = quad[1]; tri[1] = quad[3]; tri[2] = quad[2]; if (coords_triangle(pt, tri, tweights)) { if (weights) { weights[0] = 0.0; weights[1] = tweights[0]; weights[2] = tweights[2]; weights[3] = tweights[1]; } return 1; } return 0; } static float LineLength(float **pts, int pDim) { int i; float d, l; float *p0 = pts[0]; float *p1 = pts[1]; l = 0.0; for (i = 0; i < pDim; i++) { d = p0[i] - p1[i]; l += d*d; } return sqrt(l); } #define INTERPOLATE(type, round) \ { \ Interp *iPtr; \ char *sData = (char *)DXGetArrayData(sA); \ type *dPtr = (type *)DXGetArrayData(dA); \ \ DXInitGetNextSegListItem(iList); \ while (NULL != (iPtr = (Interp *)DXGetNextSegListItem(iList))) \ { \ int j, k, *elt; \ type *sPtr; \ \ elt = (int *)DXGetArrayEntry(elements, iPtr->index, eBuf); \ \ sPtr = (type *)(sData + elt[0]*itemSize); \ \ for (j = 0; j < itemValues; j++) \ fBuf[j] = iPtr->weights[0]*sPtr[j]; \ \ for (j = 1; j < nV; j++) \ { \ sPtr = (type *)(sData + elt[j]*itemSize); \ \ for (k = 0; k < itemValues; k++) \ fBuf[k] += iPtr->weights[j]*sPtr[k]; \ } \ \ for (j = 0; j < itemValues; j++) \ dPtr[j] = fBuf[j] + round; \ \ dPtr += itemValues; \ } \ } static Field CreateOutputIrreg(Field f, SegList *pList, SegList *iList) { Array cA = NULL, sA, dA = NULL; char *dptr; SegListSegment *slist; char *name, *toDelete[100], nDelete = 0; Object at; int nV; int i, nSamples = DXGetSegListItemCount(pList); float *fBuf = NULL; ArrayHandle elements = NULL; int eBuf[32]; cA = (Array)DXGetComponentValue(f, "connections"); if (! cA) { DXSetError(ERROR_INTERNAL, "missing connections"); goto error; } DXReference((Object)cA); DXDeleteComponent(f, "connections"); elements = DXCreateArrayHandle(cA); if (! elements) goto error; DXGetArrayInfo(cA, NULL, NULL, NULL, NULL, &nV); i = 0; while(NULL != (sA = (Array)DXGetEnumeratedComponentValue(f, i++, &name))) { Object at; Type t; Category c; int r, shape[32]; int size; int itemSize, itemValues; int nItems; if (! strcmp(name, "positions")) continue; if (DXGetComponentAttribute(f, name, "ref") || DXGetComponentAttribute(f, name, "der") || !strncmp(name, "invalid", 8) || NULL == (at = DXGetComponentAttribute(f, name, "dep"))) { toDelete[nDelete++] = name; continue; } if (DXGetObjectClass(at) != CLASS_STRING) { DXSetError(ERROR_BAD_CLASS, "%s dependency attribute", name); goto error; } DXGetArrayInfo(sA, NULL, &t, &c, &r, shape); itemSize = DXGetItemSize(sA); itemValues = itemSize / DXTypeSize(t); if (DXQueryConstantArray(sA, NULL, NULL)) { dA = (Array)DXNewConstantArrayV(DXGetSegListItemCount(pList), DXGetConstantArrayData(sA), t, c, r, shape); if (! dA) goto error; } else { nItems = DXGetItemSize(sA) / DXTypeSize(t); dA = DXNewArrayV(t, c, r, shape); if (! dA) goto error; if (! DXAddArrayData(dA, 0, nSamples, NULL)) goto error; if (! strcmp(DXGetString((String)at), "positions")) { fBuf = (float *)DXAllocate(nItems * sizeof(float)); if (fBuf == NULL) goto error; switch(t) { case TYPE_FLOAT: INTERPOLATE(float, 0.0); break; case TYPE_DOUBLE: INTERPOLATE(double, 0.0); break; case TYPE_SHORT: INTERPOLATE(short, 0.5); break; case TYPE_USHORT: INTERPOLATE(ushort, 0.5); break; case TYPE_INT: INTERPOLATE(int, 0.5); break; case TYPE_UINT: INTERPOLATE(uint, 0.5); break; case TYPE_UBYTE: INTERPOLATE(ubyte, 0.5); break; case TYPE_BYTE: INTERPOLATE(byte, 0.5); break; } DXFree((Pointer)fBuf); fBuf = NULL; } else if (! strcmp(DXGetString((String)at), "connections")) { char *src, *dst; Interp *iPtr; src = (char *)DXGetArrayData(sA); dst = (char *)DXGetArrayData(dA); DXInitGetNextSegListItem(iList); while (NULL != (iPtr = (Interp *)DXGetNextSegListItem(iList))) { memcpy(dst, src+(iPtr->index*itemSize), itemSize); dst += itemSize; } } else { DXSetError(ERROR_BAD_CLASS, "%s dependency attribute", name); goto error; } } if (! DXSetComponentValue(f, name, (Object)dA)) goto error; if (! DXSetComponentAttribute(f, name, "dep", (Object)DXNewString("positions"))) goto error; } for (i = 0; i < nDelete; i++) DXDeleteComponent(f, toDelete[i]); sA = (Array)DXGetComponentValue(f, "positions"); if (! sA) goto error; if (! DXGetArrayInfo(sA, NULL, NULL, NULL, NULL, &nV)) goto error; dA = DXNewArray(TYPE_FLOAT, CATEGORY_REAL, 1, nV); if (! dA) goto error; if (! DXAddArrayData(dA, 0, DXGetSegListItemCount(pList), NULL)) goto error; dptr = (char *)DXGetArrayData(dA); DXInitGetNextSegListSegment(pList); while (NULL != (slist = DXGetNextSegListSegment(pList))) { int size = DXGetSegListSegmentItemCount(slist) * nV * sizeof(float); memcpy(dptr, DXGetSegListSegmentPointer(slist), size); dptr += size; } if (! DXSetComponentValue(f, "positions", (Object)dA)) goto error; DXFreeArrayHandle(elements); DXDelete((Object)cA); return f; error: DXFreeArrayHandle(elements); DXDelete((Object)cA); DXFree((Pointer)fBuf); return NULL; } #define INTERPOLATE_POLYLINE(type, round) \ { \ PolylineInterp *iPtr; \ type *dPtr = (type *)DXGetArrayData(dA); \ \ DXInitGetNextSegListItem(iList); \ while (NULL != (iPtr = (PolylineInterp *)DXGetNextSegListItem(iList))) \ { \ int j; \ type *s0; \ type *s1; \ \ s0 = (type *)DXGetArrayEntry(ah, iPtr->i0, buf0); \ s1 = (type *)DXGetArrayEntry(ah, iPtr->i1, buf1); \ \ for (j = 0; j < itemValues; j++) \ *dPtr++ = (type)(iPtr->weights[0]*s0[j]+iPtr->weights[1]*s1[j]+round); \ } \ } static Field CreateOutputPolylines(Field f, SegList *pList, SegList *iList) { Array sA, dA = NULL; char *dptr; SegListSegment *slist; char *name, *toDelete[100], nDelete = 0; Object at; int nV = 2; int i, nSamples = DXGetSegListItemCount(pList); ArrayHandle ah = NULL; Pointer buf0 = NULL, buf1 = NULL; DXDeleteComponent(f, "polylines"); DXDeleteComponent(f, "edges"); i = 0; while(NULL != (sA = (Array)DXGetEnumeratedComponentValue(f, i++, &name))) { Object at; Type t; Category c; int r, shape[32]; int size; int itemSize, itemValues; int nItems; if (! strcmp(name, "positions")) continue; if (DXGetComponentAttribute(f, name, "ref") || DXGetComponentAttribute(f, name, "der") || !strncmp(name, "invalid", 8) || NULL == (at = DXGetComponentAttribute(f, name, "dep"))) { toDelete[nDelete++] = name; continue; } if (DXGetObjectClass(at) != CLASS_STRING) { DXSetError(ERROR_BAD_CLASS, "%s dependency attribute", name); goto error; } DXGetArrayInfo(sA, NULL, &t, &c, &r, shape); itemSize = DXGetItemSize(sA); itemValues = itemSize / DXTypeSize(t); if (DXQueryConstantArray(sA, NULL, NULL)) { dA = (Array)DXNewConstantArrayV(DXGetSegListItemCount(pList), DXGetConstantArrayData(sA), t, c, r, shape); if (! dA) goto error; } else { nItems = DXGetItemSize(sA) / DXTypeSize(t); dA = DXNewArrayV(t, c, r, shape); if (! dA) goto error; if (! DXAddArrayData(dA, 0, nSamples, NULL)) goto error; ah = DXCreateArrayHandle(sA); if (! ah) goto error; buf0 = DXAllocate(DXGetItemSize(sA)); buf1 = DXAllocate(DXGetItemSize(sA)); if (! buf0 || ! buf1) goto error; if (! strcmp(DXGetString((String)at), "positions")) { switch(t) { case TYPE_FLOAT: INTERPOLATE_POLYLINE(float, 0.0); break; case TYPE_DOUBLE: INTERPOLATE_POLYLINE(double, 0.0); break; case TYPE_SHORT: INTERPOLATE_POLYLINE(short, 0.5); break; case TYPE_USHORT: INTERPOLATE_POLYLINE(ushort, 0.5); break; case TYPE_INT: INTERPOLATE_POLYLINE(int, 0.5); break; case TYPE_UINT: INTERPOLATE_POLYLINE(uint, 0.5); break; case TYPE_UBYTE: INTERPOLATE_POLYLINE(ubyte, 0.5); break; case TYPE_BYTE: INTERPOLATE_POLYLINE(byte, 0.5); break; } } else if (! strcmp(DXGetString((String)at), "polylines")) { char *src, *dst; PolylineInterp *iPtr; dst = (char *)DXGetArrayData(dA); DXInitGetNextSegListItem(iList); while (NULL != (iPtr = (PolylineInterp *)DXGetNextSegListItem(iList))) { memcpy(dst, DXGetArrayEntry(ah, iPtr->index, buf0), itemSize); dst += itemSize; } } else { DXSetError(ERROR_BAD_CLASS, "%s dependency attribute", name); goto error; } } DXFree(buf0); buf0 = NULL; DXFree(buf1); buf1 = NULL; DXFreeArrayHandle(ah); ah = NULL; if (! DXSetComponentValue(f, name, (Object)dA)) goto error; if (! DXSetComponentAttribute(f, name, "dep", (Object)DXNewString("positions"))) goto error; } for (i = 0; i < nDelete; i++) DXDeleteComponent(f, toDelete[i]); sA = (Array)DXGetComponentValue(f, "positions"); if (! sA) goto error; if (! DXGetArrayInfo(sA, NULL, NULL, NULL, NULL, &nV)) goto error; dA = DXNewArray(TYPE_FLOAT, CATEGORY_REAL, 1, nV); if (! dA) goto error; if (! DXAddArrayData(dA, 0, DXGetSegListItemCount(pList), NULL)) goto error; dptr = (char *)DXGetArrayData(dA); DXInitGetNextSegListSegment(pList); while (NULL != (slist = DXGetNextSegListSegment(pList))) { int size = DXGetSegListSegmentItemCount(slist) * nV * sizeof(float); memcpy(dptr, DXGetSegListSegmentPointer(slist), size); dptr += size; } if (! DXSetComponentValue(f, "positions", (Object)dA)) goto error; return f; error: DXFree(buf0); DXFree(buf1); DXFree(ah); return NULL; } #define SQUEEZE(knts, del, pDim, eDim) \ { \ int i, j, k; \ \ /* \ * Squeeze out unnecessary axes \ */ \ for (i = j = 0; i < pDim; i++) \ { \ if (knts[i] != 1) \ { \ if (i != j) \ { \ knts[j] = knts[i]; \ for (k = 0; k < pDim; k++) \ del[j*pDim+k] = del[i*pDim+k]; \ } \ j++; \ } \ } \ if (j != eDim) \ { \ DXSetError(ERROR_INVALID_DATA, "dimension mismatch"); \ goto error; \ } \ for (j = j; j < pDim; j++) \ { \ knts[j] = 1; \ for (k = 0; k < pDim; k++) \ del[j*sd+k] = 0.0; \ } \ } static void _Get_Counts(int ndim, int nSamples, float *lengths, int *samples) { switch(ndim) { float A, B, a; case 1: samples[0] = nSamples; break; case 2: A = lengths[1] / lengths[0]; a = sqrt(nSamples / A); samples[0] = (int)(a + 0.5); samples[1] = (int)((a * A) + 0.5); if (! samples[0]) { samples[0] = nSamples; samples[1] = 1; } else if (! samples[1]) { samples[0] = 1; samples[1] = nSamples; } break; case 3: A = lengths[1] / lengths[0]; B = lengths[2] / lengths[0]; a = pow(nSamples / (A*B), 0.333333333); samples[0] = (int)(a + 0.5); samples[1] = (int)((a * A) + 0.5); samples[2] = (int)((a * B) + 0.5); if (! samples[0] || ! samples[1] || !samples[2]) { int i, j, tmp_samples[3]; float tmp_lengths[3]; for (i = j = 0; i < 3; i++) if (samples[i]) tmp_lengths[j++] = lengths[i]; if (j > 0) _Get_Counts(j, nSamples, tmp_lengths, tmp_samples); for (i = j = 0; i < 3; i++) if (samples[i]) samples[i] = tmp_samples[j++]; else samples[i] = 1; } break; } } static Array _Sample_Grid_Regular(Object o, int nSamples, int eltDim) { Field field; int i, j, k, sd, ed, gridConnections; float deltas[9], *d; int counts[3], samples[3]; float origin[3]; float lengths[3]; Class class; Array grid; class = DXGetObjectClass(o); if (class == CLASS_GROUP) class = DXGetGroupClass((Group)o); ed = eltDim; if (class == CLASS_COMPOSITEFIELD) { Array cA, pA; int maxGrid[3], minGrid[3], pc[3], part; float m[9], offset[3]; /* * We need the overall composite field origin and counts. */ for (part = 0, field = NULL; ;part++) { Field f = (Field)DXGetEnumeratedMember((Group)o, part, NULL); if (! f) break; if (DXEmptyField(f)) continue; if (field == NULL) { int pc[3]; int cc[3]; field = f; cA = (Array)DXGetComponentValue(field, "connections"); if (! cA) goto error; pA = (Array)DXGetComponentValue(field, "positions"); if (! pA) goto error; if (! DXQueryGridPositions(pA, &sd, pc, origin, deltas)) goto error; SQUEEZE(pc, deltas, sd, ed); /* * If the connections are grid, we can get this info from * the mesh offset info. Initialize for a search for the * maximum count, and determine the overall origin. */ gridConnections = (NULL != DXQueryGridConnections(cA, &ed, cc)); if (gridConnections) { int mo[3]; for (i = 0; i < ed; i++) minGrid[i] = 0; if (! DXGetMeshOffsets((MeshArray)cA, mo)) goto error; /* * determine overall offset by subtracting off the * offset to the mesh origin. Also initialize the max * grid search. */ for (i = 0; i < ed; i++) { for (k = 0; k < sd; k++) origin[k] -= mo[i]*deltas[i*sd+k]; maxGrid[i] = cc[i] + mo[i] - 1; } } else { /* * if not, then we use the origin of the first field * as the temporary origin. For each field, we get the * grid coordinates of its origin relative to the * (0,0,0) origin of the first field and use that * as mesh offset. This requires determining a back * transformation from (xyz) space to the (ijk) space * of the first field. * * Determine the inverse deltas matrix. The result * will map points in (xyz) space to (ijk) coordinates * in the first field's grid. */ if (! invert(sd, deltas, m)) goto error; /* * For each member of the composite field, we will * determine the pseudo-mesh offsets relative to the * grid of this field. */ for (i = 0; i < sd; i++) { minGrid[i] = 0; maxGrid[i] = pc[i] - 1; } } } else { if (gridConnections) { int mo[3], cnts[3]; cA = (Array)DXGetComponentValue(f, "connections"); if (! cA) goto error; if (! DXQueryGridConnections(cA, NULL, cnts)) goto error; if (! DXGetMeshOffsets((MeshArray)cA, mo)) goto error; for (i = 0; i < sd; i++) { j = mo[i] + cnts[i] - 1; if (j > maxGrid[i]) maxGrid[i] = j; } } else { int mo[3]; int pc[3]; float po[3]; float pdel[9]; pA = (Array)DXGetComponentValue(f, "positions"); if (! pA) goto error; if (! DXQueryGridPositions(pA, NULL, pc, po, pdel)) goto error; SQUEEZE(pc, pdel, sd, ed); /* * back-transform to get the origin in the first field's * grid coordinates */ for (j = 0; j < sd; j++) po[j] -= origin[j]; for (i = 0; i < sd; i++) { float f; f = 0; for (j = 0; j < sd; j++) f += po[j]*m[j*sd+i]; mo[i] = (int)(f + 0.5); } for (i = 0; i < sd; i++) { if (mo[i] < minGrid[i]) minGrid[i] = mo[i]; j = mo[i] + pc[i] - 1; if (j > maxGrid[i]) maxGrid[i] = j; } } } } if (gridConnections) { for (i = 0; i < ed; i++) counts[i] = maxGrid[i] + 1; } else { Array pA; float po[3]; pA = (Array)DXGetComponentValue(field, "positions"); if (! pA) goto error; if (! DXQueryGridPositions(pA, &sd, pc, po, deltas)) goto error; SQUEEZE(pc, deltas, sd, ed); /* * max and min grids are relative to this grid's origin. Overall * origin is offset from here by minGrid offsets. */ for (i = 0; i < sd; i++) { origin[i] = po[i]; for (j = 0; j < ed; j++) origin[i] += minGrid[j]*deltas[j*sd+i]; } /* * For the meaningful axes, counts are max - min, else 1. */ for (i = 0; i < ed; i++) counts[i] = (maxGrid[i] - minGrid[i]) + 1; } } else if (class == CLASS_FIELD) { Array pA; field = (Field)o; pA = (Array)DXGetComponentValue(field, "positions"); if (! pA) goto error; if (! DXQueryGridPositions(pA, &sd, counts, origin, deltas)) goto error; /* * Squeeze out unnecessary axes */ for (i = j = 0; i < sd; i++) { if (counts[i] != 1) { if (i != j) { counts[j] = counts[i]; for (k = 0; k < sd; k++) deltas[j*sd+k] = deltas[i*sd+k]; } j++; } } for (j = j; j < sd; j++) { counts[j] = 1; for (k = 0; k < sd; k++) deltas[j*sd+k] = 0.0; } } /* * Determine lengths of meaningful deltas */ for (i = 0, d = deltas; i < ed; i++) { lengths[i] = 0.0; for (j = 0; j < sd; j++, d++) lengths[i] += (*d)*(*d); lengths[i] = (counts[i]-1)*sqrt(lengths[i]); } if (nSamples > 1) { /* * Determine allocation of samples along meaningful edges */ _Get_Counts(ed, nSamples, lengths, samples); for (i; i < sd; i++) samples[i] = 1; } else for (i=0; i < sd; i++) samples[i] = 1; /* * Scale deltas by ratio of number of original steps to number of * grid steps. */ for (i = 0; i < ed; i++) { float s = ((float)(counts[i] - 1)) / samples[i]; for (j = 0; j < sd; j++) deltas[i*sd+j] *= s; } /* * Add half a delta to the origin */ for (i = 0; i < ed; i++) { for (j = 0; j < sd; j++) origin[j] += (deltas[i*sd+j] / 2.0); } /* * flesh the deltas out so that the deltas matrix is square */ for (i = ed; i < sd; i++) for (j = 0; j < sd; j++) deltas[i*sd+j] = 0.0; grid = DXMakeGridPositionsV(sd, samples, origin, deltas); return grid; error: return NULL; } static Error invert(int n, float *in, float *out) { switch(n) { case 3: { if (! invert33(in, out)) return ERROR; break; } case 2: { float m; m = in[0]*in[3] - in[1]*in[2]; if (m == 0) return ERROR; m = 1.0 / m; out[0] = in[3]*m; out[1] = -in[1]*m; out[2] = -in[2]*m; out[3] = in[0]*m; break; } case 1: { if (in[0] == 0) { DXSetError(ERROR_INVALID_DATA, "bad delta vectors"); return ERROR; } out[0] = 1.0 / in[0]; break; } } return OK; } static float rand1() { return ((float)(rand()&0x7fff)) / 0x7fff; } #define FIELD_TRIES 1024 static int FieldVolumeTask(Pointer ptr) { Field field = (Field)ptr; if (-1.0 == _dxfFieldVolume(field)) return ERROR; return OK; } static float StandardFieldVolume(Field); static float PolylineFieldVolume(Field); static float FaceFieldVolume(Field); static float _dxfFieldVolume(Field field) { Array pArray, cArray, vArray = NULL; float fvol = 0.0; if (DXEmptyField(field)) goto done; pArray = (Array)DXGetComponentValue(field, "positions"); if (! pArray) { DXSetError(ERROR_MISSING_DATA, "missing positions"); goto error; } cArray = (Array)DXGetComponentValue(field, "connections"); if (cArray) { if (DXQueryGridPositions(pArray, NULL, NULL, NULL, NULL) && DXQueryGridConnections(cArray, NULL, NULL) && ! DXGetComponentValue(field, "invalid connections")) { fvol = _dxfBoxVolume((Object)field, NULL); } else { fvol = StandardFieldVolume(field); } } else if (DXGetComponentValue(field, "polylines")) { fvol = PolylineFieldVolume(field); } else if (DXGetComponentValue(field, "faces")) { fvol = PolylineFieldVolume(field); } else { InvalidComponentHandle invalids; invalids = DXCreateInvalidComponentHandle((Object)field, NULL, "positions"); if (invalids == NULL) goto error; fvol = DXGetValidCount(invalids); DXFreeInvalidComponentHandle(invalids); invalids = NULL; goto done; } done: vArray = (Array)DXNewConstantArray(1, (Pointer)&fvol, TYPE_FLOAT, CATEGORY_REAL, 0); if (! vArray) goto error; if (! DXSetAttribute((Object)field, "volume", (Object)vArray)) goto error; return fvol; error: DXDelete((Object)vArray); return -1; } static float StandardFieldVolume(Field field) { Array pArray, cArray, vArray = NULL; InvalidComponentHandle invalids = NULL; ArrayHandle pHandle = NULL, cHandle = NULL; float fvol = 0.0; int pDim, nCon, nSamples, *element; int i, j, n; Object attr; char *str; static float (*metric)(float **, int); int eltSize; if (DXEmptyField(field)) goto done; pArray = (Array)DXGetComponentValue(field, "positions"); if (! pArray) { DXSetError(ERROR_MISSING_DATA, "missing positions"); goto error; } cArray = (Array)DXGetComponentValue(field, "connections"); DXGetArrayInfo(pArray, NULL, NULL, NULL, NULL, &pDim); DXGetArrayInfo(cArray, &nCon, NULL, NULL, NULL, &eltSize); pHandle = DXCreateArrayHandle(pArray); cHandle = DXCreateArrayHandle(cArray); if (! pHandle || !cHandle) goto error; if (DXGetComponentValue(field, "invalid connections")) { invalids = DXCreateInvalidComponentHandle((Object)field, NULL, "connections"); if (! invalids) goto error; } metric = _dxfGetVolumeMethod(field); if (! metric) goto error; fvol = 0.0; n = 0; for (i = 0; i < nCon; i++) if (!invalids || DXIsElementValid(invalids, i)) { int *elt, eBuf[32]; float *pts[32], pBuf[32][3]; n ++; elt = (int *)DXGetArrayEntry(cHandle, i, (Pointer)eBuf); for (j = 0; j < eltSize; j++) pts[j] = (float *)DXGetArrayEntry(pHandle, elt[j], (Pointer)&pBuf[j][0]); fvol += (*metric)(pts, pDim); } DXFreeArrayHandle(pHandle); pHandle = NULL; DXFreeArrayHandle(cHandle); cHandle = NULL; DXFreeInvalidComponentHandle(invalids); invalids = NULL; done: return fvol; error: if (! pHandle) DXFreeArrayHandle(pHandle); if (! cHandle) DXFreeArrayHandle(cHandle); if (! invalids) DXFreeInvalidComponentHandle(invalids); DXDelete((Object)vArray); return -1; } static float PolylineFieldVolume(Field field) { Array pArray, plArray, eArray; int i, j, pDim; float buf0[32], buf1[32]; float *ptrs[2]; int nPolylines, nEdges, *polylines, *edges; ArrayHandle pHandle; float fvol; pArray = (Array)DXGetComponentValue(field, "positions"); if (! pArray) { DXSetError(ERROR_MISSING_DATA, "missing positions"); goto error; } plArray = (Array)DXGetComponentValue(field, "polylines"); if (! plArray) { DXSetError(ERROR_MISSING_DATA, "missing polylines"); goto error; } eArray = (Array)DXGetComponentValue(field, "edges"); if (! eArray) { DXSetError(ERROR_MISSING_DATA, "polylines but no edges"); goto error; } DXGetArrayInfo(plArray, &nPolylines, NULL, NULL, NULL, NULL); DXGetArrayInfo(eArray, &nEdges, NULL, NULL, NULL, NULL); DXGetArrayInfo(pArray, NULL, NULL, NULL, NULL, &pDim); polylines = (int *)DXGetArrayData(plArray); edges = (int *)DXGetArrayData(eArray); pHandle = DXCreateArrayHandle(pArray); for (i = 0, fvol = 0.0; i < nPolylines; i++) { int start = polylines[i]; int end = (i == nPolylines-1) ? nEdges : polylines[i+1]; int knt = (end - start) - 1; int *seg = edges + start; ptrs[1] = (float *)DXGetArrayEntry(pHandle, *seg, buf1); for (j = 0, seg++; j < knt; j++, seg++) { ptrs[0] = ptrs[1]; ptrs[1] = (float *)DXGetArrayEntry(pHandle, *seg, (j & 0x1) ? buf1 : buf0); fvol += LineLength(ptrs, pDim); } } pHandle = NULL; return fvol; error: return -1; } extern Error _dxfLoopNormal(int *, ArrayHandle, float *, int, int, Vector *); static float FaceFieldVolume(Field field) { Array faces, loops, edges, positions; int *faceData, *loopData, *edgeData; ArrayHandle phandle = NULL; int nFaces, nLoops, nEdges, nDim; float area = 0; Type t; Category c; int rank, shape[32]; int i, face; float *p0 = NULL, *p1 = NULL, *p2 = NULL; InvalidComponentHandle ich = NULL; Array normals = NULL, areas = NULL; Vector *normData, *fn; float *areaData, *fa; if (DXEmptyField(field)) goto done; faces = (Array)DXGetComponentValue(field, "faces"); if (! faces) { DXSetError(ERROR_MISSING_DATA, "faces"); goto error; } DXGetArrayInfo(faces, &nFaces, &t, &c, &rank, shape); if (nFaces == 0) goto done; if (t != TYPE_INT || c != CATEGORY_REAL || !(rank == 0 || (rank == 1 && shape[0] == 1))) { DXSetError(ERROR_INVALID_DATA, "faces must be real integer scalars"); goto error; } loops = (Array)DXGetComponentValue(field, "loops"); if (! loops) { DXSetError(ERROR_MISSING_DATA, "loops"); goto error; } DXGetArrayInfo(loops, &nLoops, &t, &c, &rank, shape); if (nLoops == 0) goto done; if (t != TYPE_INT || c != CATEGORY_REAL || !(rank == 0 || (rank == 1 && shape[0] == 1))) { DXSetError(ERROR_INVALID_DATA, "loops must be real integer scalars"); goto error; } edges = (Array)DXGetComponentValue(field, "edges"); if (! edges) { DXSetError(ERROR_MISSING_DATA, "edges"); goto error; } DXGetArrayInfo(edges, &nEdges, &t, &c, &rank, shape); if (nEdges == 0) goto done; if (t != TYPE_INT || c != CATEGORY_REAL || !(rank == 0 || (rank == 1 && shape[0] == 1))) { DXSetError(ERROR_INVALID_DATA, "edges must be real integer scalars"); goto error; } positions = (Array)DXGetComponentValue(field, "positions"); if (! positions) { DXSetError(ERROR_MISSING_DATA, "positions"); goto error; } DXGetArrayInfo(positions, NULL, &t, &c, &rank, &nDim); if (t != TYPE_FLOAT || c != CATEGORY_REAL || rank != 1 || (nDim != 2 && nDim != 3)) { DXSetError(ERROR_INVALID_DATA, "positions must be real float 2- or 3-vectors"); goto error; } phandle = DXCreateArrayHandle(positions); if (! phandle) goto error; normals = DXNewArray(TYPE_FLOAT, CATEGORY_REAL, 1, 3); if (! normals) goto error; if (! DXAddArrayData(normals, 0, nFaces, NULL)) goto error; areas = DXNewArray(TYPE_FLOAT, CATEGORY_REAL, 0); if (! areas) goto error; if (! DXAddArrayData(areas, 0, nFaces, NULL)) goto error; faceData = (int *)DXGetArrayData(faces); loopData = (int *)DXGetArrayData(loops); edgeData = (int *)DXGetArrayData(edges); normData = (Vector *)DXGetArrayData(normals); areaData = (float *)DXGetArrayData(areas); if (! DXSetComponentValue(field, "face normals", (Object)normals)) goto error; normals = NULL; if (! DXSetComponentValue(field, "face areas", (Object)areas)) goto error; areas = NULL; if (DXGetComponentValue(field, "invalid faces")) { ich = DXCreateInvalidComponentHandle((Object)field, NULL, "faces"); if (ich == NULL) goto error; } fn = normData; fa = areaData; for (face = 0; face < nFaces; face++, fn++, fa++) { int fLoop = faceData[face]; int lLoop = (face == nFaces-1) ? nLoops : faceData[face+1]; int loop; float l; if (ich && DXIsElementInvalid(ich, face)) continue; for (loop = fLoop; loop < lLoop; loop++) { int fEdge = loopData[loop]; int lEdge = (loop == nLoops-1) ? nEdges : loopData[loop+1]; if (loop == fLoop) { if (! _dxfLoopNormal(edgeData+fEdge, phandle, NULL, (lEdge-fEdge), nDim, fn)) goto error; l = *fa = vector_length(fn); } else { Vector ln; if (! _dxfLoopNormal(edgeData+fEdge, phandle, NULL, (lEdge-fEdge), nDim, &ln)) goto error; if ((fn->x*ln.x + fn->y*ln.y + fn->z*ln.z) > 0) *fa += vector_length(&ln); else *fa -= vector_length(&ln); } } *fa = fabs(*fa / 2); area += *fa; if (l > 0) vector_scale((Vector *)fn, 1.0/l, (Vector *)fn); } done: if (ich) DXFreeInvalidComponentHandle(ich); DXFreeArrayHandle(phandle); return area; error: if (ich) DXFreeInvalidComponentHandle(ich); DXFreeArrayHandle(phandle); DXFree((Object)normals); DXFree((Object)areas); return -1; } static Error Cleanup_Sample_Field(Field field) { DXDeleteComponent(field, "connections"); DXSetAttribute((Object)field, "volume", NULL); return OK; } static float _dxfBoxVolume(Object object, int *nd) { Point boxPts[8], *b; int i; Point min, max; float ext, vol; if (! DXBoundingBox(object, boxPts)) { DXSetError(ERROR_INVALID_DATA, "unable to compute bounding box"); goto error; } min.x = min.y = min.z = DXD_MAX_FLOAT; max.x = max.y = max.z = -DXD_MAX_FLOAT; for (b = boxPts, i = 0; i < 8; i++, b++) { if (min.x > b->x) min.x = b->x; if (max.x < b->x) max.x = b->x; if (min.y > b->y) min.y = b->y; if (max.y < b->y) max.y = b->y; if (min.z > b->z) min.z = b->z; if (max.z < b->z) max.z = b->z; } vol = 1.0; i = 0; ext = max.x - min.x; if (ext > 0) { vol *= ext; i++; } ext = max.y - min.y; if (ext > 0) { vol *= ext; i++; } ext = max.z - min.z; if (ext > 0) { vol *= ext; i++; } if (nd) *nd = i; return vol; error: return -1.0; } typedef struct { float x, X, y, Y; } MinMax; static int InsideFace(float x, float y, int face, float *positions, int *faces, int nF, int *loops, int nL, int *edges, int nE, MinMax *loopMM) { int l; int ls = faces[face]; int le = (face == nF-1) ? nL : faces[face+1]; MinMax *lmm; int knt = 0; for (l = ls, lmm = loopMM + ls; l < le; l++, lmm++) { int e; int es = loops[l]; int ee = (l == nL-1) ? nE : loops[l+1]; if (x >= lmm->x && x <= lmm->X && y >= lmm->y && y <= lmm->Y) { int p, q = edges[es]; float *ppos, *qpos = positions+(q<<1); float *t, *b, intercept, d; for (e = es; e < ee; e++) { p = q; ppos = qpos; q = (e == ee-1) ? edges[es] : edges[e+1]; qpos = positions+(q<<1); if (qpos[1] > ppos[1]) t = qpos, b = ppos; else if (qpos[1] < ppos[1]) t = ppos, b = qpos; else continue; if (t[1] < y) continue; if (b[1] >= y) continue; if (t[0] < x && b[0] < x) continue; if (t[0] > x && b[0] > x) { knt++; continue; } d = (y - t[1]) / (b[1] - t[1]); intercept = t[0] + d*(b[0] - t[0]); if (intercept > x) knt++; } } } return knt & 0x1; } static Error FacesSample(Field field, int n, Array grid) { Array pArray, fArray, lArray, eArray, nArray, aArray; int nP, nF, nL, nE, nD, planar; float *positions2D = NULL, *positions3D; float *positions, *fareas; int *faces, *loops, *edges; Vector *fnorms; int i, f, l, e, nDim; ubyte *flags = NULL; float A, B, C, D; SegList *pList = NULL, *iList = NULL; int axis; InvalidComponentHandle ich = NULL; float totArea; MinMax *loopMM = NULL; MinMax *faceMM = NULL; MinMax minmax; if (DXEmptyField(field)) goto empty; if (! DXGetComponentValue(field, "face normals")) if ((totArea = FaceFieldVolume(field)) < 0.0) goto error; pArray = (Array)DXGetComponentValue(field, "positions"); fArray = (Array)DXGetComponentValue(field, "faces"); lArray = (Array)DXGetComponentValue(field, "loops"); eArray = (Array)DXGetComponentValue(field, "edges"); nArray = (Array)DXGetComponentValue(field, "face normals"); aArray = (Array)DXGetComponentValue(field, "face areas"); if (!pArray || !fArray || !lArray || !eArray || !nArray || !aArray) goto error; DXGetArrayInfo(pArray, &nP, NULL, NULL, NULL, &nDim); DXGetArrayInfo(fArray, &nF, NULL, NULL, NULL, NULL); DXGetArrayInfo(lArray, &nL, NULL, NULL, NULL, NULL); DXGetArrayInfo(eArray, &nE, NULL, NULL, NULL, NULL); positions = (float *)DXGetArrayData(pArray); faces = (int *)DXGetArrayData(fArray); loops = (int *)DXGetArrayData(lArray); edges = (int *)DXGetArrayData(eArray); fnorms = (Vector *)DXGetArrayData(nArray); fareas = (float *)DXGetArrayData(aArray); loopMM = (MinMax *)DXAllocate(nL*sizeof(MinMax)); faceMM = (MinMax *)DXAllocate(nF*sizeof(MinMax)); if (! loopMM || ! faceMM) goto error; if (DXGetComponentValue(field, "invalid faces")) { ich = DXCreateInvalidComponentHandle((Object)field, NULL, "faces"); if (! ich) goto error; } /* * If the dimensionality is three, we look to see if the * faces all lie in the sample plane, since if using a * single grid of prospective samples will result in a nice * even distribution of samples. */ if (nDim == 3) { Vector *f0 = NULL; Vector *fn = fnorms; int firstValid; /* * They are all planar if the face normals are very close * to collinear. */ for (f = 0, planar = 1; f < nF && planar; f++, fn++) { if (ich && DXIsElementInvalid(ich, f)) continue; if (f0) { if ((f0->x*fn->x + f0->y*fn->y + f0->z*fn->z) < 0.95) planar = 0; } else { f0 = fn; firstValid = f; } } if (! f0) goto empty; /* * If they are coplanar, then project all the positions into the * most perpendicular 2-D plane. */ if (planar) { float *p2d, *p3d; positions2D = (float *)DXAllocate(2 * nP * sizeof(float)); if (! positions2D) goto error; p2d = positions2D; p3d = (float *)positions; if (fabs(f0->x) > fabs(f0->y) && fabs(f0->x) > fabs(f0->z)) { axis = 0; for (i = 0; i < nP; i++) { ++p3d; *p2d++ = *p3d++; *p2d++ = *p3d++; } } else if (fabs(f0->y) > fabs(f0->x) && fabs(f0->y) > fabs(f0->z)) { axis = 1; for (i = 0; i < nP; i++) { *p2d++ = *p3d++; ++p3d; *p2d++ = *p3d++; } } else { axis = 2; for (i = 0; i < nP; i++) { *p2d++ = *p3d++; *p2d++ = *p3d++; ++p3d; } } positions3D = positions; positions = positions2D; /* * Determine the planar equation of the 3-D plane * using the normal to the first valid plane and * one point on that plane. */ p3d = positions3D + 3*edges[loops[faces[firstValid]]]; f0 = fnorms + firstValid; A = f0->x; B = f0->y; C = f0->z; D = -(A*p3d[0] + B*p3d[1] + C*p3d[2]); } } else if (nDim == 2) { planar = 1; } else { DXSetError(ERROR_NOT_IMPLEMENTED, "Sample on faces/loops/edges data that is %d-D", nDim); goto error; } pList = DXNewSegList(nDim*sizeof(float), n, n); iList = DXNewSegList(sizeof(int), n, n); if (! pList || ! iList) goto error; /* * If the faces are coplanar, we use a single grid of prospective * sample points for all the faces */ if (planar) { float xSize, ySize, dx, dy, ox, oy, x, y; int nx, ny, ix, iy, xstart, ystart, xend, yend; MinMax *fmm, *lmm; minmax.x = DXD_MAX_FLOAT, minmax.X = -DXD_MAX_FLOAT; minmax.y = DXD_MAX_FLOAT, minmax.Y = -DXD_MAX_FLOAT; /* * Determine the minmax of the set of faces, each face, * and each loop. */ for (f = 0, fmm = faceMM; f < nF; f++, fmm++) { int lstart = faces[f]; int lend = (f == nF-1) ? nL : faces[f+1]; fmm->x = DXD_MAX_FLOAT, fmm->X = -DXD_MAX_FLOAT; fmm->y = DXD_MAX_FLOAT, fmm->Y = -DXD_MAX_FLOAT; if (ich && DXIsElementInvalid(ich, f)) continue; for (l = lstart, lmm = loopMM + lstart; l < lend; l++, lmm++) { int estart = loops[l]; int eend = (l == nL-1) ? nE: loops[l+1]; lmm->x = DXD_MAX_FLOAT, lmm->X = -DXD_MAX_FLOAT; lmm->y = DXD_MAX_FLOAT, lmm->Y = -DXD_MAX_FLOAT; for (e = estart; e < eend; e++) { float *p = positions + 2*edges[e]; if (*p < lmm->x) lmm->x = *p; if (*p > lmm->X) lmm->X = *p; p++; if (*p < lmm->y) lmm->y = *p; if (*p > lmm->Y) lmm->Y = *p; } if (lmm->x < fmm->x) fmm->x = lmm->x; if (lmm->X > fmm->X) fmm->X = lmm->X; if (lmm->y < fmm->y) fmm->y = lmm->y; if (lmm->Y > fmm->Y) fmm->Y = lmm->Y; } if (fmm->x < minmax.x) minmax.x = fmm->x; if (fmm->X > minmax.X) minmax.X = fmm->X; if (fmm->y < minmax.y) minmax.y = fmm->y; if (fmm->Y > minmax.Y) minmax.Y = fmm->Y; } /* * Determine the x and y extents of the bounding box */ xSize = minmax.X - minmax.x; ySize = minmax.Y - minmax.y; /* * scale the number of samples by the ratio of the area of * the bounding box to that covered by the faces */ n *= ((xSize * ySize) / totArea); /* * allocate samples in rows and columns based on the aspect * ratio of the bounding box. */ nx = sqrt(n*xSize/ySize); ny = n / nx; /* * determine the grid origin and steps */ if (nx <= 0) dx = xSize, nx = 1; else dx = xSize / nx; if (ny <= 0) dy = ySize, ny = 1; else dy = ySize / ny; ox = minmax.x + dx/2; oy = minmax.y + dx/2; /* * for each face, determine which samples it covers */ for (f = 0, fmm = faceMM; f < nF; f++, fmm++) { float fx, fy; int lstart = faces[f]; int lend = (f == nF-1) ? nL : faces[f+1]; if (ich && DXIsElementInvalid(ich, f)) continue; /* * Get the range of prospective sample points that * are covered by the bounding box of the face */ xstart = (fmm->x - ox) / dx; if (ox+xstart*dx < fmm->x) xstart ++; xend = (fmm->X - ox) / dx; ystart = (fmm->y - oy) / dy; if (oy+ystart*dy < fmm->y) ystart ++; yend = (fmm->Y - oy) / dy; if (xstart < 0) xstart = 0; if (ystart < 0) ystart = 0; if (xend >= nx) xend = nx-1; if (yend >= ny) yend = ny-1; fx = ox + xstart*dx; fy = oy + ystart*dy; /* * For each prospective point covered by the bounding box * of the face, see if its inside the face. */ for (x = fx, ix = xstart; ix <= xend; ix++, x += dx) for (y = fy, iy = ystart; iy <= yend; iy++, y += dy) if (InsideFace(x, y, f, positions, faces, nF, loops, nL, edges, nE, loopMM)) { float *pptr = (float *)DXNewSegListItem(pList); Interp *iptr = (Interp *)DXNewSegListItem(iList); if (nDim == 2) { *pptr++ = x; *pptr++ = y; } else { float x3d, y3d, z3d; switch(axis) { case 0: y3d = x; z3d = y; x3d = -(B*y3d + C*z3d + D) / A; break; case 1: x3d = x; z3d = y; y3d = -(A*x3d + C*z3d + D) / B; break; case 2: x3d = x; y3d = y; z3d = -(A*x3d + B*y3d + D) / C; break; } *pptr++ = x3d; *pptr++ = y3d; *pptr++ = z3d; } iptr->index = f; } } } else { /* * Non-planar case. Determine a grid for each face * independently. */ float xSize, ySize, dx, dy, ox, oy, x, y; int nx, ny, ix, iy; MinMax *fmm, *lmm; positions3D = positions; positions2D = (float *)DXAllocate(2 * nP * sizeof(float)); if (! positions2D) goto error; for (f = 0, fmm = faceMM; f < nF; f++, fmm++) { int lstart = faces[f]; int lend = (f == nF-1) ? nL : faces[f+1]; int nface; Vector *fn = fnorms + f; float *p0; if (ich && DXIsElementInvalid(ich, f)) continue; /* * Determine the best axis of projection for the face */ if (fabs(fn->x) > fabs(fn->y) && fabs(fn->x) > fabs(fn->z)) axis = 0; else if (fabs(fn->y) > fabs(fn->x) && fabs(fn->y) > fabs(fn->z)) axis = 1; else axis = 2; /* * project the points used by this face into 2-D. As * this is being done, acquire the loop and face 2-D * minmax. */ fmm->x = DXD_MAX_FLOAT, fmm->X = -DXD_MAX_FLOAT; fmm->y = DXD_MAX_FLOAT, fmm->Y = -DXD_MAX_FLOAT; for (l = lstart, lmm = loopMM + lstart; l < lend; l++, lmm++) { int estart = loops[l]; int eend = (l == nL-1) ? nE : loops[l+1]; int e; float *p2d, *p3d; lmm->x = DXD_MAX_FLOAT, lmm->X = -DXD_MAX_FLOAT; lmm->y = DXD_MAX_FLOAT, lmm->Y = -DXD_MAX_FLOAT; for (e = estart; e < eend; e++) { p2d = positions2D + 2*edges[e]; p3d = positions3D + 3*edges[e]; switch (axis) { case 0: p2d[0] = p3d[1]; p2d[1] = p3d[2]; break; case 1: p2d[0] = p3d[0]; p2d[1] = p3d[2]; break; case 2: p2d[0] = p3d[0]; p2d[1] = p3d[1]; break; } if (p2d[0] < lmm->x) lmm->x = p2d[0]; if (p2d[0] > lmm->X) lmm->X = p2d[0]; if (p3d[1] < lmm->y) lmm->y = p3d[1]; if (p3d[1] > lmm->Y) lmm->Y = p3d[1]; } if (lmm->x < fmm->x) fmm->x = lmm->x; if (lmm->X > fmm->X) fmm->X = lmm->X; if (lmm->y < fmm->y) fmm->y = lmm->y; if (lmm->Y > fmm->Y) fmm->Y = lmm->Y; } /* * Resolve the planar equation based on the face normal * and one vertex of the face. */ p0 = positions3D + 3*edges[loops[faces[f]]]; A = fn->x; B = fn->y; C = fn->z; D = -(A*p0[0] + B*p0[1] + C*p0[2]); /* * Determine the proportion of the samples that * should lie in this face. */ nface = n * (fareas[f] / totArea); if (nface <= 0) continue; /* * determine the 2-D extents of the bounding box of the * 2-D projected face. */ xSize = fmm->X - fmm->x; ySize = fmm->Y - fmm->y; /* * scale the sample count by the ratio of the total area * of the bounding box to the area of the face. */ nface *= ((xSize * ySize) / fareas[f]); if (nface == 0) continue; /* * allocate the prospective samples in rows and columns. */ nx = sqrt(nface*xSize/ySize); ny = nface / nx; /* * get the origin and deltas of the prospective sample grid */ if (nx <= 0) dx = xSize, nx = 1; else dx = xSize / nx; if (ny <= 0) dy = ySize, ny = 1; else dy = ySize / ny; ox = fmm->x + dx/2; oy = fmm->y + dx/2; for (x = ox, ix = 0; ix < nx; ix++, x += dx) for (y = oy, iy = 0; iy < ny; iy++, y += dy) if (InsideFace(x, y, f, positions2D, faces, nF, loops, nL, edges, nE, loopMM)) { float *pptr = (float *)DXNewSegListItem(pList); Interp *iptr = (Interp *)DXNewSegListItem(iList); if (nDim == 2) { *pptr++ = x; *pptr++ = y; } else { float x3d, y3d, z3d; switch(axis) { case 0: y3d = x; z3d = y; x3d = -(B*y3d + C*z3d + D) / A; break; case 1: x3d = x; z3d = y; y3d = -(A*x3d + C*z3d + D) / B; break; case 2: x3d = x; y3d = y; z3d = -(A*x3d + B*y3d + D) / C; break; } *pptr++ = x3d; *pptr++ = y3d; *pptr++ = z3d; } iptr->index = f; } } } if (pList && DXGetSegListItemCount(pList) != 0) { if (! CreateOutputFaces(field, pList, iList)) goto error; goto done; } empty: { char *name; while (DXGetEnumeratedComponentValue(field, 0, &name)) DXDeleteComponent(field, name); } done: if (ich) DXFreeInvalidComponentHandle(ich); if (iList) DXDeleteSegList(iList); if (pList) DXDeleteSegList(pList); if (loopMM) DXFree((Pointer)loopMM); if (faceMM) DXFree((Pointer)faceMM); if (positions2D) DXFree((Pointer)positions2D); return OK; error: if (ich) DXFreeInvalidComponentHandle(ich); if (iList) DXDeleteSegList(iList); if (pList) DXDeleteSegList(pList); if (loopMM) DXFree((Pointer)loopMM); if (faceMM) DXFree((Pointer)faceMM); if (positions2D) DXFree((Pointer)positions2D); return ERROR; } static Field CreateOutputFaces(Field f, SegList *pList, SegList *iList) { Array sA, dA = NULL; char *dptr; SegListSegment *slist; char *name, *toDelete[100], nDelete = 0; Object at; int nV = 2; int i, nSamples = DXGetSegListItemCount(pList); ArrayHandle ah = NULL; Pointer buf0 = NULL, buf1 = NULL; DXDeleteComponent(f, "faces"); DXDeleteComponent(f, "loops"); DXDeleteComponent(f, "edges"); DXDeleteComponent(f, "face normals"); DXDeleteComponent(f, "face areas"); i = 0; while(NULL != (sA = (Array)DXGetEnumeratedComponentValue(f, i++, &name))) { Object at; Type t; Category c; int r, shape[32]; int size; int itemSize, itemValues; int nItems; if (! strcmp(name, "positions")) continue; if (DXGetComponentAttribute(f, name, "ref") || DXGetComponentAttribute(f, name, "der") || !strncmp(name, "invalid", 8) || NULL == (at = DXGetComponentAttribute(f, name, "dep"))) { toDelete[nDelete++] = name; continue; } if (DXGetObjectClass(at) != CLASS_STRING) { DXSetError(ERROR_BAD_CLASS, "%s dependency attribute", name); goto error; } DXGetArrayInfo(sA, NULL, &t, &c, &r, shape); itemSize = DXGetItemSize(sA); itemValues = itemSize / DXTypeSize(t); if (DXQueryConstantArray(sA, NULL, NULL)) { dA = (Array)DXNewConstantArrayV(DXGetSegListItemCount(pList), DXGetConstantArrayData(sA), t, c, r, shape); if (! dA) goto error; } else { nItems = DXGetItemSize(sA) / DXTypeSize(t); dA = DXNewArrayV(t, c, r, shape); if (! dA) goto error; if (! DXAddArrayData(dA, 0, nSamples, NULL)) goto error; ah = DXCreateArrayHandle(sA); if (! ah) goto error; buf0 = DXAllocate(DXGetItemSize(sA)); buf1 = DXAllocate(DXGetItemSize(sA)); if (! buf0 || ! buf1) goto error; if (! strcmp(DXGetString((String)at), "faces")) { char *src, *dst; PolylineInterp *iPtr; dst = (char *)DXGetArrayData(dA); DXInitGetNextSegListItem(iList); while (NULL != (iPtr = (PolylineInterp *)DXGetNextSegListItem(iList))) { memcpy(dst, DXGetArrayEntry(ah, iPtr->index, buf0), itemSize); dst += itemSize; } } else { DXSetError(ERROR_BAD_CLASS, "%s dependency attribute", name); goto error; } } DXFree(buf0); buf0 = NULL; DXFree(buf1); buf1 = NULL; DXFreeArrayHandle(ah); ah = NULL; if (! DXSetComponentValue(f, name, (Object)dA)) goto error; if (! DXSetComponentAttribute(f, name, "dep", (Object)DXNewString("positions"))) goto error; } for (i = 0; i < nDelete; i++) DXDeleteComponent(f, toDelete[i]); sA = (Array)DXGetComponentValue(f, "positions"); if (! sA) goto error; if (! DXGetArrayInfo(sA, NULL, NULL, NULL, NULL, &nV)) goto error; dA = DXNewArray(TYPE_FLOAT, CATEGORY_REAL, 1, nV); if (! dA) goto error; if (! DXAddArrayData(dA, 0, DXGetSegListItemCount(pList), NULL)) goto error; dptr = (char *)DXGetArrayData(dA); DXInitGetNextSegListSegment(pList); while (NULL != (slist = DXGetNextSegListSegment(pList))) { int size = DXGetSegListSegmentItemCount(slist) * nV * sizeof(float); memcpy(dptr, DXGetSegListSegmentPointer(slist), size); dptr += size; } if (! DXSetComponentValue(f, "positions", (Object)dA)) goto error; return f; error: DXFree(buf0); DXFree(buf1); DXFree(ah); return NULL; } /* for a field or multigrid/composite field, if there are no * connections then Sample returns a subset of the existing * points. generate a warning if the number of samples requested * exceeds the number of positions in the input. */ static void _Check_Sample_Count(Object o, int n) { Field field; Object con; int vCount = 0; int i; Class objclass; InvalidComponentHandle ich = NULL; objclass = DXGetObjectClass(o); if (objclass == CLASS_GROUP) objclass = DXGetGroupClass((Group)o); switch (objclass) { case CLASS_FIELD: field = (Field)o; if (DXEmptyField(field)) return; /* if there are connections, with interpolation we can * generate as many points as they want. */ con = DXGetComponentValue(field, "connections"); if (! con) con = DXGetComponentValue(field, "faces"); if (! con) con = DXGetComponentValue(field, "polylines"); if (con) return; ich = DXCreateInvalidComponentHandle((Object)field, NULL, "positions"); if (! ich) return; vCount = DXGetValidCount(ich); if (n > vCount) DXWarning("cannot generate %d sample points; input contains %d valid positions", n, vCount); return; case CLASS_MULTIGRID: case CLASS_COMPOSITEFIELD: for (i=0; field = (Field)DXGetEnumeratedMember((Group)o, i, NULL); i++) { if (DXEmptyField(field)) continue; /* if there are connections, with interpolation we can * generate as many points as they want. */ con = DXGetComponentValue(field, "connections"); if (! con) con = DXGetComponentValue(field, "faces"); if (! con) con = DXGetComponentValue(field, "polylines"); if (con) return; ich = DXCreateInvalidComponentHandle((Object)field, NULL, "positions"); if (! ich) return; vCount += DXGetValidCount(ich); } if (n > vCount) DXWarning("cannot generate %d sample points; input contains %d valid positions", n, vCount); return; } }