1
/*
2
    RawSpeed - RAW file decoder.
3

4
    Copyright (C) 2018 Stefan Löffler
5
    Copyright (C) 2018-2019 Roman Lebedev
6

7
    This library is free software; you can redistribute it and/or
8
    modify it under the terms of the GNU Lesser General Public
9
    License as published by the Free Software Foundation; either
10
    version 2 of the License, or (at your option) any later version.
11

12
    This library is distributed in the hope that it will be useful,
13
    but WITHOUT ANY WARRANTY; without even the implied warranty of
14
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
15
    Lesser General Public License for more details.
16

17
    You should have received a copy of the GNU Lesser General Public
18
    License along with this library; if not, write to the Free Software
19
    Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
20
*/
21

22
/*
23
  This is a decompressor for VC-5 raw compression algo, as used in GoPro raws.
24
  This implementation is similar to that one of the official reference
25
  implementation of the https://github.com/gopro/gpr project, and is producing
26
  bitwise-identical output as compared with the Adobe DNG Converter
27
  implementation.
28
 */
29

30
#include "rawspeedconfig.h" // for HAVE_OPENMP
31
#include "decompressors/VC5Decompressor.h"
32
#include "common/Array2DRef.h"            // for Array2DRef
33
#include "common/Common.h"                // for clampBits, roundUpDivision
34
#include "common/Optional.h"              // for Optional
35
#include "common/Point.h"                 // for iPoint2D
36
#include "common/RawspeedException.h"     // for RawspeedException
37
#include "common/SimpleLUT.h"             // for SimpleLUT, SimpleLUT<>::va...
38
#include "decoders/RawDecoderException.h" // for ThrowRDE
39
#include "io/Endianness.h"                // for Endianness, Endianness::big
40
#include <cassert>                        // for assert
41
#include <cmath>                          // for pow
42
#include <initializer_list>               // for initializer_list
43
#include <limits>                         // for numeric_limits
44
#include <string>                         // for string
45
#include <utility>                        // for move
46

47
namespace {
48

49
// Definitions needed by table17.inc
50
// Taken from
51
// https://github.com/gopro/gpr/blob/a513701afce7b03173213a2f67dfd9dd28fa1868/source/lib/vc5_decoder/vlc.h
52
struct RLV {
53
  uint_fast8_t size; //!< Size of code word in bits
54
  uint32_t bits;     //!< Code word bits right justified
55
  uint16_t count;    //!< Run length
56
  uint16_t value;    //!< Run value (unsigned)
57
};
58
#define RLVTABLE(n)                                                            \
59
  struct {                                                                     \
60
    const uint32_t length;                                                     \
61
    const RLV entries[n];                                                      \
62
  } constexpr
63
#include "gopro/vc5/table17.inc"
64

65 2
constexpr int16_t decompand(int16_t val) {
66 2
  double c = val;
67
  // Invert companding curve
68 2
  c += (c * c * c * 768) / (255. * 255. * 255.);
69 2
  if (c > std::numeric_limits<int16_t>::max())
70 0
    return std::numeric_limits<int16_t>::max();
71 2
  if (c < std::numeric_limits<int16_t>::min())
72 0
    return std::numeric_limits<int16_t>::min();
73 2
  return c;
74
}
75

76
#ifndef NDEBUG
77 2
const int ignore = []() {
78 2
  for (const RLV& entry : table17.entries) {
79
    assert(((-decompand(entry.value)) == decompand(-int16_t(entry.value))) &&
80
           "negation of decompanded value is the same as decompanding of "
81
           "negated value");
82
  }
83 2
  return 0;
84
}();
85
#endif
86

87 2
const std::array<RLV, table17.length> decompandedTable17 = []() {
88
  std::array<RLV, table17.length> d;
89 2
  for (auto i = 0U; i < table17.length; i++) {
90 2
    d[i] = table17.entries[i];
91 2
    d[i].value = decompand(table17.entries[i].value);
92
  }
93 2
  return d;
94
}();
95

96
} // namespace
97

98
#define PRECISION_MIN 8
99
#define PRECISION_MAX 16
100

101
#define MARKER_BAND_END 1
102

