Biomedical Imaging Group Department of Microengineering

One-dimensional B-spline coefficients

void    ConvertToInterpolationCoefficients
        (
            double  c[],        /* input samples --> output coefficients */
            long    DataLength, /* number of samples or coefficients */
            double  z[],        /* poles */
            long    NbPoles,    /* number of poles */
            double  Tolerance   /* admissible relative error */
        )

{ /* begin ConvertToInterpolationCoefficients */

    double  Lambda = 1.0;
    long    n, k;

    /* special case required by mirror boundaries */
    if (DataLength == 1L) return;
    /* compute the overall gain */
    for (k = 0L; k < NbPoles; k++) Lambda = Lambda * (1.0 - z[k]) * (1.0 - 1.0 / z[k]);
    /* apply the gain */
    for (n = 0L; n < DataLength; n++) c[n] *= Lambda;
    /* loop over all poles */
    for (k = 0L; k < NbPoles; k++) {
        /* causal initialization */
        c[0] = InitialCausalCoefficient(c, DataLength, z[k], Tolerance);
        /* causal recursion */
        for (n = 1L; n < DataLength; n++) c[n] += z[k] * c[n - 1L];
        /* anticausal initialization */
        c[DataLength - 1L] = InitialAntiCausalCoefficient(c, DataLength, z[k]);
        /* anticausal recursion */
        for (n = DataLength - 2L; 0 <= n; n--) c[n] = z[k] * (c[n + 1L] - c[n]);
    }
} /* end ConvertToInterpolationCoefficients */

double  InitialCausalCoefficient
        (
            double  c[],        /* coefficients */
            long    DataLength, /* number of coefficients */
            double  z,          /* actual pole */
            double  Tolerance   /* admissible relative error */
        )

{ /* begin InitialCausalCoefficient */

    double  Sum, zn, z2n, iz;
    long    n, Horizon;

    /* this initialization corresponds to mirror boundaries */
    Horizon = DataLength;
    if (Tolerance > 0.0) Horizon = (long)ceil(log(Tolerance) / log(fabs(z)));
    if (Horizon < DataLength) {
        /* accelerated loop */
        zn = z;
        Sum = c[0];
        for (n = 1L; n < Horizon; n++) {
            Sum += zn * c[n];
            zn *= z;
        }
        return(Sum);
    }
    else {
        /* full loop */
        zn = z;
        iz = 1.0 / z;
        z2n = pow(z, (double)(DataLength - 1L));
        Sum = c[0] + z2n * c[DataLength - 1L];
        z2n *= z2n * iz;
        for (n = 1L; n <= DataLength - 2L; n++) {
            Sum += (zn + z2n) * c[n];
            zn *= z;
            z2n *= iz;
        }
        return(Sum / (1.0 - zn * zn));
    }
} /* end InitialCausalCoefficient */

double  InitialAntiCausalCoefficient
        (
            double  c[],        /* coefficients */
            long    DataLength, /* number of samples or coefficients */
            double  z           /* actual pole */
        )

{ /* begin InitialAntiCausalCoefficient */

    /* this initialization corresponds to mirror boundaries */
    return((z / (z * z - 1.0)) * (z * c[DataLength - 2L] + c[DataLength - 1L]));
} /* end InitialAntiCausalCoefficient */

Comments and feedback to: philippe.thevenaz@epfl.ch