在GitHub: Bokjan/LabDIP的Release中可以找到生成好的二进制。

突如其来的stress

这个学期开了一门课叫数字媒体处理技术，课程内容倒是非常丰富，图像、音频、视频都讲了个遍，涉及的内容也非常广。虽说课是这么一直这么上着，但是给人的感觉是听了也就听了，不知道有什么用，怎么用。这倒是不要紧，10月13号（第六周）开始实验课了。原以为像往常的实验一样，这实验也不打紧，看到任务书倒是目瞪口呆。

修改示例程序，从一个图像显示框，改造成两个显示框，并增加一个文本参数输出框，实现类似下图的基本程序界面（可在此基础上进一步优化）。

1. 功能区可以分tab页，按照后续功能添加；
2. 图像显示区域需考虑图像的缩放与自适应显示；
3. 参数输出区，用以显示过程，以及相关统计数据和调试信息，可滚动，可选择，可复制，可清除。
   阅读程序框架，继续采用Windows多线程和OpenMP两种方式，补充实现下述功能。算法需自行实现，不能直接使用OpenCV函数。
4. 采用三阶插值的图像任意角度旋转与缩放
5. 图像的傅立叶变换，并与功能1联动，输入图像经过旋转、缩放后的傅立叶变换结果可在右侧显示
6. 给图像添加高斯噪声
7. 采用采用平滑线性滤波、高斯滤波、维纳滤波三种方法过滤不同参数的高斯噪声

我还什么都不会呢怎么上来就写一大堆东西？所有算法还必须自己实现？有点stress啊。另外为什么钦定了使用MFC？这下倒好，macOS这段时间也别用了，正好体验一下8102年的Windows 10和Visual Studio到底是什么感觉。

MFC相关

CImage

CImage是ATL/MFC中的一个图像类，可以用它来读取图片并且进行处理。考虑到之后处理图片，写代码的效率和执行效率两个方面，可以考虑将CImage包起来。

class CImageWrapper
{
private:
    CImage *img;
public:
    byte *mem;
    bool IsGray;
    int Width, Height, Pitch, BytePP;
    inline CImageWrapper(CImage *img) :
        img(img), Pitch(img->GetPitch()),
        Width(img->GetWidth()), Height(img->GetHeight()), 
        BytePP(img->GetBPP() / 8), IsGray(BytePP == 1),
        mem((byte*)img->GetBits()) { }
    inline void SetPixel(int x, int y, const byte value)
    {
        *(mem + Pitch * y + BytePP * x) = value;
    }
    inline void SetPixel(int x, int y, const byte r, const byte g, const byte b)
    {
        *(mem + Pitch * y + BytePP * x + 0) = r;
        *(mem + Pitch * y + BytePP * x + 1) = g;
        *(mem + Pitch * y + BytePP * x + 2) = b;
    }
    inline void SetPixel(int x, int y, const byte *ptr)
    {
        *(mem + Pitch * y + BytePP * x + 0) = ptr[0];
        *(mem + Pitch * y + BytePP * x + 1) = ptr[1];
        *(mem + Pitch * y + BytePP * x + 2) = ptr[2];
    }
    inline byte* GetPixel(int x, int y)
    {
        return (byte*)(mem + Pitch * y + BytePP * x);
    }
};

在微软的这一套理论中，byte是一个无符号的8位整数。虽然整个实验中并没有写处理灰度图片的代码，也不存在用灰度图片，但是这里还是有一些关于灰度图的方法，懒得改了。那么这样，将一个CImage的指针丢进去就能构造出来一个操作起来比较简便的东西了。

CImage::GetBits得到的是逻辑上一张图片最左上方的像素的起始地址；CImage::GetPitch得到的某处到下一行的同一列的偏移量（每行的长度并不严格等于每一行像素所占的空间，还加了私货）；而CImage::GetBPP则得到的是每个像素所占的二进制位数。其实这个实验有不少保证的，比如都是彩色图像，每个像素均是24比特大小等。那么上面频繁出现的mem + Pitch * y + BytePP * x + N其实就是对应$(x,\ y)$像素$N$通道的内存地址了。多写几个重载呢，使用起来比较简便。

AFX消息机制下的多线程实现

