Coverage Report

Created: 2025-06-13 06:37

/src/libjxl/lib/jxl/splines.cc
Line
Count
Source (jump to first uncovered line)
1
// Copyright (c) the JPEG XL Project Authors. All rights reserved.
2
//
3
// Use of this source code is governed by a BSD-style
4
// license that can be found in the LICENSE file.
5
6
#include "lib/jxl/splines.h"
7
8
#include <jxl/memory_manager.h>
9
10
#include <algorithm>
11
#include <cinttypes>  // PRIu64
12
#include <cmath>
13
#include <cstddef>
14
#include <cstdint>
15
#include <cstring>
16
#include <limits>
17
#include <utility>
18
#include <vector>
19
20
#include "lib/jxl/base/bits.h"
21
#include "lib/jxl/base/common.h"
22
#include "lib/jxl/base/compiler_specific.h"  // ssize_t
23
#include "lib/jxl/base/printf_macros.h"
24
#include "lib/jxl/base/rect.h"
25
#include "lib/jxl/base/status.h"
26
#include "lib/jxl/chroma_from_luma.h"
27
#include "lib/jxl/common.h"  // JXL_HIGH_PRECISION
28
#include "lib/jxl/dct_scales.h"
29
#include "lib/jxl/dec_ans.h"
30
#include "lib/jxl/dec_bit_reader.h"
31
#include "lib/jxl/image.h"
32
#include "lib/jxl/pack_signed.h"
33
34
#undef HWY_TARGET_INCLUDE
35
#define HWY_TARGET_INCLUDE "lib/jxl/splines.cc"
36
#include <hwy/foreach_target.h>
37
#include <hwy/highway.h>
38
39
#include "lib/jxl/base/fast_math-inl.h"
40
HWY_BEFORE_NAMESPACE();
41
namespace jxl {
42
namespace HWY_NAMESPACE {
43
namespace {
44
45
// These templates are not found via ADL.
46
using hwy::HWY_NAMESPACE::Mul;
47
using hwy::HWY_NAMESPACE::MulAdd;
48
using hwy::HWY_NAMESPACE::MulSub;
49
using hwy::HWY_NAMESPACE::Sqrt;
50
using hwy::HWY_NAMESPACE::Sub;
51
52
// Given a set of DCT coefficients, this returns the result of performing cosine
53
// interpolation on the original samples.
54
0
float ContinuousIDCT(const Dct32& dct, const float t) {
55
  // We compute here the DCT-3 of the `dct` vector, rescaled by a factor of
56
  // sqrt(32). This is such that an input vector vector {x, 0, ..., 0} produces
57
  // a constant result of x. dct[0] was scaled in Dequantize() to allow uniform
58
  // treatment of all the coefficients.
59
0
  constexpr float kMultipliers[32] = {
60
0
      kPi / 32 * 0,  kPi / 32 * 1,  kPi / 32 * 2,  kPi / 32 * 3,  kPi / 32 * 4,
61
0
      kPi / 32 * 5,  kPi / 32 * 6,  kPi / 32 * 7,  kPi / 32 * 8,  kPi / 32 * 9,
62
0
      kPi / 32 * 10, kPi / 32 * 11, kPi / 32 * 12, kPi / 32 * 13, kPi / 32 * 14,
63
0
      kPi / 32 * 15, kPi / 32 * 16, kPi / 32 * 17, kPi / 32 * 18, kPi / 32 * 19,
64
0
      kPi / 32 * 20, kPi / 32 * 21, kPi / 32 * 22, kPi / 32 * 23, kPi / 32 * 24,
65
0
      kPi / 32 * 25, kPi / 32 * 26, kPi / 32 * 27, kPi / 32 * 28, kPi / 32 * 29,
66
0
      kPi / 32 * 30, kPi / 32 * 31,
67
0
  };
68
0
  HWY_CAPPED(float, 32) df;
69
0
  auto result = Zero(df);
70
0
  const auto tandhalf = Set(df, t + 0.5f);
71
0
  for (int i = 0; i < 32; i += Lanes(df)) {
72
0
    auto cos_arg = Mul(LoadU(df, kMultipliers + i), tandhalf);
73
0
    auto cos = FastCosf(df, cos_arg);
74
0
    auto local_res = Mul(LoadU(df, dct.data() + i), cos);
75
0
    result = MulAdd(Set(df, kSqrt2), local_res, result);
76
0
  }
77
0
  return GetLane(SumOfLanes(df, result));
78
0
}
Unexecuted instantiation: splines.cc:jxl::N_AVX2::(anonymous namespace)::ContinuousIDCT(std::__1::array<float, 32ul> const&, float)
Unexecuted instantiation: splines.cc:jxl::N_SSE4::(anonymous namespace)::ContinuousIDCT(std::__1::array<float, 32ul> const&, float)
Unexecuted instantiation: splines.cc:jxl::N_SSE2::(anonymous namespace)::ContinuousIDCT(std::__1::array<float, 32ul> const&, float)
79
80
template <typename DF>
81
void DrawSegment(DF df, const SplineSegment& segment, const bool add,
82
                 const size_t y, const size_t x, const size_t x0,
83
0
                 float* JXL_RESTRICT rows[3]) {
84
0
  Rebind<int32_t, DF> di;
85
0
  const auto inv_sigma = Set(df, segment.inv_sigma);
86
0
  const auto half = Set(df, 0.5f);
87
0
  const auto one_over_2s2 = Set(df, 0.353553391f);
88
0
  const auto sigma_over_4_times_intensity =
89
0
      Set(df, segment.sigma_over_4_times_intensity);
90
0
  const auto dx =
91
0
      Sub(ConvertTo(df, Iota(di, x + x0)), Set(df, segment.center_x));
92
0
  const auto dy = Set(df, y - segment.center_y);
93
0
  const auto sqd = MulAdd(dx, dx, Mul(dy, dy));
94
0
  const auto distance = Sqrt(sqd);
95
0
  const auto one_dimensional_factor =
96
0
      Sub(FastErff(df, Mul(MulAdd(distance, half, one_over_2s2), inv_sigma)),
97
0
          FastErff(df, Mul(MulSub(distance, half, one_over_2s2), inv_sigma)));
98
0
  auto local_intensity =
99
0
      Mul(sigma_over_4_times_intensity,
100
0
          Mul(one_dimensional_factor, one_dimensional_factor));
101
0
  for (size_t c = 0; c < 3; ++c) {
102
0
    const auto cm = Set(df, add ? segment.color[c] : -segment.color[c]);
103
0
    const auto in = LoadU(df, rows[c] + x);
104
0
    StoreU(MulAdd(cm, local_intensity, in), df, rows[c] + x);
105
0
  }
106
0
}
Unexecuted instantiation: splines.cc:void jxl::N_AVX2::(anonymous namespace)::DrawSegment<hwy::N_AVX2::Simd<float, 8ul, 0> >(hwy::N_AVX2::Simd<float, 8ul, 0>, jxl::SplineSegment const&, bool, unsigned long, unsigned long, unsigned long, float* restrict*)
Unexecuted instantiation: splines.cc:void jxl::N_AVX2::(anonymous namespace)::DrawSegment<hwy::N_AVX2::Simd<float, 1ul, 0> >(hwy::N_AVX2::Simd<float, 1ul, 0>, jxl::SplineSegment const&, bool, unsigned long, unsigned long, unsigned long, float* restrict*)
Unexecuted instantiation: splines.cc:void jxl::N_SSE4::(anonymous namespace)::DrawSegment<hwy::N_SSE4::Simd<float, 4ul, 0> >(hwy::N_SSE4::Simd<float, 4ul, 0>, jxl::SplineSegment const&, bool, unsigned long, unsigned long, unsigned long, float* restrict*)
Unexecuted instantiation: splines.cc:void jxl::N_SSE4::(anonymous namespace)::DrawSegment<hwy::N_SSE4::Simd<float, 1ul, 0> >(hwy::N_SSE4::Simd<float, 1ul, 0>, jxl::SplineSegment const&, bool, unsigned long, unsigned long, unsigned long, float* restrict*)
Unexecuted instantiation: splines.cc:void jxl::N_SSE2::(anonymous namespace)::DrawSegment<hwy::N_SSE2::Simd<float, 4ul, 0> >(hwy::N_SSE2::Simd<float, 4ul, 0>, jxl::SplineSegment const&, bool, unsigned long, unsigned long, unsigned long, float* restrict*)
Unexecuted instantiation: splines.cc:void jxl::N_SSE2::(anonymous namespace)::DrawSegment<hwy::N_SSE2::Simd<float, 1ul, 0> >(hwy::N_SSE2::Simd<float, 1ul, 0>, jxl::SplineSegment const&, bool, unsigned long, unsigned long, unsigned long, float* restrict*)
107
108
void DrawSegment(const SplineSegment& segment, const bool add, const size_t y,
109
                 const size_t x0, const size_t x1,
110
0
                 float* JXL_RESTRICT rows[3]) {
111
0
  ssize_t start = std::llround(segment.center_x - segment.maximum_distance);
112
0
  ssize_t end = std::llround(segment.center_x + segment.maximum_distance);
113
0
  if (end < static_cast<ssize_t>(x0) || start >= static_cast<ssize_t>(x1)) {
114
0
    return;  // span does not intersect scan
115
0
  }
116
0
  size_t span_x0 = std::max<ssize_t>(x0, start) - x0;
117
0
  size_t span_x1 = std::min<ssize_t>(x1, end + 1) - x0;  // exclusive
118
0
  HWY_FULL(float) df;
119
0
  size_t x = span_x0;
120
0
  for (; x + Lanes(df) <= span_x1; x += Lanes(df)) {
121
0
    DrawSegment(df, segment, add, y, x, x0, rows);
122
0
  }
123
0
  for (; x < span_x1; ++x) {
124
0
    DrawSegment(HWY_CAPPED(float, 1)(), segment, add, y, x, x0, rows);
125
0
  }
126
0
}
Unexecuted instantiation: splines.cc:jxl::N_AVX2::(anonymous namespace)::DrawSegment(jxl::SplineSegment const&, bool, unsigned long, unsigned long, unsigned long, float* restrict*)
Unexecuted instantiation: splines.cc:jxl::N_SSE4::(anonymous namespace)::DrawSegment(jxl::SplineSegment const&, bool, unsigned long, unsigned long, unsigned long, float* restrict*)
Unexecuted instantiation: splines.cc:jxl::N_SSE2::(anonymous namespace)::DrawSegment(jxl::SplineSegment const&, bool, unsigned long, unsigned long, unsigned long, float* restrict*)
127
128
void ComputeSegments(const Spline::Point& center, const float intensity,
129
                     const float color[3], const float sigma,
130
                     std::vector<SplineSegment>& segments,
131
0
                     std::vector<std::pair<size_t, size_t>>& segments_by_y) {
132
  // Sanity check sigma, inverse sigma and intensity
133
0
  if (!(std::isfinite(sigma) && sigma != 0.0f && std::isfinite(1.0f / sigma) &&
134
0
        std::isfinite(intensity))) {
135
0
    return;
136
0
  }
137
0
#if JXL_HIGH_PRECISION
138
0
  constexpr float kDistanceExp = 5;
139
#else
140
  // About 30% faster.
141
  constexpr float kDistanceExp = 3;
142
#endif
143
  // We cap from below colors to at least 0.01.
144
0
  float max_color = 0.01f;
145
0
  for (size_t c = 0; c < 3; c++) {
146
0
    max_color = std::max(max_color, std::abs(color[c] * intensity));
147
0
  }
148
  // Distance beyond which max_color*intensity*exp(-d^2 / (2 * sigma^2)) drops
149
  // below 10^-kDistanceExp.
150
0
  const float maximum_distance =
151
0
      std::sqrt(-2 * sigma * sigma *
152
0
                (std::log(0.1) * kDistanceExp - std::log(max_color)));
153
0
  SplineSegment segment;
154
0
  segment.center_y = center.y;
155
0
  segment.center_x = center.x;
156
0
  memcpy(segment.color, color, sizeof(segment.color));
157
0
  segment.inv_sigma = 1.0f / sigma;
158
0
  segment.sigma_over_4_times_intensity = .25f * sigma * intensity;
159
0
  segment.maximum_distance = maximum_distance;
160
0
  ssize_t y0 = std::llround(center.y - maximum_distance);
161
0
  ssize_t y1 =
162
0
      std::llround(center.y + maximum_distance) + 1;  // one-past-the-end
163
0
  for (ssize_t y = std::max<ssize_t>(y0, 0); y < y1; y++) {
164
0
    segments_by_y.emplace_back(y, segments.size());
165
0
  }
166
0
  segments.push_back(segment);
167
0
}
Unexecuted instantiation: splines.cc:jxl::N_AVX2::(anonymous namespace)::ComputeSegments(jxl::Spline::Point const&, float, float const*, float, std::__1::vector<jxl::SplineSegment, std::__1::allocator<jxl::SplineSegment> >&, std::__1::vector<std::__1::pair<unsigned long, unsigned long>, std::__1::allocator<std::__1::pair<unsigned long, unsigned long> > >&)
Unexecuted instantiation: splines.cc:jxl::N_SSE4::(anonymous namespace)::ComputeSegments(jxl::Spline::Point const&, float, float const*, float, std::__1::vector<jxl::SplineSegment, std::__1::allocator<jxl::SplineSegment> >&, std::__1::vector<std::__1::pair<unsigned long, unsigned long>, std::__1::allocator<std::__1::pair<unsigned long, unsigned long> > >&)
Unexecuted instantiation: splines.cc:jxl::N_SSE2::(anonymous namespace)::ComputeSegments(jxl::Spline::Point const&, float, float const*, float, std::__1::vector<jxl::SplineSegment, std::__1::allocator<jxl::SplineSegment> >&, std::__1::vector<std::__1::pair<unsigned long, unsigned long>, std::__1::allocator<std::__1::pair<unsigned long, unsigned long> > >&)
168
169
void DrawSegments(float* JXL_RESTRICT row_x, float* JXL_RESTRICT row_y,
170
                  float* JXL_RESTRICT row_b, size_t y, size_t x0, size_t x1,
171
                  const bool add, const SplineSegment* segments,
172
                  const size_t* segment_indices,
173
0
                  const size_t* segment_y_start) {
174
0
  float* JXL_RESTRICT rows[3] = {row_x, row_y, row_b};
175
0
  for (size_t i = segment_y_start[y]; i < segment_y_start[y + 1]; i++) {
176
0
    DrawSegment(segments[segment_indices[i]], add, y, x0, x1, rows);
177
0
  }
178
0
}
Unexecuted instantiation: splines.cc:jxl::N_AVX2::(anonymous namespace)::DrawSegments(float*, float*, float*, unsigned long, unsigned long, unsigned long, bool, jxl::SplineSegment const*, unsigned long const*, unsigned long const*)
Unexecuted instantiation: splines.cc:jxl::N_SSE4::(anonymous namespace)::DrawSegments(float*, float*, float*, unsigned long, unsigned long, unsigned long, bool, jxl::SplineSegment const*, unsigned long const*, unsigned long const*)
Unexecuted instantiation: splines.cc:jxl::N_SSE2::(anonymous namespace)::DrawSegments(float*, float*, float*, unsigned long, unsigned long, unsigned long, bool, jxl::SplineSegment const*, unsigned long const*, unsigned long const*)
179
180
void SegmentsFromPoints(
181
    const Spline& spline,
182
    const std::vector<std::pair<Spline::Point, float>>& points_to_draw,
183
    const float arc_length, std::vector<SplineSegment>& segments,
184
0
    std::vector<std::pair<size_t, size_t>>& segments_by_y) {
185
0
  const float inv_arc_length = 1.0f / arc_length;
186
0
  int k = 0;
187
0
  for (const auto& point_to_draw : points_to_draw) {
188
0
    const Spline::Point& point = point_to_draw.first;
189
0
    const float multiplier = point_to_draw.second;
190
0
    const float progress_along_arc =
191
0
        std::min(1.f, (k * kDesiredRenderingDistance) * inv_arc_length);
192
0
    ++k;
193
0
    float color[3];
194
0
    for (size_t c = 0; c < 3; ++c) {
195
0
      color[c] =
196
0
          ContinuousIDCT(spline.color_dct[c], (32 - 1) * progress_along_arc);
197
0
    }
198
0
    const float sigma =
199
0
        ContinuousIDCT(spline.sigma_dct, (32 - 1) * progress_along_arc);
200
0
    ComputeSegments(point, multiplier, color, sigma, segments, segments_by_y);
201
0
  }
202
0
}
Unexecuted instantiation: splines.cc:jxl::N_AVX2::(anonymous namespace)::SegmentsFromPoints(jxl::Spline const&, std::__1::vector<std::__1::pair<jxl::Spline::Point, float>, std::__1::allocator<std::__1::pair<jxl::Spline::Point, float> > > const&, float, std::__1::vector<jxl::SplineSegment, std::__1::allocator<jxl::SplineSegment> >&, std::__1::vector<std::__1::pair<unsigned long, unsigned long>, std::__1::allocator<std::__1::pair<unsigned long, unsigned long> > >&)
Unexecuted instantiation: splines.cc:jxl::N_SSE4::(anonymous namespace)::SegmentsFromPoints(jxl::Spline const&, std::__1::vector<std::__1::pair<jxl::Spline::Point, float>, std::__1::allocator<std::__1::pair<jxl::Spline::Point, float> > > const&, float, std::__1::vector<jxl::SplineSegment, std::__1::allocator<jxl::SplineSegment> >&, std::__1::vector<std::__1::pair<unsigned long, unsigned long>, std::__1::allocator<std::__1::pair<unsigned long, unsigned long> > >&)
Unexecuted instantiation: splines.cc:jxl::N_SSE2::(anonymous namespace)::SegmentsFromPoints(jxl::Spline const&, std::__1::vector<std::__1::pair<jxl::Spline::Point, float>, std::__1::allocator<std::__1::pair<jxl::Spline::Point, float> > > const&, float, std::__1::vector<jxl::SplineSegment, std::__1::allocator<jxl::SplineSegment> >&, std::__1::vector<std::__1::pair<unsigned long, unsigned long>, std::__1::allocator<std::__1::pair<unsigned long, unsigned long> > >&)
203
}  // namespace
204
// NOLINTNEXTLINE(google-readability-namespace-comments)
205
}  // namespace HWY_NAMESPACE
206
}  // namespace jxl
207
HWY_AFTER_NAMESPACE();
208
209
#if HWY_ONCE
210
namespace jxl {
211
HWY_EXPORT(SegmentsFromPoints);
212
HWY_EXPORT(DrawSegments);
213
214
namespace {
215
216
// It is not in spec, but reasonable limit to avoid overflows.
217
template <typename T>
218
0
Status ValidateSplinePointPos(const T& x, const T& y) {
219
0
  constexpr T kSplinePosLimit = 1u << 23;
220
0
  if ((x >= kSplinePosLimit) || (x <= -kSplinePosLimit) ||
221
0
      (y >= kSplinePosLimit) || (y <= -kSplinePosLimit)) {
222
0
    return JXL_FAILURE("Spline coordinates out of bounds");
223
0
  }
224
0
  return true;
225
0
}
Unexecuted instantiation: splines.cc:jxl::Status jxl::(anonymous namespace)::ValidateSplinePointPos<long>(long const&, long const&)
Unexecuted instantiation: splines.cc:jxl::Status jxl::(anonymous namespace)::ValidateSplinePointPos<float>(float const&, float const&)
Unexecuted instantiation: splines.cc:jxl::Status jxl::(anonymous namespace)::ValidateSplinePointPos<int>(int const&, int const&)
226
227
// Maximum number of spline control points per frame is
228
//   std::min(kMaxNumControlPoints, xsize * ysize / 2)
229
constexpr size_t kMaxNumControlPoints = 1u << 20u;
230
constexpr size_t kMaxNumControlPointsPerPixelRatio = 2;
231
232
0
float AdjustedQuant(const int32_t adjustment) {
233
0
  return (adjustment >= 0) ? (1.f + .125f * adjustment)
234
0
                           : 1.f / (1.f - .125f * adjustment);
235
0
}
236
237
0
float InvAdjustedQuant(const int32_t adjustment) {
238
0
  return (adjustment >= 0) ? 1.f / (1.f + .125f * adjustment)
239
0
                           : (1.f - .125f * adjustment);
240
0
}
241
242
// X, Y, B, sigma.
243
constexpr float kChannelWeight[] = {0.0042f, 0.075f, 0.07f, .3333f};
244
245
Status DecodeAllStartingPoints(std::vector<Spline::Point>* const points,
246
                               BitReader* const br, ANSSymbolReader* reader,
247
                               const std::vector<uint8_t>& context_map,
248
0
                               const size_t num_splines) {
249
0
  points->clear();
250
0
  points->reserve(num_splines);
251
0
  int64_t last_x = 0;
252
0
  int64_t last_y = 0;
253
0
  for (size_t i = 0; i < num_splines; i++) {
254
0
    size_t dx =
255
0
        reader->ReadHybridUint(kStartingPositionContext, br, context_map);
256
0
    size_t dy =
257
0
        reader->ReadHybridUint(kStartingPositionContext, br, context_map);
258
0
    int64_t x;
259
0
    int64_t y;
260
0
    if (i != 0) {
261
0
      x = UnpackSigned(dx) + last_x;
262
0
      y = UnpackSigned(dy) + last_y;
263
0
    } else {
264
0
      x = dx;
265
0
      y = dy;
266
0
    }
267
0
    JXL_RETURN_IF_ERROR(ValidateSplinePointPos(x, y));
268
0
    points->emplace_back(static_cast<float>(x), static_cast<float>(y));
269
0
    last_x = x;
270
0
    last_y = y;
271
0
  }
272
0
  return true;
273
0
}
274
275
struct Vector {
276
  float x, y;
277
0
  Vector operator-() const { return {-x, -y}; }
278
0
  Vector operator+(const Vector& other) const {
279
0
    return {x + other.x, y + other.y};
280
0
  }
281
0
  float SquaredNorm() const { return x * x + y * y; }
282
};
283
0
Vector operator*(const float k, const Vector& vec) {
284
0
  return {k * vec.x, k * vec.y};
285
0
}
286
287
0
Spline::Point operator+(const Spline::Point& p, const Vector& vec) {
288
0
  return {p.x + vec.x, p.y + vec.y};
289
0
}
290
0
Vector operator-(const Spline::Point& a, const Spline::Point& b) {
291
0
  return {a.x - b.x, a.y - b.y};
292
0
}
293
294
// TODO(eustas): avoid making a copy of "points".
295
void DrawCentripetalCatmullRomSpline(std::vector<Spline::Point> points,
296
0
                                     std::vector<Spline::Point>& result) {
297
0
  if (points.empty()) return;
298
0
  if (points.size() == 1) {
299
0
    result.push_back(points[0]);
300
0
    return;
301
0
  }
302
  // Number of points to compute between each control point.
303
0
  static constexpr int kNumPoints = 16;
304
0
  result.reserve((points.size() - 1) * kNumPoints + 1);
305
0
  points.insert(points.begin(), points[0] + (points[0] - points[1]));
306
0
  points.push_back(points[points.size() - 1] +
307
0
                   (points[points.size() - 1] - points[points.size() - 2]));
308
  // points has at least 4 elements at this point.
309
0
  for (size_t start = 0; start < points.size() - 3; ++start) {
310
    // 4 of them are used, and we draw from p[1] to p[2].
311
0
    const Spline::Point* const p = &points[start];
312
0
    result.push_back(p[1]);
313
0
    float d[3];
314
0
    float t[4];
315
0
    t[0] = 0;
316
0
    for (int k = 0; k < 3; ++k) {
317
      // TODO(eustas): for each segment delta is calculated 3 times...
318
      // TODO(eustas): restrict d[k] with reasonable limit and spec it.
319
0
      d[k] = std::sqrt(hypotf(p[k + 1].x - p[k].x, p[k + 1].y - p[k].y));
320
0
      t[k + 1] = t[k] + d[k];
321
0
    }
322
0
    for (int i = 1; i < kNumPoints; ++i) {
323
0
      const float tt = d[0] + (static_cast<float>(i) / kNumPoints) * d[1];
324
0
      Spline::Point a[3];
325
0
      for (int k = 0; k < 3; ++k) {
326
        // TODO(eustas): reciprocal multiplication would be faster.
327
0
        a[k] = p[k] + ((tt - t[k]) / d[k]) * (p[k + 1] - p[k]);
328
0
      }
329
0
      Spline::Point b[2];
330
0
      for (int k = 0; k < 2; ++k) {
331
0
        b[k] = a[k] + ((tt - t[k]) / (d[k] + d[k + 1])) * (a[k + 1] - a[k]);
332
0
      }
333
0
      result.push_back(b[0] + ((tt - t[1]) / d[1]) * (b[1] - b[0]));
334
0
    }
335
0
  }
336
0
  result.push_back(points[points.size() - 2]);
337
0
}
338
339
// Move along the line segments defined by `points`, `kDesiredRenderingDistance`
340
// pixels at a time, and call `functor` with each point and the actual distance
341
// to the previous point (which will always be kDesiredRenderingDistance except
342
// possibly for the very last point).
343
// TODO(eustas): this method always adds the last point, but never the first
344
//               (unless those are one); I believe both ends matter.
345
template <typename Points, typename Functor>
346
0
Status ForEachEquallySpacedPoint(const Points& points, const Functor& functor) {
347
0
  JXL_ENSURE(!points.empty());
348
0
  Spline::Point current = points.front();
349
0
  functor(current, kDesiredRenderingDistance);
350
0
  auto next = points.begin();
351
0
  while (next != points.end()) {
352
0
    const Spline::Point* previous = &current;
353
0
    float arclength_from_previous = 0.f;
354
0
    for (;;) {
355
0
      if (next == points.end()) {
356
0
        functor(*previous, arclength_from_previous);
357
0
        return true;
358
0
      }
359
0
      const float arclength_to_next =
360
0
          std::sqrt((*next - *previous).SquaredNorm());
361
0
      if (arclength_from_previous + arclength_to_next >=
362
0
          kDesiredRenderingDistance) {
363
0
        current =
364
0
            *previous + ((kDesiredRenderingDistance - arclength_from_previous) /
365
0
                         arclength_to_next) *
366
0
                            (*next - *previous);
367
0
        functor(current, kDesiredRenderingDistance);
368
0
        break;
369
0
      }
370
0
      arclength_from_previous += arclength_to_next;
371
0
      previous = &*next;
372
0
      ++next;
373
0
    }
374
0
  }
375
0
  return true;
376
0
}
377
378
}  // namespace
379
380
StatusOr<QuantizedSpline> QuantizedSpline::Create(
381
    const Spline& original, const int32_t quantization_adjustment,
382
0
    const float y_to_x, const float y_to_b) {
383
0
  JXL_ENSURE(!original.control_points.empty());
384
0
  QuantizedSpline result;
385
0
  result.control_points_.reserve(original.control_points.size() - 1);
386
0
  const Spline::Point& starting_point = original.control_points.front();
387
0
  int previous_x = static_cast<int>(std::round(starting_point.x));
388
0
  int previous_y = static_cast<int>(std::round(starting_point.y));
389
0
  int previous_delta_x = 0;
390
0
  int previous_delta_y = 0;
391
0
  for (auto it = original.control_points.begin() + 1;
392
0
       it != original.control_points.end(); ++it) {
393
0
    const int new_x = static_cast<int>(std::round(it->x));
394
0
    const int new_y = static_cast<int>(std::round(it->y));
395
0
    const int new_delta_x = new_x - previous_x;
396
0
    const int new_delta_y = new_y - previous_y;
397
0
    result.control_points_.emplace_back(new_delta_x - previous_delta_x,
398
0
                                        new_delta_y - previous_delta_y);
399
0
    previous_delta_x = new_delta_x;
400
0
    previous_delta_y = new_delta_y;
401
0
    previous_x = new_x;
402
0
    previous_y = new_y;
403
0
  }
404
405
0
  const auto to_int = [](float v) -> int {
406
    // Maximal int representable with float.
407
0
    constexpr float kMax = std::numeric_limits<int>::max() - 127;
408
0
    constexpr float kMin = -kMax;
409
0
    return static_cast<int>(std::round(Clamp1(v, kMin, kMax)));
410
0
  };
411
412
0
  const auto quant = AdjustedQuant(quantization_adjustment);
413
0
  const auto inv_quant = InvAdjustedQuant(quantization_adjustment);
414
0
  for (int c : {1, 0, 2}) {
415
0
    float factor = (c == 0) ? y_to_x : (c == 1) ? 0 : y_to_b;
416
0
    for (int i = 0; i < 32; ++i) {
417
0
      const float dct_factor = (i == 0) ? kSqrt2 : 1.0f;
418
0
      const float inv_dct_factor = (i == 0) ? kSqrt0_5 : 1.0f;
419
0
      auto restored_y = result.color_dct_[1][i] * inv_dct_factor *
420
0
                        kChannelWeight[1] * inv_quant;
421
0
      auto decorrelated = original.color_dct[c][i] - factor * restored_y;
422
0
      result.color_dct_[c][i] =
423
0
          to_int(decorrelated * dct_factor * quant / kChannelWeight[c]);
424
0
    }
425
0
  }
426
0
  for (int i = 0; i < 32; ++i) {
427
0
    const float dct_factor = (i == 0) ? kSqrt2 : 1.0f;
428
0
    result.sigma_dct_[i] =
429
0
        to_int(original.sigma_dct[i] * dct_factor * quant / kChannelWeight[3]);
430
0
  }
431
0
  return result;
432
0
}
433
434
Status QuantizedSpline::Dequantize(const Spline::Point& starting_point,
435
                                   const int32_t quantization_adjustment,
436
                                   const float y_to_x, const float y_to_b,
437
                                   const uint64_t image_size,
438
                                   uint64_t* total_estimated_area_reached,
439
0
                                   Spline& result) const {
440
0
  constexpr uint64_t kOne = static_cast<uint64_t>(1);
441
0
  const uint64_t area_limit =
442
0
      std::min(1024 * image_size + (kOne << 32), kOne << 42);
443
444
0
  result.control_points.clear();
445
0
  result.control_points.reserve(control_points_.size() + 1);
446
0
  float px = std::round(starting_point.x);
447
0
  float py = std::round(starting_point.y);
448
0
  JXL_RETURN_IF_ERROR(ValidateSplinePointPos(px, py));
449
0
  int current_x = static_cast<int>(px);
450
0
  int current_y = static_cast<int>(py);
451
0
  result.control_points.emplace_back(static_cast<float>(current_x),
452
0
                                     static_cast<float>(current_y));
453
0
  int current_delta_x = 0;
454
0
  int current_delta_y = 0;
455
0
  uint64_t manhattan_distance = 0;
456
0
  for (const auto& point : control_points_) {
457
0
    current_delta_x += point.first;
458
0
    current_delta_y += point.second;
459
0
    manhattan_distance += std::abs(current_delta_x) + std::abs(current_delta_y);
460
0
    if (manhattan_distance > area_limit) {
461
0
      return JXL_FAILURE("Too large manhattan_distance reached: %" PRIu64,
462
0
                         manhattan_distance);
463
0
    }
464
0
    JXL_RETURN_IF_ERROR(
465
0
        ValidateSplinePointPos(current_delta_x, current_delta_y));
466
0
    current_x += current_delta_x;
467
0
    current_y += current_delta_y;
468
0
    JXL_RETURN_IF_ERROR(ValidateSplinePointPos(current_x, current_y));
469
0
    result.control_points.emplace_back(static_cast<float>(current_x),
470
0
                                       static_cast<float>(current_y));
471
0
  }
472
473
0
  const auto inv_quant = InvAdjustedQuant(quantization_adjustment);
474
0
  for (int c = 0; c < 3; ++c) {
475
0
    for (int i = 0; i < 32; ++i) {
476
0
      const float inv_dct_factor = (i == 0) ? kSqrt0_5 : 1.0f;
477
0
      result.color_dct[c][i] =
478
0
          color_dct_[c][i] * inv_dct_factor * kChannelWeight[c] * inv_quant;
479
0
    }
480
0
  }
481
0
  for (int i = 0; i < 32; ++i) {
482
0
    result.color_dct[0][i] += y_to_x * result.color_dct[1][i];
483
0
    result.color_dct[2][i] += y_to_b * result.color_dct[1][i];
484
0
  }
485
0
  uint64_t width_estimate = 0;
486
487
0
  uint64_t color[3] = {};
488
0
  for (int c = 0; c < 3; ++c) {
489
0
    for (int i = 0; i < 32; ++i) {
490
0
      color[c] += static_cast<uint64_t>(
491
0
          std::ceil(inv_quant * std::abs(color_dct_[c][i])));
492
0
    }
493
0
  }
494
0
  color[0] += static_cast<uint64_t>(std::ceil(std::abs(y_to_x))) * color[1];
495
0
  color[2] += static_cast<uint64_t>(std::ceil(std::abs(y_to_b))) * color[1];
496
  // This is not taking kChannelWeight into account, but up to constant factors
497
  // it gives an indication of the influence of the color values on the area
498
  // that will need to be rendered.
499
0
  const uint64_t max_color = std::max({color[1], color[0], color[2]});
500
0
  uint64_t logcolor =
501
0
      std::max(kOne, static_cast<uint64_t>(CeilLog2Nonzero(kOne + max_color)));
502
503
0
  const float weight_limit =
504
0
      std::ceil(std::sqrt((static_cast<float>(area_limit) / logcolor) /
505
0
                          std::max<size_t>(1, manhattan_distance)));
506
507
0
  for (int i = 0; i < 32; ++i) {
508
0
    const float inv_dct_factor = (i == 0) ? kSqrt0_5 : 1.0f;
509
0
    result.sigma_dct[i] =
510
0
        sigma_dct_[i] * inv_dct_factor * kChannelWeight[3] * inv_quant;
511
    // If we include the factor kChannelWeight[3]=.3333f here, we get a
512
    // realistic area estimate. We leave it out to simplify the calculations,
513
    // and understand that this way we underestimate the area by a factor of
514
    // 1/(0.3333*0.3333). This is taken into account in the limits below.
515
0
    float weight_f = std::ceil(inv_quant * std::abs(sigma_dct_[i]));
516
0
    uint64_t weight =
517
0
        static_cast<uint64_t>(std::min(weight_limit, std::max(1.0f, weight_f)));
518
0
    width_estimate += weight * weight * logcolor;
519
0
  }
520
0
  *total_estimated_area_reached += (width_estimate * manhattan_distance);
521
0
  if (*total_estimated_area_reached > area_limit) {
522
0
    return JXL_FAILURE("Too large total_estimated_area eached: %" PRIu64,
523
0
                       *total_estimated_area_reached);
524
0
  }
525
526
0
  return true;
527
0
}
528
529
Status QuantizedSpline::Decode(const std::vector<uint8_t>& context_map,
530
                               ANSSymbolReader* const decoder,
531
                               BitReader* const br,
532
                               const size_t max_control_points,
533
0
                               size_t* total_num_control_points) {
534
0
  const size_t num_control_points =
535
0
      decoder->ReadHybridUint(kNumControlPointsContext, br, context_map);
536
0
  if (num_control_points > max_control_points) {
537
0
    return JXL_FAILURE("Too many control points: %" PRIuS, num_control_points);
538
0
  }
539
0
  *total_num_control_points += num_control_points;
540
0
  if (*total_num_control_points > max_control_points) {
541
0
    return JXL_FAILURE("Too many control points: %" PRIuS,
542
0
                       *total_num_control_points);
543
0
  }
544
0
  control_points_.resize(num_control_points);
545
  // Maximal image dimension.
546
0
  constexpr int64_t kDeltaLimit = 1u << 30;
547
0
  for (std::pair<int64_t, int64_t>& control_point : control_points_) {
548
0
    control_point.first = UnpackSigned(
549
0
        decoder->ReadHybridUint(kControlPointsContext, br, context_map));
550
0
    control_point.second = UnpackSigned(
551
0
        decoder->ReadHybridUint(kControlPointsContext, br, context_map));
552
    // Check delta-deltas are not outrageous; it is not in spec, but there is
553
    // no reason to allow larger values.
554
0
    if ((control_point.first >= kDeltaLimit) ||
555
0
        (control_point.first <= -kDeltaLimit) ||
556
0
        (control_point.second >= kDeltaLimit) ||
557
0
        (control_point.second <= -kDeltaLimit)) {
558
0
      return JXL_FAILURE("Spline delta-delta is out of bounds");
559
0
    }
560
0
  }
561
562
0
  const auto decode_dct = [decoder, br, &context_map](int dct[32]) -> Status {
563
0
    constexpr int kWeirdNumber = std::numeric_limits<int>::min();
564
0
    for (int i = 0; i < 32; ++i) {
565
0
      dct[i] =
566
0
          UnpackSigned(decoder->ReadHybridUint(kDCTContext, br, context_map));
567
0
      if (dct[i] == kWeirdNumber) {
568
0
        return JXL_FAILURE("The weird number in spline DCT");
569
0
      }
570
0
    }
571
0
    return true;
572
0
  };
573
0
  for (auto& dct : color_dct_) {
574
0
    JXL_RETURN_IF_ERROR(decode_dct(dct));
575
0
  }
576
0
  JXL_RETURN_IF_ERROR(decode_dct(sigma_dct_));
577
0
  return true;
578
0
}
579
580
0
void Splines::Clear() {
581
0
  quantization_adjustment_ = 0;
582
0
  splines_.clear();
583
0
  starting_points_.clear();
584
0
  segments_.clear();
585
0
  segment_indices_.clear();
586
0
  segment_y_start_.clear();
587
0
}
588
589
Status Splines::Decode(JxlMemoryManager* memory_manager, jxl::BitReader* br,
590
0
                       const size_t num_pixels) {
591
0
  std::vector<uint8_t> context_map;
592
0
  ANSCode code;
593
0
  JXL_RETURN_IF_ERROR(DecodeHistograms(memory_manager, br, kNumSplineContexts,
594
0
                                       &code, &context_map));
595
0
  JXL_ASSIGN_OR_RETURN(ANSSymbolReader decoder,
596
0
                       ANSSymbolReader::Create(&code, br));
597
0
  size_t num_splines =
598
0
      decoder.ReadHybridUint(kNumSplinesContext, br, context_map);
599
0
  size_t max_control_points = std::min(
600
0
      kMaxNumControlPoints, num_pixels / kMaxNumControlPointsPerPixelRatio);
601
0
  if (num_splines > max_control_points ||
602
0
      num_splines + 1 > max_control_points) {
603
0
    return JXL_FAILURE("Too many splines: %" PRIuS, num_splines);
604
0
  }
605
0
  num_splines++;
606
0
  JXL_RETURN_IF_ERROR(DecodeAllStartingPoints(&starting_points_, br, &decoder,
607
0
                                              context_map, num_splines));
608
609
0
  quantization_adjustment_ = UnpackSigned(
610
0
      decoder.ReadHybridUint(kQuantizationAdjustmentContext, br, context_map));
611
612
0
  splines_.clear();
613
0
  splines_.reserve(num_splines);
614
0
  size_t num_control_points = num_splines;
615
0
  for (size_t i = 0; i < num_splines; ++i) {
616
0
    QuantizedSpline spline;
617
0
    JXL_RETURN_IF_ERROR(spline.Decode(context_map, &decoder, br,
618
0
                                      max_control_points, &num_control_points));
619
0
    splines_.push_back(std::move(spline));
620
0
  }
621
622
0
  JXL_RETURN_IF_ERROR(decoder.CheckANSFinalState());
623
624
0
  if (!HasAny()) {
625
0
    return JXL_FAILURE("Decoded splines but got none");
626
0
  }
627
628
0
  return true;
629
0
}
630
631
0
void Splines::AddTo(Image3F* const opsin, const Rect& opsin_rect) const {
632
0
  Apply</*add=*/true>(opsin, opsin_rect);
633
0
}
634
void Splines::AddToRow(float* JXL_RESTRICT row_x, float* JXL_RESTRICT row_y,
635
                       float* JXL_RESTRICT row_b, size_t y, size_t x0,
636
0
                       size_t x1) const {
637
0
  ApplyToRow</*add=*/true>(row_x, row_y, row_b, y, x0, x1);
638
0
}
639
640
268
void Splines::SubtractFrom(Image3F* const opsin) const {
641
268
  Apply</*add=*/false>(opsin, Rect(*opsin));
642
268
}
643
644
Status Splines::InitializeDrawCache(const size_t image_xsize,
645
                                    const size_t image_ysize,
646
268
                                    const ColorCorrelation& color_correlation) {
647
  // TODO(veluca): avoid storing segments that are entirely outside image
648
  // boundaries.
649
268
  segments_.clear();
650
268
  segment_indices_.clear();
651
268
  segment_y_start_.clear();
652
268
  std::vector<std::pair<size_t, size_t>> segments_by_y;
653
268
  std::vector<Spline::Point> intermediate_points;
654
268
  uint64_t total_estimated_area_reached = 0;
655
268
  std::vector<Spline> splines;
656
268
  for (size_t i = 0; i < splines_.size(); ++i) {
657
0
    Spline spline;
658
0
    JXL_RETURN_IF_ERROR(splines_[i].Dequantize(
659
0
        starting_points_[i], quantization_adjustment_,
660
0
        color_correlation.YtoXRatio(0), color_correlation.YtoBRatio(0),
661
0
        image_xsize * image_ysize, &total_estimated_area_reached, spline));
662
0
    if (std::adjacent_find(spline.control_points.begin(),
663
0
                           spline.control_points.end()) !=
664
0
        spline.control_points.end()) {
665
      // Otherwise division by zero might occur. Once control points coincide,
666
      // the direction of curve is undefined...
667
0
      return JXL_FAILURE(
668
0
          "identical successive control points in spline %" PRIuS, i);
669
0
    }
670
0
    splines.push_back(spline);
671
0
  }
672
  // TODO(firsching) Change this into a JXL_FAILURE for level 5 codestreams.
673
268
  if (total_estimated_area_reached >
674
268
      std::min(
675
268
          (8 * image_xsize * image_ysize + (static_cast<uint64_t>(1) << 25)),
676
268
          (static_cast<uint64_t>(1) << 30))) {
677
0
    JXL_WARNING(
678
0
        "Large total_estimated_area_reached, expect slower decoding: %" PRIu64,
679
0
        total_estimated_area_reached);
680
0
#ifdef FUZZING_BUILD_MODE_UNSAFE_FOR_PRODUCTION
681
0
    return JXL_FAILURE("Total spline area is too large");
682
0
#endif
683
0
  }
684
685
268
  for (Spline& spline : splines) {
686
0
    std::vector<std::pair<Spline::Point, float>> points_to_draw;
687
0
    auto add_point = [&](const Spline::Point& point, const float multiplier) {
688
0
      points_to_draw.emplace_back(point, multiplier);
689
0
    };
690
0
    intermediate_points.clear();
691
0
    DrawCentripetalCatmullRomSpline(spline.control_points, intermediate_points);
692
0
    JXL_RETURN_IF_ERROR(
693
0
        ForEachEquallySpacedPoint(intermediate_points, add_point));
694
0
    const float arc_length =
695
0
        (points_to_draw.size() - 2) * kDesiredRenderingDistance +
696
0
        points_to_draw.back().second;
697
0
    if (arc_length <= 0.f) {
698
      // This spline wouldn't have any effect.
699
0
      continue;
700
0
    }
701
0
    HWY_DYNAMIC_DISPATCH(SegmentsFromPoints)
702
0
    (spline, points_to_draw, arc_length, segments_, segments_by_y);
703
0
  }
704
705
  // TODO(eustas): consider linear sorting here.
706
268
  std::sort(segments_by_y.begin(), segments_by_y.end());
707
268
  segment_indices_.resize(segments_by_y.size());
708
268
  segment_y_start_.resize(image_ysize + 1);
709
268
  for (size_t i = 0; i < segments_by_y.size(); i++) {
710
0
    segment_indices_[i] = segments_by_y[i].second;
711
0
    size_t y = segments_by_y[i].first;
712
0
    if (y < image_ysize) {
713
0
      segment_y_start_[y + 1]++;
714
0
    }
715
0
  }
716
2.41k
  for (size_t y = 0; y < image_ysize; y++) {
717
2.14k
    segment_y_start_[y + 1] += segment_y_start_[y];
718
2.14k
  }
719
268
  return true;
720
268
}
721
722
template <bool add>
723
void Splines::ApplyToRow(float* JXL_RESTRICT row_x, float* JXL_RESTRICT row_y,
724
                         float* JXL_RESTRICT row_b, size_t y, size_t x0,
725
0
                         size_t x1) const {
726
0
  if (segments_.empty()) return;
727
0
  HWY_DYNAMIC_DISPATCH(DrawSegments)
728
0
  (row_x, row_y, row_b, y, x0, x1, add, segments_.data(),
729
0
   segment_indices_.data(), segment_y_start_.data());
730
0
}
Unexecuted instantiation: void jxl::Splines::ApplyToRow<true>(float*, float*, float*, unsigned long, unsigned long, unsigned long) const
Unexecuted instantiation: void jxl::Splines::ApplyToRow<false>(float*, float*, float*, unsigned long, unsigned long, unsigned long) const
731
732
template <bool add>
733
268
void Splines::Apply(Image3F* const opsin, const Rect& opsin_rect) const {
734
268
  if (segments_.empty()) return;
735
0
  const size_t y0 = opsin_rect.y0();
736
0
  const size_t x0 = opsin_rect.x0();
737
0
  const size_t x1 = opsin_rect.x1();
738
0
  for (size_t y = 0; y < opsin_rect.ysize(); y++) {
739
0
    ApplyToRow<add>(opsin->PlaneRow(0, y0 + y) + x0,
740
0
                    opsin->PlaneRow(1, y0 + y) + x0,
741
0
                    opsin->PlaneRow(2, y0 + y) + x0, y0 + y, x0, x1);
742
0
  }
743
0
}
Unexecuted instantiation: void jxl::Splines::Apply<true>(jxl::Image3<float>*, jxl::RectT<unsigned long> const&) const
void jxl::Splines::Apply<false>(jxl::Image3<float>*, jxl::RectT<unsigned long> const&) const
Line
Count
Source
733
268
void Splines::Apply(Image3F* const opsin, const Rect& opsin_rect) const {
734
268
  if (segments_.empty()) return;
735
0
  const size_t y0 = opsin_rect.y0();
736
0
  const size_t x0 = opsin_rect.x0();
737
0
  const size_t x1 = opsin_rect.x1();
738
0
  for (size_t y = 0; y < opsin_rect.ysize(); y++) {
739
0
    ApplyToRow<add>(opsin->PlaneRow(0, y0 + y) + x0,
740
0
                    opsin->PlaneRow(1, y0 + y) + x0,
741
0
                    opsin->PlaneRow(2, y0 + y) + x0, y0 + y, x0, x1);
742
0
  }
743
0
}
744
745
}  // namespace jxl
746
#endif  // HWY_ONCE