基于主轴变换的医学图像倾斜校正

问题描述

在MRI扫描的时候，可能会由于病人体位的关系，出现扫描图像是倾斜的状况，如下图所示。

在特定的情况下，病人的倾斜状态会影响到后面我们数据的处理。比如在MRI头部自动定位通常是通过AC-PC连线来确定的，如果病人头部存在倾斜，则AC-PC点可能并不在同一个矢状位面上，给定位增加了难度。

如何进行倾斜校正

切斜角度计算的方法很多，常见的有基于主轴变换的方法[1]和基于傅里叶变换的方法[2]。我这里主要介绍使用主轴变换来计算偏转角度。

图像矩

图像矩主要表征了图像区域的几何特征，又称为几何矩。由于其具有不变特征，又称其为不变矩，它的主要思想是使用对变换不敏感的基于区域的几个矩作为形状特征。在医学图像处理经常被应用来作为预配准步骤。

在图像处理中，我们把灰度图像看作为一块平板，每个像素点的值为该点处的密度。对某点求期望就是该图像在该点处的矩。我们一般讲的图像矩，就是原点矩。一阶矩和零阶矩就可以计算某个形状的中心，而二阶矩就可以计算形状的方向。

零阶矩：

M00=∑x∑yV(x,y)M_{00}=\sum_{x} \sum_{y} V(x,y)M00=∑x∑yV(x,y)

一阶矩：

M10=∑x∑yx⋅V(x,y)M_{10}=\sum_{x} \sum_{y}x \cdot V(x,y)M10=∑x∑yx⋅V(x,y)
M01=∑x∑yy⋅V(x,y)M_{01}=\sum_{x} \sum_{y}y \cdot V(x,y)M01=∑x∑yy⋅V(x,y)

重心坐标：

xc=M10M00x_{c}=\frac{M_{10}}{M_{00}}xc=M00M10
yc=M01M00y_{c}=\frac{M_{01}}{M_{00}}yc=M00M01

二阶矩：

M20=∑x∑yx2⋅V(x,y)M_{20}=\sum_{x} \sum_{y} x^{2} \cdot V(x, y)M20=∑x∑yx2⋅V(x,y)
M02=∑x∑yy2⋅V(x,y)M_{02}=\sum_{x} \sum_{y} y^{2} \cdot V(x, y)M02=∑x∑yy2⋅V(x,y)
M11=∑x∑yxy⋅V(x,y)M_{11}=\sum_{x} \sum_{y} xy \cdot V(x, y)M11=∑x∑yxy⋅V(x,y)

物体形状方向计算：

θ=12arctan⁡(2ba−c),θ∈[−90∘,90∘]\theta=\frac{1}{2} \arctan \left(\frac{2 b}{a-c}\right), \theta \in\left[-90^{\circ}, 90^{\circ}\right]θ=21arctan(a−c2b),θ∈[−90∘,90∘]

其中, a=M20M00−xc2,b=M11M00−xcyc,c=M02M00−yc2\mathrm{a}=\frac{M_{20}}{M_{00}}-x_{c}^{2}, b=\frac{M_{11}}{M_{00}}-x_{c} \mathrm{y}_{\mathrm{c}}, c=\frac{M_{02}}{M_{00}}-y_{c}^{2}a=M00M20−xc2,b=M00M11−xcyc,c=M00M02−yc2

ITK代码实现

创建测试图像

LabelType::Pointer CreateSampleImage()
{LabelType::Pointer image = LabelType::New();LabelType::IndexType start;start[0] = 0;start[1] = 0;LabelType::SizeType size;size[0] = 200;size[1] = 200;LabelType::RegionType region;region.SetSize(size);region.SetIndex(start);image->SetRegions(region);image->Allocate();// 椭圆方程// (x-xc)^2/a^2 + (y-yc)^2/b^2 = 1itk::ImageRegionIterator< LabelType >  it(image, image->GetRequestedRegion());it.Begin();while (!it.IsAtEnd()){if (pow((int)it.GetIndex()[0] - 90, 2) / (57 * 57) + pow((int)it.GetIndex()[1] - 110, 2) / (78 * 78)  < 1)it.Set(1);elseit.Set(0);++it;}// 加一点旋转using MatrixType = itk::Matrix<double, 2, 2>;MatrixType matrix;matrix[0][0] = cos(3.14 / 20);matrix[0][1] = sin(3.14 / 20);matrix[1][0] = -sin(3.14 / 20);matrix[1][1] = cos(3.14 / 20);// 输出下这个旋转矩阵cout << "Init Transform Matrix : " << matrix <<endl;using PointType = itk::Point<double, 2>;PointType center;center[0] = 100;center[1] = 100;typedef itk::CenteredRigid2DTransform< double > TransformType;TransformType::Pointer transform = TransformType::New();transform->SetMatrix(matrix);transform->SetCenter(center);// 应用上面的变换重采样图像typedef itk::ResampleImageFilter<LabelType, LabelType > ResampleFilterType;ResampleFilterType::Pointer resample = ResampleFilterType::New();resample->SetTransform(transform);resample->SetInput(image);resample->SetSize(image->GetLargestPossibleRegion().GetSize());resample->SetOutputOrigin(image->GetOrigin());resample->SetOutputSpacing(image->GetSpacing());resample->SetOutputDirection(image->GetDirection());resample->SetDefaultPixelValue(0);resample->Update();// 写出我们生成的倾斜图像using WriterType = itk::ImageFileWriter<LabelType>;WriterType::Pointer writer = WriterType::New();writer->SetFileName("sample.mha");writer->SetInput(resample->GetOutput());writer->SetImageIO(itk::MetaImageIO::New());writer->Write();return image;
}