对于很多图像处理算法来说，它们都是可并行性很好的。实验要求中也明确规定需要实现多线程。

创建一个线程

在AFX中，提供了AfxBeginThread函数用来创建一个线程：

CWinThread* AFXAPI AfxBeginThread(AFX_THREADPROC pfnThreadProc, LPVOID pParam,
    int nPriority = THREAD_PRIORITY_NORMAL, UINT nStackSize = 0,
    DWORD dwCreateFlags = 0, LPSECURITY_ATTRIBUTES lpSecurityAttrs = NULL);

参数众多，但是事实上只需要提供前两个没有默认值的参数。第一个参数pfnThreadProc可以猜测一下意思为pointer of function类型的thread procedure，也就是这个线程要执行什么。

typedef UINT (AFX_CDECL *AFX_THREADPROC)(LPVOID);

找到AFX_THREADPROC的定义，嗯，我们需要提供返回类型为UINT，接受一个LPVOID类型参数的函数给AfxBeginThread，而AfxBeginThread需要的pParam也是LPVOID。其实LPVOID就是void*，你需要的东西都统统弄成指针自己扔来扔去的就好了。

那么我们考虑一下将一张图片的像素数量均匀分给各个线程来执行，除了需要知道每个线程处理的起始、结束像素的索引和原始图像的CImage指针以外，还需要什么？没错，有的操作还需要更多的参数，比如旋转需要知道角度，添加高斯噪声需要知道均值和标准差。有时不能in-situ处理，比如那些基于模板的滤波器，那么还需要一个目标图像的CImage。够了吗？不够，还需要一个回调函数，这样通用性更好——我们并不知道新加入的下一个算法需要做什么它想做的事情。那么我们得到了这样的一个东西：

struct ParallelParams
{
    CImage *img;
    int wParam;
    int begin, end;
    void *ctx, *thctx;
    void(*cb)(ParallelParams*);
    ParallelParams(void) :
        wParam(1),
        cb(nullptr),
        img(nullptr), 
        ctx(nullptr), 
        thctx(nullptr) { }
};

约定img为处理完的图片。ctx为某个算法需要的更多参数，其实这样很蠢——继承多好。但是当时脑筋一是没转过来，就不想再改了。cb是回调函数；thctx目前被用来存放ParallelParams数组的指针以便全部处理完毕之后delete掉，在开始多线程处理前它们当然是被new出来而不是在栈上分配的，不然主线程跑完那一段不就没了吗？

之前提到的AFX_THREADPROC，现在有了自己接受的指针对应的东西。那么分配了一个新的线程来执行它，执行完毕之后干嘛呢？

线程函数的结构

UINT Algo::SaltAndPepperNoise(LPVOID _params)
{
    srand(static_cast<unsigned int>(time(nullptr)));
    auto params = (ParallelParams*)_params;
    CImageWrapper img(params->img);
    for (int i = params->begin; i < params->end; ++i)
    {
        int x = i % img.Width;
        int y = i / img.Width;
        if ((rand() / (double)RAND_MAX) <= NOISE_FACTOR)
        {
            byte val = (rand() & 0x1) ? 0 : 255;
            img.SetPixel(x, y, val, val, val);
        }
    }
    PostMessageW(DA->HWnd, WM_USER_EXECUTE_FINISHED, params->wParam, (LPARAM)params);
    return 0;
}

以椒盐噪声为例，我们使用AFX提供的PostMessageW向某个Dialog发送消息。第一个参数为接受消息的对话框的句柄，第二个参数为消息的类型（用户定义的），而后两个参数是消息处理函数所使用的参数，前者为一个整数，后者为一个指针。

注册消息处理函数

首先要定义一个整数作为消息的ID，下面以WM_USER_EXECUTE_FINISHED为例。之后到欲处理该消息的对话框的message map处添加一条ON_MESSAGE，例如：

ON_MESSAGE(WM_USER_EXECUTE_FINISHED, &CDigitalImageProcessingDlg::OnExecuteFinished)

然后我们就可以来编写CDigitalImageProcessingDlg::OnExecuteFinished了。

编写消息处理函数

