/***********************************************************************/ /* 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 #define OUTPUT_TYPE_FLOAT 0 /* set to 1 to make float */ /* notes to me: * in making 2d, 3d and Nd histograms, the following changes are needed. * the num bins, the min, max and ?in/out? flags need to be per dim now * what is the median for a 2d histogram? median per dim? then it too * must be a vector now. * * requested new function: average bin value. automatically created? * flag to ask for? add a separate component named "bin average"? * * related function - user request to bin data values based on one * data value and get the average of a second data value - the shape * of the histogram will be based on the original data component, but * the other value can be marked and used to color the data or display * a graph, etc. can this be done with a 2d input data, where the nbins * is normal for X and 1 for Y? * * more distantly related function - sum across rows and columns of * Nd arrays. stats as well? can this be done with a macro which * slices data first then stats module? * */ /* just too long for me to type, sorry */ #define ICH InvalidComponentHandle /* for floating point number compares - this is a small number, * but large enough to be a normalized floating point number. * (denormalized numbers cause a floating point trap on the i860.) */ #define FUZZ ((float) 1e-37) /* max supported is 32 dimensional */ #define MAXRANK 32 /* number of parallel tasks * (number of processors - 1) */ #define PFACTOR 4 /* only used when computing histogram groups in parallel */ struct parinfo { Object h; /* the histogram itself */ char *name; /* the group member name */ int didwork; /* global copy of didwork flag */ /* slicecnt here? */ }; /* values for the 'XXXset' flags below - * indicates whether or how the values were determined. */ #define NOTSET 0 /* no value yet */ #define EXPLICIT 1 /* user-given input parameter to module */ #define BYOBJECT 2 /* statistics of input object used */ /* parm list */ struct histinfo { Object o; /* use data component of this object */ int dim; /* dimensionality of data/histogram */ int bins; /* number of histogram bins */ int *md_bins; /* list, used instead of bins if dim > 1 */ int binset; /* see flag values above */ float min; /* value of minimum of first bin */ float *md_min; /* list for dim > 1 */ int minset; /* flag */ float max; /* value of maximum of last bin */ float *md_max; /* list for dim > 1 */ int maxset; /* flag */ int inoutlo; /* what to do with data below first bin */ int *md_inoutlo; /* list for dim > 1 */ int inoutloset; /* flag */ int inouthi; /* what to do with data above last bin */ int *md_inouthi; /* list for dim > 1 */ int inouthiset; /* flag */ int didwork; /* whether any points fell into any bins */ float median; /* the approximate data median */ float *md_median; /* vector for dim > 1 */ int maxparallel; /* max number of parallel tasks; based on nproc */ int goneparallel; /* number of parallel tasks already spawned */ ICH invalid; /* if set, there are invalid data values */ struct parinfo *p; /* if computing in parallel, where to put the output */ }; /* things which need a list of values, one per dim */ struct perdim { int *binsplus1; /* dim */ float *binsize; int *thisbin; float *dminlist; float *dmaxlist; int **slicecnt; float *deltas; /* dim * dim */ }; /* this routine could be in libdx */ Object _dxfHistogram(Object o, int bins, float *min, float *max, int inout); /* private routines */ static Error Histogram_Wrapper(Pointer); static Error Histogram_Wrapper2(Pointer); static Object Object_Histogram(struct histinfo *hp); static Group Group_Histogram(struct histinfo *hp); static Field Field_Histogram(struct histinfo *hp); static Field Array_Histogram(struct histinfo *hp); static Field EmptyHistogram(struct histinfo *hp); static Field EmptyHistogram1D(struct histinfo *hp); static Field HistogramND(struct histinfo *hp, int empty); static Error HasInvalid(Field f, char *, ICH *ih); static int GetVectorShape(int rank, int *shape); static Error VectorLength(Object o, int *len, int onlyone); static Error VectorExtract(Object o, Type expected_t, int len, Pointer ptr); static Error VectorAllocAndExtract(Object o, Type expected_t, int len, Pointer *ptr); static Error ScalarAlloc(Pointer inval, Type t, int len, Pointer *retval); static Error CopyHistinfo(struct histinfo *src, struct histinfo *dst); static void FreeHistinfo(struct histinfo *hp); static Error ObjectVectorStats(Object o, int *vlen, float **minlist, float **maxlist); static Error VectorStats(Object o, int knownlength, int *vlen, int initialize, float **minlist, float **maxlist); static Error init_perdim(struct perdim **pdp, int dim); static Error bins_perdim(struct perdim *pdp, int dim, int *binsize); static void free_perdim(struct perdim *pdp, int dim); static Error FullAllocHistinfo(struct histinfo *hp, int datadim, int setdef); #define INTERNALERROR DXSetError(ERROR_INTERNAL, \ "unexpected error, file %s line %d", __FILE__, __LINE__) #define DEFAULTBINS 100 #ifndef MAX #define MAX(a,b) ((a) > (b) ? (a) : (b)) #endif /* 0 = input * 1 = bins * 2 = min * 3 = max * 4 = in/out * * 0 = output * 1 = median */ int m_Histogram(Object *in, Object *out) { struct histinfo h; float tmin, tmax; int multidim; int prevdim; int i, flag; /* the idea here is to look at each of the input parms. * if they are all integers, this is a simple, 1d histogram. * if any are vector (number lists), then all the ones which * are vectors must match the vector length of the data component. * any which are still given as scalars are ok and apply to all * the dimensions. */ /* initialize histinfo struct to zero, assume 1d and set defaults. */ memset((char *)&h, '\0', sizeof(h)); h.dim = 1; h.bins = DEFAULTBINS; h.min = 0.0; h.max = 1.0; h.inoutlo = 0; h.inouthi = 0; h.median = 0.0; h.maxparallel = DXProcessors(0); if (h.maxparallel > 1) h.maxparallel = (h.maxparallel-1) * PFACTOR; prevdim = 1; /* input object to histogram */ if (!in[0]) { DXSetError(ERROR_BAD_PARAMETER, "#10000", "input"); return ERROR; } else h.o = in[0]; /* number of bins */ if (in[1]) { /* simple integer */ if (DXExtractInteger(in[1], &h.bins)) { if (h.bins <= 0) { DXSetError(ERROR_BAD_PARAMETER, "#10525", "number of bins"); return ERROR; } h.binset = EXPLICIT; } else if (VectorLength(in[1], &multidim, 1)) { /* vector list? */ if (!VectorAllocAndExtract(in[1], TYPE_INT, multidim, (Pointer *)&h.md_bins)) goto error; for (i=0; i 1 && prevdim != multidim) { DXSetError(ERROR_BAD_PARAMETER, "#10085", "bin minimum", prevdim); goto error; } if (prevdim == 1) prevdim = multidim; h.minset = EXPLICIT; } else if (ObjectVectorStats(in[2], &multidim, &h.md_min, &h.md_max)) { if (prevdim > 1 && prevdim != multidim) { DXSetError(ERROR_BAD_PARAMETER, "#10085", "bin minimum", prevdim); goto error; } if (prevdim == 1) prevdim = multidim; h.minset = BYOBJECT; } else { DXSetError(ERROR_BAD_PARAMETER, "#10562", "bin mininum"); goto error; } } /* max bin value. either an explicit number, an list, or object whose * data maximum will be used, or the data maximum from the previous * parameter if parm 2 was an object and parm 3 is NULL. */ if (in[3]) { if (DXExtractFloat(in[3], &h.max)) { h.maxset = EXPLICIT; } else if (VectorLength(in[3], &multidim, 1)) { /* vector list? */ if (!VectorAllocAndExtract(in[3], TYPE_FLOAT, multidim, (Pointer *)&h.md_max)) goto error; /* if an earlier parm has already set the length of the list, * make sure this one matches in size. if not, set the length * here and take care of replicating the values for earlier * parms at the end. */ if (prevdim > 1 && prevdim != multidim) { DXSetError(ERROR_BAD_PARAMETER, "#10085", "bin maximum", prevdim); goto error; } if (prevdim == 1) prevdim = multidim; h.maxset = EXPLICIT; } else if (ObjectVectorStats(in[3], &multidim, NULL, &h.md_max)) { if (prevdim > 1 && prevdim != multidim) { DXSetError(ERROR_BAD_PARAMETER, "#10085", "bin maximum", prevdim); goto error; } if (prevdim == 1) prevdim = multidim; h.maxset = BYOBJECT; } else { DXSetError(ERROR_BAD_PARAMETER, "#10562", "bin maximum"); goto error; } } else if (h.minset == BYOBJECT) { h.max = tmax; h.maxset = BYOBJECT; } /* whether to ignore counts outside the bins, or add them to the * first and last bins. if using the data statistics for the * bin max, set inouthi flag to 1, so data equal to the data max * will be included in the last bin. */ if (in[4]) { if (DXExtractInteger(in[4], &h.inoutlo)) { if (h.inoutlo != 0 && h.inoutlo != 1) { DXSetError(ERROR_BAD_PARAMETER, "#10070", "outlier flag"); goto error; } h.inouthi = h.inoutlo; h.inoutloset = EXPLICIT; h.inouthiset = EXPLICIT; } else if (VectorLength(in[4], &multidim, 1)) { /* vector list? */ if (!VectorAllocAndExtract(in[4], TYPE_INT, multidim, (Pointer *)&h.md_inoutlo)) goto error; for (i=0; i 1 && prevdim != multidim) { DXSetError(ERROR_BAD_PARAMETER, "#10076", "outlier flag", prevdim); goto error; } if (prevdim == 1) prevdim = multidim; h.inoutloset = EXPLICIT; h.inouthiset = EXPLICIT; } else { DXSetError(ERROR_BAD_PARAMETER, "#10075", "outlier flag"); goto error; } } else if (!h.maxset || h.maxset == BYOBJECT) { h.inouthi = 1; h.inouthiset = BYOBJECT; } multidim = prevdim; h.dim = multidim; if (!FullAllocHistinfo(&h, multidim, 1)) goto error; /* one last error check - is max smaller or equal to min? */ if (h.minset && h.maxset) { if (h.dim == 1) { if (h.md_min[0] >= h.md_max[0]) { DXSetError(ERROR_BAD_PARAMETER, "#10185", "min", "max"); goto error; } } else { for (i=0; i= h.md_max[i]) { DXSetError(ERROR_BAD_PARAMETER, "#10187", "min", i, "max", i); goto error; } } } /* do the work here */ out[0] = Object_Histogram(&h); /* if no error and the output is one histogram (as opposed to a group * of histograms), set the second output to the estimated median. * * it is inaccurate in proportion to the bin width. */ if (out[0] && DXGetObjectClass(out[0]) == CLASS_FIELD) { if (h.dim == 1) out[1] = (Object)DXNewArray(TYPE_FLOAT, CATEGORY_REAL, 0); else out[1] = (Object)DXNewArray(TYPE_FLOAT, CATEGORY_REAL, 1, h.dim); DXAddArrayData((Array)out[1], 0, 1, (Pointer)h.md_median); } else out[1] = NULL; /* it would be nice to know if the second output is being used, * because then i'd issue a warning if the input was a group and * the output had no median because of that. but a warning every * time seems annoying. */ /* is it an error if the data is all out of range? and if not, * then what is the median? i'm currently returning NULL for * median, but a field for histogram, and not setting error. */ if (out[0] && h.didwork == 0) DXWarning("no data values in histogram range"); FreeHistinfo(&h); return (out[0] ? OK : ERROR); error: FreeHistinfo(&h); out[0] = out[1] = NULL; return(ERROR); } /* general histogram routine. not used in m_Histogram; it just packages * the parms into a block and calls the internal routine to make it * easier for other libDX users to call the routine. */ Object _dxfHistogram(Object o, int bins, float *min, float *max, int inout) { struct histinfo h; Object retval = NULL; memset((char *)&h, '\0', sizeof(h)); h.o = o; h.dim = 1; if (bins > 0) { h.bins = bins; h.binset = EXPLICIT; } if (min) { h.min = *min; h.minset = EXPLICIT; } if (max) { h.max = *max; h.maxset = EXPLICIT; } h.inoutlo = inout; h.inoutloset = EXPLICIT; /* what about !max here ? */ h.inouthi = inout; h.inouthiset = EXPLICIT; if (!FullAllocHistinfo(&h, 1, 1)) goto done; retval = Object_Histogram(&h); done: FreeHistinfo(&h); return retval; } /* if the object is a Series, CompositeField or Field, make one histogram * from it. if it is a general Group, make one histogram for each member. */ static Object Object_Histogram(struct histinfo *hp) { Object subo, old, cpy; switch(DXGetObjectClass(hp->o)) { case CLASS_ARRAY: return (Object)Array_Histogram(hp); case CLASS_FIELD: return (Object)Field_Histogram(hp); case CLASS_GROUP: switch(DXGetGroupClass((Group)hp->o)) { case CLASS_COMPOSITEFIELD: /* have to make a copy here */ cpy = DXCopy(hp->o, COPY_STRUCTURE); if (!cpy) return NULL; old = hp->o; hp->o = cpy; if (!DXInvalidateDupBoundary(hp->o)) { /* delete copy and restore original object */ DXDelete(hp->o); hp->o = old; return NULL; } subo = (Object)Field_Histogram(hp); DXDelete(cpy); hp->o = old; return subo; case CLASS_MULTIGRID: case CLASS_SERIES: return (Object)Field_Histogram(hp); default: return (Object)Group_Histogram(hp); } case CLASS_SCREEN: if (!DXGetScreenInfo((Screen)hp->o, &subo, NULL, NULL)) return NULL; hp->o = subo; return Object_Histogram(hp); case CLASS_XFORM: if (!DXGetXformInfo((Xform)hp->o, &subo, NULL)) return NULL; hp->o = subo; return Object_Histogram(hp); case CLASS_CLIPPED: if (!DXGetClippedInfo((Clipped)hp->o, &subo, NULL)) return NULL; hp->o = subo; return Object_Histogram(hp); default: DXSetError(ERROR_BAD_PARAMETER, "#10190", "input"); return NULL; } /* not reached */ } /* the input is a generic group; call Histogram for each member. */ static Group Group_Histogram(struct histinfo *hp) { Group g; Group newg = NULL; Object subo; Object h; char *name; int i, members; struct parinfo *pinfo = NULL; /* handle empty groups, and groups of only 1 member */ g = (Group)hp->o; for (members=0; subo = DXGetEnumeratedMember(g, members, &name); members++) ; if (members == 0) return (Group)EmptyHistogram(hp); newg = (Group)DXCopy(hp->o, COPY_HEADER); if (!newg) return NULL; if (members == 1) { subo = DXGetEnumeratedMember(g, 0, &name); if (!subo) goto error; hp->o = subo; h = Object_Histogram(hp); if (!h) goto error; /* if the member is named, add it by name, else add by number */ if (! (name ? DXSetMember(newg, name, h) : DXSetEnumeratedMember(newg, 0, h))) goto error; return newg; } /* if we have groups nested in groups, make sure we don't start * too many parallel tasks. */ if (hp->goneparallel > hp->maxparallel) { for (i=0; subo = DXGetEnumeratedMember(g, i, &name); i++) { hp->o = subo; h = Object_Histogram(hp); if (!h) goto error; /* if the member is named, add it by name, else add by number */ if (! (name ? DXSetMember(newg, name, h) : DXSetEnumeratedMember(newg, i, h))) goto error; } return newg; } /* for each member, construct a histogram in parallel, and * add them to the group later. */ pinfo = (struct parinfo *)DXAllocateZero(members * sizeof(struct parinfo)); if (!pinfo) goto error; if (!DXCreateTaskGroup()) goto error; for (i=0; subo = DXGetEnumeratedMember(g, i, &name); i++) { hp->o = subo; hp->p = &pinfo[i]; hp->goneparallel++; if (!DXAddTask(Histogram_Wrapper, (Pointer)hp, sizeof(struct histinfo), 0.0)) goto error; } if (!DXExecuteTaskGroup()) goto error; hp->goneparallel -= members; /* now add accumulated histograms to the parent group */ for (i=0; i < members; i++) { /* if the member is named, add it by name, else add by number */ if (! (pinfo[i].name ? DXSetMember(newg, pinfo[i].name, pinfo[i].h) : DXSetEnumeratedMember(newg, i, pinfo[i].h))) goto error; hp->didwork += pinfo[i].didwork; } DXFree((Pointer)pinfo); return newg; error: DXDelete((Object)newg); DXFree((Pointer)pinfo); return NULL; } /* for going parallel. the wrapper routine which calls Object_Histogram * and puts the results in the right place. */ static Error Histogram_Wrapper(Pointer a) { struct histinfo *hp = (struct histinfo *)a; struct histinfo h2; if (!CopyHistinfo(hp, &h2)) return ERROR; hp->p->h = Object_Histogram(&h2); if (!hp->p->h) goto error; hp->p->didwork += h2.didwork; /* DXFreeInvalidComponentHandle(hp->invalid); */ /* not needed here? */ FreeHistinfo(&h2); return OK; error: FreeHistinfo(&h2); return ERROR; } /* for going parallel. the wrapper routine which calls Array_Histogram * and puts the results in the right place. */ static Error Histogram_Wrapper2(Pointer a) { struct histinfo *hp = (struct histinfo *)a; struct histinfo h2; if (!CopyHistinfo(hp, &h2)) return ERROR; hp->p->h = (Object)Array_Histogram(&h2); if (!hp->p->h) goto error; hp->p->didwork += h2.didwork; DXFreeInvalidComponentHandle(hp->invalid); FreeHistinfo(&h2); return OK; error: FreeHistinfo(&h2); return ERROR; } /* this routine takes case of simple fields, and composite fields. in * both cases, the output is 1 histogram in a simple field. * * the code goes parallel on partitioned data. */ static Field Field_Histogram(struct histinfo *hp) { Object subo; Array ao, ao1; Type type; Group g; struct parinfo *pinfo = NULL; Field *f = NULL; Field h = NULL; int isparallel = 0; int items; int multidim; int members, i, j; int totalcnt; #if OUTPUT_TYPE_FLOAT float *accum, *add; #else int *accum, *add; #endif float binsize, tmin, tmax; float gather, median; /* simple fields - histogram the data component */ if (DXGetObjectClass(hp->o) == CLASS_FIELD) { /* ignore empty fields/partitions */ if (DXEmptyField((Field)hp->o)) return EmptyHistogram(hp); ao = (Array)DXGetComponentValue((Field)hp->o, "data"); if (!ao) { /* no data component */ DXSetError(ERROR_MISSING_DATA, "#10240", "data"); return NULL; } /* this has to happen here while o is still the field */ if (!HasInvalid((Field)hp->o, "data", &hp->invalid)) return NULL; hp->o = (Object)ao; subo = (Object)Array_Histogram(hp); DXFreeInvalidComponentHandle(hp->invalid); return (Field)subo; } /* not a simple field - a composite field or series. process each * member in parallel and sum them at the end. */ /* count the members, and return if none */ if (!DXGetMemberCount((Group)hp->o, &members) || (members == 0)) return EmptyHistogram(hp); /* the default min and max bin values must be for the entire group. * set them here if they aren't already set. */ if (!hp->minset || !hp->maxset) { if (!hp->minset && !hp->maxset) { if (!ObjectVectorStats(hp->o, &multidim, &hp->md_min, &hp->md_max)) return NULL; hp->minset = BYOBJECT; hp->maxset = BYOBJECT; } if (!hp->minset) { if (!ObjectVectorStats(hp->o, &multidim, &hp->md_min, NULL)) return NULL; hp->minset = BYOBJECT; } if (!hp->maxset) { if (!ObjectVectorStats(hp->o, &multidim, NULL, &hp->md_max)) return NULL; hp->maxset = BYOBJECT; } if ((hp->dim != multidim) && !FullAllocHistinfo(hp, multidim, 0)) return NULL; } /* if the binsize hasn't been set yet, base it on the data type. * byte data default is data range (max 256); everything else is 100 bins. */ if (!hp->binset) { if (DXGetType(hp->o, &type, NULL, NULL, NULL) == NULL) return NULL; for (i=0; idim; i++) hp->md_bins[i] = (type != TYPE_UBYTE) ? DEFAULTBINS : hp->md_max[i] - hp->md_min[i] + 1; hp->binset = BYOBJECT; } /* allocate space to hold the array of histogram fields. * we are going to use the first field to return the data, and * extract the data from the others and add it in and then * discard the rest of the fields. */ pinfo = (struct parinfo *)DXAllocateZero(members * sizeof(struct parinfo)); if (!pinfo) return NULL; /* below here, goto error instead of just returning */ f = (Field *)DXAllocateZero(members * sizeof(Field)); if (!f) goto error; /* iterate thru the members - this can happen in parallel. */ if (hp->goneparallel < hp->maxparallel) { if (!DXCreateTaskGroup()) goto error; else isparallel++; } g = (Group)hp->o; for (i=0; subo = DXGetEnumeratedMember(g, i, NULL); i++) { hp->o = subo; hp->p = &pinfo[i]; /* if not a field, histogram it with the higher level routine. * what if it's a group???? */ if (DXGetObjectClass(subo) != CLASS_FIELD) { if (isparallel) { if (!DXAddTask(Histogram_Wrapper, (Pointer)hp, sizeof(struct histinfo), 0.0)) goto error; hp->goneparallel++; } else { f[i] = (Field)Object_Histogram(hp); if (DXGetError() != ERROR_NONE) goto error; } continue; } /* ignore empty fields */ if (DXEmptyField((Field)subo)) continue; ao = (Array)DXGetComponentValue((Field)subo, "data"); if (!ao) { DXSetError(ERROR_MISSING_DATA, "#10240", "data"); /* no data component */ goto error; } if (!HasInvalid((Field)subo, "data", &hp->invalid)) goto error; hp->o = (Object)ao; if (isparallel) { if (!DXAddTask(Histogram_Wrapper2, (Pointer)hp, sizeof(struct histinfo), 0.0)) goto error; hp->goneparallel++; } else { f[i] = Array_Histogram(hp); if (DXGetError() != ERROR_NONE) goto error; DXFreeInvalidComponentHandle(hp->invalid); } } if (isparallel) { if (!DXExecuteTaskGroup()) goto error; hp->goneparallel -= members; /* put histograms where they are expected to be */ for (i=0; i < members; i++) { f[i] = (Field)pinfo[i].h; hp->didwork += pinfo[i].didwork; } } /* coalesce the bins into one histogram */ accum = NULL; totalcnt = 0; for (i=0; i < members; i++) { if (!f[i] || DXEmptyField(f[i])) continue; ao1 = (Array)DXGetComponentValue(f[i], "data"); if (!ao1) { INTERNALERROR; goto error; } /* use the data array from the first field to accumulate into */ if (!accum) { DXGetArrayInfo(ao1, &items, NULL, NULL, NULL, NULL); #if OUTPUT_TYPE_FLOAT accum = (float *)DXGetArrayData(ao1); #else accum = (int *)DXGetArrayData(ao1); #endif if (!accum) { INTERNALERROR; goto error; } /* we will accumulate the total histogram into this field. * continue here so we don't accumulate this field twice. */ h = f[i]; for (j=0; j < items; j++) totalcnt += accum[j]; /* XXX md - can't recompute totalcnt like this. need slicecnt * for Nd median. */ continue; } #if OUTPUT_TYPE_FLOAT add = (float *)DXGetArrayData(ao1); #else add = (int *)DXGetArrayData(ao1); #endif if (!add) { INTERNALERROR; goto error; } for (j=0; j < items; j++) { accum[j] += add[j]; totalcnt += add[j]; /* XXX md - ditto totalcnt comment. accum is ok */ } } /* delete all but the accumulator field, and recompute the median */ for (i=0; i < members; i++) { if (f[i] && f[i] != h) { DXDelete((Object)f[i]); f[i] = NULL; } } if (!h) { DXSetError(ERROR_MISSING_DATA, "no histogram data accumulated"); goto error; } /* if no counts in any bin, no median */ if (totalcnt == 0) goto done; /* else, set it */ if (hp->bins == 0) { INTERNALERROR; goto error; } /* XXX start median recalc */ /* the code in this area is trying to set the binsize, and not * dividing by zero during the process */ /* XXX md */ if (fabs(hp->max - hp->min) < FUZZ) binsize = 1.0; else binsize = (hp->max - hp->min) / hp->bins; if (binsize == 0) { INTERNALERROR; goto error; } /* XXX md - have to do this per axis */ gather = 0; for (i=0; ibins; i++) { if ((gather+accum[i]) >= (totalcnt/2.)) break; gather += accum[i]; } if (accum[i] == 0) { INTERNALERROR; goto error; } else median = (i*binsize)+hp->min + ((((totalcnt/2.) - gather) / accum[i]) * binsize); hp->median = median; /* XXX end median calc */ done: if (!DXEndField(h)) goto error; DXFree((Pointer)pinfo); DXFree((Pointer)f); return h; error: for (i=0; i < members; i++) { if (f && f[i] && f[i] != h) DXDelete((Object)f[i]); } DXDelete((Object)h); DXFree((Pointer)pinfo); DXFree((Pointer)f); return NULL; } static Field Array_Histogram(struct histinfo *hp) { return HistogramND(hp, 0); } static Field EmptyHistogram(struct histinfo *hp) { return HistogramND(hp, 1); } /* allocate and set up space for per/dim items */ static Error init_perdim(struct perdim **pdp, int dim) { /* allocate zero so uninitialized pointers are NULL. */ *pdp = (struct perdim *)DXAllocateZero(sizeof(struct perdim)); if (! *pdp) return ERROR; if (((*pdp)->binsplus1 = (int *)DXAllocate(dim * sizeof(int))) == NULL) goto error; if (((*pdp)->binsize = (float *)DXAllocate(dim * sizeof(float))) == NULL) goto error; if (((*pdp)->thisbin = (int *)DXAllocate(dim * sizeof(int))) == NULL) goto error; if (((*pdp)->dminlist = (float *)DXAllocate(dim * sizeof(float))) == NULL) goto error; if (((*pdp)->dmaxlist = (float *)DXAllocate(dim * sizeof(float))) == NULL) goto error; if (((*pdp)->slicecnt = (int **)DXAllocateZero(dim * sizeof(int *))) == NULL) goto error; if (((*pdp)->deltas = (float *)DXAllocate(dim * dim * sizeof(float))) == NULL) goto error; return OK; error: free_perdim(*pdp, dim); return ERROR; } static Error bins_perdim(struct perdim *pdp, int dim, int *binsize) { int i; for (i=0; islicecnt[i] = (int *)DXAllocateZero(binsize[i] * sizeof(int))) == NULL) return ERROR; return OK; } static void free_perdim(struct perdim *pdp, int dim) { int i; for (i=0; islicecnt[i] != NULL) DXFree((Pointer)pdp->slicecnt[i]); DXFree((Pointer)pdp->slicecnt); DXFree((Pointer)pdp->binsplus1); DXFree((Pointer)pdp->binsize); DXFree((Pointer)pdp->thisbin); DXFree((Pointer)pdp->dminlist); DXFree((Pointer)pdp->dmaxlist); DXFree((Pointer)pdp->deltas); DXFree((Pointer)pdp); } /* make sure all lists in histinfo struct are allocated to * the proper length, and if the previous dim was 1 and datadim > 1, * reallocate and replicate the existing values. if allocating * this for the first time, fill in the values which work if this * is going to be an empty field. */ static Error FullAllocHistinfo(struct histinfo *hp, int datadim, int setdef) { int dimchanged = 0; int i; if (hp->dim != datadim) { if (hp->dim == 1 && datadim > 1) { hp->dim = datadim; dimchanged++; } else { INTERNALERROR; return ERROR; } } #define DOALLOC(var, varlist, type) \ if (!varlist) { \ varlist = (type *)DXAllocate(hp->dim * sizeof(type)); \ if (!varlist) \ return ERROR; \ if (setdef) \ for (i=0; idim; i++) \ varlist[i] = var; \ } else if (dimchanged) { \ varlist = (type *)DXReAllocate(varlist, hp->dim * sizeof(type)); \ if (!varlist) \ return ERROR; \ for (i=0; idim; i++) \ varlist[i] = varlist[0]; \ } DOALLOC(hp->bins, hp->md_bins, int); DOALLOC(hp->min, hp->md_min, float); DOALLOC(hp->max, hp->md_max, float); DOALLOC(hp->inoutlo, hp->md_inoutlo, int); DOALLOC(hp->inouthi, hp->md_inouthi, int); DOALLOC(hp->median, hp->md_median, float); return OK; } static Field HistogramND(struct histinfo *hp, int empty) { Field f = NULL; Array a = NULL; Array sa = NULL; Array data; float zero = 0.0; int i, j, k; int bintotal, totalcnt; int defbins = DEFAULTBINS; Type type; Category cat; int rank; int shape[MAXRANK]; int vshape; int nitems; int exitearly; int index; #if OUTPUT_TYPE_FLOAT float *hd = NULL; /* histogram data, float */ #else int *hd = NULL; /* histogram data, was float before */ #endif float *hsum = NULL; /* sum of each histogram data bin */ Pointer od = NULL; /* original data */ int accum; int outofrange; struct perdim *pd; /* set up parms */ if (!init_perdim(&pd, hp->dim)) goto error; isempty: if (empty) { for (i=0; idim; i++) { if (!hp->binset) hp->md_bins[i] = defbins; if (!hp->minset) hp->md_min[i] = 0.0; if (!hp->maxset) hp->md_max[i] = 1.0; } } else { if (DXGetObjectClass(hp->o) != CLASS_ARRAY) { INTERNALERROR; goto error; } data = (Array)hp->o; if (!DXGetArrayInfo(data, &nitems, &type, &cat, &rank, shape)) goto error; vshape = GetVectorShape(rank, shape); if (cat != CATEGORY_REAL || vshape < 0) { DXSetError(ERROR_INVALID_DATA, "#10511", "data"); goto error; } /* if all user parms were scalar but the data turns out * to be vector, expand struct here. */ if (hp->dim == 1 && vshape != 1) { if (!FullAllocHistinfo(hp, vshape, 1)) goto error; free_perdim(pd, 1); if (!init_perdim(&pd, hp->dim)) goto error; } if (nitems <= 0) { empty++; if (type == TYPE_UBYTE) defbins = 256; goto isempty; } /* save these to allow early exit, and set parms if still not set */ if (hp->dim == 1) { if (!DXStatistics(hp->o, "data", pd->dminlist, pd->dmaxlist, NULL, NULL)) goto error; } else { if (!VectorStats((Object)data, 1, &hp->dim, 1, &pd->dminlist, &pd->dmaxlist)) goto error; } /* if the binsize hasn't been set yet, base it on the data type. for * byte data, default is data range (up to a max of 256); everything * else is 100 bins. */ for (i=0; idim; i++) { if (!hp->minset) hp->md_min[i] = pd->dminlist[i]; if (!hp->maxset) hp->md_max[i] = pd->dmaxlist[i]; if (!hp->binset) { hp->md_bins[i] = (type != TYPE_UBYTE) ? DEFAULTBINS : hp->md_max[i] - hp->md_min[i] + 1; } } if (!bins_perdim(pd, hp->dim, hp->md_bins)) goto error; } bintotal = 1; for (i=0; idim; i++) { pd->binsplus1[i] = hp->md_bins[i] + 1; bintotal *= hp->md_bins[i]; if (fabs(hp->md_max[i] - hp->md_min[i]) < FUZZ) pd->binsize[i] = 1.0; /* arbitrary number */ else pd->binsize[i] = (hp->md_max[i] - hp->md_min[i]) / hp->md_bins[i]; for (j=0; jdim; j++) pd->deltas[i * hp->dim + j] = (i == j) ? pd->binsize[i] : 0.0; } /* create output field. */ f = DXNewField(); if (!f) goto error; /* array = makegridpos(dims, counts, origins, deltas); */ if ((a = (Array)DXMakeGridPositionsV(hp->dim, pd->binsplus1, (Pointer)hp->md_min, (Pointer)pd->deltas)) == NULL) goto error; if (!DXSetComponentValue(f, "positions", (Object)a)) goto error; a = NULL; /* array = makegridconn(dims, counts); */ if ((a = (Array)DXMakeGridConnectionsV(hp->dim, pd->binsplus1)) == NULL) goto error; if (!DXSetComponentValue(f, "connections", (Object)a)) goto error; a = NULL; if (empty) { if (((a = (Array)DXNewConstantArray(bintotal, (Pointer)&zero, TYPE_FLOAT, CATEGORY_REAL, 0)) == NULL) || !DXSetStringAttribute((Object)a, "dep", "connections")) goto error; if (!DXSetComponentValue(f, "data", (Object)a)) goto error; if (!DXSetComponentValue(f, "sum", (Object)a)) goto error; a = NULL; goto done; } /* do the data bins here, plus median */ #if OUTPUT_TYPE_FLOAT if (((a = (Array)DXNewArray(TYPE_FLOAT, CATEGORY_REAL, 0)) == NULL) #else if (((a = (Array)DXNewArray(TYPE_INT, CATEGORY_REAL, 0)) == NULL) #endif || !DXSetStringAttribute((Object)a, "dep", "connections")) goto error; if (!DXAddArrayData(a, 0, bintotal, NULL)) goto error; #if OUTPUT_TYPE_FLOAT if ((hd = (float *)DXGetArrayData(a)) == NULL) #else if ((hd = (int *)DXGetArrayData(a)) == NULL) #endif goto error; for (i=0; idim)) == NULL) || !DXSetStringAttribute((Object)sa, "dep", "connections")) goto error; if (!DXAddArrayData(sa, 0, bintotal, NULL)) goto error; if ((hsum = (float *)DXGetArrayData(sa)) == NULL) goto error; for (i=0; idim; i++) hsum[i] = 0.0; if (!DXSetComponentValue(f, "sum", (Object)sa)) goto error; sa = NULL; /* things that allow the accumulation process to be shortcut: * all the data is out of range of the bins * all the data is the same value * * for an Nd histogram, this has to be a separate loop from * the one which is used if any of the data is in range, because * this loop is per dimension, and the other loop is per item. * i could make a flag which says skip this dimension, and use * it in the per/item loop, but it seems like a small gain for * a lot of code complexity. */ #define INDEX(hp, maxdim, thisindex) \ index = 0; \ for (j=1; jmd_bins[j]; \ index += thisindex; \ } #define ADDSUM(type) \ hsum[index*hp->dim + j] += (float)(((type *)od)[i*hp->dim + j]) /* XXX md - fill in slicecnt for median calc if field is partitioned. */ exitearly = 1; outofrange = 0; if (!hp->invalid) { for (i=0; idim; i++) { if (pd->dmaxlist[i] < hp->md_min[i]) { /* if all below min, put them in bin 0 if i/o flag set */ if (hp->md_inoutlo[i]) { INDEX(hp, i, 0); hd[index] = nitems; hp->md_median[i] = hp->md_min[i] + (0.5 * hp->md_bins[i]); totalcnt = nitems; pd->thisbin[i] = 0; } else outofrange++; } else if (pd->dminlist[i] > hp->md_max[i]) { /* if all above max, put in maxbin if i/o flag set */ if (hp->md_inouthi[i]) { INDEX(hp, i, hp->md_bins[i]-1); hd[index] = nitems; hp->md_median[i] = hp->md_max[i] - (0.5 * hp->md_bins[i]); totalcnt = nitems; pd->thisbin[i] = hp->md_bins[i]-1; } else outofrange++; } else if (pd->dminlist[i] == pd->dmaxlist[i]) { /* if all are the same, put them all into one bin enmass */ pd->thisbin[i] = (pd->dmaxlist[i] - hp->md_min[i]) / hp->md_bins[i]; INDEX(hp, i, pd->thisbin[i]); hd[index] = nitems; hp->md_median[i] = pd->dminlist[i]; totalcnt = nitems; } else { /* can't use this code, drop down and loop per item */ exitearly = 0; break; } } /* if all the data is in one bin or out of range, compute sun * and then we're done. */ if (exitearly) { if (outofrange == hp->dim) goto done; else if (outofrange == 0) { index = pd->thisbin[0]; for (k=1; kdim; k++) { index *= hp->md_bins[k]; index += pd->thisbin[k]; } od = DXGetArrayData(data); for (i=0; idim; j++) ADDSUM(ubyte); break; case TYPE_BYTE: for (j=0; jdim; j++) ADDSUM(byte); break; case TYPE_USHORT: for (j=0; jdim; j++) ADDSUM(ushort); break; case TYPE_SHORT: for (j=0; jdim; j++) ADDSUM(short); break; case TYPE_UINT: for (j=0; jdim; j++) ADDSUM(uint); break; case TYPE_INT: for (j=0; jdim; j++) ADDSUM(int); break; case TYPE_FLOAT: for (j=0; jdim; j++) ADDSUM(float); break; case TYPE_DOUBLE: for (j=0; jdim; j++) ADDSUM(double); break; default: DXSetError(ERROR_INVALID_DATA, "#10320", "histogram data"); goto error; } } goto done; } /* else fall through into the following code */ } } #define THISVAL(type) (((type *)od)[i*hp->dim + j]) #define THISBIN(type) ((THISVAL(type) - hp->md_min[j]) / pd->binsize[j]) totalcnt = 0; od = DXGetArrayData(data); for (i=0; iinvalid && DXIsElementInvalidSequential(hp->invalid, i)) continue; for (j=0; jdim; j++) { /* accumulate. if bin value is outside data range, check inout * flag. if not set, ignore those points. if set, add points * to the first or last bin. keep track of total pointcount * for later computation of median. */ switch(type) { case TYPE_UBYTE: pd->thisbin[j] = THISBIN(ubyte); break; case TYPE_BYTE: pd->thisbin[j] = THISBIN(byte); break; case TYPE_USHORT: pd->thisbin[j] = THISBIN(ushort); break; case TYPE_SHORT: pd->thisbin[j] = THISBIN(short); break; case TYPE_UINT: pd->thisbin[j] = THISBIN(uint); break; case TYPE_INT: pd->thisbin[j] = THISBIN(int); break; case TYPE_FLOAT: pd->thisbin[j] = THISBIN(float); break; case TYPE_DOUBLE: pd->thisbin[j] = THISBIN(double); break; default: DXSetError(ERROR_INVALID_DATA, "#10320", "histogram data"); goto error; } /* three cases: below, above, or within range */ if (pd->thisbin[j] < 0) { if (hp->md_inoutlo[j]) { pd->thisbin[j] = 0; } else break; } else if (pd->thisbin[j] >= hp->md_bins[j]) { if (hp->md_inouthi[j]) { pd->thisbin[j] = hp->md_bins[j]-1; } else break; } } /* did we exit early because at least one axis was out of range? */ if (j != hp->dim) continue; /* keep track of total counts per slice for each axis; this * is needed for median calc later. */ for (k=0; kdim; k++) pd->slicecnt[k][pd->thisbin[k]]++; index = pd->thisbin[0]; for (k=1; kdim; k++) { index *= hp->md_bins[k]; index += pd->thisbin[k]; } hd[index]++; totalcnt++; /* if the data is a scalar great, but otherwise it has to be * per component of the vector, which is more than a simple sum. * the sum has to be indexed one more level, shape[0]. */ switch(type) { case TYPE_UBYTE: for (j=0; jdim; j++) ADDSUM(ubyte); break; case TYPE_BYTE: for (j=0; jdim; j++) ADDSUM(byte); break; case TYPE_USHORT: for (j=0; jdim; j++) ADDSUM(ushort); break; case TYPE_SHORT: for (j=0; jdim; j++) ADDSUM(short); break; case TYPE_UINT: for (j=0; jdim; j++) ADDSUM(uint); break; case TYPE_INT: for (j=0; jdim; j++) ADDSUM(int); break; case TYPE_FLOAT: for (j=0; jdim; j++) ADDSUM(float); break; case TYPE_DOUBLE: for (j=0; jdim; j++) ADDSUM(double); break; default: DXSetError(ERROR_INVALID_DATA, "#10320", "histogram data"); goto error; } } /* compute median. * find bin which spans center value & interpolate */ for (i=0; idim; i++) { accum = 0; if (totalcnt == 0) { hp->md_median[i] = 0.0; continue; } for (j=0; jmd_bins[i]; j++) { if (accum + pd->slicecnt[i][j] >= totalcnt/2.) break; accum += pd->slicecnt[i][j]; } if (pd->slicecnt[i][j] == 0) { INTERNALERROR; goto error; } else hp->md_median[i] = (j*pd->binsize[i]) + hp->md_min[i] + (((totalcnt/2.) - accum) / pd->slicecnt[i][j] * pd->binsize[i]); } done: if (!DXEndField(f)) goto error; hp->didwork += totalcnt; free_perdim(pd, hp->dim); return f; error: DXDelete((Object)f); DXDelete((Object)a); DXDelete((Object)sa); free_perdim(pd, hp->dim); return NULL; } static Error HasInvalid(Field f, char *component, ICH *ih) { Error rc = OK; char *dep = NULL; char *invalid = NULL; *ih = NULL; if (!DXGetStringAttribute(DXGetComponentValue(f, component), "dep", &dep)) return OK; #define INVLEN 10 /* strlen("invalid "); */ if (!(invalid = (char *)DXAllocate(strlen(dep) + INVLEN))) return ERROR; strcpy(invalid, "invalid "); strcat(invalid, dep); /* if no component, return NULL */ if (!DXGetComponentValue(f, invalid)) goto done; *ih = DXCreateInvalidComponentHandle((Object)f, NULL, dep); if (!*ih) { rc = ERROR; goto done; } done: DXFree((Pointer)invalid); return rc; } static Error VectorLength(Object o, int *len, int onlyone) { int i, nitems, value = 0; int rank, shape[100]; Type t; Category c; if (DXGetObjectClass(o) != CLASS_ARRAY) return ERROR; if (!DXGetArrayInfo((Array)o, &nitems, &t, &c, &rank, shape)) return ERROR; if (onlyone && nitems > 1) return ERROR; if (c != CATEGORY_REAL) return ERROR; /* see comment below about what this does */ value = GetVectorShape(rank, shape); if (value < 0) return ERROR; if (len) *len = value; return OK; } /* if rank == 0, scalar. vector length = 1. * if rank == 1, vector length = shape[0]. * if rank > 1, error unless only one of the ranks has * a value > 1. * * returns -1 if error. */ static int GetVectorShape(int rank, int *shape) { int i; int value = -1; switch (rank) { case 0: value = 1; break; case 1: value = shape[0]; break; default: for (i=0; i 1) { if (value == -1) value = shape[i]; else return -1; } if (value == -1) value = 1; break; } return value; } static Error VectorExtract(Object o, Type expected_t, int len, Pointer d_ptr) { Array a = NULL; int i; int datasize; Pointer s_ptr; if (!(a = DXArrayConvert((Array)o, expected_t, CATEGORY_REAL, 1, len))) return ERROR; datasize = DXTypeSize(expected_t); s_ptr = DXGetArrayData(a); if (!s_ptr) goto error; memcpy(d_ptr, s_ptr, datasize * len); DXDelete((Object)a); return OK; error: DXDelete((Object)a); return ERROR; } static Error VectorAllocAndExtract(Object o, Type expected_t, int len, Pointer *d_ptr) { Array a = NULL; int i; int datasize; Pointer s_ptr; if (!(a = DXArrayConvert((Array)o, expected_t, CATEGORY_REAL, 1, len))) return ERROR; datasize = DXTypeSize(expected_t); s_ptr = DXGetArrayData(a); if (!s_ptr) goto error; if ((*d_ptr = DXAllocate(datasize * len)) == NULL) goto error; memcpy(*d_ptr, s_ptr, datasize * len); DXDelete((Object)a); return OK; error: DXDelete((Object)a); return ERROR; } /* currently unused */ static Error ScalarAlloc(Pointer inval, Type t, int len, Pointer *retval) { int i; int datasize; Pointer p; datasize = DXTypeSize(t); if ((*retval = DXAllocate(datasize * len)) == NULL) return ERROR; p = *retval; for (i=0; io) -- that is the input object which * we read but did not allocate and so should not free. * also, do not delete the contents of hp->p nor hp->invalid. * they aren't copied by the copy routine are are freed elsewhere. */ if (hp->md_bins) DXFree((Pointer)hp->md_bins); if (hp->md_min) DXFree((Pointer)hp->md_min); if (hp->md_max) DXFree((Pointer)hp->md_max); if (hp->md_inoutlo) DXFree((Pointer)hp->md_inoutlo); if (hp->md_inouthi) DXFree((Pointer)hp->md_inouthi); if (hp->md_median) DXFree((Pointer)hp->md_median); } static Error CopyHistinfo(struct histinfo *src, struct histinfo *dst) { /* first do a straight copy of the struct. * then zero all the pointers so they don't point into the src lists. * then allocate new space and copy any lists. * (if you don't zero the pointers first, you can't do meaningful * cleanup on error.) */ memcpy((char *)dst, (char *)src, sizeof(struct histinfo)); dst->md_bins = NULL; dst->md_min = NULL; dst->md_max = NULL; dst->md_inoutlo = NULL; dst->md_inouthi = NULL; dst->md_median = NULL; if (src->md_bins) { dst->md_bins = (int *)DXAllocate(src->dim * sizeof(int)); if (!dst->md_bins) goto error; memcpy(dst->md_bins, src->md_bins, src->dim * sizeof(int)); } if (src->md_min) { dst->md_min = (float *)DXAllocate(src->dim * sizeof(float)); if (!dst->md_min) goto error; memcpy(dst->md_min, src->md_min, src->dim * sizeof(float)); } if (src->md_max) { dst->md_max = (float *)DXAllocate(src->dim * sizeof(float)); if (!dst->md_max) goto error; memcpy(dst->md_max, src->md_max, src->dim * sizeof(float)); } if (src->md_inoutlo) { dst->md_inoutlo = (int *)DXAllocate(src->dim * sizeof(int)); if (!dst->md_inoutlo) goto error; memcpy(dst->md_inoutlo, src->md_inoutlo, src->dim * sizeof(int)); } if (src->md_inouthi) { dst->md_inouthi = (int *)DXAllocate(src->dim * sizeof(int)); if (!dst->md_inouthi) goto error; memcpy(dst->md_inouthi, src->md_inouthi, src->dim * sizeof(int)); } if (src->md_median) { dst->md_median = (float *)DXAllocate(src->dim * sizeof(float)); if (!dst->md_median) goto error; memcpy(dst->md_median, src->md_median, src->dim * sizeof(float)); } return OK; error: FreeHistinfo(dst); return ERROR; } static Error ObjectVectorStats(Object o, int *vlen, float **minlist, float **maxlist) { int i; Object subo, old, cpy; switch (DXGetObjectClass(o)) { case CLASS_ARRAY: return VectorStats(o, 0, vlen, 1, minlist, maxlist); case CLASS_FIELD: return VectorStats(o, 0, vlen, 1, minlist, maxlist); case CLASS_GROUP: switch (DXGetGroupClass((Group)o)) { case CLASS_MULTIGRID: case CLASS_SERIES: case CLASS_COMPOSITEFIELD: /* i'm just going for min/max on my stats here, * i can do that w/o worrying about boundary pts. */ subo = DXGetEnumeratedMember((Group)o, 0, NULL); if (VectorStats(subo, 0, vlen, 1, minlist, maxlist) == ERROR) return ERROR; for (i=1; subo = DXGetEnumeratedMember((Group)o, i, NULL); i++) { if (VectorStats(subo, 1, vlen, 0, minlist, maxlist) == ERROR) return ERROR; } return OK; default: return ERROR; } case CLASS_SCREEN: if (!DXGetScreenInfo((Screen)o, &subo, NULL, NULL)) return ERROR; return ObjectVectorStats(subo, vlen, minlist, maxlist); case CLASS_XFORM: if (!DXGetXformInfo((Xform)o, &subo, NULL)) return ERROR; return ObjectVectorStats(subo, vlen, minlist, maxlist); case CLASS_CLIPPED: if (!DXGetClippedInfo((Clipped)o, &subo, NULL)) return ERROR; return ObjectVectorStats(subo, vlen, minlist, maxlist); default: DXSetError(ERROR_BAD_PARAMETER, "#10190", "input"); return NULL; } /* not reached */ } /* this routine is getting more complicated than it was originally * intended to be. the first object can either be a simple field or * an array. if a field, the routine operates on the data component. * if the knownlength flag is set, the *vlen must match the * array or it is an error. if the flag isn't set, return the vector length * in *vlen. if initialize is set, the min/max lists are initialized to * max/min floats, otherwise they must contain previous results and the * stats from this field are combined with them. if one of minlist/maxlist * comes in as NULL, don't accumulate that value. one must be non-null. * if either of *minlist or *maxlist are NULL, allocate space for stats. * if *minlist or *maxlist are already allocated, use the space which is there. */ static Error VectorStats(Object o, int knownlength, int *vlen, int initialize, float **minlist, float **maxlist) { Array a = NULL, na = NULL; int i, j; int datasize; int nitems; int domin = 1, domax = 1; int allocmin = 0, allocmax = 0; float *d_ptr; if (DXGetObjectClass(o) == CLASS_FIELD) { a = (Array)DXGetComponentValue((Field)o, "data"); if (!a) goto error; } else if (DXGetObjectClass(o) == CLASS_ARRAY) { a = (Array)o; } else { INTERNALERROR; goto error; } if (!knownlength && !VectorLength((Object)a, vlen, 0)) goto error; if (minlist == NULL) domin = 0; if (maxlist == NULL) domax = 0; if (!domin && !domax) goto error; if (domin && (*minlist == NULL)) { if ((*minlist = (float *)DXAllocate(sizeof(float) * (*vlen))) == NULL) goto error; allocmin++; } if (domax && (*maxlist == NULL)) { if ((*maxlist = (float *)DXAllocate(sizeof(float) * (*vlen))) == NULL) goto error; allocmax++; } if (!(na = DXArrayConvert((Array)a, TYPE_FLOAT, CATEGORY_REAL, 1, *vlen))) return ERROR; d_ptr = (float *)DXGetArrayData(na); if (!d_ptr) goto error; if (!DXGetArrayInfo(na, &nitems, NULL, NULL, NULL, NULL)) goto error; if (initialize) for (i=0; i<*vlen; i++) { if (domin) (*minlist)[i] = DXD_MAX_FLOAT; if (domax) (*maxlist)[i] = DXD_MIN_FLOAT; } for (i = 0; i < nitems; i++) for (j = 0; j < *vlen; j++) { if (domin && (d_ptr[i * (*vlen) + j] < (*minlist)[j])) (*minlist)[j] = d_ptr[i * (*vlen) + j]; if (domax && (d_ptr[i * (*vlen) + j] > (*maxlist)[j])) (*maxlist)[j] = d_ptr[i * (*vlen) + j]; } DXDelete((Object)na); return OK; error: for (i = 0; i < *vlen; i++) { if (domin) (*minlist)[i] = DXD_MIN_FLOAT; if (domax) (*maxlist)[i] = DXD_MAX_FLOAT; } if (allocmin) DXFree((Pointer)minlist); if (allocmax) DXFree((Pointer)maxlist); DXDelete((Object)na); return ERROR; } #if 0 /* old routine */ static Field Array_Histogram1(struct histinfo *hp) { Field f = NULL; Array ao = (Array)hp->o; Array ap = NULL; /* positions, connections, data, median */ Array ac = NULL; Array ad = NULL; float *hd = NULL, *od = NULL; /* histogram data, original data */ int osbins; /* oversample the bins? */ int i, thisbin, nitems; int totalcnt; float accum; /* this could be an integer */ float binsize, median = 0; float binmin, binmax; float dmin, dmax; Point origin, delta; Type type; Category cat; int rank, shape[MAXRANK]; if (!DXGetArrayInfo(ao, &nitems, &type, &cat, &rank, shape)) goto error; if (cat != CATEGORY_REAL || (rank != 0 && (rank != 1 || shape[0] != 1))) { DXSetError(ERROR_INVALID_DATA, "#10080", "data"); goto error; } if (nitems <= 0) goto error; /* save these to allow early exit, and set parms if still not set */ if (!DXStatistics((Object)ao, "data", &dmin, &dmax, NULL, NULL)) goto error; /* set bin range, and number of bins */ binmin = hp->minset ? hp->min : dmin; binmax = hp->maxset ? hp->max : dmax; if (binmin > binmax) { DXSetError(ERROR_BAD_PARAMETER, "#10185", "min", "max"); goto error; } /* if the binsize hasn't been set yet, base it on the data type. * byte data default is data range (max 256); everything else is 100 bins. */ if (!hp->binset) { if ((type != TYPE_UBYTE) || (rank != 0 && (rank != 1 || shape[0] != 1))) { hp->bins = DEFAULTBINS; hp->binset = BYOBJECT; } else { hp->bins = binmax - binmin + 1; hp->binset = BYOBJECT; } } osbins = hp->bins; /* could be multiplier here to make smaller bins */ if (osbins <= 0) { DXSetError(ERROR_INVALID_DATA, "#10020", "number of bins"); goto error; } hd = (float *)DXAllocateLocalZero(osbins * sizeof(float)); if (!hd) goto error; if (fabs(binmax - binmin) < FUZZ) binsize = 1.0; /* arbitrary number */ else binsize = (binmax - binmin) / osbins; if (binsize <= 0) goto error; origin.x = binmin; origin.y = origin.z = 0.0; delta.x = binsize; delta.y = delta.z = 0.0; /* things that allow the accumulation process to be shortcut: * all the data is out of range of the bins * all the data is the same value */ totalcnt = 0; if (dmax < binmin) { if (hp->inoutlo) { hd[0] = nitems; totalcnt = nitems; median = binmin + (0.5 * binsize); } } else if (dmin > binmax) { if (hp->inouthi) { hd[osbins-1] = nitems; totalcnt = nitems; median = binmax - (0.5 * binsize); } } else if (dmin == dmax) { thisbin = (dmax - binmin) / binsize; if (thisbin >= 0 && thisbin < osbins) { hd[thisbin] = nitems; totalcnt = nitems; median = dmin; } else if (thisbin < 0 && hp->inoutlo) { hd[0] = nitems; totalcnt = nitems; median = binmin + (0.5 * binsize); } else if (thisbin >= osbins && hp->inouthi) { hd[osbins-1] = nitems; totalcnt = nitems; median = binmax - (0.5 - binsize); } } else { od = (float *)DXGetArrayData(ao); if (!od) goto error; /* accumulate. if bin value is outside data range, check inout * flag. if not set, ignore those points. if set, add points * to the first or last bin. keep track of total pointcount * for later computation of median. */ totalcnt = 0; #define THISBIN1(type) (((type *)od)[i] - binmin) / binsize for (i=0; iinvalid && DXIsElementInvalidSequential(hp->invalid, i)) continue; switch(type) { case TYPE_UBYTE: thisbin = THISBIN1(ubyte); break; case TYPE_BYTE: thisbin = THISBIN1(byte); break; case TYPE_USHORT: thisbin = THISBIN1(ushort); break; case TYPE_SHORT: thisbin = THISBIN1(short); break; case TYPE_UINT: thisbin = THISBIN1(uint); break; case TYPE_INT: thisbin = THISBIN1(int); break; case TYPE_FLOAT: thisbin = THISBIN1(float); break; case TYPE_DOUBLE: thisbin = THISBIN1(double); break; default: DXSetError(ERROR_INVALID_DATA, "#10320", "histogram data"); goto error; } /* three cases: below, above, or within range */ if (thisbin < 0) { if (hp->inoutlo) { hd[0]++; totalcnt++; } /* else ignore */ } else if (thisbin >= osbins) { if (hp->inouthi) { hd[osbins-1]++; totalcnt++; } /* else ignore */ } else { hd[thisbin]++; totalcnt++; } } /* compute median. * find bin which spans center value & interpolate */ if (totalcnt != 0) { accum = 0; for (i=0; i= totalcnt/2.) break; accum += hd[i]; } if (hd[i] == 0) { INTERNALERROR; goto error; } else median = (i*binsize)+binmin + ((((totalcnt/2.) - accum) / hd[i]) * binsize); } } /* create output field */ f = DXNewField(); if (!f) goto error; if (((ad = DXNewArray(TYPE_FLOAT, CATEGORY_REAL, 0)) == NULL) || !DXAddArrayData(ad, 0, osbins, (Pointer)hd) || !DXSetStringAttribute((Object)ad, "dep", "connections")) goto error; if (totalcnt > 0) hp->median = median; if (((ac = (Array)DXNewPathArray(osbins+1)) == NULL) || !DXSetStringAttribute((Object)ac, "ref", "positions") || !DXSetStringAttribute((Object)ac, "element type", "lines")) goto error; if ((ap = (Array)DXNewRegularArray(TYPE_FLOAT, 1, osbins+1, (Pointer)&origin, (Pointer)&delta)) == NULL) goto error; /* reset each array ptr, so it won't get deleted twice for errors * which happen after it becomes a member of the field. */ if (!DXSetComponentValue(f, "data", (Object)ad)) goto error; ad = NULL; if (!DXSetComponentValue(f, "connections", (Object)ac)) goto error; ac = NULL; if (!DXSetComponentValue(f, "positions", (Object)ap)) goto error; ap = NULL; if (!DXEndField(f)) goto error; DXFree((Pointer)hd); DXFreeInvalidComponentHandle(hp->invalid); hp->didwork += totalcnt; return f; error: DXDelete((Object)f); DXDelete((Object)ac); DXDelete((Object)ap); DXDelete((Object)ad); DXFree((Pointer)hd); DXFreeInvalidComponentHandle(hp->invalid); return NULL; } #endif /* old routine, unused now */