ANTsX
get_center_of_mass() blows up, if image intensity mean is close to zero.
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.
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 )
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
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
'```
probably need to share the data
brian
On Wed, May 13, 2020 at 8:50 AM Maximilian Hoffmann < [email protected]> 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 .
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 )
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.
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?
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.
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.
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()
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()
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 .....
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)
...
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 ...
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)
...
if the input intensity is bound in zero to one and everything is set to identity then 100
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?
see https://itk.org/Doxygen/html/classitk_1_1ImageMomentsCalculator.html for details
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.
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.
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
Alright...I meant the non-negativity part... Feel free to change the title.
@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
@MaximilianHoffmann yes - that is an issue if intensity is close to zero mean
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?