LRESULT CDigitalImageProcessingDlg::OnExecuteFinished(WPARAM wParam, LPARAM lParam)
{
    static int cnt = 0;
    ++cnt;
    if (DisplayAgent::GetInstance()->GetThreadOption().count == cnt || wParam == 0)
    {
        cnt = 0;
        auto p = (ParallelParams*)lParam;
        if (p->cb != nullptr)
            p->cb(p);
        DA->PrintTimeElapsed();
        DA->PaintCImageToCStatic(((ParallelParams*)lParam)->img, &mPictureControlRight);
        if (p->ctx != nullptr)
            delete p->ctx;
        if (p->thctx != nullptr)
            delete[] p->thctx;
    }
    this->mButtonExecute.EnableWindow(true);
    return LRESULT();
}

基本道理就是每次一个线程结束都将计数器加一，直到到达设定的线程数（全部完成）。此时若有算法定义的回调就执行，接着在消息区打出时间统计信息、在结果区绘制结果，然后将需要释放的内存释放掉。

其他

除了上述两点以外，MFC当然还有很多需要揣摩的地方——各个控件怎么使用？不同的Dialog之间如何优雅地通信？界面的绘制有什么坑点？高分屏支持怎么样？这些都不涉及“设计”相关的东西，也就不再多说了。

算法

加噪

椒盐噪声

椒盐噪声是最简单的。我们设定一个阈值，再生成一个均匀分布的随机数，如果超过阈值则对某像素加噪。加噪时再生成一个随机数，用它的奇偶来决定是加胡椒（黑）还是盐粒（白）。其实具体的代码都已经贴过了，就在上面线程函数的结构部分。

高斯噪声

其实高斯噪声也是利用随机数，只不过是将一个符合用户给定的均值$\mu$和标准差$\sigma$的正态分布的随机数值加到像素上。那么主要难点也只是如何得到一个正态分布的随机数。

如果要自己实现的话，我们可以采用Box-Muller变换。在之前的《Radiom，来自电磁波的熵》一文中其实已经提到过了。Box-Muller变换可将一对在$[0,\ 1]$上均匀分布的随机数转换为标准正态分布：
$$
\begin{aligned}
Z_0&=\sqrt{-2\ln{U_1}}\cos(2\pi U_2) \
Z_1&=\sqrt{-2\ln{U_1}}\sin(2\pi U_2)
\end{aligned}
$$
要转换到$N(\mu,\ \sigma)$，只需$\sigma Z+\mu$即可。

双三次插值

旋转和缩放都要用到插值。上图来自Wikipedia，我觉得非常的形象易懂。在一维上来说，最邻近插值就是选择映射到原图中的点更靠近的某整数点的值作为插值值；线性插值则是取两个点来拟合一个一次多项式；三次插值拟合的就是三次多项式。而拓展到二维空间，则是先在某方向上处理然后将结果套用在另一方向上处理得到最终结果。

