search

ANTsX / ANTsPy

A fast medical imaging analysis library in Python with algorithms for registration, segmentation, and more.
https://antspyx.readthedocs.io
Apache License 2.0
634 stars 162 forks source link

get_center_of_mass() blows up, if image intensity mean is close to zero. #176

Closed MaximilianHoffmann closed 8 months ago

MaximilianHoffmann commented 4 years ago

Describe the bug When I pass images (3D) containing negative values into the antspy registration using the Mattes metric the "joint PDF summed to zero" error is thrown.

To Reproduce For example:

pars_ants = {"type_of_transform": 'Rigid', "aff_metric": 'mattes', "aff_sampling": 32, 
             "aff_random_sampling_rate": 0.01 }
ants.registration(im1,im2,pars_ants)

Expected behavior If I correct the images by upshifting them to positive values it works, and I would have expected that this might be done automatically, if it's a problem, or that a warning should be issued, that negative values are not allowed.

stnava commented 4 years ago

that's not really how we expect the code to be called.

try: 

import ants
fi = ants.image_read( ants.get_data( "r16" ) ) - 128
mi = ants.image_read( ants.get_data( "r64" ) )
reg = ants.registration( fi, mi, 'Rigid', aff_random_sampling_rate=1, random_seed=1, verbose=1 )
MaximilianHoffmann commented 4 years ago 
fi = ants.from_numpy(tmpl.astype('float32'))
mi = ants.from_numpy(vol1_shift_c.astype('float32'))
reg = ants.registration( fi, mi, 'Rigid', aff_random_sampling_rate=0.1, random_seed=1, verbose=1 )

Leads to 

*** Running Euler3DTransform registration ***

Exception caught: 
itk::ExceptionObject (0x55d7cda9e110)
Location: "unknown" 
File: /ANTsPy/itksource/Modules/Registration/Metricsv4/include/itkMattesMutualInformationImageToImageMetricv4.hxx
Line: 312
Description: itk::ERROR: MattesMutualInformationImageToImageMetricv4(0x55d7cda811f0): Joint PDF summed to zero

  Elapsed time (stage 0): 8.5036e-01

Total elapsed time: 8.5187e-01

![Removed]()

Volumes are overlapping (3 projections)

tmpl32=tmpl.astype('float32')
print(tmpl32.shape,tmpl32.max(),tmpl32.min(),tmpl32.dtype)
vol32=vol1_shift_c.astype('float32')
print(vol32.shape,vol32.max(),vol32.min(),vol32.dtype)

gives

(120, 200, 380) 789.5251 -94.60905 float32
(120, 200, 380) 552.9299 -92.44084 float32

The whole output:


b'All_Command_lines_OK
Using single precision for computations.
=============================================================================
The composite transform comprises the following transforms (in order): 
  1. Center of mass alignment using fixed image: 0x55d7cba1e2b0 and moving image: 0x55d7cf52c640 (type = Euler3DTransform)
=============================================================================
[...]

  Shrink factors (level 1 out of 4): [6, 6, 6]
  Shrink factors (level 2 out of 4): [4, 4, 4]
  Shrink factors (level 3 out of 4): [2, 2, 2]
  Shrink factors (level 4 out of 4): [1, 1, 1]
  smoothing sigmas per level: [3, 2, 1, 0]
  regular sampling (percentage = 1.0000e+00)

*** Running Euler3DTransform registration ***

Exception caught: 
itk::ExceptionObject (0x55d7cdb22ce0)
Location: "unknown" 
File: /ANTsPy/itksource/Modules/Registration/Metricsv4/include/itkMattesMutualInformationImageToImageMetricv4.hxx
Line: 312
Description: itk::ERROR: MattesMutualInformationImageToImageMetricv4(0x55d7cdaf4f00): Joint PDF summed to zero

  Elapsed time (stage 0): 8.9281e-01

Total elapsed time: 8.9459e-01
'```
stnava commented 4 years ago

probably need to share the data

brian

On Wed, May 13, 2020 at 8:50 AM Maximilian Hoffmann < notifications@github.com> wrote:

fi = ants.from_numpy(tmpl.astype('float32')) mi = ants.from_numpy(vol1_shift_c.astype('float32')) reg = ants.registration( fi, mi, 'Rigid', aff_random_sampling_rate=0.1, random_seed=1, verbose=1 )

