Matlab 进行两个体数据的配准 并显示配准误差 彩色化

进行两个体数据间的配准，并且显示配准后的误差：
http://cn.mathworks.com/help/images/ref/imregister.html?requestedDomain=cn.mathworks.com
这里采用的图片是matlab子带的两张MR膝盖图，“knee1.dcm” 作为参考图像，"knee2.dcm"为浮动图像！

Plain Text code

?
可能接下来大家关注的问题就是这两幅图像到底有什么区别，这种区别又有多大的可视化程度，下面就为推荐两个比较好用的函数用于观测两幅图像的区别。

Plain Text code

?
imshowpair函数就是指以成双成对的形式显示图片，其中一个重要的参数就是‘method’，他又4个选择

（1）‘falsecolor’ 字面意思理解就是伪彩色的意思了，其实就是把两幅图像的差异用色彩来表示，这个是默认的参数。

（2）‘blend’ 这是一种混合透明处理类型，技术文档的翻译是alpha blending，大家自己理解吧。

（3）‘diff’ 这是用灰度信息来表示亮度图像之间的差异，这是对应‘falsecolor’的一种方式。

（4）参数‘monotage’可以理解成‘蒙太奇’，这是一种视频剪辑的艺术手法，其实在这里我们理解成拼接的方法就可以了。

为什么在这里罗里吧嗦的说这么多的显示呢，大家知道"人靠衣装，美靠...."(就不多说了吧)，总之就是一个好的视觉效果能给人以耳目一新的效果。



嗯嗯，这个就是蒙太奇的效果了，





这两个则分别是伪彩色，混合透明处理了，至于大家接受那个就要看自己的爱好了

说到了这里终于结束了这关没有意义的读图环节，请大家原谅我的无耻吧。

二，初始配准（粗配准）

初始配准就是大致的使图像对其，使其差别不要太明显，以方便下一步的精细配准环节。

函数imregconfig这在个环节可是主角，从名字上看就知道他要设置一些参数和方法了，其实他真正的作用是配置优化器和度量准则，

Plain Text code

?
参数modality指定fixed image, moving image之间的关系，有两种选择‘monomodal’, 'multimodal'两种，分别质量两幅图像是单一模态还是多模态，根据需要自己选择。

返回的参数optimizer是用于优化度量准则的优化算法，这里有

registration.optimizer.RegularStepGradientDescent 或者 registration.optimizer.OnePlusOneEvolutionary两种可供选择。

输出参数metric则是注明了度量两幅图片相似度的方法，提供了均方误差（registration.metric.MeanSquares）和互信息（registration.metric.MattesMutualInformation）两种供选择。

当然大家也可以根据结构扩充这两个参数。

到这里优化器和度量准别就已将做好了，是不是简单到没朋友。

要上大菜了，配准代码

Plain Text code

?
imregister函数根据取得的optimizer,metric参数对2D,3D参考图像做变换（transform）目的是fixed，moving image对其，大家关注到有一个参数‘affine’，他是指该变化是仿射变换，同样该参数还可以选为

‘translation’ （x,y）坐标平移变换，不牵涉到旋转个尺度变换

‘rigid’ 刚性变换（平移和旋转）

‘similarity’ 改变换包括了平移，旋转和尺度变换

‘affine’ 在similarity的基础上加入了shear（图像的剪辑）



该图片就是粗配准的结果了，大家可以在右上角看到明显的不重合现象。

三，提高配准精度

粗配准的结果一般情况下达不到实际应用的要求，为此很有必要进一步提高精度，如果有对精度要求不高的朋友看到这里就可以结束了。

Plain Text code

?
这两条指令可以看到默认生成的优化器和度量函数参数，当然这里提高精度的途径就是通过修改这两个参数了！



在这里我们通过修改两个参数，观察对配准效果的改进：

（1）改变优化器的步长已达到对更加精细的变换。

Plain Text code

?
把原步长缩小为原来的3.5倍，结果如下



貌似效果还是有点的啊，大家在看右上角的阴影好像不见了

（2）在（1）基础上改变最大迭代次数

Plain Text code

?
效果如下：正上的阴影好像也减小了



四，改变初始条件提高精度

这里的思想就像我们在做雕塑一样，假如我们要用石头雕一个人，首先我们可以大刀阔斧的把头部留出来，然后把脖子留的比头部更细，简单的说就是美女留出S轮廓，或者o型的（哈哈，对号入座就可以了），下一步精雕细琢的时候就会轻松很多，这里的初始条件就是先用简单的变换做出一个初始配准图像，然后以初始配准的结果作为输入做精细配准。

