zsolve_inc.c

Go to the documentation of this file.

static void ZSolveL1(tvregsolver *S);
static void ZSolveL2(tvregsolver *S);
static void ZSolvePoisson(tvregsolver *S);

#ifndef DOXYGEN
#ifndef _FIDELITY
/* Recursively include file 3 times to define different z solvers */
#define _FIDELITY  1
#include __FILE__           /* Define ZSolveL1 */
#define _FIDELITY  2
#include __FILE__           /* Define ZSolveL2 */
#define _FIDELITY  3
#include __FILE__           /* Define ZSolvePoisson */
#else /* if _FIDELITY is defined */

#include "tvregopt.h"

#if _FIDELITY == 1
static void ZSolveL1(tvregsolver *S)
#elif _FIDELITY == 2
static void ZSolveL2(tvregsolver *S)
#elif _FIDELITY == 3
static void ZSolvePoisson(tvregsolver *S)
#endif
{
    num *z = S->z;
    num *ztilde = S->ztilde;
    const num *Ku = S->Ku;
    const num *f = S->f;
    const num *VaryingLambda = S->Opt.VaryingLambda;
    const num Gamma2 = S->Opt.Gamma2;
    const int Width = S->Width;
    const int Height = S->Height;
    const int NumChannels = S->NumChannels;
    const int PadWidth = S->PadWidth;
    const int PadHeight = S->PadHeight;
    
    const long PadJump = ((long)PadWidth)*(PadHeight - Height);
    num znew;
    long k;
    int x, y;

    if(!VaryingLambda) /* Constant fidelity weight */
    {
        const num Beta = S->Opt.Lambda / Gamma2;
        
        for(k = 0; k < NumChannels; k++, Ku += PadJump)
            for(y = 0; y < Height; y++,
                z += Width, ztilde += Width, f += Width, Ku += PadWidth)
            {
                for(x = 0; x < Width; x++)
                {
#                   if _FIDELITY == 1      /* L1 fidelity      */
                        znew = Ku[x] - f[x] + z[x] - ztilde[x];
                
                        if(znew > Beta)
                            znew += f[x] - Beta;
                        else if(znew < -Beta)
                            znew += f[x] + Beta;
                        else
                            znew = f[x];
#                   elif _FIDELITY == 2    /* L2 fidelity      */
                        znew = (Ku[x] + z[x] - ztilde[x] + Beta*f[x])
                            / (1 + Beta);
#                   elif _FIDELITY == 3    /* Poisson fidelity */
                        znew = (Ku[x] + z[x] - ztilde[x] - Beta)/2;
                        znew = znew + (num)sqrt(znew*znew + Beta*f[x]);
#                   endif
                    
                    ztilde[x] += 2*znew - z[x] - Ku[x];
                    z[x] = znew;
                }
            }
    }
    else    /* Spatially varying fidelity weight */
    {
        const num *LambdaPtr;
        num Beta;
        
        for(k = 0; k < NumChannels; k++, Ku += PadJump)
            for(y = 0, LambdaPtr = VaryingLambda; y < Height; y++,
                z += Width, ztilde += Width, f += Width, Ku += PadWidth)
            {
                for(x = 0; x < Width; x++)
                {
                    Beta = LambdaPtr[x] / Gamma2;
                    
#                   if _FIDELITY == 1      /* L1 fidelity      */
                        znew = Ku[x] - f[x] + z[x] - ztilde[x];
                
                        if(znew > Beta)
                            znew += f[x] - Beta;
                        else if(znew < -Beta)
                            znew += f[x] + Beta;
                        else
                            znew = f[x];
#                   elif _FIDELITY == 2    /* L2 fidelity      */
                        znew = (Ku[x] + z[x] - ztilde[x] + Beta*f[x])
                            / (1 + Beta);
#                   elif _FIDELITY == 3    /* Poisson fidelity */
                        znew = (Ku[x] + z[x] - ztilde[x] - Beta)/2;
                        znew = znew + (num)sqrt(znew*znew + Beta*f[x]);
#                   endif
                    
                    ztilde[x] += 2*znew - z[x] - Ku[x];
                    z[x] = znew;
                }
                
                LambdaPtr += Width;
            }
    }
}
#undef _FIDELITY
#endif /* _FIDELITY */
#endif /* DOXYGEN */