@PertuyF Looks good from here except I think vector.PushBack() should be vector.push_back() :)

from simpleelastix.

Hi sunhanjie. I guess you are registering 2D images? This problem occurs when a 2D image is passed to elastix, but elastix interpret the image as 3D and tries to smooth the image along the third dimension (dimension 2 mentioned in the log) which has a length of 1. You have to options.

1. Use the SimpleITK's ExtractImageFilter() to extract the 2D slice from your 3D image and pass that to elastix, or
2. use "FixedShrinkingImagePyramid" and "MovingShrinkingImagePyramid" when can handle single-slice 3D images.

Let me know how it goes!

from simpleelastix.

Hi kaspermarstal. Thanks for your response.
Yes, you are right, I test to register 2D images, using two slices of DICOM image as fixed and moving image, as follow,

fixedImage = sitk.ReadImage('/home/s/Data_img/DICOM1/Recon0210.dcm')
movingImage = sitk.ReadImage('/home/s/Data_img/DICOM2/Recon0210.dcm')

I am very sorry, this is my first time to use it, and I just follow the introduction to test. I am wondering what kind of format of data can be read by sitk.ReadImage(). And I still have no idea how to solve it with your answer. Could you please explain more details and post some codes of demo.
Thanks a lot.

from simpleelastix.

sitk.ReadImage() can read most file formats you need and indeed DICOM images. That part of your code is fine. ReadImage() interprets the dicom slice as 3d because probably it's part of a 3d series. Right? In this case the header will contain 3d information. To use the ExtractImagerFilter, you can use the following code (python)

fixedImage3d = sitk.ReadImage('/home/s/Data_img/DICOM1/Recon0210.dcm')
fixedImage2d = sitk.Extract(fixedImage3d, (fixedImage3d.GetWidth(), fixedImage3d.GetHeight(), 0), (0, 0, 0))

The first tuple gives the size of the region you want to extract (width, height). The zero denotes that the 3rd dimension should be discarded. The second tuple (0,0,0) gives the starting index of the bounding box you wish to extract. In this case these settings are easy because you want to extract the whole image. Try this first. This will also tell elastix to do 2d registration.

If this doesnt work, you can try to to use shrinking pyramid. Simply change these settings in the parameter map before you start elastix.

parameterMap = sitk.GetDefaultParameterMap("affine")
parameterMap["FixedImagePyramid"] = ["FixedShrinkingImagePyramid"] 
parameterMap["MovingImagePyramid"] = ["MovingShrinkingImagePyramid"] 

SimpleElastix = sitk.SimpleElastix()
SimpleElastix.FixedImage(fixedImage)
SimpleElastix.MovingImage(movingImage)
SimpleElastix.SetParameterMap(parameterMap)
SimpleElastix.Execute()
sitk.WriteImage(SimpleElastix.GetResultImage(), "resultImage.dcm")

While the image readers and writers can work with dicom data, it is worth mentioning that any meta info in the original dicom header is lost. Only spatial information (origin, orientation, pixel spacing) is saved in the result image.

from simpleelastix.

Thanks a lot for your considerate explain and forward. I get your points, but it doesn't work either. The code of first solution is followed, is it right?

fixedImage3d = sitk.ReadImage('/home/s/Data_img/DICOM/Recon0210.dcm')
fixedImage2d = sitk.Extract(fixedImage3d, (fixedImage3d.GetWidth(), fixedImage3d.GetHeight(), 0), (0, 0, 0))
movingImage3d = sitk.ReadImage('/home/s/Data_img/DICOM/Recon0215.dcm')
movingImage2d = sitk.Extract(movingImage3d, (movingImage3d.GetWidth(), movingImage3d.GetHeight(), 0), (0, 0, 0))

parameterMap = sitk.GetDefaultParameterMap('translation')

elastix = sitk.SimpleElastix()
elastix.SetFixedImage(fixedImage2d)
elastix.SetMovingImage(movingImage2d)
elastix.SetParameterMap(parameterMap)

elastix.Execute()

However, both of two ways lead to the same error. It seems I need to follow as the post two ways, but how?

