conv.c

Go to the documentation of this file.

#include <string.h>
#include "conv.h"


#define CLAMP(X,A,B)    (((X) < (A)) ? (A) : (((X) > (B)) ? (B) : (X)))


const filter NullFilter = {NULL, 0, 0};


void SampledConv1D(float *Dest, int DestStride, const float *Src,
    int SrcStride, filter Filter, boundaryext Boundary, int N, 
    int nStart, int nStep, int nEnd)
{    
    const int SrcStrideStep = SrcStride*nStep;
    const int LeftmostTap = 1 - Filter.Delay - Filter.Length;    
    const int StartInterior = CLAMP(-LeftmostTap, 0, N - 1);
    const int EndInterior = (Filter.Delay <  0) ? 
                            (N + Filter.Delay - 1) : (N - 1);    
    const float *SrcN, *SrcK;
    float Accum;
    int n, k;

    
    if(nEnd < nStart || nStep <= 0 || N <= 0)
        return;
    
    /* Handle the left boundary */
    for(n = nStart; n < StartInterior; n += nStep, Dest += DestStride)
    {
        for(k = 0, Accum = 0; k < Filter.Length; k++)
            Accum += Filter.Coeff[k] 
                * Boundary(Src, SrcStride, N, n - Filter.Delay - k);
        
        *Dest = Accum;
    }
    
    /* Compute the convolution on the interior of the signal:
    
       In the inner accumulation loop
       SrcK = &inputdata[n - FilterDelay - k],  k = FilterLength-1, ..., 0.
       
       The SrcN pointer is adjusted such that
       SrcN = &inputdata[n + LeftmostTap]. 
       
       If n == StartInterior, then the loop starts with
          n = -LeftmostTap, SrcN = &inputdata[0]  if LeftmostTap <= 0
          n = 0, SrcN = &inputdata[LeftmostTap]   if LeftmostTap >= 0. */
    SrcN = (LeftmostTap <= 0) ? Src : (Src + SrcStride*LeftmostTap);
    
    /* Adjust if n > StartInterior */
    SrcN += SrcStride*(n - StartInterior);
    
    for(; n <= EndInterior; n += nStep, SrcN += SrcStrideStep, Dest += DestStride)
    {
        Accum = 0;
        SrcK = SrcN;
        k = Filter.Length;
        
        while(k)
        {
            Accum += Filter.Coeff[--k] * (*SrcK);
            SrcK += SrcStride;
        }
        
        *Dest = Accum;
    }
    
    /* Handle the right boundary */
    for(; n <= nEnd; n += nStep, Dest += DestStride)
    {
        for(k = 0, Accum = 0; k < Filter.Length; k++)
            Accum += Filter.Coeff[k]
                * Boundary(Src, SrcStride, N, n - Filter.Delay - k);
        
        *Dest = Accum;
    }
}


void SeparableConv2D(float *Dest, float *Buffer, const float *Src,
    filter FilterX, filter FilterY, boundaryext Boundary, 
    int Width, int Height, int NumChannels)
{
    const int NumPixels = Width*Height;
    int i, Channel;
    
    for(Channel = 0; Channel < NumChannels; Channel++)
    {
        /* Filter Src horizontally and store the result in Buffer */
        for(i = 0; i < Height; i++)
            Conv1D(Buffer + Width*i, 1, Src + Width*i, 1, 
                FilterX, Boundary, Width);
        
        /* Filter Buffer vertically and store the result in Dest */
        for(i = 0; i < Width; i++)
            Conv1D(Dest + i, Width, Buffer + i, Width, FilterY, 
                Boundary, Height);
            
        Src += NumPixels;
        Dest += NumPixels;
    }
}


filter MakeFilter(float *Coeff, int Delay, int Length)
{
    filter Filter;
    
    Filter.Coeff = Coeff;
    Filter.Delay = Delay;
    Filter.Length = Length;
    return Filter;
}


filter AllocFilter(int Delay, int Length)
{
    float *Coeff;
    
    if(Length > 0 && (Coeff = (float *)Malloc(sizeof(float)*Length)))
        return MakeFilter(Coeff, Delay, Length);
    else
        return NullFilter;
}


int IsNullFilter(filter Filter)
{
    return (Filter.Coeff == NULL) ? 1:0;
}


filter GaussianFilter(double Sigma, int R)
{
    filter Filter = AllocFilter(-R, 2*R+1);
    

    if(!IsNullFilter(Filter))
    {
        if(Sigma == 0)
            Filter.Coeff[0] = 1;
        else
        {
            float Sum;
            int r;
            
            for(r = -R, Sum = 0; r <= R; r++)
            {
                Filter.Coeff[R + r] = (float)exp(-r*r/(2*Sigma*Sigma));
                Sum += Filter.Coeff[R + r];
            }
            
            for(r = -R; r <= R; r++)
                Filter.Coeff[R + r] /= Sum;
        }
    }
    
    return Filter;
}


static float ZeroPaddedExtension(const float *Src, int Stride, int N, int n)
{
    return (0 <= n && n < N) ? Src[Stride*n] : 0;
}


