[Insight-developers] Re: Bug somewhere in ResampleImageFilter?
Lydia Ng
lng at insightful . com
Wed, 25 Jun 2003 21:49:49 -0700
This is a multi-part message in MIME format.
------_=_NextPart_001_01C33B9E.5F8755F5
Content-Type: text/plain;
charset="us-ascii"
Content-Transfer-Encoding: quoted-printable
Hi Mark and Luis,
A couple of years back I did compare the current version with a version
which explicitly reduces it to 1D lerp (hard coded for 3D) and as I
recall I didn't get much of a speedup.
I did play around with the interpolator this afternoon - using both
reduction to 1D lerp and Stephen's trick. It seems to be actually slower
- but it could be my implementation. I have attached the code for those
who would like to play with it :-)
Cheers,
Lydia
> -----Original Message-----
> From: Luis Ibanez [mailto:luis . ibanez at kitware . com]
> Sent: Wednesday, June 25, 2003 8:42 PM
> To: Lydia Ng
> Cc: Stephen R. Aylward; Mark Foskey; Derek W Cool;
jisung at email . unc . edu;
> piloo at email . unc . edu; Insight-Developers (E-mail)
> Subject: Re: [Insight-developers] Re: Bug somewhere in
> ResampleImageFilter?
>=20
>=20
> Hi Lydia,
>=20
> I agree with you in that interpolation in a binary image doesn't make
> much sense. However the same rounding error will degrade the
> interpolation of any image with integer pixel type (e.g. char, short
> long, signed or unsigned).
>=20
>=20
> Extending the trick to N-D should be possible, however the
> code in itkLinearInterpolateImageFunction.txx will have
> a very different look.
>=20
> Currently in itkLinearInterpolateImageFunction.txx we
> are computing the full contribution factor for each neighbor
> of the continuous index to be evaluated. We end up computing
> somehting like
>=20
> interpolatedValue =3D 0
>=20
> for all neighbors
> {
> interpolatedValue +=3D K1 * K2 * K3 * K4.... * KN * =
neighborPixelValue
> }
>=20
>=20
> Where "N" is the image dimension and the terms Ki are
> "(overlap)" or "(1-overlap)" along dimension "i".
> The number of neighbors in the delaunay region used for linear
> interpolation is 2^N.
>=20
> The drawback of the current computation is that we are not
> taking advantage of the potential factorization in the products
> of the Ki terms,
>=20
> So in an N-D cell we are computing
>=20
> N * (2 ^ N) Multiplications
> N * (2^ (N-1) ) Sums
>=20
> (e.g. in 3D =3D=3D 24 mult and 12 sums )
>=20
> when in principle the same task could be done with
>=20
> 2^N-1 Multiplications
> 2^(N+1) -2 Sums
>=20
> (e.g. in 3D =3D=3D 7 mult and 14 sums )
>=20
> when the terms are factorized.
>=20
> The implementation will require a bit of auxiliary
> memory in which the interpolations will be done
> along a single direction at a time, and each time
> the dimension of the problem will be reduced by
> one.
>=20
>=20
> Since this is the most commonly used interpolator
> in the registration framework, it is worth to give
> it a try to implement the factorization.
>=20
>=20
>=20
> Luis
>=20
>=20
> ----------------------
> Lydia Ng wrote:
>=20
> >Hi Stephen and Luis,
> >
> >I am not totally convinced that we should classify the "rounding
error"
> >as a bug.
> >
> >I am not quite sure that using a linear interpolation for a binary
image
> >makes sense if what you really want back is a warp/transform binary
> >image. One way to this more robustly is to create a signed distance
map
> >from the binary image, then warp/transform the map and threshold to
get
> >the warped binary image - ala
> >S. P. Raya and J. K. Udupa,
> >"Shape-based interpolation of multidimensional objects",
> >IEEE Trans. on Medical Imaging, vol. 9, no. 1, pp 32-42, 1990.
> >
> >
> >
> >
> >>>If instead you do
> >>>v =3D pixel1 + overlap * (pixel2-pixel1)
> >>>there are fewer multiplications (faster?) and the rounding error
> >>>
> >>>
> >will
> >
> >
> >>>not occur.
> >>>
> >>>
> >
> >I agree this would help with rounding error. Does this trick extend
to
> >N-D?
> >
> >In any case, if we are to replace/modify the existing linear
> >interpolator we really should do some performance testing so that we
can
> >systematically compare before and after.
> >
> >My 2 cents worth ....
> >
> >- Lydia
> >
> >
> >
> >
> >
>=20
>=20
------_=_NextPart_001_01C33B9E.5F8755F5
Content-Type: application/x-zip-compressed;
name="Lerp.zip"
Content-Transfer-Encoding: base64
Content-Description: Lerp.zip
Content-Disposition: attachment;
filename="Lerp.zip"
UEsDBBQAAAAIAPB62S7yEDrWmgIAAHUFAAATAAAATGVycC9DTWFrZUxpc3RzLnR4dKVUXWvbMBR9
L/Q/XLJC3TE89rA9jYFry5kax84Ue+xpRrUVR1SRgj/ajtH/vis5SdstjELzYKyTq6N7zzlyOA9m
pJzTlM6LecnIt4IyEnnfCVvSLIUP/qeL05MFy65ImHuJaLe4pGmYFBHB6mmRBKwkPxaMLG29N/np
vz2bYM3pyRvIhVIQzvmNOO8gpmlULoJwFkwJVGaz4bqGu7VoBfQGlDE3sDItmKGFWOr68/am+WJJ
NqYelOj805Mlyb1wbDeLioQgW/4Vzn7/gz0gSHUnm3UfbLdKVryXRnflMitYSMqIsoddh5HR5z1c
D1LVwPWvfi11A4PG8zro1wJuRdvhVjCrcQ6QHayRF4Q2Q7P23ZRYt+F95faSZEnekzSiMXRrMyDt
tXBMinc9KKlFB1I7ZCWVQAIae5PDDMFVxsqd9g/+AabpE3gCU0aCnDA056ObA7uwjTjdgOYzsNJW
Rq9kM4zyDh223vuu8KkRHlZfuB7wpYyzIo1wCbBz2EMdES+WpIxpQqxqdsC/iufoveWKgzxISsJY
xiy8/01CrrXZq4wSd8AfXYE72aNQve3bB1gowbHXTjjAWuXbNDlFnx07+hfUtdNS3Itq6Pm1Eo7P
pajDR3WQOYgizCkJizy4TIgn+5sEzeAt1T2G2ijeC7rhjYgHXdm+coF+vajKr+7vXUN5wKYkLxOa
zvBxyQJGydJ7GYkTPZ8FqjEtDrDpYIdc8k5WsVS49YClw0a0shrXSH+r1f4vmu3fQrxiKK9dHeSi
q+OpvuMYd2PAqPqdvZtbxTGkNkWjaWifTfc+7tag/4YbLd7J7uLyunw/BmyJMTiSr9yed+ym7kay
F2Cc48hHATM3eca2S+CwbVpeu8sz8uGX8EkUXzfSH1BLAwQUAAAACADqfdkuNRbuvbQFAABoDwAA
KQAAAExlcnAvaXRrTGluZWFySW50ZXJwb2xhdGVJbWFnZUZ1bmN0aW9uMi5orVdtT9tIEP6OxH8Y
tRVKKDiFU3VqaCulFNSoEFBirlepEtrYk2QPe9faXQdQxf32m9m1E4eCQtHlA6zXM8+8v7iz/eH/
+m1ubG4AnBs9NSLvAkBfWTmdORjhNEflhJNawRYMcSqtM+Ex1jq7ko45T3VaZsiM8Gp4OJpIfpDu
6kQqFKavHJpCZ8JhPxdTPC5Vwgj70WxnDq8Y4ESoaUmvCOLw9Wu++UzUHpAgw3n/zZv9zpt3nb09
2HvXfftnd++PwPwXGktwXS8d5zI87EX7b/k9Uxzq4tZ4g1pJ2wMtLDzUymrjZJlH0Msy8GQWDFo0
c0wjZh8hQj/+ukCJ3I0DbWDmXNHtdK6vryMyNtJm2vkSn550loQzl8OEKFN0QmY2CurQL55JC1ZP
3LUwCHRO2bFyXDpM4Vs//nJ2EUNv8B2+9YbD3iD+fgDX0s106QDnqKCCcTNizotMEhchGaHcLegJ
nB4ND78QW+9T/6Qff2dlj/vx4Gg0guOzIfTgvDeM+4cXJ71hDXV+MTw/Gx1FwVwGFmM9R0gWvlPa
yQStNyjXrLaiY+6zwVu2PtGe+NvubG68lBOV4gQu1+fR5YzIiZaonkrO+CrJyhThBXE8RhvNXjCt
EjnaQiTISb258ZPvOtvb8CPJhLWwRh45mEjHRpI1gTS7Jd8tqEEoCiJxgHBgC0zkhONZaCuZn7MG
tj3IGkGQPYBuK2y+opyn/BDOgwkKqNrl6ykaKOQNZguZUUjQYB4dHFKSEVjqOSktTEg9VVBCBnx3
W7AlqX+RaG1Sqdg4gwUXU91EPJlHaWE0jWCSaeGLKdXlOMP20livwaSyDa61uQqpN9hNZc6maCWy
ILzhoh9UBkqq6QPsoBW5hiECk68osInIRGW/R2ANbQTHRDfHxC3JS4vUa/hmXV9bamPF01iWHKS7
0WUBK+8tv6OiqOMA70No4j5HwJPuVNGKD9n3QyzgQ3Dux82N8IZa2OXR3+dnw3htInV9VygoIjKB
x8jeN6UvxH4EXyCBt8swXCkjR6khTFopyT6mgiU3s1VQP6/Va4TZ5KDJ8Lu6jcoCjddhBWaUC+PO
tS+b9yzkI1RPj5MlFBYHgdiPEddkqQ0flmrXUb6GxG90TGhxsRj0ZQU5UnNPbbt2CPWZmBhORWJ0
a41Xdh71QvsgDBxW5NRL8OmfGKyKcUapNp35ktXjfyhRYSI4W28begzwOqjBprYP2LQa86x05GRW
dOEiWxYFzdP7geX/3EUbAeh2G+zLY0PlZQifLeIexOpjQ9QQRfZsIQvm+tAA/lz3KihpnFFv5t7U
6JsNR4+4RSY+k4LDvaILgB1CoJ1FUbZQmu2sWLlC2F5xYYo3z/QccVZOq04NXNLSSVXq0j5fwj0M
L+uBu4bUo7nISm5/nLCL3k7zQ/yiUD3JmNW3VvpL8XGlUdazPzQr0wdmJWFX3M3RTBSNYTnQMNYU
X5qYM0yufIhppdMKo4o3JomBi16Q9SV3BE1KoJ9CUlWrHEsfl5MJmmhF85WypuBY3l9T/ORJW21I
aIcYI4+olHG9GrVOY6SixwqIpl3G+rG40HTqWM2lcSVN1EZV1v7uuXvubVVbY+iBDwVyi5zIcWgH
Gh/EwmhHPQZTPxbWdLVW2/fRf9eS/bzzhHMtU/qEIRdzn2pZl3a7mnZqFPkWaLvj05j97/811aLi
U5wN1QoCJySpVZlWXW0BewL3dh683v/lWmTFTNS2B1f9pF7P2QcVFLyuqGgTqmBgt35HbfYueEzO
+ePnKf6C5UDaYoBOpyiJ0iLlOK3u/jsB+WvOb3HBX5qqUVC///B73HU9Dsp8TKsgfXAopE+EsTY2
pGCdzwtluU5Dmlnf5Spxi46WacpJyC8HNU41Z+58hO5IHbi3hC+/EHivOe0NLnonl/3BKKbPnn4v
7p8N7u3465Y2d3NDy/5LVKmcePjq9B9QSwMEFAAAAAgA+IDZLmkzqhiYBQAARA8AACsAAABMZXJw
L2l0a0xpbmVhckludGVycG9sYXRlSW1hZ2VGdW5jdGlvbjIudHh4rVd/b9s2EP0/QL7DNRsKOU4k
xxuG1U4CeG6yGkuTwHZXFEUQyBJts5ZITaTsZEW+++5I/XSQpWsmILBE3j2+e3c8Mt7+yf/17O7s
7gBcp3KR+nEPAEZC8cVSw4QtYia0r7kU8BrGbMGVTu3nVMpoxTV5vpdhFjFyhB/Hw8mc0wfXqwsu
mJ+OhGZpIiNfs1HsL9h5JgJC6Lr67u5gDT8SxIUvFhlOIsiw3aaRt2hvIBHUvnc7nZ+8zi9etwtH
b3o/v+kddazznyxVCNgz67M1tx9HbvdXmieLoUzuUxOSE7QIqFvGOJRCyVTzLHZhEEVgzBSkTLF0
zUKX3CeMwWj6R4mCzDXIFJZaJz3P22w2LobrynThvZu+v/Aqw6WOYY6WIdM+j5Rr6eAzXXIFSs71
xk8Z4HtI0vJZplkIH0fTd1cfpjC4/AQfB+Px4HL6qQ8brpcy08DWTEAOo5foHCcRRy9ESn2h70HO
4f3ZePgO3Qa/jS5G009E9nw0vTybTOD8agwDuB6Mp6Phh4vBuIC6/jC+vpqcuTZcAvZncs0gKLUT
UvOAKRNQLIm2wNfY1IOJ7PlS+8Zn39vd+YHPRcjmcPt8Jd1iJaEDWqPdf3DgIoiykMHeNxTrco8i
rFzWIvLw7xbDX5rJBppxpv0iBRWYHiGsr2Wawwg/ZirxA0bbZHfnK415+/uYiX14a6OgBIgsnrGU
8ikYZmAmU0UmKI5mmHQkeRxEvlIwHYkk02bRA8iHhlKm4Zglp7s7AVGAjCpeYKFEUix2d54J+LgJ
WsIB4vV68e1lwQhO4AiOG+a9nvl5y7F70GbsU3xVhEaRNAtQj+8K50XMn3F2Wnk6Hpqcr1Mu9IRF
8+9ivJY8fCHvkoCjdNjrSVSQ+fFrkOoAe1mIbRq3I/20wKTbhEENgqvD00mWsNTQqgOha+7SfxTv
2dqPMqQIPuISIYN+B4lUnNh+lwz6PmFU+vAiJa4yjXNTxHqhpEWMA40VqbnIZKZIyjuHlLObZmuG
Vn1tpWjqXG4uFBd7edwH8LxcM+yiKY3ZzWCPAKszmN0QI0MGM18VIp+gdlIxnX/OWCQ3qDxCb/nQ
oeELbCPzVMbWArSsQVl7z1RCIb4RgwLBbVrEZFzMV59sQ5nNogr+c3M73xgb6iKgCx8awcPAoTiR
f6dvXo5hqxHQaLsNLXvmfLU/KNQUu91cRhgnR9Q9Pt+DWSSDFZ2M7K+MY6aYjU03LCNqlfnBWIIt
/TXNoXHg44k+jzDrbjmLTxnsZ6Rzg3QdCqZlLR2rnJ1q9a0fn4NTGz7FCN1OHkUZx9PAlWu/MH2w
LyxS7BGMVfVf3YmQzZJjrVvw6qRmWHKrwRbAh4f9aujhSe5lbhuGZU3kVjVVDvO6cbaxChUNhKl+
Dy8ZeCIl4IOKKUmpOSvxOLTn5FMFO+F/m7fcnj4NdvXpnvMocqBrFn0Kxp7MNaAaiIvMCMipgeYR
VAa2TVRxNlajW8tTF4DjksUpFGMGu/goJXAgvs37Vi6OJUF/jW4jMAv189iDrjE058Sa0TF7bBNz
StfaLNLKgQKtgRTIjDqp2b82nM2SUz5fFaTckRroMxE6rcd72E+S6B72ImzFe7T5mB8sIfF5ajYo
LZ0gmWIPdw5xrNYUTW3ZtkNtmR0hC0X/egS3ga90EYFT6uP+zjTRyGur3eZ1OaEI9XMeVLttatoc
hxdI0cmXOSjuv9+22kG1ATo3Ty/+kHd5D2kEGf57smaoTUMhqnYiFrOQ00Gb8zW3a5XNFLY9Zs8S
KxFd/cykUyXtC12++vjzuNF+edRmAZPuneTFQU8z3XbM9PASf2Vb+QrxEXJVQTZgK61X2xpvJeHm
4FFajm5qkn6pJK0TbNdZl52o7CYp01kqSuTOTb+8wD1QCrZu23hVZyLkeJP7B1BLAwQUAAAACABU
gtkuqozA5U4IAADGHgAALgAAAExlcnAvaXRrTGluZWFySW50ZXJwb2xhdGVJbWFnZUZ1bmN0aW9u
VGVzdC5jeHjtWf9P20oS/x2J/2HKPVUGgvMF8npNAifKo9foKEVJ+p6q3lNl7E2yxfFG3jWEq/jf
b2Z3/TUGSo8frrozIrG9M5+dLzuzM5vmzuFzXZsbmxsAF7GYxd6iBwDDSPLZXMGYzRYsUp7iIoKX
MGIzLlVsHidChFdcEed7ESQhI0b4ZXQynnJ64OrqjEfMi4eRYvFShJ5iw4U3Y2+TyCeECZPK9Ver
xjX8QihnXjRLcBxxTnZ36c1vyKJREdfcd1qt/War0+zsQ6fTOzjodVuG+XcWS8TsaRHYNTcPbbfz
isaJ4kQsb2OtleNvE1AnU/NERFLEiicLF47DEDSZhJhJFl+zwCX2MWMwnPwjQ3HVSoGIYa7Ustds
3tzcuKixK+JZ893k/VkzJ5yrBUyRMmDK46F0jTh4TeZcghRTdePFDPA+IOvyy0SxAP4YTt59+DiB
4/NP8MfxaHR8PvnUhxuu5iJRwK5ZBBZGzZF5sQw5ciFS7EXqFsQU3p+OTt4h2/Gb4dlw8omEfTuc
nJ+Ox/D2wwiO4eJ4NBmefDw7HqVQFx9HFx/Gp65Rl4C9S3HNwM9sFwnFfSa1QgtBYkd4u9BLQmv2
+Gr7zmunSXB/4ZEfJgGDARdoHOYtjkqvt9DoelW5863K6zH/V83bh9dkZ53hHMlxpZ4z1P9S3Mup
GSusE75gGFWXTJ6IMGS+EvEbT1aFuo7CJv5/kR46ka3hpMNo47kd3NxQt0sWsCkFWa9Hig72j+Dh
i6gmyNavcGstBgmFQoQrSM4xEhpQhdNUFfZAJJchq5+tcp0IEQcjtqwT4BGPDLKpG0WUI8pSKYuI
K8DVIVQzCtiKbtdky0Ye5L8QPFK1/NnIg/yYYxSPEpHIqiQ1I33yceQtmFx6Ptvc+EbxubnR3IHf
2BTNZWMe7QIS3QpeFMByfiu574UYrGglHqEhJVAQpY43pIfw7Vun1YAD/P9r6+4Op7J+FBjhPILP
+38SEbTc7lkDxSt+F6hJMh7Ncuq2pWp1z6BB3x1Dv7lxRx/NnR1UYgco64MHMyYWDJOdD0syn4t7
i0riSIKKE1RtCoro5p6EpSclCxqaObZEUy+UKDBaIb7hktEYKqoYBhCqDQM/RCaYFL2AaeMSdyw9
/9/TybXrnM0NHzcAVabH2bh+bKTDmZ9fGpkbFpFL2kcC1shsQzr87oWoyLb2HXlPqqDX8yl3Dwaw
ZcB6eINPGq1vyDTi5bVmPrQS7B0N7RRvkumUxY7hgO1+DbAhtMgGyELzqZMivzjMpCYRaRV+M19l
tJ2dHTiNYxHjdm6oMT0kYQCXzEyQoeC9ZmVREPYNlHGW8ZV9dZdLcs/81oJVA5zSC3RtWfU15bXR
re5F1c2UaRb94l1Kx06xV3QWHEGb7b1ORcqEesAsBqVilRyyxiy1htGmKVioNF8ZwXLTHPdGl5+l
FLQg5pRnDa/HoquSz743vGrS4Esj/X8UaBrMLgmN9qRA0xz/i4F2vObGki3+H3lp5G1uUD5aeDxy
rL74SGWjUB8xsBwvnvnbDTPkzz298Euj19vw+U+7flP+aejN0CmtPuQXbf48wJ0fGwUdvUvTr8GN
iK/+BqZcXrcYRSTu07Za4FkIUguHm+BcBLL3z2grcxdOgy2Q8Gkf9bBkoJrUzGcQaB6iy6oyWxqx
2BIcFofO2Y2TrpjCa2oksQGkoiTWt9lKpgd3zBQVLQ4VLCm3Bt87wqEzNBrKcyGk5LiEDZZjWNep
TSSz4CGyVGEtamYGZE5rIqqvbMFDndW9pVZ56g+a2TEY65KNDaBjgUtzDyOuuBdSwZaXepTVsRvP
p8oryTSz0VtqzRxaRTLEVs0sI3M7wIrP3u/ubucxBob/c4dqOT3cz4cyuFjcGDC6GWABqe8qQDlY
m8CQoF8ezOB8zL8ajm4G2JDruxq4HLJFkEjUryHIzHrBVyx0zMZhdSUxkW27wneXP5ZzgHbBCTab
ypTWvOKN0j6W+aK2Y6Cw0CMUF2skxfBIMzCqMIyWidLh4mi9qiQXMZWthUAvrR1dAzCTxynKPVWp
B3QDT3r5IsbWdinwVTTLjDFjxco8X23lvJJuFAhv03413dZ1PH5xoeYtVVoB01u9MZuqJNeqifJq
9Wdo0KWQXGtmt8osQgyx3dkCOgq5bacNShv7HdP33NmJjCw4WCOoY7lTyxt5kLamwBlU/XrkpLWN
naOh94wGvGoVPEX74guLu51m/Ladr5QrJrEXSTplqQoqLmwiMn1MNlu5Ri7JXm59HhLd9Dg/KDk6
LPMSrr88h11iwmRxjaM61lE/4qbO87rp1+9S9r/NRU+RuugeOny53zlsKXmIVMjtttheu9Wv8d2+
9V37NVZ1luMH3Lj/zNH2+md041OkLroR8/KjyfDA+qnzI0F28Fze0dV3A36yIHu61LRp0YF5tMer
G1eNb7rWN123021Au+Pi50HHfYqHus8bP93X7qvuzxhCTxWcDqHvOaqnatM8pW7QZxXZaTmJfk5o
LXNZqvw4/wv9iISNHnp0/6D7a6an9b3n+8kC+VtuK5Mnm9IdKy9WztaZEMv2Vq4PWtcpi3AFupK+
wjr6HL92d2sa6bpybJWWYnWgXw3oVxjs42cRE4ol+mr1+Sut3ILOScRpBTiAOYZ6OCLYg3amPF1Z
8f30g4jVqoDUbJbqUhzT5xGHhaOI+rMAY/rdw+y4IpNr3Q1iWfRCeXTElgK9lLunJI+ZpUYALNTN
mRz9zJiEWGcLeyqXNZJTcPRifYF+KHVF62ciYwwUwz31eMiC2oKcrvQEBHH7qQ8YnfbdC08/iBKy
BJn42DnIaRKGtzZQH5un1TcGtUcl63//BlBLAwQKAAAAAABUgtkuAAAAAAAAAAAAAAAABQAAAExl
cnAvUEsBAhQAFAAAAAgA8HrZLvIQOtaaAgAAdQUAABMAAAAAAAAAAQAgALaBAAAAAExlcnAvQ01h
a2VMaXN0cy50eHRQSwECFAAUAAAACADqfdkuNRbuvbQFAABoDwAAKQAAAAAAAAABACAAtoHLAgAA
TGVycC9pdGtMaW5lYXJJbnRlcnBvbGF0ZUltYWdlRnVuY3Rpb24yLmhQSwECFAAUAAAACAD4gNku
aTOqGJgFAABEDwAAKwAAAAAAAAABACAAtoHGCAAATGVycC9pdGtMaW5lYXJJbnRlcnBvbGF0ZUlt
YWdlRnVuY3Rpb24yLnR4eFBLAQIUABQAAAAIAFSC2S6qjMDlTggAAMYeAAAuAAAAAAAAAAEAIAC2
gacOAABMZXJwL2l0a0xpbmVhckludGVycG9sYXRlSW1hZ2VGdW5jdGlvblRlc3QuY3h4UEsBAhQA
CgAAAAAAVILZLgAAAAAAAAAAAAAAAAUAAAAAAAAAAAAQAP9BQRcAAExlcnAvUEsFBgAAAAAFAAUA
gAEAAGQXAAAAAA==
------_=_NextPart_001_01C33B9E.5F8755F5--