本文主要介绍“图像的操作插值算法”.
本文为了简明体现整个操作过程,因此程序框架及软件设计并不规范,所以代码仅供参考。
上一篇介绍了图像坐标系的转换,见:连载:徒手写一个DICOM阅图软件(004)之图像的操作-CSDN博客
0. 临近插值的问题
上篇中,我们完成了对图像的操作(移动、缩放、旋转)。尤其在图像的缩放和旋转中,图像会表现出马赛克状,如下图:
原因是在图像显示时,咱们采用的是临近插值算法,因此,在图像显示时,采用双线性插值即可改善此问题。
1.双线性插值的代码
代码如下:
void ShowWork::DrawInPixel( )
{
short* pPixelBuffer = mData.GetPixelBuffer( ); // DICOM原始像素
if ( pPixelBuffer == NULL )
return;
int pixelWidth = mData.GetWidth( ); // DICOM图像的宽度
int pixelHeight = mData.GetHeight( ); // DICOM图像的高度
float a, b;
mData.GetWLFactor( a, b ); // 窗宽窗位转换系数
int w = mShowSize[0]; // 显示尺寸
int h = mShowSize[1];
M3DMatrix44f ScreenToDicom;
m3dInvertMatrix44( ScreenToDicom, mDicomToScreen );
M3DVector3f v1, v2;
UCHAR c = 0;
for ( int y = 0; y < h; y++ )
{
for ( int x = 0; x < w; x++ )
{
m3dLoadVector3( v1, float( x ), float( y ), 0 ); // 当前窗口位置
m3dTransformVector3( v2, v1, ScreenToDicom ); // 计算出原始图像中位置
// 测试是否越界
if ( v2[0] < 1 || v2[0] >= pixelWidth-1 ||
v2[1] < 1 || v2[1] >= pixelHeight-1 )
{
mShowPixel[x + y * w ] = 0; // 如果越界,则置0
continue;
}
// 未越界,则取出CT值
/* // 旧算法:临近插值
int posX = int( v2[0] + 0.5f );
int posY = int( v2[1] + 0.5f );
short pixel = pPixelBuffer[ posX + posY * pixelWidth];
float v = pixel * a + b; // 将原始数据进行转换,y=x*a+b
*/
// 以下是双线性插值的算法
int intX = int( v2[0] );
int intY = int( v2[1] );
float fx = v2[0] - intX;
float fy = v2[1] - intY;
short* tempP = pPixelBuffer + intX + intY*pixelWidth;
// 取出周围4个点
short v00 = *tempP;
short v01 = *(tempP+1);
short v10 = *(tempP+pixelWidth);
short v11 = *(tempP+pixelWidth+1);
float fxInv = 1.0f - fx;
float fyInv = 1.0f - fy;
float v0001 = v00*fxInv + v01*fx;
float v1011 = v10*fxInv + v11*fx;
float v00011011 = v0001 * fyInv + v1011 * fy;
float v = v00011011 * a + b; // 计算显示的像素值
if ( v < 0.0f )
c = 0;
else if ( v > 255.0f )
c = 255;
else
c = UCHAR( v );
mShowPixel[x + y * w] = RGB( c, c, c ); // 显示像素为32bit,因为显示灰度,所以3通道颜色一样。
}
}
}
2. 双线性插值的效果
效果如下:
这样就没什么马赛克了。
当然,当双线性插值的图像显示大约5倍以上,还是能肉眼看到马赛克,此时可以采用立方插值来改善(医学图像下双线性插值基本上够了)