usolve_dct_inc.c

Go to the documentation of this file.

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


static void AdjBlurDct(num *ATrans, FFT(plan) TransformA,
    const num *KernelTrans,
    int Width, int Height, int NumChannels, num Alpha)
{
    const long NumPixels = ((long)Width) * ((long)Height);
    long i;
    int k;
    
    /* Compute ATrans = DCT[A] */
    FFT(execute)(TransformA);
    
    /* Compute ATrans = Alpha . KernelTrans . ATrans */
    for(k = 0; k < NumChannels; k++, ATrans += NumPixels)
        for(i = 0; i < NumPixels; i++)
            ATrans[i] = Alpha * KernelTrans[i] * ATrans[i];
}


static int InitDeconvDct(tvregsolver *S)
{
    num *KernelTrans = S->KernelTrans;
    num *DenomTrans = S->DenomTrans;
    num *B = S->B;
    const num *Kernel = S->Opt.Kernel;
    const int KernelWidth = S->Opt.KernelWidth;
    const int KernelHeight = S->Opt.KernelHeight;
    const int Width = S->Width;
    const int Height = S->Height;
    const num Alpha = S->Alpha;
    const long NumPixels = ((long)Width) * ((long)Height);
    const long PadNumPixels = ((long)Width + 1) * ((long)Height + 1);
    FFT(plan) Plan = NULL;
    FFT(r2r_kind) Kind[2];
    long i;
    int x0, y0, x, y, xi, yi, Size[2];
    
    for(i = 0; i < PadNumPixels; i++)
        B[i] = 0;
    
    x0 = -KernelWidth/2;
    y0 = -KernelHeight/2;
    
    /* Pad Kernel to size Width by Height.  If Kernel
       happens to be larger, it is folded. */
    for(y = 0; y < y0 + KernelHeight; y++)
    {
        yi = WSymExtension(Height + 1, y);
    
        for(x = 0; x < x0 + KernelWidth; x++)
        {
            xi = WSymExtension(Width + 1, x);
            B[xi + (Width + 1)*yi] += Kernel[(x - x0) + KernelWidth*(y - y0)];
        }
    }
    
    /* Compute the DCT-I transform of the padded Kernel */
    if(!(Plan = FFT(plan_r2r_2d)(Height + 1, Width + 1, B, KernelTrans,
        FFTW_REDFT00, FFTW_REDFT00, FFTW_ESTIMATE | FFTW_DESTROY_INPUT)))
        return 0;
    
    FFT(execute)(Plan);
    FFT(destroy_plan)(Plan);
    
    /* Cut last row and column from KernelTrans */
    for(y = 1, i = Width; y < Height; y++, i += Width)
        memmove(KernelTrans + i, KernelTrans + i + y, sizeof(num)*Width);
    
    /* Precompute the denominator that will be used in the u-subproblem. */
    for(y = 0, i = 0; y < Height; y++)
        for(x = 0; x < Width; x++, i++)
            DenomTrans[i] =
                (num)(4*NumPixels*(Alpha*KernelTrans[i]*KernelTrans[i]
                + 2*(2 - cos(x*M_PI/Width) - cos(y*M_PI/Height))));

    /* Plan DCT-II transforms */
    Size[1] = Width;
    Size[0] = Height;
    Kind[0] = Kind[1] = FFTW_REDFT10;
    
    if(!(S->TransformA = FFT(plan_many_r2r)(2, Size, S->NumChannels,
        S->A, NULL, 1, NumPixels, S->ATrans, NULL, 1, NumPixels, Kind,
        FFTW_ESTIMATE | FFTW_DESTROY_INPUT))
        || !(S->TransformB = FFT(plan_many_r2r)(2, Size, S->NumChannels,
        S->B, NULL, 1, NumPixels, S->BTrans, NULL, 1, NumPixels, Kind,
        FFTW_ESTIMATE | FFTW_DESTROY_INPUT)))
        return 0;
    
    /* Plan inverse DCT-II transforms (DCT-III) */
    Kind[0] = Kind[1] = FFTW_REDFT01;
    
    if(!(S->InvTransformA = FFT(plan_many_r2r)(2, Size, S->NumChannels,
        S->ATrans, NULL, 1, NumPixels, S->A, NULL, 1, NumPixels, Kind,
        FFTW_ESTIMATE | FFTW_DESTROY_INPUT))
        || !(S->InvTransformB = FFT(plan_many_r2r)(2, Size, S->NumChannels,
        S->BTrans, NULL, 1, NumPixels, S->B, NULL, 1, NumPixels, Kind,
        FFTW_ESTIMATE | FFTW_DESTROY_INPUT)))
        return 0;
    
    /* Compute ATrans = Alpha . KernelTrans . DCT[f] */
    if(!S->UseZ)
    {
        memcpy(S->A, S->f, sizeof(num)*NumPixels*S->NumChannels);
        AdjBlurDct(S->ATrans, S->TransformA,
            KernelTrans, Width, Height, S->NumChannels, Alpha);
    }
    
    S->Ku = S->A;
    return 1;
}


