反Radon变换 C++实现

于 2021-04-29 16:38:16 发布
阅读量922 收藏 13
点赞数 2
分类专栏： 学基于MFC的图像处理时随便记记 文章标签： c++ 图像处理
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.youkuaiyun.com/weixin_44290152/article/details/116272054
版权
本文详细介绍了Radon变换的原理，包括其在图像滤波和反投影过程中的应用。首先，通过添加RampFilter对图像进行滤波处理，然后通过反投影算法将滤波后的频域数据转换回图像。在反投影过程中，采用线性插值计算每个像素的值，并对结果进行了归一化处理。最后，创建新的DIB对象并写入重建后的图像数据。代码示例展示了整个流程，包括FFT和IFFT的使用，以及图像数据的处理和存储。

摘要生成于 C知道 ，由 DeepSeek-R1 满血版支持， 前往体验 >

文章目录

Radon变换原理

百度百科里面说的蛮清楚的，可以自己看一下。

频域添加了Ramp Filter(如下)
H ( ω ) = ∣ ω ∣ b ω ( ω ) b ω ( ω ) = { 1 ∣ ω ∣ < ω 0 o t h e r s H(\omega) = |\omega|b_\omega (\omega)\\ b_\omega(\omega) = \begin{cases} 1 \quad |\omega| < \omega\\ 0 \quad others \end{cases} H(ω)=∣ω∣bω​(ω)bω​(ω)={1∣ω∣<ω0others​
正弦图如下
原图
重建后图像如下
在这里插入图片描述
里面有一点波纹可能是频谱泄露或者滤波器的问题

代码实现

对图像做滤波

int nLevel = ceil(log(1.0*m_nHeight) / log(2.0));
	int FFTLen = pow(2.0, nLevel);
	complex<double>* wFData = new complex<double>[FFTLen];
	complex<double>* wTData = new complex<double>[FFTLen];
	
	double* rawDataBuffer = new double[m_nWidth*m_nHeight];
	for (int i = 0; i < m_nWidth*m_nHeight; i++)
		rawDataBuffer[i] = 0;

	//将数据添加滤波器
	for (int lineX = 0; lineX < m_nWidth; lineX++)
	{	
		unsigned char tempPixel;
		for (int lineY = 0; lineY < FFTLen; lineY++)
		{
			tempPixel = lineY < m_nHeight ? m_pDibBits[lineX + lineY * m_nWidth] : 0;
			wTData[lineY] = complex<double>(tempPixel, 0);
			wFData[lineY] = complex<double>(0, 0);
		}
		//转到频域
		FFT_1D(wTData, wFData, nLevel);
		
		double filterFactor;
		// 添加滤波器
		for (int w = 0; w < FFTLen/2; w++)
		{
			filterFactor = (w *1.0 / FFTLen * 2);
			wFData[w] = complex<double>(wFData[w].real()*filterFactor, wFData[w].imag()*filterFactor);
			wFData[FFTLen - w - 1] = complex<double>\
				(wFData[FFTLen - w - 1].real()*filterFactor, wFData[FFTLen - w - 1].imag()*filterFactor);
		}

		IFFT_1D(wFData, wTData, nLevel);

		//加过滤波器后的数据放回buffer的原位置
		for (int lineY = 0; lineY < m_nHeight; lineY++)
		{
			rawDataBuffer[lineY*m_nWidth + lineX] = wTData[lineY].real();
		}
	}

对图像做反投影

	//反投影
	int t1 ,t2 , newHeight = ceil(sqrt(m_nHeight * m_nHeight / 2.0));
	double t;
	//maxData 反投影后像素最大值 minData 反投影后像素最小值
	//tempPixel 算每个Pixel的中间缓存 pixelWeight 做插值时的权重
	double maxData,minData,tempPixel,pixelWeight; 
	maxData = minData = rawDataBuffer[0];
	tempPixel = pixelWeight = 0;

	double* dataBuffer = new double[newHeight*newHeight];
	
	for (int lineX = 0; lineX < newHeight; lineX++)
	{
		for (int lineY = 0; lineY < newHeight; lineY++)
		{
			tempPixel = 0;
			for (int theta = 0; theta < m_nWidth; theta++)
			{
				t = (lineX - newHeight / 2) * cos(theta*Pi / 180) \
					+ (lineY - newHeight / 2) * sin(theta*Pi / 180)\
					+ m_nHeight / 2;

				t1 = floor(t);
				t1 = t1 % m_nHeight;//防止越界
				t2 = (t1 == m_nHeight - 1)? t1 : ceil(t); //超限的话就一起是t了

				pixelWeight = t - t1;

				tempPixel += pixelWeight * rawDataBuffer[t1*m_nWidth + theta] / 180\
					+ (1 - pixelWeight) * rawDataBuffer[t2*m_nWidth + theta] / 180;
			}
			dataBuffer[lineY*newHeight + lineX] = tempPixel;
			maxData = maxData > tempPixel ? maxData : tempPixel;
			minData = minData < tempPixel ? minData : tempPixel;
		}
	}

创建新的dib对象并写入数据

//销毁原数据
	Destroy();
	//创建新的数据区
	int BPP = 8;
	Create(newHeight, newHeight, BPP, 0);
	m_nWidth = newHeight;
	m_nHeight = newHeight;
	//调色板赋值
	if (IsIndexed())
	{
		int nColors = 256;
		if (nColors > 0)
		{
			RGBQUAD *pal = new RGBQUAD[nColors];
			for (int i = 0; i < nColors; i++)
			{
				pal[i].rgbRed = i;
				pal[i].rgbBlue = i;
				pal[i].rgbGreen = i;
				pal[i].rgbReserved = BYTE(0);
			}
			SetColorTable(0, nColors, pal);
			delete[] pal;
		}
	}
	//数据传递
	m_pDibBits = (unsigned char*)GetBits() + (m_nHeight - 1)*GetPitch();
	for (int i = 0; i < newHeight*newHeight; i++)
	{
		//可能出现了越界 排一下bug
		m_pDibBits[i] =(BYTE) ((dataBuffer[i] - minData) / (maxData - minData) * 255);
	}

	delete[] rawDataBuffer;
	delete[] dataBuffer;
	delete[] wFData;
	delete[] wTData;

完整代码

void CDib::reConstruction()
{
	if (!this) return;
	// 准备部分
	int nLevel = ceil(log(1.0*m_nHeight) / log(2.0));
	int FFTLen = pow(2.0, nLevel);
	complex<double>* wFData = new complex<double>[FFTLen];
	complex<double>* wTData = new complex<double>[FFTLen];
	
	double* rawDataBuffer = new double[m_nWidth*m_nHeight];
	for (int i = 0; i < m_nWidth*m_nHeight; i++)
		rawDataBuffer[i] = 0;

	//将数据添加滤波器
	for (int lineX = 0; lineX < m_nWidth; lineX++)
	{	
		unsigned char tempPixel;
		for (int lineY = 0; lineY < FFTLen; lineY++)
		{
			tempPixel = lineY < m_nHeight ? m_pDibBits[lineX + lineY * m_nWidth] : 0;
			wTData[lineY] = complex<double>(tempPixel, 0);
			wFData[lineY] = complex<double>(0, 0);
		}
		//转到频域
		FFT_1D(wTData, wFData, nLevel);
		
		double filterFactor;
		// 添加滤波器
		for (int w = 0; w < FFTLen/2; w++)
		{
			filterFactor = (w *1.0 / FFTLen * 2);
			wFData[w] = complex<double>(wFData[w].real()*filterFactor, wFData[w].imag()*filterFactor);
			wFData[FFTLen - w - 1] = complex<double>\
				(wFData[FFTLen - w - 1].real()*filterFactor, wFData[FFTLen - w - 1].imag()*filterFactor);
		}

		IFFT_1D(wFData, wTData, nLevel);

		//加过滤波器后的数据放回buffer的原位置
		for (int lineY = 0; lineY < m_nHeight; lineY++)
		{
			rawDataBuffer[lineY*m_nWidth + lineX] = wTData[lineY].real();
		}
	}

	//反投影
	int t1 ,t2 , newHeight = ceil(sqrt(m_nHeight * m_nHeight / 2.0));
	double t;
	//maxData 反投影后像素最大值 minData 反投影后像素最小值
	//tempPixel 算每个Pixel的中间缓存 pixelWeight 做插值时的权重
	double maxData,minData,tempPixel,pixelWeight; 
	maxData = minData = rawDataBuffer[0];
	tempPixel = pixelWeight = 0;

	double* dataBuffer = new double[newHeight*newHeight];
	
	for (int lineX = 0; lineX < newHeight; lineX++)
	{
		for (int lineY = 0; lineY < newHeight; lineY++)
		{
			tempPixel = 0;
			for (int theta = 0; theta < m_nWidth; theta++)
			{
				t = (lineX - newHeight / 2) * cos(theta*Pi / 180) \
					+ (lineY - newHeight / 2) * sin(theta*Pi / 180)\
					+ m_nHeight / 2;

				t1 = floor(t);
				t1 = t1 % m_nHeight;//防止越界
				t2 = (t1 == m_nHeight - 1)? t1 : ceil(t); //超限的话就一起是t了

				pixelWeight = t - t1;

				tempPixel += pixelWeight * rawDataBuffer[t1*m_nWidth + theta] / 180\
					+ (1 - pixelWeight) * rawDataBuffer[t2*m_nWidth + theta] / 180;
			}
			dataBuffer[lineY*newHeight + lineX] = tempPixel;
			maxData = maxData > tempPixel ? maxData : tempPixel;
			minData = minData < tempPixel ? minData : tempPixel;
		}
	}

	//销毁原数据
	Destroy();
	//创建新的数据区
	int BPP = 8;
	Create(newHeight, newHeight, BPP, 0);
	m_nWidth = newHeight;
	m_nHeight = newHeight;
	//调色板赋值
	if (IsIndexed())
	{
		int nColors = 256;
		if (nColors > 0)
		{
			RGBQUAD *pal = new RGBQUAD[nColors];
			for (int i = 0; i < nColors; i++)
			{
				pal[i].rgbRed = i;
				pal[i].rgbBlue = i;
				pal[i].rgbGreen = i;
				pal[i].rgbReserved = BYTE(0);
			}
			SetColorTable(0, nColors, pal);
			delete[] pal;
		}
	}
	//数据传递
	m_pDibBits = (unsigned char*)GetBits() + (m_nHeight - 1)*GetPitch();
	for (int i = 0; i < newHeight*newHeight; i++)
	{
		//可能出现了越界 排一下bug
		m_pDibBits[i] =(BYTE) ((dataBuffer[i] - minData) / (maxData - minData) * 255);
	}

	delete[] rawDataBuffer;
	delete[] dataBuffer;
	delete[] wFData;
	delete[] wTData;
}
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值