Short version

If you have a asinh compressed stack (C) and want the linear version of it (L):

 L = BOFFSET + BSOFTEN * (exp(C / alpha) - exp(-C / alpha))

where BOFFSET and BSOFTEN are supplied in the header, and alpha = log10(exp(2.5)) = 1.0857362.

If you want to compress something with asinh compression:

  BOFFSET = mean(L)
  BSOFTEN = sqrt(alpha) * stddev(L)
  C = alpha * asinh( (L - BOFFSET) / (2.0 * BSOFTEN) )

Long version

Due to the expected large range of data values in the final stacked image, saving them as compressed 16-bit integer images with linear BSCALE and BZERO scaling values is likely to offer poor reconstructions of the stacked image. This will lead either to truncation of the extrema of the image, or quantized values that are poorly spaced for the image histogram. Saving the images as 32-bit floating point values would alleviate this quantization issue, at the cost of a large increase in the disk space required for the stacked images.

Transforming the data prior to writing to disk by taking the logarithm of the pixel values can resolve this, with the complication that all data values must first be made positive, which then sets the highest quantization sampling near the lowest values in the image. Following techniques used by SDSS, we have instead opted to use the inverse hyperbolic sine function to transform the data. This domain of this function allows any input value to be converted. In addition, the quantization sampling can be tuned by placing the zero of the inverse hyperbolic sine function at the most interesting pixel values.

Formally, prior to being written to disk, the pixel values are transformed by C = alpha * asinh( (L - BOFFSET) / (2.0 * BSOFTEN) ), where L is the linear input pixel values, C the transformed values, alpha = 2.5 log10(e), and BOFFSET = \langle L \rangle and BSOFTEN = sqrt(alpha) * sigma_{L} are saved to the header as the scaling parameters used. The image is then passed to the standard BSCALE and BZERO calculation and saved to disk.

To reverse this process (on subsequent reads of the image, for example in warp-stack difference calculations), the BOFFSET and BSOFTEN parameters are read from the header and the transformation inverted, such that: L = BOFFSET + BSOFTEN * (exp(C / alpha) - exp(-C / alpha) ).