util_deconv.h

Go to the documentation of this file.

#ifndef _UTIL_DECONV_H_
#define _UTIL_DECONV_H_

#include "tvregopt.h"


static void Divergence(num *DivV, int DivWidth, int DivHeight,
    const numvec2 *V, int Width, int Height, int NumChannels)
{
    int x, y, k;
    
    for(k = 0; k < NumChannels; k++)
    {
        /* Top-left corner */
        DivV[0] = V[0].x + V[0].y;
        
        /* Top row, x = 1, ..., Width - 2 */
        for(x = 1; x < Width - 1; x++)
            DivV[x] = V[x].x - V[x - 1].x + V[x].y;
        
        /* Top-right corner */
        DivV[x] = V[x].y;
        DivV += DivWidth;
        V += Width;
        
        for(y = 1; y < Height - 1; y++, DivV += DivWidth, V += Width)
        {
            /* Left edge */
            DivV[0] = V[0].x + V[0].y - V[-Width].y;
            
            /* Interior */
            for(x = 1; x < Width - 1; x++)
                DivV[x] = V[x].x - V[x - 1].x
                    + V[x].y - V[x - Width].y;
            
            /* Top-right corner */
            DivV[x] = V[x].y - V[x - Width].y;
        }
        
        /* Bottom-reft corner */
        DivV[0] = V[0].x;
        
        /* Bottom row, x = 1, ..., Width - 2 */
        for(x = 1; x < Width - 1; x++)
            DivV[x] = V[x].x - V[x - 1].x;
        
        /* Bottom-right corner */
        DivV[x] = 0;
        DivV += ((long)DivWidth)*((long)DivHeight - Height + 1);
        V += Width;
    }
}


static num UUpdate(tvregsolver *S)
{
    num *u = S->u;
    const num *B = S->B;
    const int Width = S->Width;
    const int Height = S->Height;
    const int PadWidth = S->PadWidth;
    const int PadHeight = S->PadHeight;
    const long PadJump = ((long)PadWidth) * (PadHeight - Height);
    num Norm = 0;
    int x, y, k;
    
    for(k = 0; k < S->NumChannels; k++, B += PadJump)
        for(y = 0; y < Height; y++, u += Width, B += PadWidth)
            for(x = 0; x < Width; x++)
            {
                num unew = B[x];
                num Diff = unew - u[x];
                Norm += Diff * Diff;
                u[x] = unew;
            }
    
    return (num)sqrt(Norm) / S->fNorm;
}


static ATTRIBUTE_ALWAYSINLINE int WSymExtension(int N, int i)
{
    while(1)
    {
        if(i < 0)
            i = -i;
        else if(i >= N)
            i = (2*N - 2) - i;
        else
            return i;
    }
}


static ATTRIBUTE_ALWAYSINLINE int PeriodicExtension(int N, int i)
{
    while(1)
    {
        if(i < 0)
            i += N;
        else if(i >= N)
            i -= N;
        else
            return i;
    }
}

#endif /* _UTIL_DECONV_H_ */