How to improve imregister for rigid 3D image registration of very small images?

I have 3D image stacks of a dendrite over time and want to register every image stack to the first one of the series. I found the function imregister and want to use it, but the resulting alignment is not very well. Please find two example images and a registration attempt attached. The images are created using two-photon microscopy.

I played around with the parameters but I feel like I can do this forever and I still have Mean-Squared Errors around 380. The images are also quite small, voxel size of 0.007 microns in X&Y and 0.6 microns in Z; image size = 512x512x22. The 3D array has values from 0-255.

So I wanted to ask for some experience/tips in choosing parameters especially for such an image size, or if you spot any error in my code. How should I proceed?

Here is my code

fixed = get3DImage(<fixedImage>);   %returns 3D array containing my image, uint8
moving = get3DImage(<movingImage>);
Rfixed = imref3d(size(fixed),0.000007,0.000007,0.0006);   %Pixel extent in world parameters
Rmoving = imref3d(size(moving),0.000007,0.000007,0.0006);

[optimizer, metric] = imregconfig('monomodal');
optimizer.MaximumStepLength = 6.25e-05; 
optimizer.MinimumStepLength = 0.9e-05;
optimizer.MaximumIterations = 1000;

movingRegistered = imregister(moving,Rmoving, fixed, Rfixed, 'rigid', optimizer, metric);

Any help appreciated! Thanks

Answer 1

[Matt J]

I recommend that you use imregtform instead, since it will output the transformation as a result, instead of just the registered image. I also recommend that you use the 'InitialTransformation' option to initialize the registration with the immediately preceding registration result in the series. Assuming the motion is gradual over time, this will make it easier for the solver to deal with very large displacements relative to the fixed image.

Finally, are you sure that a 'rigid' motion model is a good one for dendrites? I thought the only biological tissue that moved rigidly was bone, but then again I'm not a biologist.

[Hello Matlab]

Ok, but the Algorithm is the same right? Imregister calls imregtform.

The imaged object gets removed and put back in the camera after some days. So the motion is not gradual but independent and due to the fact that the object is always slightly different oriented in the camera. So I cannot conclude the next orientation based from the previous image.

The rigid model is a specification from the lab, I am not a biologist either.

[Matt J]

The algorithm is the same....

Well, they didn't make it easy on you. I would try the multimodal regconfig option, at least as a way of getting a decent initial guess of the transform. You could also try changing just the optimizer to the Evolutionary algorithm.

[Hello Matlab]

The combination Evolutionary optimizer plus Mean-Square metric looks a bit better. Then I tried to run imregtform a second time with the result from the first registration as initial transformation. Like so:

optimizer = registration.optimizer.OnePlusOneEvolutionary;
optimizer.InitialRadius = 6.25e-05;
optimizer.MaximumIterations = 300;
metric = registration.metric.MeanSquares;

T = imregtform(moving,Rmoving, fixed, Rfixed, 'rigid', optimizer, metric);
T2 = imregtform(moving,Rmoving,fixed,Rfixed,'rigid',optimizer,metric,'DisplayOptimization',true,'InitialTransformation',T);
movingRegistered = imwarp(moving,T2);

But this didn't bring a good result, error around 380. Any idea what to do next?

[Matt J]

It's not possible to look at the regmex implementation. It's a dll. But looking at your images, it is far more likely that the problem is in the initial tform guess, or the assumption of a rigid object, rather than the algorithm implementation. Registration algorithms typically depend on partial overlap between structures in the fixed and moving images to get a foothold in the search. Considering how spatially varying your image structures are, I'm betting you have very little overlap (I see none in the images you've posted). Without a good initial guess, this makes the solution very hard for imregtform to find.

I think the place to start is to try to manually obtain a good initial guess of the rototranslation and use it to test the assumption of rigid motion.

