psFitsImageWrite BSCALE/BZERO needs some work

[eugene]

There are a few issues going on in psFitsImage.c that need to be cleaned up. I have disabled a bit of code with #if(0) until these are resolved.

Currently, it looks like the steps for CFITSIO output are as follows:

• the file is opened for write
• compression options are applied to the fits pointer with psFitsCompressionApply
• the ps-specific float type is written if selected, otherwise:

(I get a little lost in the options at this point : are these next two steps only for float types?)

• a new value for BZERO and BSCALE are determined based on the data range (confusingly, this is applied for write, but not update?)

or,

• the old values of BZERO and BSCALE (from the header?) are applied.

• the image has random fuzz added to it (if input is float and output is int type)

The problems:

1) the user needs the option of control over BZERO and BSCALE. for example, you may know that the realistic dynamic range of the output image is (say, 0 to 100k). The autoscaling trick can be very dangerous because of division by near-zero values.

2) the user needs more control over the application of the fuzz when converting from float to int. It is necessary (or useful at least) to apply the fuzz when writing a float-valued image out as an int, but in some cases, you have an image being carried as a float which is actually an int: using ppImage or ppMerge to manipulate masks, for example. In that case, you want the hard edge transition from one int value to the next.

Also, I find this code very confusing. I think it does not help that the steps are nested rather than sequential. So, for example rather than define a function 'convertImageWrite' that calls two possible write functions, one of which calls the other, I would rather see these as isolated steps, eg:

scaleImage
convertImage
writeImage

the way these are laid out, you cannot tell by looking at the top level what is going on: you need to keep too many of the separate functions in your head at once.

[Paul Price]

0. I think that the confusion in the work flow has resulted from poorly chosen function names. For example, "convertImageWrite" did not actually write an image, but merely converted the image *for* writing to disk (as opposed to converting the image after reading from disk). To reduce confusion, I've renamed several of the functions and added to the comments in the prefaces.

BSCALE and BZERO are generated from the image data only for floating-point types. However, BZERO may be set if the image in memory is an unsigned integer; see bug 1023.

In answer to the question about the difference in setting BSCALE/BZERO between writing and update, I figured that we do NOT want to generate a new BSCALE and BZERO when updating, since it's very possible that you're updating only a piece of the image (as we do, e.g., in ppMerge), and generating a new BSCALE and BZERO would mean having to touch all of the pixels already written to disk.

1. There are several ways that we could choose BSCALE and BZERO. At the moment, the code attempts to preserve the dynamic range. We could choose to preserve some minimum scale based on image statistics (e.g., 0.1 * standard deviation) or some user-defined value, at the (possible) cost of killing the dynamic range. We might also desire to have the user set BSCALE and/or BZERO explicitly. It would be nice to have all of these options, but I fear making the interface (at the psModules level) too messy. Which of these options should we support? Which should be the default mode?

2. Integer images (in memory) are never fuzzed. I'm particularly confused about your statement about fuzzing mask images. They should always be U8 in memory, so there would never be any need to fuzz them on write. You mention using ppImage and ppMerge to manipulate mask images. ppImage should never be used to 'manipulate' a mask image (if by that you mean specifying the mask image as a "-file") because it converts it to F32 on read (you've told it that it's an image, so it will treat it as an "image", not as a "mask"); I expect that ppArith with the "-mask" flag will accomplish what you desire here. ppMerge generates masks properly, at least under the modes that are for generating masks, otherwise all bets are off. Can you please give me a concrete example of what you want to do?

[eugene]

regarding (2), there actually two separate questions:

1) the fuzzing of int-valued data (masks). detrend_norm_apply.pl uses ppImage to construct the normalized image. This treats the mask as an input -file image, calculates stats, etc, and writes out an image with a normalization of 1.0. I agree that the mask should not be used as a -file input image. The correct fix here is to modify detrend_norm_apply so it does not apply ppImage to the masks. This could be to copy the image, or if we want to inject information in the header, then we need to use another tool to perform that operation.

2) should the 'fuzz' be optional? It is certainly necessary in some cases, but there are cases where you would not want it. The easiest example is in creating a tool that reads in and writes out an identical copy of the image, perhaps with header modifications. For example, a compression program should not alter the real data values.

[eugene]

regarding (1), what options we should be supporting. I want to see the following:

psLib: at the psLib level, we need the following features:

• user-defined BZERO,BSCALE. the user should be able to specify the output BZERO and BSCALE values as desired, accepting that this will potentially sacrifice dynamic range.

