******************************************************************************** * resample.c -- Routines to resample a statistical series objects * * User routines * Pro_ResampledStat() creates a resampled statistical series object. * Pro_ResampledSize() returns the size of the needed series buffer. * * Mark Showalter, PDS Ring-Moon Systems Node, March 1998 *******************************************************************************/ #include <stdio.h> #include <string.h> #include <math.h> #include "profile.h" #include "fortran.h" /******************** * Type definitions * ********************/ typedef struct ZPRO_RESAMPLED_STRUCT { PRO_OBJECT *stat, *oldpsf, *newpsf, *curve; RL_INT4 segment; RL_INT4 oldk1, oldk2, newk1, newk2; RL_FLT8 oldx1, olddx, newx1, newdx; RL_INT4 oldtype, newtype; RL_FLT8 oldhw, oldpar, newhw, newpar; RL_FLT8 expand0; RL_FLT8 segment_min, segment_max; } ZPRO_RESAMPLED; /******************************** * Internal function prototypes * ********************************/ static void ZPro_ResampledKRange RL_PROTO((RL_INT4 knew, RL_INT4 *k1, RL_INT4 *k2, RL_VOID *params)); static RL_FLT8 ZPro_ResampledWeights RL_PROTO((RL_INT4 knew, RL_INT4 k, RL_VOID *params)); static void ZPro_PrintResampled RL_PROTO((RL_VOID *params)); /* ******************************************************************************** * EXPORTED USER ROUTINES ******************************************************************************** *$ Component_name: * Pro_ResampledStat (resample.c) *$ Abstract: * This routine resamples a statistical series so that it is uniformly * spaced in a new variable and based on a new point spread function (PSF). *$ Keywords: * PROFILE, SERIES * C, PUBLIC, SUBROUTINE *$ Declarations: * PRO_OBJECT *Pro_ResampledStat(stat, oldpsf, newpsf, curve, segment, * x1, dx, k1, k2) * PRO_OBJECT *stat, *oldpsf, *newpsf, *curve; * RL_INT4 segment, k1, k2; * RL_FLT8 x1, dx; *$ Inputs: * stat pointer to a statistical series to resample. * oldpsf pointer to the PSF of the original series, possibly * unnormalized, in units of series samples. * newpsf pointer to the new PSF desired, possibly unnormalized, * in units of new series samples. * curve pointer to the curve that evaluates new coordinate from * old. * segment segment index of the curve to use for the resampling. * x1, dx new x coordinate starting value and interval. * k1, k2 new range of indices. *$ Outputs: * none *$ Returns: * pointer to a new weighted statistical series; NULL on non-fatal error. *$ Detailed_description: * This routine resamples a statistical series so that it is uniformly * spaced in a new variable and based on a new point spread function (PSF). * * The returned statistical series is uniformly sampled in a new * coordinate, such that the give curve function converts from sampling * parameter of the original series to that of the new. For example, if * you have a ring occultation profile sampled uniformly in time plus a * curve that converts from time to ring-plane intercept radius, you can * quickly generate a revised profile that is uniformly sampled in radius. * You can also specify the desired PSF of the new profile. Since the * quantity returned is a statistical series, it remains possible to * relate the uncertainty in new samples to the uncertainties in the * original, and also to evaluate the covariances between new samples. *$ External_references: * Profile toolkit *$ Side_effects: * none *$ Examples: * Suppose we have a statistical series called stat containing a time * series of photon counts from a star during a ring occultation. There is * no dead time between samples so the PSF is effectively a boxcar of unit * width. We wish to resample the series into 1-km bins in the ring plane, * where each new sample is a uniform average of ring region falling within * a given bin. The new radius limits will be 139000. and 141000. km. We * already have a monotonic curve called radius that maps from the series' * sampling parameter to radius projected into the ring plane. * * PRO_OBJECT *oldpsf, *newpsf, *resamp; * * oldpsf = Pro_BoxcarFunc(1., 0.5); * newpsf = Pro_BoxcarFunc(1., 0.5); * resamp = Pro_ResampledStat(oldpsf, newpsf, radius, 1, 139000., 1., * 1, 2001); * * At this point, resamp is a statistical series with the desired * properties. Sample #1 gives the mean photon counts at 139000 km; * sample #2001 gives the mean counts at 141000 km. Because it is a * statistical series, it is also possible to use Pro_StatCovar to * request the variances of samples and the covariances between samples. *$ Error_handling: * Profile library error handling is in effect. * * Conditions raised: * RL_MEMORY_ERROR on memory allocation failure. * PRO_CLASS_ERROR if stat is NULL or is not a stat series; * if curve is NULL or not a curve; * if either psf is NULL or is not a supported * function. * PRO_COORD_MISMATCH if curve xname does not match stat xname. * PRO_DOMAIN_ERROR if segment is out of limits for curve. *$ Limitations: * none *$ Author_and_institution: * Mark R. Showalter * NASA/Ames Research Center *$ Version_and_date: * 1.0: March 1998 *$ Change_history: * none *******************************************************************************/ PRO_OBJECT *Pro_ResampledStat(stat, oldpsf, newpsf, curve, segment, x1, dx, k1, k2) PRO_OBJECT *stat, *oldpsf, *newpsf, *curve; RL_INT4 segment, k1, k2; RL_FLT8 x1, dx; { ZPRO_RESAMPLED *params; PRO_OBJECT *new; RL_FLT8 temp, height1, height2; /* Verify object types */ if (XPro_StatPtr(stat) == NULL || XPro_CurvePtr(curve) == NULL) return NULL; /* Check coordinate names */ if (strcmp(Pro_ObjectName(stat,1), Pro_ObjectName(curve,1)) != 0) { XPro_CoordMismatch("resampled series", stat, curve, Pro_ObjectName(stat,1), Pro_ObjectName(curve,1)); } /* Check segment index */ if (segment < 1 || segment > Pro_CurveSegments(curve)) { XPro_IDomainError("segment index", curve, 1, Pro_CurveSegments(curve), segment); return NULL; } /* Create new resampled statistical series structure */ params = (ZPRO_RESAMPLED *) XRL_Malloc(sizeof(ZPRO_RESAMPLED)); if (params == NULL) return NULL; /* Fill structure */ params->stat = stat; params->oldpsf = oldpsf; params->newpsf = newpsf; params->curve = curve; params->segment = segment; params->newx1 = x1; params->newdx = dx; params->newk1 = k1; params->newk2 = k2; params->oldtype = XPro_PSFInfo(oldpsf, &(params->oldhw), &(params->oldpar)); params->newtype = XPro_PSFInfo(newpsf, &(params->newhw), &(params->newpar)); if (params->oldtype == 0 || params->newtype == 0) { XRL_Free((RL_VOID *) params); return NULL; } (void) Pro_SeriesIndices(stat, &(params->oldk1), &(params->oldk2)); (void) Pro_SeriesSampling(stat, &(params->oldx1), NULL, &(params->olddx)); params->expand0 = params->newdx / params->olddx; (void) Pro_CurveExtremum(curve, segment-1, &(params->segment_min), NULL); (void) Pro_CurveExtremum(curve, segment, &(params->segment_max), NULL); if (params->segment_min > params->segment_max) { temp = params->segment_min; params->segment_min = params->segment_max; params->segment_max = temp; } /* Create new weighted statistical series object */ new = Pro_WeightedStat(stat, k1, k2, x1, dx, ZPro_ResampledKRange, ZPro_ResampledWeights, params, sizeof(ZPRO_RESAMPLED), TRUE, ZPro_PrintResampled); XRL_Free((RL_VOID *) params); /* Enslave old objects and check for errors */ if (new != NULL && (XPro_EnslaveObject(new, curve) || XPro_EnslaveObject(new, newpsf) || XPro_EnslaveObject(new, oldpsf))) { Pro_FreeObject(new); new = NULL; } /* NOTE: It is not necessary to enslave stat because Pro_WeightedStat() * does that for us! */ return new; } /* ******************************************************************************** *$ Component_name: * Pro_ResampledSize (resample.c) *$ Abstract: * This routine returns the size of a buffer needed for the series being * resampled. *$ Keywords: * PROFILE, SERIES * C, PUBLIC, SUBROUTINE *$ Declarations: * RL_INT4 Pro_ResampledSize(series, oldpsf, newpsf, curve, * segment, r1, r2, dr); * PRO_OBJECT *series, *oldpsf, *newpsf, *curve; * RL_INT4 segment; * RL_FLT8 r1, r2, dr; *$ Inputs: * series pointer to a series to resample. * oldpsf pointer to the PSF of the original series, possibly * unnormalized, in units of series samples. * newpsf pointer to the new PSF desired, possibly unnormalized, * in units of new series samples. * curve pointer to the curve that evaluates new coordinate from * old. * segment segment index of the curve to use for the resampling. * r1, r2, dr new x coordinate limits and interval. *$ Outputs: * none *$ Returns: * the needed size of a series buffer, or 0 on error. *$ Detailed_description: * This routine returns the size of a buffer needed for a series being * resampled. The estimate is somewhat conservative to ensure the value * returned is adequate. * * A buffer defines the finite subset of series that must be kept in memory * at any one time. Buffers are used by buffered series (see buffer.c) and * by column series (see column.c). It is very inefficient to filter or * resample a series when the size of the filter exceeds the segment of the * series kept in memory. Hence, this routine is useful when one wishes to * determine how big a buffer to use before calling Pro_ResampledStat() to * create the resampled series. *$ External_references: * Profile toolkit *$ Side_effects: * none *$ Examples: * Suppose we have a statistical series called stat containing a time * series of photon counts from a star during a ring occultation. There is * no dead time between samples so the PSF is effectively a boxcar of unit * width. We wish to resample the series into 1-km bins in the ring plane, * where each new sample is a uniform average of ring region falling within * a given bin. The new radius limits will be 139000. and 141000. km. We * already have a monotonic curve called radius that maps from the series' * sampling parameter to radius projected into the ring plane. * * integer*4 oldpsf, newpsf, resamp * * oldpsf = FPro_BoxcarFunc(1.d0, 0.5d0) * newpsf = FPro_BoxcarFunc(1.d0, 0.5d0) * resamp = FPro_ResampledStat(oldpsf, newpsf, radius, 1, 139000.d0, 1.d0, * & 1, 2001) * * At this point, resamp is a statistical series with the desired * properties. Sample #1 gives the mean photon counts at 139000 km; * sample #2001 gives the mean counts at 141000 km. Because it is a * statistical series, it is also possible to use Pro_StatCovar to * request the variances of samples and the covariances between samples. *$ Error_handling: * Profile library error handling is in effect. * * Conditions raised: * PRO_CLASS_ERROR if series is NULL or is not a series; * if curve is NULL or not a curve; * if either psf is NULL or is not a function. * PRO_DOMAIN_ERROR if segment is out of limits for curve. *$ Limitations: * none *$ Author_and_institution: * Mark R. Showalter * NASA/Ames Research Center *$ Version_and_date: * 1.0: March 1998 *$ Change_history: * none *******************************************************************************/ /* Number of steps to take when searching interval for minimum slope */ #define NSTEPS 50 RL_INT4 Pro_ResampledSize(series, oldpsf, newpsf, curve, segment, r1, r2, dr) PRO_OBJECT *series, *oldpsf, *newpsf, *curve; RL_INT4 segment; RL_FLT8 r1, r2, dr; { RL_FLT8 xr1, xr2, xtemp, min_drdx, drdx, x, xstep, size, dx, expand, oldwidth, newwidth; /* Find limits in old series' sampling parameter */ xr1 = Pro_CurveInverse(curve, segment, r1); xr2 = Pro_CurveInverse(curve, segment, r2); if (xr1 == xr2) return 0; if (xr1 > xr2) {xtemp = xr1; xr1 = xr2; xr2 = xtemp;} /* Find minimum absolute value of slope within domain * * Note that there is no Profile Tool to find an extremum of a slope. * Instead, we simply search the interval and assume the function is * not too perverse. */ Pro_CurveValue(curve, xr1, &min_drdx); if (min_drdx < 0.) min_drdx = -min_drdx; xstep = (xr2 - xr1) / NSTEPS; for (x = xr2; x > xr1; x -= xstep) { Pro_CurveValue(curve, x, &drdx); if (drdx < 0.) drdx = -drdx; if (drdx < min_drdx) min_drdx = drdx; } /* Now solve for the maximum expansion factor */ size = Pro_SeriesSampling(series, NULL, NULL, &dx); if (size == 0.) return 0; expand = fabs(dr / (min_drdx * dx)); /* Get the widths of the filter functions */ newwidth = Pro_ObjectDomain(newpsf, NULL, NULL); if (newwidth == 0) return 0; oldwidth = Pro_ObjectDomain(oldpsf, NULL, NULL); if (oldwidth == 0) return 0; /* Return a conservative estimate of the needed samples */ return 1 + (RL_INT4) (1.5 * (oldwidth + newwidth*expand)); } /* ******************************************************************************** * INTERNAL ROUTINES ******************************************************************************** * ZPro_ResampledKRange(knew, k1, k2, params) * * This internal routine returns the range of nonzero weights in the old series * corresponding to a given index in the new series. * * Inputs: * newk new index. * params pointer to the ZPRO_RESAMPLED structure. * * Outputs: * k1, k2 range of indices in old series for which weight is * nonzero. *******************************************************************************/ static void ZPro_ResampledKRange(newk, k1, k2, params) RL_INT4 newk, *k1, *k2; RL_VOID *params; { ZPRO_RESAMPLED *resamp; RL_FLT8 newx, oldx, dydx, center, expand, sumhw; resamp = (ZPRO_RESAMPLED *) params; /* Map new psf center to X-coordinate in new series */ newx = resamp->newx1 + (newk - resamp->newk1) * resamp->newdx; if (newx < resamp->segment_min || newx > resamp->segment_max) { *k1 = 1; *k2 = 0; return; } /* Map new psf center to X-coordinate in old series */ oldx = Pro_CurveInverse(resamp->curve, resamp->segment, newx); /* Map new psf center to sample index in old series */ center = resamp->oldk1 + (oldx - resamp->oldx1) / resamp->olddx; /* Determine expansion factor between sample grids */ (void) Pro_CurveValue(resamp->curve, oldx, &dydx); expand = fabs(resamp->expand0 / dydx); /* Determine halfwidth of convolved psfs */ sumhw = resamp->oldhw + expand*resamp->newhw - 1.; /* Return rounded values */ *k1 = (RL_INT4) floor(center - sumhw); *k2 = (RL_INT4) ceil(center + sumhw); } /* ******************************************************************************** * ZPro_ResampledWeight(knew, k, params) * * This internal routine returns the weight a sample in the old series has in the * new series. * * Inputs: * newk index in new series. * k index in old series. * params pointer to the ZPRO_RESAMPLED structure. * * Return: weight value. *******************************************************************************/ static RL_FLT8 ZPro_ResampledWeights(newk, k, params) RL_INT4 newk, k; RL_VOID *params; { ZPRO_RESAMPLED *resamp; RL_FLT8 newx, oldx, dydx, center, expand; resamp = (ZPRO_RESAMPLED *) params; /* Map new psf center to X-coordinate in new series */ newx = resamp->newx1 + (newk - resamp->newk1) * resamp->newdx; if (newx < resamp->segment_min || newx > resamp->segment_max) { return 0.; } /* Map new psf center to X-coordinate in old series */ oldx = Pro_CurveInverse(resamp->curve, resamp->segment, newx); /* Map new psf center to sample index in old series */ center = resamp->oldk1 + (oldx - resamp->oldx1) / resamp->olddx; /* Determine expansion factor between sample grids */ (void) Pro_CurveValue(resamp->curve, oldx, &dydx); expand = fabs(resamp->expand0 / dydx); /* * The integral we wish to return is * f[k] = Integral[ oldpsf(x-k) newpsf((x-center)/expand) dx ] * * The integral rearranges to * f[k] = Integral[ oldpsf(x) newpsf((x+(k-center))/expand) dx ] */ return XPro_PSFInteg(resamp->oldtype, resamp->newtype, 1./expand, (k-center)/expand, resamp->oldpar, resamp->newpar); } /* ******************************************************************************** * ZPro_PrintResampled(params) * * This internal routine prints information about a resampled stat series. * * Inputs: * params pointer to the ZPRO_RESAMPLED data structure. *******************************************************************************/ static void ZPro_PrintResampled(params) RL_VOID *params; { ZPRO_RESAMPLED *resamp; RL_INT4 i; resamp = (ZPRO_RESAMPLED *) params; /* Print resampled statistical series info... */ printf("\nResampled statistical series parameters...\n"); printf(" series = "); XPro_PrintInfo(resamp->stat); printf("original PSF = "); XPro_PrintInfo(resamp->oldpsf); printf(" new PSF = "); XPro_PrintInfo(resamp->newpsf); printf(" curve = "); XPro_PrintInfo(resamp->curve); printf(" segment = %d\n", resamp->segment); } /* ******************************************************************************** * FORTRAN INTERFACE ROUTINES ******************************************************************************** *$ Component_name: * FPro_ResampledStat (resample.c) *$ Abstract: * This routine resamples a statistical series so that it is uniformly * spaced in a new variable and based on a new point spread function (PSF). *$ Keywords: * PROFILE, SERIES * FORTRAN, PUBLIC, SUBROUTINE *$ Declarations: * integer*4 FPro_ResampledStat(stat, oldpsf, newpsf, curve, segment, * x1, dx, k1, k2) * integer*4 stat, oldpsf, newpsf, curve * integer*4 segment, k1, k2 * real*8 x1, dx *$ Inputs: * stat FORTRAN pointer to a statistical series to resample. * oldpsf FORTRAN pointer to the PSF of the original series, * possibly unnormalized, in units of series samples. * newpsf FORTRAN pointer to the new PSF desired, possibly * unnormalized, in units of new series samples. * curve FORTRAN pointer to the curve that evaluates new * coordinate from old. * segment segment index of the curve to use for the resampling. * x1, dx new x coordinate starting value and interval. * k1, k2 new range of indices. *$ Outputs: * none *$ Returns: * pointer to a new weighted statistical series; NULL on non-fatal error. *$ Detailed_description: * This routine resamples a statistical series so that it is uniformly * spaced in a new variable and based on a new point spread function (PSF). * * The returned statistical series is uniformly sampled in a new * coordinate, such that the give curve function converts from sampling * parameter of the original series to that of the new. For example, if * you have a ring occultation profile sampled uniformly in time plus a * curve that converts from time to ring-plane intercept radius, you can * quickly generate a revised profile that is uniformly sampled in radius. * You can also specify the desired PSF of the new profile. Since the * quantity returned is a statistical series, it remains possible to * relate the uncertainty in new samples to the uncertainties in the * original, and also to evaluate the covariances between new samples. *$ External_references: * Profile toolkit *$ Side_effects: * none *$ Examples: * none *$ Error_handling: * Profile library error handling is in effect. * * Conditions raised: * RL_MEMORY_ERROR on memory allocation failure. * PRO_CLASS_ERROR if stat is NULL or is not a stat series; * if curve is NULL or not a curve; * if either psf is NULL or is not a supported * function. * PRO_COORD_MISMATCH if curve xname does not match stat xname. * PRO_DOMAIN_ERROR if segment is out of limits for curve. * FORTRAN_POINTER_ERROR if any object is not a valid FORTRAN pointer. *$ Limitations: * none *$ Author_and_institution: * Mark R. Showalter * NASA/Ames Research Center *$ Version_and_date: * 1.0: March 1998 *$ Change_history: * none *******************************************************************************/ RL_INT4 FORTRAN_NAME(fpro_resampledstat) (stat, oldpsf, newpsf, curve, segment, x1, dx, k1, k2) RL_INT4 *stat, *oldpsf, *newpsf, *curve; RL_INT4 *segment, *k1, *k2; RL_FLT8 *x1, *dx; { RL_INT4 index; RL_VOID *statptr, *oldpsfptr, *newpsfptr, *curveptr, *newptr; /* Look up pointers */ statptr = FORT_GetPointer(*stat); if (statptr == NULL) return 0; oldpsfptr = FORT_GetPointer(*oldpsf); if (oldpsfptr == NULL) return 0; newpsfptr = FORT_GetPointer(*newpsf); if (newpsfptr == NULL) return 0; curveptr = FORT_GetPointer(*curve); if (curveptr == NULL) return 0; /* Call function */ newptr = (RL_VOID *) Pro_ResampledStat((PRO_OBJECT *) statptr, (PRO_OBJECT *) oldpsfptr, (PRO_OBJECT *) newpsfptr, (PRO_OBJECT *) curveptr, *segment, *x1, *dx, *k1, *k2); if (newptr == NULL) return 0; /* Save new pointer */ index = FORT_AddPointer(newptr); if (index == 0) Pro_FreeObject((PRO_OBJECT *) newptr); return index; } /* ******************************************************************************** *$ Component_name: * FPro_ResampledSize (resample.c) *$ Abstract: * This routine returns the size of a buffer needed for the series being * resampled. *$ Keywords: * PROFILE, SERIES * FORTRAN, PUBLIC, SUBROUTINE *$ Declarations: * integer function FPro_ResampledSize(series, oldpsf, newpsf, curve, * segment, r1, r2, dr) * integer*4 series, oldpsf, newpsf, curve, segment * real*8 r1, r2, dr *$ Inputs: * series FORTRAN pointer to the series to resample. * oldpsf FORTRAN pointer to the PSF of the original series, * possibly unnormalized, in units of series samples. * newpsf FORTRAN pointer to the new PSF desired, possibly * unnormalized, in units of new series samples. * curve FORTRAN pointer to the curve that evaluates new * coordinate from old. * segment segment index of the curve to use for the resampling. * r1, r2, dr new x coordinate limits and interval. *$ Outputs: * none *$ Returns: * the needed size of a series buffer, or 0 on error. *$ Detailed_description: * This routine returns the size of a buffer needed for a series being * resampled. The estimate is somewhat conservative to ensure the value * returned is adequate. * * A buffer defines the finite subset of series that must be kept in memory * at any one time. Buffers are used by buffered series (see buffer.c) and * by column series (see column.c). It is very inefficient to filter or * resample a series when the size of the filter exceeds the segment of the * series kept in memory. Hence, this routine is useful when one wishes to * determine how big a buffer to use before calling Pro_ResampledStat() to * create the resampled series. *$ External_references: * Profile toolkit *$ Side_effects: * none *$ Examples: * none *$ Error_handling: * Profile library error handling is in effect. * * Conditions raised: * PRO_CLASS_ERROR if series is NULL or is not a series; * if curve is NULL or not a curve; * if either psf is NULL or is not a function. * PRO_DOMAIN_ERROR if segment is out of limits for curve. * FORTRAN_POINTER_ERROR if any object is not a valid FORTRAN pointer. *$ Limitations: * none *$ Author_and_institution: * Mark R. Showalter * NASA/Ames Research Center *$ Version_and_date: * 1.0: March 1998 *$ Change_history: * none *******************************************************************************/ RL_INT4 FORTRAN_NAME(fpro_resampledsize) (series, oldpsf, newpsf, curve, segment, r1, r2, dr) RL_INT4 *series, *oldpsf, *newpsf, *curve, *segment; RL_FLT8 *r1, *r2, *dr; { RL_VOID *seriesptr, *oldpsfptr, *newpsfptr, *curveptr, *newptr; /* Look up pointers */ seriesptr = FORT_GetPointer(*series); if (seriesptr == NULL) return 0; oldpsfptr = FORT_GetPointer(*oldpsf); if (oldpsfptr == NULL) return 0; newpsfptr = FORT_GetPointer(*newpsf); if (newpsfptr == NULL) return 0; curveptr = FORT_GetPointer(*curve); if (curveptr == NULL) return 0; /* Call function */ return Pro_ResampledSize((PRO_OBJECT *) seriesptr, (PRO_OBJECT *) oldpsfptr, (PRO_OBJECT *) newpsfptr, (PRO_OBJECT *) curveptr, *segment, *r1, *r2, *dr); } /********************************************************************************/