[MINC-users] mincresample modifying data
Andrew Wood
andrew at biospective.com
Thu Aug 15 10:41:37 EDT 2013
Hi all,
I looked in the linear interpolation code, as Vladamir suggested, and it
looks to me like there are edge cases that actually get extrapolated.
Values that fall off the edge of the input volume, but inside a small
epsilon, are extrapolated. I tried removing the epsilon from linear
interpolation code and it seems to fix this case.
I tried another couple of things comparing the official mincresample with
my patched version, using short and ushort input
( data available https://www.dropbox.com/sh/kzziqemv3ro855q/0hkYHLKmzF )
$ mincresample in_short.mnc res_short.mnc
$ mincresample -clobber in_ushort.mnc res_ushort.mnc
$ mincresample_patched in_short.mnc res_short_PATCH.mnc
$ mincresample_patched -clobber in_ushort.mnc res_ushort_PATCH.mnc
$ mincstats -min res_short.mnc
Min: -5.36601874e-11
$ mincstats -min res_ushort.mnc
Min: -5.36601874e-11
$ mincstats -min res_short_PATCH.mnc
Min: 0
$ mincstats -min res_ushort_PATCH.mnc
Min: 0
With both data types, the official mincresample results in small negative
mins, and my patched version results in zero mins. I think this is
encouraging, and it makes me wonder if the epsilon in the out-of-volume
test condition should be removed. I noticed that the corresponding test
condition in the tricubic interpolation code doesn't have any epsilon
wiggle room.
To continue the experiment, I added a small linear transform and things get
weirder.
$ mincresample -transormation linear.xfm -like like.mnc in_short.mnc
res_short_xfm.mnc
$ mincresample -transormation linear.xfm -like like.mnc in_ushort.mnc
res_ushort_xfm.mnc
$ mincresample_patched -transormation linear.xfm -like like.mnc
in_short.mnc res_short_xfm_PATCH.mnc
$ mincresample_patched -transormation linear.xfm -like like.mnc
in_ushort.mnc res_ushort_xfm_PATCH.mnc
$ mincstats -min res_short_xfm.mnc
*** mincstats - reported min (-0.010563) doesn't equal header (-0.010563)
Min: -0.01056300123
$ mincstats -min res_ushort_xfm.mnc
Min: -0.01056300123
$ mincstats -min res_short_xfm_PATCH.mnc
Min: -3.637978807e-12
$ mincstats -min res_ushort_xfm_PATCH.mnc
Min: 0
It looks like both runs using the official mincresample ended up
extrapolating a minimum of -0.01056300123. I'm not sure what to make of the
difference in the patched runs. It looks like the numerical stability of
some expressions in the linear interpolation code might depend on the
magnitude of voxel values they're given. The the following lines, the v***
variables are voxel values from the input image.
*result =
r0 * (volume->scale[slcind] *
(r1r2 * v000 +
r1f2 * v001 +
f1r2 * v010 +
f1f2 * v011) + volume->offset[slcind]);
*result +=
f0 * (volume->scale[slcind+1] *
(r1r2 * v100 +
r1f2 * v101 +
f1r2 * v110 +
f1f2 * v111) + volume->offset[slcind+1]);
- Andrew
On Wed, Aug 14, 2013 at 8:20 PM, Andrew Janke <a.janke at gmail.com> wrote:
> 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
> _______________________________________________
> MINC-users at bic.mni.mcgill.ca
> http://www.bic.mni.mcgill.ca/mailman/listinfo/minc-users
>
More information about the MINC-users
mailing list