/* Univariate polynomial splines defined on irregular dyadic grids */ /* Last edited on 2004-08-24 02:12:12 by stolfi */ #ifndef dgpulse_H #define dgpulse_H #include #include #include /* UNIVARIATE POLYNOMIAL SPLINES A /univariate polynomial spline/ is a spline whose domain {D} is a subset of the real line. FINITE PULSES A univariate finite element is a /finite pulse/. STANDARD DYADIC PULSES For unidimensional dyadic splines, the bases that we consider here consist of /standard dyadic finite pulses/, or /dyadic pulses/ for short, which are characterized by the following parameters: * an enumerated /kind/, e.g. {DG_PK_H}, {DG_PK_B}, etc. * its order of continuity {c\geq -1} * its degree {g \geq 0} * a /pulse index/ {pix \geq 0} * a /rank/ {r\geq 0} * a /shift/ {j} in the range {0..2^r-1} Any dyadic pulse {b} is a polynomial spline of continuity {c} and degree {g}, whose mesh is level {r} of the unidimensional dyadic multigrid {G*} -- namely on a partition of the circular unit interval {T = [0_1)} into {2^r} equal cells (intervals) and {2^r} corners. STANDARD PULSE KINDS This module presently implements two kinds of finite pulses: /H-pulses/ ({kind = DG_PK_H}) are useful for Hermite-style interpolation of data given at grid corners. /B-pulses/ ({kind = DG_PK_B}) have the property that they are a partition of unity (non-negative and add to 1 everywhere). When {c = 0, g = 1}, the two kinds coincide with the /piecewise linear/ splines -- more properly called /piecewise affine/. */ typedef enum { DG_PK_B, DG_PK_H } dg_PulseKind; typedef int32 dg_PulseIndex; /* Identifies a finite pulse among those with same kind, degree, rank, and shift. */ double dg_pulse_eval ( dg_PulseKind kind, dg_Cont c, dg_Degree g, dg_PulseIndex pix, dg_GridSize sz, dg_GridPos j, double x ); /* Evaluates the dyadic pulse of kind {kind}, continuity {c}, degree {g}, pulse index {pix}, rank {r} such that {sz=2^r}, and shift {j} in {0..sz-1}, at the argument {x}. Takes wrap-around into account when appropriate (see below). */ /* STANDARD DYADIC PULSE SUPPORT The support of a standard dyadic pulse of rank {r}, on level {r} of the univariate multigrid {G*}, consists of some number {wd} -- the pulse's /support count/ -- of consecutive cells, possibly including their lower vertices. By convention, the shift {j} specifies the grid position of the *last* cell in the pulse's support. That is, the support is the interval {[(j-wd+1)/2^r _ (j+1)/2^r)}, taken modulo 1, except possibly for some vertices contained in it. For large enough {r}, the support count {wd} depends only on {kind}, {c}, {g}, and {pix}. For small {r}, the count {wd} is reduced by self-overlap (see below). */ dg_GridSize dg_pulse_supp_count ( dg_PulseKind kind, dg_Cont c, dg_Degree g, dg_PulseIndex pix, dg_GridSize sz ); /* Returns the number {wd} of consecutive cells in the support of the standard dyadic pulse with given {kind}, continuity {c}, degree {g}, pulse index {pix}, and rank {r} such that {sz = 2^r}. Takes self-overlap into account, so that the result never exceeds {sz}. */ /* SELF-SIMILARITY AND MOTHER PULSES Dyadic pulses that differ only in rank and shift are related to each other by argument scaling and translation (like wavelets). The other four parameters -- {kind}, continuity {c}, degree {g}, and pulse index {pix} -- define the basic shape of the pulse, and thus will be called the /shape parameters/. All dyadic pulses with the same shape parameters are derived from a single /mother pulse/ {w(x)}. The mother pulse is a polynomial spline of continuity order {c\geq -1} and degree {g > c}, defined on the non-periodic integer grid -- the uniform infinite grid on the real line with integer vertices. The dyadic pulse {b} of rank {r} and shift {j} with those shape parameters is given by {b[r,j](x) = SUM { w(z) : z \eqv 2^r*x - j (mod 2^r) } } The mother pulse {w} has always bounded support, which by convention is the interval {[-(m-1)_1]}, where {m} is the /natural support count/ for pulses of that shape. Thus, in the above definition, there are at most {ceil(m/2^r)} arguments {z} such that {w(z)} is non-zero, so the above sum is finite. In fact, when the rank is large enough (namely {2^r \geq m}), there is at most one {z} in that interval. Then every dyadic pulse {b} with those shape parameters looks like a scaled and translated version of {w}, and its support consists of {m} consecutive cells, whose positions are congruent to {j-m+1} through {j}, modulo {2^r}. Note also that while the mother pulse {w} is not periodic, the derived dyadic pulses are periodic with period 1. SELF-OVERLAP When the rank {r} is too low (namely {2^r < m}), the reduction of the argument {2^r*x - j} modulo {2^r} yields several {z} within the mother pulse's support. Intuitively, the mother pulse gets folded over itself when the grid is too small; in that case, the overlapped parts are added together to give the dyadic pulse {b}. */ double dg_mother_pulse_eval ( dg_PulseKind kind, dg_Cont c, dg_Degree g, dg_PulseIndex pix, double z ); /* Evaluates the mother pulse of kind {kind}, continuity {c}, degree {g}, and pulse index {pix}, at argument {z}. The result is zero if {z} is outside the interval {[-m+1_1)} where {m = dg_mother_pulse_supp_count(kind,c,g,pix)}. */ dg_GridSize dg_mother_pulse_supp_count ( dg_PulseKind kind, dg_Cont c, dg_Degree g, dg_PulseIndex pix ); /* Computes the natural support count {m} of the mother spline of the given {kind}, with continuity {c}, degree {g}, and pulse index {pix}. */ void dg_mother_pulse_to_bezier ( dg_PulseKind kind, dg_Cont c, dg_Degree g, dg_PulseIndex pix, double **bp ); /* Splits the {d}-dimensional mother pulse {t} of the given {kind}, continuity order {c}, degree {d}, and pulse index {pix} into its {NP} constituent polynomial patches {pp[0..NP-1]}, computes the vector of {g+1} Bézier coeffs of each patch {tp[k]}, and stores that array into the given array {bp[k][0..g]}. If {bp[k]} is NULL, an appropriate coeff array is automatically allocated and initialized to zeros. */ dg_Degree dg_mother_pulse_min_degree ( dg_PulseKind kind, dg_Cont c); dg_Degree dg_mother_pulse_max_degree ( dg_PulseKind kind, dg_Cont c); /* Computes the natural support count {m} of the mother spline of the given {kind}, with continuity {c}, degree {g}, and pulse index {pix}. */ double dg_mother_pulse_sum ( dg_PulseKind kind, dg_Cont c, dg_Degree g, dg_PulseIndex pix, dg_GridSize sz, double z ); /* Evaluates the mother pulse of kind {kind}, continuity {c}, degree {g}, and index {pix}, summed over all arguments that are congruent to {z} modulo {sz}. */ dg_Interval dg_mother_pulse_sum_range ( dg_PulseKind kind, dg_Cont c, dg_Degree g, dg_PulseIndex pix, dg_GridSize sz ); /* The range of values of {dg_mother_pulse_sum(kind,c,g,pix,sz,z)}, for any argument {z}. Useful for plotting, but beware: due to roudoff errors, the returned range may be slightly smaller than the true one. */ /* H-PULSES The dyadic H-pulses are defined only for {c\geq 0}, {g = 2*c+1}. There are exactly {c+1} H-pulses for each rank {r} and shift {j}, distinguished by their index {pix} in {0..c}. The mother H-pulses are all supported by the open interval {(-1_+1)}. At -1, 0, and +1, the mother pulse {w^g_c[pix]} of index {pix} has all derivatives up to order {c} equal to zero at integer {z} values, except for the derivative of order {pix} which is equal to {choose(g,pix)}. The mother pulse is symmetric or anti-symmetric, depending on the parity of {pix}; so, for {z} in {(-1_0]}, { w^g_c[pix](z) = (-1)^pix w^g_c[pix](-z) } It follows that a dyadic H-pulse of rank {r} and shift {j} is zero for any argument {x} that lies outside the interval {((j-1)/2^r_(j+1)/2^r)}, when both are taken modulo 1. Thus, for {r\geq 1} its support consists of the grid vertex {j/2^r} and the two adjacent grid cells, with positions {j-1} and {j}. For {r = 0} the grid has only one cell and the pulse overlaps itself -- that is, the two halves of the mother pulse are added to give the H-pulse. The {c+1} shapes for a given {c} remain distinct in spite of the self-overlap. */ double dg_mother_H_pulse_eval ( dg_Cont c, dg_PulseIndex pix, double z ); /* Computes the mother H-pulse of continuity {c \geq 0}, degree {g = 2*c+1}), and index {pix} (which must lie in {0..c}) for argument {z}. */ /* B-PULSES The dyadic B-pulses are defined for any {c \geq -1} and any non-negative degree {g \geq 2*c+1}. For a given {g} and {c}, there are {g-c} distinct mother pulses, with pulse index {pix} ranging from 0 to {g-c-1}. They fall into two sub-classes: * When {pix} is in {c+1..g-c-1}, the mother spline is supported on {(0_1)} only, and coincides with the Bernstein-Bézier polynomial {BB^g_pix}: { w^g_c[pix](z) = BB^g_pix(z) = choose(g,pix) z^{g-pix} (1-z)^pix } In this case, the dyadic B-pulses of rank {r} and shift {j} are supported by the interval {[j/2^r_(j+1)/2^r)}. * When {pix} is in {0..c}, the mother spline {w^g_c[pix]} is supported on the interval {(-1_+1)}. In the positive half {[0_1)}, the mother pulse can be expressed as a combination of Bernstein-Bézier polynomials: { w^g_c[pix](z) = SUM { choose(c-k,pix-k)/2^{c-k} BB^g_{k}(z) : k = 0..pix } } For {z} in {[-1_0]}, the following relation holds { w^g_c[pix](z) = w^g_c[c-pix](-z) } In this case, the dyadic B-pulses of rank {r} and shift {j} are supported by the interval {[(j-1)/2^r_ (j+1)/2^r)}, reduced modulo 1. */ /* Maximum degree for a B-pulse, for safety: */ #define DG_MAX_B_PULSE_DEGREE 20 double dg_mother_B_pulse_eval ( dg_Cont c, dg_Degree g, dg_PulseIndex pix, double z ); /* Computes the univariate mother B-pulse of continuity {c}, degree {g} (at least {2*c + 1}), and index {pix} (which must lie in {0..g-c-1}) for the argument {z}. */ #endif