[gdal-dev] Pansharpening produces grayscale image?

[Enrico Mattea]

Hello all,
I'm trying to perform pansharpening of Landsat 8 images using the 
gdal_pansharpen.py script from 
https://svn.osgeo.org/gdal/trunk/gdal/swig/python/scripts/gdal_pansharpen.py, 
with mixed luck.

The command
gdal_pansharpen.py 8.tif 4.tif 3.tif 2.tif ps.tif,
where the tifs are the 16-bits original landsat bands, runs with no 
errors but the resulting image is - according to both imagemagick and 
the gimp - a 16-bits *grayscale* image instead of the expected 16-bits 
rgb image.

The tile I'm specifically working on is LC81950282016101LGN00.

What am I missing here?
Thank you very much,
Enrico