You might be interested in our 3dfractionize program that does a voting method for computing something like the nearest neighbor interpolation. The voting method looks at the fractional overlap of different voxel values when going from a fine grid to a coarse grid usually and picks the largest fraction. See the help for 3dfractionize for more information on that.
I'm not sure any of these interpolation methods will provide exactly what you want and preserve small nuclei. Another program to try is 3dAllineate to apply the transformation matrix with "-final wsinc5" to use an alternative weighted sinc interpolation method instead.
Of course, you need not transform to a standard template at all, but instead use your data in its original space by using manually drawn or transformed standard ROIs.