New to ITK - moving from ImageJ macros - SimpleITK arithmetic Q

richardM · July 27, 2019, 7:59pm

Hi all,

I’ve started a project to move a bunch of very long complicated ImageJ macros over to using SimpleITK.

I’m new to both Python and ITK so please be nice

I’m working to get some simple registration going for some DICOM files and was using a notebook to play with the code. I wanted to use interact() to change alpha display between two images.

My original files were CT dicom, so were signed int 32’s.

I was calculating a new alpha blended image like this:

blended = (1-alpha)*original + alpha*moved

However this didn’t work, and I figured out it’s because multiplying int32’s by a float like the above results in an array full of zero int32’s.

So I changed things to

blended = (1.0 - alpha)*sitk.Cast(original, sitk.sitkFloat32) + alpha*sitk.Cast(moved, sitk.sitkFloat32)

and I got a new image with alpha applied like I wanted.

So my question is, if you wanted to do arbitrary arithmetic on a sitk image of int32’s, I presume you have to first cast to float, then finally cast back to int32. Is that right?

If so, what is the most efficient (fastest, least memory, functional) way of doing that? I would presume the sitk cast method above is, …is that right?

Thanks,
Richard

blowekamp · July 28, 2019, 1:33pm

Welcome to ITK!

It looks like you had some mark up issues with your post. I corrected it above so that the * operator is there.

The question I have is what is the range and value of alpha?

You can do arithmetic on just about any pixel type in SimpleITK. However, there is no implicit coercion of images types and pixel types. Meaning the types of the images and scalars need to be the same type. So in your case the 1-aplha is convert to int32, which will be 0 in many cases which is why you are getting a black image.

Look into using the ShiftScaleImageFilter to multiply by the float alpha component without having to cast to float.

richardM · July 28, 2019, 11:00pm

Great - thanks for the info Bradley!

So I’ll always remember to cast

I’ve got a follow on question:
I’ve got my initial registration wrapper function working, and it seems to be fine for what I’ll need.
I’m using simple mean squares metric as the similarity, as I’ll always be using images with similar greyscale values/morphology.
So I’m testing things, using a sample image, where I simply translate/rotate the image to test to register back on itself.

So in that case, for Mean Squares, I would have (naively?) thought that the mean square value when the images are registered would be ~0.

However, when I run my code, I get a great fit, in a reasonably fast time, that is robust for the translations/rotations I’m trying - however the Mean Squared final values never approach 0.

I don’t understand why not In my mind they should.

…is that because there is a region where the fixed and moving images don’t overlap, but the registration still try to calculate a mean square value for those pixels, but then gets some non-zero value???

The code I’m using is below. When I run it and view the metric value history it never approaches zero, even for pretty much perfect registrations

import SimpleITK as sitk