If you manually pick out some corresponding landmarks in each image, you can register their coordinates using ABSOR. You can then convert the output of absor to an initial tform guess by doing

    tform = affine3d(reg.M.');

If this doesn't get you into a regime where imregtform performs well, the chances are pretty good that a rigid registration model doesn't fit very well. Incidentally, if the registration error given by absor

[reg,~,errors]=absor(A,B)

is large, that would also be a sign that a rigid model is bad.

[Hello Matlab]

Yeah, I figured that.. Didn't want to accuse the algorithm but wanted to look how it is implement out of interest.

Cool, I tried this and it looks promising! But weirdly imregister always diverged with this as InitialTransformation, no matter how small the parameters. But when I imwarp() the image with it, then the rotation is very accurate, the translation not yet optimal. Those are the errors: errlsq: 15.8247, errmax: 12.6509. I have no clue, why the registration still won't work..

Also, picking the landmark points is really hard and impractical the way I did it. Is there a smoother way for finding 3d coordinates? I tried ginput and rather estimated the z-page. Maybe this is also why the translation is not 100%.

Do you have a hint here?

[Matt J]

errlsq: 15.8247, errmax: 12.6509

What units are these? Millimeters? Voxels? And are your voxels isotropic (equal lengths in x,y, and z)? ABSOR assumes you have common units for the x,y, and z coordinates.

Also, it would be informative to examine the distances between your landmarks in both images. If the motion is rigid, the inter-landmark distances should be preserved.

Likewise with angles. If the triangle formed by any three landmarks have nearly the same angles in both images, it is evidence that the transformation is rigid. And if not, that it is not rigid.

Also, picking the landmark points is really hard

I didn't intend for the final solution to be a landmark-based registration. As far as I was concerned, this was just for the purposes of analysis and setting performance baselines.

There is no generic way of doing automatic landmark detection that I know of, but if you think you can design something based on your insights into this class of objects, it could be something to pursue. regionprops3() gives a variety of segmentation tools.

[Hello Matlab]

One voxel is of 0.007x0.007x0.6 microns, or 0.000007x0.000007x0.0006 mm. So no, they are anisotropic.

Okay, I will check this.

About the landmarks, I didn't mean a method for automatic detection, but simple clicking. Something like ginput, but for 3d or so. Because I had to estimate the z-plane of the voxels, while x and y were precise thanks to ginput.

I am thinking, isn't it possible to use results from regionprops as initial Transformation parameters? For example the difference of the centroids in x/y/z as Init for the translation etc.?

[Matt J]

Something like ginput, but for 3d or so.

You could try computing maximum intensity projections

figure(1);
imshow(squeeze(   max(yourImage,[],1)   )); %MIP in y,z

figure(2);
imshow(squeeze(   max(yourImage,[],2)   )); %MIP in x,z

and use ginput to explore these.

I am thinking, isn't it possible to use results from regionprops as initial Transformation parameters? For example the difference of the centroids in x/y/z as Init for the translation etc.?

Yes, but I think the WeightedCentroid would be better. And you could similarly analyze the EigenVectors to get the rotation, although there again I think a weighted eigenvector computation would be better, maybe using pca().

[Hello Matlab]

Cool, I will try this. But: "And you could similarly analyze the EigenVectors to get the rotation" - How do I get the 3 rotation angles of those? And how do they differ from the 'Orientation' option?

[Matt J]

The three Eigenvectors form the rows of a rotation matrix whose Euler angle (or Tait-Bryan angle?) decomposition gives Orientation. How that decomposition works is something you probably understand better than anyone else, since you explained how it works here.

But I don't see why you would need to analyze things in terms of angles. If E1 and E2 are the 3x3 eigenvector matrices of Image 1 and Image 2, then they will be related by a rotation matrix R

     E1*R=E2

which you can just solve for by matrix algebra,

   R=E1.'*E2

[Hello Matlab]

Aaah, now I see that connection! One matrix multiplication is indeed easier than rotating manually by the angles. Especially when axes are unclear.

Now I use a transformation object (rotation matrix from eigenvectors and weighted centroids diff for translation) as Initial Transformation. Mean-Square Errors are still pretty high, though, and the registered images not euqally aligned. This will take a little longer

You helped me a lot, Matt, thank you.