Differential D25318 Diff 69786 kstars/fitsviewer/stretch.cpp

Changeset View

Standalone View

View Options

kstars/fitsviewer/stretch.cpp

  • This file was added.
1​#include "stretch.h"
2​
3​/* Stretch
4​
5​ This application is free software; you can redistribute it and/or
6​ modify it under the terms of the GNU General Public
7​ License as published by the Free Software Foundation; either
8​ version 2 of the License, or (at your option) any later version.
9​*/
10​
11​#include "stretch.h"
12​
13​#include <math.h>
14​#include <QtConcurrent>
15​
16​namespace {
17​
18​// Returns the median v of the vector.
19​// The vector is modified in an undefined way.
20​template <typename T>
21​T median(std::vector<T>& values)
22​{
23​ const int middle = values.size() / 2;
24​ std::nth_element(values.begin(), values.begin() + middle, values.end());
25​ return values[middle];
26​}
27​
28​// Returns the median of the sample values.
29​// The values are not modified.
30​template <typename T>
31​T median(T *values, int size, int sampleBy)
32​{
33​ const int downsampled_size = size / sampleBy;
34​ std::vector<T> samples(downsampled_size);
35​ for (int index = 0, i = 0; index < size; ++i, index += sampleBy)
36​ samples[i] = values[index];
37​ return median(samples);
38​}
39​
40​// This stretches one channel given the input parameters.
41​// Based on the spec in section 8.5.6
42​// http://pixinsight.com/doc/docs/XISF-1.0-spec/XISF-1.0-spec.html
43​// Uses multiple threads, blocks until done.
44​// The extension parameters are not used.
45​template <typename T>
46​void stretchOneChannel(T *input_buffer, QImage *output_image,
47​ const StretchParams& stretch_params,
48​ int input_range, int image_height, int image_width)
49​{
50​ QVector<QFuture<void>> futures;
51​
52​ // We're outputting uint8, so the max output is 255.
53​ constexpr int maxOutput = 255;
54​
55​ // Maximum possible input value (e.g. 1024*64 - 1 for a 16 bit unsigned int).
56​ const float maxInput = input_range > 1 ? input_range - 1 : input_range;
57​
58​ const float midtones = stretch_params.grey_red.midtones;
59​ const float highlights = stretch_params.grey_red.highlights;
60​ const float shadows = stretch_params.grey_red.shadows;
61​
62​ // Precomputed expressions moved out of the loop.
63​ // hightlights - shadows, protecting for divide-by-0, in a 0->1.0 scale.
64​ const float hsRangeFactor = highlights == shadows ? 1.0 : 1.0 / (highlights - shadows);
65​ // Shadow and highlight values translated to the ADU scale.
66​ const T nativeShadows = shadows * maxInput;
67​ const T nativeHighlights = highlights * maxInput;
68​ // Constants based on above needed for the stretch calculations.
69​ const float k1 = (midtones - 1) * hsRangeFactor * maxOutput / maxInput;
70​ const float k2 = ((2 * midtones) - 1) * hsRangeFactor / maxInput;
71​
72​ for (int j = 0; j < image_height; j++)
73​ {
74​ futures.append(QtConcurrent::run([ = ]()
75​ {
76​ T * inputLine = input_buffer + j * image_width;
77​ auto * scanLine = output_image->scanLine(j);
78​
79​ for (int i = 0; i < image_width; i++)
80​ {
81​ const T input = inputLine[i];
82​ if (input < nativeShadows) scanLine[i] = 0;
83​ else if (input >= nativeHighlights) scanLine[i] = maxOutput;
84​ else
85​ {
86​ const T inputFloored = (input - nativeShadows);
87​ scanLine[i] = (inputFloored * k1) / (inputFloored * k2 - midtones);
88​ }
89​ }
90​ }));
91​ }
92​ for(QFuture<void> future : futures)
93​ future.waitForFinished();
94​}
95​
96​// This is like the above 1-channel stretch, but extended for 3 channels.
97​// This could have been more modular, but the three channels are combined
98​// into a single qRgb value at the end, so it seems the simplest thing is to
99​// replicate the code. It is assume the colors are not interleaved--the red image
100​// is stored fully, then the green, then the blue.
101​template <typename T>
102​void stretchThreeChannels(T *inputBuffer, QImage *outputImage,
103​ const StretchParams& stretchParams,
104​ int inputRange, int imageHeight, int imageWidth)
105​{
106​ QVector<QFuture<void>> futures;
107​
108​ // We're outputting uint8, so the max output is 255.
109​ constexpr int maxOutput = 255;
110​
111​ // Maximum possible input value (e.g. 1024*64 - 1 for a 16 bit unsigned int).
112​ const float maxInput = inputRange > 1 ? inputRange - 1 : inputRange;
113​
114​ const float midtonesR = stretchParams.grey_red.midtones;
115​ const float highlightsR = stretchParams.grey_red.highlights;
116​ const float shadowsR = stretchParams.grey_red.shadows;
117​ const float midtonesG = stretchParams.green.midtones;
118​ const float highlightsG = stretchParams.green.highlights;
119​ const float shadowsG = stretchParams.green.shadows;
120​ const float midtonesB = stretchParams.blue.midtones;
121​ const float highlightsB = stretchParams.blue.highlights;
122​ const float shadowsB = stretchParams.blue.shadows;
123​
124​ // Precomputed expressions moved out of the loop.
125​ // hightlights - shadows, protecting for divide-by-0, in a 0->1.0 scale.
126​ const float hsRangeFactorR = highlightsR == shadowsR ? 1.0 : 1.0 / (highlightsR - shadowsR);
127​ const float hsRangeFactorG = highlightsG == shadowsG ? 1.0 : 1.0 / (highlightsG - shadowsG);
128​ const float hsRangeFactorB = highlightsB == shadowsB ? 1.0 : 1.0 / (highlightsB - shadowsB);
129​ // Shadow and highlight values translated to the ADU scale.
130​ const T nativeShadowsR = shadowsR * maxInput;
131​ const T nativeShadowsG = shadowsG * maxInput;
132​ const T nativeShadowsB = shadowsB * maxInput;
133​ const T nativeHighlightsR = highlightsR * maxInput;
134​ const T nativeHighlightsG = highlightsG * maxInput;
135​ const T nativeHighlightsB = highlightsB * maxInput;
136​ // Constants based on above needed for the stretch calculations.
137​ const float k1R = (midtonesR - 1) * hsRangeFactorR * maxOutput / maxInput;
138​ const float k1G = (midtonesG - 1) * hsRangeFactorG * maxOutput / maxInput;
139​ const float k1B = (midtonesB - 1) * hsRangeFactorB * maxOutput / maxInput;
140​ const float k2R = ((2 * midtonesR) - 1) * hsRangeFactorR / maxInput;
141​ const float k2G = ((2 * midtonesG) - 1) * hsRangeFactorG / maxInput;
142​ const float k2B = ((2 * midtonesB) - 1) * hsRangeFactorB / maxInput;
143​
144​ const int size = imageWidth * imageHeight;
145​
146​ for (int j = 0; j < imageHeight; j++)
147​ {
148​ futures.append(QtConcurrent::run([ = ]()
149​ {
150​ // R, G, B input images are stored one after another.
151​ T * inputLineR = inputBuffer + j * imageWidth;
152​ T * inputLineG = inputLineR + size;
153​ T * inputLineB = inputLineG + size;
154​
155​ auto * scanLine = reinterpret_cast<QRgb*>(outputImage->scanLine(j));
156​
157​ for (int i = 0; i < imageWidth; i++)
158​ {
159​ const T inputR = inputLineR[i];
160​ const T inputG = inputLineG[i];
161​ const T inputB = inputLineB[i];
162​
163​ uint8_t red, green, blue;
164​
165​ if (inputR < nativeShadowsR) red = 0;
166​ else if (inputR >= nativeHighlightsR) red = maxOutput;
167​ else
168​ {
169​ const T inputFloored = (inputR - nativeShadowsR);
170​ red = (inputFloored * k1R) / (inputFloored * k2R - midtonesR);
171​ }
172​
173​ if (inputG < nativeShadowsG) green = 0;
174​ else if (inputG >= nativeHighlightsG) green = maxOutput;
175​ else
176​ {
177​ const T inputFloored = (inputG - nativeShadowsG);
178​ green = (inputFloored * k1G) / (inputFloored * k2G - midtonesG);
179​ }
180​
181​ if (inputB < nativeShadowsB) blue = 0;
182​ else if (inputB >= nativeHighlightsB) blue = maxOutput;
183​ else
184​ {
185​ const T inputFloored = (inputB - nativeShadowsB);
186​ blue = (inputFloored * k1B) / (inputFloored * k2B - midtonesB);
187​ }
188​ scanLine[i] = qRgb(red, green, blue);
189​ }
190​ }));
191​ }
192​ for(QFuture<void> future : futures)
193​ future.waitForFinished();
194​}
195​
196​// See section 8.5.7 in above link http://pixinsight.com/doc/docs/XISF-1.0-spec/XISF-1.0-spec.html
197​template <typename T>
198​void computeParamsOneChannel(T *buffer, StretchParams1Channel *params,
199​ int inputRange, int height, int width)
200​{
201​ // Find the median sample.
202​ constexpr int maxSamples = 500000;
203​ const int sampleBy = width * height < maxSamples ? 1 : width * height / maxSamples;
204​ const int size = width * height;
205​ T medianSample = median(buffer, width * height, sampleBy);
206​
207​ // Find the Median deviation: 1.4826 * median of abs(sample[i] - median).
208​ const int numSamples = width * height / sampleBy;
209​ std::vector<T> deviations(numSamples);
210​ for (int index = 0, i = 0; index < size; ++i, index += sampleBy)
211​ {
212​ if (medianSample > buffer[index])
213​ deviations[i] = medianSample - buffer[index];
214​ else
215​ deviations[i] = buffer[index] - medianSample;
216​ }
217​
218​ // Shift everything to 0 -> 1.0.
219​ const float medDev = median(deviations);
220​ const float normalizedMedian = medianSample / static_cast<float>(inputRange);
221​ const float MADN = 1.4826 * medDev / static_cast<float>(inputRange);
222​
223​ const bool upperHalf = normalizedMedian > 0.5;
224​
225​ const float shadows = (upperHalf || MADN == 0) ? 0.0 :
226​ fmin(1.0, fmax(0.0, (normalizedMedian + -2.8 * MADN)));
227​
228​ const float highlights = (!upperHalf || MADN == 0) ? 1.0 :
229​ fmin(1.0, fmax(0.0, (normalizedMedian - -2.8 * MADN)));
230​
231​ float X, M;
232​ constexpr float B = 0.25;
233​ if (!upperHalf) {
234​ X = normalizedMedian - shadows;
235​ M = B;
236​ } else {
237​ X = B;
238​ M = highlights - normalizedMedian;
239​ }
240​ float midtones;
241​ if (X == 0) midtones = 0;
242​ else if (X == M) midtones = 0.5;
243​ else if (X == 1) midtones = 1.0;
244​ else midtones = ((M - 1) * X) / ((2 * M - 1) * X - M);
245​
246​ // Store the params.
247​ params->shadows = shadows;
248​ params->highlights = highlights;
249​ params->midtones = midtones;
250​ params->shadows_expansion = 0.0;
251​ params->highlights_expansion = 1.0;
252​}
253​
254​// Need to know the possible range of input values.
255​// Using the type of the sample and guessing.
256​// Perhaps we should examine the contents for the file
257​// (e.g. look at maximum value and extrapolate from that).
258​int getRange(uint8_t sample) { return 256; }
259​int getRange(char sample) { return 256; }
260​int getRange(unsigned short sample) { return 64 * 1024; }
261​int getRange(short sample) { return 64 * 1024; }
262​int getRange(int sample) { return 64 * 1024; }
263​int getRange(unsigned int sample) { return 64 * 1024; }
264​int getRange(long sample) { return 64 * 1024; }
265​int getRange(long long sample) { return 64 * 1024; }
266​int getRange(float sample) { return 64 * 1024; }
267​int getRange(double sample) { return 64 * 1024; }
268​
269​} // namespace
270​
271​template <typename T>
272​Stretch<T>::Stretch(T * image_buffer, int width, int height, int channels)
273​{
274​ buffer = image_buffer;
275​ image_width = width;
276​ image_height = height;
277​ image_channels = channels;
278​
279​ T sample;
280​ input_range = getRange(sample);
281​}
282​
283​template <typename T>
284​void Stretch<T>::run(QImage *outputImage)
285​{
286​ if (image_channels == 1)
287​ stretchOneChannel(buffer, outputImage, params, input_range, image_height, image_width);
288​ else if (image_channels == 3)
289​ stretchThreeChannels(buffer, outputImage, params, input_range, image_height, image_width);
290​}
291​
292​template <typename T>
293​StretchParams Stretch<T>::computeParams()
294​{
295​ StretchParams result;
296​ for (int channel = 0; channel < image_channels; ++channel)
297​ {
298​ int offset = channel * image_width * image_height;
299​ StretchParams1Channel *params = channel == 0 ? &result.grey_red :
300​ (channel == 1 ? &result.green : &result.blue);
301​ computeParamsOneChannel(buffer + offset, params, input_range, image_height, image_width);
302​ }
303​ return result;
304​}
305​
306​// These declarations force the compiler to instantiate these templates,
307​// allowing us to place the template definitions in the .cpp instead of the .h file.
308​template Stretch<unsigned char>::Stretch(unsigned char *, int , int, int);
309​template void Stretch<unsigned char>::run(QImage *);
310​template StretchParams Stretch<unsigned char>::computeParams();
311​
312​template Stretch<short>::Stretch(short *, int , int, int);
313​template void Stretch<short>::run(QImage *);
314​template StretchParams Stretch<short>::computeParams();
315​
316​template Stretch<unsigned short>::Stretch(unsigned short *, int , int, int);
317​template void Stretch<unsigned short>::run(QImage *);
318​template StretchParams Stretch<unsigned short>::computeParams();
319​
320​template Stretch<int>::Stretch(int *, int , int, int);
321​template void Stretch<int>::run(QImage *);
322​template StretchParams Stretch<int>::computeParams();
323​
324​template Stretch<unsigned int>::Stretch(unsigned int *, int , int, int);
325​template void Stretch<unsigned int>::run(QImage *);
326​template StretchParams Stretch<unsigned int>::computeParams();
327​
328​template Stretch<float>::Stretch(float *, int , int, int);
329​template void Stretch<float>::run(QImage *);
330​template StretchParams Stretch<float>::computeParams();
331​
332​template Stretch<long long>::Stretch(long long *, int , int, int);
333​template void Stretch<long long>::run(QImage *);
334​template StretchParams Stretch<long long>::computeParams();
335​
336​template Stretch<double>::Stretch(double *, int , int, int);
337​template void Stretch<double>::run(QImage *);
338​template StretchParams Stretch<double>::computeParams();