本文首发于公众号【Opencv视觉实践】，各种有趣的小操作带你进阶opencv
之前写过一篇
【opencv】带你再学一遍直方图​mp.weixin.qq.com
里面的内容可以看下图。所以今天还要再再再学一个直方图的API：直方图反投影。
直方图反向投影是干啥的呢，它用于图像分割或在图像中查找感兴趣的对象。用什么来查找呢？用直方图。
直方图在一定程度上可以反应图像的特征，我们截取一个有固定特征的样例，比如草地，然后计算该块草地的直方图，然后用这个直方图去和整幅图像的直方图做对比，根据一定的判断条件，就能得出相似的即为草地。
看着就像是语义分割，其实一定意义上这就是语义分割，这不过直方图反向分割的依据是人为计算的（直方图）,后者分割的依据是靠在神经网络中学习得来的。
我们先看一下opencv直方图反向计算得API：
void cv::calcBackProject(
  const Mat *  images, //要进行投影的输入图像的地址，注意该API要求输入的是地址
  int          nimages,//输入图像的数目
  const int *  channels,//要进行投影的通道数
  InputArray   hist,//样本得直方图
  OutputArray  backProject,//输出得反向投影，为Mat类型
  const float ** ranges, //输入直方图得特征空间的取值范围
  double         scale = 1,
  bool          uniform = true
)
一：该API得实现原理是什么呢
假设我们现在有一个四行四列得灰度图，它得灰度值如下图：
说这幅图有什么特征呢？直观上看类似于一个边角，但这是直观上，怎么表示出来呢？深度学习是靠神经网络黑箱计算出来得，我们可以用直方图。
那我们就计算这幅灰度图得直方图，如果以组距为1计算直方图并反向投影到原图，得到得为下图：
可以大概表述一下边角得特征：左下角有6个像素值相同得三角形区域，中间斜向下有四个像素值相同得边界线，以此类推。这就是用直方图得到得边角得特征。
那如果以组距为2计算直方图呢？反向投影后为：
可以看到特征描述得更为广泛了，就像深度学习里，提取更高层次得特征，虽然更为普适，但也会忽略掉一些细节特征。
我们就是拿这个反向投影所表达得特征信息，去和整幅图做对比，来得到特征相似得部分，达到分割得效果。
二：利用反向投影进行语义分割
先看一下我们今天要处理得图片：
我们今天要做得就是将这条公路给提取出来。
1，先读取原图以及样本图，并转换为HSV格式。
//【1】读取图片
Mat srcImage = imread("风景图.jpg");
Mat RoiImage = imread("公路2.png"); 
//【2】转换为HSV图像
Mat HsvImage, RoiImage_HSV;
cvtColor(srcImage, HsvImage, COLOR_BGR2HSV);
cvtColor(RoiImage, RoiImage_HSV, COLOR_BGR2HSV);
为什么转HSV呢？因为HSV表达颜色更为方便区分，我们今天用到得只有前两个通道：H（色调）和S（饱和度），不用V（亮度）。
来看一下我们截取得样本图：
这是我们在公路上随便截取得一块儿样本，可以看到整条公路大概都是这个样子。
2，计算样本图得直方图并进行归一化
//【3】计算公路的直方图
MatND roiHist; //直方图对象
int dims = 2;  //特征数目（直方图维度）
float hranges[] = { 0,180 }; //特征空间的取值范围
float Sranges[] = { 0,256 };
const float *ranges[] = { hranges,Sranges };
int size[] = { 20,32 };  //存放每个维度的直方图的尺寸的数组
int channels[] = {0,1};  //通道数
calcHist(&RoiImage_HSV,1, channels, Mat(), roiHist, dims, size, ranges);
//【4】直方图归一化
normalize(roiHist, roiHist, 0, 255, NORM_MINMAX);
代码第7行，我们计算H和S通道直方图的组距分别为20和32，我们采用更小的组距来抓取更大的特征区域。
为什么要归一化呢，直方图反向投影到原图后，原图各位置表示的是整幅图中等于该点像素值的数量，归一化后就变成概率了。
3，将计算的归一化后的直方图进行反向投影
//【5】反向投影
Mat proImage; //投影输出图像
calcBackProject(&HsvImage, 1, channels, roiHist, proImage, ranges);
反向投影后的原图：
可以大概看出公路的轮廓了吧？
图像腌膜Mask的常规操作你真的信手拈来吗？​mp.weixin.qq.com
4，上一篇介绍了掩码操作，这里我们就要用掩码将公路给抠出来显示：
//图像掩码Mask操作
threshold(proImage, proImage, 50, 255, THRESH_BINARY);//对mask进行二值化，将mask进一步处理
Mat dstImage = Mat::zeros(srcImage.size(),CV_8UC3);
srcImage.copyTo(dstImage, proImage);
imshow("掩码操作", dstImage);
显示结果：
效果并不是想象中那么好是吧哈，感觉还是语义分割要更准确些。但准确能当女朋友吗？
后期再加上些边缘检测和霍夫直线变换，和一些其他骚操作，就可以提取公路啊，车道线啥的了。
欢迎关注本知乎专栏！
如果觉的很干货，不要光收藏，点个赞呀~
点个赞祝你收藏的东西永不吃土哦~

