[MINC-users] mincresample modifying data

Wed Aug 14 20:20:05 EDT 2013

Hi Andrew,

On 15 August 2013 05:18, Andrew Wood <andrew at biospective.com> wrote:
> I've been experiencing some behaviour coming out of mincresample that I
> can't explain. I've been trolling through the source code and wondering
> about a number of things that might contribute to my problem, but I've
> boiled it down to one simple example:
>
> Here is a volume in.mnc (a signed float volume):
> https://www.dropbox.com/s/cyfujyenc4j7h8b/in.mnc
>
> $ mincresample in.mnc rsl.mnc
>
> Then compare with mincstats:
>
> $ mincstats -min in.mnc
> Min:               0
> $ mincstats -min rsl.mnc
> Min:               -5.345658291e-11

My guess is that there is something afoot (NaN perhaps?) in zslice 12.

   $ mincheader rsl.mnc
   ...
   data:

   image-min = 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -5.34565829102274e-11, 0, 0,
       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,

If you then do this:

   mincreshape -dimrange zspace=0,10 in.mnc z.mnc
   mincresample z.mnc rslz.mnc
   mincdiff z.mnc rslz.mnc

There are some rounding errors (e-36) but the image data itself is
identical.

A bit of digging reveals this:

   mincreshape -dimrange zspace=12,0 in.mnc z.mnc
   mincresample z.mnc rslz.mnc -clobber

   mincdump z.mnc > z.dat
   mincdump rslz.mnc > rslz.dat

   diff z.dat rslz.dat