大致做法如下：

Plain Text code

?
用similarity的变换方式做初始配准，或者你还可以用rigid，transform的方式都可以

Plain Text code

?
在这里imregtform把变化矩阵输出；

然后用imref2d限制变换后的图像与参考图像有相同的坐标分布

Plain Text code

?
imwarp函数执行几何变换，当然依据则是tformSimilarity的变换矩阵了。

Plain Text code

?
得到的tformsimilarity.T就是传说中的变换矩阵了

tformSimilarity.T= 1.0331 -0.1110 0

0.1110 1.0331 0

-51.1491 6.9891 1.0000

下面就是精配准的部分了：

Plain Text code

?
初始配准结果：


进一步精细配准：


五，到这里就是你说了算了Deciding When Enough is Enough

Plain Text code

?
选择一个合适的，理想的你想要的结果，去飞去装逼吧。

代码全文如下：

%% Registering Multimodal MRI Images
% This example shows how you can use |imregister| to automatically
% align two magnetic resonance images (MRI) to a common coordinate
% system using intensity-based image registration.  Unlike some other
% techniques, it does not find features or use control points.
% Intensity-based registration is often well-suited for medical and
% remotely sensed imagery.

% Copyright 2011-2013 The MathWorks, Inc.

%% Step 1: Load Images
% This example uses two magnetic resonance (MRI) images of a knee.
% The fixed image is a spin echo image, while the moving image is a
% spin echo image with inversion recovery.  The two sagittal slices
% were acquired at the same time but are slightly out of alignment.

fixed = dicomread('knee1.dcm');
moving = dicomread('knee2.dcm');

%%
% The |imshowpair| function is a useful function for visualizing
% images during every part of the registration process.  Use it to see
% the two images individually in a montage fashion or display them
% stacked to show the amount of misregistration.

figure, imshowpair(moving, fixed, 'montage')
title('Unregistered')

%%
% In the overlapping image from |imshowpair|, gray areas correspond to
% areas that have similar intensities, while magenta and green areas
% show places where one image is brighter than the other.  In some
% image pairs, green and magenta areas don't always indicate
% misregistration, but in this example it's easy to use the color
% information to see where they do.

figure, imshowpair(moving, fixed)
title('Unregistered')

%% Step 2: Set up the Initial Registration
% The |imregconfig| function makes it easy to pick the correct
% optimizer and metric configuration to use with |imregister|. These
% two images have different intensity distributions, which suggests a
% multimodal configuration.

[optimizer,metric] = imregconfig('multimodal');

%%
% The distortion between the two images includes scaling, rotation,
% and (possibly) shear.  Use an affine transformation to register the
% images.
%
% It's very, very rare that |imregister| will align images perfectly
% with the default settings.  Nevertheless, using them is a useful way
% to decide which properties to tune first.

movingRegisteredDefault = imregister(moving, fixed, 'affine', optimizer, metric);

figure, imshowpair(movingRegisteredDefault, fixed)
title('A: Default registration')

%% Step 3: Improve the Registration
% The initial registration is not very good. There are still significant
% regions of poor alignment, particularly along the right edge.  Try to
% improve the registration by adjusting the optimizer and metric
% configuration properties.
%
% The optimizer and metric variables are objects whose properties
% control the registration.

disp(optimizer)
disp(metric)

%%
% The InitialRadius property of the optimizer controls the initial step
% size used in parameter space to refine the geometric transformation. When
% multi-modal registration problems do not converge with the default
% parameters, the InitialRadius is a good first parameter to adjust. Start
% by reducing the default value of InitialRadius by a scale factor of 3.

optimizer.InitialRadius = optimizer.InitialRadius/3.5;

movingRegisteredAdjustedInitialRadius = imregister(moving, fixed, 'affine', optimizer, metric);
figure, imshowpair(movingRegisteredAdjustedInitialRadius, fixed)
title('Adjusted InitialRadius')

%%
% Adjusting the InitialRadius had a positive impact. There is a noticeable
% improvement in the alignment of the images at the top and right edges.

%%
% The MaximumIterations property of the optimizer controls the maximum
% number of iterations that the optimizer will be allowed to take.
% Increasing the MaximumIterations allows the registration search to run
% longer and potentially find better registration results. Does the
% registration continue to improve if the InitialRadius from the last step
% is used with a large number of interations?