103
namespace rawspeed {
104

105 1
void VC5Decompressor::Wavelet::setBandValid(const int band) {
106 1
  mDecodedBandMask |= (1 << band);
107
}
108

109 1
bool VC5Decompressor::Wavelet::isBandValid(const int band) const {
110 1
  return mDecodedBandMask & (1 << band);
111
}
112

113 1
bool VC5Decompressor::Wavelet::allBandsValid() const {
114 1
  return mDecodedBandMask == static_cast<uint32_t>((1 << numBands) - 1);
115
}
116

117
Array2DRef<const int16_t>
118 1
VC5Decompressor::Wavelet::bandAsArray2DRef(const unsigned int iBand) const {
119 1
  return {bands[iBand]->data.data(), width, height};
120
}
121

122
namespace {
123 1
const auto convolute = [](int row, int col, std::array<int, 4> muls,
124
                          const Array2DRef<const int16_t> high, auto lowGetter,
125
                          int DescaleShift = 0) {
126 1
  auto highCombined = muls[0] * high(row, col);
127 1
  auto lowsCombined = [muls, lowGetter]() {
128 1
    int lows = 0;
129 1
    for (int i = 0; i < 3; i++)
130 1
      lows += muls[1 + i] * lowGetter(i);
131 1
    return lows;
132 1
  }();
133
  // Round up 'lows' up
134 1
  lowsCombined += 4;
135
  // And finally 'average' them.
136 1
  auto lowsRounded = lowsCombined >> 3;
137 1
  auto total = highCombined + lowsRounded;
138
  // Descale it.
139 1
  total <<= DescaleShift;
140
  // And average it.
141 1
  total >>= 1;
142 1
  return total;
143
};
144

145
struct ConvolutionParams {
146
  struct First {
147
    static constexpr std::array<int, 4> mul_even = {+1, +11, -4, +1};
148
    static constexpr std::array<int, 4> mul_odd = {-1, +5, +4, -1};
149
    static constexpr int coord_shift = 0;
150
  };
151
  static constexpr First First{};
152

153
  struct Middle {
154
    static constexpr std::array<int, 4> mul_even = {+1, +1, +8, -1};
155
    static constexpr std::array<int, 4> mul_odd = {-1, -1, +8, +1};
156
    static constexpr int coord_shift = -1;
157
  };
158
  static constexpr Middle Middle{};
159

160
  struct Last {
161
    static constexpr std::array<int, 4> mul_even = {+1, -1, +4, +5};
162
    static constexpr std::array<int, 4> mul_odd = {-1, +1, -4, +11};
163
    static constexpr int coord_shift = -2;
164
  };
165
  static constexpr Last Last{};
166
};
167

168
constexpr std::array<int, 4> ConvolutionParams::First::mul_even;
169
constexpr std::array<int, 4> ConvolutionParams::First::mul_odd;
170

171
constexpr std::array<int, 4> ConvolutionParams::Middle::mul_even;
172
constexpr std::array<int, 4> ConvolutionParams::Middle::mul_odd;
173

174
constexpr std::array<int, 4> ConvolutionParams::Last::mul_even;
175
constexpr std::array<int, 4> ConvolutionParams::Last::mul_odd;
176

177
} // namespace
178

179 1
void VC5Decompressor::Wavelet::reconstructPass(
180
    const Array2DRef<int16_t> dst, const Array2DRef<const int16_t> high,
181
    const Array2DRef<const int16_t> low) const noexcept {
182 1
  auto process = [low, high, dst](auto segment, int row, int col) {
183 1
    auto lowGetter = [&row, &col, low](int delta) {
184 1
      return low(row + decltype(segment)::coord_shift + delta, col);
185
    };
186 1
    auto convolution = [&row, &col, high, lowGetter](std::array<int, 4> muls) {
187 1
      return convolute(row, col, muls, high, lowGetter, /*DescaleShift*/ 0);
188
    };
189

190 1
    int even = convolution(decltype(segment)::mul_even);
191 1
    int odd = convolution(decltype(segment)::mul_odd);
192

193 1
    dst(2 * row, col) = static_cast<int16_t>(even);
194 1
    dst(2 * row + 1, col) = static_cast<int16_t>(odd);
195 1
  };
196

197
  // Vertical reconstruction
198
#ifdef HAVE_OPENMP
199 1
#pragma omp for schedule(static)
200
#endif
201 1
  for (int row = 0; row < height; ++row) {
202 1
    if (row == 0) {
203
      // 1st row
204 1
      for (int col = 0; col < width; ++col)
205 1
        process(ConvolutionParams::First, row, col);
206 1
    } else if (row + 1 < height) {
207
      // middle rows
208 1
      for (int col = 0; col < width; ++col)
209 1
        process(ConvolutionParams::Middle, row, col);
210
    } else {
211
      // last row
212 1
      for (int col = 0; col < width; ++col)
213 1
        process(ConvolutionParams::Last, row, col);
214
    }
215
  }
216
}
217

218 1
void VC5Decompressor::Wavelet::combineLowHighPass(
219
    const Array2DRef<int16_t> dst, const Array2DRef<const int16_t> low,
220
    const Array2DRef<const int16_t> high, int descaleShift,
221
    bool clampUint = false) const noexcept {
222 1
  auto process = [low, high, descaleShift, clampUint, dst](auto segment,
223 1
                                                           int row, int col) {
224 1
    auto lowGetter = [&row, &col, low](int delta) {
225 1
      return low(row, col + decltype(segment)::coord_shift + delta);
226
    };
227 1
    auto convolution = [&row, &col, high, lowGetter,
228
                        descaleShift](std::array<int, 4> muls) {
229 1
      return convolute(row, col, muls, high, lowGetter, descaleShift);
230
    };
231

232 1
    int even = convolution(decltype(segment)::mul_even);
233 1
    int odd = convolution(decltype(segment)::mul_odd);
234

235 1
    if (clampUint) {
236 1
      even = clampBits(even, 14);
237 1
      odd = clampBits(odd, 14);
238
    }
239 1
    dst(row, 2 * col) = static_cast<int16_t>(even);
240 1
    dst(row, 2 * col + 1) = static_cast<int16_t>(odd);
241 1
  };
242

243
  // Horizontal reconstruction
244
#ifdef HAVE_OPENMP
245 1
#pragma omp for schedule(static)
246
#endif
247 1
  for (int row = 0; row < dst.height; ++row) {
248
    // First col
249 1
    int col = 0;
250 1
    process(ConvolutionParams::First, row, col);
251
    // middle cols
252 1
    for (col = 1; col + 1 < width; ++col) {
253 1
      process(ConvolutionParams::Middle, row, col);
254
    }
255
    // last col
256 1
    process(ConvolutionParams::Last, row, col);
257
  }
258
}
259

260 1
void VC5Decompressor::Wavelet::ReconstructableBand::processLow(
261
    const Wavelet& wavelet) noexcept {
262 1
  Array2DRef<int16_t> lowpass;
263
#ifdef HAVE_OPENMP
264 0
#pragma omp single copyprivate(lowpass)
265
#endif
266 1
  lowpass = Array2DRef<int16_t>::create(&lowpass_storage, wavelet.width,
267 1
                                        2 * wavelet.height);
268

269 1
  const Array2DRef<const int16_t> highlow = wavelet.bandAsArray2DRef(2);
270 1
  const Array2DRef<const int16_t> lowlow = wavelet.bandAsArray2DRef(0);
271

272
  // Reconstruct the "immediates", the actual low pass ...
273 1
  wavelet.reconstructPass(lowpass, highlow, lowlow);
274
}
275

276 1
void VC5Decompressor::Wavelet::ReconstructableBand::processHigh(
277
    const Wavelet& wavelet) noexcept {
278 1
  Array2DRef<int16_t> highpass;
279
#ifdef HAVE_OPENMP
280 0
#pragma omp single copyprivate(highpass)
281
#endif
282 1
  highpass = Array2DRef<int16_t>::create(&highpass_storage, wavelet.width,
283 1
                                         2 * wavelet.height);
284

285 1
  const Array2DRef<const int16_t> highhigh = wavelet.bandAsArray2DRef(3);
286 1
  const Array2DRef<const int16_t> lowhigh = wavelet.bandAsArray2DRef(1);
287

288 1
  wavelet.reconstructPass(highpass, highhigh, lowhigh);
289
}
290

291 1
void VC5Decompressor::Wavelet::ReconstructableBand::combine(
292
    const Wavelet& wavelet) noexcept {
293 1
  int16_t descaleShift = (wavelet.prescale == 2 ? 2 : 0);
294

295 1
  Array2DRef<int16_t> dest;
296
#ifdef HAVE_OPENMP
297 0
#pragma omp single copyprivate(dest)
298
#endif
299 1
  dest =
300 1
      Array2DRef<int16_t>::create(&data, 2 * wavelet.width, 2 * wavelet.height);
301

302 1
  const Array2DRef<int16_t> lowpass(lowpass_storage.data(), wavelet.width,
303 1
                                    2 * wavelet.height);
304 1
  const Array2DRef<int16_t> highpass(highpass_storage.data(), wavelet.width,
305 1
                                     2 * wavelet.height);
306

307
  // And finally, combine the low pass, and high pass.
308 1
  wavelet.combineLowHighPass(dest, lowpass, highpass, descaleShift, clampUint);
309
}
310

311 1
void VC5Decompressor::Wavelet::ReconstructableBand::decode(
312
    const Wavelet& wavelet) noexcept {
313 0
  assert(wavelet.allBandsValid());
314 0
  assert(data.empty());
315 1
  processLow(wavelet);
316 1
  processHigh(wavelet);
317 1
  combine(wavelet);
318
}
319

320 1
VC5Decompressor::VC5Decompressor(ByteStream bs, const RawImage& img)
321 1
    : mRaw(img), mBs(std::move(bs)) {
322 1
  if (!mRaw->dim.hasPositiveArea())
323 0
    ThrowRDE("Bad image dimensions.");
324

325 1
  if (mRaw->dim.x % mVC5.patternWidth != 0)
326 0
    ThrowRDE("Width %u is not a multiple of %u", mRaw->dim.x,
327
             mVC5.patternWidth);
328

329 1
  if (mRaw->dim.y % mVC5.patternHeight != 0)
330 0
    ThrowRDE("Height %u is not a multiple of %u", mRaw->dim.y,
331
             mVC5.patternHeight);
332

333
  // Initialize wavelet sizes.
334 1
  for (Channel& channel : channels) {
335 1
    channel.width = mRaw->dim.x / mVC5.patternWidth;
336 1
    channel.height = mRaw->dim.y / mVC5.patternHeight;
337

338 1
    uint16_t waveletWidth = channel.width;
339 1
    uint16_t waveletHeight = channel.height;
340 1
    for (Wavelet& wavelet : channel.wavelets) {
341
      // Pad dimensions as necessary and divide them by two for the next wavelet
342 1
      for (auto* dimension : {&waveletWidth, &waveletHeight})
343 1
        *dimension = roundUpDivision(*dimension, 2);
344 1
      wavelet.width = waveletWidth;
345 1
      wavelet.height = waveletHeight;
346
    }
347
  }
348

349 0
  if (img->whitePoint <= 0 || img->whitePoint > int(((1U << 16U) - 1U)))
350 0
    ThrowRDE("Bad white level %i", img->whitePoint);
351

352 1
  outputBits = 0;
353 1
  for (int wp = img->whitePoint; wp != 0; wp >>= 1)
354 1
    ++outputBits;
355 0
  assert(outputBits <= 16);
356

357 1
  parseVC5();
358
}
359

360 1
void VC5Decompressor::initVC5LogTable() {
361 1
  mVC5LogTable = decltype(mVC5LogTable)(
362 1
      [outputBits = outputBits](unsigned i, unsigned tableSize) {
363
        // The vanilla "inverse log" curve for decoding.
364 1
        auto normalizedCurve = [](auto normalizedI) {
365 1
          return (std::pow(113.0, normalizedI) - 1) / 112.0;
366
        };
367

368 1
        auto normalizeI = [tableSize](auto x) { return x / (tableSize - 1.0); };
369 1
        auto denormalizeY = [maxVal = std::numeric_limits<uint16_t>::max()](
370 1
                                auto y) { return maxVal * y; };
371
        // Adjust for output whitelevel bitdepth.
372 1
        auto rescaleY = [outputBits](auto y) {
373 1
          auto scale = 16 - outputBits;
374 1
          return y >> scale;
375 1
        };
376

377 1
        const auto naiveY = denormalizeY(normalizedCurve(normalizeI(i)));
378 1
        const auto intY = static_cast<unsigned int>(naiveY);
379 1
        const auto rescaledY = rescaleY(intY);
380 1
        return rescaledY;
381 1
      });
382
}
383

384 1
void VC5Decompressor::parseVC5() {
385 1
  mBs.setByteOrder(Endianness::big);
386

387 0
  assert(mRaw->dim.x > 0);
388 0
  assert(mRaw->dim.y > 0);
389

390
  // All VC-5 data must start with "VC-%" (0x56432d35)
391 1
  if (mBs.getU32() != 0x56432d35)
392 0
    ThrowRDE("not a valid VC-5 datablock");
393

394 1
  bool done = false;
395 1
  while (!done) {
396 1
    auto tag = static_cast<VC5Tag>(mBs.getU16());
397 1
    uint16_t val = mBs.getU16();
398

399 1
    bool optional = matches(tag, VC5Tag::Optional);
400 1
    if (optional)
401 1
      tag = -tag;
402

403 1
    switch (tag) {
404 1
    case VC5Tag::ChannelCount:
405 1
      if (val != numChannels)
406 0
        ThrowRDE("Bad channel count %u, expected %u", val, numChannels);
407 1
      break;
408 1
    case VC5Tag::ImageWidth:
409 1
      if (val != mRaw->dim.x)
410 0
        ThrowRDE("Image width mismatch: %u vs %u", val, mRaw->dim.x);
411 1
      break;
412 1
    case VC5Tag::ImageHeight:
413 1
      if (val != mRaw->dim.y)
414 0
        ThrowRDE("Image height mismatch: %u vs %u", val, mRaw->dim.y);
415 1
      break;
416 1
    case VC5Tag::LowpassPrecision:
417 1
      if (val < PRECISION_MIN || val > PRECISION_MAX)
418 0
        ThrowRDE("Invalid precision %i", val);
419 1
      mVC5.lowpassPrecision = val;
420 1
      break;
421 1
    case VC5Tag::ChannelNumber:
422 1
      if (val >= numChannels)
423 0
        ThrowRDE("Bad channel number (%u)", val);
424 1
      mVC5.iChannel = val;
425 1
      break;
426 1
    case VC5Tag::ImageFormat:
427 1
      if (val != mVC5.imgFormat)
428 0
        ThrowRDE("Image format %i is not 4(RAW)", val);
429 1
      break;
430 0
    case VC5Tag::SubbandCount:
431 0
      if (val != numSubbands)
432 0
        ThrowRDE("Unexpected subband count %u, expected %u", val, numSubbands);
433 0
      break;
434 1
    case VC5Tag::MaxBitsPerComponent:
435 1
      if (val != VC5_LOG_TABLE_BITWIDTH) {
436 0
        ThrowRDE("Bad bits per componend %u, not %u", val,
437
                 VC5_LOG_TABLE_BITWIDTH);
438
      }
439 1
      break;
440 1
    case VC5Tag::PatternWidth:
441 1
      if (val != mVC5.patternWidth)
442 0
        ThrowRDE("Bad pattern width %u, not %u", val, mVC5.patternWidth);
443 1
      break;
444 1
    case VC5Tag::PatternHeight:
445 1
      if (val != mVC5.patternHeight)
446 0
        ThrowRDE("Bad pattern height %u, not %u", val, mVC5.patternHeight);
447 1
      break;
448 1
    case VC5Tag::SubbandNumber:
449 1
      if (val >= numSubbands)
450 0
        ThrowRDE("Bad subband number %u", val);
451 1
      mVC5.iSubband = val;
452 1
      break;
453 1
    case VC5Tag::Quantization:
454 1
      mVC5.quantization = static_cast<int16_t>(val);
455 1
      break;
456 1
    case VC5Tag::ComponentsPerSample:
457 1
      if (val != mVC5.cps)
458 0
        ThrowRDE("Bad component per sample count %u, not %u", val, mVC5.cps);
459 1
      break;
460 1
    case VC5Tag::PrescaleShift:
461
      // FIXME: something is wrong. We get this before VC5Tag::ChannelNumber.
462
      // Defaulting to 'mVC5.iChannel=0' seems to work *for existing samples*.
463 1
      for (int iWavelet = 0; iWavelet < numWaveletLevels; ++iWavelet) {
464 1
        auto& channel = channels[mVC5.iChannel];
465 1
        auto& wavelet = channel.wavelets[iWavelet];
466 1
        wavelet.prescale =
467 1
            extractHighBits(val, 2 * iWavelet, /*effectiveBitwidth=*/14) & 0x03;
468
      }
469 1
      break;
470 1
    default: { // A chunk.
471 1
      unsigned int chunkSize = 0;
472 1
      if (matches(tag, VC5Tag::LARGE_CHUNK)) {
473 1
        chunkSize = static_cast<unsigned int>(
474 1
            ((static_cast<std::underlying_type<VC5Tag>::type>(tag) & 0xff)
475 1
             << 16) |
476 1
            (val & 0xffff));
477 1
      } else if (matches(tag, VC5Tag::SMALL_CHUNK)) {
478 1
        chunkSize = (val & 0xffff);
479
      }
480

481 1
      if (is(tag, VC5Tag::LargeCodeblock)) {
482 1
        parseLargeCodeblock(mBs.getStream(chunkSize, 4));
483 1
        break;
484
      }
485

486
      // And finally, we got here if we didn't handle this tag/maybe-chunk.
487

488
      // Magic, all the other 'large' chunks are actually optional,
489
      // and don't specify any chunk bytes-to-be-skipped.
490 1
      if (matches(tag, VC5Tag::LARGE_CHUNK)) {
491 1
        optional = true;
492 1
        chunkSize = 0;
493
      }
494

495 1
      if (!optional) {
496 0
        ThrowRDE("Unknown (unhandled) non-optional Tag 0x%04hx",
497
                 static_cast<std::underlying_type<VC5Tag>::type>(tag));
498
      }
499

500 1
      if (chunkSize)
501 1
        mBs.skipBytes(chunkSize, 4);
502

503 1
      break;
504
    }
505
    }
506

507 1
    done = true;
508 1
    for (int iChannel = 0; iChannel < numChannels && done; ++iChannel) {
509 1
      Wavelet& wavelet = channels[iChannel].wavelets[0];
510 1
      if (!wavelet.allBandsValid())
511 1
        done = false;
512
    }
513
  }
514
}
515

516 1
VC5Decompressor::Wavelet::LowPassBand::LowPassBand(const Wavelet& wavelet,
517
                                                   ByteStream bs_,
518 1
                                                   uint16_t lowpassPrecision_)
519 1
    : AbstractDecodeableBand(std::move(bs_)),
520 1
      lowpassPrecision(lowpassPrecision_) {
521
  // Low-pass band is a uncompressed version of the image, hugely downscaled.
522
  // It consists of width * height pixels, `lowpassPrecision` each.
523
  // We can easily check that we have sufficient amount of bits to decode it.
524 1
  const auto waveletArea = iPoint2D(wavelet.width, wavelet.height).area();
525 1
  const auto bitsTotal = waveletArea * lowpassPrecision;
526 1
  const auto bytesTotal = roundUpDivision(bitsTotal, 8);
527 1
  bs = bs.getStream(bytesTotal); // And clamp the size while we are at it.
528
}
529

530 1
void VC5Decompressor::Wavelet::LowPassBand::decode(const Wavelet& wavelet) {
531
  const auto dst =
532 1
      Array2DRef<int16_t>::create(&data, wavelet.width, wavelet.height);
533

534 1
  BitPumpMSB bits(bs);
535 1
  for (auto row = 0; row < dst.height; ++row) {
536 1
    for (auto col = 0; col < dst.width; ++col)
537 1
      dst(row, col) = static_cast<int16_t>(bits.getBits(lowpassPrecision));
538
  }
539
}
540

541 1
void VC5Decompressor::Wavelet::HighPassBand::decode(const Wavelet& wavelet) {
542 1
  auto dequantize = [quant = quant](int16_t val) -> int16_t {
543 1
    return val * quant;
544 1
  };
545

546 1
  Array2DRef<int16_t>::create(&data, wavelet.width, wavelet.height);
547

548 1
  BitPumpMSB bits(bs);
549
  // decode highpass band
550 1
  int pixelValue = 0;
551 1
  unsigned int count = 0;
552 1
  int nPixels = wavelet.width * wavelet.height;
553 1
  for (int iPixel = 0; iPixel < nPixels;) {
554 1
    getRLV(&bits, &pixelValue, &count);
555 1
    for (; count > 0; --count) {
556 1
      if (iPixel >= nPixels)
557 0
        ThrowRDE("Buffer overflow");
558 1
      data[iPixel] = dequantize(pixelValue);
559 1
      ++iPixel;
560
    }
561
  }
562 1
  getRLV(&bits, &pixelValue, &count);
563
  static_assert(decompand(MARKER_BAND_END) == MARKER_BAND_END, "passthrough");
564 1
  if (pixelValue != MARKER_BAND_END || count != 0)
565 0
    ThrowRDE("EndOfBand marker not found");
566
}
567

568 1
void VC5Decompressor::parseLargeCodeblock(const ByteStream& bs) {
569 1
  static const auto subband_wavelet_index = []() {
570
    std::array<int, numSubbands> wavelets;
571 1
    int wavelet = 0;
572 1
    for (auto i = wavelets.size() - 1; i > 0;) {
573 1
      for (auto t = 0; t < numWaveletLevels; t++) {
574 1
        wavelets[i] = wavelet;
575 1
        i--;
576
      }
577 1
      if (i > 0)
578 1
        wavelet++;
579
    }
580 1
    wavelets.front() = wavelet;
581 1
    return wavelets;
582 1
  }();
583 1
  static const auto subband_band_index = []() {
584
    std::array<int, numSubbands> bands;
585 1
    bands.front() = 0;
586 1
    for (auto i = 1U; i < bands.size();) {
587 1
      for (int t = 1; t <= numWaveletLevels;) {
588 1
        bands[i] = t;
589 1
        t++;
590 1
        i++;
591
      }
592
    }
593 1
    return bands;
594 1
  }();
595

596 1
  if (!mVC5.iSubband.hasValue())
597 0
    ThrowRDE("Did not see VC5Tag::SubbandNumber yet");
598

599 1
  const int idx = subband_wavelet_index[mVC5.iSubband.getValue()];
600 1
  const int band = subband_band_index[mVC5.iSubband.getValue()];
601

602 1
  auto& wavelets = channels[mVC5.iChannel].wavelets;
603

604 1
  Wavelet& wavelet = wavelets[idx];
605 1
  if (wavelet.isBandValid(band)) {
606 0
    ThrowRDE("Band %u for wavelet %u on channel %u was already seen", band, idx,
607
             mVC5.iChannel);
608
  }
609

610 1
  std::unique_ptr<Wavelet::AbstractBand>& dstBand = wavelet.bands[band];
611 1
  if (mVC5.iSubband.getValue() == 0) {
612 0
    assert(band == 0);
613
    // low-pass band, only one, for the smallest wavelet, per channel per image
614 1
    if (!mVC5.lowpassPrecision.hasValue())
615 0
      ThrowRDE("Did not see VC5Tag::LowpassPrecision yet");
616 1
    dstBand = std::make_unique<Wavelet::LowPassBand>(
617 1
        wavelet, bs, mVC5.lowpassPrecision.getValue());
618 1
    mVC5.lowpassPrecision.reset();
619
  } else {
620 1
    if (!mVC5.quantization.hasValue())
621 0
      ThrowRDE("Did not see VC5Tag::Quantization yet");
622 1
    dstBand = std::make_unique<Wavelet::HighPassBand>(
623 1
        bs, mVC5.quantization.getValue());
624 1
    mVC5.quantization.reset();
625
  }
626 1
  wavelet.setBandValid(band);
627

628
  // If this wavelet is fully specified, mark the low-pass band of the
629
  // next lower wavelet as specified.
630 1
  if (idx > 0 && wavelet.allBandsValid()) {
631 1
    Wavelet& nextWavelet = wavelets[idx - 1];
632 0
    assert(!nextWavelet.isBandValid(0));
633 1
    nextWavelet.bands[0] = std::make_unique<Wavelet::ReconstructableBand>();
634 1
    nextWavelet.setBandValid(0);
635
  }
636

637 1
  mVC5.iSubband.reset();
638
}
639

640 1
void VC5Decompressor::prepareBandDecodingPlan() {
641 0
  assert(allDecodeableBands.empty());
642 1
  allDecodeableBands.reserve(numSubbandsTotal);
643
  // All the high-pass bands for all wavelets,
644
  // in this specific order of decreasing worksize.
645 1
  for (int waveletLevel = 0; waveletLevel < numWaveletLevels; waveletLevel++) {
646 1
    for (auto channelId = 0; channelId < numChannels; channelId++) {
647 1
      for (int bandId = 1; bandId <= numHighPassBands; bandId++) {
648 1
        auto& channel = channels[channelId];
649 1
        auto& wavelet = channel.wavelets[waveletLevel];
650 1
        auto* band = wavelet.bands[bandId].get();
651
        auto* decodeableHighPassBand =
652 0
            dynamic_cast<Wavelet::HighPassBand*>(band);
653 1
        allDecodeableBands.emplace_back(decodeableHighPassBand, wavelet);
654
      }
655
    }
656
  }
657
  // The low-pass bands at the end. I'm guessing they should be fast to
658
  // decode.
659 1
  for (Channel& channel : channels) {
660
    // Low-pass band of the smallest wavelet.
661 1
    Wavelet& smallestWavelet = channel.wavelets.back();
662
    auto* decodeableLowPassBand =
663 0
        dynamic_cast<Wavelet::LowPassBand*>(smallestWavelet.bands[0].get());
664 1
    allDecodeableBands.emplace_back(decodeableLowPassBand, smallestWavelet);
665
  }
666 0
  assert(allDecodeableBands.size() == numSubbandsTotal);
667
}
668

669 1
void VC5Decompressor::prepareBandReconstruction() {
670 0
  assert(reconstructionSteps.empty());
671 1
  reconstructionSteps.reserve(numLowPassBandsTotal);
672
  // For every channel, recursively reconstruct the low-pass bands.
673 1
  for (auto& channel : channels) {
674
    // Reconstruct the intermediate lowpass bands.
675 1
    for (int waveletLevel = numWaveletLevels - 1; waveletLevel > 0;
676
         waveletLevel--) {
677 1
      Wavelet* wavelet = &(channel.wavelets[waveletLevel]);
678 1
      Wavelet& nextWavelet = channel.wavelets[waveletLevel - 1];
679

680 0
      auto* band = dynamic_cast<Wavelet::ReconstructableBand*>(
681 1
          nextWavelet.bands[0].get());
682 1
      reconstructionSteps.emplace_back(wavelet, band);
683
    }
684
    // Finally, reconstruct the final lowpass band.
685 1
    Wavelet* wavelet = &(channel.wavelets.front());
686 1
    reconstructionSteps.emplace_back(wavelet, &(channel.band));
687
  }
688 0
  assert(reconstructionSteps.size() == numLowPassBandsTotal);
689
}
690

691 1
void VC5Decompressor::prepareDecodingPlan() {
692 1
  prepareBandDecodingPlan();
693 1
  prepareBandReconstruction();
694
}
695

696 1
void VC5Decompressor::decodeThread(bool* exceptionThrown) const noexcept {
697
  // Decode all the existing bands. May fail.
698 1
  decodeBands(exceptionThrown);
699

700
  // Proceed only if decoding did not fail.
701 1
  if (*exceptionThrown)
702 0
    return;
703

704
  // And now, reconstruct the low-pass bands.
705 1
  reconstructLowpassBands();
706

707
  // And finally!
708 1
  combineFinalLowpassBands();
709
}
710

711 1
void VC5Decompressor::decode(unsigned int offsetX, unsigned int offsetY,
712
                             unsigned int width, unsigned int height) {
713 0
  if (offsetX || offsetY || mRaw->dim != iPoint2D(width, height))
714 0
    ThrowRDE("VC5Decompressor expects to fill the whole image, not some tile.");
715

716 1
  initVC5LogTable();
717

718 1
  prepareDecodingPlan();
719

720 1
  bool exceptionThrown = false;
721
#ifdef HAVE_OPENMP
722 1
#pragma omp parallel default(none) shared(exceptionThrown)                     \
723 1
    num_threads(rawspeed_get_number_of_processor_cores())
724
#endif
725
  decodeThread(&exceptionThrown);
726

727 1
  std::string firstErr;
728 1
  if (mRaw->isTooManyErrors(1, &firstErr)) {
729 0
    assert(exceptionThrown);
730 0
    ThrowRDE("Too many errors encountered. Giving up. First Error:\n%s",
731
             firstErr.c_str());
732
  } else {
733 0
    assert(!exceptionThrown);
734
  }
735
}
736

737 1
void VC5Decompressor::decodeBands(bool* exceptionThrown) const noexcept {
738
#ifdef HAVE_OPENMP
739 1
#pragma omp for schedule(dynamic, 1)
740
#endif
741 1
  for (auto decodeableBand = allDecodeableBands.begin();
742 1
       decodeableBand < allDecodeableBands.end(); ++decodeableBand) {
743
    try {
744 1
      decodeableBand->band->decode(decodeableBand->wavelet);
745 0
    } catch (RawspeedException& err) {
746
      // Propagate the exception out of OpenMP magic.
747 0
      mRaw->setError(err.what());
748
#ifdef HAVE_OPENMP
749 0
#pragma omp atomic write
750
#endif
751
      *exceptionThrown = true;
752
#ifdef HAVE_OPENMP
753 0
#pragma omp cancel for
754
#endif
755
    }
756
  }
757
}
758

759 1
void VC5Decompressor::reconstructLowpassBands() const noexcept {
760 1
  for (const ReconstructionStep& step : reconstructionSteps) {
761 1
    step.band.decode(step.wavelet);
762

763
#ifdef HAVE_OPENMP
764
#pragma omp single nowait
765
#endif
766 1
    step.wavelet.clear(); // we no longer need it.
767
  }
768
}
769

770 1
void VC5Decompressor::combineFinalLowpassBands() const noexcept {
771 1
  const Array2DRef<uint16_t> out(mRaw->getU16DataAsUncroppedArray2DRef());
772

773 1
  const int width = out.width / 2;
774 1
  const int height = out.height / 2;
775

776
  const Array2DRef<const int16_t> lowbands0 = Array2DRef<const int16_t>(
777 1
      channels[0].band.data.data(), channels[0].width, channels[0].height);
778
  const Array2DRef<const int16_t> lowbands1 = Array2DRef<const int16_t>(
779 1
      channels[1].band.data.data(), channels[1].width, channels[1].height);
780
  const Array2DRef<const int16_t> lowbands2 = Array2DRef<const int16_t>(
781 1
      channels[2].band.data.data(), channels[2].width, channels[2].height);
782
  const Array2DRef<const int16_t> lowbands3 = Array2DRef<const int16_t>(
783 1
      channels[3].band.data.data(), channels[3].width, channels[3].height);
784

785
  // Convert to RGGB output
786
#ifdef HAVE_OPENMP
787 1
#pragma omp for schedule(static) collapse(2)
788
#endif
789
  for (int row = 0; row < height; ++row) {
790
    for (int col = 0; col < width; ++col) {
791 1
      const int mid = 2048;
792

793 1
      int gs = lowbands0(row, col);
794 1
      int rg = lowbands1(row, col) - mid;
795 1
      int bg = lowbands2(row, col) - mid;
796 1
      int gd = lowbands3(row, col) - mid;
797

798 1
      int r = gs + 2 * rg;
799 1
      int b = gs + 2 * bg;
800 1
      int g1 = gs + gd;
801 1
      int g2 = gs - gd;
802

803 1
      out(2 * row + 0, 2 * col + 0) = static_cast<uint16_t>(mVC5LogTable[r]);
804 1
      out(2 * row + 0, 2 * col + 1) = static_cast<uint16_t>(mVC5LogTable[g1]);
805 1
      out(2 * row + 1, 2 * col + 0) = static_cast<uint16_t>(mVC5LogTable[g2]);
806 1
      out(2 * row + 1, 2 * col + 1) = static_cast<uint16_t>(mVC5LogTable[b]);
807
    }
808
  }
809
}
810

811 1
inline void VC5Decompressor::getRLV(BitPumpMSB* bits, int* value,
812
                                    unsigned int* count) {
813
  unsigned int iTab;
814

815
  static constexpr auto maxBits = 1 + table17.entries[table17.length - 1].size;
816

817
  // Ensure the maximum number of bits are cached to make peekBits() as fast as
818
  // possible.
819 1
  bits->fill(maxBits);
820 1
  for (iTab = 0; iTab < table17.length; ++iTab) {
821 1
    if (decompandedTable17[iTab].bits ==
822 1
        bits->peekBitsNoFill(decompandedTable17[iTab].size))
823 1
      break;
824
  }
825 1
  if (iTab >= table17.length)
826 0
    ThrowRDE("Code not found in codebook");
827

828 1
  bits->skipBitsNoFill(decompandedTable17[iTab].size);
829 1
  *value = decompandedTable17[iTab].value;
830 1
  *count = decompandedTable17[iTab].count;
831 1
  if (*value != 0) {
832 1
    if (bits->getBitsNoFill(1))
833 1
      *value = -(*value);
834
  }
835
}
836

837
} // namespace rawspeed

Read our documentation on viewing source code .

Loading