STEP BY STEP

发表于 分类于 阅读次数： Valine：
本文字数： 825 阅读时长 ≈ 1 分钟

首先，说明下为什么要用MATLAB。原因是，支持稀疏矩阵的矩阵库比较少，即使支持去学习使用时间也更多。所以，就用MATLAB了。目的是使用MATLAB计算稀疏矩阵的特征向量，再用C++调用MATLAB代码。

问题1：如何定义稀疏矩阵。MATLAB使用sparse函数创建稀疏矩阵。如下调用：S = sparse(i,j,s,m,n)，根据参数创建稀疏矩阵。i和j代表不为0元素的行和列值，s则是对应的值，m和n是矩阵的大小。当然还有其它接口，具体查阅MATLAB文档。

问题2：怎么计算稀疏矩阵的特征向量？MATLAB提供了eigs函数。其实，还有个eig函数，经本人实验，在使用C++调用MATLAB引擎的过程中，该函数失败了。那么，就只能使用eigs函数了。eigs函数有几种形式。下面说下我用到的几种形式。

形式1：d = eigs(A)，返回的是稀疏矩阵A的特征值（所有的特征值构成一个向量d，本人的理解是没有保证d是有序的）。

形式2：d = eigs(A,k)，返回的是最大的k个特征值，文档也没有说是有序的，如果设置k为1，那么返回的是最大的特征值。

形式3：[V,D] = eigs(A)，V是特征向量构成的矩阵，每列是一个特征矩阵，D则是一个对角矩阵，对角线上的值是和V对应的特征值。不过，V和D的大小都固定为6。

形式4：[V,D] = eigs(A，k)，则可以返回最大的k个特征向量和对应的特征值，但是也没有排序，需要根据D排序V，以方便使用。

问题3：如何实现返回k个根据特征值排序后的特征向量？这是最有价值的问题了。直接提供代码吧。

1
2
3
[V,D] = eigs(sm,k);
[newD,IX] = sort(diag(D)));;
newV=V(:,IX);

newV就是根据特征值排序后的特征向量，sm当然需要先建立了。
写在最后的话，之所以使用MATLAB计算稀疏矩阵的特征向量，是因为可以用C++调用。所以，关于这篇文章这些知识就够了，也算是我为自己做的整理吧。

发表于 分类于 ， ， 阅读次数： Valine：
本文字数： 4.9k 阅读时长 ≈ 4 分钟

前面提到priority_queue的内部使用了stl的泛型heap算法。现在让我们看看怎么来实现heap。至于heap，就不作仔细介绍了，本质上是用一个数组表示的完全二叉树，并且父节点总是大于（或者小于）子节点的值。我们先来看push_heap。

下图是push_heap的示意图。

从图中可以看到算法的过程是将新加入堆的值（50），层层上挪，直到正确的位置。下面来看，摘录出来的代码。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
template <class RandomAccessIterator>
inline void push_heap(RandomAccessIterator first, RandomAccessIterator last) {
    // 注意，调用该函数时候，新元素位于最后一个位置(last-1)。
    __push_heap_aux(first, last, distance_type(first), value_type(first));
}

template <class RandomAccessIterator, class Distance, class T>
inline void __push_heap_aux(RandomAccessIterator first,
    RandomAccessIterator last, Distance*, T*) {
    __push_heap(first, Distance((last - first) - 1), Distance(0),
        T(*(last - 1)));
    // (last-first)–1代表新元素的索引，0是堆首的索引，*(last - 1)是新加入的值
}

template <class RandomAccessIterator, class Distance, class T>
void __push_heap(RandomAccessIterator first, Distance holeIndex,
    Distance topIndex, T value) {
    Distance parent = (holeIndex - 1) / 2;  // 找出父節點
    while (holeIndex > topIndex  *(first + parent) < value) {
        // 尚未到达顶端，且父节点小于新值
        // 由于以上使用 operator<，可知 STL heap 是max-heap
        *(first + holeIndex) = *(first + parent);   // 令洞值为父值
        holeIndex = parent; // percolate up：调整洞号，向上提升至父节点。
        parent = (holeIndex - 1) / 2;   // 新洞的父节点
    }    // 持续至顶端，或满足 heap 的次序特性为止。
    *(first + holeIndex) = value;   // 令洞值为新值。
}