Leads to

Running Euler3DTransform registration

Exception caught: itk::ExceptionObject (0x55d7cda9e110) Location: "unknown" File: /ANTsPy/itksource/Modules/Registration/Metricsv4/include/itkMattesMutualInformationImageToImageMetricv4.hxx Line: 312 Description: itk::ERROR: MattesMutualInformationImageToImageMetricv4(0x55d7cda811f0): Joint PDF summed to zero

Elapsed time (stage 0): 8.5036e-01

Total elapsed time: 8.5187e-01

[image: Figure 51] https://user-images.githubusercontent.com/8191558/81814345-f6332980-9528-11ea-9e13-12acc714b9f3.png

Volumes are overlapping (3 projections)

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/ANTsX/ANTsPy/issues/176#issuecomment-627961232, or unsubscribe https://github.com/notifications/unsubscribe-auth/AACPE7SNRMGVBPUCOJG3D2DRRKJRXANCNFSM4M7QP5MA .

MaximilianHoffmann commented 4 years ago

A cropout is here

stnava commented 4 years ago

can you provide a reproducible example based on this h5 data?

MaximilianHoffmann commented 4 years ago 
import h5py as h5
import ants
import os
import io
import tempfile

path_to_data=pt+'/antspy_test.h5'

os.environ["ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS"] = "1"
original_stdout_fd = sys.__stdout__.fileno()
# Duplicate FD to __stdout__
saved_stdout_fd = os.dup(original_stdout_fd)
# create tempfile and redirect __stdout__
tfile = tempfile.TemporaryFile(mode='w+b')
os.dup2(tfile.fileno(), original_stdout_fd)

f_share=h5.File(path_to_data)
fixed=np.copy(f_share['fixed'])
moving=np.copy(f_share['moving'])

fi = ants.from_numpy(fixed.astype('float32'))
mi = ants.from_numpy(moving.astype('float32'))
reg = ants.registration( fi, mi, 'Rigid', aff_random_sampling_rate=1, random_seed=1, verbose=1 )

tfile.flush()
tfile.seek(0, io.SEEK_SET)
ants_out = tfile.read()
tfile.close()
# reset __stdout__
os.dup2(saved_stdout_fd, original_stdout_fd)
print(str(ants_out).replace('\\n','\n'))

or minimal without console output in notebook

import h5py as h5
import ants
import os

path_to_data=pt+'/antspy_test.h5'

os.environ["ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS"] = "1"

f_share=h5.File(path_to_data)
fixed=np.copy(f_share['fixed'])
moving=np.copy(f_share['moving'])
fi = ants.from_numpy(fixed.astype('float32'))
mi = ants.from_numpy(moving.astype('float32'))
reg = ants.registration( fi, mi, 'Rigid', aff_random_sampling_rate=1, random_seed=1, verbose=1 )
stnava commented 4 years ago

not reproducible - lots of undefined stuff etc

stnava commented 4 years ago

ok - your update helps

stnava commented 4 years ago

here is the issue:

fi2=ants.image_math( fi, "Normalize" )
mi2=ants.image_math( mi, "Normalize" )

>>> ants.get_center_of_mass( fi )
(29994810.950913575, 54067577.80331416, 93929963.00075227)
>>> ants.get_center_of_mass( fi2 )
(34.487496393787595, 49.47746143828112, 39.460844423879415)
reg = ants.registration( fi2, mi2, 'Rigid', aff_random_sampling_rate=1, random_seed=1, verbose=1 )
# also works
reg = ants.registration( fi-1, mi-1, 'Rigid', aff_random_sampling_rate=1, random_seed=1, verbose=1 )

center of mass calculation - used to initialize the registration - gives nonsense ..... this has something to do with your conversion from numpy and/or your read from the h5 file. not sure worth debugging as likely to be rare / dataset specific. if you can find out why, we would appreciate the assistance/enlightenment. or if you demonstrate that it's not specific to this I/O setup.

MaximilianHoffmann commented 4 years ago

Thanks! at least a point to start.

MaximilianHoffmann commented 4 years ago

But allow me one question: after conversion to the antspy image the date is fixed in float 32 format. How can any center of mass calculation give something that is outside of the image coordinates given arbitrary inputs?

stnava commented 4 years ago

it looks to be a type issue or maybe pointer issue.

