blob: 4f21ac637246128218cdc891b8597b01673a6ec8 [file] [log] [blame]
Nick Deakinf6bca5a2022-11-04 10:43:43 -04001/*
2 * Copyright 2022 The Android Open Source Project
3 *
4 * Licensed under the Apache License, Version 2.0 (the "License");
5 * you may not use this file except in compliance with the License.
6 * You may obtain a copy of the License at
7 *
8 * http://www.apache.org/licenses/LICENSE-2.0
9 *
10 * Unless required by applicable law or agreed to in writing, software
11 * distributed under the License is distributed on an "AS IS" BASIS,
12 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13 * See the License for the specific language governing permissions and
14 * limitations under the License.
15 */
16
17#include <cmath>
Harish Mahendrakar555a06b2022-12-14 09:37:27 -080018#include <vector>
Nick Deakinf6bca5a2022-11-04 10:43:43 -040019#include <jpegrecoverymap/recoverymapmath.h>
20
21namespace android::recoverymap {
22
Harish Mahendrakar555a06b2022-12-14 09:37:27 -080023constexpr size_t kPqOETFPrecision = 10;
24constexpr size_t kPqOETFNumEntries = 1 << kPqOETFPrecision;
25
26static const std::vector<float> kPqOETF = [] {
27 std::vector<float> result;
28 float increment = 1.0 / kPqOETFNumEntries;
29 float value = 0.0f;
30 for (int idx = 0; idx < kPqOETFNumEntries; idx++, value += increment) {
31 result.push_back(pqOetf(value));
32 }
33 return result;
34}();
35
36constexpr size_t kPqInvOETFPrecision = 10;
37constexpr size_t kPqInvOETFNumEntries = 1 << kPqInvOETFPrecision;
38
39static const std::vector<float> kPqInvOETF = [] {
40 std::vector<float> result;
41 float increment = 1.0 / kPqInvOETFNumEntries;
42 float value = 0.0f;
43 for (int idx = 0; idx < kPqInvOETFNumEntries; idx++, value += increment) {
44 result.push_back(pqInvOetf(value));
45 }
46 return result;
47}();
48
49constexpr size_t kHlgOETFPrecision = 10;
50constexpr size_t kHlgOETFNumEntries = 1 << kHlgOETFPrecision;
51
52static const std::vector<float> kHlgOETF = [] {
53 std::vector<float> result;
54 float increment = 1.0 / kHlgOETFNumEntries;
55 float value = 0.0f;
56 for (int idx = 0; idx < kHlgOETFNumEntries; idx++, value += increment) {
57 result.push_back(hlgOetf(value));
58 }
59 return result;
60}();
61
62constexpr size_t kHlgInvOETFPrecision = 10;
63constexpr size_t kHlgInvOETFNumEntries = 1 << kHlgInvOETFPrecision;
64
65static const std::vector<float> kHlgInvOETF = [] {
66 std::vector<float> result;
67 float increment = 1.0 / kHlgInvOETFNumEntries;
68 float value = 0.0f;
69 for (int idx = 0; idx < kHlgInvOETFNumEntries; idx++, value += increment) {
70 result.push_back(hlgInvOetf(value));
71 }
72 return result;
73}();
74
75constexpr size_t kSRGBInvOETFPrecision = 10;
76constexpr size_t kSRGBInvOETFNumEntries = 1 << kSRGBInvOETFPrecision;
77static const std::vector<float> kSRGBInvOETF = [] {
78 std::vector<float> result;
79 float increment = 1.0 / kSRGBInvOETFNumEntries;
80 float value = 0.0f;
81 for (int idx = 0; idx < kSRGBInvOETFNumEntries; idx++, value += increment) {
82 result.push_back(srgbInvOetf(value));
83 }
84 return result;
85}();
Ram Mohanfe723d62022-12-15 00:59:11 +053086
87// Use Shepard's method for inverse distance weighting. For more information:
88// en.wikipedia.org/wiki/Inverse_distance_weighting#Shepard's_method
89
90float ShepardsIDW::euclideanDistance(float x1, float x2, float y1, float y2) {
91 return sqrt(((y2 - y1) * (y2 - y1)) + (x2 - x1) * (x2 - x1));
92}
93
94void ShepardsIDW::fillShepardsIDW(float *weights, int incR, int incB) {
95 for (int y = 0; y < mMapScaleFactor; y++) {
96 for (int x = 0; x < mMapScaleFactor; x++) {
97 float pos_x = ((float)x) / mMapScaleFactor;
98 float pos_y = ((float)y) / mMapScaleFactor;
99 int curr_x = floor(pos_x);
100 int curr_y = floor(pos_y);
101 int next_x = curr_x + incR;
102 int next_y = curr_y + incB;
103 float e1_distance = euclideanDistance(pos_x, curr_x, pos_y, curr_y);
104 int index = y * mMapScaleFactor * 4 + x * 4;
105 if (e1_distance == 0) {
106 weights[index++] = 1.f;
107 weights[index++] = 0.f;
108 weights[index++] = 0.f;
109 weights[index++] = 0.f;
110 } else {
111 float e1_weight = 1.f / e1_distance;
112
113 float e2_distance = euclideanDistance(pos_x, curr_x, pos_y, next_y);
114 float e2_weight = 1.f / e2_distance;
115
116 float e3_distance = euclideanDistance(pos_x, next_x, pos_y, curr_y);
117 float e3_weight = 1.f / e3_distance;
118
119 float e4_distance = euclideanDistance(pos_x, next_x, pos_y, next_y);
120 float e4_weight = 1.f / e4_distance;
121
122 float total_weight = e1_weight + e2_weight + e3_weight + e4_weight;
123
124 weights[index++] = e1_weight / total_weight;
125 weights[index++] = e2_weight / total_weight;
126 weights[index++] = e3_weight / total_weight;
127 weights[index++] = e4_weight / total_weight;
128 }
129 }
130 }
131}
132
Nick Deakin594a4ca2022-11-16 20:57:42 -0500133////////////////////////////////////////////////////////////////////////////////
134// sRGB transformations
Nick Deakinf6bca5a2022-11-04 10:43:43 -0400135
Harish Mahendrakar1107ff32022-12-07 17:24:35 -0800136static const float kMaxPixelFloat = 1.0f;
137static float clampPixelFloat(float value) {
138 return (value < 0.0f) ? 0.0f : (value > kMaxPixelFloat) ? kMaxPixelFloat : value;
139}
140
Nick Deakin65f492a2022-11-29 22:47:40 -0500141// See IEC 61966-2-1, Equation F.7.
142static const float kSrgbR = 0.2126f, kSrgbG = 0.7152f, kSrgbB = 0.0722f;
Nick Deakin594a4ca2022-11-16 20:57:42 -0500143
144float srgbLuminance(Color e) {
145 return kSrgbR * e.r + kSrgbG * e.g + kSrgbB * e.b;
Nick Deakinf6bca5a2022-11-04 10:43:43 -0400146}
147
Nick Deakin65f492a2022-11-29 22:47:40 -0500148// See ECMA TR/98, Section 7.
Nick Deakinf6bca5a2022-11-04 10:43:43 -0400149static const float kSrgbRCr = 1.402f, kSrgbGCb = 0.34414f, kSrgbGCr = 0.71414f, kSrgbBCb = 1.772f;
150
Nick Deakin594a4ca2022-11-16 20:57:42 -0500151Color srgbYuvToRgb(Color e_gamma) {
Harish Mahendrakar1107ff32022-12-07 17:24:35 -0800152 return {{{ clampPixelFloat(e_gamma.y + kSrgbRCr * e_gamma.v),
153 clampPixelFloat(e_gamma.y - kSrgbGCb * e_gamma.u - kSrgbGCr * e_gamma.v),
154 clampPixelFloat(e_gamma.y + kSrgbBCb * e_gamma.u) }}};
Nick Deakinf6bca5a2022-11-04 10:43:43 -0400155}
156
Nick Deakin65f492a2022-11-29 22:47:40 -0500157// See ECMA TR/98, Section 7.
158static const float kSrgbYR = 0.299f, kSrgbYG = 0.587f, kSrgbYB = 0.114f;
Nick Deakin594a4ca2022-11-16 20:57:42 -0500159static const float kSrgbUR = -0.1687f, kSrgbUG = -0.3313f, kSrgbUB = 0.5f;
160static const float kSrgbVR = 0.5f, kSrgbVG = -0.4187f, kSrgbVB = -0.0813f;
161
162Color srgbRgbToYuv(Color e_gamma) {
Nick Deakin65f492a2022-11-29 22:47:40 -0500163 return {{{ kSrgbYR * e_gamma.r + kSrgbYG * e_gamma.g + kSrgbYB * e_gamma.b,
Nick Deakin594a4ca2022-11-16 20:57:42 -0500164 kSrgbUR * e_gamma.r + kSrgbUG * e_gamma.g + kSrgbUB * e_gamma.b,
165 kSrgbVR * e_gamma.r + kSrgbVG * e_gamma.g + kSrgbVB * e_gamma.b }}};
166}
167
Nick Deakin65f492a2022-11-29 22:47:40 -0500168// See IEC 61966-2-1, Equations F.5 and F.6.
Nick Deakin594a4ca2022-11-16 20:57:42 -0500169float srgbInvOetf(float e_gamma) {
170 if (e_gamma <= 0.04045f) {
171 return e_gamma / 12.92f;
Nick Deakinf6bca5a2022-11-04 10:43:43 -0400172 } else {
Nick Deakin594a4ca2022-11-16 20:57:42 -0500173 return pow((e_gamma + 0.055f) / 1.055f, 2.4);
Nick Deakinf6bca5a2022-11-04 10:43:43 -0400174 }
175}
176
Nick Deakin594a4ca2022-11-16 20:57:42 -0500177Color srgbInvOetf(Color e_gamma) {
178 return {{{ srgbInvOetf(e_gamma.r),
179 srgbInvOetf(e_gamma.g),
180 srgbInvOetf(e_gamma.b) }}};
181}
182
Harish Mahendrakar555a06b2022-12-14 09:37:27 -0800183// See IEC 61966-2-1, Equations F.5 and F.6.
184float srgbInvOetfLUT(float e_gamma) {
185 uint32_t value = static_cast<uint32_t>(e_gamma * kSRGBInvOETFNumEntries);
186 //TODO() : Remove once conversion modules have appropriate clamping in place
187 value = CLIP3(value, 0, kSRGBInvOETFNumEntries - 1);
188 return kSRGBInvOETF[value];
189}
190
191Color srgbInvOetfLUT(Color e_gamma) {
192 return {{{ srgbInvOetfLUT(e_gamma.r),
193 srgbInvOetfLUT(e_gamma.g),
194 srgbInvOetfLUT(e_gamma.b) }}};
195}
Nick Deakin594a4ca2022-11-16 20:57:42 -0500196
197////////////////////////////////////////////////////////////////////////////////
198// Display-P3 transformations
199
Nick Deakin65f492a2022-11-29 22:47:40 -0500200// See SMPTE EG 432-1, Table 7-2.
201static const float kP3R = 0.20949f, kP3G = 0.72160f, kP3B = 0.06891f;
Nick Deakin6bd90432022-11-20 16:26:37 -0500202
203float p3Luminance(Color e) {
204 return kP3R * e.r + kP3G * e.g + kP3B * e.b;
205}
Nick Deakin594a4ca2022-11-16 20:57:42 -0500206
207
208////////////////////////////////////////////////////////////////////////////////
209// BT.2100 transformations - according to ITU-R BT.2100-2
210
Nick Deakin65f492a2022-11-29 22:47:40 -0500211// See ITU-R BT.2100-2, Table 5, HLG Reference OOTF
Nick Deakin594a4ca2022-11-16 20:57:42 -0500212static const float kBt2100R = 0.2627f, kBt2100G = 0.6780f, kBt2100B = 0.0593f;
213
214float bt2100Luminance(Color e) {
215 return kBt2100R * e.r + kBt2100G * e.g + kBt2100B * e.b;
216}
217
Nick Deakin65f492a2022-11-29 22:47:40 -0500218// See ITU-R BT.2100-2, Table 6, Derivation of colour difference signals.
Nick Deakin594a4ca2022-11-16 20:57:42 -0500219static const float kBt2100Cb = 1.8814f, kBt2100Cr = 1.4746f;
220
221Color bt2100RgbToYuv(Color e_gamma) {
222 float y_gamma = bt2100Luminance(e_gamma);
223 return {{{ y_gamma,
224 (e_gamma.b - y_gamma) / kBt2100Cb,
225 (e_gamma.r - y_gamma) / kBt2100Cr }}};
226}
227
Nick Deakin65f492a2022-11-29 22:47:40 -0500228// Derived by inversing bt2100RgbToYuv. The derivation for R and B are pretty
229// straight forward; we just invert the formulas for U and V above. But deriving
230// the formula for G is a bit more complicated:
Nick Deakin594a4ca2022-11-16 20:57:42 -0500231//
232// Start with equation for luminance:
233// Y = kBt2100R * R + kBt2100G * G + kBt2100B * B
234// Solve for G:
235// G = (Y - kBt2100R * R - kBt2100B * B) / kBt2100B
236// Substitute equations for R and B in terms YUV:
237// G = (Y - kBt2100R * (Y + kBt2100Cr * V) - kBt2100B * (Y + kBt2100Cb * U)) / kBt2100B
238// Simplify:
239// G = Y * ((1 - kBt2100R - kBt2100B) / kBt2100G)
240// + U * (kBt2100B * kBt2100Cb / kBt2100G)
241// + V * (kBt2100R * kBt2100Cr / kBt2100G)
242//
243// We then get the following coeficients for calculating G from YUV:
244//
245// Coef for Y = (1 - kBt2100R - kBt2100B) / kBt2100G = 1
246// Coef for U = kBt2100B * kBt2100Cb / kBt2100G = kBt2100GCb = ~0.1645
247// Coef for V = kBt2100R * kBt2100Cr / kBt2100G = kBt2100GCr = ~0.5713
248
249static const float kBt2100GCb = kBt2100B * kBt2100Cb / kBt2100G;
250static const float kBt2100GCr = kBt2100R * kBt2100Cr / kBt2100G;
251
252Color bt2100YuvToRgb(Color e_gamma) {
Harish Mahendrakar1107ff32022-12-07 17:24:35 -0800253 return {{{ clampPixelFloat(e_gamma.y + kBt2100Cr * e_gamma.v),
254 clampPixelFloat(e_gamma.y - kBt2100GCb * e_gamma.u - kBt2100GCr * e_gamma.v),
255 clampPixelFloat(e_gamma.y + kBt2100Cb * e_gamma.u) }}};
Nick Deakinf6bca5a2022-11-04 10:43:43 -0400256}
257
Nick Deakin65f492a2022-11-29 22:47:40 -0500258// See ITU-R BT.2100-2, Table 5, HLG Reference OETF.
Nick Deakinf6bca5a2022-11-04 10:43:43 -0400259static const float kHlgA = 0.17883277f, kHlgB = 0.28466892f, kHlgC = 0.55991073;
260
Nick Deakin65f492a2022-11-29 22:47:40 -0500261float hlgOetf(float e) {
Nick Deakinf6bca5a2022-11-04 10:43:43 -0400262 if (e <= 1.0f/12.0f) {
263 return sqrt(3.0f * e);
264 } else {
265 return kHlgA * log(12.0f * e - kHlgB) + kHlgC;
266 }
267}
268
269Color hlgOetf(Color e) {
270 return {{{ hlgOetf(e.r), hlgOetf(e.g), hlgOetf(e.b) }}};
271}
272
Harish Mahendrakar555a06b2022-12-14 09:37:27 -0800273float hlgOetfLUT(float e) {
274 uint32_t value = static_cast<uint32_t>(e * kHlgOETFNumEntries);
275 //TODO() : Remove once conversion modules have appropriate clamping in place
276 value = CLIP3(value, 0, kHlgOETFNumEntries - 1);
277
278 return kHlgOETF[value];
279}
280
281Color hlgOetfLUT(Color e) {
282 return {{{ hlgOetfLUT(e.r), hlgOetfLUT(e.g), hlgOetfLUT(e.b) }}};
283}
284
Nick Deakin65f492a2022-11-29 22:47:40 -0500285// See ITU-R BT.2100-2, Table 5, HLG Reference EOTF.
286float hlgInvOetf(float e_gamma) {
Nick Deakin594a4ca2022-11-16 20:57:42 -0500287 if (e_gamma <= 0.5f) {
288 return pow(e_gamma, 2.0f) / 3.0f;
289 } else {
290 return (exp((e_gamma - kHlgC) / kHlgA) + kHlgB) / 12.0f;
291 }
292}
293
294Color hlgInvOetf(Color e_gamma) {
295 return {{{ hlgInvOetf(e_gamma.r),
296 hlgInvOetf(e_gamma.g),
297 hlgInvOetf(e_gamma.b) }}};
298}
299
Harish Mahendrakar555a06b2022-12-14 09:37:27 -0800300float hlgInvOetfLUT(float e_gamma) {
301 uint32_t value = static_cast<uint32_t>(e_gamma * kHlgInvOETFNumEntries);
302 //TODO() : Remove once conversion modules have appropriate clamping in place
303 value = CLIP3(value, 0, kHlgInvOETFNumEntries - 1);
304
305 return kHlgInvOETF[value];
306}
307
308Color hlgInvOetfLUT(Color e_gamma) {
309 return {{{ hlgInvOetfLUT(e_gamma.r),
310 hlgInvOetfLUT(e_gamma.g),
311 hlgInvOetfLUT(e_gamma.b) }}};
312}
313
Nick Deakin65f492a2022-11-29 22:47:40 -0500314// See ITU-R BT.2100-2, Table 4, Reference PQ OETF.
Nick Deakin6bd90432022-11-20 16:26:37 -0500315static const float kPqM1 = 2610.0f / 16384.0f, kPqM2 = 2523.0f / 4096.0f * 128.0f;
316static const float kPqC1 = 3424.0f / 4096.0f, kPqC2 = 2413.0f / 4096.0f * 32.0f,
317 kPqC3 = 2392.0f / 4096.0f * 32.0f;
318
Nick Deakin65f492a2022-11-29 22:47:40 -0500319float pqOetf(float e) {
320 if (e <= 0.0f) return 0.0f;
321 return pow((kPqC1 + kPqC2 * pow(e, kPqM1)) / (1 + kPqC3 * pow(e, kPqM1)),
Nick Deakin6bd90432022-11-20 16:26:37 -0500322 kPqM2);
323}
324
325Color pqOetf(Color e) {
326 return {{{ pqOetf(e.r), pqOetf(e.g), pqOetf(e.b) }}};
327}
328
Harish Mahendrakar555a06b2022-12-14 09:37:27 -0800329float pqOetfLUT(float e) {
330 uint32_t value = static_cast<uint32_t>(e * kPqOETFNumEntries);
331 //TODO() : Remove once conversion modules have appropriate clamping in place
332 value = CLIP3(value, 0, kPqOETFNumEntries - 1);
333
334 return kPqOETF[value];
335}
336
337Color pqOetfLUT(Color e) {
338 return {{{ pqOetfLUT(e.r), pqOetfLUT(e.g), pqOetfLUT(e.b) }}};
339}
340
Nick Deakin65f492a2022-11-29 22:47:40 -0500341// Derived from the inverse of the Reference PQ OETF.
342static const float kPqInvA = 128.0f, kPqInvB = 107.0f, kPqInvC = 2413.0f, kPqInvD = 2392.0f,
343 kPqInvE = 6.2773946361f, kPqInvF = 0.0126833f;
344
345float pqInvOetf(float e_gamma) {
346 // This equation blows up if e_gamma is 0.0, and checking on <= 0.0 doesn't
347 // always catch 0.0. So, check on 0.0001, since anything this small will
348 // effectively be crushed to zero anyways.
349 if (e_gamma <= 0.0001f) return 0.0f;
350 return pow((kPqInvA * pow(e_gamma, kPqInvF) - kPqInvB)
351 / (kPqInvC - kPqInvD * pow(e_gamma, kPqInvF)),
352 kPqInvE);
Nick Deakin6bd90432022-11-20 16:26:37 -0500353}
354
355Color pqInvOetf(Color e_gamma) {
356 return {{{ pqInvOetf(e_gamma.r),
357 pqInvOetf(e_gamma.g),
358 pqInvOetf(e_gamma.b) }}};
359}
360
Harish Mahendrakar555a06b2022-12-14 09:37:27 -0800361float pqInvOetfLUT(float e_gamma) {
362 uint32_t value = static_cast<uint32_t>(e_gamma * kPqInvOETFNumEntries);
363 //TODO() : Remove once conversion modules have appropriate clamping in place
364 value = CLIP3(value, 0, kPqInvOETFNumEntries - 1);
365
366 return kPqInvOETF[value];
367}
368
369Color pqInvOetfLUT(Color e_gamma) {
370 return {{{ pqInvOetfLUT(e_gamma.r),
371 pqInvOetfLUT(e_gamma.g),
372 pqInvOetfLUT(e_gamma.b) }}};
373}
374
Nick Deakin594a4ca2022-11-16 20:57:42 -0500375
376////////////////////////////////////////////////////////////////////////////////
377// Color conversions
378
Nick Deakin6bd90432022-11-20 16:26:37 -0500379Color bt709ToP3(Color e) {
380 return {{{ 0.82254f * e.r + 0.17755f * e.g + 0.00006f * e.b,
381 0.03312f * e.r + 0.96684f * e.g + -0.00001f * e.b,
382 0.01706f * e.r + 0.07240f * e.g + 0.91049f * e.b }}};
383}
384
385Color bt709ToBt2100(Color e) {
386 return {{{ 0.62740f * e.r + 0.32930f * e.g + 0.04332f * e.b,
387 0.06904f * e.r + 0.91958f * e.g + 0.01138f * e.b,
388 0.01636f * e.r + 0.08799f * e.g + 0.89555f * e.b }}};
389}
390
391Color p3ToBt709(Color e) {
392 return {{{ 1.22482f * e.r + -0.22490f * e.g + -0.00007f * e.b,
393 -0.04196f * e.r + 1.04199f * e.g + 0.00001f * e.b,
394 -0.01961f * e.r + -0.07865f * e.g + 1.09831f * e.b }}};
395}
396
397Color p3ToBt2100(Color e) {
398 return {{{ 0.75378f * e.r + 0.19862f * e.g + 0.04754f * e.b,
399 0.04576f * e.r + 0.94177f * e.g + 0.01250f * e.b,
400 -0.00121f * e.r + 0.01757f * e.g + 0.98359f * e.b }}};
401}
402
403Color bt2100ToBt709(Color e) {
404 return {{{ 1.66045f * e.r + -0.58764f * e.g + -0.07286f * e.b,
405 -0.12445f * e.r + 1.13282f * e.g + -0.00837f * e.b,
406 -0.01811f * e.r + -0.10057f * e.g + 1.11878f * e.b }}};
407}
408
409Color bt2100ToP3(Color e) {
410 return {{{ 1.34369f * e.r + -0.28223f * e.g + -0.06135f * e.b,
411 -0.06533f * e.r + 1.07580f * e.g + -0.01051f * e.b,
412 0.00283f * e.r + -0.01957f * e.g + 1.01679f * e.b
413 }}};
414}
415
416// TODO: confirm we always want to convert like this before calculating
417// luminance.
418ColorTransformFn getHdrConversionFn(jpegr_color_gamut sdr_gamut, jpegr_color_gamut hdr_gamut) {
Nick Deakin65f492a2022-11-29 22:47:40 -0500419 switch (sdr_gamut) {
Nick Deakin6bd90432022-11-20 16:26:37 -0500420 case JPEGR_COLORGAMUT_BT709:
421 switch (hdr_gamut) {
422 case JPEGR_COLORGAMUT_BT709:
423 return identityConversion;
424 case JPEGR_COLORGAMUT_P3:
425 return p3ToBt709;
426 case JPEGR_COLORGAMUT_BT2100:
427 return bt2100ToBt709;
428 case JPEGR_COLORGAMUT_UNSPECIFIED:
429 return nullptr;
430 }
431 break;
432 case JPEGR_COLORGAMUT_P3:
433 switch (hdr_gamut) {
434 case JPEGR_COLORGAMUT_BT709:
435 return bt709ToP3;
436 case JPEGR_COLORGAMUT_P3:
437 return identityConversion;
438 case JPEGR_COLORGAMUT_BT2100:
439 return bt2100ToP3;
440 case JPEGR_COLORGAMUT_UNSPECIFIED:
441 return nullptr;
442 }
443 break;
444 case JPEGR_COLORGAMUT_BT2100:
445 switch (hdr_gamut) {
446 case JPEGR_COLORGAMUT_BT709:
447 return bt709ToBt2100;
448 case JPEGR_COLORGAMUT_P3:
449 return p3ToBt2100;
450 case JPEGR_COLORGAMUT_BT2100:
451 return identityConversion;
452 case JPEGR_COLORGAMUT_UNSPECIFIED:
453 return nullptr;
454 }
455 break;
456 case JPEGR_COLORGAMUT_UNSPECIFIED:
457 return nullptr;
458 }
459}
460
Nick Deakin594a4ca2022-11-16 20:57:42 -0500461
462////////////////////////////////////////////////////////////////////////////////
463// Recovery map calculations
464
Nick Deakin5c20b9e2022-11-15 17:39:24 -0500465uint8_t encodeRecovery(float y_sdr, float y_hdr, float hdr_ratio) {
Nick Deakinf6bca5a2022-11-04 10:43:43 -0400466 float gain = 1.0f;
467 if (y_sdr > 0.0f) {
468 gain = y_hdr / y_sdr;
469 }
470
Nick Deakin65f492a2022-11-29 22:47:40 -0500471 if (gain < (1.0f / hdr_ratio)) gain = 1.0f / hdr_ratio;
Nick Deakinf6bca5a2022-11-04 10:43:43 -0400472 if (gain > hdr_ratio) gain = hdr_ratio;
473
474 return static_cast<uint8_t>(log2(gain) / log2(hdr_ratio) * 127.5f + 127.5f);
475}
476
Nick Deakin5c20b9e2022-11-15 17:39:24 -0500477Color applyRecovery(Color e, float recovery, float hdr_ratio) {
Harish Mahendrakara5ddcc22022-12-13 12:45:23 -0800478 float recoveryFactor = pow(hdr_ratio, recovery);
479 return e * recoveryFactor;
Nick Deakinf6bca5a2022-11-04 10:43:43 -0400480}
481
Harish Mahendrakarf25991f2022-12-16 11:57:44 -0800482Color applyRecoveryLUT(Color e, float recovery, RecoveryLUT& recoveryLUT) {
483 float recoveryFactor = recoveryLUT.getRecoveryFactor(recovery);
484 return e * recoveryFactor;
485}
486
Nick Deakinf6bca5a2022-11-04 10:43:43 -0400487Color getYuv420Pixel(jr_uncompressed_ptr image, size_t x, size_t y) {
488 size_t pixel_count = image->width * image->height;
489
490 size_t pixel_y_idx = x + y * image->width;
491 size_t pixel_uv_idx = x / 2 + (y / 2) * (image->width / 2);
492
493 uint8_t y_uint = reinterpret_cast<uint8_t*>(image->data)[pixel_y_idx];
494 uint8_t u_uint = reinterpret_cast<uint8_t*>(image->data)[pixel_count + pixel_uv_idx];
495 uint8_t v_uint = reinterpret_cast<uint8_t*>(image->data)[pixel_count * 5 / 4 + pixel_uv_idx];
496
497 // 128 bias for UV given we are using jpeglib; see:
498 // https://github.com/kornelski/libjpeg/blob/master/structure.doc
499 return {{{ static_cast<float>(y_uint) / 255.0f,
500 (static_cast<float>(u_uint) - 128.0f) / 255.0f,
501 (static_cast<float>(v_uint) - 128.0f) / 255.0f }}};
502}
503
Nick Deakin594a4ca2022-11-16 20:57:42 -0500504Color getP010Pixel(jr_uncompressed_ptr image, size_t x, size_t y) {
505 size_t pixel_count = image->width * image->height;
Nick Deakinf6bca5a2022-11-04 10:43:43 -0400506
Nick Deakin594a4ca2022-11-16 20:57:42 -0500507 size_t pixel_y_idx = x + y * image->width;
508 size_t pixel_uv_idx = x / 2 + (y / 2) * (image->width / 2);
509
Nick Deakin86207ba2022-11-21 16:07:36 -0500510 uint16_t y_uint = reinterpret_cast<uint16_t*>(image->data)[pixel_y_idx]
511 >> 6;
512 uint16_t u_uint = reinterpret_cast<uint16_t*>(image->data)[pixel_count + pixel_uv_idx * 2]
513 >> 6;
514 uint16_t v_uint = reinterpret_cast<uint16_t*>(image->data)[pixel_count + pixel_uv_idx * 2 + 1]
515 >> 6;
Nick Deakin594a4ca2022-11-16 20:57:42 -0500516
517 // Conversions include taking narrow-range into account.
Nick Deakin38125332022-12-12 15:48:24 -0500518 return {{{ (static_cast<float>(y_uint) - 64.0f) / 876.0f,
519 (static_cast<float>(u_uint) - 64.0f) / 896.0f - 0.5f,
520 (static_cast<float>(v_uint) - 64.0f) / 896.0f - 0.5f }}};
Nick Deakin594a4ca2022-11-16 20:57:42 -0500521}
522
523typedef Color (*getPixelFn)(jr_uncompressed_ptr, size_t, size_t);
524
525static Color samplePixels(jr_uncompressed_ptr image, size_t map_scale_factor, size_t x, size_t y,
526 getPixelFn get_pixel_fn) {
527 Color e = {{{ 0.0f, 0.0f, 0.0f }}};
Nick Deakinf6bca5a2022-11-04 10:43:43 -0400528 for (size_t dy = 0; dy < map_scale_factor; ++dy) {
529 for (size_t dx = 0; dx < map_scale_factor; ++dx) {
Nick Deakin594a4ca2022-11-16 20:57:42 -0500530 e += get_pixel_fn(image, x * map_scale_factor + dx, y * map_scale_factor + dy);
Nick Deakinf6bca5a2022-11-04 10:43:43 -0400531 }
532 }
533
534 return e / static_cast<float>(map_scale_factor * map_scale_factor);
535}
536
Nick Deakin594a4ca2022-11-16 20:57:42 -0500537Color sampleYuv420(jr_uncompressed_ptr image, size_t map_scale_factor, size_t x, size_t y) {
538 return samplePixels(image, map_scale_factor, x, y, getYuv420Pixel);
Nick Deakinf6bca5a2022-11-04 10:43:43 -0400539}
540
Nick Deakin594a4ca2022-11-16 20:57:42 -0500541Color sampleP010(jr_uncompressed_ptr image, size_t map_scale_factor, size_t x, size_t y) {
542 return samplePixels(image, map_scale_factor, x, y, getP010Pixel);
Nick Deakinf6bca5a2022-11-04 10:43:43 -0400543}
Nick Deakin5c20b9e2022-11-15 17:39:24 -0500544
Nick Deakin65f492a2022-11-29 22:47:40 -0500545// TODO: do we need something more clever for filtering either the map or images
546// to generate the map?
547
548static size_t clamp(const size_t& val, const size_t& low, const size_t& high) {
549 return val < low ? low : (high < val ? high : val);
550}
551
552static float mapUintToFloat(uint8_t map_uint) {
553 return (static_cast<float>(map_uint) - 127.5f) / 127.5f;
554}
555
556static float pythDistance(float x_diff, float y_diff) {
557 return sqrt(pow(x_diff, 2.0f) + pow(y_diff, 2.0f));
558}
559
Ram Mohanfe723d62022-12-15 00:59:11 +0530560// TODO: If map_scale_factor is guaranteed to be an integer, then remove the following.
Nick Deakin65f492a2022-11-29 22:47:40 -0500561float sampleMap(jr_uncompressed_ptr map, size_t map_scale_factor, size_t x, size_t y) {
562 float x_map = static_cast<float>(x) / static_cast<float>(map_scale_factor);
563 float y_map = static_cast<float>(y) / static_cast<float>(map_scale_factor);
564
565 size_t x_lower = static_cast<size_t>(floor(x_map));
566 size_t x_upper = x_lower + 1;
567 size_t y_lower = static_cast<size_t>(floor(y_map));
568 size_t y_upper = y_lower + 1;
569
570 x_lower = clamp(x_lower, 0, map->width - 1);
571 x_upper = clamp(x_upper, 0, map->width - 1);
572 y_lower = clamp(y_lower, 0, map->height - 1);
573 y_upper = clamp(y_upper, 0, map->height - 1);
574
575 // Use Shepard's method for inverse distance weighting. For more information:
576 // en.wikipedia.org/wiki/Inverse_distance_weighting#Shepard's_method
577
578 float e1 = mapUintToFloat(reinterpret_cast<uint8_t*>(map->data)[x_lower + y_lower * map->width]);
579 float e1_dist = pythDistance(x_map - static_cast<float>(x_lower),
580 y_map - static_cast<float>(y_lower));
581 if (e1_dist == 0.0f) return e1;
582
583 float e2 = mapUintToFloat(reinterpret_cast<uint8_t*>(map->data)[x_lower + y_upper * map->width]);
584 float e2_dist = pythDistance(x_map - static_cast<float>(x_lower),
585 y_map - static_cast<float>(y_upper));
586 if (e2_dist == 0.0f) return e2;
587
588 float e3 = mapUintToFloat(reinterpret_cast<uint8_t*>(map->data)[x_upper + y_lower * map->width]);
589 float e3_dist = pythDistance(x_map - static_cast<float>(x_upper),
590 y_map - static_cast<float>(y_lower));
591 if (e3_dist == 0.0f) return e3;
592
593 float e4 = mapUintToFloat(reinterpret_cast<uint8_t*>(map->data)[x_upper + y_upper * map->width]);
594 float e4_dist = pythDistance(x_map - static_cast<float>(x_upper),
595 y_map - static_cast<float>(y_upper));
596 if (e4_dist == 0.0f) return e2;
597
598 float e1_weight = 1.0f / e1_dist;
599 float e2_weight = 1.0f / e2_dist;
600 float e3_weight = 1.0f / e3_dist;
601 float e4_weight = 1.0f / e4_dist;
602 float total_weight = e1_weight + e2_weight + e3_weight + e4_weight;
603
604 return e1 * (e1_weight / total_weight)
605 + e2 * (e2_weight / total_weight)
606 + e3 * (e3_weight / total_weight)
607 + e4 * (e4_weight / total_weight);
608}
609
Ram Mohanfe723d62022-12-15 00:59:11 +0530610float sampleMap(jr_uncompressed_ptr map, size_t map_scale_factor, size_t x, size_t y,
611 ShepardsIDW& weightTables) {
612 // TODO: If map_scale_factor is guaranteed to be an integer power of 2, then optimize the
613 // following by computing log2(map_scale_factor) once and then using >> log2(map_scale_factor)
614 int x_lower = x / map_scale_factor;
615 int x_upper = x_lower + 1;
616 int y_lower = y / map_scale_factor;
617 int y_upper = y_lower + 1;
618
619 x_lower = std::min(x_lower, map->width - 1);
620 x_upper = std::min(x_upper, map->width - 1);
621 y_lower = std::min(y_lower, map->height - 1);
622 y_upper = std::min(y_upper, map->height - 1);
623
624 float e1 = mapUintToFloat(reinterpret_cast<uint8_t*>(map->data)[x_lower + y_lower * map->width]);
625 float e2 = mapUintToFloat(reinterpret_cast<uint8_t*>(map->data)[x_lower + y_upper * map->width]);
626 float e3 = mapUintToFloat(reinterpret_cast<uint8_t*>(map->data)[x_upper + y_lower * map->width]);
627 float e4 = mapUintToFloat(reinterpret_cast<uint8_t*>(map->data)[x_upper + y_upper * map->width]);
628
629 // TODO: If map_scale_factor is guaranteed to be an integer power of 2, then optimize the
630 // following by using & (map_scale_factor - 1)
631 int offset_x = x % map_scale_factor;
632 int offset_y = y % map_scale_factor;
633
634 float* weights = weightTables.mWeights;
635 if (x_lower == x_upper && y_lower == y_upper) weights = weightTables.mWeightsC;
636 else if (x_lower == x_upper) weights = weightTables.mWeightsNR;
637 else if (y_lower == y_upper) weights = weightTables.mWeightsNB;
638 weights += offset_y * map_scale_factor * 4 + offset_x * 4;
639
640 return e1 * weights[0] + e2 * weights[1] + e3 * weights[2] + e4 * weights[3];
641}
642
Nick Deakin6bd90432022-11-20 16:26:37 -0500643uint32_t colorToRgba1010102(Color e_gamma) {
644 return (0x3ff & static_cast<uint32_t>(e_gamma.r * 1023.0f))
645 | ((0x3ff & static_cast<uint32_t>(e_gamma.g * 1023.0f)) << 10)
646 | ((0x3ff & static_cast<uint32_t>(e_gamma.b * 1023.0f)) << 20)
647 | (0x3 << 30); // Set alpha to 1.0
648}
649
Nick Deakinf6bca5a2022-11-04 10:43:43 -0400650} // namespace android::recoverymap