push_heap的用法是输入迭代器对，并且保证[first,last-1)是最大堆，*(last-1)是新加入的元素。push_heap调用辅助函数push_heap_aux。至于为什么需要这个辅助函数了？应该是为了提取出distance_type和value_type吧，这两个内联函数的定义，可以参考stl源码剖析迭代器的那章。下面来思考真正的实现函数push_heap。这个函数需要新加入元素位置holeIndex和堆首位置topIndex，另外还有保存好的新加入值。算法的过程很简单，就是上溯holeIndex，找到真正满足条件的位置（无法继续上回溯），然后把value放入该位置即可。
pop_heap实际上是一个相反的过程。实现思路是将堆大小加一后，再找出最后一个元素应该放入的位置holeIndex，最后再加入这个值。示意图如下：


下面看看摘录出来的代码，思想类似于push_heap，只需要求出最终的holeIndex。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
template <class RandomAccessIterator>
inline void pop_heap(RandomAccessIterator first, RandomAccessIterator last) {
    __pop_heap_aux(first, last, value_type(first));
}

template <class RandomAccessIterator, class T>
inline void __pop_heap_aux(RandomAccessIterator first,
    RandomAccessIterator last, T*) {
    __pop_heap(first, last - 1, last - 1, T(*(last - 1)), distance_type(first));
    // pop动作的結果为底层容器的第一個元素。因此，首先设定欲调整值为尾值，然后將首值調至 
    // 尾节点（所以以上將迭代器result设为last-1）。然后重整 [first, last-1)，
    // 使之重新成一個合格的 heap。
}

template <class RandomAccessIterator, class T, class Distance>
inline void __pop_heap(RandomAccessIterator first, RandomAccessIterator last,
    RandomAccessIterator result, T value, Distance*) {
    *result = *first; // 設定尾值为首值，于是尾值即是結果，
    // 可由调用底层容器之 pop_back() 取出尾值。
    __adjust_heap(first, Distance(0), Distance(last - first), value);
    // 以上欲重新調整 heap，洞号为 0，欲調整值为value。
}

template <class RandomAccessIterator, class Distance, class T>
void __adjust_heap(RandomAccessIterator first, Distance holeIndex,
    Distance len, T value) {
    Distance topIndex = holeIndex;
    Distance secondChild = 2 * holeIndex + 2;   // 洞节点之右子节点
    while (secondChild < len) {
        // 比较洞节点之左右兩个子值，然后以 secondChild 代表较大子节点。
        if (*(first + secondChild) < *(first + (secondChild - 1)))
            secondChild--;
        // Percolate down：令较大大子值为洞值，再令洞号下移至较大子节点处。
        *(first + holeIndex) = *(first + secondChild);
        holeIndex = secondChild;
        // 找出新洞节点的右子节点
        secondChild = 2 * (secondChild + 1);
    }

    if (secondChild == len) { // 沒有右子节点，只有左子节点
        // Percolate down：令左子值为洞值，再令洞号下移至左子节点处。
        *(first + holeIndex) = *(first + (secondChild - 1));
        holeIndex = secondChild - 1;
    }

    // 將欲调整值填入目前的洞号內。注意，此時肯定滿足次序特性。
    // 依侯捷之见，下面直接改為 *(first + holeIndex) = value; 应该可以。
    __push_heap(first, holeIndex, topIndex, value);
}

类似于push_heap，pop_heap也是调用辅助函数pop_heap_aux。pop_heap_aux调用pop_heap。pop_heap调用adjust_heap调整holeIndex，最终在holeIndex处放入value（原最后一个的值）。关键代码是adjust_heap中的循环。循环的主要意思是将holeIndex不断下放，直到最底层。最后的if语句的意思是，如果最底层有左子节点，而没有右子节点，那么最终位置肯定是这个左子节点。最后一句代码的意思是加入value到holeIndex，由于已经调整完毕，所以一个赋值操作也可以达到要求，参见侯捷注释。
sort_heap就比较简单了，不断将极值移动到末尾，不断pop_heap。

1
2
3
4
5
6
7
8
9
// 以下這個 sort_heap() 不允許指定「大小比較標準」
template <class RandomAccessIterator>
void sort_heap(RandomAccessIterator first, RandomAccessIterator last) {
    // 以下，每執行一次 pop_heap()，極值（在STL heap中為極大值）即被放在尾端。
    // 扣除尾端再執行一次 pop_heap()，次極值又被放在新尾端。一直下去，最後即得
    // 排序結果。
    while (last - first > 1)
        pop_heap(first, last--); // 每執行 pop_heap() 一次，操作範圍即退縮一格。
}