Running elastix with parameter map 0
Current time: Tue Feb 23 11:01:01 2016.
ERROR:
FixedImageType: 2D short
MovingImageType: 2D short
elastix was not compiled with this combination of ImageTypes!
You have two options to solve this:

1. Add the combination to the CMake parameters ELASTIX_IMAGE_nD_PIXELTYPES and ELASTIX_IMAGE_DIMENSIONS, re-cmake and re-compile.
2. Change the parameters FixedInternalImagePixelType and/or MovingInternalImagePixelType in the elastix parameter file.

ERROR:
Something went wrong in the ComponentDatabase
Errors occurred!

from simpleelastix.

Elastix is only compiled for float in 2d by default so you will have to cast your image to this pixel type. You can cast your image to float like this

fixedimage2dfloat = sitk.Cast(fixedimage2d, sitk.sitkFloat32)

This is inded not very user friendly and will be changed in the next release :)

from simpleelastix.

Thanks a lot. It finally can work.
And so does saving Resultimage after registration, I have to change pixel type to int in order to sitk.WriteImage(), is it?

from simpleelastix.

Something like this?

resultFloatImage = SimpleElastix.GetResultImage()
resultIntImage = sitk.Cast(resultFloatImage, sitk.sitkInt32) # Or 8, 16, 64, whatever you like ;)
sitk.WriteImage(resultIntImage, "resultImage.dcm")

from simpleelastix.

Yes, I have already handled it ,^_^ ~ Thanks all the same.
I am very grateful for your time and help.

Best Regards.

在 2016年2月23日，下午12:03，Jorrit Posthuma [email protected] 写道：

Something like this?

resultFloatImage = SimpleElastix.GetResultImage()
resultIntImage = sitk.Cast(resultFloatImage, sitk.sitkInt32) # Or 8, 16, 64, whatever you like ;)
sitk.WriteImage(resultIntImage, "resultImage.dcm")
—
Reply to this email directly or view it on GitHub #10 (comment).

from simpleelastix.

@kaspermarstal Sorry for reopening the thread. I was trying to save 3d image in .dcm using
sitk.WriteImage(resultImage, "resultImage.dcm"). It gave me an error that GDCM doesn't support that type.
But as I changed it to
finalImage = sitk.Cast(resultImage, sitk.sitkUInt32)

It successfully saved into DICOM image. But as we move from Float to int ot Unit. Informations are lost and the image looks different. Is there any other way to save DICOM images?
Again when I moved to .mha, it worked. But ImageJ failed to open that. Anyway my preference is to save output image in .dcm without information loss.

Thanks

from simpleelastix.

You should use the ImageSeriesWriter() for writing dicom series :)

from simpleelastix.

@subhajitchatu Beware of ImageJ as by default it will adjust the dynamic range displayed for each image. I would advise to set the same range if you want to compare two images by eye (Open your two images then open the B&C settings with Image>Adjust>Brightness/Contrast... or Ctrl + C, then Set and enter min and max value, check Propagate to all other open images then OK).

However in your case I suspect that casting your image back to sitkUInt32 may play a role, as CT image usually makes use of signed integers you may want to cast to sitkInt32 instead.

Note that I personally disliked using DICOM format for images manipulations. In my experience transfer rates is often lower than for other formats, and you duplicate a lot of useless metadata each time you write a modified image.
I usually keep the raw image as DICOM somewhere (then I can get metadata if needed), then I create a copy in a different format (choose one that is simpler to handle) and only use this format for all modified images.

from simpleelastix.

@kaspermarstal

Thanks a lot for ur suggestion. The issue or rather overhear I found that, while the registeredImage is a 3D image with 60 slice, while I use ImageSeriesWriter(), I have to specify 60 different names, thats also fine. But to make it a single dicom image()combing 60 slices) as registeredImage is again not possible. :P
Any suggestion?

@PertuyF
The Image is a MRI image. sitkInt32 is not worked so I have used sitkUInt32. But it is visually clear that the image is losing some information, it looks like some morphological operations are done on the image.
I saved in .mha format, do u know any good software to display that type, ImageJ is not able to display .mha in Linux.

Another doubt is about reading a 4D image, like DCE-MRI
It has total 4800, (1.dcm to 4800.dcm), these are total 100 stack, each has 48 slices (100*48=4800). One way is to use ImageSeriesReader() in a loop for 48 slices and save 100 stack in a 4D array , but if image Registration or smoothing functions are not supporting 4D image, then issue is still there.