that's why the fi -1 mi - 1 conversion has an impact ---- the internal conversions we do to perform math on the images seems to "fix" whatever issue is going on .... it was just a guess that led me to discover this.

MaximilianHoffmann commented 4 years ago 
import ants
import numpy as np
np.random.seed(21)
ants.from_numpy(np.random.randn(100,100,100).astype('float32')).get_center_of_mass()

gives

(12086.293415810545, -107.79403426152939, 2221.4008447057618)

Numpy 1.18.1 but you can also execute it a couple of times with a random seed and it will yield something similar.

stnava commented 4 years ago

try 

import ants
import numpy as np
np.random.seed(21)
ants.from_numpy(np.random.randn(100,100,100).astype('float64')).get_center_of_mass()
MaximilianHoffmann commented 4 years ago

Identical. Also

import ants
import numpy as np
np.random.seed(21)
ants.from_numpy(np.random.randn(100,100,100)).get_center_of_mass()
stnava commented 4 years ago

hmm - I see repeated calls go from nonsense ( 1st call ) to something expected 2nd call and so on ...

import ants
import numpy as np
np.random.seed(21)
ants.from_numpy(np.random.randn(100,100,100)).get_center_of_mass()
ants.from_numpy(np.random.randn(100,100,100)).get_center_of_mass()
ants.from_numpy(np.random.randn(100,100,100)).get_center_of_mass()

gives

(12086.293415810545, -107.79403426152939, 2221.4008447057618)
(65.76023786090398, 45.39262471805497, 66.97123286337336)
(58.08016628651665, 43.97373502386715, 57.241743248005584)
(24.93346743221081, 105.05868214978264, 96.25459333865997)

first one is wrong - the others are correct .....

MaximilianHoffmann commented 4 years ago

I think it's (pseudo)random, if i try 100 runs...

(12086.293415810545, -107.79403426152939, 2221.4008447057618)
(65.76023786090398, 45.39262471805497, 66.97123286337336)
(58.08016628651665, 43.97373502386715, 57.241743248005584)
(24.93346743221081, 105.05868214978264, 96.25459333865997)
(37.51878263120262, 37.97909291308354, 35.07218977822534)
(-26.43213420101111, -78.66991920410356, 118.21344568918884)
(-7.5657297061659365, 28.68783884145177, 83.81479392087338)
(9.114983768329353, 22.557688173311405, 57.84102932654726)
(76.41799911396069, 71.84808402413948, 60.508508665312505)
(78.36049935083078, 100.56934582371844, 53.125090126883)
(57.510135623927354, 53.21558306281195, 62.23516082889067)
(136.1119847899379, 332.0815631902609, -307.1114304849334)
(15.709290226620517, 45.33531006434985, 89.29800818652997)
(249.21274780792217, 89.38111331967765, 183.01316433456785)
...
(545.7921723066764, -277.6391008412984, -1637.0821599519966)
...
(553.2534689040832, 172.56216743122363, -218.0160528193693)
...
(-161.29343213100105, 5.7332478989255335, 198.38503488514007)
...
(1493.0586377559048, 1248.94025028152, -4492.808397278152)
...
(922.9460252268083, 719.2123546863406, -409.57641697390443)
...
stnava commented 4 years ago

those all look correct to me (except the first one). the center of mass is just the intensity multiplied by the coordinate system and then averaged ....

also note:

import ants
import numpy as np
np.random.seed(21)
x=ants.from_numpy(np.random.randn(100,100,100))
( x - 1e-2).get_center_of_mass()

appears to "fix" things again. just trying to develop some idea here about what is going on ...

MaximilianHoffmann commented 4 years ago

Ok, maybe I don't get the standard coordinate system, what would be the maximum value for 100x100x100 pixel?

import ants
import numpy as np
np.random.seed(21)
for ii in range(10000):
    ct=ants.from_numpy(np.random.randn(100,100,100).astype('float32')).get_center_of_mass()
    if ct[0]>1000:
        print(ct)