最后要将的是make_heap，即将一个迭代器对里面的内容构造为最大堆。代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
// 將 [first,last) 排列為一個 heap。
template <class RandomAccessIterator>
inline void make_heap(RandomAccessIterator first, RandomAccessIterator last) {
    __make_heap(first, last, value_type(first), distance_type(first));
}

// 以下這組 make_heap() 不允許指定「大小比較標準」。
template <class RandomAccessIterator, class T, class Distance>
void __make_heap(RandomAccessIterator first, RandomAccessIterator last, T*,
    Distance*) {
    if (last - first < 2) return;   // 如果長度為 0 或 1，不必重新排列。
    Distance len = last - first;
    // 找出第一個需要重排的子樹頭部，以 parent 標示出。由於任何葉節點都不需執行 
    // perlocate down，所以有以下計算。parent 命名不佳，名為 holeIndex 更好。
    Distance parent = (len - 2) / 2;

    while (true) {
        // 重排以 parent 為首的子樹。len 是為了讓 __adjust_heap() 判斷操作範圍
        __adjust_heap(first, parent, len, T(*(first + parent)));
        if (parent == 0) return;    // 走完根節點，就結束。
        parent--;                   // （即將重排之子樹的）頭部向前一個節點
    }
}

make_heap中代码的思路也很简单。从原序列的中间位置开始不断调整（调用adjust_heap），每次调整的目的是以当前位置为根的构建一个子堆。至于为什么从中间位置开始就可以了？原因很简单，最底层元素的数目大致就会占了一半了。

发表于 分类于 ， ， 阅读次数： Valine：
本文字数： 1.3k 阅读时长 ≈ 1 分钟

类似于queue，stack也是个简单的适配器。在实现上，stack与queue非常类似。底层可以使用同样的容器，比如默认都采用deque。因此，stack的源码和queue的非常类似。代码过于简单，都可以自解释了。下面是摘录出来的代码，可以顺便编译通过。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
#ifndef _STL_STACK_H
#define _STL_STACK_H

namespace std
{
    template <class T, class Sequence = deque<T> >
    class stack {
        friend bool operator==(const stack, const stack);
        friend bool operator<(const stack, const stack);

    public:
        typedef typename Sequence::value_type value_type;
        typedef typename Sequence::size_type size_type;
        typedef typename Sequence::reference reference;
        typedef typename Sequence::const_reference const_reference;
    protected:
        Sequence c; // 底层容器

    public:
        // 以下完全利用 Sequence c 的操作，完成 stack 的操作。
        bool empty() const { return c.empty(); }
        size_type size() const { return c.size(); }
        reference top() { return c.back(); }
        const_reference top() const { return c.back(); }

        // deque 是两头可进出，stack 是末端进，末端出（所以后进者先出）。
        void push(const value_type x) { c.push_back(x); }
        void pop() { c.pop_back(); }
    };

    template <class T, class Sequence>
    bool operator==(const stack<T, Sequence> x, const stack<T, Sequence> y) {
        return x.c == y.c;
    }

    template <class T, class Sequence>
    bool operator<(const stack<T, Sequence> x, const stack<T, Sequence> y) {
        return x.c < y.c;
    }

};

#endif

代码非常类似于queue，除了pop弹出的是末尾元素，top也是返回末尾元素之外。

发表于 分类于 ， ， 阅读次数： Valine：
本文字数： 1.5k 阅读时长 ≈ 1 分钟

接着上面一篇，这次要讲的是另一个适配器容器queue。queue的操作很简单，先进先出。现在能够使用底层容器来实现，所以也非常简单。总之，这就是个简单的适配器吧。queue需要两个模版参数，类型和底层容器，默认的底层容器是deque。摘录出的代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
#ifndef _STL_QUEUE_H
#define _STL_QUEUE_H
#include <deque>

namespace std
{
    template <class T, class Sequence = deque<T> >
    class queue {
        friend bool operator== (const queue x, const queue y);
        friend bool operator< (const queue x, const queue y);

    public:
        typedef typename Sequence::value_type value_type;
        typedef typename Sequence::size_type size_type;
        typedef typename Sequence::reference reference;
        typedef typename Sequence::const_reference const_reference;

    protected:
        Sequence c; // 底层容器

