Closed subhajitchatu closed 7 years ago
Hi subhajitchatu,
As you seem to be starting to use Python I would strongly advise that you use Python3 (v.3.6 is latest) instead of Python2. And especially you can start with a "scientific" distribution, I personally like Anaconda Python very much, as it bundles a lot of packages from scratch together with a very nice package manager (conda). The Spyder IDE can also be of an help, with the extremely useful ctrl+i shortcut to render any object docstring inplace (and nicely).
Regarding this part, I believe @kaspermarstal will be the best to answer.
Regarding point 2. as it seems you're migrating from Matlab to Python you should probably have a look at matplotlib, scipy, scikit-image and numpy, which are libraries often cited for providing some functions and objects with a syntax close to Matlab. All these libraries have wonderfull documentations so you will find your way using it easily. There is a lot of articles/blog posts to help this transition, you can start with this one or look yourself for others.
Regarding your specific question about imshowpair(), I did not found a similar function, but it happens that I coded one for my personal use, so you can have a look if it helps.
import numpy as np
import matplotlib.pyplot as plt
import SimpleITK as sitk
def overlay(img1, img2, title=None, interpolation=None, sizeThreshold=128):
"""
Description
-----------
Displays an overlay of two 2D images using matplotlib.pyplot.imshow().
Args
----
img1 : np.array or sitk.Image
image to be displayed
img2 : np.array or sitk.Image
image to be displayed
Optional
--------
title : string
string used as image title
(default None)
interpolation : string
desired option for interpolation among matplotlib.pyplot.imshow options
(default nearest)
sizeThreshold : integer
size under which interpolation is automatically set to 'nearest' if all
dimensions of img1 and img2 are below
(default 128)
"""
# Check for type of images and convert to np.array
if isinstance(img1, sitk.Image):
img1 = sitk.GetArrayFromImage(img1)
if isinstance(img2, sitk.Image):
img2 = sitk.GetArrayFromImage(img2)
if type(img1) is not type(img2) is not np.ndarray:
raise NotImplementedError('Please provide images as np.array or '
'sitk.Image.')
# Check for size of images
if not img1.ndim == img2.ndim == 2:
raise NotImplementedError('Only supports 2D images.')
if interpolation:
plt.imshow(img1, cmap='summer', interpolation=interpolation)
plt.imshow(img2, cmap='autumn', alpha=0.5, interpolation=interpolation)
elif max(max(img1.shape), max(img2.shape)) > sizeThreshold:
plt.imshow(img1, cmap='summer')
plt.imshow(img2, cmap='autumn', alpha=0.5)
else:
plt.imshow(img1, cmap='summer', interpolation='nearest')
plt.imshow(img2, cmap='autumn', alpha=0.5, interpolation='nearest')
plt.title(title)
plt.axis('off')
plt.show()
You can feed it either SimpleITK.Image() or numpy.array(), so it is quicker to use it with SimpleElastiX. It only supports a "blend" mode, but if you want a difference image it is very easy to get using numpy.
# Something like that
diffImage = img1 - img2
Considering this part, I believe ITK has features like Hausdroff distance, so SimpleITK should have it too, just have a look in their doc, or explore SimpleITK from Python.
Hope it helps.
@PertuyF
Hi,
Thanks a lot for your valuable reply.
I have coded using matplotlib for displaying images. But I thought SimpleITK may provide this too.
I have searched for Huaudroff Distance in python SimpleITK codes, didn't find anything.
There are several plotting libraries in Python, with different aims, but matplotlib
certainly is the most used for simple local display or plotting.
Actually SimpleITK is not meant for visualisation. Although it does provide a show() function, it requires external software for the display (ImageJ by default). See there for resources about SimpleITK.
You might be especially interested into notebooks 01, 02, 03 and 10 for images insights (subtilities between sitk.Image and np.array, visualisation). Be aware that there is a mistake in these documents, the author sometimes calls imshow(), which is either pylab's or matplotlib's (doesn't matter), but does not import it properly. You should do either:
import pylab
pylab.imshow(yourImage)
or
from pylab import imshow()
imshow(yourImage)
About Hausdroff distance and such tools, I found this tutorial. I find the code repelling (too much text) but it follows the ITK "explicit everything has much as you can" verbose style, and you should find relevant information there.
Hi @subhajitchatu,
Parameter maps are only available for the elastix-based registration methods. The ITK methods are also available but use a different approach to setup the registration. Check out http://insightsoftwareconsortium.github.io/SimpleITK-Notebooks/60_Registration_Introduction.html for an introduction. You can find examples for specific methods at http://insightsoftwareconsortium.github.io/SimpleITK-Notebooks/60_Registration_Introduction.html.
Looks like @PertuyF has answered you other questions :D
@kaspermarstal First of all a big Thanks for building such a useful platform.
A. I was looking for deformable registrations using SimpleITK. The link u shared are not telling about any deforamable registration, even sitkImageRegistrationMethod.h, there is no FEM or thin plate spline related functions in that header fine functions.
Anyway in RegistrationITKV4 there are FEM based ( itkFEMRegistrationFilter.h) deformable registration and thin plate spline code is available, I was trying to figure out how to start using them in SimpleITK or python based platform. Please throw some lights.
Demon registration is available in SimpleITK.
B. Anyway** there are no hauzdroff distance API for accuracy estimation for Registration. Is there any other way to look into these things. Please suggest.
@PertuyF Thanks for sharing the link for accuracy estimation. http://insightsoftwareconsortium.github.io/SimpleITK-Notebooks/34_Segmentation_Evaluation.html helped me for accuracy estimation.
But I am now looking for way to use deformable registration methods.
Thanks
Then you should be looking at ElastiX registrations using BSplineTransform and SplineKernelTransform (thin plate), these allow to perform deformable registration. If needed, have a look at ElastiX manual for details.
Regarding registration accuracy, could you explain what you aim to do? In my experience, assessing the quality of a deformable registration if quite subjective, and however scoring is somewhat doable for binary images, it is no very reliable for plain images.
@PertuyF Oh Thanks. As @kaspermarstal suggested, all Elastix Registration methods are included in SimpleITK , so there is no issue to use BSplineTransform and SplineKernelTransform.
But I was thinking how one can use FEM or meshing using some element(triangle,quadrilateral shape) in this settings. Right now there are no open source stable software (commercial s/w: ScanIP) which can do meshing using particular geometric shaped element on DCM image stacks and using some deformation model(Neo Hookean or viscous model etc) analysis the displacements. (ANSYS).
Requirement: There are several nonrigid registration methods which are mainly used to correct local nonrigid elastix motions. Now FEM or any deformable registration method can simulate displacement on each node of a mesh and by warping them back one can align moving image to fixed image . (another advantages are these are not intensity based, so it can be used in multimodal settings). Now after registering two image, If I want to check how good it corrects local motions,(for eg. in a 2-3mm tumor alignment), we need some mechanism for visualization and quantification. So I was searching for some methods in SimpleITK.
P.S. all these requirements are for DCM images.
@subhajitchatu I would suggest that you go ask these questions to the ElastiX mailing list, where you may find more complete answers. You'll be able to search the mailing for previous discussions, and I do not doubt you will find plenty of usefull information about ElastiX capabilities there.
@PertuyF Thanks, but ElastiX is not SimpleElastix right? so suppose even if Elastix supports one method , does that mean SimpleElastix will do the same?
@subhajitchatu I am pretty sure SimpleElastiX does support most of ElastiX methods, if not all. Anyway, active people on the mailing list will likely be able to tell you whether a feature you're asking for is missing or not.
Am I right @kaspermarstal ?
@pertuyF as you suggested and I have also found SplineKernelTransfrom in elastix, can you please elaborate how to use it in Simplelastix in Python settings. I was unable to find it in simpleITK class , I might be wrong, but I need your suggestions . I did go through the code and found bspline but SplineKernelTransfrom was only in elastix , how can I use or access API for thin plate. Thanks
@subhajitchatu with ElastiX, and SimpleElastiX alike, you must tune your registration settings by modifying a parameterMap
. You can read about the basics in SimpleElastiX doc.
When you build SimpleElastiX, it will wrap all features the reference ElastiX install was built with, and that you selected in CMake. The fastest way to see if it is available in your build is to try a registration using it.
Basically, my approach is that if I see something I need/want in ElastiX documentation, I just try it in SimpleElastiX and see if it works. Otherwise I ask Kasper or the ElastiX mailing list.
For your case the first step would be to load or create a parameterMap and to change the transform to SplineKernelTransfrom in order to test it. As SplineKernelTransform seems to require specific parameters, I recommand that you carefully look at this document.
import SimpleITK as sitk
selx = sitk.SimpleElastix()
selx.SetLogToConsoleOn()
P = selx.GetDefaultParameterMap("bspline")
P['Transform'] = ['SplineKernelTransform'] # Edit PMap by key accessing
# Change SplineKernel, SplineRelaxationFactor or SplinePoissonRatio if needed
selx.SetParameterMap(P)
selx.SetFixedImage(yourFixedImage)
selx.SetFixedPointSet(yourFixedLandMarks) # Mandatory with SplineKernelTransform
selx.SetMovingImage(yourMovingImage)
selx.Execute()
At this step you should see directly into your python console if anything goes wrong. As I never tried this specific transform, and I cannot run selx right now, I have no idea if it will work, or if my piece of code is correct or not.
As I feel that you do not seem to like the BSplineTransform for now, I would just indicate that by tweaking the parameterMap I was able to compute very nice population averages from µCT datasets just by using EulerTransform, AffineTransform and BSplineTransfom. I am not very well aware of all technologies available in this field, but in my experience ElastiX's BSpline should be convenient for most registrations.
@PertuyF yes you are correct :) SimpleElastix supports all the features of the elastix C++ api (@subhajitchatu that means some component-specific settings, like the options to TransformRigidityPenalty, are not used, but all registrations methods are in there).
@subhajitchatu There is extensive documentation on how to use SimpleElastix but it will admittedly require a little effort on your part to read up on it. You can see how to setup and use parameter maps at https://simpleelastix.readthedocs.io/ParameterMaps.html, you can find information on the SplineKernelTransform in the manual http://elastix.isi.uu.nl/download/elastix_manual_v4.8.pdf on page 10 and you can find info on the SpineKernelTransform paramerters in the api reference documenation here http://elastix.isi.uu.nl/doxygen/classelastix_1_1SplineKernelTransform.html.
@PertuyF Thank you so much. If any feature is in Elastix but not in simpleElastix, can it be possible to add in Make file to get the feature or have to code it separately and build new SimpleITK.so. I mean is this platform generic with Elastix? You used EulerTransform, AffineTransform and BSplineTransfom altogether right? (rigid followed by non-rigid). So nice to hear from you @kaspermarstal . I have gone through the code, if I am not wrong GetDefaultParameterMap() must have an Enum list for which transformation it supports. Right? parameterMap['Transform'] = ['BSplineTransform'], But I didn't find SplineKernelTransform or FEM based algorithms, that why I asked my doubts. But as the code is embedded , So I think I am missing some FLAGS in MAKE file or just mis-conception.
@subhajitchatu Yes I usually perform Rigid registration, then Affine, then several rounds of BSpline with decreasing values for FinalGridSpacingInVoxels. As Kasper said, it will take you some time to read the doc, but I think the results are worth it.
@subhajitchatu your welcome!
You can always change the default parameter maps to suit your needs, so you can do fx
parameterMap = sitk.GetDefaultParameterMap('bspline')
parameterMap['Transform'] = ['SplineKernelTransform']
That being said, there should be a SplineKernelTransform default parameter map in there, as seen here https://simpleelastix.readthedocs.io/PointBasedRegistration.html.
Note that you have to provide fixed point set as described in the manual. This is done with the SetFixedPointSetFileName()
in SimpleElastix. The format is described here https://simpleelastix.readthedocs.io/PointBasedRegistration.html.
@PertuyF Just an out of context question, In your functions of overlay, you have used sequential monochromatic colormaps, right?
plt.imshow(img1, cmap='summer')
plt.imshow(img2, cmap='autumn', alpha=0.5)
This brings color to background also, to make it transparent or tune accordingly, changing cmap or alpha didn't work. How we can bring that feature, it will show the color overlay more accurately, one can easily identify how well registration works. Now it is overshadowed by the background color itself.
@subhajitchatu The colormaps aren't monochromatic, they are sequential. I think your problem comes from the fact that 'autumn' and 'summer' do not faint to black for the lower intensity levels. You can change them by choosing among matplotlib's colormap collection. My idea there was to use two colormaps that do no share any color. Maybe the new 'viridis' can be used with 'magma' to correct this problem a bit. Just give it a try.
Note that if you're willing to, you can create a colormap yourself. Just look at the details of the matplotlib.colors.Colormap
class for this.
@PertuyF @kaspermarstal There are two different kinds of APIs available for registration in SimpleElastix and SimpleITK.
simpleITK: R = sitk.ImageRegistrationMethod() R.SetMetricAsCorrelation()
and then setting Optimizer , we can achieve registration results.
If we use SimpleElatix, SimpleElastix = sitk.SimpleElastix()
again using setParameterMap() API we can implement different registration method.
In my understanding SimpleElastix is another higher level of API abstraction to simpleITK.
are they same? which one is better to use? Is there any guidelines, which one should we use. In most of the examples in SimpleElastix, SimpleITK APIs are used.
@subhajitchatu ITK is a C++ toolbox for image processing, and SimpleITK offers an interface to ITK in several languages. ElastiX itself is an independent toolbox for image registration, which is based on ITK. What @kaspermarstal did with SimpleElastiX is a fork of SimpleITK that integrates an interface to ElastiX. So basically in SimpleElastiX you can access ITK's registration tools as well as ElastiX's tools, which are different. However, SimpleElastiX will use some code from ITK as ElastiX is based on it. That's for who's who.
Now for who should you use, it's up to you, try it both and choose.
@PertuyF Thanks a lot! SimpleElastix = union(simpleITK,ElastiX) right? So functionality wise SimpleElastix can provide a lot depending upon the functions ElastiX is providing which are better and extra than simpleITK. But if simpleElastix is using SimpleITK, then in this case, it might end up with some extra function calls. That's it!(simpleElastix-->simpleITK--->ITK)
@PertuyF
Hi, I have registered two images, let's say fixed and moving are Registered. After registration I want to measure overlap ratio etc. The overlap_measures_filter.Execute(fixed, moving) and hausdroff_measures_filter.Execute() returns huge false negatives and positives. It seems it is not working. If you have used it before please let me know how exactly can I use these functions for registration evaluation?
@subhajitchatu What are you comparing exactly? These methods are meant to compare segmented images if I am not mistaken. I never used it (only Dice score one or twice) as I rarely registered segmented images. If I really needed to evaluate registration I'd use my eyes to judge, or maybe a difference image sometimes.
You'll have to dig ITK's doc and SimpleITK's notebooks for more information I'm afraid.
@PertuyF I was comparing registration methods accuracy for a type of image. Diff image is not sufficient.anyway thanks. Are there any way to show 3D/4D image in ur display function.
@subhajitchatu Unfortunately this function is only intended for facilitating 2D display. For 3D display from code I think it would be better to use another library dedicated to 3D, VTK is one (C++ and wrapped to several languages, among which python), MayaVi or VisPy are others (only for Python, and OpenCL may be required at some point).
Anyway, manipulating 3D rendering from code will most likely require a lot of time and knowledge. Especially for CT-like data, which is most often better visualised using volume renders, and not surface renders. I prefer using some external viewer for this part. Among commercial solutions Amira/Avizo is very efficient, and Osirix is renowned among mac users. Among free solutions 3DSlicer is open-source and quite good, MeVisLab well known but I used it very rarely, and I recently discovered TomViz which seems to be quite interesting (visualisation only).
But for now I may stop talking about unrelated materials in a SimpleElastiX issue, which is not the place for this.
@kaspermarstal
I have question regarding SimpleElastix API. After using registration like rigid, affine or non-rigid with some similarity metric and optimizer. I want to fetch GetMetricValue() as it is in SimpleITK, ierations it took to terminate or time taken to register. It would be nice to get all information after registering two images using some method. It will help us to see which similarity measure works well as well as which optimizer is working best. Many of these are available in SimpleITK but can't find it in SimpleElastix.
Any Suggestion?
Thanks
@subhajitchatu things like the metric value are not available in code from elastix unfortunately. But perhaps you can SetLogToFile('filename') and parse the log after registration has completed? The log will have all the information about metric values.
However, allow me to humbly suggest to use LabelShapeStatisticsImageFilter(), LabelOverlapMeasuresImageFilter() or just plain visual inspection to assess registration quality. The raw metric value is not well suited for this in my experience.
@subhajitchatu @kaspermarstal
Concerning final metric value, I was able to extract it from elastix.log
using this regex: 'Final metric value = (?P<value>[+-.0-9]{9})'
. Maybe it will help.
Hi,
How can I use "MultiMetricMultiResolutionRegistration" with some weights in SimpleElastix?
How can I search initial guess randomly for an optimization method? As we know our objective function can be non-convex also(I don't know when or why? may be Mutual information is a non-convex? SSD is convex always right? ), which might have local minima, I just want to avoid this using some exhaustive or random search of initial point? Is there any other efficient way to get rid of this local minima? I have tried with some different guess but couldn't able to automate that.
Please let me know if Multiple metric with some weight function is allowed in SimpleElastix or not? @kaspermarstal @PertuyF
Thanks
@PertuyF @kaspermarstal
As you people suggested , I tried to find out metric value from simple elastix logs and it looks like
=============== end of TransformParameterFile ===============
Creating the TransformParameterFile took 0.45s
Registration result checksum: 268701954
Final metric value = -1.392004
The metric was stochastic gradient descent and the negative reminds of a minimization function. Is that right? Is it looking for a minimized value? I mean more less the value more better registration, can we conclude that?
Thanks
@subhajitchatu
Stochastic gradient descent is the optimizer, the type of metric should be "AdvancedMattesMutualInformation", or "AdvancedMeanSquares", or any other metric available in ElastiX.
Yes, the idea is that the transform parameters are optimized according to the metric value, but I am not sure that the metric value can be used to asses the quality of the registration because:
Note that with some metric, optimization can go toward a higher value (saw that in ElastiX but I can't remember which one was concerned)
@PertuyF
Yes it will change from one image pair to other. But lets suppose we fix image pair and metric, now we want to see which optimizer gives better result. Then I think metric value is somehow helpful. But it is not a minimizing function as I have seen the final metric value ,closer to zero gives better result. (-0.067 is better than - 1.66).
@PertuyF Can you please throw some light on my previous question on MultiMetricMultiResolutionRegistration, if anyone wants to use multiple metric together with some weightage, How to achieve it in SimpleElastix? Actually it is in Elastix documentation advance. So I thought it may possible here also.
and also on searching the initial guess efficiently. Please help if you have idea on these.
Thanks
@subhajitchatu This is getting too technical for me, I'm afraid I can't give you answers here. I do not understand these aspects well enough, sorry.
Regarding the weights and MultiMetricMultiResolutionRegistration, what parameter do you want to use the weights on?
@PertuyF
Ok, Thanks so much. @kaspermarstal might help on metric value interpretation question. Please help!
MultiMetricMultiResolutionRegistration 's usage is not clear to me. It is said that multiple metric can be used for different resolution. By default weight is (1/#metric), but we can change it to like for two metric like MI and CC (0.8 , 0.2). Have you used that before? I don't know how it helps in registration!
According to my understanding if we use MI for lower resolution, then it's result might be a good starting for the next resolution where we will use both CC and MI. Again it is not clear why it is helpful, may be get rid of local minima or saddle point? OR some other reason? People also use multi-image , multi-metric multi resolution registration , but do you know how will it help in getting good results?
In my cases DCE-MRI has 40 different stack of images, that means if we take one fixed image then 39 moving images, so multi-image multi-metric might be helpful for me. That's why I just want try them in my case. if SimpleElastix supports that will be great.
Thanks
@subhajitchatu It is supported, just add the correct parameters to your parameterFile as described in ElastiX doc. Here you should add something like:
P = youParameterFileObject
P['Metric0Weight'] = ['1', '0.5', '0.8']
P['Metric1Weight'] = ['0', '0.5', '0.8']
Note that the actual weight will be normalized between the metrics, so 0.5 / 0.5 is the same as 0.8 / 0.8 or 1 / 1. You can either adjust the weights for fine tuning, or use the boolean "Metric0Use" to deactivate a metric for a specific resolution.
It can definitely be helpful in some cases, but it is totally case-dependent so you will have to try on yourself. Expect to spend a lot of time on it, it is not easy to optimize such parameters.
Hi @PertuyF
Thanks so much. But it is not clear what will be the advantages of using it? The purpose of using it is getting rid of local minima problem or something else? Just a out of context question, Do you have idea on SimpleITK registration framework? SimpleElastix might also have features to use demon registrations (diffeomorphic), but just in case you have used SimpleITK python framework, please throw some lights on the below question https://github.com/InsightSoftwareConsortium/SimpleITK-Notebooks/issues/132
Thanks a lot
Hi @subhajitchatu,
On the advantages of weighting metrics, I did see some improvements in previous works, but the main one was using RigidityPenalty
, which has the specific purpose of limiting deformations during bspline registration. It will probably not help you there.
A little off-topic digression regarding your specific problem of local minima. A local minimum is a value (for simplicity, lets assume we only optimize one value) that is not the global minimum, but which "traps" the optimizer. These minima depend on the loss function used in the optimization (in the case of ElastiX, the metric). In this specific case the loss is a representation of the difference between moving and fixed images, and the optimizer tries to reduce this difference as much as possible. The process is iterative and often uses the Standard Gradient Descent (SGD) algorithm and variations around SGD.
To avoid a local minimum you can either:
Using multiple metrics obviously refers to option 1, and using weighted metrics will allow you to finely compose the contribution of each metric to your final loss function. However, the idea here is that using multiple metrics that work in a similar way will not help much, as they may "share" a local minimum. So a little bit of thinking can help to be efficient here. Regarding option 2, there is a nice review of variations of SGD optimizers (with some recent development linked to the deeplearning field) which can help you understand how it works. However I am not aware of such variants implemented for ElastiX, so this is only for your knowledge.
About your other question, I have no real experience with ITK's internal registration framework, so I won't be able to help there, sorry for that.
Note that I am not an expert in SGD, and I don't know which implementation of SGD is AdaptativeGradientDescentOptimizer
used in ElastiX. These informations are just to help you figure out some options to look in if you want to solve your problem.
Thanks so much @PertuyF
Can you tell me one thing, while I tried to change spacing in Bspline mesh, I tried to print default parameters in simpleElastix, while I did this with console Log on and also tried with write to log file, it both failed to write anything both in console or file.
sitk.PrintParameterMaps(SimpleElastix.GetDefaultParameterMaps("bspline") )
What could be the possible reason for not displaying anything in console? Thanks
@subhajitchatu Are you using IPython? or a standard python console? I know that using Spyder as an IDE, calling this command from IPython console ends printing the output in the standard python console.
@PertuyF
Normal python 2.7 console. Is there any other way to print?
Hi, I am using Python based framework of SimpleElastix. I want to use different non-rigid registration method. I have some issues while using the APIs.
In Transformation enumeration ( GetDefaultParameterMap("bspline")) of namespace Simple, I was able to find only very few methods like bspline,daemons etc. But as in the code base of ITK, there are thin spline, elastic methods, FEM also. How can I use them in transformation methods of python API?
I am also having a problem to display images, like in Matlab there is imshowpair(), which will display overlap between fixed and moving images. Is there any APIs available to achieve the same.
In matlab there are different validation or accuracy estimation methods for rigid and nonrigid registration like Hausdroff distance, SSD etc. Is there any APIs available to achieve the same.
Thanks