快速Gabor滤波(~1 ms) & python ctype 踩坑经验

原创已于 2024-09-20 14:20:58 修改 · 667 阅读

10 ·

CC 4.0 BY-SA版权

文章标签：

#python #开发语言 #c++ #图像处理 #信号处理

于 2024-09-19 15:14:05 首次发布

最重要的更新：CSDN资源已经上传，想跑一下的可以直接下载！

链接：https://download.csdn.net/download/Zhang__zj__/89771265?spm=1001.2014.3001.5501

然后可能之后会小修小补一点代码的bug（有时候表现不是很好），不过可能不会及时在CSDN上更新（但是github大概率会），所以如果你有issue或者pull request，欢迎来github，我不太常用CSDN，只是ctype吃了很多史才想着来CSDN发一个（还有github我的代码也没人看来引流）。

# 一些更新内容

update: CSDN上同步了和github上一样的代码，之后估计不会有较大的变化，所以就这样了。主要是更新了虚部的计算（不然出的图不好）顺便修了一些小bug。

艹，所有的kernel都写成kernal了……fix typo以后再说吧

update:更新了一点算三角函数的小bug，顺便更新了虚部的计算（因为不算虚部效果确实不太好）

虚部的计算方法放在最后了，在资源给的代码里稍微改改就行。因为可能还有bug于是没有更新资源的内容，（不然还要再打包再改）想要开箱即用在github上直接下载吧。

建议用Frequency = np.float32(np.pi**2/2)，同时注意传到C++代码里面的元素最大为1（我也不知道为啥），否则会出问题！

# 代码介绍

自己的CV代码想要做到realtime，需要一个巨快的Gabor滤波（而且我们需要好多次，效率一定要非常高），但是网上没有现成的轮子（大部分是基于卷积的，速度比较慢），于是下定决心自己写一个（而且也不难）。

实现基于J. Kim, S. Um and D. Min, "Fast 2D Complex Gabor Filter With Kernel Decomposition," in IEEE Transactions on Image Processing, vol. 27, no. 4, pp. 1713-1722, April 2018。他们自称是速度最快的方法，不过我们只算了实数部分，虚数部分没算，因为大部分纯卷积的也就只算实数部分。（update：我们后来加入了虚数部分因为纯实数效果不太好。）

# 代码可用性

代码开源在github上，链接放在下面。

如果看的人多了之后开一个gitee或者什么国内的应该也行？CSDN这个B资源还要审核，总之求大家给点star吧（你GitHub - Zhang-Zhaoji/FastGabor_Python: Based on J. Kim, S. Um and D. Min, "Fast 2D Complex Gabor Filter With Kernel Decomposition," in IEEE Transactions on Image Processing, vol. 27, no. 4, pp. 1713-1722, April 2018,速度大概是：

对于一个 120x160 的矩阵，频率为 8 像素，核大小为 17x17（等效），我们的代码的运行时间为0.0010004043579101562 秒 （只看gabor滤波，不看什么plt之类的）

# 怎么用

python提供了原生的ctypes，基本上就是个C的接口，可以把我的python代码文件打包成一个函数直接调用，或者稍微改改也行。只要dll在旁边就能直接调，不用什么pybind之类的。

# 踩坑经验

不管是中英互联网大部分人都没有搞清楚ctypes还是比较恶心的，不支持C++的大部分内容，包括iostream/string/vector/...，所以，如果你想进一步优化我的代码，注意不要想着找C++的工具。

并行计算听起来很美妙，但是实操貌似会稍微慢一点。

-O3编译器优化非常有效，大概提高了20倍运行速度，赞美机神！

截图留档

# （可能）常见的问题

运行python文件，但是没有plt的输出：

需要输入的image格式大概是np.float32(image)，还必须是二维数组。

明明就在同一个文件夹下，却说找不到dll

这个问题可能是你改了c++的代码？我遇到的大部分情况都是加入了C不兼容的（因此ctype也不兼容的）C++特性代码。如果不是，也许可以试试把写绝对的地址。

## 虚部的计算方法（已经在新代码中更改，如果用的是更新之后的可以不管）

// define IF
        float** IF = (float**)malloc(M * sizeof(float*));
        for (int i = 0; i < M; ++i) {
            IF[i] = (float*)malloc(N * sizeof(float));
        }

// ... 中间其他的计算代码

IF[_y][x] += sin_ws_y[_y] * S_sigma[std::abs(__kernel)] *\
                            (cos_s[std::abs(__kernel)] * RJ[mirror_index(_y,__kernel,M)][x] \
                           + sin_s[std::abs(__kernel)] * IJ[mirror_index(_y,__kernel,M)][x] * (1?__kernel>0:-1)) \
                              - cos_ws_y[_y] * S_sigma[std::abs(__kernel)] *\
                            (sin_s[std::abs(__kernel)] * RJ[mirror_index(_y,__kernel,M)][x] * (1?__kernel>0:-1)\
                           - cos_s[std::abs(__kernel)] * IJ[mirror_index(_y,__kernel,M)][x])  ;

// 其他计算代码

// write answer
        for (int i = 0; i < M; ++i) {
            for (int j = 0; j < N; ++j) {
                // send results towards original matrix
                image[i * N + j] = std::sqrt(F[i][j] * F[i][j] + IF[i][j] * IF[i][j]);
            }
        }
# 释放内存
       for (int i = 0; i < M; ++i) {
        free(IF[i]);
        }
        free(IF);
    }