/* internal ana subroutines and functions */ /* contains changes needed for 64 bit machines */ /*this is file fun3.c */ #include #include #include #include #include "ana_structures.h" #include "defs.h" /* for SGI only (?) */ #if __sgi #include #include #endif extern struct sym_desc sym[]; extern struct sym_list *subr_sym_list[]; extern int temp_base; extern int edb_context; extern int ana_float(); extern char *find_sym_name(); extern int ana_type_size[]; extern int vfix_top, num_ana_subr, next_user_subr_num; extern int lastmin_sym, lastmax_sym, max_temp; extern float float_arg(); extern double double_arg(); extern int get_p_c(), get_p_c_r(); union types_ptr { byte *b; short *w; int *l; float *f; double *d;}; int scalemax = 255, scalemin= 0, lastmaxloc, lastminloc; int maxhistsize = 500000, histmin, histmax; float trim_min=.001, trim_max=.999; union scalar lastmin, lastmax; int callback_count = 0, fade_xsize, fade_ysize, fade_dim[2]; short *fade1, *fade2; /* used for fade intermediates */ unsigned int bitmask[] = { 0x00, 0x01, 0x03, 0x07, 0x0f, 0x1f, 0x3f, 0x7f, 0xff, 0x1ff,0x3ff,0x7ff,0xfff,0x1fff,0x3fff,0x7fff,0xffff, 0x0001ffff,0x0003ffff,0x0007ffff,0x000fffff,0x001fffff,0x003fffff, 0x007fffff,0x00ffffff,0x01ffffff,0x03ffffff,0x07ffffff,0x0fffffff, 0x1fffffff,0x3fffffff,0x7fffffff,0xffffffff}; static unsigned int bits[16]={1,2,4,8,16,32,64,128,256,512,1024,2048,4096, 8192,16384,32768}; int fgccdshiftDebug; /*------------------------------------------------------------------------ */ int ana_extract_bits_f(narg,ps) /* get a bit field from data */ /* call is xint = extract_bits(array, start, width) */ /* if width is < 0, then abs(width) is used for width and the sign is extended */ int narg, ps[]; { int result_sym, value; if (extract_bits(narg, ps, &value) != 1) return -1; result_sym = scalar_scratch(1); sym[result_sym].spec.scalar.w = (short) value; return result_sym; } /*------------------------------------------------------------------------ */ int ana_extract_bits(int narg, int ps[]) /* get a bit field from data */ /* call is extract_bits, result, array, start, width */ /* if width is < 0, then abs(width) is used for width and the sign is extended */ { int result_sym, value; short v; result_sym = ps[0]; if (extract_bits(narg-1, &ps[1], &value) != 1) return -1; v = value; redef_scalar(result_sym, 1, &v); return 1; } /*------------------------------------------------------------------------ */ int extract_bits(int narg, int ps[],int *value) /* internal routine */ { int n, start, width, iq, j, sign_flag, type; short *q; iq = ps[0]; n = get_p_c(&iq, &q); if (n <= 0) return -1; type = sym[iq].type; n = n * ana_type_size[type]; /* now the # of bytes */ if (int_arg_stat(ps[1], &start) != 1) return -1; if (int_arg_stat(ps[2], &width) != 1) return -1; /* check for sign extension */ if (width < 0) { width = - width; sign_flag = 1; } else sign_flag = 0; /* start is a bit count, get the I*2 index */ j = start/16; iq = (int) *(q+j); j = start % 16; /* our shift */ iq = (iq >> j) & bitmask[width]; if (sign_flag) { /* messy, but probably not used too often */ if (iq & bits[width-1]) { /* the top bit of the bit field is set, so need to extend it */ iq = iq | ( 0xffffffff - bitmask[width]); }} *value = iq; return 1; } /*------------------------------------------------------------------------ */ int ana_hkextract(int narg, int ps[]) /* get a bit field from HK data using EGSE info*/ /* call is xint = hkextract(array, startbyte, startbit, sizebit, mode), return value in xint is I*4 */ /* the modes supported are: 0 SB 1 UB 3 IU1 5 IU2 7 UL1 */ { int result_sym, value, bigendian, eflag, margin; int startbyte, startbit, sizebit, mode, iq, n, type; char b1, b2, *q; { int one = 1; bigendian = (*(char *)&one == 0); } iq = ps[0]; n = get_p_c(&iq, &q); if (n <= 0) return -1; type = sym[iq].type; n = n * ana_type_size[type]; /* now the # of bytes */ if (int_arg_stat(ps[1], &startbyte) != 1) return -1; if (int_arg_stat(ps[2], &startbit) != 1) return -1; if (int_arg_stat(ps[3], &sizebit) != 1) return -1; if (int_arg_stat(ps[4], &mode) != 1) return -1; value = eflag = 0; margin = 1; if (mode > 1) { if (mode > 5) margin = 4; else margin = 2; } if (startbyte > (n - margin)) { fprintf(stderr, "HKEXTRACT - start byte out of range, %d >= %d\n",startbyte, n); return -1; } b1 = *(q+startbyte); switch (mode) { case 0: case 1: { int shift; value = (int) b1; /* should work for both endians */ if (sizebit > 8 || sizebit <= 0) { eflag = 2; break; } shift = 8 - sizebit - startbit; if (shift < 0 || shift > 7) { eflag = 1; break; } value = value >> shift; value = value & bitmask[sizebit]; /* and for case 0 mode, we have to do a sign extension */ if (mode == 0 && (value & bits[sizebit-1])) { /* the top bit of the bit field is set, so need to extend it */ value = value | ( 0xffffffff - bitmask[sizebit]); } break; } case 2: case 3: /* the IU1 mode */ { int jq, shift; iq = (int) b1 & 0xff; /* this starts as bits 8-15 */ jq = (int) *(q+startbyte+1) & 0xff; /* this starts as bits 0-7 */ value = (iq << 8) + jq; shift = 16 - sizebit - startbit; if (shift < 0 || shift > 15) { eflag = 1; break; } if (shift) value = value >> shift; value = value & bitmask[sizebit]; /* and for case 2 mode, we have to do a sign extension */ if (mode == 2 && (value & bits[sizebit-1])) { /* the top bit of the bit field is set, so need to extend it */ value = value | ( 0xffffffff - bitmask[sizebit]); } break; } case 4: case 5: /* the IU2 mode, swap first */ { int jq, shift; iq = (int) b1 & 0xff; /* this starts as bits 0-7 */ jq = (int) *(q+startbyte+1) & 0xff; /* this starts as bits 8-15 */ value = (jq << 8) + iq; shift = 16 - sizebit - startbit; if (shift < 0 || shift > 15) { eflag = 1; break; } if (shift) value = value >> shift; value = value & bitmask[sizebit]; /* and for case 4 mode, we have to do a sign extension */ if (mode == 2 && (value & bits[sizebit-1])) { /* the top bit of the bit field is set, so need to extend it */ value = value | ( 0xffffffff - bitmask[sizebit]); } break; } case 7: /* do in an endian independent way (but takes longer) */ { int i1, i2, i3, i4, shift; i1 = (int) b1 & 0xff; /* this starts as bits 24-31 */ i2 = (int) *(q+startbyte+1) & 0xff; /* this starts as bits 16-23 */ i3 = (int) *(q+startbyte+2) & 0xff; /* this starts as bits 8-15 */ i4 = (int) *(q+startbyte+3) & 0xff; /* this starts as bits 1-7 */ value = (i1<< 24) + (i2<< 16) + (i3<< 8) +i4; shift = 32 - sizebit - startbit; if (shift < 0 || shift > 31) { eflag = 1; break; } if (shift) value = value >> shift; value = value & bitmask[sizebit]; break; } case 6: case 8: default: /* these are unsupported modes */ eflag = 99; break; } if (eflag) { fprintf(stderr, "HKEXTRACT - "); switch (eflag) { case 1: fprintf(stderr, "start bit %d and width %d bad for mode %d\n", startbit, sizebit, mode); break; case 2: fprintf(stderr, "width %d bad for mode %d\n", sizebit, mode); break; case 99: fprintf(stderr, "unsupported mode %d\n", mode); break; } return -1; } result_sym = scalar_scratch(2); sym[result_sym].spec.scalar.l = value; return result_sym; } /*------------------------------------------------------------------------- */ int ana_strtol(int narg, int ps[]) /* convert string to integer with optional base */ { char *p1; int base; int result_sym; if ( sym[ ps[0] ].class != 2 ) return execute_error(70); p1 = (char *) sym[ps[0] ].spec.array.ptr; result_sym = scalar_scratch(2); base = 10; if (narg > 1) { if (int_arg_stat(ps[1], &base) != 1) return -1; } sym[result_sym].spec.scalar.l = strtoul( p1, (char **)NULL, base); return result_sym; } /*------------------------------------------------------------------------- */ int ana_fgccdshift(int narg, int ps[]) /* specialized for processing FG images from Solar B (aka Hinode). Does an inplace pixel shift based on the parameter list usually taken from FITS keywords. This is a subr. The call is: fgccdshift, x, roistop, fg_ix1, fg_iy1, sum, bin, genid */ { /* this doesn't yet handle the half pixel shift for 2x2's */ int iq, roistop, fg_ix1, fg_iy1, ix0, iy0, sum, bin, genid; int rhs, lhs, ix1, iy1, nx, ny, extras, nxp, lhsmid, rhsmid; short *ptr, *pin, *pout, *pbase; int nym, nxm, cstride, nq, nsides, nhalf, partial; struct ahead *h; iq = ps[0]; /* only I*2 input accepted */ if ( sym[iq].class != 4 || sym[iq].type != 1) return execute_error(111); /* get the size of the image(s) */ h = (struct ahead *) sym[iq].spec.array.ptr; if ( h->ndim < 2 ) return execute_error(101); /* note - we should modify to do all images in astack when ndim = 3 */ nx = h->dims[0]; ny = h->dims[1]; ptr = (short *) ( (char *) h + sizeof( struct ahead ) ); pbase = ptr; /* rest of the list is presently a set of scalars: roistop, ix0, iy0, sum, bin */ if (int_arg_stat(ps[1], &roistop) != 1) return -1; if (int_arg_stat(ps[2], &fg_ix1) != 1) return -1; if (int_arg_stat(ps[3], &fg_iy1) != 1) return -1; if (int_arg_stat(ps[4], &sum) != 1) return -1; if (int_arg_stat(ps[5], &bin) != 1) return -1; if (int_arg_stat(ps[6], &genid) != 1) return -1; /* genid not used yet */ if (bin != 1) { printf("FGCCDSHIFT doesn't like bin > 1 yet\n"); return -1; } rhs = lhs = lhsmid = rhsmid = 0; /* partial = 0 is the full camera readout case */ if (roistop == 1025) partial = 0; else partial = 1; /* note that ix0 and iy0 are assumed inverted from FITS keywords (see solarb.ana), this puts them in correct direction for the FITS images */ ix0 = 4095 - fg_ix1; iy0 = 2047 - fg_iy1; switch (sum) { case 1: ix1 = ix0 + nx - 1; iy1 = iy0 + ny - 1; extras = 1; nhalf = 2048; /* for the middle of the CCD */ /* but if not full height, no rollover */ if (iy1 != 2047) extras = 0; /* note that this is equivalent to fg_ix0 = 0 */ break; case 2: ix1 = ix0 + 2*nx -1; /* note that ix0 and ix1 are in unsummed CCD pixels */ iy1 = iy0 + 2*ny -1; extras = 3; /* preliminary, 3 is in double pixels, I think the value is actually 5 singles */ nhalf = 1024; /* for the middle of the CCD */ /* but if not full height, no rollover */ if (iy1 != 2047) extras = 0; break; case 4: ix1 = ix0 + 4*nx -1; iy1 = iy0 + 4*ny -1; nhalf = 512; /* for the middle of the CCD */ printf("FGCCDSHIFT doesn't handle sum = 4 yet\n"); return -1; /* but if not full height, no rollover */ if (iy1 != 511) extras = 0; break; default: if (fgccdshiftDebug) printf("FGCCDSHIFT passed an illegal sum = %d\n", sum); return -1; break; } if (fgccdshiftDebug) printf("initial ix0, ix1, iy0, iy1 = %d, %d, %d, %d\n", ix0, ix1, iy0, iy1); /* adjust ix0 and ix1 for summed pixels */ if (sum > 1) { ix0 = ix0/sum; ix1 = (ix1+sum-1)/sum -1; iy0 = ix0/sum; iy1 = (iy1+sum-1)/sum -1; } if (fgccdshiftDebug) printf("ix0, ix1 = %d, %d\n", ix0, ix1); /* do the pixel shift to account for extra pixels read in by the camera electronics */ /* check if we have a rhs and/or a lhs and the middle column situations */ if (ix0 < nhalf) { lhs = 1; if (ix1 >= (nhalf-1)) lhsmid =1 ; } if (ix1 >= nhalf) { rhs = 1; if (ix0 <= nhalf) rhsmid =1 ; } if (fgccdshiftDebug) printf("partial, lhs, rhs, extras = %d, %d, %d, %d\n", partial, lhs, rhs, extras); if (fgccdshiftDebug) printf("lhsmid, rhsmid = %d, %d\n", lhsmid, rhsmid); nsides = lhs + rhs; /* we don't need to do anything shift-wise if a full read (not partial) and the top is out of the fov. The image location center is changed however. But that has to be handled elsewhere. The top edge is outside the fov if extras = 0. */ //if (partial || extras) { /* always run through the loop even if no change since we are planning to do dark and flat removal in this loop eventually */ while (nsides--) { /* clumsy but allows common code to be used for the 2 sides, do lhs first and then zero lhs for second pass for rhs (if any) */ if (lhs) { /* we have at least one N/S column on lhs */ nxp = nhalf - ix0; if (nxp > nx) nxp = nx; /* we have a split array to shift in the original CCD sense, this is different for the lhs and rhs, note that the array in memory is nx (not nxp) but we process just the nxp leftmost x values. */ pout = pbase + nx * ny - nx; pin = pout - nx * extras; if (fgccdshiftDebug) { printf("lhs\npout-pbase = %d\n", pout-pbase); printf("pin -pbase = %d\n", pin -pbase); } cstride = nx * ny + 1; nxm = nxp - 1; /* here we use nxp */ /* partials drop the right most column of the lhs if it is the middle column of the CCD, this has 2 effects on the shifts */ if (partial && lhsmid) { pin += 1; nxm--; } if (fgccdshiftDebug) printf("lhs nxm = %d\n", nxm); lhs = 0; } else if (rhs) { pout = pbase + nx * ny - 1; pin = pout - nx * extras; cstride = nx * ny - 1; nxm = ix1 - nhalf; /* columns (minus 1) in rhs */ /* partials drop the left most column of the rhs if it is the middle column of the CCD, this has 2 effects on the shifts */ if (partial && rhsmid) { pin -= 1; nxm--; } if (fgccdshiftDebug) { printf("rhs\npout-pbase = %d\n", pout-pbase); printf("pin -pbase = %d\n", pin -pbase); printf("rhs nxm = %d\n", nxm); } rhs = 0; } if (nxm >0 ) while (nxm--) { nym = ny - extras; /* do the column except for the wrap around */ while (nym--) { *pout = *pin; pout -= nx; pin -= nx; } pin += cstride; nq = extras; while (nq--) { *pout = *pin; pout -= nx; pin -= nx; } /* this should leave pout at the bottom of the column, so we need to add cstride to it now */ pout += cstride; } /* and the last column to finish, like the others but "extras" at the end are zeroed */ nym = ny - extras; if (fgccdshiftDebug) { printf("nym = %d\n", nym); printf("pout-pbase = %d\n", pout-pbase); printf("pin -pbase = %d\n", pin -pbase); } while (nym--) { *pout = *pin; pout -= nx; pin -= nx; } if (extras) { nq = extras; while (nq--) { *pout = 0; pout -= nx; } } } //} /* end of if (partial || extras) conditional */ /* check if the 2047 or 2048 column is in the image */ if (partial && (rhsmid || lhsmid)) { short *pa, *pb; int ix; float wa, wb, zq; if (fgccdshiftDebug) printf("middle column processing\n"); pa = 0; pb = 0; wa = wb = 0.33333333; /* usually we get both in an image but lots of special cases possible */ if (lhsmid) { ix = nhalf - ix0 - 1; pout = ptr + ix; /* this is the middle column */ if (ix > 0) { pa = pb = pout - 1; /* this is column to the left of it, also set pb in case */ } else { /* check if good rhs available */ if ((ix+2) >= nx) { /* if neither, we have a problem */ printf("image too narrow? nx = %d, ix0, ix1 = %d, %d\n", nx, ix0, ix1); return -1; } pa = pb = pout + 2; /* only one, so set both pa and pb */ } /* find the rhs column for the lhs middle */ ix = ix + 2; if (ix < nx) pb = ptr + ix; /* otherwise we have pb = pa */ nym = ny; /* a full column here */ if (fgccdshiftDebug) printf("lhs pa-ptr, pb-ptr = %d, %d\n", pa-ptr, pb-ptr); wa = 0.6666666667; wb = 0.33333333; while (nym--) { zq = wa * (float) *pa + wb * (float) *pb; *pout = (short) zq; pout += nx; pa += nx; pb += nx; } } if (rhsmid) { ix = nhalf - ix0; pout = ptr + ix; /* this is the middle column */ if (ix1 > nhalf) { pa = pb = pout + 1; /* this is column to the right of it, also set pb in case */ } else { /* we should have a good lhs available in this case or we would have bailed out in the lhsmid section above */ if ((ix-2) < 0) { /* if neither, we have a problem */ printf("image too narrow? nx = %d, ix0, ix1 = %d, %d\n", nx, ix0, ix1); return -1; } pa = pb = pout - 2; /* only one, so set both pa and pb */ } /* find the rhs column for the lhs middle */ ix = ix - 2; if (ix > 0) pa = ptr + ix; /* otherwise we have pb = pa */ nym = ny; /* a full column here */ if (fgccdshiftDebug) printf("rhs pa-ptr, pb-ptr = %d, %d\n", pa-ptr, pb-ptr); wb = 0.6666666667; wa = 0.33333333; while (nym--) { zq = wa * (float) *pa + wb * (float) *pb; *pout = (short) zq; pout += nx; pa += nx; pb += nx; } } } return 1; } /*------------------------------------------------------------------------- */ int ana_ccdflatmtrim(int narg, int ps[]) /* specialized flat field function */ /* subroutine version, call is CCDFLATMTRIM, XIN, DARK, GAIN, IMAGE where XIN is the image to be processed, it is returned modified according to XIN = ( XIN - DARK)*GAIN and has the same data type and dimension as the input 9/11/89 - but we assume that the input is in raw format if the dimensions fit certain critera; i.e., if the data is 1-D and of certain lengths IMAGE is loaded with a trim byte image suitable for display 9/17/93 - designed to work with both I*2 and I*1 arrays, the special raw formats are different for the I*2 case and the I*1 case */ { int i, j, iq, n, result_sym, nf, rawflag, nx, ny, ns, nd, dim[2], *hp; unsigned int kq; int i2flag, delta, jq; short *p1, *p2, *p3, *pr, *psave; byte *p1b, *p2b, *prb, *psaveb, jqb; struct ahead *h; byte *pb; int htop, hbot; float xq, yq; extern int scrat[]; /* scratch storage, also used in anatest.c */ /* must be 8192 I*4 for purposes here */ //extern int histmin, histmax; /* 4 arguments, the first is an I*2 image */ iq = ps[0]; if ( sym[iq].class != 4 ) return execute_error(66); switch ( sym[iq].type ) { case 0: i2flag = 0; break; case 1: i2flag = 1; break; default: return execute_error(119); } h = (struct ahead *) sym[iq].spec.array.ptr; if (i2flag == 1) { p1 = (short *) ((char *)h + sizeof(struct ahead)); pr = psave = p1; /* check if raw, must be 1-D and one of certain sizes */ rawflag = 0; ny = 0; if ( h->ndim == 1) { /* maybe */ nx = h->dims[0]; delta = 26; switch (nx ) { case 276506: h->ndim = 2; h->dims[0] = h->dims[1] = 512; rawflag = 512; p1 += 1076; break; /* 512 x 512 case */ case 1076250: h->ndim = 2; h->dims[0] = h->dims[1] = 1024; rawflag = 1024; p1 += 1076; break; /* 1024 x 1024 case */ case 263194: h->ndim = 2; h->dims[0] = h->dims[1] = 512; rawflag = 0; p1 += 1050; break; /* m1c case */ } } } else { p1b = (byte *) ((char *)h + sizeof(struct ahead)); prb = psaveb = p1b; rawflag = 0; ny = 0; if ( h->ndim == 2) { /* maybe */ nx = h->dims[0]; ny = h->dims[1]; if ( nx == 1360 && ny == 1036) { /* yes */ p1b += 8200; h->dims[0] = 1280; h->dims[1] = 1024; rawflag = 1280; delta = 80; } } } /* whatever, nx and ny should now agree with array dimensions */ nx = h->dims[0]; ny = h->dims[1]; ns = nx * ny; /* dark counts array */ iq = ps[1]; if ( sym[iq].class != 4 ) return execute_error(66); h = (struct ahead *) sym[iq].spec.array.ptr; if ( i2flag ) { if ( sym[iq].type != 1 ) return execute_error(111); p2 = (short *) ((char *)h + sizeof(struct ahead)); } else { if ( sym[iq].type != 0 ) return execute_error(118); p2b = (byte *) ((char *)h + sizeof(struct ahead)); } n = 1; /* # of elements must agree with data array (as possibly modified) don't really mind if dimension details are right */ nd = h->ndim; for (j=0;jdims[j]; if ( n != ns ) { printf("data and dark arrays are incompatiable\n"); return -1; } /* gain array, must be I*2 always, a scaled type of data */ iq = ps[2]; if ( sym[iq].class != 4 ) return execute_error(66); if ( sym[iq].type != 1 ) return execute_error(111); h = (struct ahead *) sym[iq].spec.array.ptr; p3 = (short *) ((char *)h + sizeof(struct ahead)); n = 1; /* # of elements must array with data array (as possibly modified) don't really mind if dimension details are right */ nd = h->ndim; for (j=0;jdims[j]; if ( n != ns ) { printf("data and gain arrays are incompatiable\n"); return -1; } /* set up the output byte array for display */ iq = ps[3]; dim[0] = nx; dim[1] = ny; redef_array(iq, 0, 2, dim); /* assumed to be 2-D */ h = (struct ahead *) sym[iq].spec.array.ptr; pb = (byte *) ((char *)h + sizeof(struct ahead)); /* use scrat for the histogram */ hp = scrat; for (j=0;j<4096;j++) *hp++ = 0; /*zero histogram */ hp = scrat; if (rawflag == 0) { nx = ns/4; ny = 1; } else { nx = nx/4; } /* top of flat field loop, also does histogram */ if (i2flag) { /* note that for non-raw, the array is contiguous and we set ny=1 with nx the element count/4, the 26 added doesn't happen unless it is raw which has the normal ny, nx and therefore loops over ny */ while (ny--) { n = nx; while (n--) { jq = *p1++ - *p2++; if (jq < 0) jq = 0; iq = ((int) *p3++ * (int) jq) >> 14; *pr++ = iq; /* histogram entry, every fourth one */ *(hp + (iq & 0xfff)) += 1; /* now 3 more pixels */ jq = *p1++ - *p2++; if (jq < 0) jq = 0; iq = ((int) *p3++ * (int) jq) >> 14; *pr++ = iq; jq = *p1++ - *p2++; if (jq < 0) jq = 0; iq = ((int) *p3++ * (int) jq) >> 14; *pr++ = iq; jq = *p1++ - *p2++; if (jq < 0) jq = 0; iq = ((int) *p3++ * (int) jq) >> 14; *pr++ = iq; } p1 += delta; /* effective only for full frame raws */ } } else { /* note that for non-raw, the array is contiguous and we set ny=1 with nx the element count/4, the 26 added doesn't happen unless it is raw which has the normal ny, nx and therefore loops over ny */ while (ny--) { n = nx; while (n--) { jq = (int) *p1b++ - (int) *p2b++; if (jq < 0) jq = 0; kq = ((unsigned int) *p3++ * (unsigned int) jq) >> 14; kq = MIN(kq, 255); *prb++ = (byte) kq; /* histogram entry, every fourth one */ *(hp + (kq & 0xff)) += 1; /* now 3 more pixels */ jq = (int) *p1b++ - (int) *p2b++; if (jq < 0) jq = 0; kq = ((unsigned int) *p3++ * (unsigned int) jq) >> 14; kq = MIN(kq, 255); *prb++ = (byte) kq; jq = (int) *p1b++ - (int) *p2b++; if (jq < 0) jq = 0; kq = ((unsigned int) *p3++ * (unsigned int) jq) >> 14; kq = MIN(kq, 255); *prb++ = (byte) kq; jq = (int) *p1b++ - (int) *p2b++; if (jq < 0) jq = 0; kq = ((unsigned int) *p3++ * (unsigned int) jq) >> 14; kq = MIN(kq, 255); *prb++ = (byte) kq; } p1b += delta; /* effective only for full frame raws */ } } /* done with flat field and we have a partial histogram (every fourth) */ /* convert the fractional histograms targets into # of points */ hbot = (int) ( (float) (ns/4) * trim_min); htop = (int) ( (float) (ns/4) * trim_max); /* printf("hbot, htop = %d %d\n",hbot, htop); */ /* scan histogram */ n = 4096; iq = 0; hp = scrat; while (n--) { iq += *hp++; if (iq >= hbot) break; } hbot = hp - scrat; while (n--) { iq += *hp++; if (iq >= htop) break; } htop = hp - scrat; /* printf("bottom and top indices = %d %d\n",hbot, htop); */ /* set up the lookup table, also use scratch memory */ hp = scrat; if (htop <= hbot) { /* not a good omen, make everybody gray */ for (j=0;j<8192;j++) *hp++ = 127; } else { for (j=0;j0) { *p++ = *(ind + *q++); *p++ = *(ind + *q++); *p++ = *(ind + *q++); *p++ = *(ind + *q++); } /* the rest */ while (nr-- >0) { *p++ = *(ind + *q++); } } return 1; } /*------------------------------------------------------------------------- */ int ana_factors(int narg, int ps[]) /* factors function */ /* return a vector of factors of an integer */ { static int primes[]= {2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,57,61 ,67,71,73,79,83,89,97,101,103,107,109,113,127,131,137,139,149,151,157 ,159,161,163,167,173,179,181,191,193,197,199,211,223,227,229,239,241 ,251,257,263,269,271,277,281,283,293,307,311,313,317,331,337,347,349 ,353,359,367,373,379,383,389,397,401,409,419,421,431,433,439,443 ,449,457,461,463,467,479,487,491,499,503,509,521,523,541,547,557,563, 569,571,577,587,593,599,601,607,613,617,619,631,641,643,647,653,659, 661,673,677,683,691,701,709,719,727,733,739,743,751,757,761,769,773, 787,797,809,811,821,823,827,829,839,853,857,859,863,877,881,883,887,907,911, 919,929,937,941,947,953,967,971,977,983,991,997,1009, 1013,1019,1021,1031,1033}; extern int scrat[]; /* scratch storage, also used in anatest.c */ #define NPRIMES (sizeof( primes) / sizeof(int)) int iq, lim, i, nfac, *p, *q, result_sym, fac; struct ahead *h; iq = int_arg( ps[0]); /* need a scalar integer */ if ( iq == 0 ) return 4; /* a scalar zero if input is 0 */ if ( iq < 0 ) iq = -iq; /* lim = (int) sqrt( (float) iq) + 1; */ nfac = 0; p = scrat; i = 0; fac = primes[i]; while ( fac < iq) { while ( iq%fac == 0 ) { /* got a factor ? */ iq = iq/fac; *p++ = fac; nfac++; } i++; if (i >= NPRIMES) { printf( "FACTORS, exceeded prime table, only checks for prime factors le 599\n"); break; } lim = (int) sqrt( (float) iq) + 1; fac = primes[i]; if (fac > lim) break; } /* done, last factor may be in iq if it is not 1 */ if (iq > 1) { *p = iq; nfac++; } /* factors in scrat, nfac of them, copy into a long array */ result_sym = array_scratch(2, 1, &nfac); h = (struct ahead *) sym[result_sym].spec.array.ptr; q = (int *) ((char *)h + sizeof(struct ahead)); bcopy( (char *) scrat, (char *) q, nfac*sizeof(int)); return result_sym; } /*------------------------------------------------------------------------- */ int ana_execute_string(int narg, int ps[]) /* parse and execute a string */ { int iq, nsym; char *s; iq = ps[0]; if (sym[iq].class != 2) { return execute_error(70); } s = (char *) sym[iq].spec.array.ptr; nsym = ana_execute_symbol(s,0); if (nsym > 0) return execute(nsym); printf("EXECUTE: parsed something non-executable\n"); return -1; } /*------------------------------------------------------------------------- */ int ana_execute_symbol(s,context) /* parse a string and return nsym */ char *s; int context; { extern single_string_parse_flag, single_string_parse_nest, nest; extern byte *input_compile_string; /* used for compiling strings */ extern byte line2[]; extern struct sdesc cur_line; extern char *strsavsd(); struct sdesc save_cur_line; int nsym, iq, save_edb_context; char call_back_name[20]; /* check if something in the parser buffer and save if so */ /* we'll put it back at the end */ if (save_cur_line.n=cur_line.n) save_cur_line.p=(byte *) strsavsd(&cur_line); cur_line.n = 0; input_compile_string = (byte *) s; single_string_parse_flag = 1; single_string_parse_nest = nest; /* set to current nest level */ save_edb_context = edb_context; if (context) { /* the context arg is currently used to make edb's that have the context of code blocks and are therefore permanent, they have unique block names with lc letters (so user names can't match) */ /* this is not a real code block because the block structures won't know where the compiled code is, we call it directly by the nsym returned by parser, not by the block name */ sprintf(call_back_name,"callback%d", callback_count); callback_count++; if ( (iq = setup_block(call_back_name, call_back_name)) < 0 ) { printf("error while processing callback execution string/n"); } edb_context = iq + 1 + MAX_USER_SUBRS + MAX_USER_FUNCS; } nsym = parser(); edb_context = save_edb_context; /* printf("nsym = %d\n", nsym); */ single_string_parse_flag = 0; /* important to do this before an execute since the execute may parse something */ /* in case there was a partial line at the start, we copy the original line back into line2 (our preprocessed buffer) */ if (cur_line.n = save_cur_line.n) { strncpy((char *) line2, (char *) save_cur_line.p, save_cur_line.n); cur_line.p = line2; free(save_cur_line.p); } return nsym; } /*end of return */ /*------------------------------------------------------------------------- */ int ana_eval(int narg, int ps[]) /* eval function */ /* compiles and evaluates a string, interpreting it as code */ /* non-strings are just evaluated */ { extern byte line[],line2[]; extern struct sdesc cur_line, token; extern char *strsavsd(); /*extern int symbol_context;*/ struct sdesc save_cur_line; int iq; byte *p, *pq; iq = ps[0]; switch (sym[iq].class) { case 1: case 4: default: return iq; /* anything but a string is just passed */ case 2: /* compile the string as an expression, run through conversion table and put in line2 (after saving it in case someone was using it)*/ /*printf("eval: symbol_context = %d\n", symbol_context);*/ if (save_cur_line.n=cur_line.n) save_cur_line.p=(byte *) strsavsd(&cur_line); cur_line.n = 0; pq=line2; p = (byte *) sym[iq].spec.array.ptr; iq = linevert(p, pq); if (iq <= 0) { execute_error(108); iq=-1; } else { token.n = cur_line.n; token.p = cur_line.p; cur_line.n = 0; /* reset or parser gets confused */ if ( (iq = compile_tok()) < 0) { execute_error(108); iq=-1; } } /* because we don't want to keep the malloc'ed string around forever, we copy the original line back into line2 (our preprocessed buffer) */ if (cur_line.n = save_cur_line.n) { strncpy((char *) line2, (char *) save_cur_line.p, save_cur_line.n); cur_line.p = line2; free(save_cur_line.p); } if (iq >= 0) return eval(iq); else return iq; } } /*------------------------------------------------------------------------- */ int ana_tense(int narg, int ps[]) /* tense function */ /* splines under tension */ /* THE CALL IS YF=TENSE(X,Y,XF,[SIGMA],[DLHS,DRHS]) */ { int i, iq, len[3], n, result_sym, nf; double sigma, der_left, der_right, *st, *yf; union types_ptr p[3]; struct ahead *h; /* first 3 args must be 1-D vectors */ for (i=0;i<3;i++) { iq = ps[i]; if ( sym[iq].class != 4 ) return execute_error(66); iq = ana_double(1, &iq); if ( (len[i] = get_p_c(&iq, &p[i])) <= 0) return -1; } /* take smaller of X and Y size */ n = MIN( len[0], len[1]); sigma = -1.0; der_left = 0.0; der_right = 0.0; /* defaults */ if (narg > 3) sigma = -double_arg( ps[3] ); /* sign of sigma is a flag, + => that slopes (der's) are input */ if (narg > 4) { der_left = double_arg( ps[4] ); sigma = - sigma; } if (narg > 5) der_right = double_arg( ps[5] ); result_sym = array_clone(ps[2], 4); h = (struct ahead *) sym[result_sym].spec.array.ptr; yf = (double *) ((char *)h + sizeof(struct ahead)); nf = len[2]; /* scratch storage for 2 double arrays */ st = (double *) malloc( 16 * n ); /* printf("n = %d\n", n); */ #ifdef NeXT iq = curv1(&n, p[0].d, p[1].d, &der_left, &der_right, st, (st+n), &sigma, p[2].d, yf, &nf); #else iq = curv1_(&n, p[0].d, p[1].d, &der_left, &der_right, st, (st+n), &sigma, p[2].d, yf, &nf); #endif /* printf("returned iq = %d\n", iq); */ free( st); if ( iq == 1) return result_sym; else return -1; } /*------------------------------------------------------------------------- */ int ana_tense_curve(int narg, int ps[]) /* tense_curve function */ /* splines under tension */ /* THE CALL IS XY=TENSE_CURVE(X,Y,XS,[SIGMA],[DLHS,DRHS]) */ { /* this differs from plain TENSE in that it allows any open ended curve to be defined, it does not have to be a single valued function of x the specification of points along the curve is complicated because of this, the XS input is in units of the polygon length of the curve and can range from 0 to 1.0 the output is a 2-D array containing matching x and y coordinates on the curve */ int i, iq, len[3], n, result_sym, dim[2], nf; double sigma, der_left, der_right, *st, *yf, *xf; union types_ptr p[3]; struct ahead *h; /* first 3 args must be 1-D vectors */ for (i=0;i<3;i++) { iq = ps[i]; if ( sym[iq].class != 4 ) return execute_error(66); iq = ana_double(1, &iq); if ( (len[i] = get_p_c(&iq, &p[i])) <= 0) return -1; } /* take smaller of X and Y size */ n = MIN( len[0], len[1]); sigma = 1.0; der_left = 0.0; der_right = 0.0; /* defaults */ if (narg > 3) sigma = double_arg( ps[3] ); if (narg > 4) der_left = double_arg( ps[4] ); if (narg > 5) der_right = double_arg( ps[5] ); sigma = - sigma; dim[0] = nf = len[2]; dim[1] = 2; result_sym = array_scratch(4, 2, dim); h = (struct ahead *) sym[result_sym].spec.array.ptr; xf = (double *) ((char *)h + sizeof(struct ahead)); yf = xf + nf; /* scratch storage for 3 double arrays */ st = (void *) malloc( 24 * n ); #ifdef NeXT iq = kurv1( &n, p[0].d, p[1].d, &der_left, &der_right, st,(st+n),(st+n+n), &sigma, p[2].d, xf, yf, &nf); #else iq = kurv1_( &n, p[0].d, p[1].d, &der_left, &der_right, st,(st+n),(st+n+n), &sigma, p[2].d, xf, yf, &nf); #endif free( st); if ( iq == 1) return result_sym; else return -1; } /*------------------------------------------------------------------------- */ int ana_tense_loop(int narg, int ps[]) /* tense_loop function */ /* generates spline under tension for any closed (x,y) curve */ /* THE CALL IS XY=TENSE_LOOP(X,Y,XS,[SIGMA]) */ { /* this differs from TENSE_CURVE in that it allows a closed curve to be defined, it does not have to be a single valued function of x the specification of points along the curve is complicated because of this, the XS input is in units of the polygon length of the curve and can range from 0 to 1.0 the output is a 2-D array containing matching x and y coordinates on the curve */ int i, iq, len[3], n, result_sym, dim[2], nf; double sigma, der_left, der_right, *st, *yf, *xf; union types_ptr p[3]; struct ahead *h; /* first 3 args must be 1-D vectors */ for (i=0;i<3;i++) { iq = ps[i]; if ( sym[iq].class != 4 ) return execute_error(66); iq = ana_double(1, &iq); if ( (len[i] = get_p_c(&iq, &p[i])) <= 0) return -1; } /* take smaller of X and Y size */ n = MIN( len[0], len[1]); sigma = 1.0; /* defaults */ if (narg > 3) sigma = double_arg( ps[3] ); sigma = - sigma; dim[0] = nf = len[2]; dim[1] = 2; result_sym = array_scratch(4, 2, dim); h = (struct ahead *) sym[result_sym].spec.array.ptr; xf = (double *) ((char *)h + sizeof(struct ahead)); yf = xf + nf; /* scratch storage for 3 double arrays */ /* one of which is 2*nf long */ st = (void *) malloc( 32 * n ); #if NeXT iq = kurvp1( &n, p[0].d, p[1].d, st,(st+n),(st+n+n), &sigma, p[2].d, xf, yf, &nf); #else iq = kurvp1_( &n, p[0].d, p[1].d, st,(st+n),(st+n+n), &sigma, p[2].d, xf, yf, &nf); #endif free( st); if ( iq == 1) return result_sym; else return -1; } /*------------------------------------------------------------------------- */ int ana_histr(int narg, int ps[]) /* histr function */ /* running sum histogram, normalized to 1.0 */ { int iq, n, nd, j, type; float sum, fac; struct ahead *h; union types_ptr q1, q2; /* first get a normal histogram */ if ( (iq = ana_hist(narg,ps) ) <= 0 ) return iq; /* the returned symbol should be a long array, convert in place to a float, this depends on float and long (int) being 32 bits! */ /* check for some impossible errors */ if ( sym[iq].class != 4 ) return execute_error(3); type = sym[iq].type; if (type != 2 ) return execute_error(3); n = get_p_c(&iq, &q1); if (n <= 0) return -1; sym[iq].type = 3; /* change to float */ sum = 0.0; j = n; q2.f = q1.f; while (j--) { sum += (float) *q2.l; *q2.f++ = sum; } fac = 1.0 / sum; j = n; q2.f = q1.f; while (j--) { *q2.f = fac * *q2.f; q2.f++; } return iq; } /*------------------------------------------------------------------------- */ int ana_hist(int narg, int ps[]) /* histogram function (frequency distrib.) */ /* general histogram function */ { struct ahead *h; int iq, i, n, nd, j, range, type, result_sym; union types_ptr q1, q2; iq = ps[0]; /* arg must be an array */ if ( sym[iq].class != 4 ) return execute_error(66); type = sym[iq].type; n = get_p_c(&iq, &q1); if (n <= 0) return -1; /* always need the range */ minmax( q1.l, n, type); /* get long (int) versions of min and max */ ana_convert(&lastmin, type, 2); ana_convert(&lastmax, type, 2); /* create a long array for results */ histmin = lastmin.l; histmax = lastmax.l; /* make the min 0 if it is greater */ histmin = histmin > 0 ? 0 : histmin; //if (histmin < 0 ) {f // printf("WARNING - histogram result contains negative entries, use\n"); // printf("!histmin and !histmax to find range included\n"); } range = histmax - histmin +1; if (range > maxhistsize) { printf("range (%d) larger than !maxhistsize (%d)\n",range, maxhistsize); return -1; } result_sym = array_scratch(2, 1, &range); h = (struct ahead *) sym[result_sym].spec.array.ptr; q2.l = (int *) ((char *)h + sizeof(struct ahead)); ana_zero( 1, &result_sym); /* need to zero initially */ /* now accumulate the distribution */ switch (type) { case 0: while (n--) { i = *q1.b++ - histmin; *(q2.l + i) += 1; } break; case 1: while (n--) { i = *q1.w++ - histmin; *(q2.l + i) += 1; } break; case 2: while (n--) { i = *q1.l++ - histmin; *(q2.l + i) += 1; } break; case 3: while (n--) { i = *q1.f++ - histmin; *(q2.l + i) += 1; } break; case 4: while (n--) { i = *q1.d++ - histmin; *(q2.l + i) += 1; } break; } return result_sym; } /*------------------------------------------------------------------------- */ int ana_key_locations(int narg, int ps[]) /* locate integer keys */ /* for integer arrays only, finds unique values and their locations, it is similar in some ways to where_unique but also returns all the indices for each unique value (key) in a sorted result array */ /* call is key_locations, array, values, counts, indexsymarr array is input, others are output */ { struct ahead *h; int iq, i, n, nd, j, range, type, m, nuniq, k, mq; int *values, *counts, *p, *q, val_sym, cnt_sym, indexsymarr; int *countp, **indexptrs, **qp, *valueoffsets, *pv; union types_ptr q1, q2; struct sym_desc *pbase, *pdesc; iq = ps[0]; val_sym = ps[1]; cnt_sym = ps[2]; indexsymarr = ps[3]; /* arg must be an array */ if ( sym[iq].class != 4 ) return execute_error(66); type = sym[iq].type; if (type >2 ) { printf("non-integer array for key_locations\n"); return -1; } n = get_p_c(&iq, &q1); if (n <= 0) return -1; /* always need the range */ minmax( q1.l, n, type); /* get long (int) versions of min and max */ ana_convert(&lastmin, type, 2); ana_convert(&lastmax, type, 2); /* create a long array for results */ histmin = lastmin.l; histmax = lastmax.l; range = histmax - histmin +1; if (range > maxhistsize) { printf("range (%d) larger than !maxhistsize (%d) in key_locations\n",range, maxhistsize); return -1; } /* need some scratch to keep the temporary values and indices */ values = (int *) malloc(range * sizeof (int)); valueoffsets = (int *) malloc(range * sizeof (int)); bzero( (void *) values, range * sizeof (int)); bzero( (void *) valueoffsets, range * sizeof (int)); nuniq = 0; q2 = q1; /* save for later */ /* first a pass to get the count of each unique value and the # of uniques */ switch (type) { case 0: for (j=0;jclass = ARRAY; pdesc->type = LONG; pdesc->xx = - 1; /* this is not named */ mq = 4 * cq + sizeof(struct ahead); pdesc->spec.array.bstore = mq; if ( (pdesc->spec.array.ptr = (int *) malloc(mq) ) == NULL ) return ana_malloc_error(mq); h = (struct ahead *) pdesc->spec.array.ptr; h->ndim = 1; h->c1 = 0; h->c2 = 0; h->dims[0] = cq; h->facts = NULL; /*no known facts */ *qp++ = (int *) ((char *)pdesc->spec.array.ptr + sizeof(struct ahead)); *pv++ = i; i++; pdesc++; if (i > nuniq) { printf("internal error 1 in ana_key_locations\n"); return -1; } } else *pv++ = -1; /* we need a pv for each value within "range" */ } /* now we have to scan the original array again to collect all the locations and populate the index arrays */ q1 = q2; switch (type) { case 0: for (j=0;j2 ) { printf("non-integer array for where_unique\n"); return -1; } n = get_p_c(&iq, &q1); if (n <= 0) return -1; /* always need the range */ minmax( q1.l, n, type); /* get long (int) versions of min and max */ ana_convert(&lastmin, type, 2); ana_convert(&lastmax, type, 2); /* create a long array for results */ histmin = lastmin.l; histmax = lastmax.l; /* make the min 0 if it is greater */ histmin = histmin > 0 ? 0 : histmin; range = histmax - histmin +1; if (range > maxhistsize) { printf("range (%d) larger than !maxhistsize (%d) in where_unique\n",range, maxhistsize); return -1; } /* need some scratch to keep the temporary values and indices */ values = (int *) malloc(range * sizeof (int)); indices = (int *) malloc(range * sizeof (int)); bzero( (void *) values, range * sizeof (int)); switch (type) { case 0: for (j=0;j does not contain non-zero elements, then a scalar -1 is returned. alternate use: X=SIEVE(condition) yields the same as X=SIEVE(LONG(INDGEN(condition)),condition) */ { int iq, i, j, n, nd, nc, cnt, *p, type, typec, result_sym; struct ahead *h; union types_ptr q1, q2, q3; iq = ps[0]; if (narg == 2) { /* two-argument case */ /* first arg must be an array */ if ( sym[iq].class != 4 ) return execute_error(66); type = sym[iq].type; h = (struct ahead *) sym[iq].spec.array.ptr; q1.l = (int *) ((char *)h + sizeof(struct ahead)); nd = h->ndim; n = 1; for(j=0;jdims[j]; /* second argument is the conditional */ iq = ps[1]; } nc = get_p_c(&iq, &p); if (nc <= 0) return -1; typec = sym[iq].type; /* unlike the VMS code, we do a double pass, first to get the count */ /* use the min of size of array and conditional if 2 args */ if (narg == 1) {n = nc; type = LONG;} else n = n <= nc ? n : nc; nc = n; cnt = 0; q3.l = p; switch ( typec) { case 0: while (nc--) if (*q3.b++ != 0) cnt += 1; break; case 1: while (nc--) if (*q3.w++ != 0) cnt += 1; break; case 2: while (nc--) if (*q3.l++ != 0) cnt += 1; break; case 3: while (nc--) if (*q3.f++ != 0) cnt += 1; break; case 4: while (nc--) if (*q3.d++ != 0) cnt += 1; break; } if (cnt <= 0) return 4; /* this returns a scalar 0 */ q3.l = p; result_sym = array_scratch(type,1, &cnt); /*for the result */ h = (struct ahead *) sym[result_sym].spec.array.ptr; q2.l = (int *) ((char *)h + sizeof(struct ahead)); nc = n; if (narg == 1) { /* one-argument case LS 12may92 */ n = 0; /* use as indgen counter */ switch (typec) { case BYTE: while (nc--) {if (*q3.b++ != 0) *q2.l++ = n; n++;} break; case WORD: while (nc--) {if (*q3.w++ != 0) *q2.l++ = n; n++;} break; case LONG: while (nc--) {if (*q3.l++ != 0) *q2.l++ = n; n++;} break; case FLOAT: while (nc--) {if (*q3.f++ != 0) *q2.l++ = n; n++;} break; case DOUBLE: while (nc--) {if (*q3.d++ != 0) *q2.l++ = n; n++;} break; } } else { /* two-argument case */ /* this is a 25 case situation */ switch ( typec) { case 0: switch (type) { case 0: while (nc--) {if (*q3.b++ != 0) *q2.b++ = *q1.b; q1.b++;} break; case 1: while (nc--) {if (*q3.b++ != 0) *q2.w++ = *q1.w; q1.w++;} break; case 2: while (nc--) {if (*q3.b++ != 0) *q2.l++ = *q1.l; q1.l++;} break; case 3: while (nc--) {if (*q3.b++ != 0) *q2.f++ = *q1.f; q1.f++;} break; case 4: while (nc--) {if (*q3.b++ != 0) *q2.d++ = *q1.d; q1.d++;} break; } break; case 1: switch (type) { case 0: while (nc--) {if (*q3.w++ != 0) *q2.b++ = *q1.b; q1.b++;} break; case 1: while (nc--) {if (*q3.w++ != 0) *q2.w++ = *q1.w; q1.w++;} break; case 2: while (nc--) {if (*q3.w++ != 0) *q2.l++ = *q1.l; q1.l++;} break; case 3: while (nc--) {if (*q3.w++ != 0) *q2.f++ = *q1.f; q1.f++;} break; case 4: while (nc--) {if (*q3.w++ != 0) *q2.d++ = *q1.d; q1.d++;} break; } break; case 2: switch (type) { case 0: while (nc--) {if (*q3.l++ != 0) *q2.b++ = *q1.b; q1.b++;} break; case 1: while (nc--) {if (*q3.l++ != 0) *q2.w++ = *q1.w; q1.w++;} break; case 2: while (nc--) {if (*q3.l++ != 0) *q2.l++ = *q1.l; q1.l++;} break; case 3: while (nc--) {if (*q3.l++ != 0) *q2.f++ = *q1.f; q1.f++;} break; case 4: while (nc--) {if (*q3.l++ != 0) *q2.d++ = *q1.d; q1.d++;} break; } break; case 3: switch (type) { case 0: while (nc--) {if (*q3.f++ != 0) *q2.b++ = *q1.b; q1.b++;} break; case 1: while (nc--) {if (*q3.f++ != 0) *q2.w++ = *q1.w; q1.w++;} break; case 2: while (nc--) {if (*q3.f++ != 0) *q2.l++ = *q1.l; q1.l++;} break; case 3: while (nc--) {if (*q3.f++ != 0) *q2.f++ = *q1.f; q1.f++;} break; case 4: while (nc--) {if (*q3.f++ != 0) *q2.d++ = *q1.d; q1.d++;} break; } break; case 4: switch (type) { case 0: while (nc--) {if (*q3.d++ != 0) *q2.b++ = *q1.b; q1.b++;} break; case 1: while (nc--) {if (*q3.d++ != 0) *q2.w++ = *q1.w; q1.w++;} break; case 2: while (nc--) {if (*q3.d++ != 0) *q2.l++ = *q1.l; q1.l++;} break; case 3: while (nc--) {if (*q3.d++ != 0) *q2.f++ = *q1.f; q1.f++;} break; case 4: while (nc--) {if (*q3.d++ != 0) *q2.d++ = *q1.d; q1.d++;} break; } break; } } return result_sym; } /*------------------------------------------------------------------------- */ int ana_nwhere(int narg, int ps[]) /* nwhere or wherecount function */ /* based on sieve but returns just the number of incidents for which the condition is true, often this is all we want and the full sieve just wastes time */ { int iq, n, cnt, *p, type, result_sym; union types_ptr q1; iq = ps[0]; n = get_p_c(&iq, &q1); if (n <= 0) return -1; type = sym[iq].type; cnt = 0; switch (type) { case 0: while (n--) if (*q1.b++ != 0) cnt += 1; break; case 1: while (n--) if (*q1.w++ != 0) cnt += 1; break; case 2: while (n--) if (*q1.l++ != 0) cnt += 1; break; case 3: while (n--) if (*q1.f++ != 0) cnt += 1; break; case 4: while (n--) if (*q1.d++ != 0) cnt += 1; break; } result_sym = scalar_scratch(type); /*for the result */ sym[result_sym].spec.scalar.l = cnt; return result_sym; } /*------------------------------------------------------------------------- */ int ana_maxf(int narg, int ps[]) /* finds the max element in an array */ { return maxormin(narg, ps, 0); } /*------------------------------------------------------------------------- */ int maxormin(int narg, int ps[],int code) /* used by max, min, maxloc, minloc */ { struct ahead *h, *h2; int iq, type, j, result_sym, nd, n; union types_ptr q1; iq = ps[0]; switch (sym[iq].class) { default: return execute_error(32); case 8: iq = class8_to_1(iq); /*scalar ptr case */ case 1: /*scalar case */ /* if min or max return the symbols, if loc, then return a scalar 0 */ if (code <2 ) return iq; else return 4; case 4: /*array case */ type = sym[iq].type; h = (struct ahead *) sym[iq].spec.array.ptr; q1.l = (int *) ((char *)h + sizeof(struct ahead)); nd = h->ndim; n = 1; for(j=0;jdims[j]; minmax( q1.l, n, type); /* puts results in globals */ /* now it depends on code */ if (code < 2) { result_sym = scalar_scratch(type); switch (code) { case 0: bcopy( (char *) &lastmax.l, (char *) &sym[result_sym].spec.scalar.l, 8); break; case 1: bcopy( (char *) &lastmin.l, (char *) &sym[result_sym].spec.scalar.l, 8); break; } return result_sym; } /* a location needed */ result_sym = array_scratch(2, 1, &nd); h2 = (struct ahead *) sym[result_sym].spec.array.ptr; q1.l = (int *) ((char *)h2 + sizeof(struct ahead)); if (code == 2) iq = lastmaxloc; else iq = lastminloc; for (j=0;jdims[j]; iq = iq/h->dims[j]; } return result_sym; } } /*------------------------------------------------------------------------- */ int ana_minf(int narg, int ps[]) /* finds the min element in an array */ { return maxormin(narg, ps, 1); } /*------------------------------------------------------------------------- */ int ana_minloc(int narg, int ps[]) /* finds the min location in an array */ { return maxormin(narg, ps, 3); } /*------------------------------------------------------------------------- */ int ana_maxloc(int narg, int ps[]) /* finds the max location in an array */ { return maxormin(narg, ps, 2); } /*------------------------------------------------------------------------- */ int ana_scale(int narg, int ps[]) /* scale an array to the range (!scalemin,!scalemax) and put into a byte array !scalemax must be <= 255 */ { struct ahead *h; int iq, type, j, result_sym, nd, n; union scalar min, max; union types_ptr q1, q2; float fq, range, sf, x, qf; double dq, drange, sd, d, qd; iq = ps[0]; if ( sym[iq].class != 4 ) return execute_error(66); type = sym[iq].type; h = (struct ahead *) sym[iq].spec.array.ptr; q1.l = (int *) ((char *)h + sizeof(struct ahead)); nd = h->ndim; n = 1; for(j=0;jdims[j]; result_sym = array_clone(iq, 0); h = (struct ahead *) sym[result_sym].spec.array.ptr; q2.l = (int *) ((char *)h + sizeof(struct ahead)); /* if only one arg., then we scale between the min and max of array */ if (narg == 1) { minmax( q1.l, n, type); /* get the min and max */ simple_scale( q1.l, n, type, q2.b ); /* scale and put in output array */ return result_sym; /* done for this case */ } else { /* more than one arg., so we expect the min and max to be specified */ if (narg != 3 ) return execute_error(55); /* this is messier than the simple scale because of over/under flows */ if (type == 4) { /*array is double, get input as double */ min.d = double_arg(ps[1]); max.d = double_arg(ps[2]); drange = max.d - min.d; if (drange == 0.0) { neutral( q2.b, n); return result_sym; } sd =(double) scalemax; qd = (double) scalemin; dq = (sd - qd) / drange; } else { min.f = float_arg(ps[1]); max.f = float_arg(ps[2]); range = max.f - min.f; if (range == 0.0) { neutral( q2.b, n); return result_sym; } sf = (float) scalemax; qf= (float) scalemin; fq = (sf-qf)/ range; } switch (type) { case 0: while (n--) { x = fq * (*q1.b++ - min.f) + qf; x = x > qf ? x : qf; x = x < sf ? x : sf; *q2.b++ = (byte) x; } return result_sym; case 1: while (n--) { x = fq * (*q1.w++ - min.f) + qf; x = x > qf ? x : qf; x = x < sf ? x : sf; *q2.b++ = (byte) x; } return result_sym; case 2: while (n--) { x = fq * (*q1.l++ - min.f) + qf; x = x > qf ? x : qf; x = x < sf ? x : sf; *q2.b++ = (byte) x; } return result_sym; case 3: while (n--) { x = fq * (*q1.f++ - min.f) + qf; x = x > qf ? x : qf; x = x < sf ? x : sf; *q2.b++ = (byte) x; } return result_sym; case 4: while (n--) { d = dq * (*q1.d++ - min.d) + qd; d = d > qd ? d : qd; d = d < sd ? d : sd; *q2.b++ = (byte) d; } return result_sym; } } } /*------------------------------------------------------------------------- */ int minmax(int *p, int n, int type) /* finds the min and max (and their locations) for an array */ { union types_ptr q; union scalar min, max; /* 3/25/94, make maxloc and minloc pointers to work on 64 bit machines */ int nc; byte *minloc, *maxloc; q.l = p; nc = n - 1; minloc = maxloc = (byte *) p; switch (type) { case 0: min.b = max.b = *q.b++; while (nc--) { if (*q.b > max.b) { max.b = *q.b; maxloc = (byte *) q.b; } if (*q.b < min.b) { min.b = *q.b; minloc = (byte *) q.b; } q.b++; } lastmin.b = min.b; lastmax.b = max.b; sym[lastmin_sym].type = 0; sym[lastmax_sym].type = 0; break; case 1: min.w = max.w = *q.w++; while (nc--) { if (*q.w > max.w) { max.w = *q.w; maxloc = (byte *) q.w; } if (*q.w < min.w) { min.w = *q.w; minloc = (byte *) q.w; } q.w++; } lastmin.w = min.w; lastmax.w = max.w; sym[lastmin_sym].type = 1; sym[lastmax_sym].type = 1; break; case 2: min.l = max.l = *q.l++; while (nc--) { if (*q.l > max.l) { max.l = *q.l; maxloc = (byte *) q.l; } if (*q.l < min.l) { min.l = *q.l; minloc = (byte *) q.l; } q.l++; } lastmin.l = min.l; lastmax.l = max.l; sym[lastmin_sym].type = 2; sym[lastmax_sym].type = 2; break; case 3: min.f = max.f = *q.f++; while (nc--) { if (*q.f > max.f) { max.f = *q.f; maxloc = (byte *) q.f; } if (*q.f < min.f) { min.f = *q.f; minloc = (byte *) q.f; } q.f++; } lastmin.f = min.f; lastmax.f = max.f; sym[lastmin_sym].type = 3; sym[lastmax_sym].type = 3; break; case 4: min.d = max.d = *q.d++; while (nc--) { if (*q.d > max.d) { max.d = *q.d; maxloc = (byte *) q.d; } if (*q.d < min.d) { min.d = *q.d; minloc = (byte *) q.d; } q.d++; } lastmin.d = min.d; lastmax.d = max.d; sym[lastmin_sym].type = 4; sym[lastmax_sym].type = 4; break; } /* end of type switch */ lastmaxloc = (maxloc - (byte *) p) / ana_type_size[type]; lastminloc = (minloc - (byte *) p) / ana_type_size[type]; return 1; } /*------------------------------------------------------------------------- */ int simple_scale(int *p1, int n, int type, byte *p2) /* scale n elements starting at p1 into byte array at p2 */ /* uses the min and max stashed in lastmin and lastmax */ /* assumes that the input array values are all inbetween the min and max */ /* this means we don't need to check for over and under flow */ { union types_ptr q; int xq; float fq; double dq; q.l = p1; /* 12/18/97 - change the ranges to be the delta + 1, this gives a better mapping after truncation when range gt the scale difference, only for I*1 and I*2. The I*4 and F*4 and F*8 use rounding */ xq = (int) (scalemax - scalemin); switch (type) { case 0: { register int min, range, iscalemin; iscalemin = scalemin; min = lastmin.b; range = lastmax.b - min; if ( range == 0 ) return neutral(p2, n ); if (range > xq) { range++; xq++; } while (n--) *p2++ = (byte) (( xq * (*q.b++ - min) )/range + iscalemin); } return; case 1: { register int min, range; min = lastmin.w; range = lastmax.w - min; if ( range == 0 ) return neutral(p2, n ); if (range > xq) { range++; xq++; } while (n--) *p2++ = (byte) (( xq * (*q.w++ - min) )/ range)+scalemin; } return; case 2: /* to avoid possible range problems, calculation is done in fp */ { register int min, range; min = lastmin.l; range = lastmax.l - min; if ( range == 0 ) return neutral(p2, n ); fq = (float) xq/ (float) range; while (n--) *p2++ = (byte) (fq * (*q.l++ - min)+0.5) + scalemin; }return; case 3: { register float min, range; min = lastmin.f; range = lastmax.f - min; if ( range == 0 ) return neutral(p2, n ); fq = (float) xq/ range; while (n--) *p2++ = (byte) ( fq * (*q.f++ - min)+0.5 ) + scalemin; } return; case 4: { register double min, range; min = lastmin.d; range = lastmax.d - min; if ( range == 0 ) return neutral(p2, n ); dq = (double) xq/ range; while (n--) *p2++ = (byte) ( dq * (*q.d++ - min)+0.5 )+scalemin; } return 1; } } /*------------------------------------------------------------------------- */ int neutral(byte *p, int n) /* fills a byte array with (scalemax-scalemin)/2 */ { byte xq; xq = (scalemax-scalemin)/2; while (n--) *p++ = xq; return 1; } /*------------------------------------------------------------------------- */ int ana_fade_init(int narg, int ps[]) /* initializes a fade (or dissolve) between 2 byte arrays, the actual fade result is obtained with the subroutine fade call is: fade_init, x1, x2 */ { byte *p1, *p2; short *q1, *q2; int n1, n2, n, nd, i, iq; struct ahead *h; iq = ps[0]; if ( sym[iq].class != 4 ) return execute_error(66); if (sym[iq].type != 0) return execute_error(118); h = (struct ahead *) sym[iq].spec.array.ptr; if ( h->ndim != 2 ) return execute_error(101); p1 = (byte *) ((char *)h + sizeof(struct ahead)); fade_xsize = h->dims[0]; fade_ysize = h->dims[1]; fade_dim[0] = fade_xsize; fade_dim[1] = fade_ysize; iq = ps[1]; if ( sym[iq].class != 4 ) return execute_error(66); if (sym[iq].type != 0) return execute_error(118); h = (struct ahead *) sym[iq].spec.array.ptr; if ( h->ndim != 2 ) return execute_error(101); if (h->dims[0] != fade_xsize || h->dims[1] != fade_ysize) { printf("FADE_INIT: arguments must have the same dimensions\n"); } p2 = (byte *) ((char *)h + sizeof(struct ahead)); /* set up the I*2 arrays for quick fades */ if (fade1) { free(fade1); fade1 = NULL; } if (fade2) { free(fade2); fade2 = NULL; } /* since the inputs are byte arrays, n1 is the size in bytes, we need 2*n1 for the fade arrays */ n = 2 * fade_xsize * fade_ysize; fade1 = (short *) malloc(n); fade2 = (short *) malloc(n); n = fade_xsize * fade_ysize; q1 = fade1; q2 = fade2; while (n--) { *q1++ = ( (short) *p2) << 8 ; *q2++ = (( (short) *p1++) - ( (short) *p2++)); } /* while (n--) { *q1++ = ( (short) *p1++); *q2++ = ( (short) *p1++); } */ return 1; } /*------------------------------------------------------------------------- */ int ana_fade(int narg, int ps[]) /* does a fade (or dissolve) between 2 byte arrays, must be initialized with fade_init call is: fade, weight, result - where n is from 0 to 256 and represents the weight of the first array that was passed to fade_init result is a subroutine argument to avoid mallocs */ { byte *p1; short *q1, *q2, wt; int weight, n, iq; struct ahead *h; if (int_arg_stat(ps[0], &weight) != 1) return -1; wt = (short) weight; iq = ps[1]; if (redef_array(iq, 0, 2, fade_dim) != 1) return -1; h = (struct ahead *) sym[iq].spec.array.ptr; p1 = (byte *) ((char *)h + sizeof(struct ahead)); n = fade_xsize * fade_ysize; q1 = fade1; q2 = fade2; while (n--) { *p1++ = (byte) (( wt * *q2++ + *q1++) >> 8); } return 1; } /*------------------------------------------------------------------------- */ int ana_rscale(int narg, int ps[]) /* re-scale a byte array to a specified range from the array's min and max; i.e., scales min/max range to input range; input must be a byte array and is returned modified; this is a subroutine call is: rscale, array, new_min, new_max */ { extern int scrat[]; /* scratch storage */ int iq, type, i, n, nsave, *hist; byte *pq, *lut; float new_min, new_max, old_min, old_max; float fq; iq = ps[0]; if ( sym[iq].class != 4 ) return execute_error(66); type = sym[iq].type; if (type != 0) return execute_error(118); n = nsave = get_p_c(&iq, &pq); if (n <= 0) return -1; if (float_arg_stat(ps[1], &new_min) != 1) return -1; if (float_arg_stat(ps[2], &new_max) != 1) return -1; /* get the max/min via a histogram, use hist pointer and scrat */ hist = (int *) scrat; bzero((char *) hist, 256*sizeof(int)); { register byte *p = pq; int nd4=n/4, nr; nr = n -4*nd4; while (nd4-- >0) { (hist[*p++])++; (hist[*p++])++; (hist[*p++])++; (hist[*p++])++; } /* the rest */ while (nr-- >0) { (hist[*p++])++; } } n = nsave; /* the first non-zero entry in lut is the min */ for (i=0;i<255;i++) { if (hist[i]) break; } old_min = (float) i; /* highest non-zero is the max */ for (i=255;i>0;i--) { if (hist[i]) break; } old_max = (float) i; /* printf("max, min = %f, %f\n", old_max, old_min); */ fq = (new_max - new_min + 1)/( old_max - old_min + 1); /* 12/24/96 modify to use lookup table for conversion, should speed it up */ lut = (byte *) scrat; for (i=0;i<256;i++) *lut++ = (byte) ( ((float) i - old_min) * fq + new_min ); /* clamp just to be sure */ lut = (byte *) scrat; for (i=0;i<256;i++) { if (*lut < (byte) new_min) *lut = (byte) new_min; else if (*lut > (byte) new_max) *lut = (byte) new_max; lut++; } lut = (byte *) scrat; /* while (n--) { *pq = lut[ *pq]; pq++; } */ { register byte *p = pq; register byte *ind = lut; int nd4=n/4, nr; nr = n -4*nd4; while (nd4-- >0) { *p++ = *(ind + *p); *p++ = *(ind + *p); *p++ = *(ind + *p); *p++ = *(ind + *p); } /* the rest */ while (nr-- >0) { *p++ = *(ind + *p); } } return 1; } /*------------------------------------------------------------------------- */ int ana_circle_fit_lsq(int narg, int ps[]) /* fit a bunch of points to a circle circle_fit_lsq, x, y, w, xc, yc, r, chsq */ { float xc, yc, r, chi, *x, *y, *w; int n, nd, j, nx, iq, wf; /* get pointer to x and size, and float it */ iq = ps[0]; if ( sym[iq].class != 4 ) return execute_error(66); iq = ana_float(1, &iq); nx = get_p_c(&iq, &x); if (nx <= 0) return -1; /* similarly for y, check against x */ iq = ps[1]; if ( sym[iq].class != 4 ) return execute_error(66); iq = ana_float(1, &iq); n = get_p_c(&iq, &y); if (n <= 0) return -1; if (nx != n) {printf("CIRCLE_FIT_LSQ: input arrays must match\n"); return -1;} /* now the weight array */ iq = ps[2]; if ( sym[iq].class != 4 ) return execute_error(66); iq = ana_float(1, &iq); n = get_p_c(&iq, &w); if (n <= 0) return -1; if (nx != n) { printf("CIRCLE_FIT_LSQ: weight array must match (x,y)\n"); return -1; } /* plan to upgrade to allow a scalar in w which would imply all weights are unity */ wf = 0; circle_fit(x, y, w, n, wf, &xc, &yc, &r, &chi); if (set_scalar_value( ps[3], 3, &xc) != 1) return execute_error(13); if (set_scalar_value( ps[4], 3, &yc) != 1) return execute_error(13); if (set_scalar_value( ps[5], 3, &r) != 1) return execute_error(13); if (narg > 6) { if (set_scalar_value( ps[6], 3, &chi) != 1) return execute_error(13); } return 1; } /*------------------------------------------------------------------------- */ #define EPSILON 1.e-10 /*------------------------------------------------------------------------- */ int circle_fit(x, y, w, n, wf, xc_out, yc_out, r_out, chi_out) float *x, *y, *w, *xc_out, *yc_out, *r_out, *chi_out; int n, wf; { /* 4/11/96 - adapted from the C++ code in the Gism package which in turn was taken from something that was taken from PREPMORT file by CVERT on 1 Mar 1991 which is after the method of chernov and ososkov in computer physics communications 33 (1984) 329-333. */ /* return convention same as Gismo version: 0 OK 1 < 3 points (but not checked here) 2 no convergence, returns non-iterative solution 3 singular determinant 4 rsquared < 0 (imaginary radius!) */ double xxyy, sx,sy,sw,swinv, /*sums*/ sxx8,syy8,sxy8,sxxy8,syxy8,sxxyy8, /*more sums*/ xt,xtp,top,bot, x8,y8,w8,x08,y08,rsq,chisq8, f,g,h,p,q,t, /*see the paper*/ a0,b0,c0,d0,det, gamma, gamma0,g0inv; float *xp, *yp, *wp; int i, nptw, iflag = 0; int j, k, nn; /* To improve accuracy, we do all calculation in a shifted coordinate*/ /* system where the weighted sum of the xi and the yi are both zero */ /* this seems to work even in single precision, but I leave it in */ /* double just for clean living. */ /* Form the sums, which are used to shift the coordinate system */ sx = 0.0; sy = 0.0; sw = 0.0; nptw = 0; nn = n; xp = x; yp = y; wp = w; sx = sy = sw = 0.0; while (nn--) { w8 = *wp++; sw += w8; sx += w8 * *xp++; sy += w8 * *yp++; } swinv = 1./sw; sx = sx*swinv; sy = sy*swinv; /* normalize */ /* These are the sums needed in the calculation */ sxx8 = syy8 = sxy8 = sxxy8 = syxy8 = sxxyy8 = 0.0; nn = n; xp = x; yp = y; wp = w; while (nn--) { w8 = *wp++; x8 = *xp++ - sx; y8 = *yp++ - sy; xxyy = x8*x8 + y8*y8; sxx8 += x8*x8*w8; syy8 += y8*y8*w8; sxy8 += x8*y8*w8; sxxy8 += x8*xxyy*w8; syxy8 += y8*xxyy*w8; sxxyy8 += xxyy*xxyy*w8; } f = swinv*(3.*sxx8 + syy8); g = swinv*(sxx8 + 3.*syy8); h = 2.*swinv*sxy8; p = swinv*sxxy8; q = swinv*syxy8; t = swinv*sxxyy8; /* The initial guess gamma comes from the non-iterative fast 'fit' */ gamma0 = swinv*(sxx8 + syy8); g0inv = 1./gamma0; a0 = -4.; b0 = (f*g - t - h*h)*g0inv*g0inv; c0 = (t*(f + g) - 2.*(p*p + q*q))*g0inv*g0inv*g0inv; d0 = (t*(h*h - f*g) + 2.*(p*p*g + q*q*f) - 4.*p*q*h)*g0inv*g0inv*g0inv*g0inv; /* Here we get the root by newton's method */ xt = 1.; for(i= 0; i < 10; i++) {top = xt*(xt*(xt*(xt + a0) + b0) + c0) + d0; bot = xt*(xt*(4.*xt + 3.*a0) + 2.*b0) + c0; xtp = xt -top/bot; if(fabs(xtp-xt) < EPSILON) break; xt = xtp; } if (i >= 10) { iflag = 2; xtp = 1.0; } /* no convergence */ /* for the no convergence case, we use the original gamma0 above */ gamma = gamma0*xtp; det = h*h - (f-gamma)*(g-gamma); /* Check on singularity of determinant */ if(fabs(det) <= 1.0e-16) {return 3;} x08 = (q*h - p*(g-gamma))/det; y08 = (p*h - q*(f-gamma))/det; rsq = x08*x08 + y08*y08 + gamma; if (rsq < 0.0) { return 4; } /* a problem, r^2 < 0 */ *r_out = (float) sqrt (rsq); /* Shift the coordinate system back */ *xc_out = (float) (x08 + sx); *yc_out = (float) (y08 + sy); /* Calculate pseudo-chi-squared/dof */ if(n == 3) { *chi_out = 0.0; } else { chisq8 = 0.0; nn = n; wp = w; xp = x; yp = y; while (nn--) { w8 = *wp++; x8 = *xp++ - sx; y8 = *yp++ - sy; chisq8 += w8*((x08+sx-x8)*(x08+sx-x8) + (y08+sy-y8)*(y08+sy-y8) - rsq)* ((x08+sx-x8)*(x08+sx-x8) + (y08+sy-y8)*(y08+sy-y8) -rsq); } *chi_out = (float) (chisq8/(4.0*rsq)/(n-3)); } return iflag; } /*------------------------------------------------------------------------- */