然后我们就得到这样一个二值图像

这个变换的旋转矩阵为：

Init Transform Matrix : 0.987701 0.156356 -0.156356 0.987701

计算主轴变换矩阵并重采样

#include <itkImageFileReader.h>
#include <itkImageFileWriter.h>
#include <itkGDCMImageIO.h>
#include <itkMetaImageIO.h>
#include <itkImageMomentsCalculator.h>
#include <itkOtsuThresholdImageFilter.h>
#include <itkBinaryFillholeImageFilter.h>
#include <itkCenteredRigid2DTransform.h>
#include <itkResampleImageFilter.h>int main(int argc, char * argv[])
{using LabelType = itk::Image<unsigned char, 2>;// 读取图像using ReaderType = itk::ImageFileReader<LabelType>;ReaderType::Pointer reader = ReaderType::New();reader->SetFileName("sample.mha");reader->SetImageIO(itk::MetaImageIO::New());reader->Update();// 计算图像矩using ImageMomentsFilter = itk::ImageMomentsCalculator<LabelType>;ImageMomentsFilter::Pointer moments = ImageMomentsFilter::New();moments->SetImage(reader->GetOutput());moments->Compute();moments->Print(cout);// 这里我们要将偏转的图像重新摆正，// 所以旋转矩阵使用主轴变换的逆矩阵using MatrixType = itk::Matrix<double, 2, 2>;MatrixType trans(moments->GetPrincipalAxes().GetInverse());typedef itk::CenteredRigid2DTransform< double > TransformType;TransformType::Pointer transform = TransformType::New();transform->SetMatrix(trans);typedef itk::ResampleImageFilter<LabelType,LabelType > ResampleFilterType;ResampleFilterType::Pointer resample = ResampleFilterType::New();resample->SetTransform(transform);resample->SetInput(reader->GetOutput());resample->SetSize(reader->GetOutput()->GetLargestPossibleRegion().GetSize());resample->SetOutputOrigin(reader->GetOutput()->GetOrigin());resample->SetOutputSpacing(reader->GetOutput()->GetSpacing());resample->SetOutputDirection(reader->GetOutput()->GetDirection());resample->SetDefaultPixelValue(0);try{resample->Update();}catch (itk::ExceptionObject & excp){std::cerr << excp << std::endl;}using WriterType = itk::ImageFileWriter<LabelType>;WriterType::Pointer writer = WriterType::New();writer->SetFileName("resample.mha");writer->SetInput(resample->GetOutput());writer->SetImageIO(itk::MetaImageIO::New());try{writer->Write();}catch (itk::ExceptionObject & excp){std::cerr << excp << std::endl;}//moments->Print(std::cout, 0);return EXIT_SUCCESS;
}

我们的校正的倾斜图像

主轴变换的相关计算结果如下：

ImageMomentsCalculator (000001C489893E70)RTTI typeinfo:   class itk::ImageMomentsCalculator<class itk::Image<unsigned char,2> >Reference Count: 1Modified Time: 47Debug: OffObject Name:Observers:noneImage: 000001C4898B1EB0Valid: 1Zeroth Moment about origin: 13680First Moment about origin: [88.5567, 108.308]Second Moment about origin: 809.857 -108.129
-108.129 1477.97Center of Gravity: [88.5567, 108.308]Second central moments: 809.857 -108.129
-108.129 1477.97Principal Moments: [1.08454e+07, 2.0452e+07]Principal axes: 0.987776 0.155883
-0.155883 0.987776

参考文献

[1]. 潘梅森. 医学图像倾斜校正方法与应用研究[D].中南大学,2011.
[2]. https://blog.csdn.net/linqianbi/article/details/78863839