tvregopt.h

Go to the documentation of this file.

#ifndef _TVREGOPT_H_
#define _TVREGOPT_H_

#if defined(TVREG_NONGAUSSIAN) || (defined(TVREG_INPAINT) && defined(TVREG_DECONV))
#define TVREG_USEZ
#endif

#ifdef TVREG_DECONV
#include <fftw3.h>
#endif
#include "tvreg.h"

#define ALGSTRING_SIZE  128

#define _TVREG_CONCAT(A,B)    A ## B

#ifdef NUM_SINGLE
#define FFT(S)      _TVREG_CONCAT(fftwf_,S)
#else
#define FFT(S)      _TVREG_CONCAT(fftw_,S)
#endif


/* Internal type definitions */

typedef struct
{
    num x;      
    num y;      
} numvec2;

typedef num numcomplex[2];

typedef enum {
    NOISEMODEL_L2,
    NOISEMODEL_L1,
    NOISEMODEL_POISSON
} noisemodel;

struct tag_tvregopt
{
    num Lambda;
    const num *VaryingLambda;
    int LambdaWidth;
    int LambdaHeight;
    const num *Kernel;
    int KernelWidth;
    int KernelHeight;
    num Tol;
    num Gamma1;
    num Gamma2;
    int MaxIter;
    noisemodel NoiseModel;
    int (*PlotFun)(int, int, num, const num*, int, int, int, void*);
    void *PlotParam;
    char *AlgString;
};

typedef struct tag_tvregsolver
{
    num *u;                     
    const num *f;               
    numvec2 *d;                 
    numvec2 *dtilde;            
    num *Ku;                    
    num fNorm;                  
    num Alpha;                  
    int Width;                  
    int Height;                 
    int PadWidth;               
    int PadHeight;              
    int NumChannels;            
    tvregopt Opt;               
    int UseZ;                   
#ifdef TVREG_USEZ
    num *z;                     
    num *ztilde;                
#endif
    
#ifdef TVREG_DECONV
    num *A, *B;                 
    num *ATrans, *BTrans;       
    num *DenomTrans;            
    num *KernelTrans;           
    FFT(plan) TransformA;       
    FFT(plan) TransformB;       
    FFT(plan) InvTransformA;    
    FFT(plan) InvTransformB;    
#endif
} tvregsolver;

typedef num (*usolver)(tvregsolver*);
typedef void (*zsolver)(tvregsolver*);

tvregopt TvRegDefaultOpt = {TVREGOPT_DEFAULT_LAMBDA, NULL, 0, 0, NULL, 0, 0,
    (num)(TVREGOPT_DEFAULT_TOL), TVREGOPT_DEFAULT_GAMMA1,
    TVREGOPT_DEFAULT_GAMMA2, TVREGOPT_DEFAULT_MAXITER, NOISEMODEL_L2,
    TvRestoreSimplePlot, NULL, NULL};

static int TvRestoreChooseAlgorithm(int *UseZ, int *DeconvFlag, int *DctFlag,
    usolver *USolveFun, zsolver *ZSolveFun, const tvregopt *Opt);


/* If GNU C language extensions are available, apply the "unused" attribute
   to avoid warnings.  TvRestoreSimplePlot is a plotting callback function
   for TvRestore, so the unused arguments are indeed required. */
int TvRestoreSimplePlot(int State, int Iter, num Delta,
    ATTRIBUTE_UNUSED const num *u,
    ATTRIBUTE_UNUSED int Width,
    ATTRIBUTE_UNUSED int Height,
    ATTRIBUTE_UNUSED int NumChannels,
    ATTRIBUTE_UNUSED void *Param)
{
    switch(State)
    {
    case 0: /* TvRestore is running */
        /* We print to stderr so that messages are displayed on the console
           immediately, during the TvRestore computation.  If we use stdout,
           messages might be buffered and not displayed until after TvRestore
           completes, which would defeat the point of having this real-time
           plot callback. */
        fprintf(stderr, "   Iteration %4d     Delta %7.4f\r", Iter, Delta);
        break;
    case 1: /* Converged successfully */
        fprintf(stderr, "Converged in %d iterations.           \n", Iter);
        break;
    case 2: /* Maximum iterations exceeded */
        fprintf(stderr, "Maximum number of iterations exceeded.\n");
        break;
    }
    return 1;
}


tvregopt *TvRegNewOpt()
{
    tvregopt *Opt;
        
    if((Opt = (tvregopt *)Malloc(sizeof(tvregopt))))
        *Opt = TvRegDefaultOpt;
    
    if(!(Opt->AlgString = (char *)Malloc(sizeof(char)*ALGSTRING_SIZE)))
    {
        Free(Opt);
        return NULL;
    }
    
    return Opt;
}


void TvRegFreeOpt(tvregopt *Opt)
{
    if(Opt)
    {
        if(Opt->AlgString)
            Free(Opt->AlgString);
        Free(Opt);
    }
}


void TvRegSetLambda(tvregopt *Opt, num Lambda)
{
    if(Opt)
        Opt->Lambda = Lambda;
}


