The Fourier transform as a diagonalization

One of the benefits of using the Fourier transform of a function is that convolutions become multiplications. This is important when solving a differential equation with its Green's function. If the Green's function comes from a differential operator $D^*D$, where $D$ is a differential operator and $D^*$ is its adjunct, then the Green's function is not singular at the origin, and is continuous. It expands a function space called a reproducing kernel hilbert space, RKHS, and all functions in this space can be written as linear combinations of the Green's function evaluated on one argument, and the solution to the differential equation $D^*D u = y$ would be of that form. OK, don't digrees anymore... to the cheese...

In the Fourier domain we operate on frequencies $\omega$. For example, to attenuate the noise, we decrease the power in the high omegas, which accounts for a convolution (with a Gaussian, for example). If we see this linear operation as a matrix, the convolution operator that has one (let's say) dimensional Gaussians in its rows (in the time/space domain) becomes a diagonal in the Fourier domain.

The page popped up with much to follow on. In particular, I liked this paragraph

The moral of the story is that the Fourier Transform may be thought of as a change of basis.  The Fourier integral projects a function onto the basis functions of a new coordinate system whose basis functions are the complex exponentials.  In this new basis, the convolution operator is diagonal and everything is simple.  The convolution operator acts on each Fourier component independently by multiplying the component by an associated magnitude and phase.

In Matlab

C=[4 1 2 3; 3 4 1 2; 2 3 4 1; 1 2 3 4]

C =

     4     1     2     3
     3     4     1     2
     2     3     4     1
     1     2     3     4

 F=fft(C)

F =

  10.0000            10.0000            10.0000            10.0000        
   2.0000 - 2.0000i  -2.0000 - 2.0000i  -2.0000 + 2.0000i   2.0000 + 2.0000i
   2.0000            -2.0000             2.0000            -2.0000        
   2.0000 + 2.0000i  -2.0000 + 2.0000i  -2.0000 - 2.0000i   2.0000 - 2.0000i

  F*C*F'

ans =

  1.0e+003 *

   4.0000                  0                  0                  0          
        0             0.0640 - 0.0640i        0                  0          
        0                  0             0.0320                  0          
        0                  0                  0             0.0640 + 0.0640i

I am an Analytic Bastard and I will fight Delusional Geometers to death

This post is dedicated to AMG: friend and enemy, mentor and destroyer, wise and fool.

AMG is obsessed to equate generalized functions (as in Swartz distribution theory) to probability distributions (as in measure theory). According to AMG I am an Analytic Bastard, I agree. And this was the least thing I could take from a Delusional Geometer. Therefore I left AMG.

Schwarz distributions are NOT probability distributions

I can say it louder but not clearer.

It doesn't matter that they are both called distributions, sometimes it does happen in mathematics that two different things are similarly called. It doesn't matter how hard you try to make them the same, it doesn't matter how proud you are and how little you think of the people that surround you and that are not Field medalists.

The fact that Strichartz's book shows a bell-like $C^{\infty}$, compactly supported function does not imply it is a distribution. What is more, this bell-like $C^{\infty}$, compactly supported function is clearly stated to belong to the set $\mathcal{D}$, the set of test functions. Therefore it is not a distribution itself, but the objects to which the linear functional (the Swartz distribution) is applied to. If $\varphi \in \mathcal{D}$ then there is a test function. It looks similar to a DENSITY function (the Gaussian) but the density is not the distribution, nor the test function is the (other kind of) distribution.

Furthermore, forcing my brains so as to accept $\varphi \in \mathcal{D}'$ and call it a (Swartz) distribution, then you can't write $\varphi(x)$ outside the integral symbol. It is a functional, which means that it is applied to some $\phi \in \mathcal{D}$, so what makes sense is $\varphi(\phi)$, don't get angry with me because of this, this is a fact. If $\varphi(x)$ were a linear functional, it could be written as $\int \varphi(x) \phi(x) dx$, why on Earth do you say $\varphi(x)$ anyway?

At this point, why the hell do you use $\varphi(x)$ to name a "bell-like" function with $x= \arg \max_y \varphi(y)$?

Then it remains going full retarded and try to apply density estimation methods to image analysis following this logic:

1. David Mumford develops an axiomatic theory that describes images as generalized functions (check)
2. We have methods that work in the density estimation field fairly well (check)
3. Since YOU (and only you) say generalized functions = probability distributions, then our methods must be very powerful in image analysis (FAIL)

FAIL! Because the only supporting argument you have is your pride.

So, let me out!

Denoising by convolution using CuFFT

In this post we will create a CUDA program that removes some noise with simple isotropic diffusion. We are going to diffuse the function (image) according to the heat equation
$$u_t = \Delta u$$
Where $\Delta$ is the laplacian operator and $u_t$ is the derivative with respect time. The solution at time $t$ is
$$u(x) = \int G_t (x-y) u_o(y) dy$$
Where $x,y\in \mathbb{R}^2$, $u_0$ is the initial noisy image and $G_t$ is the bivariate (normalized) gaussian kernel with standard deviation $t$.

Since a convolution is equivalent to a multiplication of the Fourier transforms as below (by the convolution theorem, and transforming back), we can use the CuFFT library to perform these steps.
$$\hat{u}(\omega) = \hat{G_t}(\omega)\hat{u_o}(\omega)$$

With this project, I want to hunt down two birds with one shot. I want to make a managed code application with Windows forms that reads and visualizes a JPG file, adds noise and then blurs it to remove some noise.

Now, you must know that you cannot mix managed code and CUDA. The solution to this is to separate the CUDA computation and the presentation application in two different parts. We will encapsulate the CUDA computations into a DLL and the Windows program in a Windows Forms application. Both will be part of the same Visual Studio solution so everything will work in the final step. The Windows Forms will be the main project within the VS solution.

Under the DLL project, you must set up CUDA libraries. For the compiler to look for CUDA and CUDAutils header files, we must add

$(NVSDKCOMPUTE_ROOT)\C\common\inc;"$(CUDA_PATH)/include";./;../../common/inc;../../../shared/inc

to  project properties -> C/C++ -> General. Then we must add the directory where additional support header files can be found, in this case,

$(NVSDKCOMPUTE_ROOT)\C\common\inc

. This must be done under project properties -> CUDA Runtime API -> Additional Include Directories. In project properties ->Linker -> Additional library directories you must add

 $(CUDA_PATH)/lib/$(PlatformName);../../common/lib/$(PlatformName);$(NVSDKCOMPUTE_ROOT)\C\common\lib\Win32

Then our DLL function takes a matrix of values and performs the Fourier transform of both (previously padded) the input data and the kernel. We can precompute the kernel with Matlab (with fspecial) and copy the values. In this case, for a kernel with $t=1$, we store the values in h_kernel. Then we perform the multiplication and compute the inverse Fourier transform. This code excerpt contains the important lines:

DLL int blur(float *h_Data, int dataH, int dataW, float **d, int *dH, int *dW) {
    float
        *d_Data,
        *d_PaddedData,
        *d_Kernel,
        *d_PaddedKernel;

    fComplex
        *d_DataSpectrum,
        *d_KernelSpectrum;

    cufftHandle
        fftPlanFwd,
        fftPlanInv;

    const int kernelH = 7;
    const int kernelW = 6;
    ...

    float *h_ResultGPU = (float *)malloc(fftH    * fftW * sizeof(float));

    cutilSafeCall( cudaMalloc((void **)&d_Data,   dataH   * dataW   * sizeof(float)) );
    ...

    float h_Kernel[] = {0.0001F,0.0006F,0.0016F,0.0016F,0.0006F,0.0001F,
        0.0009F,0.0070F,0.0190F,0.0190F,0.0070F,0.0009F,
        0.0043F,0.0314F,0.0854F,0.0854F,0.0314F,0.0043F,
        0.0070F,0.0518F,0.1407F,0.1407F,0.0518F,0.0070F,
        0.0043F,0.0314F,0.0854F,0.0854F,0.0314F,0.0043F,
        0.0009F,0.0070F,0.0190F,0.0190F,0.0070F,0.0009F,
        0.0001F,0.0006F,0.0016F,0.0016F,0.0006F,0.0001F};

    ...

    cufftSafeCall( cufftExecR2C(fftPlanFwd, 
         (cufftReal *)d_PaddedKernel, (cufftComplex *)d_KernelSpectrum) );

    cutilSafeCall( cutilDeviceSynchronize() );
   
    cufftSafeCall( cufftExecR2C(fftPlanFwd, 
         (cufftReal *)d_PaddedData, (cufftComplex *)d_DataSpectrum) );
    modulateAndNormalize(d_DataSpectrum, d_KernelSpectrum, fftH, fftW, 1);
    cufftSafeCall( cufftExecC2R(fftPlanInv, 
         (cufftComplex *)d_DataSpectrum, (cufftReal *)d_PaddedData) );

    cutilSafeCall( cutilDeviceSynchronize() );

    cutilSafeCall( cudaMemcpy(h_ResultGPU, 
         d_PaddedData, fftH * fftW * sizeof(float), cudaMemcpyDeviceToHost) );
    ...

    *d=h_ResultGPU;
    *dH=fftH;
    *dW=fftW;

    return 1;
}

I have highlighted the lines where the Fourier transform, the multiplication and the inverse Fourier transform with the induced blur take place.

In the application, we must have a way to call the DLL function. This code loads the DLL, creates a function pointer and casts it to the DLL function.
private: int CallMyDLL(float *f, int h, int w, float **d, int *dh, int *dw)
             {
                 HINSTANCE hGetProcIDDLL = LoadLibrary(L"cudaint.dll");
                 FARPROC lpfnGetProcessID = 
                      GetProcAddress(HMODULE (hGetProcIDDLL),"blur");
                 typedef int (__stdcall * pICFUNC)
                          (float *, int, int, float**, int *, int*);

                 pICFUNC blur;
                 blur = pICFUNC(lpfnGetProcessID);

                 int MyReturnVal = blur(f, h, w, d, dh, dw);

                 FreeLibrary(hGetProcIDDLL);

                 return MyReturnVal;
             }

Then we can call with the pixels as a linearized matrix and the computations will be made within the GPU. Conversions between .Net objects and float matrices are straighforward.

  

I made the source code available here.