(12086.293415810545, -107.79403426152939, 2221.4008447057618)
(1493.0586377559048, 1248.94025028152, -4492.808397278152)
(3022.2361681970824, 57.52413359525336, -3747.033571040534)
(1197.9346702658174, 1081.3120974717056, -434.2635009868153)
(4928.967510084977, 3708.299860828199, 1273.2053591009965)
(15061.958531328044, -3089.8641795535177, 2506.3287609435283)
(4745.902291833249, -1139.5378107221363, 3874.5457226527647)
(1120.3095204663543, 303.97339436128357, 67.37104296708097)
(5923.067654958986, 2116.759613789938, -9006.318428378694)
(1308.4338114947448, -146.17851725843462, -514.9508493131104)
(3924.276457587353, 8327.882051282497, -634.5892166825039)
(1998.9225874503834, 2366.8752179412554, 982.4465048077085)
(1023.4391452055147, -539.618788661797, -468.05093883555116)
(2244.1273565383617, -981.7173734849376, 3246.020934147846)
(9790.857828605005, 1081.05869244328, -471.0502058501832)
(1555.7318569556517, -343.90314222323417, 291.9174101017077)
(1426.3579268964638, 1397.2584926578418, 1049.595807627448)
(1357.3845677541715, -5195.014709839194, 12743.92136976191)
(5005.834420744816, 2914.6394615228296, 3858.6426469265216)
(6550.641978373897, 11210.882152036234, -8017.70162527153)
(174308.95881019876, 28834.33186581375, -104604.03019130287)
(2503.1509991212524, 2179.934210821127, 2465.701280461111)
...
stnava commented 4 years ago

if the input intensity is bound in zero to one and everything is set to identity then 100

MaximilianHoffmann commented 4 years ago

So it's not necessarily a real center of mass, in the sense of three coordinates, if the image intensity is not bounded by zero and one?

stnava commented 4 years ago

see https://itk.org/Doxygen/html/classitk_1_1ImageMomentsCalculator.html for details

MaximilianHoffmann commented 4 years ago

center of gravity ? is it like M_01/M_00 . Do you happen to know, whether there is a mathematical description of what the moment calculator is implementing, it expect also negative pixels, yes? Because otherwise M_00 could get very close to zero.

stnava commented 4 years ago

I ran the same experiment in ANTsR and the behavior is similar .... 

getCenterOfMass( antsImageRead( '/tmp/temp.nii.gz' )
[1] 12086.293  -107.794  2221.401
> getCenterOfMass(antsImageRead( '/tmp/temp.nii.gz' ) + 0.01)
[1] 53.42208 49.44875 50.20769

suggests that the issue is inside the underlying C++ to which I linked.

MaximilianHoffmann commented 4 years ago

I think, that's what's going on.

import ants
import numpy as np

np.random.seed(21)
x=ants.from_numpy(np.random.randn(100,100,100))
x=x+x.min()
print(( x - 1e-5).get_center_of_mass())
print(( x - 1e-7).get_center_of_mass())
print(( x - 1e-4).get_center_of_mass())

Gives

(49.49165402783811, 49.50010906338797, 49.49849406443995)
(49.491654010072445, 49.500109063675566, 49.49849406135992)
(49.491654187952314, 49.500109061939696, 49.498494093187055)

Which is what you would expect for a noisy constant density. So the issue seems to carry a proper title after all

stnava commented 4 years ago

It has nothing to do with the metric - it's the initialization.

MaximilianHoffmann commented 4 years ago

Alright...I meant the non-negativity part... Feel free to change the title.

stnava commented 4 years ago

@ntustison let me know what you think about this ... short story is the itk moments calculator easily produces nonsense if the intensity is a zero centered gaussian .... did we encounter / discuss this before? seems familiar

stnava commented 4 years ago

@MaximilianHoffmann yes - that is an issue if intensity is close to zero mean

MaximilianHoffmann commented 4 years ago

EDITED AFTER CODING MISTAKE

import ants
import numpy as np
xx=np.arange(-100,100,0.1)
X,Y=np.meshgrid(xx,xx)
im=ants.from_numpy(Y)
im.get_center_of_mass()
(im).get_center_of_mass()
(-665667.0000128918, 999.4999999999973)

It would be good if a warning is issued, no?

stnava commented 4 years ago 
(im*0+1).get_center_of_mass()

is the com of the coordinate system

ntustison commented 4 years ago

I can't say I remember a discussion but there's a check for exactly non-zero mass due to normalization so I'm guessing that's where you're running into problems.

stnava commented 4 years ago

Will add TODO: document behavior of get_center_of_mass in ANTsR, ANTsPy and possibly add a warning based on the mean of the image intensity.