OpenJPH
Open-source implementation of JPEG2000 Part-15
ojph_colour_avx2.cpp
Go to the documentation of this file.
1//***************************************************************************/
2// This software is released under the 2-Clause BSD license, included
3// below.
4//
5// Copyright (c) 2019, Aous Naman
6// Copyright (c) 2019, Kakadu Software Pty Ltd, Australia
7// Copyright (c) 2019, The University of New South Wales, Australia
8//
9// Redistribution and use in source and binary forms, with or without
10// modification, are permitted provided that the following conditions are
11// met:
12//
13// 1. Redistributions of source code must retain the above copyright
14// notice, this list of conditions and the following disclaimer.
15//
16// 2. Redistributions in binary form must reproduce the above copyright
17// notice, this list of conditions and the following disclaimer in the
18// documentation and/or other materials provided with the distribution.
19//
20// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
21// IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
22// TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
23// PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
24// HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
25// SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
26// TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
27// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
28// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
29// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
30// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
31//***************************************************************************/
32// This file is part of the OpenJPH software implementation.
33// File: ojph_colour_avx2.cpp
34// Author: Aous Naman
35// Date: 11 October 2019
36//***************************************************************************/
37
38#include <climits>
39#include <cmath>
40
41#include "ojph_defs.h"
42#include "ojph_arch.h"
43#include "ojph_mem.h"
44#include "ojph_colour.h"
45
46#include <immintrin.h>
47
48namespace ojph {
49 namespace local {
50
52 // https://github.com/seung-lab/dijkstra3d/blob/master/libdivide.h
53 static inline
54 __m256i avx2_mm256_srai_epi64(__m256i a, int amt, __m256i m)
55 {
56 // note than m must be obtained using
57 // __m256i m = _mm256_set1_epi64x(1ULL << (63 - amt));
58 __m256i x = _mm256_srli_epi64(a, amt);
59 x = _mm256_xor_si256(x, m);
60 __m256i result = _mm256_sub_epi64(x, m);
61 return result;
62 }
63
65 void avx2_rev_convert(const line_buf *src_line,
66 const ui32 src_line_offset,
67 line_buf *dst_line,
68 const ui32 dst_line_offset,
69 si64 shift, ui32 width)
70 {
71 if (src_line->flags & line_buf::LFT_32BIT)
72 {
73 if (dst_line->flags & line_buf::LFT_32BIT)
74 {
75 const si32 *sp = src_line->i32 + src_line_offset;
76 si32 *dp = dst_line->i32 + dst_line_offset;
77 __m256i sh = _mm256_set1_epi32((si32)shift);
78 for (int i = (width + 7) >> 3; i > 0; --i, sp+=8, dp+=8)
79 {
80 __m256i s = _mm256_loadu_si256((__m256i*)sp);
81 s = _mm256_add_epi32(s, sh);
82 _mm256_storeu_si256((__m256i*)dp, s);
83 }
84 }
85 else
86 {
87 const si32 *sp = src_line->i32 + src_line_offset;
88 si64 *dp = dst_line->i64 + dst_line_offset;
89 __m256i sh = _mm256_set1_epi64x(shift);
90 for (int i = (width + 7) >> 3; i > 0; --i, sp+=8, dp+=8)
91 {
92 __m256i s, t;
93 s = _mm256_loadu_si256((__m256i*)sp);
94
95 t = _mm256_cvtepi32_epi64(_mm256_extracti128_si256(s, 0));
96 t = _mm256_add_epi64(t, sh);
97 _mm256_storeu_si256((__m256i*)dp, t);
98
99 t = _mm256_cvtepi32_epi64(_mm256_extracti128_si256(s, 1));
100 t = _mm256_add_epi64(t, sh);
101 _mm256_storeu_si256((__m256i*)dp + 1, t);
102 }
103 }
104 }
105 else
106 {
107 assert(src_line->flags | line_buf::LFT_64BIT);
108 assert(dst_line->flags | line_buf::LFT_32BIT);
109 const si64 *sp = src_line->i64 + src_line_offset;
110 si32 *dp = dst_line->i32 + dst_line_offset;
111 __m256i low_bits = _mm256_set_epi64x(0, (si64)ULLONG_MAX,
112 0, (si64)ULLONG_MAX);
113 __m256i sh = _mm256_set1_epi64x(shift);
114 for (int i = (width + 7) >> 3; i > 0; --i, sp+=8, dp+=8)
115 {
116 __m256i s, t;
117 s = _mm256_loadu_si256((__m256i*)sp);
118 s = _mm256_add_epi64(s, sh);
119
120 t = _mm256_shuffle_epi32(s, _MM_SHUFFLE(0, 0, 2, 0));
121 t = _mm256_and_si256(low_bits, t);
122
123 s = _mm256_loadu_si256((__m256i*)sp + 1);
124 s = _mm256_add_epi64(s, sh);
125
126 s = _mm256_shuffle_epi32(s, _MM_SHUFFLE(2, 0, 0, 0));
127 s = _mm256_andnot_si256(low_bits, s);
128
129 t = _mm256_or_si256(s, t);
130 t = _mm256_permute4x64_epi64(t, _MM_SHUFFLE(3, 1, 2, 0));
131 _mm256_storeu_si256((__m256i*)dp, t);
132 }
133 }
134 }
135
138 const ui32 src_line_offset,
139 line_buf *dst_line,
140 const ui32 dst_line_offset,
141 si64 shift, ui32 width)
142 {
143 if (src_line->flags & line_buf::LFT_32BIT)
144 {
145 if (dst_line->flags & line_buf::LFT_32BIT)
146 {
147 const si32 *sp = src_line->i32 + src_line_offset;
148 si32 *dp = dst_line->i32 + dst_line_offset;
149 __m256i sh = _mm256_set1_epi32((si32)(-shift));
150 __m256i zero = _mm256_setzero_si256();
151 for (int i = (width + 7) >> 3; i > 0; --i, sp += 8, dp += 8)
152 {
153 __m256i s = _mm256_loadu_si256((__m256i*)sp);
154 __m256i c = _mm256_cmpgt_epi32(zero, s); // 0xFFFFFFFF for -ve val
155 __m256i v_m_sh = _mm256_sub_epi32(sh, s); // - shift - value
156 v_m_sh = _mm256_and_si256(c, v_m_sh); // keep only -shift-val
157 s = _mm256_andnot_si256(c, s); // keep only +ve or 0
158 s = _mm256_or_si256(s, v_m_sh); // combine
159 _mm256_storeu_si256((__m256i*)dp, s);
160 }
161 }
162 else
163 {
164 const si32 *sp = src_line->i32 + src_line_offset;
165 si64 *dp = dst_line->i64 + dst_line_offset;
166 __m256i sh = _mm256_set1_epi64x(-shift);
167 __m256i zero = _mm256_setzero_si256();
168 for (int i = (width + 7) >> 3; i > 0; --i, sp += 8, dp += 8)
169 {
170 __m256i s, t, u0, u1, c, v_m_sh;
171 s = _mm256_loadu_si256((__m256i*)sp);
172
173 t = _mm256_cmpgt_epi32(zero, s); // find -ve 32bit -1
174 u0 = _mm256_unpacklo_epi32(s, t); // correct 64bit data
175 c = _mm256_unpacklo_epi32(t, t); // 64bit -1 for -ve value
176
177 v_m_sh = _mm256_sub_epi64(sh, u0); // - shift - value
178 v_m_sh = _mm256_and_si256(c, v_m_sh); // keep only - shift - value
179 u0 = _mm256_andnot_si256(c, u0); // keep only +ve or 0
180 u0 = _mm256_or_si256(u0, v_m_sh); // combine
181
182 u1 = _mm256_unpackhi_epi32(s, t); // correct 64bit data
183 c = _mm256_unpackhi_epi32(t, t); // 64bit -1 for -ve value
184
185 v_m_sh = _mm256_sub_epi64(sh, u1); // - shift - value
186 v_m_sh = _mm256_and_si256(c, v_m_sh); // keep only - shift - value
187 u1 = _mm256_andnot_si256(c, u1); // keep only +ve or 0
188 u1 = _mm256_or_si256(u1, v_m_sh); // combine
189
190 t = _mm256_permute2x128_si256(u0, u1, (2 << 4) | 0);
191 _mm256_storeu_si256((__m256i*)dp, t);
192
193 t = _mm256_permute2x128_si256(u0, u1, (3 << 4) | 1);
194 _mm256_storeu_si256((__m256i*)dp + 1, t);
195 }
196 }
197 }
198 else
199 {
200 assert(src_line->flags | line_buf::LFT_64BIT);
201 assert(dst_line->flags | line_buf::LFT_32BIT);
202 const si64 *sp = src_line->i64 + src_line_offset;
203 si32 *dp = dst_line->i32 + dst_line_offset;
204 __m256i sh = _mm256_set1_epi64x(-shift);
205 __m256i zero = _mm256_setzero_si256();
206 __m256i half_mask = _mm256_set_epi64x(0, (si64)ULLONG_MAX,
207 0, (si64)ULLONG_MAX);
208 for (int i = (width + 7) >> 3; i > 0; --i, sp += 8, dp += 8)
209 {
210 // s for source, t for target, p for positive, n for negative,
211 // m for mask, and tm for temp
212 __m256i s, t, p, n, m, tm;
213 s = _mm256_loadu_si256((__m256i*)sp);
214
215 m = _mm256_cmpgt_epi64(zero, s); // 64b -1 for -ve value
216 tm = _mm256_sub_epi64(sh, s); // - shift - value
217 n = _mm256_and_si256(m, tm); // -ve
218 p = _mm256_andnot_si256(m, s); // +ve
219 tm = _mm256_or_si256(n, p);
220 tm = _mm256_shuffle_epi32(tm, _MM_SHUFFLE(0, 0, 2, 0));
221 t = _mm256_and_si256(half_mask, tm);
222
223 s = _mm256_loadu_si256((__m256i*)sp + 1);
224 m = _mm256_cmpgt_epi64(zero, s); // 64b -1 for -ve value
225 tm = _mm256_sub_epi64(sh, s); // - shift - value
226 n = _mm256_and_si256(m, tm); // -ve
227 p = _mm256_andnot_si256(m, s); // +ve
228 tm = _mm256_or_si256(n, p);
229 tm = _mm256_shuffle_epi32(tm, _MM_SHUFFLE(2, 0, 0, 0));
230 tm = _mm256_andnot_si256(half_mask, tm);
231
232 t = _mm256_or_si256(t, tm);
233 t = _mm256_permute4x64_epi64(t, _MM_SHUFFLE(3, 1, 2, 0));
234 _mm256_storeu_si256((__m256i*)dp, t);
235 }
236 }
237 }
238
241 const line_buf *g,
242 const line_buf *b,
243 line_buf *y, line_buf *cb, line_buf *cr,
244 ui32 repeat)
245 {
246 assert((y->flags & line_buf::LFT_REVERSIBLE) &&
252
253 if (y->flags & line_buf::LFT_32BIT)
254 {
255 assert((y->flags & line_buf::LFT_32BIT) &&
256 (cb->flags & line_buf::LFT_32BIT) &&
257 (cr->flags & line_buf::LFT_32BIT) &&
258 (r->flags & line_buf::LFT_32BIT) &&
259 (g->flags & line_buf::LFT_32BIT) &&
260 (b->flags & line_buf::LFT_32BIT));
261 const si32 *rp = r->i32, * gp = g->i32, * bp = b->i32;
262 si32 *yp = y->i32, * cbp = cb->i32, * crp = cr->i32;
263 for (int i = (repeat + 7) >> 3; i > 0; --i)
264 {
265 __m256i mr = _mm256_load_si256((__m256i*)rp);
266 __m256i mg = _mm256_load_si256((__m256i*)gp);
267 __m256i mb = _mm256_load_si256((__m256i*)bp);
268 __m256i t = _mm256_add_epi32(mr, mb);
269 t = _mm256_add_epi32(t, _mm256_slli_epi32(mg, 1));
270 _mm256_store_si256((__m256i*)yp, _mm256_srai_epi32(t, 2));
271 t = _mm256_sub_epi32(mb, mg);
272 _mm256_store_si256((__m256i*)cbp, t);
273 t = _mm256_sub_epi32(mr, mg);
274 _mm256_store_si256((__m256i*)crp, t);
275
276 rp += 8; gp += 8; bp += 8;
277 yp += 8; cbp += 8; crp += 8;
278 }
279 }
280 else
281 {
282 assert((y->flags & line_buf::LFT_64BIT) &&
283 (cb->flags & line_buf::LFT_64BIT) &&
284 (cr->flags & line_buf::LFT_64BIT) &&
285 (r->flags & line_buf::LFT_32BIT) &&
286 (g->flags & line_buf::LFT_32BIT) &&
288 __m256i v2 = _mm256_set1_epi64x(1ULL << (63 - 2));
289 const si32 *rp = r->i32, *gp = g->i32, *bp = b->i32;
290 si64 *yp = y->i64, *cbp = cb->i64, *crp = cr->i64;
291 for (int i = (repeat + 7) >> 3; i > 0; --i)
292 {
293 __m256i mr32 = _mm256_load_si256((__m256i*)rp);
294 __m256i mg32 = _mm256_load_si256((__m256i*)gp);
295 __m256i mb32 = _mm256_load_si256((__m256i*)bp);
296 __m256i mr, mg, mb, t;
297 mr = _mm256_cvtepi32_epi64(_mm256_extracti128_si256(mr32, 0));
298 mg = _mm256_cvtepi32_epi64(_mm256_extracti128_si256(mg32, 0));
299 mb = _mm256_cvtepi32_epi64(_mm256_extracti128_si256(mb32, 0));
300
301 t = _mm256_add_epi64(mr, mb);
302 t = _mm256_add_epi64(t, _mm256_slli_epi64(mg, 1));
303 _mm256_store_si256((__m256i*)yp, avx2_mm256_srai_epi64(t, 2, v2));
304 t = _mm256_sub_epi64(mb, mg);
305 _mm256_store_si256((__m256i*)cbp, t);
306 t = _mm256_sub_epi64(mr, mg);
307 _mm256_store_si256((__m256i*)crp, t);
308
309 yp += 4; cbp += 4; crp += 4;
310
311 mr = _mm256_cvtepi32_epi64(_mm256_extracti128_si256(mr32, 1));
312 mg = _mm256_cvtepi32_epi64(_mm256_extracti128_si256(mg32, 1));
313 mb = _mm256_cvtepi32_epi64(_mm256_extracti128_si256(mb32, 1));
314
315 t = _mm256_add_epi64(mr, mb);
316 t = _mm256_add_epi64(t, _mm256_slli_epi64(mg, 1));
317 _mm256_store_si256((__m256i*)yp, avx2_mm256_srai_epi64(t, 2, v2));
318 t = _mm256_sub_epi64(mb, mg);
319 _mm256_store_si256((__m256i*)cbp, t);
320 t = _mm256_sub_epi64(mr, mg);
321 _mm256_store_si256((__m256i*)crp, t);
322
323 rp += 8; gp += 8; bp += 8;
324 yp += 4; cbp += 4; crp += 4;
325 }
326 }
327 }
328
331 const line_buf *cb,
332 const line_buf *cr,
333 line_buf *r, line_buf *g, line_buf *b,
334 ui32 repeat)
335 {
336 assert((y->flags & line_buf::LFT_REVERSIBLE) &&
342
343 if (y->flags & line_buf::LFT_32BIT)
344 {
345 assert((y->flags & line_buf::LFT_32BIT) &&
346 (cb->flags & line_buf::LFT_32BIT) &&
347 (cr->flags & line_buf::LFT_32BIT) &&
348 (r->flags & line_buf::LFT_32BIT) &&
349 (g->flags & line_buf::LFT_32BIT) &&
351 const si32 *yp = y->i32, *cbp = cb->i32, *crp = cr->i32;
352 si32 *rp = r->i32, *gp = g->i32, *bp = b->i32;
353 for (int i = (repeat + 7) >> 3; i > 0; --i)
354 {
355 __m256i my = _mm256_load_si256((__m256i*)yp);
356 __m256i mcb = _mm256_load_si256((__m256i*)cbp);
357 __m256i mcr = _mm256_load_si256((__m256i*)crp);
358
359 __m256i t = _mm256_add_epi32(mcb, mcr);
360 t = _mm256_sub_epi32(my, _mm256_srai_epi32(t, 2));
361 _mm256_store_si256((__m256i*)gp, t);
362 __m256i u = _mm256_add_epi32(mcb, t);
363 _mm256_store_si256((__m256i*)bp, u);
364 u = _mm256_add_epi32(mcr, t);
365 _mm256_store_si256((__m256i*)rp, u);
366
367 yp += 8; cbp += 8; crp += 8;
368 rp += 8; gp += 8; bp += 8;
369 }
370 }
371 else
372 {
373 assert((y->flags & line_buf::LFT_64BIT) &&
374 (cb->flags & line_buf::LFT_64BIT) &&
375 (cr->flags & line_buf::LFT_64BIT) &&
376 (r->flags & line_buf::LFT_32BIT) &&
377 (g->flags & line_buf::LFT_32BIT) &&
379 __m256i v2 = _mm256_set1_epi64x(1ULL << (63 - 2));
380 __m256i low_bits = _mm256_set_epi64x(0, (si64)ULLONG_MAX,
381 0, (si64)ULLONG_MAX);
382 const si64 *yp = y->i64, *cbp = cb->i64, *crp = cr->i64;
383 si32 *rp = r->i32, *gp = g->i32, *bp = b->i32;
384 for (int i = (repeat + 7) >> 3; i > 0; --i)
385 {
386 __m256i my, mcb, mcr, tr, tg, tb;
387 my = _mm256_load_si256((__m256i*)yp);
388 mcb = _mm256_load_si256((__m256i*)cbp);
389 mcr = _mm256_load_si256((__m256i*)crp);
390
391 tg = _mm256_add_epi64(mcb, mcr);
392 tg = _mm256_sub_epi64(my, avx2_mm256_srai_epi64(tg, 2, v2));
393 tb = _mm256_add_epi64(mcb, tg);
394 tr = _mm256_add_epi64(mcr, tg);
395
396 __m256i mr, mg, mb;
397 mr = _mm256_shuffle_epi32(tr, _MM_SHUFFLE(0, 0, 2, 0));
398 mr = _mm256_and_si256(low_bits, mr);
399 mg = _mm256_shuffle_epi32(tg, _MM_SHUFFLE(0, 0, 2, 0));
400 mg = _mm256_and_si256(low_bits, mg);
401 mb = _mm256_shuffle_epi32(tb, _MM_SHUFFLE(0, 0, 2, 0));
402 mb = _mm256_and_si256(low_bits, mb);
403
404 yp += 4; cbp += 4; crp += 4;
405
406 my = _mm256_load_si256((__m256i*)yp);
407 mcb = _mm256_load_si256((__m256i*)cbp);
408 mcr = _mm256_load_si256((__m256i*)crp);
409
410 tg = _mm256_add_epi64(mcb, mcr);
411 tg = _mm256_sub_epi64(my, avx2_mm256_srai_epi64(tg, 2, v2));
412 tb = _mm256_add_epi64(mcb, tg);
413 tr = _mm256_add_epi64(mcr, tg);
414
415 tr = _mm256_shuffle_epi32(tr, _MM_SHUFFLE(2, 0, 0, 0));
416 tr = _mm256_andnot_si256(low_bits, tr);
417 mr = _mm256_or_si256(mr, tr);
418 mr = _mm256_permute4x64_epi64(mr, _MM_SHUFFLE(3, 1, 2, 0));
419
420 tg = _mm256_shuffle_epi32(tg, _MM_SHUFFLE(2, 0, 0, 0));
421 tg = _mm256_andnot_si256(low_bits, tg);
422 mg = _mm256_or_si256(mg, tg);
423 mg = _mm256_permute4x64_epi64(mg, _MM_SHUFFLE(3, 1, 2, 0));
424
425 tb = _mm256_shuffle_epi32(tb, _MM_SHUFFLE(2, 0, 0, 0));
426 tb = _mm256_andnot_si256(low_bits, tb);
427 mb = _mm256_or_si256(mb, tb);
428 mb = _mm256_permute4x64_epi64(mb, _MM_SHUFFLE(3, 1, 2, 0));
429
430 _mm256_store_si256((__m256i*)rp, mr);
431 _mm256_store_si256((__m256i*)gp, mg);
432 _mm256_store_si256((__m256i*)bp, mb);
433
434 yp += 4; cbp += 4; crp += 4;
435 rp += 8; gp += 8; bp += 8;
436 }
437 }
438 }
439
440 }
441}
si64 * i64
Definition: ojph_mem.h:173
si32 * i32
Definition: ojph_mem.h:172
void avx2_rct_forward(const line_buf *r, const line_buf *g, const line_buf *b, line_buf *y, line_buf *cb, line_buf *cr, ui32 repeat)
void avx2_rct_backward(const line_buf *y, const line_buf *cb, const line_buf *cr, line_buf *r, line_buf *g, line_buf *b, ui32 repeat)
void avx2_rev_convert(const line_buf *src_line, const ui32 src_line_offset, line_buf *dst_line, const ui32 dst_line_offset, si64 shift, ui32 width)
void avx2_rev_convert_nlt_type3(const line_buf *src_line, const ui32 src_line_offset, line_buf *dst_line, const ui32 dst_line_offset, si64 shift, ui32 width)
static __m256i avx2_mm256_srai_epi64(__m256i a, int amt, __m256i m)
int64_t si64
Definition: ojph_defs.h:57
int32_t si32
Definition: ojph_defs.h:55
uint32_t ui32
Definition: ojph_defs.h:54