本文首发于公众号【Opencv视觉实践】，各种有趣的小操作带你进阶opencv
之前写过一篇
【opencv】带你再学一遍直方图​mp.weixin.qq.com
里面的内容可以看下图。所以今天还要再再再学一个直方图的API：直方图反投影。
直方图反向投影是干啥的呢，它用于图像分割或在图像中查找感兴趣的对象。用什么来查找呢？用直方图。
直方图在一定程度上可以反应图像的特征，我们截取一个有固定特征的样例，比如草地，然后计算该块草地的直方图，然后用这个直方图去和整幅图像的直方图做对比，根据一定的判断条件，就能得出相似的即为草地。
看着就像是语义分割，其实一定意义上这就是语义分割，这不过直方图反向分割的依据是人为计算的（直方图）,后者分割的依据是靠在神经网络中学习得来的。
我们先看一下opencv直方图反向计算得API：
void 
一：该API得实现原理是什么呢
假设我们现在有一个四行四列得灰度图，它得灰度值如下图：
说这幅图有什么特征呢？直观上看类似于一个边角，但这是直观上，怎么表示出来呢？深度学习是靠神经网络黑箱计算出来得，我们可以用直方图。
那我们就计算这幅灰度图得直方图，如果以组距为1计算直方图并反向投影到原图，得到得为下图：
可以大概表述一下边角得特征：左下角有6个像素值相同得三角形区域，中间斜向下有四个像素值相同得边界线，以此类推。这就是用直方图得到得边角得特征。
那如果以组距为2计算直方图呢？反向投影后为：
可以看到特征描述得更为广泛了，就像深度学习里，提取更高层次得特征，虽然更为普适，但也会忽略掉一些细节特征。
我们就是拿这个反向投影所表达得特征信息，去和整幅图做对比，来得到特征相似得部分，达到分割得效果。
二：利用反向投影进行语义分割
先看一下我们今天要处理得图片：
我们今天要做得就是将这条公路给提取出来。
1，先读取原图以及样本图，并转换为HSV格式。
//【1】读取图片
为什么转HSV呢？因为HSV表达颜色更为方便区分，我们今天用到得只有前两个通道：H（色调）和S（饱和度），不用V（亮度）。
来看一下我们截取得样本图：
这是我们在公路上随便截取得一块儿样本，可以看到整条公路大概都是这个样子。
2，计算样本图得直方图并进行归一化
//【3】计算公路的直方图
代码第7行，我们计算H和S通道直方图的组距分别为20和32，我们采用更小的组距来抓取更大的特征区域。
为什么要归一化呢，直方图反向投影到原图后，原图各位置表示的是整幅图中等于该点像素值的数量，归一化后就变成概率了。
3，将计算的归一化后的直方图进行反向投影
//【5】反向投影
反向投影后的原图：
可以大概看出公路的轮廓了吧？
图像腌膜Mask的常规操作你真的信手拈来吗？​mp.weixin.qq.com
4，上一篇介绍了掩码操作，这里我们就要用掩码将公路给抠出来显示：
//图像掩码Mask操作
显示结果：
效果并不是想象中那么好是吧哈，感觉还是语义分割要更准确些。但准确能当女朋友吗？
后期再加上些边缘检测和霍夫直线变换，和一些其他骚操作，就可以提取公路啊，车道线啥的了。
欢迎关注本知乎专栏！
如果觉的很干货，不要光收藏，点个赞呀~
点个赞祝你收藏的东西永不吃土哦~

滤波反投影重建算法实现及应用（matlab）
1. 滤波反投影重建算法原理

滤波反投影重建算法常用在CT成像重建中，背后的数学原理是傅立叶变换：对投影的一维傅立叶变换等效于对原图像进行二维的傅立叶变换。（傅立叶中心切片定理）

