lprefilt.c

Go to the documentation of this file.

#include <stdlib.h>
#include <math.h>
#include <fftw3.h>
#include "adaptlob.h"
#include "lprefilt.h"


static void PrefilterScan(float *Data, int Stride, int N, float alpha,
    boundaryhandling Boundary)
{
    const float Eps = 1e-4;
    float Sum, Weight, Last;
    int i, iEnd, n0;

    n0 = (int)ceil(log(Eps)/log(fabs(alpha)));

    if(n0 > N)
        n0 = N;

    switch(Boundary)
    {
        case BOUNDARY_CONSTANT:
            Sum = Data[0]/(1 - alpha);
            break;
        case BOUNDARY_WSYMMETRIC:
            Sum = Data[0];
            Weight = 1;
            iEnd = n0*Stride;

            for(i = Stride; i < iEnd; i += Stride)
            {
                Weight *= alpha;
                Sum += Data[i]*Weight;
            }
            break;
        default: /* BOUNDARY_HSYMMETRIC */
            Sum = Data[0]*(1 + alpha);
            Weight = alpha;
            iEnd = n0*Stride;

            for(i = Stride; i < iEnd; i += Stride)
            {
                Weight *= alpha;
                Sum += Data[i]*Weight;
            }
            break;
    }

    Last = Data[0] = Sum;
    iEnd = (N - 1)*Stride;

    for(i = Stride; i < iEnd; i += Stride)
    {
        Data[i] += alpha*Last;
        Last = Data[i];
    }

    switch(Boundary)
    {
        case BOUNDARY_CONSTANT:
            Last = Data[iEnd] = (alpha*(-Data[iEnd] + (alpha - 1)*alpha*Last))
                /((alpha - 1)*(alpha*alpha - 1));
            break;
        case BOUNDARY_HSYMMETRIC:
            Data[iEnd] += alpha*Last;
            Last = Data[iEnd] *= alpha/(alpha - 1);
            break;
        case BOUNDARY_WSYMMETRIC:
            Data[iEnd] += alpha*Last;
            Last = Data[iEnd] = (alpha/(alpha*alpha - 1))
                * ( Data[iEnd] + alpha*Data[iEnd - Stride] );
            break;
    }

    for(i = iEnd - Stride; i >= 0; i -= Stride)
    {
        Data[i] = alpha*(Last - Data[i]);
        Last = Data[i];
    }
}


void PrefilterImage(float *Data, int Width, int Height, int NumChannels,
    const float *alpha, int NumFilterPairs, float ConstantFactor,
    boundaryhandling Boundary)
{
    const int NumPixels = Width*Height;
    int k, x, y, Channel;


    /* Square the ConstantFactor for two spatial dimensions */
    ConstantFactor = ConstantFactor*ConstantFactor;

    for(Channel = 0; Channel < NumChannels; Channel++)
    {
        for(x = 0; x < Width; x++)
            for(k = 0; k < NumFilterPairs; k++)
                PrefilterScan(Data + x, Width, Height, alpha[k], Boundary);

        for(y = 0; y < Height; y++)
            for(k = 0; k < NumFilterPairs; k++)
                PrefilterScan(Data + Width*y, 1, Width, alpha[k], Boundary);

        for(k = 0; k < NumPixels; k++)
            Data[k] *= ConstantFactor;

        Data += NumPixels;
    }
}


typedef struct {
    float (*Kernel)(float);
    float KernelRadius;
    float (*Psf)(float, const void*);
    const void *PsfParams;
    float x;
} convolutionparams;


static float ConvIntegrand(float t, const void *Params)
{
    float (*Kernel)(float) = ((convolutionparams*)Params)->Kernel;
    float (*Psf)(float, const void*) = ((convolutionparams*)Params)->Psf;
    const void *PsfParams = ((convolutionparams*)Params)->PsfParams;
    float x = ((convolutionparams*)Params)->x;

    return Kernel(t) * Psf(x - t, PsfParams);
}


static float ConvIntegrandNormalized(float t, const void *Params)
{
    float (*Kernel)(float) = ((convolutionparams*)Params)->Kernel;
    const int R = 2*ceil(((convolutionparams*)Params)->KernelRadius);
    float (*Psf)(float, const void*) = ((convolutionparams*)Params)->Psf;
    const void *PsfParams = ((convolutionparams*)Params)->PsfParams;
    float Sum, x = ((convolutionparams*)Params)->x;
    int r;

    for(r = -R, Sum = 0; r <= R; r++)
        Sum += Kernel(t - r);

    return (Kernel(t)/Sum) * Psf(x - t, PsfParams);
}


void PsfConvCoeff(float *Coeff, int NumCoeffs,
    float (*Psf)(float, const void*), const void *PsfParams,
    float (*Kernel)(float), float KernelRadius, int KernelNormalize)
{
    float (*ConvFun)(float, const void*) = ConvIntegrand;
    convolutionparams ConvParams;
    int m;

    ConvParams.Kernel = Kernel;
    ConvParams.KernelRadius = KernelRadius;
    ConvParams.Psf = Psf;
    ConvParams.PsfParams = PsfParams;

    if(KernelNormalize)
        ConvFun = ConvIntegrandNormalized;

    for(m = 0; m < NumCoeffs; m++)
    {
        ConvParams.x = m;
        Coeff[m] = AdaptLob(ConvFun, -KernelRadius, KernelRadius,
            1e-7f, (const void *)&ConvParams);
    }
}