    public:
        // 以下完全利用 Sequence c 的操作，完成 queue 的操作。
        bool empty() const { return c.empty(); }
        size_type size() const { return c.size(); }
        reference front() { return c.front(); }
        const_reference front() const { return c.front(); }
        reference back() { return c.back(); }
        const_reference back() const { return c.back(); }

        // deque 是两头可进出，queue 是末端进，前端出（所以先进者先出）。
        void push(const value_type x) { c.push_back(x); }
        void pop() { c.pop_front(); }
    };

    template <class T, class Sequence>
    bool operator==(const queue<T, Sequence> x, const queue<T, Sequence> y) {
        return x.c == y.c;
    }

    template <class T, class Sequence>
    bool operator<(const queue<T, Sequence> x, const queue<T, Sequence> y) {
        return x.c < y.c;
    }

};

#endif

总之，代码非常简洁。分析下代码，typedef typename Sequence::value_type value_type;这其实，就是声明嵌套模版类型。还有必须注意的是，声明嵌套模版类型，必须加typename关键字，否则会将其视作定义。

发表于 分类于 ， ， 阅读次数： Valine：
本文字数： 2.2k 阅读时长 ≈ 2 分钟

暑假时候好好阅读过，stl源码剖析这本书，之后又看了c++标准程序库，effective c++，more
effecttive c++。虽然每次都如醍醐灌顶般，明白了很多知识概念，有种拨开迷雾的感觉，但是终究抵不过
遗忘，所以了还是随便写写吧。
首先从实现代码差不多最短的优先队列开始，虽然简单的代码后面依赖很庞大的东西，比如泛型
迭代器，还有泛型二叉堆算法。我是这样理解优先队列的。首先，数据是存在一个泛型序列容器(线性容器)，
这个可以指定，默认是vector，当然指定list和deque也是可以的。另外还需要给数据类型提供个小于比较仿
函数，有些类型有默认的less泛型仿函数。所以了，对模版比较熟悉的同学，很自然就猜到了，priority_queue
的外观->一个模版（类型，线性容器，less仿函数）。
优先队列支持哪些操作了？什么empty，size，top，直接转交给底层容器就行了，所谓的适配器模式。
只需要思考，如何实现push和pop操作。对算法熟悉的同学，都知道底层数据的组织是用二叉堆的，而且是最大
堆。问题是，我们还需不需要实现个堆操作了。不需要了，stl做了所有的事情。push_heap就是将输入范围的最
后一个元素作为在原堆后面加入的新元素，再重构堆。pop_heap则是将堆首元素交换到堆末，再重构缩小之后的堆。
至于二叉堆的原理，就不仔细介绍了。其实也比较简单，如果不去思考一些复杂度的证明。
到最后，priority_queue的代码可以很清爽的出来了。尤其注意下面的那些typedef，可以说是模版
trick的常见手法，习惯了就好。这种写法还有个用法是为了使内嵌模版声明可见，记得effect c++提到过。以下
源码基本来自stl源码剖析随书奉送代码，做了些修改。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
#include <algorithm>

template <typename T, typename Sequence = vector<T>,
    typename Compare = less<typename Sequence::value_type >>
class  priority_queue
{
public:
    typedef typename Sequence::value_type value_type;
    typedef typename Sequence::size_type size_type;
    typedef typename Sequence::reference reference;
    typedef typename Sequence::const_reference const_reference; //为了使内嵌模版声明可见

protected:
    Sequence c; // 底层容器
    Compare comp;   // 元素大小比较仿函数

    priority_queue() : c() {}
    explicit priority_queue(const Compare x) : c(), comp(x) {}

    template <class InputIterator>
    priority_queue(InputIterator first, InputIterator last, const Compare x)
        : c(first, last), comp(x) {
        make_heap(c.begin(), c.end(), comp);
    }

    template <class InputIterator>
    priority_queue(InputIterator first, InputIterator last)
        : c(first, last) {
        make_heap(c.begin(), c.end(), comp);
    }

    bool empty() const { return c.empty(); }
    size_type size() const { return c.size(); }
    const_reference top() const { return c.front(); }

    void push(const value_type x)
    {
        try
        {
            c.push_back(x);
            push_heap(c.begin(), c.end(), comp);
        }
        catch (...)
        {
            c.clear();
            throw;
        }
    }

    void pop()
    {
        try
        {
            pop_heap(c.begin(), c.end(), comp);
            c.pop_back();
        }
        catch (...)
        {
            c.clear();
            throw;
        }
    }
};

