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 }