Is there any method to handle(register 99 stack with 1 stack) this type of 4D Image stack?(In SimpleElastix or SimpleITK)

from simpleelastix.

@subhajitchatu

I would say the loss of information you are observing sounds surprising to me, but I am not sure how ITK Cast() deals with float64 to int32 conversion, so it may be responsible for it.
Are you sure it is not linked to the registration itself ? You could try casting the image to float, then back to int just see if you observe the same kind of thing. Or just try saving the result image in a format that support float64, and see if these alterations are present.

Open .mha images

There is probably an ImageJ plugin to open .mha images, but I don't know which one. Anyway, it seems that Fiji deals with this format out of the box, so you may want to use this distribution instead of the "vanilla" ImageJ.

Open 4D .dcm series

When working with DICOM images I found PyDicom module easier than SimpleITK's ImageSeriesReader(). And if you are interested there is a solution here to deal with 4D images.

Register 4D images

I do not really get what are your 100 stacks, time series? different contrasts?
What would be your fixed and moving images in this case?

If you want to register all stacks from the series together you could try the groupwise registration.

If you want to register stack0 to a reference image, then transform all stacks according to registered stack0, IF the stacks are aligned to one another, then you can:

• perform registration of stack0 to your reference
• get the transformParameterMap
• loop over the other stacks and apply the transform using TransformiX

from simpleelastix.

@PertuyF

Thanks for your reply!

Image Type Cast :

I don't think registration has an impact on the result Image. Because when I use matplotlib just to display 2D slices, it shows perfectly. The problem arises when I try to write the result image in .dcm
resultImage = SimpleElastix.GetResultImage();

It fails with a message that GDCM doesn't support this type.
Now I tried with different type like int16, int32 but it was not working. Then I changed it to UInt32 and it successfully executed. The image intensity range of the original image was 0 - 1000 , so I thought there shouldn't be any problem. But as I found the registered Image is not like the original one, it seemed like some pixels were lost.
Anyway I will try your suggestions and will let you know.

Open .mha images

Thanks! I will try to install Fiji. Is there any preferable image type extension which is popularly used in medical imaging?

Open 4D .dcm series

I have seen PyDicom before, but I am not sure after reading 4D image using PyDicom how can I pass it to SimpleITK or SimpleElastix framework. Because to use SimpleITK functionality we need to change it to SimpleITK image object (may be via numpy).
Again the same issue will happen while writing
resultImage = SimpleElastix.GetResultImage();
ds.dicom_write(resultImage) ??
There is no clear documentation regarding this.

Register 4D images

DCE-MRI is a 4D time series image. With different contrast.
Stacks means
[stack0 , stack1 , ..... , stack100]
and each stack has 48 slices
[{1,2,3,....,48};{1,2,3,...,48}; ...... {1,2,3,...,48}]
So, slice 1 will repeat after 48 th slice but with different contrast. These are subsequent stacks after injecting contrast agent.
Lets say we have a 3D fixed Image. First we can register stack0 with the fixed image and then
to register other 99 stack, what we need is to register stack1 , stack 2 , ... , stack 100 with Stack0 individually.
Once all the stacks are aligned w.r.t Stack0. It will automatically aligned with the fixed Image.

What you suggest is absolutely correct in this scenario. But can you elaborate what will be the difference if we do GroupWise Registration(first register stack0 with reference and then groupwise registration among the other stacks) and the method you suggested later(register with reference and get parameter map and loop through the other stacks).

from simpleelastix.

@subhajitchatu

Image file formats

I would say that whereas DICOM as spread well and is now the standard format for raw data that output from medical imaging devices, there is no such consensus for any other file format.
Some communities or institutes have set there how formats, but the perfect one does not exists yet.