CT重建算法大致分为解析重建算法和迭代重建算法，随着CT技术的发展，重建算法也变得多种多样，各有各的有特点。本文使用目前应用最广泛的重建算法——滤波反投影算法（FBP）作为模型的基础算法。FBP算法是在傅立叶变换理论基础之上的一种空域处理技术。它的特点是在反投影前将每一个采集投影角度下的投影进行卷积处理，从而改善点扩散函数引起的形状伪影，重建的图像质量较好。

上图应可以清晰的描述傅立叶中心切片定理的过程：对投影的一维傅立叶变换等效于对原图像进行二维的傅立叶变换
傅立叶切片定理的意义在于，通过投影上执行傅立叶变换，可以从每个投影中得到二维傅立叶变换。从而投影图像重建的问题，可以按以下方法进行求解：采集不同时间下足够多的投影（一般为180次采集），求解各个投影的一维傅立叶变换，将上述切片汇集成图像的二维傅立叶变换，再利用傅立叶反变换求得重建图像。

投影相关知识请参考fbp的matlab实现

2. 滤波反投影重建算法过程（以平行束为例）
投影重建的过程是，先把投影由线阵探测器上获得的投影数据进行一次一维傅立叶变换，再与滤波器函数进行卷积运算，得到各个方向卷积滤波后的投影数据；然后把它们沿各个方向进行反投影，即按其原路径平均分配到每一矩阵单元上，进行重叠后得到每一矩阵单元的CT值；再经过适当处理后得到被扫描物体的断层图像

算法步骤如下：

1. 将原始投影进行一次一维傅立叶变换

2. 设计合适的滤波器，在φ_i的角度下将得到原始投影p(x_r,φ_i)进行卷积滤波，得到滤波后的投影。

3. 将滤波后的投影进行反投影，得到满足x_r=r cos⁡（(θ - φ_i)）方向上的原图像的密度。

4. 将所有反投影进行叠加，得到重建后的投影。
3. 滤波器(滤波函数)和内插函数的选取
由于直接使用反投影算法会存在两个对实验结果影响很不好的因素：
不准确的数据重建图像就会产生各种伪影。
投影的数据是天然离散的，处理不当的话会产生很大的误差。
常见的滤波器有R-S滤波函数和S-L滤波函数。R-L滤波函数滤波计算简单，避免了大量的正弦、余弦计算，得到的采样序列分是分段现行的，并没有明显的降低图像质量，所以重建图像轮廓清楚，空间分辨率高。
常见的插值方法有最近邻插值和双线插值，最近邻插值即将离散点中间的缺失值用离它最近的整数处的投影值来替代。
4. FBP的matlab实现
使用R-L滤波器和最近邻插值方法
clc,clear;
%% 各个参数信息
%重建后图片像素个数
M=512;%建议先和接收传感器的个数一样

%旋转的角度 180次旋转
theta=xlsread('angles_180_3.xlsx','Sheet1','a2:a181');

%投影为512行，180列的数据，列的数据对应每一次旋转，512个接收传感器
R=xlsread('Accessory_A.xls','附件3');

%由投影进行反变换
size(R,1);%512

% 设置快速傅里叶变换的宽度
width = 2^nextpow2(size(R,1));  

%% 对投影做快速傅里叶变换并滤波
%傅立叶变换
proj_fft = fft(R, width);

% filter 滤波
% R-L是一种基础的滤波算法
filter = 2*[0:(width/2-1), width/2:-1:1]'/width;

% 滤波后结果 proj_filtered
proj_filtered = zeros(width,180);
for i = 1:180
    proj_filtered(:,i) = proj_fft(:,i).*filter;
end
figure
subplot(1,2,1),imshow(proj_fft),title('傅立叶变换')
subplot(1,2,2),imshow(proj_filtered),title('傅立叶变换+滤波')

%% 逆快速傅里叶变换并反投影
% 逆快速傅里叶变换 proj_ifft
proj_ifft = real(ifft(proj_filtered)); 
figure,imshow(proj_ifft),title('逆傅立叶变换')

%反投影到x轴，y轴
fbp = zeros(M); % 假设初值为0
for i = 1:180
    rad = theta(i);%弧度， %这个rad 是投影角，不是投影线与x轴夹角，他们之间相差 pi/2
    for x = 1:M
        for y = 1:M
            %{
            %最近邻插值法
            t = round((x-M/2)*cos(rad)-(y-M/2)*sin(rad));%将每个元素舍X入到最接近的整数。
            if t<size(R,1)/2 && t>-size(R,1)/2
                fbp(x,y)=fbp(x,y)+proj_ifft(round(t+size(R,1)/2),i);
            end
            %}
            t_temp = (x-M/2) * cos(rad) - (y-M/2) * sin(rad)+M/2  ;
             %最近邻插值法
            t = round(t_temp) ;
            if t>0 && t<=512
                fbp(x,y)=fbp(x,y)+proj_ifft(t,i);
            end
        end
    end