optimizer.MaximumIterations = 300;
movingRegisteredAdjustedInitialRadius300 = imregister(moving, fixed, 'affine', optimizer, metric);

figure, imshowpair(movingRegisteredAdjustedInitialRadius300, fixed)
title('B: Adjusted InitialRadius, MaximumIterations = 300, Adjusted InitialRadius.')

%%
% Further improvement in registration were achieved by reusing the
% InitialRadius optimizer setting from the previous registration and
% allowing the optimizer to take a large number of iterations.

%% Step 4: Use Initial Conditions to Improve Registration
% Optimization based registration works best when a good initial condition
% can be given for the registration that relates the moving and fixed
% images. A useful technique for getting improved registration results is
% to start with more simple transformation types like 'rigid', and then use
% the resulting transformation as an initial condition for more complicated
% transformation types like 'affine'.
%
% The function |imregtform| uses the same algorithm as imregister, but
% returns a geometric transformation object as output instead of a
% registered output image. Use |imregtform| to get an initial
% transformation estimate based on a 'similarity' model
% (translation,rotation, and scale).
%
% The previous registration results showed in improvement after modifying
% the MaximumIterations and InitialRadius properties of the optimizer.
% Keep these optimizer settings while using initial conditions while
% attempting to refine the registration further.

tformSimilarity = imregtform(moving,fixed,'similarity',optimizer,metric);

%%
% Because the registration is being solved in the default MATLAB coordinate
% system, also known as the intrinsic coordinate system, obtain the default
% spatial referencing object that defines the location and resolution of
% the fixed image.

Rfixed = imref2d(size(fixed));

%%
% Use |imwarp| to apply the geometric transformation output from
% |imregtform| to the moving image to align it with the fixed image. Use
% the 'OutputView' option in |imwarp| to specify the world limits and
% resolution of the output resampled image. Specifying Rfixed as the
% 'OutputView' forces the resampled moving image to have the same
% resolution and world limits as the fixed image.

movingRegisteredRigid = imwarp(moving,tformSimilarity,'OutputView',Rfixed);
figure, imshowpair(movingRegisteredRigid, fixed);
title('C: Registration based on similarity transformation model.');

%%
% The "T" property of the output geometric transformation defines the
% transformation matrix that maps points in moving to corresponding
% points in fixed.

tformSimilarity.T

%%
% Use the 'InitialTransformation' Name/Value in imregister to refine this
% registration by using an 'affine' transformation model with the 'similarity'
% results used as an initial condition for the geometric transformation.
% This refined estimate for the registration includes the possibility of
% shear.

movingRegisteredAffineWithIC = imregister(moving,fixed,'affine',optimizer,metric,...
'InitialTransformation',tformSimilarity);
figure, imshowpair(movingRegisteredAffineWithIC,fixed);
title('D: Registration from affine model based on similarity initial condition.');

%%
% Using the 'InitialTransformation' to refine the 'similarity' result of
% |imregtform| with a full affine model has also yielded a nice
% registration result.

%% Step 5: Deciding When Enough is Enough
% Comparing the results of running |imregister| with different
% configurations and initial conditions, it becomes apparent that there are
% a large number of input parameters that can be varied in imregister, each
% of which may lead to different registration results.

figure
imshowpair(movingRegisteredDefault, fixed)
title('A - Default settings.');

figure
imshowpair(movingRegisteredAdjustedInitialRadius, fixed)
title('B - Adjusted InitialRadius, 100 Iterations.');

figure
imshowpair(movingRegisteredAdjustedInitialRadius300, fixed)
title('C - Adjusted InitialRadius, 300 Iterations.');

figure
imshowpair(movingRegisteredAffineWithIC, fixed)
title('D - Registration from affine model based on similarity initial condition.');

%%
% It can be difficult to quantitatively compare registration results
% because there is no one quality metric that accurately describes the
% alignment of two images. Often, registration results must be judged
% qualitatively by visualizing the results. In The results above, the
% registration results in C) and D) are both very good and are difficult to
% tell apart visually.

%% Step 6: Alternate Visualizations
% Often as the quality of multimodal registrations improve it becomes more
% difficult to judge the quality of registration visually.  This is because
% the intensity differences can obscure areas of misalignment.  Sometimes
% switching to a different display mode for |imshowpair| exposes hidden
% details.  (This is not always the case.)

displayEndOfDemoMessage(mfilename)