There are many options yet it is highly dependent on the properties of the image being resampled. Images downsampled with bicubic interpolation are regarded as smoother and with fewer interpolation artifacts. You may also try nearest neighbor interpolation and bilinear interpolation. Make sure a low pass filter is applied before further processing with downsampling strategies, otherwise the aliasingaliasing artifacts may be induced.
Besides, you may also consider vectorization algorithm in which a resolution independent vector representation is created before rendered to the image at a desired resolution. it is mainly used in up-sampling but the vectorization method is helpful in preserving the feature connectivity in down-sampling as well.
Another interesting blog on comparison between different image re-size approaches.
EDIT
Regarding the design of low pass filter served as a pre-processing method before downsampling, you can consider this way:
Suppose a low pass filter $g$ works on the image $x$, and the filtered image is $u = x*g$. After that, simply decimate $u$ by 2, you get $v$. Now you want to recover the original image. You can zero-upsample the $v$ into $w$ where $w_{2n} = v_{n}$, and $w_{2n-1} = 0$. Then another filter $h$ is applied, trying to obtain $\hat{x} = w * h$. After Z-transform, $U(z) = X(z)G(z)$, $\hat{X}(z) = W(z)H(z)$.
Since $W(z) = 1/2(U(z) + U(-z))$, we have $\hat{X}(z)=1/2(H(z)G(z))X(z) + 1/2(H(z)G(-z))X(-z)$
A good downsampling is equivalent to a good recovery to some extent. In order to best recover the original image ($\hat{X}(z)=X(z)$), you are trying to make $H(z)G(z)$ as close to $2$ as possible, while $H(z)G(-z)$ as close to $0$ as possible. The design of low pass filter $G(z)$ is based on this criteria.