static int PsfPreFilterScan(float *Data, int Stride, int ScanSize,
    int ScanStride, int NumScans, int ChannelStride, int NumChannels,
    const float *Coeff, int NumCoeffs, boundaryhandling Boundary)
{
    const int PadCoeff = ScanSize + ((Boundary == BOUNDARY_HSYMMETRIC) ? 1:0);
    float *Buf = NULL, *CoeffDct = NULL;
    fftwf_plan Plan = 0;
    fftwf_iodim Dims[1], HowManyDims[2];
    fftw_r2r_kind Kind[1];
    int i, Scan, Channel, Success = 0;


    if((Boundary != BOUNDARY_HSYMMETRIC && Boundary != BOUNDARY_WSYMMETRIC)
        || !(Buf = (float *)fftwf_malloc(sizeof(float)*ScanSize*NumScans*NumChannels))
        || !(CoeffDct = (float *)fftwf_malloc(sizeof(float)*PadCoeff)))
        goto Catch;

    for(i = 0; i < NumCoeffs && i < PadCoeff; i++)
        Buf[i] = Coeff[i];
    for(; i < PadCoeff; i++)
        Buf[i] = 0;

    /* Perform DCT-I */
    if(!(Plan = fftwf_plan_r2r_1d(PadCoeff, Buf, CoeffDct,
        FFTW_REDFT00, FFTW_ESTIMATE | FFTW_DESTROY_INPUT)))
        goto Catch;

    fftwf_execute(Plan);
    fftwf_destroy_plan(Plan);

    /* Incorporate the normalization scale factor into the coefficients */
    for(i = 0; i < ScanSize; i++)
        CoeffDct[i] *= 2*(PadCoeff - 1);

    /* Forward transform of the image data */
    Dims[0].n = ScanSize;
    Dims[0].is = Stride;
    Dims[0].os = 1;
    HowManyDims[0].n = NumChannels;
    HowManyDims[0].is = ChannelStride;
    HowManyDims[0].os = ScanSize*NumScans;
    HowManyDims[1].n = NumScans;
    HowManyDims[1].is = ScanStride;
    HowManyDims[1].os = ScanSize;

    /* Use DCT-II for half-sample symmetric boundaries,
    DCT-I for whole-sample symmetric. */
    Kind[0] = (Boundary == BOUNDARY_HSYMMETRIC) ? FFTW_REDFT10 : FFTW_REDFT00;

    if(!(Plan = fftwf_plan_guru_r2r(1, Dims, 2, HowManyDims, Data,
        Buf, Kind, FFTW_ESTIMATE | FFTW_DESTROY_INPUT)))
        goto Catch;

    fftwf_execute(Plan);
    fftwf_destroy_plan(Plan);

    /* Divide */
    for(Channel = 0; Channel < NumChannels; Channel++)
    {
        for(Scan = 0; Scan < NumScans; Scan++)
        {
            for(i = 0; i < ScanSize; i++)
                Buf[i + ScanSize*(Scan + NumScans*Channel)] /= CoeffDct[i];
        }
    }

    /* Perform inverse transform */
    Dims[0].n = ScanSize;
    Dims[0].is = 1;
    Dims[0].os = Stride;
    HowManyDims[0].n = NumChannels;
    HowManyDims[0].is = ScanSize*NumScans;
    HowManyDims[0].os = ChannelStride;
    HowManyDims[1].n = NumScans;
    HowManyDims[1].is = ScanSize;
    HowManyDims[1].os = ScanStride;

    /* Use DCT-III (inverse of DCT-II) for half-sample symmetric boundaries,
    DCT-I (inverse of DCT-I) for whole-sample symmetric. */
    Kind[0] = (Boundary == BOUNDARY_HSYMMETRIC) ? FFTW_REDFT01 : FFTW_REDFT00;

    if(!(Plan = fftwf_plan_guru_r2r(1, Dims, 2, HowManyDims, Buf,
        Data, Kind, FFTW_ESTIMATE | FFTW_DESTROY_INPUT)))
        goto Catch;

    fftwf_execute(Plan);
    fftwf_destroy_plan(Plan);
    Success = 1;
Catch:
    if(CoeffDct)
        fftwf_free(CoeffDct);
    if(Buf)
        fftwf_free(Buf);
    fftwf_cleanup();
    return Success;
}


int PsfPreFilter(float *Data, int Width, int Height, int NumChannels,
    const float *Coeff, int NumCoeffs, boundaryhandling Boundary)
{
    if(!Data || !Coeff)
        return 0;

    if(!PsfPreFilterScan(Data, Width, Height, 1, Width, Width*Height,
            NumChannels, Coeff, NumCoeffs, Boundary))
        return 0;

    if(!PsfPreFilterScan(Data, 1, Width, Width, Height, Width*Height,
            NumChannels, Coeff, NumCoeffs, Boundary))
        return 0;

    return 1;
}