-
Notifications
You must be signed in to change notification settings - Fork 1
/
register_images.m
358 lines (314 loc) · 14.8 KB
/
register_images.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
function [Ireg,Bx,By,Fx,Fy,transformedMatrix] = register_images(Imoving,Istatic,Options)
% This function register_images is the most easy way to register two
% images both affine and nonrigidly.
%
% Features:
% - It can be used with images from different type of scans or modalities.
% - It uses both a rigid transform and a nonrigid registration.
% - It uses multilevel refinement
% - It can be used with images of different sizes.
% - The function will automaticaly detect single modality or multiple
% modalities, and choose the right registration method.
%
% [Ireg,Bx,By,Fx,Fy] = register_images(Imoving,Istatic,Options);
%
% Inputs,
% Imoving : The image which will be registerd
% Istatic : The image on which Imoving will be registered
% Options : Registration options, see help below
%
% Outputs,
% Ireg : The registered moving image
% Bx, By : The backwards transformation fields of the pixels in
% x and y direction seen from the static image to the moving image.
% Fx, Fy : The (approximated) forward transformation fields of the pixels in
% x and y direction seen from the moving image to the static image.
% (See the function backwards2forwards)
% Options,
% Options.SigmaFluid : The sigma of the gaussian smoothing kernel of the pixel
% velocity field / update field, this is a form of fluid
% regularization, (default 4)
% Options.SigmaDiff : The sigma for smoothing the transformation field
% is not part of the orignal demon registration, this is
% a form of diffusion regularization, (default 1)
% Options.Interpolation : Linear (default) or Cubic.
% Options.Alpha : Constant which reduces the influence of edges (and noise)
% and limits the update speed (default 4).
% Options.Similarity : Choose 'p' for single modality and 'm' for
% images of different modalities. (default autodetect)
% Options.Registration: Rigid, Affine, NonRigid , Similarity
% Options.MaxRef : Maximum number of grid refinements steps.
% Options.Verbose: Display Debug information 0,1 or 2
%
% Notes,
% In case of Multiple Modalities affine registration is done with mutual
% information. The non-rigid registration is done by first doing a
% modality transformation (paints regions in image 1 with the intensity
% pallette of those regions in image2 and visa versa), and than
% using "normal" pixel based demon registration. See MutualTransform.m
%
%
% Example,
% % Read two greyscale images of Lena
% Imoving=imread('images/lenag1.png');
% Istatic=imread('images/lenag3.png');
%
% % Register the images
% [Ireg,Bx,By,Fx,Fy] = register_images(Imoving,Istatic,struct('Similarity','p'));
%
% % Show the registration result
% figure,
% subplot(2,2,1), imshow(Imoving); title('moving image');
% subplot(2,2,2), imshow(Istatic); title('static image');
% subplot(2,2,3), imshow(Ireg); title('registerd moving image');
% % Show also the static image transformed to the moving image
% Ireg2=movepixels(Istatic,Fx,Fy);
% subplot(2,2,4), imshow(Ireg2); title('registerd static image');
%
% % Show the transformation fields
% figure,
% subplot(2,2,1), imshow(Bx,[]); title('Backward Transf. in x direction');
% subplot(2,2,2), imshow(Fx,[]); title('Forward Transf. in x direction');
% subplot(2,2,3), imshow(By,[]); title('Backward Transf. in y direction');
% subplot(2,2,4), imshow(Fy,[]); title('Forward Transf. in y direction');
%
% % Calculate strain tensors
% E = strain(Fx,Fy);
% % Show the strain tensors
% figure,
% subplot(2,2,1), imshow(E(:,:,1,1),[]); title('Strain Tensors Exx');
% subplot(2,2,2), imshow(E(:,:,1,2),[]); title('Strain Tensors Exy');
% subplot(2,2,3), imshow(E(:,:,2,1),[]); title('Strain Tensors Eyx');
% subplot(2,2,4), imshow(E(:,:,2,2),[]); title('Strain Tensors Eyy');
%
% Example Multi-Modalities
% % Read two brain images
% Imoving=im2double(imread('images/brain_T1_wave.png'));
% Istatic=im2double(imread('images/brain_T2.png'));
%
% % Register the images
% [Ireg,Bx,By] = register_images(Imoving,Istatic,struct('SigmaFluid',4));
%
% figure,
% subplot(1,3,1), imshow(Imoving); title('moving image');
% subplot(1,3,2), imshow(Istatic); title('static image');
% subplot(1,3,3), imshow(Ireg); title('registerd moving image');
%
% % Read normal T1 image and transformation field
% Inormal=im2double(imread('images/brain_T1.png'));
% load('images/wave_field.mat');
%
% % Show the difference with ideal image
% figure, imshow(Imoving-Inormal,[-0.5 0.5]); title('unregistered')
% figure, imshow(Ireg-Inormal,[-0.5 0.5]); title('registered');
% disp(['pixel abs difference : ' num2str(sum(abs(Imoving(:)-Inormal(:))))])
% disp(['pixel abs difference : ' num2str(sum(abs(Imoving(:)-Ireg(:))))])
%
% % Show Warp field
% figure,
% subplot(2,2,1), imshow(BxNormal,[-20 20]); title('Bx Normal');
% subplot(2,2,2), imshow(Bx,[-20 20]); title('Bx');
% subplot(2,2,3), imshow(ByNormal,[-20 20]); title('By Normal');
% subplot(2,2,4), imshow(By,[-20 20]); title('By');
% Function is written by D.Kroon University of Twente (March 2009)
% add all needed function paths
try
functionname='register_images.m';
functiondir=which(functionname);
functiondir=functiondir(1:end-length(functionname));
addpath([functiondir '/functions'])
addpath([functiondir '/functions_affine'])
addpath([functiondir '/functions_nonrigid'])
catch me
disp(me.message);
end
% Disable warning
warning('off', 'MATLAB:maxNumCompThreads:Deprecated')
% Process inputs
defaultoptions=struct('Similarity',[],'Registration','NonRigid','MaxRef',[],'Verbose',2,'SigmaFluid',4,'Alpha',4,'SigmaDiff',1,'Interpolation','Linear');
if(~exist('Options','var')),
Options=defaultoptions;
else
tags = fieldnames(defaultoptions);
for i=1:length(tags)
if(~isfield(Options,tags{i})), Options.(tags{i})=defaultoptions.(tags{i}); end
end
if(length(tags)~=length(fieldnames(Options))),
warning('register_images:unknownoption','unknown options found');
end
end
% Set parameters
MaxRef=Options.MaxRef;
% Start time measurement
if(Options.Verbose>0), tic; end
% Store the class of the inputs
Iclass=class(Imoving);
% Convert the inputs to double
Imoving=im2double(Imoving);
Istatic=im2double(Istatic);
% Resize the moving image to fit the static image
if(sum(size(Istatic)-size(Imoving))~=0)
Imoving = imresize(Imoving,size(Istatic),'bicubic');
end
% Make smooth images for histogram and fast affine registration
ISmoving=imgaussian(Imoving,2.5,[10 10]);
ISstatic=imgaussian(Istatic,2.5,[10 10]);
% Detect if the mutual information or pixel distance can be used as
% similarity measure. By comparing the histograms.
if(isempty(Options.Similarity))
Hmoving= hist(ISmoving(:),60)./numel(Imoving);
Hstatic = hist(ISstatic(:),60)./numel(Istatic);
Hmoving(1)=0; Hstatic(1)=0;
if(sum(log(abs(Hmoving-Hstatic)+1))>0.3),
Options.Similarity='m';
if(Options.Verbose>0), disp('Multi Modalities, Mutual information is used'); drawnow; end
else
Options.Similarity='p';
if(Options.Verbose>0), disp('Same Modalities, Pixel Distance is used'); drawnow; end
end
end
if(Options.Similarity(1)=='p'), type_affine='sd'; else type_affine='mi'; end
% Register the moving image affine to the static image
% Affine register the smoothed images to get the registration parameters
if(strcmpi(Options.Registration(1),'R'))
if(Options.Verbose>0), disp('Start Rigid registration'); drawnow; end
% Parameter scaling of the Translation and Rotation
scale=[1 1 1];
% Set initial affine parameters
x=[0 0 0];
elseif(strcmpi(Options.Registration(1),'S'))
if(Options.Verbose>0), disp('Start Similarity Registration'); drawnow; end
% Parameter scaling of Translation, Rotation, Resize
scale = [ 1 1 1 .01 ];
% Set initial similarity parameters
x = [ 0 0 0 100 ];
elseif(strcmpi(Options.Registration(1),'A'))
if(Options.Verbose>0), disp('Start Affine registration'); drawnow; end
% Parameter scaling of the Translation, Rotation, Resize and Shear
scale=[1 1 1 0.01 0.01 1e-4 1e-4];
% Set initial affine parameters
x=[0 0 0 100 100 0 0];
elseif(strcmpi(Options.Registration(1),'N'))
if(Options.Verbose>0), disp('Start Affine part of Non-Rigid registration'); drawnow; end
% Parameter scaling of the Translation, Rotation, Resize and Shear
scale=[1 1 1 0.01 0.01 1e-4 1e-4];
% Set initial affine parameters
x=[0 0 0 100 100 0 0];
else
warning('register_images:unknownoption','unknown registration method');
end
for refine_itt=1:2
if(refine_itt==2)
ISmoving=Imoving; ISstatic=Istatic;
end
% Use struct because expanded optimset is part of the Optimization Toolbox.
optim=struct('GradObj','off','GoalsExactAchieve',1,'Display','off','MaxIter',100,'MaxFunEvals',1000,'TolFun',1e-14,'DiffMinChange',1e-6);
if(Options.Verbose>0), optim.Display='iter'; end
x=fminlbfgs(@(x)affine_registration_error(x,scale,ISmoving,ISstatic,type_affine),x,optim);
end
% Scale the translation, resize and rotation parameters to the real values
x=x.*scale;
transformedMatrix = x;
if(strcmpi(Options.Registration(1),'R'))
% Make the rigid transformation matrix
M=make_transformation_matrix(x(1:2),x(3));
elseif(strcmpi(Options.Registration(1),'S'))
M = make_transformation_matrix(x(1:2),x(3),[x(4) x(4)]);
else
% Make the affine transformation matrix
M=make_transformation_matrix(x(1:2),x(3),x(4:5));
end
% Make center of the image transformation coordinates 0,0
[x,y]=ndgrid(0:(size(Imoving,1)-1),0:(size(Imoving,2)-1));
xd=x-(size(Imoving,1)/2); yd=y-(size(Imoving,2)/2);
% Calculate the backwards transformation fields
Bx = ((size(Imoving,1)/2) + M(1,1) * xd + M(1,2) *yd + M(1,3) * 1)-x;
By = ((size(Imoving,2)/2) + M(2,1) * xd + M(2,2) *yd + M(2,3) * 1)-y;
% Initialize the modality transformed image variables
M_TF=[]; F_TF=[];
% The nonrigid part of the registration
if(strcmpi(Options.Registration(1),'N'))
% Demon registration parameters
refinements=floor(log2(min(size(Imoving))/16));
if(refinements>MaxRef), refinements=MaxRef; end
parameters.sigma_diff=Options.SigmaDiff;
% Non-rigid registration
if(Options.Verbose>0), disp('Start non-rigid demon registration'); drawnow; end
% Do every refinements step twice if modality transformation enabled
if(Options.Similarity(1)=='m'), loop=2; else loop=1; end
% Loop trough all refinements steps.
for j=0:refinements
for l=1:loop
% Set scaling parameters.resizepercentageentage
resizepercentage=1/2^(refinements-j);
if(resizepercentage>1), resizepercentage=1; end
parameters.alpha=Options.Alpha*sqrt(resizepercentage);
parameters.sigma_fluid=Options.SigmaFluid;
if(Options.Verbose>0), disp(['Scaling resizepercentageentage : ' num2str(resizepercentage)]), end
% Incase of multiple modalities, transform both images to their
% opposite modalities.
if(Options.Similarity(1)=='m')
if(Options.Verbose>0), disp('Start modality transformation'); drawnow; end
Bx_large=imresize(Bx,size(Imoving),'bicubic')*(size(Imoving,1)/size(Bx,1));
By_large=imresize(By,size(Imoving),'bicubic')*(size(Imoving,2)/size(By,2));
[Imoving_TF,Istatic_TF]=MutualTransform(Imoving,Istatic,15*sqrt(1/resizepercentage),4,Bx_large,By_large);
if(Options.Verbose>0), disp('Finished modality transformation'); drawnow; end
end
sigma = 0.3/resizepercentage;
% Set and resize the moving image and static image
M=imresize(imgaussian(Imoving,sigma,[sigma*6 sigma*6]),resizepercentage,'bicubic');
F=imresize(imgaussian(Istatic,sigma,[sigma*6 sigma*6]),resizepercentage,'bicubic');
% Resize the modality transformed images
if(Options.Similarity(1)=='m')
M_TF=imresize(imgaussian(Imoving_TF,sigma,[sigma*6 sigma*6]),resizepercentage,'bicubic');
F_TF=imresize(imgaussian(Istatic_TF,sigma,[sigma*6 sigma*6]),resizepercentage,'bicubic');
end
% Resize the transformation field to current image size
Bx=imresize(Bx,size(M),'bicubic')*(size(M,1)/size(Bx,1));
By=imresize(By,size(M),'bicubic')*(size(M,2)/size(By,2));
% Put transformation fields in x and y direction in one variable
B=zeros([size(M) 2]); B(:,:,1)=Bx; B(:,:,2)=By;
% Store the dimensions of transformation field, and make a long vector from T
sizes=size(B); B=B(:);
% Parameters
options.sigma_fluid=parameters.sigma_fluid;
options.sigma_diff=parameters.sigma_diff;
options.alpha=parameters.alpha;
options.interpolation=Options.Interpolation;
% Optimizer parameters
optim=struct('Display','off','StoreN',10,'GoalsExactAchieve',0,'HessUpdate','lbfgs','GradObj','on','OutputFcn', @store_transf,'MaxIter',200,'TolFun',1e-14,'DiffMinChange',1e-5);
if(l==loop),
optim.TolX = 0.02;
else
optim.TolX = 0.1;
end
if(Options.Verbose>1), optim.Display='iter'; end
% Start the demon energy registration optimizer
B=fminlbfgs(@(x)demons_energy(M,F,M_TF,F_TF,x,sizes,options),B,optim);
% Reshape B from a vector to an x and y transformation field
B=reshape(B,sizes);
Bx=B(:,:,1); By=B(:,:,2);
end
end
% Scale everything back if not already
if(resizepercentage~=1)
Bx=imresize(Bx,size(Imoving),'bicubic')*(size(Imoving,1)/size(Bx,1));
By=imresize(By,size(Imoving),'bicubic')*(size(Imoving,2)/size(By,2));
end
end
% Transform the input image
Ireg=movepixels(Imoving,Bx,By,[], 3);
if ( nargout>3 )
% Make the forward transformation fields from the backwards
[Fx,Fy]=backwards2forwards(Bx,By);
end
% Set the class of output to input class
if(strcmpi(Iclass,'uint8')), Ireg=im2uint8(Ireg); end
if(strcmpi(Iclass,'uint16')), Ireg=im2uint16(Ireg); end
if(strcmpi(Iclass,'uint32')), Ireg=im2uint32(Ireg); end
if(strcmpi(Iclass,'int8')), Ireg=im2int8(Ireg); end
if(strcmpi(Iclass,'int16')), Ireg=im2int16(Ireg); end
if(strcmpi(Iclass,'int32')), Ireg=im2int32(Ireg); end
if(strcmpi(Iclass,'single')), Ireg=im2single(Ireg); end
% End time measurement
if(Options.Verbose>0), toc, end