static float ConstantExtension(const float *Src, int Stride, int N, int n)
{
    return Src[(n < 0) ? 0 : ((n >= N) ? Stride*(N - 1) : n)];
}


static float LinearExtension(const float *Src, int Stride, int N, int n)
{
    if(0 <= n)
    {
        if(n < N)
            return Src[Stride*n];
        else if(N == 1)
            return Src[0];
        else
        {
            Src += Stride*(N - 1);
            return Src[0] + (N - 1 - n)*(Src[-Stride] - Src[0]);
        }
    }
    else if(N == 1)
        return Src[0];
    else
        return Src[0] + n*(Src[Stride] - Src[0]);
}


static float PeriodicExtension(const float *Src, int Stride, int N, int n)
{
    if(n < 0)
    {
        do
        {
            n += N;
        }while(n < 0);
    }
    else if(n >= N)
    {
        do
        {
            n -= N;
        }while(n >= N);
    }
    
    return Src[Stride*n];
}


static float SymhExtension(const float *Src, int Stride, int N, int n)
{
    while(1)
    {
        if(n < 0)
            n = -1 - n;
        else if(n >= N)
            n = 2*N - 1 - n;
        else
            break;
    }
    
    return Src[Stride*n];
}


static float SymwExtension(const float *Src, int Stride, int N, int n)
{
    while(1)
    {
        if(n < 0)
            n = -n;
        else if(n >= N)
            n = 2*(N - 1) - n;
        else
            break;
    }
    
    return Src[Stride*n];
}


static float AsymhExtension(const float *Src, int Stride, int N, int n)
{
    float Jump, Offset;
    
    
    /* Use simple formulas for -N <= n <= 2*N - 1 */
    if(0 <= n)
    {
        if(n < N)
            return Src[Stride*n];        
        else if(n <= 2*N - 1)
            return 3*Src[Stride*(N - 1)] - Src[Stride*(N - 2)]
                - Src[Stride*(2*N - 1 - n)];
    }
    else if(-N <= n)
        return 3*Src[0] - Src[Stride] - Src[Stride*(-1 - n)];
    
    /* N == 1 is a special case */
    if(N == 1)
        return Src[0];
    
    /* General formula for extension at an arbitrary n */
    Jump = 3*(Src[Stride*(N - 1)] - Src[0]) 
        - (Src[Stride*(N - 2)] - Src[Stride]);
    Offset = 0;        
    
    if(n >= N)
    {
        do
        {
            Offset += Jump;
            n -= 2*N;
        }while(n >= N);
    }
    else
    {
        while(n < -N)
        {
            Offset -= Jump;
            n += 2*N;
        }
    }
    
    if(n >= 0)
        return Src[Stride*n] + Offset;
    else
        return 3*Src[0] - Src[Stride] - Src[Stride*(-1 - n)] + Offset;
}


static float AsymwExtension(const float *Src, int Stride, int N, int n)
{
    float Jump, Offset;
    
    
    /* Use simple formulas for -N < n < 2*N - 1 */
    if(0 <= n)
    {
        if(n < N)
            return Src[Stride*n];        
        else if(n < 2*N - 1)
            return 2*Src[Stride*(N - 1)] - Src[Stride*(2*(N - 1) - n)];
    }
    else if(-N < n)
        return 2*Src[0] - Src[Stride*(-n)];
    
    /* N == 1 is a special case */
    if(N == 1)
        return Src[0];
    
    /* General formula for extension at an arbitrary n */
    Jump = 2*(Src[Stride*(N - 1)] - Src[0]);
    Offset = 0;        
    
    if(n >= N)
    {
        do
        {
            Offset += Jump;
            n -= 2*(N - 1);
        }while(n >= N);
    }
    else
    {
        while(n <= -N)
        {
            Offset -= Jump;
            n += 2*(N - 1);
        }
    }
    
    if(n >= 0)
        return Src[Stride*n] + Offset;
    else
        return 2*Src[0] - Src[Stride*(-n)] + Offset;
}


boundaryext GetBoundaryExt(const char *Boundary)
{
    if(!strcmp(Boundary, "zpd") || !strcmp(Boundary, "zero"))
        return ZeroPaddedExtension;
    else if(!strcmp(Boundary, "sp0") || !strcmp(Boundary, "const"))
        return ConstantExtension;
    else if(!strcmp(Boundary, "sp1") || !strcmp(Boundary, "linear"))
        return LinearExtension;
    else if(!strcmp(Boundary, "per") || !strcmp(Boundary, "periodic"))
        return PeriodicExtension;
    else if(!strcmp(Boundary, "sym") 
        || !strcmp(Boundary, "symh") || !strcmp(Boundary, "hsym"))
        return SymhExtension;
    else if(!strcmp(Boundary, "symw") || !strcmp(Boundary, "wsym"))
        return SymwExtension;
    else if(!strcmp(Boundary, "asym") 
        || !strcmp(Boundary, "asymh") || !strcmp(Boundary, "hasym"))
        return AsymhExtension;
    else if(!strcmp(Boundary, "asymw") || !strcmp(Boundary, "wasym"))
        return AsymwExtension;
    else
    {
        ErrorMessage("Unknown boundary extension \"%s\".\n", Boundary);
        return NULL;
    }
}