再啰嗦几句，刚看源码剖析的时候，对某些模版写法很头疼，但是看完书之后，另外看了几本其它的书，
就对这些写法习以为常了，觉得不这么写才怪了。有点艺术的感觉。

发表于 分类于 阅读次数： Valine：
本文字数： 454 阅读时长 ≈ 1 分钟

这里统一考虑二维直线和三维平面的情况。
假设法线（二维直线和三维平面的法线）为n，其上的任意一点为p0。那么，
可以得到方程为：(p - p0) * n = 0。展开得到p * n = p0 * n。令p0 * n = d，由此可以得到
方程可以简单的表示为p * n = d，即(x, y, z) * (nx ,ny, nz) = d，二维情况直线的为(x, y) * (nx, ny) = d。
至此，我们可以方便的用(n,d)代表任意的二维直线或者三维平面了。
下面，开始应用这个表示。
比如，求点到直线或者平面的距离，然后求投影点，再求对称点。
最关键的是如何方便求出距离。假设我们要求q到直线的距离。我们现在的直线方程为
p * n = d。那么经过q与当前直线的平行直线的法线也为n。由于q在经过q的平行直线上面，有q * n = d’。
有符号距离为 d - d’ = d - q * n。通过这个式子，也可以知道d其实就是原点到直线或者平面的有符号距离。
那么投影点q’ = q + (d - d’) = q + (d - q * n)。
对称点则为 q’’ = q’ + (d - d’) = q + 2 * (d - q * n)。
利用这种表示方法，基本上可以直接给出答案。

发表于 分类于 ， 阅读次数： Valine：
本文字数： 286 阅读时长 ≈ 1 分钟

一直以来都默认删除指针前，都得判断下是否是空指针。直到最近在看More Effective C++
的时候，作者明确指出了C++在语言层次上保证删除空指针是安全的，我才意识到这个问题。
我用vs2013实践了下，发现没有运行错误。所以了，这大概是从学习C语言起遗留下来的思维
定势吧。以前用C调用malloc申请内存的时候，对应的是free释放，free肯定是不能释放NULL指针的，所以
自然而然就觉得delete也不能处理NULL指针。
虽然说，这件事情无伤大雅，但是确实没必要在delete之前判断下是否为NULL了，毕竟自己
判断和编译器判断都是一样的结果，也不存在什么效率问题。

发表于 分类于 ， 阅读次数： Valine：
本文字数： 978 阅读时长 ≈ 1 分钟

首先，当然是安装好vs2013。第二步是配置最新的cuda6.5版本，这个有相关的教程。
这里要讲的是如何生成相应的opencv的gpu版本。
首先，当然是去opencv官方网站下载opencv，然后解压到指定的目录，最好不要出现
中文路径，因为后面需要用到cmake。2.4.10的下载地址，已经不好找了。最新的是3.0beta版本，
由于模块分类完全变了，所以还是别去吃螃蟹了。
下载好解压之后，会出现两个子目录，build和source。如果，你只使用cpu版本，那么build
里面的内容就满足你的要求了。否则，打开source目录，你会发现有一个CMakeLists.txt文件。
那么，肯定是需要使用cmake生成相应的工程了。
第二步，下载cmake，安装好cmake。
第三步，用cmake生成工程。操作方法：假如opencv在f盘根目录下面，那么设置好cmake的
source目录和build目录。然后点击configure按钮，选择对应的平台如visual studio 12。注意，现在
记得选择WITH_CUDA，WITH_CUFFT，WITH_CUBLAS等相关的选项。
然后，点击generate按钮就能生成对应的工程。如果，出现配置错误，可以根据cmake的
提示寻找原因。基本上是些环境变量的设置问题，或者是平台不匹配等。
如图所示：


注意build目录不用和source目录设置为同一个，否则可能引起问题，cmake也会出警告。
第四步，就是打开工程生成文件了。用vs2013或者你的使用版本，打开build目录下的OpenCV.sln
文件。最简单的办法，就是批生成->全选->重新生成。
关键的问题出现，如果你只是这样做，可能有几个dll无法编译出来，比如说opencv_gpu2410.dll。
问题在哪里了。我试验了好几次，每次几个小时，都是这样的结果。顿时觉得很郁闷，
只能去分析vs的输出内容了。突然发现有个错误提示：大致内容是NCV.cu文件中的std::max没有
定义。这个cu文件在..\opencv\sources\modules\gpu\src\nvidia\core下面。
你最后需要做的就是打开这个文件，然后包含algorithm头文件就行了，然后再重新编译吧。
PS：我只生成过32平台下的gpu版本了，64位的就没有去尝试了。