static void UTransSolveDct(num *BTrans, num *B, FFT(plan) TransformB,
    num *ATrans, const numvec2 *dtilde,
    const num *DenomTrans, int Width, int Height, int NumChannels)
{
    const long NumPixels = ((long)Width) * ((long)Height);
    long i;
    int k;
    
    /* Compute B = div(dtilde) */
    Divergence(B, Width, Height, dtilde, Width, Height, NumChannels);
    
    /* Compute BTrans = DCT[B] */
    FFT(execute)(TransformB);

    /* Compute BTrans = ( ATrans - BTrans ) / DenomTrans */
    for(k = 0; k < NumChannels; k++, ATrans += NumPixels, BTrans += NumPixels)
        for(i = 0; i < NumPixels; i++)
            BTrans[i] = (ATrans[i] - BTrans[i]) / DenomTrans[i];
}


static num UDeconvDct(tvregsolver *S)
{
    /* BTrans = ( ATrans - DCT[div(dtilde)] ) / DenomTrans */
    UTransSolveDct(S->BTrans, S->B, S->TransformB, S->ATrans, S->dtilde,
        S->DenomTrans, S->Width, S->Height, S->NumChannels);
    /* B = IDCT[BTrans] */
    FFT(execute)(S->InvTransformB);
    /* Compute ||B - u||, and assign u = B */
    return UUpdate(S);
}


#if defined(TVREG_USEZ) || defined(DOXYGEN)

static num UDeconvDctZ(tvregsolver *S)
{
    num *ATrans = S->ATrans;
    num *BTrans = S->BTrans;
    const num *KernelTrans = S->KernelTrans;
    const int NumChannels = S->NumChannels;
    const long NumPixels = ((long)S->Width) * ((long)S->Height);
    long i;
    int k;
    
    /* Compute ATrans = Alpha . KernelTrans . DCT[ztilde] */
    memcpy(S->A, S->ztilde, sizeof(num)*NumPixels*NumChannels);
    AdjBlurDct(ATrans, S->TransformA, KernelTrans,
        S->Width, S->Height, NumChannels, S->Alpha);
    /* BTrans = ( ATrans - DCT[div(dtilde)] ) / DenomTrans */
    UTransSolveDct(BTrans, S->B, S->TransformB, ATrans, S->dtilde,
        S->DenomTrans, S->Width, S->Height, NumChannels);
    
    /* Compute ATrans = KernelTrans . BTrans */
    for(k = 0; k < NumChannels; k++, ATrans += NumPixels, BTrans += NumPixels)
        for(i = 0; i < NumPixels; i++)
            ATrans[i] = KernelTrans[i] * BTrans[i];
    
    /* A = IDCT[ATrans] = new Ku */
    FFT(execute)(S->InvTransformA);
    /* B = IDCT[BTrans] = new u */
    FFT(execute)(S->InvTransformB);
    
    /* Compute ||B - u||, and assign u = B */
    return UUpdate(S);
}
#endif