• auto-scaled BZERO,BSCALE. the psFitsImageWrite APIs should be able to choose BZERO,BSCALE in any of several ways: maintain dynamic range, sample sigma with Nbits, there may be others.

These options can be easily added by defining elements in the psFits structure equivalent to the compression options. They can be set at the psLib level with APIs equivalent to the psFitsSetCompression function.

psModules: at the psModules level, the above APIs should be used to set the desired options in pmFPAfileIO.c, just as psFitsSetCompression is called. We should extend the camera.config:FITS METADATA structure to include fields for BZERO,BSCALE, with allowed values that specify the auto-scaling mode. eg:

BZERO BSCALE
0.0 1.0
AUTO FULL_DYNAMIC
AUTO FULL_SAMPLE

[Paul Price]

Thanks for the guidance. I believe I've implemented all the desired features at the psLib level, though the exact method of specifying them at the psModules level is a little different.

I recast the psFits structure to include a pointer to some psFitsOptions which can control how psLib will read/write FITS files (if non-NULL). The idea is that at the most basic level, users will simply want to let psLib do its thing (leave the fits->options = NULL), while people who care how the data is written can go to the trouble of setting the options explicitly.

At the psModules level, I broke the old manner of specifying the FITS options using a "TYPE", and resorted to separate METADATAs. This is because not all options are required; some options are used in combination while others need not be mentioned at all. Below is an example.

FITS METADATA
# BITPIX is the bits per pixel for writing the output data
# COMP = NONE|RICE|GZIP|HCOMPRESS|PLIO is the compression algorithm
# TILE.[XYZ] are the tile sizes. 0 means entire the dimension, so (0,1,1) forms tiles from rows
# NOISE [0..16] is the number of "noise bits" to preserve when quantising floating point data; 16 for no loss
# HSCALE is the scale factor for lossy compression with HCOMPRESS; 0 or 1 for none; 2*RMS --> 10x compression
# HSMOOTH is the smoothing to apply to HCOMPRESSed data when decompressing; 0 for none

# BITPIX(S32) is the bits per pixel for writing the output data
# COMPRESSION(STR) = NONE|RICE|GZIP|HCOMPRESS|PLIO is the compression algorithm
# TILE.[XYZ](S32) are the tile sizes. 0 means entire the dimension, so (0,1,1) forms tiles from rows
# NOISE(S32) [0..16] is the number of "noise bits" to preserve when quantising floating point data
# HSCALE(S32) is the scale factor for lossy compression with HCOMPRESS; 0 or 1 for none; 2*RMS --> 10x compression
# HSMOOTH(S32) is the smoothing to apply to HCOMPRESSed data when decompressing; 0 for none
# SCALING(STR) = NONE|RANGE|STDEV_POSITIVE|STDEV_NEGATIVE|STDEV_BOTH|MANUAL is the scaling scheme
# BSCALE(F32) is the manual scaling to apply (when SCALING = MANUAL)
# BZERO(F32) is the manual zero-point to apply (when SCALING = MANUAL)
# STDEV.BITS(S32) is the number of bits to map to a standard deviation (when SCALING = STDEV_*)
# STDEV.NUM(F32) is the number of standard deviations to the edge (when SCALING = STDEV_NEGATIVE|STDEV_POSITIVE)
# FLOAT(STR) is the name of a custom floating-point type

DET_IMAGE METADATA

BITPIX S32 -32

END
DET_MASK METADATA

BITPIX S32 8

END
DET_WEIGHT METADATA

BITPIX S32 -32

END

COMPRESSED_POSITIVE METADATA

BITPIX S32 16
SCALING STR SIGMA_POSITIVE
SIGMA.BITS S32 4
SIGMA.NUM F32 10
COMPRESSSION STR RICE
TILE.X S32 0
TILE.Y S32 1
TILE.Z S32 1

END
COMPRESSED_MASK METADATA

COMPRESSION STR PLIO
TILE.X S32 0
TILE.Y S32 1
TILE.Z S32 1

END
COMPRESSED_SUBTRACTION METADATA

BITPIX S32 16
SCALING STR SIGMA_BOTH
SIGMA.BITS S32 4
SIGMA.NUM F32 5
COMPRESSION STR RICE
TILE.X S32 0
TILE.Y S32 1
TILE.Z S32 1

END

END

I have merged my branch into the mainline. I think that concludes this bug.