• OME-TIFF is an improvement over the classical TIFF format, compatible with the former and using the extension .ome.tif or .ome.tiff, but readable using .tif or .tiff too. It is spreading with the use of the BioFormat library, strongly used in ImageJ/Fiji. However, as it relies on TIFF, it may cause problems with 4D images.
• NIfTI I used a lot, it is derived from the Analyze 7.5 format that used the pair .hdr + .img files, but can also be used with one only file with .nii extension. It seems this on is used in the neuro community. Note that is uses the RAS orientation by default, which can lead several software to flip the image at opening (those following LAS orientation by default).
• NRRD is pretty handy to use too, either a .nrrd file or a pair of .nrhd + .raw files.
• MINC is used by several big centers too, with .mnc files.
• And there is plenty others...

PyDicom to SimpleITK

PyDicom uses Numpy arrays natively for image data, so once you extract the data you can convert it very simply to sitk.Image()

import dicom

ds = dicom.read_file(yourFile.dcm)  # Read file
array = ds.pixel_array  # Extract data as numpy array
# Reshape the 4D dcm data here if needed
image = sitk.GetImageFromArray(array)  # Convert to sitk.Image

Register 4D images

Then Groupwise registration won't fit in your case. It was meant for registration of 3D or 4D time series (3D+T or 2D+T) but can also be used to create population average if you create a file that is 3D + Samples for instance. It will not allow you to register your data to an external reference.

Then you should use the loop strategy

import SimpleITK as sitk
selx = sitk.SimpleElastix
tfmx = sitk.SimpleTransformix

stacks  # Assume this is you 4D data as sitk.Image in XYZT order
moving = stacks[:, :, :, 0]  # This should extract stack0 but unfortunately I cannot test it, so try it before
vector = sitk.VectorOfImage()  # This will be used to store your result images

selx.SetFixedImage(fixed)
selx.SetMovingImage(moving)
...  # Set your parameters and everything needed
selx.Execute()

transfromParameters = selx.GetTransformParameterMap()  # Get the transform parameters to apply transform to other stacks

vector.push_back(selx.GetResultImage())  # Store registered stack0

# Now iterate through your 4D stack:
for i in range(1, stacks.GetSize()[0] + 1):
    moving = stacks[:, :, :, i]
    tfmx.SetMovingImage(moving)
    tfmx.SetTransformParameterMap(transformParameters)
    tfmx.Execute()

    vector.push_back(tfmx.GetResultImage())

registeredImage = sitk.JoinSeries(vector)  # Get the result as an image
# Save it or do anything

Note that I cannot test this, so there is likely typos, or maybe I got something wrong (I usually use numpy so my SimpleITK is a little rusted), but at least you have an idea of the process.
@kaspermarstal please correct if I said something wrong.

EDIT: PushBack() was corrected to push_back() following Kasper's comment (Thank you Kasper :) )

from simpleelastix.

@PertuyF

Thanks for this great post!

Read Dicom from a directory