1c1
< hdf5 z {
---
> hdf5 rslz {
12c12
< 		image:valid_range = 0.f, 26989.4f ;
---
> 		image:valid_range = -5.345658e-11f, 9456.072f ;
24c24
< 	double xspace ;
---
> 	int xspace ;
31a32,33
> 		xspace:step = 0.06 ;
> 		xspace:start = -6.33 ;
35,37c37
< 		xspace:step = 0.06 ;
< 		xspace:start = -6.33 ;
< 	double yspace ;
---
> 	int yspace ;
44a45,46
> 		yspace:step = 0.06 ;
> 		yspace:start = -9.663 ;
48,50c50
< 		yspace:step = 0.06 ;
< 		yspace:start = -9.663 ;
< 	double zspace ;
---
> 	int zspace ;
57a58,59
> 		zspace:step = 0.06 ;
> 		zspace:start = -4.08 ;
61,62d62
< 		zspace:step = 0.06 ;
< 		zspace:start = -4.08 ;
68c68
< 		:ident = "rotor:cai-harold.cai.uq.edu.au:2013.08.15.10.06.37:22297:1" ;
---
> 		:ident = "rotor:cai-harold.cai.uq.edu.au:2013.08.15.10.07.58:23195:1" ;
70a71
> 			"Thu Aug 15 10:07:58 2013>>> mincresample z.mnc rslz.mnc -clobber\n",
1878,1907c1879,1908
<   0, 3761.673, 3887.281, 3472.979, 3282.3, 3688.367, 4092.374, 4016.597,
<     3932.583, 3795.031, 3634.005, 3617.12, 3808.21, 3992.71, 4135.204,
<     4278.11, 4130.262, 3923.935, 3872.867, 3951.939, 3999.3, 3952.763,
<     3906.637, 3985.709, 4068.488, 4096.492, 4098.139, 4105.552, 4124.085,
<     4122.438, 3948.644, 3773.616, 2756.803, 1558.372, 1236.731, 1678.215,
<     2236.659, 3150.926, 4064.781, 3755.083, 3395.554, 3867.513, 4738.95,
<     4939.512, 4081.254, 3305.775, 3250.177, 3194.168, 3290.537, 3423.559,
<     3542.99, 3647.595, 3690.014, 3533.518, 3379.905, 3450.74, 3532.282,
<     3569.759, 3585.821, 3615.061, 3665.716, 3697.839, 3541.343, 3386.082,
<     3501.395, 3683.425, 3729.138, 3646.771, 3590.351, 3619.179, 3646.36,
<     3553.698, 3454.034, 3452.387, 3503.454, 3524.046, 3492.746, 3438.796,
<     3092.858, 2746.095, 2904.239, 3197.463, 3329.661, 3303.304, 3302.068,
<     3386.906, 3469.272, 3465.154, 3454.858, 3486.157, 3541.343, 3549.991,
<     3474.626, 3430.56, 3829.213, 4224.572, 4589.867, 4949.808, 5077.064,
<     4960.104, 4894.622, 5021.878, 5146.251, 5205.967, 5258.682, 5258.27,
<     5228.206, 5159.43, 5014.054, 4871.148, 4755.835, 4640.522, 4571.334,
<     4514.913, 4481.143, 4472.906, 4364.595, 3829.213, 3294.655, 3419.44,
<     3613.413, 3629.063, 3533.106, 3496.041, 3581.702, 3662.421, 3624.533,
<     3587.056, 3603.941, 3638.947, 3580.878, 3415.734, 3268.298, 3206.111,
<     3145.16, 3293.42, 3468.037, 3582.526, 3656.244, 3714.312, 3739.022,
<     3766.203, 3888.105, 4010.007, 3852.276, 3599.823, 3569.759, 3804.915,
<     3978.296, 3854.335, 3731.609, 3865.454, 4033.07, 4059.427, 3989.828,
<     3970.471, 4066.428, 4159.502, 4147.971, 4136.852, 4094.021, 4039.659,
<     3951.115, 3819.329, 3687.131, 3555.757, 3426.03, 3770.733, 4186.683,
<     4351.828, 4337.826, 4180.506, 3680.542, 3187.579, 3587.056, 3982.003,
<     3392.671, 2423.219, 1788.998, 1593.378, 1591.318, 2714.796, 3837.862,
<     4080.019, 4178.858, 4173.093, 4088.255, 4043.366, 4097.728, 4153.325,
<     4130.674, 4107.611, 4013.302, 3890.576, 3850.217, 3923.523, 4001.359,
<     4099.787, 4196.979, 4229.514, 4250.517, 4208.51, 4115.024, 4032.658,
<     3982.003, 3932.995, 3992.299, 4051.602, 4024.421, 3962.647,
---
>   -5.345658e-11, 3761.673, 3887.281, 3472.979, 3282.3, 3688.367, 4092.374,
>     4016.597, 3932.583, 3795.031, 3634.005, 3617.12, 3808.21, 3992.71,
>     4135.204, 4278.11, 4130.262, 3923.935, 3872.867, 3951.939, 3999.3,
>     3952.763, 3906.637, 3985.709, 4068.488, 4096.492, 4098.139, 4105.552,
>     4124.085, 4122.438, 3948.644, 3773.616, 2756.803, 1558.372, 1236.731,
>     1678.215, 2236.659, 3150.926, 4064.781, 3755.083, 3395.554, 3867.513,
>     4738.95, 4939.512, 4081.254, 3305.775, 3250.177, 3194.168, 3290.537,
>     3423.559, 3542.99, 3647.595, 3690.014, 3533.518, 3379.905, 3450.74,
>     3532.282, 3569.759, 3585.821, 3615.061, 3665.716, 3697.839, 3541.343,
>     3386.082, 3501.395, 3683.425, 3729.138, 3646.771, 3590.351, 3619.179,
>     3646.36, 3553.698, 3454.034, 3452.387, 3503.454, 3524.046, 3492.746,
>     3438.796, 3092.858, 2746.095, 2904.239, 3197.463, 3329.661, 3303.304,
>     3302.068, 3386.906, 3469.272, 3465.154, 3454.858, 3486.157, 3541.343,
>     3549.991, 3474.626, 3430.56, 3829.213, 4224.572, 4589.867, 4949.808,
>     5077.064, 4960.104, 4894.622, 5021.878, 5146.251, 5205.967, 5258.682,
>     5258.27, 5228.206, 5159.43, 5014.054, 4871.148, 4755.835, 4640.522,
>     4571.334, 4514.913, 4481.143, 4472.906, 4364.595, 3829.213, 3294.655,
>     3419.44, 3613.413, 3629.063, 3533.106, 3496.041, 3581.702, 3662.421,
>     3624.533, 3587.056, 3603.941, 3638.947, 3580.878, 3415.734, 3268.298,
>     3206.111, 3145.16, 3293.42, 3468.037, 3582.526, 3656.244, 3714.312,
>     3739.022, 3766.203, 3888.105, 4010.007, 3852.276, 3599.823, 3569.759,
>     3804.915, 3978.296, 3854.335, 3731.609, 3865.454, 4033.07, 4059.427,
>     3989.828, 3970.471, 4066.428, 4159.502, 4147.971, 4136.852, 4094.021,
>     4039.659, 3951.115, 3819.329, 3687.131, 3555.757, 3426.03, 3770.733,
>     4186.683, 4351.828, 4337.826, 4180.506, 3680.542, 3187.579, 3587.056,
>     3982.003, 3392.671, 2423.219, 1788.998, 1593.378, 1591.318, 2714.796,
>     3837.862, 4080.019, 4178.858, 4173.093, 4088.255, 4043.366, 4097.728,
>     4153.325, 4130.674, 4107.611, 4013.302, 3890.576, 3850.217, 3923.523,
>     4001.359, 4099.787, 4196.979, 4229.514, 4250.517, 4208.51, 4115.024,
>     4032.658, 3982.003, 3932.995, 3992.299, 4051.602, 4024.421, 3962.647,
4789c4790
<  image-min = 0 ;
---
>  image-min = -5.34565829102274e-11 ;
4791c4792
<  image-max = 26989.4044242532 ;
---
>  image-max = 9456.07226562497 ;

So, it's either something up with that value of zero or this might
well be within the discretisation range a value of around 5000 for the
way mincresample does linear interpolation. (Macros, IIRC).

If you use -tricubic or -sinc you also get differences (the most with
sinc), given that the difference is also in the first voxel (corner)
there might also be small edge effects in some cases. As Vlad
mentions, -nearest does what you'd expect and there is no inherrent
blur (from the resampling) in this.

At this point I wouldn't call this a bug, I wouldn't expect two
volumes to be identical after resampling as there is always some form
of blurring in this unless you use -nearest_neighbour for an identity
transform.

a