发表于 分类于 ， 阅读次数： Valine：
本文字数： 450 阅读时长 ≈ 1 分钟

经测试，第一次gpu调用，无论是用opencv的gpu模块或者cuda都会比较耗时，可能将近1s。
额，这也不是我一个人碰到的问题，确实有这回事。stackoverflow上面有个帖子就是关于这个的。
帖子地址：关于速度慢，关于解决办法。
如果，看了上面的帖子话，这篇文章也不用看下去了。因为我要讲的就是这件事情。
我的观点也是，在第一次调用之前先建立cuda context。为什么需要做这样的事情了？额，假如，
你只调用一次gpu处理，总不能太慢了吧。那样就看不到速度了，所以先调用次垃圾操作初始化cuda环境。
方法是调用cudaFree(0)，在这之前最好调用先调用cudaSetDevice(0)。记住包含cuda的头文件，
如果只有opencv的头文件，这2个函数是找不到的。还有，建立的工程是cuda runtime模版。
在我的L0Smooth代码里，这样的处理之后，初次L0Smooth调用能减少1s左右，从3s变成了2s。
其余的耗时，基本都在gpu版本的dft和idft，还得继续寻找加快速度的办法。

发表于 分类于 ， 阅读次数： Valine：
本文字数： 1.6k 阅读时长 ≈ 1 分钟

最近打算使用gpu优化L0Smooth的代码，但是不熟悉opencv的gpu部分的使用，
直接用cuda觉得太麻烦，还是继续试试opencv吧。
首先，从网站上下载下来的默认版本是不支持gpu的，必须下载代码下来，用cmake生成工程，注意
选择支持cuda等选项，具体参考教程。这样生成的版本才能支持gpu运算。
那么如何测试自己的opencv是否支持gpu计算了，以及自己的硬件是否符合要求？
使用这句代码打印cuda设备个数：
printf(“Device Num:%d\n”, cv::gpu::getCudaEnabledDeviceCount());
如果，个数大于0，那么说明你的显卡是支持cuda的，并且你的opencv版本支持gpu运算了。
下面是我测试成功的gpumat实现的dft和idft函数，输入和输出的都是cpu上的mat。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
void Dft(cv::Mat in, cv::Mat out)
{
    cv::gpu::GpuMat gpuIn0(in.size(), CV_32FC1);
    gpuIn0.upload(in);

    std::vector<cv::gpu::GpuMat> planes;

    planes.push_back(gpuIn0);

    cv::Mat zero = cv::Mat::zeros(in.size(), CV_32FC1);
    cv::gpu::GpuMat gpuZero(in.size(), CV_32FC1);
    gpuZero.upload(zero);

    planes.push_back(gpuZero);

    cv::gpu::GpuMat gpuIn(in.size(), CV_32FC2);
    cv::gpu::merge(planes, gpuIn);

    cv::gpu::GpuMat gpuOut(gpuIn.size(), CV_32FC2);

    cv::gpu::dft(gpuIn, gpuOut, gpuIn.size(), 0);

    out.create(in.size(), CV_32FC2);
    gpuOut.download(out);
}

void IDft(cv::Mat in, cv::Mat out)
{
    cv::gpu::GpuMat gpuIn(in.size(), CV_32FC2);
    cv::gpu::GpuMat gpuOut(in.size(), CV_32FC2);

    gpuIn.upload(in);
    cv::gpu::dft(gpuIn, gpuOut, gpuIn.size(), cv::DFT_INVERSE);

    cv::gpu::GpuMat splitter[2];
    cv::gpu::split(gpuOut, splitter);

    out.create(in.size(), CV_32FC1);

    double minV, maxV;
    cv::gpu::minMax(splitter[0], minV, maxV);
    splitter[0].convertTo(splitter[0], CV_32F, 255.0 / maxV);
    splitter[0].download(out);
}

从代码上，可以看出，要尽可能把运算放到gpu上去，包括很多辅助运算，比如说merge,convert等可以加快速度。
这个版本的代码，dft的输入和输出都是双通道的，也就是一个通道实数，一个是复数的矩阵。通过merge生成和split分离。
这两个函数在分离通道单独处理时候非常有用，比如说可以分离彩色图像的三个通道，单独进行dft处理。最后再将结果合并。