PyDicom provides only read_file() API , which is only for a single file while DICOM images are normally stored in a directory which contains Image.dcm [eg. Image001.dcm, Image002.dcm] , so to make it read to a single stack as you suggested, one needs to have knowledge about the which file to read and what to read! If we use

    `data_directory` = os.path.dirname(inputPath)

    `# Get the list of files belonging to a specific series ID.
    `reader = sitk.ImageSeriesReader()
    `series_ID = reader.GetGDCMSeriesIDs(data_directory)[0];
    `# Use the functional interface to read the image series.
    `original_image = sitk.ReadImage(reader.GetGDCMSeriesFileNames(data_directory, series_ID))

It will automatically reads all Images from the Input Directory, Is there other way around in PyDicom.

Write PyDicom:

write_file() on numpy array normally doesn't work. As http://stackoverflow.com/questions/14350675/create-pydicom-file-from-numpy-array suggested it is really difficult to write our output to a file.

from simpleelastix.

As you pointed out, with single page DICOM folders, you have to recursively read the files anyway.
Writing a DICOM file is not hard in itself, but creating one that is compliant with the standard and will be handled properly by all softwares... ...that's why I don't keep my files in DICOM format.

The tricky part is getting the metadata right, with some unique identifiers, so creating dicom from scratch is a mess. It is easier to edit an image and rewrite it using the original metadata. If you use single-page DICOM you will need to recursively write the files with editing the metadata correctly for each. It could be easier to use multi-page DICOM, but I already tried and whereas pydicom was able to read the file back without any problem, other software weren't.

Anyway, this is off topic regarding SimpleElastiX, so you should probably talk directly on the github project of pydicom if you are interested.

Note that it seems that PyDicom is right in the middle of a transition to 1.0 version, which is a major release and comes with refactored module and methods too.

from simpleelastix.

Hi @manningzhang ,

Your message is barely comprehensible. Could you at least post the accurate error message ? without that it is not possible to help.

from simpleelastix.

@PertuyF I'd like to know that If I can read a series DICOM files and convert them into a 3D image. Just like the data type "mhd”

from simpleelastix.

@manningzhang That is actually a SimpleITK question, while this is the SimpleElastix repo. But nevertheless, this should get you going:

path = '/some/local/dicom/folder'
reader = sitk.ImageSeriesReader()
serie_ids = reader.GetGDCMSeriesIDs(path)
serie_images = []

for serie_id in serie_ids:
    dicom_filepaths = reader.GetGDCMSeriesFileNames(path, serie_id)
    reader.SetFileNames(dicom_filepaths)
    image = reader.Execute()
    serie_images.append(image)

(This is actually for reading multiple series in one folder, but it shows what you can do.)

from simpleelastix.

@JorritPosthuma Thanks a lot !!!

from simpleelastix.

@kaspermarstal Hello , I am doing 2D to 2D image registration with C++, and I meet some similiar issues discripted from @sunhanjie, the log file told me to" add the combination to the CMake parameters ELASTIX IMAGE_nD_PIXELTYPES".

I see your solution using 'extract' function on the above , but the using method seems different with python, can you teach me how to use extract function in C++? thanks so much~

fixedImage3d = sitk.ReadImage('/home/s/Data_img/DICOM1/Recon0210.dcm') fixedImage2d = sitk.Extract(fixedImage3d, (fixedImage3d.GetWidth(), fixedImage3d.GetHeight(), 0), (0, 0, 0))

from simpleelastix.

You can see the C++ signature of the function here:

Image itk::simple::Extract ( const Image & image1, 
                             const std::vector< unsigned int > &, size = std::vector< unsigned int >(4, 1), 
                             const std::vector< int > & | index = std::vector< int >(4, 0), 
                             ExtractImageFilter::DirectionCollapseToStrategyType directionCollapseToStrategy = itk::simple::ExtractImageFilter::DIRECTIONCOLLAPSETOGUESS)

The usage in C++ is identical to that of Python (the Python version is the C++ version wrapped by SWIG after all), except that you pass std::vectors instead of tuples. Happy coding!

from simpleelastix.

Hi @kaspermarstal
Right now , what I am facing is that both of my original input data are 2D images,
once I load affine and b-spline txt file inside the elastix process, it will break down
my code is as below

sitk::Image CT = sitk::ReadImage("CT.tif");
sitk::Image CTT = sitk::Cast(CT, sitk::sitkFloat32);
sitk::Image MR = sitk::ReadImage("MR.tif");
sitk::Image MRR = sitk::Cast(MR, sitk::sitkFloat32);
sitk::SimpleElastix elastixRigid;
ElastixProcess(CTT, MRR, elastixRigid, "Affine.txt");

and my sub function of ElastixProcess is as below
void ElastixProcess(sitk::Image fixed, sitk::Image moving, sitk::SimpleElastix &elastix, std::string parameter){
elastix.SetFixedImage(fixed);
elastix.SetMovingImage(moving);
cout << "1"<<endl;
elastix.SetParameterMap(sitk::ReadParameterFile(parameter));
cout << "2" << endl;
elastix.LogToFileOn();
elastix.LogToConsoleOn();
cout << "3" << endl;
elastix.Execute();
cout << "4" << endl;

}
the code will break in .Execute and the log error will be :

Description: itk::ERROR: AdaptiveStochasticGradientDescent(000000000B24F090): ERROR: Scales must have value greater than epsilon! Scale[2] = 0

and my affine parameter is as below which is setting for 3d registration and in use for 2d images :

(FixedInternalImagePixelType "float")
(MovingInternalImagePixelType "float")
(UseDirectionCosines "true")
(Registration "MultiMetricMultiResolutionRegistration")
(Interpolator "BSplineInterpolator")
(ResampleInterpolator "FinalBSplineInterpolator")
(Resampler "DefaultResampler")
(FixedImagePyramid "FixedSmoothingImagePyramid")
(MovingImagePyramid "MovingSmoothingImagePyramid")
(Optimizer "AdaptiveStochasticGradientDescent")
(Transform "AffineTransform")
(Metric "AdvancedMattesMutualInformation" "TransformBendingEnergyPenalty")
(Metric0Weight 10)
(Metric1Weight 200)
(AutomaticScalesEstimation "true")
(AutomaticTransformInitialization "true")
(HowToCombineTransforms "Compose")
(NumberOfHistogramBins 60)
(ErodeMask "false")
(NumberOfResolutions 4)
(ImagePyramidSchedule 40 40 8 20 20 4 10 10 2 5 5 1 )
(MaximumNumberOfIterations 1500)
(NumberOfSpatialSamples 5000)
(NewSamplesEveryIteration "true")
(ImageSampler "MultiInputRandomCoordinate")
(SampleRegionSize 150)
(UseRandomSampleRegion "true")
(MaximumNumberOfSamplingAttempts 5)
(RequiredRatioOfValidSamples 0.1)
(BSplineInterpolationOrder 1)
(FinalBSplineInterpolationOrder 3)
(ShowExactMetricValue "false")
(WriteTransformParametersEachResolution "true")
(DefaultPixelValue 0)
(WriteResultImage "true")
(ResultImagePixelType "short")
(ResultImageFormat "nii.gz")

Sorry for my description is so long , it my first quize, please let me know if any comment, thanks

from simpleelastix.

Elastix thinks your images are 3D, hence the message Scale[2]=0i.e. the third entry in Scale is zero. This probably because your sitk::Image is 3D with one slice. What does e. g.CT.GetSize() return? If it is something like (256,256,1) (depending on how big your images are of course) you can use sitk::Extract to extract the plane so that GetSize() will return (256,256) i.e. remove the trailing singleton dimension.

from simpleelastix.

Hi @kasper
after I change my code as below , my elastix process start to work
sitk::Image CT = sitk::ReadImage("CT.tif", sitk::sitkFloat32);
sitk::Image MR = sitk::ReadImage("MR.tif", sitk::sitkFloat32);
but When I use the original one as below , it will fail
sitk::Image CT = sitk::ReadImage("CT.tif", sitk::sitkUInt8);
sitk::Image CTT = sitk::Cast(CT, sitk::sitkFloat32);
Do you have any idea about the difference between 2 method ? thanks so much for the kindly reply

btw do you know how to \n in the code description `` ?

from simpleelastix.

What is the error message with which it fails?

from simpleelastix.

Hi kaspermarstal,
I am new to elastix, sorry if you had answered this question before, what's the best way to register RGB images? When I cast my RGB images to float, the result image after registration is also with pixel type float, is there a way to get the result image as RGB?
Thanks

from simpleelastix.

You can split the images into three grey-level images, one for each color channel, and then set three fixed images, and three moving images. Just make sure that the color channels correspond to each other, i.e. the first fixed images is red, the first moving image is red, etc etc.

from simpleelastix.

Many thanks kaspermarstal, I will try that.

from simpleelastix.

Hello @kaspermarstal~
I am curious that during the registration process of Elastix , whether the system execute the process of feature extraction like Sift or Surf, or it just aims whole fixed and moving images to calculate the transformation formula.
and one more question is that is it possible to let elastix to registrate based on some additional landmark? thanks so much ~

from simpleelastix.

@peterzzz12 It randomly selects a user-specified number of positions to compute the similarity (typically ~2000). It is possible to do registration based on other features.

You can base registration also on Landmarks and/or pack custom features into an image and use MultiResolutionRegistrationWithFeatures.

from simpleelastix.

@kaspermarstal Do you mean user set numbers of landmark in the ParameterMap, and the system randomly allocate these landmarks on the fixed and moving images to register? If so, may I know which parameter it is ?

In addition, I also want to ask that, right now i am using sitk library, but I dont see any function about MultiResolutionRegistrationWithFeatures, do you have any idea that how can I use this function with my setting landmark?
Thanks for your kindly reply :)

from simpleelastix.