def register(fixed, moving, debug=False):
    
    if debug:
        print("fixed:")
        print("  origin:", fixed.GetOrigin())
        print("  direction:", fixed.GetDirection())
        print("  width/height:", fixed.GetWidth(), fixed.GetHeight())
        print("  spacing:", fixed.GetSpacing())
        print("  type:", fixed.GetPixelIDTypeAsString())
        print("  min/max:", sitk.GetArrayViewFromImage(fixed).min(), sitk.GetArrayViewFromImage(fixed).max())
        print("moving:")
        print("  origin:", moving.GetOrigin())
        print("  direction:", moving.GetDirection())
        print("  width/height:", moving.GetWidth(), moving.GetHeight())
        print("  spacing:", moving.GetSpacing())
        print("  type:", moving.GetPixelIDTypeAsString())
        print("  min/max:", sitk.GetArrayViewFromImage(moving).min(), sitk.GetArrayViewFromImage(moving).max())

    if debug:
        metricVals = []
        
        
    # set min values to both images to be zero
    fixedP = fixed - sitk.GetArrayViewFromImage(fixed).min()
    movingP = moving - sitk.GetArrayViewFromImage(moving).min()
    
    # cast to float for registration
    fixedP = sitk.Cast(fixedP, sitk.sitkFloat32)
    movingP = sitk.Cast(movingP, sitk.sitkFloat32)
    
    if debug:
        print("zeroed values")
        print("  fixed min/max:", sitk.GetArrayViewFromImage(fixedP).min(), sitk.GetArrayViewFromImage(fixedP).max())
        print("  moving min/max:", sitk.GetArrayViewFromImage(movingP).min(), sitk.GetArrayViewFromImage(movingP).max())
    
    # set initial transform
    initial_transform = sitk.CenteredTransformInitializer(
        fixedP, 
        movingP, 
        sitk.Euler2DTransform(),
        sitk.CenteredTransformInitializerFilter.GEOMETRY)
    
    
    # setup registration
    registration = sitk.ImageRegistrationMethod()
    registration.SetInterpolator(sitk.sitkLinear)

    # similarity metric to use
    registration.SetMetricAsMeanSquares()
    
    # set optimizer
    registration.SetOptimizerAsRegularStepGradientDescent(
        learningRate = 1.0,
        minStep=1e-3,
        numberOfIterations=100,
        gradientMagnitudeTolerance=1e-6 )
    
    # set transform
    registration.SetInitialTransform(initial_transform)
    
    # setup multi-resolution
    registration.SetShrinkFactorsPerLevel(shrinkFactors = [4,2,1])
    registration.SetSmoothingSigmasPerLevel(smoothingSigmas=[2,1,0])
    registration.SmoothingSigmasAreSpecifiedInPhysicalUnitsOn()   

    if debug:
        registration.AddCommand(sitk.sitkIterationEvent, lambda: metricVals.append(registration.GetMetricValue()))
    
    
    # run registration
    transform = registration.Execute(fixedP, movingP)
    
    
    # finished
    if debug:
        return transform, registration.GetOptimizerStopConditionDescription(), metricVals
    else:
        return transform, registration.GetOptimizerStopConditionDescription()

dzenanz · July 29, 2019, 1:57pm

How close to zero do metric values get? If you try to register similar, but not identical images, how much larger are metric values then? Only perfect registration would result in zero metric, and that is essentially not achievable with reasonably small number of iterations and/or “large” minStep.

richardM · July 29, 2019, 3:50pm

Thanks for the reply Dženan,

When I take a single image and translate/rotate a copy, then run registration, the mean square reported goes from a high at the start of 650,000 to a fitted value of 300,000. You can see this in the images below, showing the metric as well as the fitted image overlaid the “original”.

So these images are identical - except for translation/rotation. So I would expect that a mean squares metric would be essentially zero. (assuming its a mean squared difference that’s being reported).

I tried setting the translation to only (2,2) and rotation to only 2 degrees, and even in that case the metric returned is still ~300,000 when fitted.

These images are 512x512 originally int32’s with a value range ~ -3000-3000. In the example below the 2nd image is a translation by (10,10) and rotation of 25 degrees.

The only reason I can think it isn’t, is that the metric is being calculated without a mask to make sure to only to take values from where the images overlap. But even then, I made sure I moved the minimum of each image to 0 prior to registration, so that if it was doing that, and using a default value of 0 for any undefined points, it wouldn’t affect the result

Unknown Unknown-2

blowekamp · July 29, 2019, 5:51pm

Hello,

One thing missing in your registration code is setting the parameters scale. You can use a method such as ImageRegistrationMethod::SetOptimizerScalesFromPhysicalShift. Setting the scales are important when the transform parameters have different ranges.

If you are uncertain about the resulting MSE, you can subtract the resulting images after the transformation have been applies. Then you can visualize the difference and and computing statistics in it.

richardM · July 29, 2019, 8:23pm

Thanks Bradley - I’ll check it out!

…btw, is there any documentation/info anywhere on the differences b/w all the different options for registration, what’s the pros/cons of setting parameters, and what the differences are? (apart from trying to read all the code docs, which can be rather terse/obscure)

zivy · July 30, 2019, 2:09pm

Hi @richardM,

I think that what you are asking for is somewhere in between the theory of registration (really a broad domain and independent of a toolkit) and the doxygen documentation describing a specific tool (e.g. SimpleITK). I’m not sure we currently have exactly what you want, possibly:

The registration notebooks in our notebook repository.
The ITK software guide’s section on registration which differs slightly from the interface exposed in SimpleITK, but is relevant.

Hopefully one of these is what you are looking for.

richardM · August 9, 2019, 8:22am

Thanks Ziv!

I hadn’t seen the ITK guide - I’ll read that tomorrow for some background