end
fbp = (fbp*pi)/180;%512x512 原图像每个像素位置的密度

xlswrite('problem2_origin.xlsx',fbp,'Sheet1');%将得到的重建后的图像数据写入
% 显示结果
figure,imshow(fbp),title('反投影变换后的图像')



其中每一个文件都有作用的说明，如需源文件请留言哈。（如有错误请指正）
fbp算法实现案例
github代码下载  别忘了给star哟~
参考文献：

【1】	范慧赟.CT 图像滤波反投影重建算法的研究[D].硕士学位论文,西北工业大学,2007.

【2】	余晓锷,龚剑,马建华等.CT 原理与技术[M].北京:科学出版社,2013,95-97.

【3】	毛小渊. 二维CT图像重建算法研究[D].南昌航空大学,2016.

【4】	洪虹. CT中金属伪影的校正研究[D].南方医科大学,2013.

【5】	范慧赟. CT图像滤波反投影重建算法的研究[D].西北工业大学,2007.

原博客地址：https://blog.csdn.net/shuangyue/article/details/5811384


                    
            


关于CT重建的算法有很多，在这里给大家介绍的是滤波反投影算法，其原理如下：
设f(x,y)表示需要重建的图像，用p(t,θ)表示在角度获取的f(x,y)的一个平行投影，t表示投影射线到对称中心（即旋转中心）的距离，设s为于角度θ下的投影X射线平行的坐标轴，它与t所在的坐标轴垂直，则：
        （1）
对p(t,θ)进行一维傅里叶变换
    （2）
将（1）代入（2）式得到
  （3）



图1 图像空间坐标系和投影空间坐标系
由上图可以看出，q 点在xoy坐标系内的坐标
x = Rcosα, y= Rsinα,  式中R为q点到原点的距离。
同样 q点在 tos中的坐标为：

 t = Rcos(a-θ)= x cosθ+y sinθ     (4)

 s = Rsin(a-θ) = ycosθ - xsinθ    (5)
将（4）（5）式代入（3）得到
         （6）
图像f(x,y)的二维傅里叶变换为F(u,v)
              (7)
令u = w cosθ, v = w sinθ,则
F(wcosθ,wsinθ)=P(w,θ)               (8)
由傅里叶变换知，图像函数f（x,y)可以通过其傅里叶反变换F(u,v)中恢复，即
           (9)
令 u = w cosθ, v = w sinθ，并根据(8)
       (10)
利用对称关系P(w,θ+ π) = P(-w,θ),(9)式变为
     (11)
如果令  
                        (12)
则(11)式重写为
            (13)
上述(12)(13)就是滤波反投影的主要公式，(12)中的|w|表示滤波函数。

滤波器的设计与选择
在滤波反投影重建算法中，滤波器的设计是关键。理想的滤波器是频带无限的V型滤波函数，在无穷积分区间上的积分发散，根据佩利一维纳准则 ，这一理想滤波器是不可实现的。但是如果结合具体的成像过程,则不但能够实现，而且可以达到足够的精度。为此，需要对理想滤波函数进行加窗处理，即只保留滤波函数的低频段。常用滤波函数如下几种
Ramp-Lak滤波器，它实际上是直接截断V型滤波器高频部分的结果，该滤波函数的特点是形式简单．重建的图像轮廓清楚。缺点是由于在频域中用矩形窗函数截断了滤波函数，在相应的空域中会造成振荡响应．即Gibbs现象。
把Ramp-Lak滤波器与sin(x)/x进行卷积，就得到Shepp-Logan滤波器。用Shepp-Logan滤波器重建的图像中振荡相应较小，对含噪声的数据重建出来的图像质量也较Ramp-Lak滤波函数要好。但是由于该滤波函数在高频段偏离了理想的滤波函数|w|，因而重建图像在高频段的响应不如R-L滤波函数。
其它常用滤波器还有Hamming滤波器，它是通过Ramp-Lak滤波器与Hamming窗进行卷积而得到。Hanning滤波器，它实际上是Ramp-Lak滤波器与Hanning窗进行卷积。

关于其他CT重建技术，请参考本人的上一篇文章
http://blog.csdn.net/shuangyue/archive/2010/08/13/5810789.aspx
            

            </div>