由于真正的双三次插值计算过程异常复杂，Wikipedia上给出了一个性质很接近的卷积核：
$$
W(x)=\left{\begin{matrix}
\begin{aligned}
&(a+2)|x|^3-(a+3)|x|^2+1 &&\text{for}\ |x|\leq 1, \
&a|x|^3-5a|x|^2+8a|x|-4a &&\text{for}\ 1<|x|<2,\
&0 &&\text{otherwise}
\end{aligned}
\end{matrix}\right.
$$
其中，常取$a=-0.5$，$x$为采样点距映射到原图的点的距离。使用该卷积核卷积则得插值结果，但是要记得clamp到合法范围避免溢出。

滤波器

实现中各个滤波器均取$3\times 3$大小的模板。

中值滤波器

对于椒盐噪声，适合使用中值滤波。从字面可知中值滤波就是取每个像素的指定大小邻域的值的中值。由于椒盐噪声要么是最小值要么是最大值，很容易就被过滤掉了。

算术均值滤波器

这应该是实现起来最简单的平滑线性滤波器，对邻域取平均值即可。

高斯滤波器

高斯滤波器实际上是使用一个给定标准差$\sigma$的二维高斯（正态）分布卷积核对模板进行卷积来得到结果。由于卷积核只受$\sigma$影响，因此只需生成一次。我们知道二维高斯分布是：
$$
G(x,\ y)=\frac{1}{2\pi\sigma^2}\exp{\left (-\frac{x^2+y^2}{2\sigma^2}\right )}
$$
那用代码写起来也很简单：

static void GetGaussianTemplate(double t[3][3], double stddev)
{
    const int center = 1; // [[0, 1, 2], [0, 1, 2], [0, 1, 2]], center is (1, 1)
    double total = 0;
    static const double PI = acos(-1);
    for (int i = 0; i < 3; ++i)
    {
        double xsq = pow(i - center, 2.0);
        for (int j = 0; j < 3; ++j)
        {
            double ysq = pow(j - center, 2.0);
            double f = 1 / (2.0 * PI * stddev);
            double e = exp(-(xsq + ysq) / (2.0 * stddev * stddev));
            t[i][j] = f * e;
            total += t[i][j];
        }
    }
    for (int i = 0; i < 3; ++i)
        for (int j = 0; j < 3; ++j)
            t[i][j] /= total;
}

要注意，为了使得结果落在$[0,\ 255]$，要进行归一化处理。

维纳滤波器

作为一个自适应的滤波，这个算是做起来最费劲的了。

首先对每个像素的邻域求均值、方差，那么我们总共会得到$w\times h$（$w$宽，$h$高）个均值$\mu$、方差$\sigma^2$。全部计算完毕后，令估计的噪声$\nu^2=\frac{1}{wh}\sum\sigma^2$。有了这些数据，就可以正式进行滤波了。像素$(i,\ j)$的滤波结果为：
$$
\mu{i,j} +\frac{\max(0,\ \sigma^2{i,j}-\nu^2)}{\max(\sigma^2{i,j},\ \nu^2)}\cdot(\text{p}{i,j}-\mu_{i,j})
$$

注意到完全并行化需要分为两个步骤，一是并行计算全部像素的邻域的均值和方差，这之后才能计算估计的噪声，然后再并行计算滤波结果。

傅里叶变换

为了使效果比较明显，我们先将彩色图片转换为灰度（0.299红、0.587绿、0.114蓝）。实现一个蛮力法的二维离散傅里叶变换（2D DFT）其实不难，其实就是将公式转换为代码:
$$
F(u,\ v)=\sum{x=0}^{M-1}\sum{y=0}^{N-1}f(x,\ y)\textrm{e}^{-\textrm{j}2\pi(\frac{ux}{M}+\frac{vy}{N})}
$$
得到的结果是一个复数，我们想得到的是傅里叶谱，即：
$$
\left | F(u,\ v) \right |=\sqrt{R^2(u,\ v)+I^2(u,\ v)}
$$
其中$R$、$I$是实部和虚部。

这样是做了傅里叶变换没错，但是并不能做到肉眼可见，因此要进一步标定处理。考虑取对数$\ln(|F|+1)$，再乘以一个系数（14左右）用来提升亮度，这样得到的结果就很清晰了。如果要进行中心化，则要将$f(x,\ y)$即原始像素的值乘以$(-1)^{x+y}$，代码实现上可以用f=((x+y)&1)?-f:f;来代替。

进行一番标定之后，用课本上的例子测试一下，效果就很好了。

基于OpenMP的并行化处理

以傅里叶变换为例，基于AFX的多线程处理部分是：

for (int i = 0; i < th.count; ++i)
{
    tasks[i].ctx = img;
    AfxBeginThread(Algo::FourierTransform, tasks + i);
}

其中tasks是th.count个元素的ParallelParams数组。

改写为OpenMP形式其实很简单，可以直接利用已有的代码。可以直接使用#pragma预处理器指令来使用OpenMP，那么在循环中直接调用对应的算法函数就可以了。

#pragma omp parallel for num_threads(th.count)
for (int i = 0; i < th.count; ++i)
{
    tasks[i].ctx = img;
    Algo::FourierTransform(tasks + i);
}

通过这种方式，还可以同时利用AFX的消息处理来解决之后的问题，可以说是修改成本很低。但是一定要注意一点，这个循环不能直接写在UI主线程中，而是要新开一个线程在其中执行。若在UI线程中直接使用该循环，OpenMP只会创建th.count-1个新线程，也就是说UI线程也要承担一份计算工作，这样一来若有AFX消息发出，则消息处理方面会出现问题，没法正确地结束一次完整的操作。

GPGPU?

第一次实验给了2周时间来完成。然而第14天还没有到，第二次实验内容就下来了——为第一次实验的算法编写CUDA和OpenCL版本。由于我只有AMD的显卡，CUDA自然是做不了的，那么只能做OpenCL的版本。

当我开始动手时，农企却给我当头一棒。当你点开这个链接来到AMD Developer Central的时候，你会发现根本没有对应的SDK下载。没关系，我们可以谷歌，然而当我点开AMD APP SDK的链接时，换来的却是一个只有标题没有内容的网页。好在万能的GitHub上面能够找到一个第三方版本：GPUOpen-LibrariesAndSDKs/OCL-SDK。没错整个SDK只有363KiB，不要怀疑，这是真的可以用。NVIDIA的开发者支持当然是要友好得多，把CUDA装好，对应的OpenCL也是一并支持了的。

OpenCL Wrapper

考虑到OpenCL使用起来非常具有Khronos特色，必须要把它包装起来。

class CLAgent
{
private:
    static CLAgent clAgent;
    CLAgent(void);
    ~CLAgent(void);

    cl_device_id*       Devices;
    cl_context          Context;
    cl_kernel           Kernel;
    cl_program          Program;
    cl_command_queue    Queue;
public:
    struct Str
    {
        char* s;
        size_t l;
        inline void Release(void)
        {
            delete[] this->s;
        }
    };
    inline static CLAgent* GetInstance(void)
    {
        return &CLAgent::clAgent;
    }
    Str ReadFile(const char *fn);
    bool LoadKernel(const char * kernel_name);
    bool LoadKernel(const char *fn, const char * kernel_name);
    cl_mem CreateMemoryBuffer(size_t size, void *pointer);
    bool ReadBuffer(cl_mem obj, size_t size, void *dst);
    bool SetKernelArg(cl_uint index, size_t size, const void *pointer);
    bool RunKernel(cl_uint work_dim, const size_t *local_work_size, const size_t *global_work_size);
    void Cleanup(void);
};

优化

实验任务书给出的sample中，采用了新分配6个长度为$w\times h$的数组的方法。也就是说，将原始图片拆分成三个通道再处理，结果也是分开的，最终再合并。看起来就像：

for (int y = 0; y < height; y++)
{
    for (int x = 0; x < width; x++)
    {
        pixelR[index] = *(pRealData + pit * y + x * bitCount + 2);
        pixelG[index] = *(pRealData + pit * y + x * bitCount + 1);
        pixelB[index] = *(pRealData + pit * y + x * bitCount);
        index++;
    }
}
MedianFilter_CL(pixelR, pixelIndex, width, height);
MedianFilter_CL(pixelG, pixelIndex, width, height);
MedianFilter_CL(pixelB, pixelIndex, width, height);
index = 0;
for (int y = 0; y < height; y++)
{
    for (int x = 0; x < width; x++)
    {
        *(pRealData + pit*y + x*bitCount + 2) = pixelR[index];
        *(pRealData + pit*y + x*bitCount + 1) = pixelG[index];
        *(pRealData + pit*y + x*bitCount) = pixelB[index];
        index++;
    }
}

内存拷贝还好，毕竟不大，开销不用太在意。但是这里茫茫多的偏移量计算，多少个乘法啊！全都要放在CPU端单线程进行。我们可以考虑像之前的CPU版本一样直接读写CImage格式的原始数据。但是要将这样的数据送入显存，有点麻烦，需要知道数据的起始地址和长度。然而CImage::GetBits得到的是逻辑上第一个像素的位置而不是物理上的第一个像素。因为查看一下CImage::GetPitch的返回值就会知道pitch是一个负数，也就是说GetBits实际上得到的是物理上最后一行的第一个像素的起始位置。那么起始位置可以推出来了，是bits+pitch*(height-1)。总长度倒是简单一些，就是-pitch+height。

当然核函数部分要用起来比较一致，最好将地址变换回逻辑起始位置。也就是p-pitch*(height-1)这里。我没做与分通道拷贝的对比测试，但有同学修改为优化后的版本后，缩放、旋转的性能提升可以达到100%（Ryzen 2700、GTX 960）。

改写的注意事项

TDR修改

Timeout Detection and Recovery是Windows上处理显卡未响应的功能。如果Windows在一定时间（默认值2s）内没有收到显卡的响应，系统将会重置显卡驱动。因此在进行蛮力傅里叶等非常耗时的操作时，很容易会触发TDR使显卡驱动重启，那么结果就有问题了。我们可以将超时时间设的更长来规避这个问题。在注册表编辑器中找到HKEY_LOCAL_MACHINE\SYSTEM\CurrentControlSet\Control\GraphicsDrivers并添加DWORD项TdrDelay设为一个你觉得合适的值即可。

注意核函数的范围

处理图片时，work size设为二维的，那么在核函数中使用get_global_id(0)和get_global_id(1)来取出像素坐标。但是要注意：由于分配 global work size 时是将图片的宽高按 local work size 的倍数向上取整，因此有一些核函数接收到的位置其实位于真正图像之外，应当先判断是否越界。

随机数？

添加噪声都需要随机数，CUDA至少还提供了随机数功能，OpenCL根本没有，怎么办？方法有二，第一个办法是在 CPU 端生成足量的随机数序列，之后拷贝到显存中供 GPU 使用；二是在 GPU 端自行实现一个线性同余法的伪随机数发生器。两个方案各有优劣，方案一实现起来比较方便，但速度较慢也有额外的拷贝开销；方案二不需要拷贝大块内存，但是需要在各个核函数间同步来保证 PRNG 序列的正确性。我是个懒人，我选方案一。

性能对比

CPU 为 Intel i7-4870HQ（4 核 8 线程），GPU 为 AMD Radeon R9 M370X。时间单位为毫秒，计时精度约为 16.67ms。旋转/缩放操作使用的图片尺寸为 2880×1800 像素；高斯滤波器使用图片尺寸为 14400×9000 像素；高斯噪声生成使用图片尺寸为 5760×3600 像素。傅里叶变换使用图片尺寸为 480×320 像素。


在好 CPU 搭配差 GPU 的情况下，测试结果依然具有代表性。当 CPU 线程数少于等于 CPU 物理核数时，多线程加速效果基本上呈线性关系；而超线程时提升比较有限，这是当然的。而采用 OpenCL 的结果，基本上均为 CPU 最好成绩的50%耗时。这说明图像处理这类可并行的计算密集型应用，使用 GPU 的成百上千个处理器处理有着很强的加速效果。至于添加高斯噪声这一操作，GPU 版本慢于 8 线程 CPU 版本的原因，上文已经提到。为使用的测试图片添加高斯噪声，需要生成 6220.8 万个随机数（数据大小为 237.3MiB）并且将它们拷贝到显存，且这个操作在 CPU 是单线程进行。即使把均匀分布转换高斯分布这一复杂运算转移至 GPU 进行，也无法显著消除随机数生成步骤的大量开销。因此这个添加噪声的 OpenCL 实现对比 CPU 实现不占优势。

后记

你们应该看出来我有些地方贴了实验报告了。

其实说实话这两次实验感觉还好，毕竟时间总共加起来也给了5个星期，还是比较充裕的。虽然都是完全没有接触过的内容，但是实际动手也没有想象的那样困难。以后图像处理估计是不会再做了，但是这次玩了这么多花样，包括GPGPU，还真是一种不一样的体验。另外关于MFC，东西是很老旧这个必须要指出，但使用起来好像还是比较直观比较易于理解的，这点确实是让我有所改观，不过之后肯定也不会用就是了。

在Visual Studio里面工作了这么一段时间，体验不能说是没有，但是总体来说我还是给不了好评。当然这和我自身估计也有点关系，毕竟不是使用这套工具链的人。微软的这套工具链，包括微软自己搞的这一大套东西，历史遗留的不是历史遗留的，都给人一种难以把握的感觉……而且很多东西的行为和另一个世界的那堆工具链实在是太不一样了，一下子很难适应。

所以，我现在当然是赶紧切回macOS去做OpenGL啊！