File indexing completed on 2024-03-24 15:17:03

0001 /*  Stretch
0002 
0003     SPDX-License-Identifier: GPL-2.0-or-later
0004 */
0005 
0006 #include "stretch.h"
0007 
0008 #include <fitsio.h>
0009 #include <math.h>
0010 #include <QtConcurrent>
0011 
0012 namespace
0013 {
0014 
0015 // Returns the median value of the vector.
0016 // The vector is modified in an undefined way.
0017 template <typename T>
0018 T median(std::vector<T> &values)
0019 {
0020     const int middle = values.size() / 2;
0021     std::nth_element(values.begin(), values.begin() + middle, values.end());
0022     return values[middle];
0023 }
0024 
0025 // Returns the rough max of the buffer.
0026 template <typename T>
0027 T sampledMax(T const *values, int size, int sampleBy)
0028 {
0029     T maxVal = 0;
0030     for (int i = 0; i < size; i += sampleBy)
0031         if (maxVal < values[i])
0032             maxVal = values[i];
0033     return  maxVal;
0034 }
0035 
0036 // Returns the median of the sample values.
0037 // The values are not modified.
0038 template <typename T>
0039 T median(T const *values, int size, int sampleBy)
0040 {
0041     const int downsampled_size = size / sampleBy;
0042     std::vector<T> samples(downsampled_size);
0043     for (int index = 0, i = 0; i < downsampled_size; ++i, index += sampleBy)
0044         samples[i] = values[index];
0045     return median(samples);
0046 }
0047 
0048 // This stretches one channel given the input parameters.
0049 // Based on the spec in section 8.5.6
0050 // https://pixinsight.com/doc/docs/XISF-1.0-spec/XISF-1.0-spec.html
0051 // Uses multiple threads, blocks until done.
0052 // The extension parameters are not used.
0053 // Sampling is applied to the output (that is, with sampling=2, we compute every other output
0054 // sample both in width and height, so the output would have about 4X fewer pixels.
0055 template <typename T>
0056 void stretchOneChannel(T *input_buffer, QImage *output_image,
0057                        const StretchParams &stretch_params,
0058                        int input_range, int image_height, int image_width, int sampling)
0059 {
0060     QVector<QFuture<void>> futures;
0061 
0062     // We're outputting uint8, so the max output is 255.
0063     constexpr int maxOutput = 255;
0064 
0065     // Maximum possible input value (e.g. 1024*64 - 1 for a 16 bit unsigned int).
0066     const float maxInput = input_range > 1 ? input_range - 1 : input_range;
0067 
0068     const float midtones   = stretch_params.grey_red.midtones;
0069     const float highlights = stretch_params.grey_red.highlights;
0070     const float shadows    = stretch_params.grey_red.shadows;
0071 
0072     // Precomputed expressions moved out of the loop.
0073     // highlights - shadows, protecting for divide-by-0, in a 0->1.0 scale.
0074     const float hsRangeFactor = highlights == shadows ? 1.0f : 1.0f / (highlights - shadows);
0075     // Shadow and highlight values translated to the ADU scale.
0076     const T nativeShadows = shadows * maxInput;
0077     const T nativeHighlights = highlights * maxInput;
0078     // Constants based on above needed for the stretch calculations.
0079     const float k1 = (midtones - 1) * hsRangeFactor * maxOutput / maxInput;
0080     const float k2 = ((2 * midtones) - 1) * hsRangeFactor / maxInput;
0081 
0082     // Increment the input index by the sampling, the output index increments by 1.
0083     for (int j = 0, jout = 0; j < image_height; j += sampling, jout++)
0084     {
0085         futures.append(QtConcurrent::run([ = ]()
0086         {
0087             T * inputLine  = input_buffer + j * image_width;
0088             auto * scanLine = output_image->scanLine(jout);
0089 
0090             for (int i = 0, iout = 0; i < image_width; i += sampling, iout++)
0091             {
0092                 const T input = inputLine[i];
0093                 if (input < nativeShadows) scanLine[iout] = 0;
0094                 else if (input >= nativeHighlights) scanLine[iout] = maxOutput;
0095                 else
0096                 {
0097                     const T inputFloored = (input - nativeShadows);
0098                     scanLine[iout] = (inputFloored * k1) / (inputFloored * k2 - midtones);
0099                 }
0100             }
0101         }));
0102     }
0103     for(QFuture<void> future : futures)
0104         future.waitForFinished();
0105 }
0106 
0107 // This is like the above 1-channel stretch, but extended for 3 channels.
0108 // This could have been more modular, but the three channels are combined
0109 // into a single qRgb value at the end, so it seems the simplest thing is to
0110 // replicate the code. It is assume the colors are not interleaved--the red image
0111 // is stored fully, then the green, then the blue.
0112 // Sampling is applied to the output (that is, with sampling=2, we compute every other output
0113 // sample both in width and height, so the output would have about 4X fewer pixels.
0114 template <typename T>
0115 void stretchThreeChannels(T *inputBuffer, QImage *outputImage,
0116                           const StretchParams &stretchParams,
0117                           int inputRange, int imageHeight, int imageWidth, int sampling)
0118 {
0119     QVector<QFuture<void>> futures;
0120 
0121     // We're outputting uint8, so the max output is 255.
0122     constexpr int maxOutput = 255;
0123 
0124     // Maximum possible input value (e.g. 1024*64 - 1 for a 16 bit unsigned int).
0125     const float maxInput = inputRange > 1 ? inputRange - 1 : inputRange;
0126 
0127     const float midtonesR   = stretchParams.grey_red.midtones;
0128     const float highlightsR = stretchParams.grey_red.highlights;
0129     const float shadowsR    = stretchParams.grey_red.shadows;
0130     const float midtonesG   = stretchParams.green.midtones;
0131     const float highlightsG = stretchParams.green.highlights;
0132     const float shadowsG    = stretchParams.green.shadows;
0133     const float midtonesB   = stretchParams.blue.midtones;
0134     const float highlightsB = stretchParams.blue.highlights;
0135     const float shadowsB    = stretchParams.blue.shadows;
0136 
0137     // Precomputed expressions moved out of the loop.
0138     // highlights - shadows, protecting for divide-by-0, in a 0->1.0 scale.
0139     const float hsRangeFactorR = highlightsR == shadowsR ? 1.0f : 1.0f / (highlightsR - shadowsR);
0140     const float hsRangeFactorG = highlightsG == shadowsG ? 1.0f : 1.0f / (highlightsG - shadowsG);
0141     const float hsRangeFactorB = highlightsB == shadowsB ? 1.0f : 1.0f / (highlightsB - shadowsB);
0142     // Shadow and highlight values translated to the ADU scale.
0143     const T nativeShadowsR = shadowsR * maxInput;
0144     const T nativeShadowsG = shadowsG * maxInput;
0145     const T nativeShadowsB = shadowsB * maxInput;
0146     const T nativeHighlightsR = highlightsR * maxInput;
0147     const T nativeHighlightsG = highlightsG * maxInput;
0148     const T nativeHighlightsB = highlightsB * maxInput;
0149     // Constants based on above needed for the stretch calculations.
0150     const float k1R = (midtonesR - 1) * hsRangeFactorR * maxOutput / maxInput;
0151     const float k1G = (midtonesG - 1) * hsRangeFactorG * maxOutput / maxInput;
0152     const float k1B = (midtonesB - 1) * hsRangeFactorB * maxOutput / maxInput;
0153     const float k2R = ((2 * midtonesR) - 1) * hsRangeFactorR / maxInput;
0154     const float k2G = ((2 * midtonesG) - 1) * hsRangeFactorG / maxInput;
0155     const float k2B = ((2 * midtonesB) - 1) * hsRangeFactorB / maxInput;
0156 
0157     const int size = imageWidth * imageHeight;
0158 
0159     for (int j = 0, jout = 0; j < imageHeight; j += sampling, jout++)
0160     {
0161         futures.append(QtConcurrent::run([ = ]()
0162         {
0163             // R, G, B input images are stored one after another.
0164             T * inputLineR  = inputBuffer + j * imageWidth;
0165             T * inputLineG  = inputLineR + size;
0166             T * inputLineB  = inputLineG + size;
0167 
0168             auto * scanLine = reinterpret_cast<QRgb*>(outputImage->scanLine(jout));
0169 
0170             for (int i = 0, iout = 0; i < imageWidth; i += sampling, iout++)
0171             {
0172                 const T inputR = inputLineR[i];
0173                 const T inputG = inputLineG[i];
0174                 const T inputB = inputLineB[i];
0175 
0176                 uint8_t red, green, blue;
0177 
0178                 if (inputR < nativeShadowsR) red = 0;
0179                 else if (inputR >= nativeHighlightsR) red = maxOutput;
0180                 else
0181                 {
0182                     const T inputFloored = (inputR - nativeShadowsR);
0183                     red = (inputFloored * k1R) / (inputFloored * k2R - midtonesR);
0184                 }
0185 
0186                 if (inputG < nativeShadowsG) green = 0;
0187                 else if (inputG >= nativeHighlightsG) green = maxOutput;
0188                 else
0189                 {
0190                     const T inputFloored = (inputG - nativeShadowsG);
0191                     green = (inputFloored * k1G) / (inputFloored * k2G - midtonesG);
0192                 }
0193 
0194                 if (inputB < nativeShadowsB) blue = 0;
0195                 else if (inputB >= nativeHighlightsB) blue = maxOutput;
0196                 else
0197                 {
0198                     const T inputFloored = (inputB - nativeShadowsB);
0199                     blue = (inputFloored * k1B) / (inputFloored * k2B - midtonesB);
0200                 }
0201                 scanLine[iout] = qRgb(red, green, blue);
0202             }
0203         }));
0204     }
0205     for(QFuture<void> future : futures)
0206         future.waitForFinished();
0207 }
0208 
0209 template <typename T>
0210 void stretchChannels(T *input_buffer, QImage *output_image,
0211                      const StretchParams &stretch_params,
0212                      int input_range, int image_height, int image_width, int num_channels, int sampling)
0213 {
0214     if (num_channels == 1)
0215         stretchOneChannel(input_buffer, output_image, stretch_params, input_range,
0216                           image_height, image_width, sampling);
0217     else if (num_channels == 3)
0218         stretchThreeChannels(input_buffer, output_image, stretch_params, input_range,
0219                              image_height, image_width, sampling);
0220 }
0221 
0222 // See section 8.5.7 in above link  https://pixinsight.com/doc/docs/XISF-1.0-spec/XISF-1.0-spec.html
0223 template <typename T>
0224 void computeParamsOneChannel(T const *buffer, StretchParams1Channel *params,
0225                              int inputRange, int height, int width)
0226 {
0227     // Find the median sample.
0228     constexpr int maxSamples = 500000;
0229     const int sampleBy = width * height < maxSamples ? 1 : width * height / maxSamples;
0230 
0231     T medianSample = median(buffer, width * height, sampleBy);
0232     // Find the Median deviation: 1.4826 * median of abs(sample[i] - median).
0233     const int numSamples = width * height / sampleBy;
0234     std::vector<T> deviations(numSamples);
0235     for (int index = 0, i = 0; i < numSamples; ++i, index += sampleBy)
0236     {
0237         if (medianSample > buffer[index])
0238             deviations[i] = medianSample - buffer[index];
0239         else
0240             deviations[i] = buffer[index] - medianSample;
0241     }
0242 
0243     // Shift everything to 0 -> 1.0.
0244     const float medDev = median(deviations);
0245     const float normalizedMedian = medianSample / static_cast<float>(inputRange);
0246     const float MADN = 1.4826 * medDev / static_cast<float>(inputRange);
0247 
0248     const bool upperHalf = normalizedMedian > 0.5;
0249 
0250     const float shadows = (upperHalf || MADN == 0) ? 0.0 :
0251                           fmin(1.0, fmax(0.0, (normalizedMedian + -2.8 * MADN)));
0252 
0253     const float highlights = (!upperHalf || MADN == 0) ? 1.0 :
0254                              fmin(1.0, fmax(0.0, (normalizedMedian - -2.8 * MADN)));
0255 
0256     float X, M;
0257     constexpr float B = 0.25;
0258     if (!upperHalf)
0259     {
0260         X = normalizedMedian - shadows;
0261         M = B;
0262     }
0263     else
0264     {
0265         X = B;
0266         M = highlights - normalizedMedian;
0267     }
0268     float midtones;
0269     if (X == 0) midtones = 0.0f;
0270     else if (X == M) midtones = 0.5f;
0271     else if (X == 1) midtones = 1.0f;
0272     else midtones = ((M - 1) * X) / ((2 * M - 1) * X - M);
0273 
0274     // Store the params.
0275     params->shadows = shadows;
0276     params->highlights = highlights;
0277     params->midtones = midtones;
0278     params->shadows_expansion = 0.0;
0279     params->highlights_expansion = 1.0;
0280 }
0281 
0282 // Need to know the possible range of input values.
0283 // Using the type of the sample and guessing.
0284 // Perhaps we should examine the contents for the file
0285 // (e.g. look at maximum value and extrapolate from that).
0286 int getRange(int data_type)
0287 {
0288     switch (data_type)
0289     {
0290         case TBYTE:
0291             return 256;
0292         case TSHORT:
0293             return 64 * 1024;
0294         case TUSHORT:
0295             return 64 * 1024;
0296         case TLONG:
0297             return 64 * 1024;
0298         case TFLOAT:
0299             return 64 * 1024;
0300         case TLONGLONG:
0301             return 64 * 1024;
0302         case TDOUBLE:
0303             return 64 * 1024;
0304         default:
0305             return 64 * 1024;
0306     }
0307 }
0308 
0309 }  // namespace
0310 
0311 Stretch::Stretch(int width, int height, int channels, int data_type)
0312 {
0313     image_width = width;
0314     image_height = height;
0315     image_channels = channels;
0316     dataType = data_type;
0317     input_range = getRange(dataType);
0318 }
0319 
0320 void Stretch::run(uint8_t const *input, QImage *outputImage, int sampling)
0321 {
0322     Q_ASSERT(outputImage->width() == (image_width + sampling - 1) / sampling);
0323     Q_ASSERT(outputImage->height() == (image_height + sampling - 1) / sampling);
0324     recalculateInputRange(input);
0325 
0326     switch (dataType)
0327     {
0328         case TBYTE:
0329             stretchChannels(reinterpret_cast<uint8_t const*>(input), outputImage, params,
0330                             input_range, image_height, image_width, image_channels, sampling);
0331             break;
0332         case TSHORT:
0333             stretchChannels(reinterpret_cast<short const*>(input), outputImage, params,
0334                             input_range, image_height, image_width, image_channels, sampling);
0335             break;
0336         case TUSHORT:
0337             stretchChannels(reinterpret_cast<unsigned short const*>(input), outputImage, params,
0338                             input_range, image_height, image_width, image_channels, sampling);
0339             break;
0340         case TLONG:
0341             stretchChannels(reinterpret_cast<long const*>(input), outputImage, params,
0342                             input_range, image_height, image_width, image_channels, sampling);
0343             break;
0344         case TFLOAT:
0345             stretchChannels(reinterpret_cast<float const*>(input), outputImage, params,
0346                             input_range, image_height, image_width, image_channels, sampling);
0347             break;
0348         case TLONGLONG:
0349             stretchChannels(reinterpret_cast<long long const*>(input), outputImage, params,
0350                             input_range, image_height, image_width, image_channels, sampling);
0351             break;
0352         case TDOUBLE:
0353             stretchChannels(reinterpret_cast<double const*>(input), outputImage, params,
0354                             input_range, image_height, image_width, image_channels, sampling);
0355             break;
0356         default:
0357             break;
0358     }
0359 }
0360 
0361 // The input range for float/double is ambiguous, and we can't tell without the buffer,
0362 // so we set it to 64K and possibly reduce it when we see the data.
0363 void Stretch::recalculateInputRange(uint8_t const *input)
0364 {
0365     if (input_range <= 1) return;
0366     if (dataType != TFLOAT && dataType != TDOUBLE) return;
0367 
0368     float mx = 0;
0369     if (dataType == TFLOAT)
0370         mx = sampledMax(reinterpret_cast<float const*>(input), image_height * image_width, 1000);
0371     else if (dataType == TDOUBLE)
0372         mx = sampledMax(reinterpret_cast<double const*>(input), image_height * image_width, 1000);
0373     if (mx <= 1.01f) input_range = 1;
0374 }
0375 
0376 StretchParams Stretch::computeParams(uint8_t const *input)
0377 {
0378     recalculateInputRange(input);
0379     StretchParams result;
0380     for (int channel = 0; channel < image_channels; ++channel)
0381     {
0382         int offset = channel * image_width * image_height;
0383         StretchParams1Channel *params = channel == 0 ? &result.grey_red :
0384                                         (channel == 1 ? &result.green : &result.blue);
0385         switch (dataType)
0386         {
0387             case TBYTE:
0388             {
0389                 auto buffer = reinterpret_cast<uint8_t const*>(input);
0390                 computeParamsOneChannel(buffer + offset, params, input_range,
0391                                         image_height, image_width);
0392                 break;
0393             }
0394             case TSHORT:
0395             {
0396                 auto buffer = reinterpret_cast<short const*>(input);
0397                 computeParamsOneChannel(buffer + offset, params, input_range,
0398                                         image_height, image_width);
0399                 break;
0400             }
0401             case TUSHORT:
0402             {
0403                 auto buffer = reinterpret_cast<unsigned short const*>(input);
0404                 computeParamsOneChannel(buffer + offset, params, input_range,
0405                                         image_height, image_width);
0406                 break;
0407             }
0408             case TLONG:
0409             {
0410                 auto buffer = reinterpret_cast<long const*>(input);
0411                 computeParamsOneChannel(buffer + offset, params, input_range,
0412                                         image_height, image_width);
0413                 break;
0414             }
0415             case TFLOAT:
0416             {
0417                 auto buffer = reinterpret_cast<float const*>(input);
0418                 computeParamsOneChannel(buffer + offset, params, input_range,
0419                                         image_height, image_width);
0420                 break;
0421             }
0422             case TLONGLONG:
0423             {
0424                 auto buffer = reinterpret_cast<long long const*>(input);
0425                 computeParamsOneChannel(buffer + offset, params, input_range,
0426                                         image_height, image_width);
0427                 break;
0428             }
0429             case TDOUBLE:
0430             {
0431                 auto buffer = reinterpret_cast<double const*>(input);
0432                 computeParamsOneChannel(buffer + offset, params, input_range,
0433                                         image_height, image_width);
0434                 break;
0435             }
0436             default:
0437                 break;
0438         }
0439     }
0440     return result;
0441 }