void TvRegSetVaryingLambda(tvregopt *Opt,
    const num *VaryingLambda, int LambdaWidth, int LambdaHeight)
{
    if(Opt)
    {
        Opt->VaryingLambda = VaryingLambda;
        Opt->LambdaWidth = LambdaWidth;
        Opt->LambdaHeight = LambdaHeight;
    }
}


void TvRegSetKernel(tvregopt *Opt,
    const num *Kernel, int KernelWidth, int KernelHeight)
{
    if(Opt)
    {
        Opt->Kernel = Kernel;
        Opt->KernelWidth = KernelWidth;
        Opt->KernelHeight = KernelHeight;
    }
}


void TvRegSetTol(tvregopt *Opt, num Tol)
{
    if(Opt)
        Opt->Tol = Tol;
}


void TvRegSetGamma1(tvregopt *Opt, num Gamma1)
{
    if(Opt)
        Opt->Gamma1 = Gamma1;
}


void TvRegSetGamma2(tvregopt *Opt, num Gamma2)
{
    if(Opt)
        Opt->Gamma2 = Gamma2;
}


void TvRegSetMaxIter(tvregopt *Opt, int MaxIter)
{
    if(Opt)
        Opt->MaxIter = MaxIter;
}


int TvRegSetNoiseModel(tvregopt *Opt, const char *NoiseModel)
{
    if(!Opt)
        return 0;
    
    if(!NoiseModel || !strcmp(NoiseModel, "L2") || !strcmp(NoiseModel, "l2")
        || !strcmp(NoiseModel, "Gaussian") || !strcmp(NoiseModel, "gaussian"))
        Opt->NoiseModel = NOISEMODEL_L2;
    else if(!strcmp(NoiseModel, "L1") || !strcmp(NoiseModel, "l1")
        || !strcmp(NoiseModel, "Laplace") || !strcmp(NoiseModel, "laplace")
        || !strcmp(NoiseModel, "Laplacian") || !strcmp(NoiseModel, "laplacian"))
        Opt->NoiseModel = NOISEMODEL_L1;
    else if(!strcmp(NoiseModel, "Poisson") || !strcmp(NoiseModel, "poisson"))
        Opt->NoiseModel = NOISEMODEL_POISSON;
    else
        return 0;
    
    return 1;
}


void TvRegSetPlotFun(tvregopt *Opt,
    int (*PlotFun)(int, int, num, const num*, int, int, int, void*),
    void *PlotParam)
{
    if(Opt)
    {
        Opt->PlotFun = PlotFun;
        Opt->PlotParam = PlotParam;
    }
}


void TvRegPrintOpt(const tvregopt *Opt)
{
    if(!Opt)
        Opt = &TvRegDefaultOpt;
    
    printf("lambda    : ");

    if(!Opt->VaryingLambda)
        printf("%g\n", Opt->Lambda);
    else
        printf("[%d x %d]\n",
            Opt->LambdaWidth, Opt->LambdaHeight);
    
    printf("K         : ");

    if(!Opt->Kernel)
        printf("(identity)\n");
    else
        printf("[%d x %d]\n", Opt->KernelWidth, Opt->KernelHeight);

    printf("tol       : %g\n", (double)Opt->Tol);
    printf("max iter  : %d\n", Opt->MaxIter);
    printf("gamma1    : %g\n", (double)Opt->Gamma1);
    printf("gamma2    : %g\n", (double)Opt->Gamma2);
    printf("noise     : ");

    switch(Opt->NoiseModel)
    {
    case NOISEMODEL_L2:
        printf("L2\n");
        break;
    case NOISEMODEL_L1:
        printf("L1\n");
        break;
    case NOISEMODEL_POISSON:
        printf("Poisson\n");
        break;
    default:
        printf("(invalid)\n");
        break;
    }

    printf("plotting  : ");

    if(Opt->PlotFun == TvRestoreSimplePlot)
        printf("default\n");
    else if(!Opt->PlotFun)
        printf("none\n");
    else
        printf("custom\n");
    
    printf("algorithm : %s\n", TvRegGetAlgorithm(Opt));
}


const char *TvRegGetAlgorithm(const tvregopt *Opt)
{
    static const char *DefaultAlgorithm =
        (char *)"split Bregman (d = grad u) Gauss-Seidel u-solver";
    static const char *Invalid = (char *)"(invalid)";
    usolver USolveFun;
    zsolver ZSolveFun;
    int UseZ, DeconvFlag, DctFlag;
    
    if(!Opt)
        return DefaultAlgorithm;
    
    if(!TvRestoreChooseAlgorithm(&UseZ, &DeconvFlag,
        &DctFlag, &USolveFun, &ZSolveFun, Opt))
        return Invalid;
    
    sprintf(Opt->AlgString, "split Bregman (%s) %s u-solver",
            (UseZ) ?
                "d = grad u, z = Ku" :
                "d = grad u",
            (!DeconvFlag) ?
                "Gauss-Seidel" :
                ((DctFlag) ?
                    "DCT" :
                    "Fourier"));
    return Opt->AlgString;
}

#endif /* _TVREGOPT_H_ */