/***********************************************************************/ /* 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 /* * $Header: /home/gda/dxcvs/dx/src/exec/dxmods/reduce.c,v 1.3 1999/05/10 15:45:29 gda Exp $: */ /*** SECTION: Rendering MODULE: Reduce SHORTDESCRIPTION: creates a set of reduced-resolution versions of an input field CATEGORY: Import and Export INPUTS: input; field; none; field to reduce factor; scalar list or vector list; all powers of two; set of reduction factors OUTPUTS: reduced; field or group; NULL; set of reduced-resolution versions FLAGS: BUGS: AUTHOR: Greg Abram END: ***/ #include #include #include #define MAX_DIMENSION 10 #define MAX_ELEMENTS 16 #define MAX_STACKBUF 15 /* 5x filter, 3-vector */ #define MAX_REDUCTION 32 #define TRUE 1 #define FALSE 0 #define NULL_DEPENDENCY 0 #define POSITIONS_DEPENDENT 1 #define CONNECTIONS_DEPENDENT 2 #define INVALID_POSITIONS 1 #define INVALID_CONNECTIONS 2 typedef struct oMapElt { int inOffset; int outOffset; } OMapElt; typedef struct originMap { int nDim; OMapElt **map; } *OriginMap; typedef struct { int left, right; } Shoulder; Error m_Reduce(Object *, Object *); static Object Reduce(Object, int, float *, int); static Object ReduceObject(Object, Private, int); static Object ReduceGroup(Group, Private, int); static Object ReduceCField(CompositeField, Private, int); static Object ReduceSeries(Series, Private, int); static Object ReduceField(Field, Private, int, Private); static Error ReduceFieldTask(Pointer); static Array ReduceArrayReg(Array, int, Shoulder *, int *); static Array ReduceArrayIrreg(Array, int, int *, Shoulder *, int *, int *, int *, int); static Error ReduceSumAxis(Pointer, Pointer, int, int, int, int *, int *, int *, int *, int, float *, Type); static Error InvalidateData(Field, Field, int *, int); static void AxesOrder(int *, int *, int); static Object GrowObject(Object, float); static Private ReducedMeshOffsets(CompositeField, Private); static Error DeleteOriginMapObject(Pointer); static Error DeleteFactorsObject(Pointer); #define GetEM(g,i,n) DXGetEnumeratedMember((g),(i),(n)) #define GetECV(g,i,n) DXGetEnumeratedComponentValue((g),(i),(n)) #define GetECA(g,i,n,a) DXGetEnumeratedComponentAttribute((g),(i),(n),(a)) Error m_Reduce(Object *in, Object *out) { Object data; int nFactors; float *factors, maxFactor; int i; int nDim; Type type; Category cat; int rank, shape[32]; if ((data = in[0]) == NULL) DXErrorReturn(ERROR_MISSING_DATA, "data argument missing"); if (in[1] == NULL) { nFactors = 1; nDim = 1; factors = NULL; } else { if (DXGetObjectClass(in[1]) != CLASS_ARRAY) DXErrorReturn(ERROR_BAD_PARAMETER, "bad reduction specification"); if (DXGetType(in[1], &type, &cat, &rank, shape)) { if ((type != TYPE_INT && type != TYPE_FLOAT) || cat != CATEGORY_REAL || (rank != 0 && rank != 1)) DXErrorReturn(ERROR_BAD_PARAMETER, "illegal reduction specification"); if (rank == 1) nDim = shape[0]; else nDim = 1; } else DXErrorReturn(ERROR_BAD_PARAMETER, "bad reduction specification"); if (nDim > MAX_DIMENSION) DXErrorReturn(ERROR_INVALID_DATA, "dimension too high"); if (! DXGetArrayInfo((Array)in[1], &nFactors, NULL, NULL, NULL, NULL)) DXErrorReturn(ERROR_BAD_PARAMETER, "illegal reduction specification"); if ((factors = (float *)DXAllocate(nDim*nFactors*sizeof(float))) == NULL) DXErrorReturn(ERROR_NO_MEMORY, "error allocating factor array"); if (!DXExtractParameter(in[1], TYPE_FLOAT, nDim, nFactors, (Pointer)factors)) { DXFree((Pointer)factors); DXErrorReturn(ERROR_INTERNAL, "error extracting factor parameter"); } maxFactor = 1.0; for (i = 0; i < nDim*nFactors; i++) if (factors[i] < 1.0) { DXSetError(ERROR_BAD_PARAMETER, "reduction factors must be >= 1.0"); DXFree((Pointer)factors); out[0] = NULL; return ERROR; } else if (factors[i] > maxFactor) maxFactor = factors[i]; } out[0] = Reduce(data, nFactors, factors, nDim); DXFree((Pointer)factors); if (! out[0]) return ERROR; else return OK; } static Object Reduce(Object input, int nFactors, float *factorList, int nDim) { int i, j; float *flist[MAX_DIMENSION], *factors, maxFactor; Object template, result, reduced; Private plist[MAX_DIMENSION]; int reduceFlags[MAX_DIMENSION]; template = NULL; reduced = NULL; result = NULL; factors = NULL; for (i = 0; i < nFactors; i++) { flist[i] = (float *)DXAllocate(MAX_DIMENSION * sizeof(float)); if (! flist[i]) { for (--i; i >= 0; i--) DXDelete((Object)plist[i]); return NULL; } plist[i] = DXNewPrivate((Pointer)flist[i], DeleteFactorsObject); if (! plist[i]) { for (--i; i >= 0; i--) DXDelete((Object)plist[i]); return NULL; } DXReference((Object)plist[i]); } if (factorList) { maxFactor = 1; if (nDim == 1) for (i = 0; i < nFactors; i++) { factors = flist[i]; reduceFlags[i] = FALSE; if (factorList[i] <= 0) factorList[i] = 1.0; for (j = 0; j < MAX_DIMENSION; j++) factors[j] = factorList[i]; if (maxFactor < factorList[i]) maxFactor = factorList[i]; if (factorList[i] > 1.0) reduceFlags[i] = TRUE; } else for (i = 0; i < nFactors; i++) { factors = flist[i]; reduceFlags[i] = FALSE; for (j = 0; j < nDim; j++) { if (factorList[i*nDim+j] <= 0) factorList[i*nDim+j] = 1.0; if ((factors[j] = factorList[i*nDim+j]) > 1.0) reduceFlags[i] = TRUE; if (maxFactor < factors[j]) maxFactor = factors[j]; } } } else { for (i = 0; i < nFactors; i++) { factors = flist[i]; for (j = 0; j < MAX_DIMENSION; j++) factors[j] = 2; /* * always reduce... */ reduceFlags[i] = TRUE; } maxFactor = 1 << nFactors; } if (nFactors == 1 && maxFactor == 1.0) return input; template = DXCopy(input, COPY_STRUCTURE); if (! template) goto error; if (!GrowObject(template, maxFactor)) goto error; if (!template || DXGetError()) goto error; DXCreateTaskGroup(); if (nFactors == 1) { result = DXCopy(template, COPY_STRUCTURE); if (! result) { DXAbortTaskGroup(); goto error; } if (!ReduceObject(result, plist[0], nDim)) { DXAbortTaskGroup(); goto error; } } else { result = (Object)DXNewGroup(); if (! result) goto error; for (i = 0; i < nFactors; i++) { if (reduceFlags[i]) { reduced = DXCopy(template, COPY_STRUCTURE); if (! reduced) { DXAbortTaskGroup(); goto error; } if (!ReduceObject(reduced, plist[i], nDim)) { DXAbortTaskGroup(); goto error; } } else reduced = input; if (! DXSetMember((Group)result, NULL, (Object)reduced)) { DXDelete((Object)result); DXDelete((Object)template); return NULL; } reduced = NULL; } if (i == 0) goto error; } if (! DXExecuteTaskGroup() || DXGetError()) goto error; if (template != input) DXDelete((Object)template); for (i = 0; i < nFactors; i++) DXDelete((Object)plist[i]); return result; error: for (i = 0; i < nFactors; i++) DXDelete((Object)plist[i]); if (template != input) DXDelete(template); DXDelete(result); DXDelete(reduced); return NULL; } static Object GrowObject(Object object, float maxFactor) { Class class; Object child; int i, j; char *name, *components[64]; Field field, field0; int overlap; Array array; float pos; overlap = (((int)ceil(maxFactor)) | 0x1) >> 1; class = DXGetObjectClass(object); if (class == CLASS_GROUP) class = DXGetGroupClass((Group)object); switch(class) { case CLASS_XFORM: if (! DXGetXformInfo((Xform)object, &child, 0)) return NULL; if (! GrowObject(child, maxFactor)) return NULL; break; case CLASS_SERIES: i = 0; while (NULL != (child=DXGetSeriesMember((Series)object, i, &pos))) { child = GrowObject(child, maxFactor); if (! child || DXGetError()) return NULL; if (! DXSetSeriesMember((Series)object, i, pos, child)) return NULL; i++; } break; case CLASS_MULTIGRID: case CLASS_GROUP: i = 0; while (NULL != (child = GetEM((Group)object, i, &name))) { child = GrowObject(child, maxFactor); if (! child || DXGetError()) return NULL; if (! DXSetEnumeratedMember((Group)object, i, child)) return NULL; i++; } break; case CLASS_COMPOSITEFIELD: field0 = NULL; for (i = 0; NULL != (field = DXGetPart(object, i)); i++) { if (DXGetObjectClass((Object)field) != CLASS_FIELD) { DXSetError(ERROR_INVALID_DATA, "invalid member of composite field"); return NULL; } if (DXEmptyField(field) || !DXGetComponentValue(field, "connections")) continue; array = (Array)DXGetComponentValue(field, "connections"); if (! array) { DXSetError(ERROR_MISSING_DATA, "missing connections"); return NULL; } if (! DXQueryGridConnections(array, NULL, NULL)) { DXSetError(ERROR_INVALID_DATA, "irregular connections"); return NULL; } field0 = field; } if (field0 == NULL) return object; i = j = 0; while (NULL != GetECV(field0, i++, &components[j])) if (strcmp(components[j], "statistics") && strcmp(components[j], "neighbors") && strcmp(components[j], "positions") && strcmp(components[j], "connections") && strcmp(components[j], "box")) j++; components[j] = NULL; if (! DXGrowV(object, overlap, NULL, components)) return NULL; break; case CLASS_FIELD: if (DXEmptyField((Field)object) || !DXGetComponentValue((Field)object, "connections")) break; i = j = 0; while (NULL != GetECV((Field)object, i++, &components[j])) if (strcmp(components[j], "statistics") && strcmp(components[j], "neighbors") && strcmp(components[j], "positions") && strcmp(components[j], "connections") && strcmp(components[j], "box")) j++; components[j] = NULL; if (! DXGrowV(object, overlap, NULL, components)) return NULL; break; default: DXSetError(ERROR_INVALID_DATA, "unsupported object class encountered"); return NULL; } return object; } static Object ReduceObject(Object object, Private factors, int nDim) { float *f; Class class; Object result; Object child; int i; f = (float *)DXGetPrivateData(factors); for (i = 0; i < nDim; i++) if (f[i] > 1.0) goto reduce; return object; reduce: class = DXGetObjectClass(object); if (class == CLASS_GROUP) class = DXGetGroupClass((Group)object); switch (class) { case CLASS_XFORM: if (! DXGetXformInfo((Xform)object, &child, 0)) return NULL; result = ReduceObject(child, factors, nDim); break; case CLASS_FIELD: result = ReduceField((Field)object, factors, nDim, NULL); break; case CLASS_COMPOSITEFIELD: result = ReduceCField((CompositeField)object, factors, nDim); break; case CLASS_SERIES: result = ReduceSeries((Series)object, factors, nDim); break; case CLASS_MULTIGRID: case CLASS_GROUP: result = ReduceGroup((Group)object, factors, nDim); break; default: DXSetError(ERROR_INVALID_DATA, "invalid data object"); result = NULL; break; } return result; } static Object ReduceGroup(Group group, Private factors, int nDim) { int i; Object child; i = 0; while ((child = DXGetEnumeratedMember(group, i++, NULL)) != NULL) if (! ReduceObject(child, factors, nDim)) return NULL; return (Object)group; } static Object ReduceCField(CompositeField cfield, Private factors, int nDim) { int i; Object child; Private omap; omap = ReducedMeshOffsets(cfield, factors); if (! omap) return NULL; i = 0; while ((child = DXGetEnumeratedMember((Group)cfield, i++, NULL)) != NULL) { if (DXGetObjectClass((Object)child) != CLASS_FIELD) { DXSetError(ERROR_INVALID_DATA, "invalid member of composite field"); return NULL; } if (! ReduceField((Field)child, factors, nDim, omap)) return NULL; } return (Object)cfield; } static Object ReduceSeries(Series series, Private factors, int nDim) { int i; Object child; i = 0; while ((child = DXGetSeriesMember((Series)series, i++, NULL)) != NULL) if (! ReduceObject(child, factors, nDim)) return NULL; return (Object)series; } typedef struct { Field field; int nDim; Private omap; Private factors; } FieldTask; static Object ReduceField(Field field, Private factors, int nDim, Private omap) { FieldTask task; if (DXEmptyField(field) || !DXGetComponentValue(field, "connections")) return (Object)field; task.field = field; task.nDim = nDim; task.factors = (Private)DXReference((Object)factors); if (omap) task.omap = (Private)DXReference((Object)omap); else task.omap = NULL; if (! DXAddTask(ReduceFieldTask, (Pointer)&task, sizeof(task), 1.0)) return NULL; else return (Object)field; } static Error ReduceFieldTask(Pointer ptr) { FieldTask *task; Field oField, rField; float *inFactors; int nDim; OriginMap omap; OMapElt **maps; Array inCons, outCons, comp; int nDimCons; int inCCounts[MAX_DIMENSION], origCCounts[MAX_DIMENSION]; int outCCounts[MAX_DIMENSION]; int origMeshOffsets[MAX_DIMENSION]; int reducedMeshOffsets[MAX_DIMENSION]; int meshOffsets[MAX_DIMENSION]; int i, j; float factors[MAX_DIMENSION]; int axes[MAX_DIMENSION]; int filter[MAX_DIMENSION]; char *name, *str; Object attr; Shoulder shoulder[MAX_DIMENSION]; int dep; task = (FieldTask *)ptr; rField = task->field; nDim = task->nDim; inFactors = (float *)DXGetPrivateData(task->factors); if (task->omap) { omap = (OriginMap)DXGetPrivateData(task->omap); maps = omap->map; } else { omap = NULL; maps = NULL; } DXChangedComponentValues(rField, "data"); if (NULL == (inCons = (Array)DXGetComponentValue(rField, "connections"))) { DXSetError(ERROR_MISSING_DATA, "connections missing"); return ERROR; } if (NULL == DXQueryGridConnections(inCons, &nDimCons, inCCounts)) { DXSetError(ERROR_INVALID_DATA, "field connections irregular"); return ERROR; } if (nDim != 1 && nDim != nDimCons) { DXSetError(ERROR_INVALID_DATA, "dimension disagreement"); return ERROR; } if (! DXGetMeshOffsets((MeshArray)inCons, meshOffsets)) { DXSetError(ERROR_INVALID_DATA, "field connections irregular"); return ERROR; } if (! DXQueryOriginalMeshExtents(rField, origMeshOffsets, origCCounts)) { for (i = 0; i < nDimCons; i++) { shoulder[i].left = shoulder[i].right = 0; } } else { for (i = 0; i < nDimCons; i++) { meshOffsets[i] += origMeshOffsets[i]; shoulder[i].left = origMeshOffsets[i]; shoulder[i].right = (inCCounts[i]-origCCounts[i])-shoulder[i].left; inCCounts[i] = origCCounts[i]; } } if (nDim != 1) for (i = 0; i < nDimCons; i++) if (inFactors[i] != 0) factors[i] = inFactors[i]; else factors[i] = 1.0; else for (i = 0; i < nDimCons; i++) if (inFactors[0] != 0) factors[i] = inFactors[0]; else factors[i] = 1.0; for (i = 0; i < nDimCons; i++) { if (inCCounts[i] < 1) { DXSetError(ERROR_INVALID_DATA, "<1 poins along %d axis", i); return ERROR; } if (factors[i] == 1 || inCCounts[i] == 1) { outCCounts[i] = inCCounts[i]; filter[i] = 1; } else { outCCounts[i] = (inCCounts[i] + 1) / factors[i]; if (outCCounts[i] < 2) outCCounts[i] = 2; filter[i] = ((int)ceil(factors[i])) | 0x01; } } /* * DXCopy input field for use as a template. Original field gets * reduced in place. */ oField = (Field)DXCopy((Object)rField, COPY_STRUCTURE); if (! oField) return ERROR; /* * Make regular output connections */ outCons = DXMakeGridConnectionsV(nDimCons, outCCounts); if (! outCons) goto error; if (! DXSetComponentValue(rField, "connections", (Object)outCons)) goto error; if (maps) { for (i = 0; i < nDimCons; i++) { for (j = 0; maps[i][j].inOffset != -1; j++) if (maps[i][j].inOffset == meshOffsets[i]) { reducedMeshOffsets[i] = maps[i][j].outOffset; break; } if (maps[i][j].inOffset == -1) { DXSetError(ERROR_INTERNAL, "mesh offset error"); goto error; } } if (! DXSetMeshOffsets((MeshArray)outCons, reducedMeshOffsets)) { DXSetError(ERROR_INTERNAL, "mesh offset error"); goto error; } } /* * Re-order axes for optimal reduction in irregular case */ AxesOrder(filter, axes, nDimCons); /* * Now reduce each other component that deps on positions */ i = -1; while ((comp=(Array)GetECV(oField, ++i, &name)) != NULL) { if (! strncmp(name, "original", 8)) { DXDeleteComponent(rField, name); continue; } if (! strncmp(name, "invalid ", 8) || ! strcmp(name, "connections")) continue; attr = GetECA(oField, i, NULL, "dep"); if (! attr) continue; str = DXGetString((String)attr); if (! str) goto error; if (! strcmp(str, "positions")) dep = POSITIONS_DEPENDENT; else if (! strcmp(str, "connections")) dep = CONNECTIONS_DEPENDENT; else continue; if (DXQueryGridPositions(comp, NULL, NULL, NULL, NULL)) comp = ReduceArrayReg(comp, nDimCons, shoulder, outCCounts); else { if (! strcmp(name, "positions")) comp = ReduceArrayIrreg(comp, nDimCons, NULL, shoulder, inCCounts, outCCounts, axes, dep); else comp = ReduceArrayIrreg(comp, nDimCons, filter, shoulder, inCCounts, outCCounts, axes, dep); } if (! comp) goto error; if (! DXSetComponentValue(rField, name, (Object)comp)) goto error; } if (DXGetComponentValue(oField, "invalid positions")) if (! InvalidateData(oField, rField, filter, INVALID_POSITIONS)) goto error; if (DXGetComponentValue(oField, "invalid connections")) if (! InvalidateData(oField, rField, filter, INVALID_CONNECTIONS)) goto error; DXDelete((Object)oField); DXDeleteComponent(rField, "box"); DXEndField(rField); DXDelete((Object)task->factors); DXDelete((Object)task->omap); return OK; error: DXDelete((Object)task->factors); DXDelete((Object)task->omap); DXDelete((Object)oField); return NULL; } static Array ReduceArrayReg(Array inArray, int nDimGrid, Shoulder *shoulder, int *outCountsGrid) { int i, j, nDimData, nextAxis; int inCountsData[MAX_DIMENSION], outCountsData[MAX_DIMENSION]; float origins[MAX_DIMENSION]; float inDeltas[MAX_DIMENSION], outDeltas[MAX_DIMENSION]; float inMax, outMax; Array outArray; float *odPtr, *idPtr; int *ocPtr, *icPtr; DXQueryGridPositions(inArray, &nDimData, inCountsData, origins, inDeltas); nextAxis = 0; for (i = 0; i < nDimData; i++) { if (inCountsData[i] != 1) { outCountsData[i] = outCountsGrid[nextAxis]; for (j = 0; j < nDimData; j++) origins[j] += shoulder[nextAxis].left * inDeltas[i*nDimData+j]; inCountsData[i] -= (shoulder[nextAxis].left + shoulder[nextAxis].right); nextAxis ++; } else outCountsData[i] = 1; } if (nextAxis != nDimGrid) { DXSetError(ERROR_INVALID_DATA, "regular data dimensions mismatch grid"); return NULL; } odPtr = outDeltas; idPtr = inDeltas; ocPtr = outCountsData; icPtr = inCountsData; for (i = 0; i < nDimData; i++) { /* * Make sure we do not reduce to less than two samples along * each axis. */ for (j = 0; j < nDimData; j++) { if (*ocPtr > 1) { *odPtr = ((double)((*icPtr) - 1) * (double)(*idPtr)) / (double)((*ocPtr) - 1); /* * A kludge to ensure that the output grid lies entirely inside * the input grid. */ inMax = (*idPtr) * ((*icPtr) - 1); if (inMax < 0.0) inMax = -inMax; outMax = (*odPtr) * ((*ocPtr) - 1); if (outMax < 0.0) outMax = -outMax; while(outMax > inMax) { *odPtr *= 0.9999999; outMax = (*odPtr) * ((*ocPtr) - 1); if (outMax < 0.0) outMax = -outMax; } } else *odPtr = 0.0; idPtr++; odPtr++; } icPtr++; ocPtr++; } outArray = DXMakeGridPositionsV(nDimData, outCountsData, origins, outDeltas); return outArray; } typedef struct { byte *inStart; byte *outStart; int inSkip; int outSkip; int increment; int limit; } FilterLoopCntl; static Array ReduceArrayIrreg(Array inDataArray, int nDim, int *filter, Shoulder *shoulders, int *inCounts, int *outCounts, int *axes, int dep) { Array outDataArray; int axis, a; int i, nElts; int tmpOutCounts[MAX_DIMENSION], tmpInCounts[MAX_DIMENSION]; int tmpLeft[MAX_DIMENSION], tmpRight[MAX_DIMENSION]; Pointer buf0, buf1, inBuf; Type type; Category cat; int rank, shape[MAX_DIMENSION]; int nItems, itemSize, bufSize; int status; int outItems; float strides[MAX_DIMENSION]; int filterLen; int outDataCounts[MAX_DIMENSION], inDataCounts[MAX_DIMENSION]; if (dep == POSITIONS_DEPENDENT) for (i = 0; i < nDim; i++) { outDataCounts[i] = outCounts[i]; inDataCounts[i] = inCounts[i]; } else for (i = 0; i < nDim; i++) { outDataCounts[i] = outCounts[i] - 1; inDataCounts[i] = inCounts[i] - 1; } outItems = 1; for (i = 0; i < nDim; i++) outItems *= outDataCounts[i]; if (!DXGetArrayInfo(inDataArray, &nItems, &type, &cat, &rank, shape)) return NULL; if (cat != CATEGORY_REAL || nItems <= 0) { DXSetError(ERROR_INVALID_DATA, "bad data component"); return NULL; } if (outItems == nItems) return inDataArray; itemSize = DXGetItemSize(inDataArray); nElts = itemSize / DXTypeSize(type); bufSize = itemSize * nItems; inBuf = DXGetArrayData(inDataArray); buf0 = inBuf; buf1 = NULL; for (i = 0; i < nDim; i++) { tmpInCounts[i] = tmpOutCounts[i] = inDataCounts[i] + shoulders[i].left + shoulders[i].right; if (outDataCounts[i] > 1) strides[i] = (inDataCounts[i]-1.0)/(outDataCounts[i]-1.0); else strides[i] = -1; tmpLeft[i] = shoulders[i].left; tmpRight[i] = shoulders[i].right; } for (a = 0; a < nDim; a++) { axis = axes[a]; tmpOutCounts[axis] = outDataCounts[axis]; if (tmpInCounts[axis] != tmpOutCounts[axis]) { if (! filter || filter[axis] <= 1) filterLen = 1; else filterLen = filter[axis]; bufSize = nElts * itemSize; for (i = 0; i < nDim; i++) bufSize *= tmpOutCounts[i]; buf1 = DXAllocate(bufSize); if (! buf1) { if (buf0 != inBuf) DXFree(buf0); return NULL; } status = ReduceSumAxis(buf0, buf1, filterLen, nDim, nElts, tmpInCounts, tmpOutCounts, tmpLeft, tmpRight, axis, strides, type); if (status != OK) { if (buf0 != inBuf) DXFree(buf0); DXFree(buf1); return NULL; } if (buf0 != inBuf) DXFree(buf0); buf0 = buf1; buf1 = NULL; } tmpInCounts[axis] = outDataCounts[axis]; tmpLeft[axis] = tmpRight[axis] = 0; } if (! (outDataArray = DXNewArrayV(type, cat, rank, shape))) return NULL; if (! DXAddArrayData(outDataArray, 0, outItems, buf0)) { DXDelete((Object)outDataArray); return NULL; } DXFree(buf0); return outDataArray; } #define AXISLOOP(type) \ { \ while(1) \ { \ byte *outPtr, *inPtr; \ float *sumElt; \ float *lastSumElt; \ type *outElt; \ type *inElt; \ float *bufElt, *bufPtr; \ \ outPtr = loop[0].outStart; \ inPtr = loop[0].inStart + startSkip; \ \ /* \ * Initialize accumulator and buffer \ */ \ memset(sumPtr, 0, nElts*sizeof(float)); \ memset(buf, 0, filter*nElts*sizeof(float)); \ \ bufI = filter; \ stride = filter; \ \ for (i = 0; i < loopLength; i++) \ { \ /* \ * Handle wrapping buf \ */ \ if (++bufI >= filter) \ { \ bufPtr = (float *)buf; \ bufI = 0; \ } \ \ sumElt = (float *)sumPtr; \ lastSumElt = (float *)lastSumPtr; \ inElt = (type *)inPtr; \ bufElt = bufPtr; \ \ for (j = 0; j < nElts; j++) \ { \ *lastSumElt = *sumElt; \ *sumElt += (((float)*inElt) - *bufElt); \ *bufElt = *inElt; \ lastSumElt++; \ sumElt++; \ inElt++; \ bufElt++; \ } \ \ /* \ * We may have to interpolate an output sample \ */ \ if (stride <= 1.0 && i != (loopLength-1)) \ { \ /* \ * Pointers to interpolate \ */ \ outElt = (type *)outPtr; \ sumElt = (float *)sumPtr; \ lastSumElt = (float *)lastSumPtr; \ \ if (stride == 1.0) \ for (j = 0; j < nElts; j++) \ *outElt++ = (type) \ (divisor * (*sumElt++) + round); \ else if (stride == 0.0) \ for (j = 0; j < nElts; j++) \ *outElt++ = (type) \ (divisor * (*lastSumElt++) + round); \ else \ for (j = 0; j < nElts; j++) \ { \ *outElt++ = (type) \ (divisor * (*lastSumElt + stride*(*sumElt-*lastSumElt)) + round); \ lastSumElt ++; \ sumElt ++; \ } \ \ stride += strideLength; \ outPtr += outAxisSkip; \ } \ \ if (strideLength > 0) \ stride -= 1.0; \ \ if (i >= incStart && i < incEnd) \ inPtr += inAxisSkip; \ \ bufPtr += nElts; \ } \ \ /* \ * For the last value along the axis. sum up the buffer \ */ \ if (++bufI >= filter) \ bufI = 0; \ \ bufPtr = buf + nElts*bufI; \ \ sumElt = sumPtr; \ bufElt = bufPtr; \ for (i = 0; i < nElts; i++) \ *sumElt++ = *bufElt++; \ \ for (i = 1; i < filter; i++) \ { \ if (++bufI >= filter) \ bufI = 0; \ \ bufPtr = buf + nElts*bufI; \ \ bufElt = bufPtr; \ sumElt = sumPtr; \ for (j = 0; j < nElts; j++) \ *sumElt++ += *bufElt++; \ } \ \ sumElt = sumPtr; \ outElt = (type *)outPtr; \ for (i = 0; i < nElts; i++) \ *outElt++ = (type)(divisor * *sumElt++ + round); \ \ if (nDim == 1) \ break; \ \ /* \ * Bump loops inner-to-outer. Terminate \ * either when out of loops or find a un-finished loop. \ */ \ l = loop + 1; \ for (i = 1; i < nDim; i++) \ { \ if (l->increment < l->limit) \ { \ l->increment++; \ break; \ } \ \ l++; \ } \ \ /* \ * If all loops complete, continue loop that iterates \ * the axis currently being filtered. Break loop that \ * handles outer loop. \ */ \ \ if (i == nDim) \ break; \ \ /* \ * Otherwise, bump the incompleted loop. \ */ \ l->inStart += l->inSkip; \ l->outStart += l->outSkip; \ \ /* \ * Now re-initialize lower axis loops. \ */ \ for (j = 0; j < i; j++) \ { \ loop[j].increment = 0; \ loop[j].inStart = l->inStart; \ loop[j].outStart = l->outStart; \ } \ \ /* \ * And continue the outer loop. \ */ \ } \ } static Error ReduceSumAxis(Pointer inBuf, Pointer outBuf, int filter, int nDim, int nElts, int *inCounts, int *outCounts, int *leftAxis, int *rightAxis, int axis, float *strides, Type type) { int inSkip[MAX_DIMENSION]; int outSkip[MAX_DIMENSION]; int overhang; int i, j; FilterLoopCntl loop[MAX_DIMENSION], *l; int limit; int inAxisSkip, outAxisSkip; float *buf, *sumPtr, *lastSumPtr; int bufI; float strideLength, stride; int left, right; int incStart, incEnd, startSkip, loopLength; int typeSize; float divisor = 1.0 / filter; float round; buf = sumPtr = lastSumPtr = NULL; typeSize = DXTypeSize(type); strideLength = strides[axis]; left = leftAxis[axis]; right = rightAxis[axis]; outSkip[nDim-1] = inSkip[nDim-1] = nElts * typeSize; for (i = (nDim-2); i >= 0; i--) { inSkip[i] = inSkip[i+1] * inCounts[i+1]; outSkip[i] = outSkip[i+1] * outCounts[i+1]; } /* * Make sure that the filter is odd in length */ filter |= 0x01; /* * Here's the setup for looping in each non-filter direction: */ l = loop; l->limit = inCounts[axis] - 1; l->inSkip = inSkip[axis]; l->outSkip = outSkip[axis]; l->inStart = (byte *)inBuf; l->outStart = (byte *)outBuf; l->increment = 0; l++; for (i = 0; i < nDim; i++) if (i != axis) { l->limit = inCounts[i] - 1; l->inSkip = inSkip[i]; l->outSkip = outSkip[i]; l->inStart = (byte *)inBuf; l->outStart = (byte *)outBuf; l->increment = 0; l++; } /* * Get buffer and accumulator that will span the * box filter in the current filter axis. */ buf = (float *)DXAllocateLocal(filter*nElts*sizeof(float)); if (! buf) { DXResetError(); buf = (float *)DXAllocate(filter*nElts*sizeof(float)); } if (! buf) goto error; sumPtr = (float *)DXAllocateLocal(nElts*sizeof(float)); if (! sumPtr) { DXResetError(); sumPtr = (float *)DXAllocate(nElts*sizeof(float)); } if (! sumPtr) goto error; lastSumPtr = (float *)DXAllocateLocal(nElts*sizeof(float)); if (! lastSumPtr) { DXResetError(); lastSumPtr = (float *)DXAllocate(nElts*sizeof(float)); } if (! lastSumPtr) goto error; limit = loop[0].limit; inAxisSkip = loop[0].inSkip; outAxisSkip = loop[0].outSkip; overhang = (filter - 1) >> 1; /* * Basic loop operates on the total number of samples accumulated * to produce the output line. This is the number of input samples * minus extra added for shoulders plus the two overhangs. */ loopLength = (inCounts[axis] - (left + right)) + (overhang + overhang); /* * If the left shoulder is wider than the overhang, we will need to skip * unused samples. Since there is at least overhang shoulder samples, * we start bumping the input pointer immediately. Otherwise, we start * operating on the first sample. In this case, we begin bumping the * pointer when we have replicated (overhang - left) samples. */ if (left > overhang) { startSkip = (left - overhang) * inAxisSkip; incStart = 0; } else { startSkip = 0; incStart = overhang - left; } /* * If the right shoulder is shorter than the overhang, we cease * bumping the input pointer when we have exhausted the shoulder. */ if (right < overhang) { incEnd = (loopLength - 1) - (overhang - right); } else { incEnd = loopLength - 1; } switch(type) { case TYPE_BYTE: round = 0.5; AXISLOOP(byte); break; case TYPE_UBYTE: round = 0.5; AXISLOOP(ubyte); break; case TYPE_SHORT: round = 0.5; AXISLOOP(short); break; case TYPE_USHORT: round = 0.5; AXISLOOP(ushort); break; case TYPE_INT: round = 0.5; AXISLOOP(int); break; case TYPE_UINT: round = 0.5; AXISLOOP(uint); break; case TYPE_FLOAT: round = 0.0; AXISLOOP(float); break; case TYPE_DOUBLE: round = 0.0; AXISLOOP(double); break; } DXFree((Pointer)buf); DXFree((Pointer)sumPtr); DXFree((Pointer)lastSumPtr); return OK; error: DXFree((Pointer)sumPtr); DXFree((Pointer)lastSumPtr); DXFree((Pointer)buf); return ERROR; } static void AxesOrder(int *filter, int *axes, int nDim) { int flags[MAX_DIMENSION]; int max, i, j; memset(flags, 1, nDim*sizeof(int)); for (i = 0; i < nDim; i++) { max = -1; for (j = 0; j < nDim; j++) if (flags[j] && max < filter[j]) { axes[i] = j; max = filter[j]; } flags[axes[i]] = 0; } } typedef struct node { struct node *next; int offset; } Node; #define MAX_NODES 2048 static Private ReducedMeshOffsets(CompositeField cfield, Private pfactors) { Field part; Node lists[MAX_DIMENSION]; int axisCounts[MAX_DIMENSION]; Node nodeStash[MAX_NODES]; int growthOffsets[MAX_DIMENSION]; int meshOffsets[MAX_DIMENSION]; int nextn, nDim; int i, j, k, l; Node *n, *last, *next; Array conn; int first, offset; OriginMap omap; OMapElt **maps; Private omapObject; float *factors; omap = NULL; if (! pfactors) return NULL; factors = (float *)DXGetPrivateData(pfactors); if (! factors) return NULL; first = 1; nextn = 0; /* * First, look through the partitions looking for those that lie * along the axes. */ i = 0; while(NULL != (part = (Field)DXGetPart((Object)cfield, i++))) { conn = (Array)DXGetComponentValue(part, "connections"); if (! conn) { DXSetError(ERROR_MISSING_DATA, "partition with no connections found"); return NULL; } if (first) { if (! DXQueryGridConnections(conn, &nDim, NULL)) { DXSetError(ERROR_INVALID_DATA, "irregular connections found"); return NULL; } /* * Set up empty linked lists for the insertion of partition * origins. */ for (j = 0; j < nDim; j++) { lists[j].next = NULL; lists[j].offset = 0; axisCounts[j] = 1; /* for 0 slot */ } first = 0; } if (! DXGetMeshOffsets((MeshArray)conn, meshOffsets)) return NULL; if (DXQueryOriginalMeshExtents(part, growthOffsets, NULL)) for (j = 0; j < nDim; j++) meshOffsets[j] += growthOffsets[j]; for (j = 0; j < nDim; j++) { offset = meshOffsets[j]; if (offset == 0) continue; last = &lists[j]; next = last->next; while (next != NULL) { if (next->offset >= offset) break; last = next; next = next->next; } if (next && next->offset == offset) continue; n = nodeStash + nextn++; if (nextn >= MAX_NODES) { DXSetError(ERROR_INTERNAL, "exhausted supply of list nodes"); return NULL; } n->offset = offset; n->next = next; last->next = n; axisCounts[j] ++; } } omap = (OriginMap)DXAllocate(sizeof(struct originMap)); if (! omap) goto error; omap->map = (OMapElt **)DXAllocateZero(nDim*sizeof(struct oMapElt *)); if (! omap->map) goto error; omap->nDim = nDim; maps = omap->map; for (i = 0; i < nDim; i++) { maps[i]=(OMapElt *)DXAllocate((axisCounts[i]+1)*sizeof(struct originMap)); if (! maps[i]) goto error; j = k = 0; next = &lists[i]; while (next) { maps[i][j].inOffset = next->offset; maps[i][j].outOffset = k; last = next; next = next->next; if (! next) break; if (factors[i] > 1.0) { /* * The input count in the last partition is * ((nextOffset - lastOffset) + 1). The output * counts for that (see ReduceFieldTask) is * (counts + 1) / factor. The offset of the next * partition is last + (outCounts-1) */ l = (((next->offset - last->offset) + 2)/factors[i]); if (l < 2) l = 2; l -= 1; } else l = next->offset - last->offset; k += l; j++; } maps[i][j+1].inOffset = -1; maps[i][j+1].outOffset = -1; } omapObject = DXNewPrivate((Pointer)omap, DeleteOriginMapObject); if (! omapObject) goto error; return omapObject; error: if (omap) DeleteOriginMapObject((Pointer) omap); return NULL; } static Error DeleteOriginMapObject(Pointer p) { int i; OriginMap o; if (! p) { DXSetError(ERROR_INTERNAL,"DeleteOriginMapObject called on NULL pointer"); return ERROR; } o = (OriginMap)p; for (i = 0; i < o->nDim; i++) DXFree((Pointer)o->map[i]); DXFree((Pointer)o->map); DXFree(p); return OK; } static Error DeleteFactorsObject(Pointer p) { DXFree(p); return OK; } #define INDICES(ref, indices, strides, ndim) \ { \ int i, j = (ref); \ for (i = 0; i < (ndim)-1; i++) { \ (indices)[i] = j / (strides)[i]; \ j = j % (strides)[i]; \ } \ (indices)[i] = j; \ } #define REFERENCE(indices, ref, strides, ndim) \ { \ int i, j = 0; \ for (i = 0; i < (ndim)-1; i++) \ j += (indices)[i]*(strides)[i]; \ (ref) = j + (indices)[(ndim)-1]; \ } static Error InvalidateData(Field original, Field reduced, int *filter, int invalid) { InvalidComponentHandle inInvalids = NULL, outInvalids = NULL; Array inArray, outArray; int i, nDim, index; int inCounts[10], inStrides[10]; int outCounts[10], outStrides[10]; float axisScale[10]; int indices[10], min[10], max[10]; inArray = (Array)DXGetComponentValue(original, "connections"); if (! inArray) { DXSetError(ERROR_INTERNAL, "original has no connections"); goto error; } if (! DXQueryGridConnections(inArray, &nDim, inCounts)) { DXSetError(ERROR_INTERNAL, "original positions are not grid"); goto error; } outArray = (Array)DXGetComponentValue(reduced, "connections"); if (! outArray) { DXSetError(ERROR_INTERNAL, "reduced has no positions"); goto error; } if (! DXQueryGridConnections(outArray, NULL, outCounts)) { DXSetError(ERROR_INTERNAL, "reduced connections are not grid"); goto error; } if (invalid == INVALID_POSITIONS) { inInvalids = DXCreateInvalidComponentHandle((Object)original, NULL, "positions"); outInvalids = DXCreateInvalidComponentHandle((Object)reduced, NULL, "positions"); } else { for (i = 0; i < nDim; i++) { inCounts[i] -= 1; outCounts[i] -= 1; } inInvalids = DXCreateInvalidComponentHandle((Object)original, NULL, "connections"); outInvalids = DXCreateInvalidComponentHandle((Object)reduced, NULL, "connections"); } inStrides[nDim-1] = outStrides[nDim-1] = 1; for (i = nDim-2; i >= 0; i--) { inStrides[i] = inStrides[i+1] * inCounts[i+1]; outStrides[i] = outStrides[i+1] * outCounts[i+1]; } for (i = 0; i < nDim; i++) if (inCounts[i] == 0) axisScale[i] = 1.0; else axisScale[i] = ((float)outCounts[i]) / ((float)inCounts[i]); if (! inInvalids || ! outInvalids) goto error; DXSetAllValid(outInvalids); DXInitGetNextInvalidElementIndex(inInvalids); while (-1 != (index = DXGetNextInvalidElementIndex(inInvalids))) { INDICES(index, indices, inStrides, nDim); for (i = 0; i < nDim; i++) { int f = filter[i] >> 1; /* * Get the range of input samples along the axis that the * invalid input entry affects. This is the range of * input points that the invalid point affects in the * filtering process. */ min[i] = indices[i] - f; if (min[i] < 0) min[i] = 0; max[i] = indices[i] + f; if (max[i] >= inCounts[i]) max[i] = inCounts[i] - 1; /* * Now transform these input indices to output indices * by scaling. This gives us the range of output points * that are affected by this input point. */ min[i] = floor(min[i] * axisScale[i]); max[i] = floor(max[i] * axisScale[i]); /* * initialize the loop */ indices[i] = min[i]; } /* * Now loop to set the appropriate output points invalid */ for ( ;; ) { int r; REFERENCE(indices, r, outStrides, nDim); if (! DXSetElementInvalid(outInvalids, r)) goto error; for (i = 0; i < nDim; i++) { indices[i] ++; if (indices[i] <= max[i]) break; indices[i] = min[i]; } if (i == nDim) break; } } if (! DXSaveInvalidComponent(reduced, outInvalids)) goto error; DXFreeInvalidComponentHandle(outInvalids); DXFreeInvalidComponentHandle(inInvalids); outInvalids = inInvalids = NULL; return OK; error: DXFreeInvalidComponentHandle(outInvalids); DXFreeInvalidComponentHandle(inInvalids); return ERROR; }