OpenJPH
Open-source implementation of JPEG2000 Part-15
Loading...
Searching...
No Matches
ojph_colour_sse2.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_sse2.cpp
34// Author: Aous Naman
35// Date: 11 October 2019
36//***************************************************************************/
37
38#include "ojph_arch.h"
39#if defined(OJPH_ARCH_I386) || defined(OJPH_ARCH_X86_64)
40
41#include <climits>
42#include <cmath>
43
44#include "ojph_defs.h"
45#include "ojph_mem.h"
46#include "ojph_colour.h"
47
48#include "ojph_params.h"
49#include "../codestream/ojph_params_local.h"
50
51#include <emmintrin.h>
52
53namespace ojph {
54 namespace local {
55
57 void sse2_cnvrt_float_to_si32_shftd(const float *sp, si32 *dp, float mul,
58 ui32 width)
59 {
60 uint32_t rounding_mode = _MM_GET_ROUNDING_MODE();
61 _MM_SET_ROUNDING_MODE(_MM_ROUND_NEAREST);
62 __m128 shift = _mm_set1_ps(0.5f);
63 __m128 m = _mm_set1_ps(mul);
64 for (int i = (width + 3) >> 2; i > 0; --i, sp+=4, dp+=4)
65 {
66 __m128 t = _mm_loadu_ps(sp);
67 __m128 s = _mm_add_ps(t, shift);
68 s = _mm_mul_ps(s, m);
69 _mm_storeu_si128((__m128i*)dp, _mm_cvtps_epi32(s));
70 }
71 _MM_SET_ROUNDING_MODE(rounding_mode);
72 }
73
75 void sse2_cnvrt_float_to_si32(const float *sp, si32 *dp, float mul,
76 ui32 width)
77 {
78 uint32_t rounding_mode = _MM_GET_ROUNDING_MODE();
79 _MM_SET_ROUNDING_MODE(_MM_ROUND_NEAREST);
80 __m128 m = _mm_set1_ps(mul);
81 for (int i = (width + 3) >> 2; i > 0; --i, sp+=4, dp+=4)
82 {
83 __m128 t = _mm_loadu_ps(sp);
84 __m128 s = _mm_mul_ps(t, m);
85 _mm_storeu_si128((__m128i*)dp, _mm_cvtps_epi32(s));
86 }
87 _MM_SET_ROUNDING_MODE(rounding_mode);
88 }
89
91 static inline
92 __m128i ojph_mm_max_ge_epi32(__m128i a, __m128i b, __m128 x, __m128 y)
93 {
94 __m128 ct = _mm_cmpge_ps(x, y); // 0xFFFFFFFF for x >= y
95 __m128i c = _mm_castps_si128(ct); // does not generate any code
96 __m128i d = _mm_and_si128(c, a); // keep only a, where x >= y
97 __m128i e = _mm_andnot_si128(c, b); // keep only b, where x < y
98 return _mm_or_si128(d, e); // combine
99 }
100
102 static inline
103 __m128i ojph_mm_min_lt_epi32(__m128i a, __m128i b, __m128 x, __m128 y)
104 {
105 __m128 ct = _mm_cmplt_ps(x, y); // 0xFFFFFFFF for x < y
106 __m128i c = _mm_castps_si128(ct); // does not generate any code
107 __m128i d = _mm_and_si128(c, a); // keep only a, where x < y
108 __m128i e = _mm_andnot_si128(c, b); // keep only b, where x >= y
109 return _mm_or_si128(d, e); // combine
110 }
111
113 template <bool NLT_TYPE3>
114 static inline
115 void local_sse2_irv_convert_to_integer(const line_buf *src_line,
116 line_buf *dst_line, ui32 dst_line_offset,
117 ui32 bit_depth, bool is_signed, ui32 width)
118 {
119 assert((src_line->flags & line_buf::LFT_32BIT) &&
120 (src_line->flags & line_buf::LFT_INTEGER) == 0 &&
121 (dst_line->flags & line_buf::LFT_32BIT) &&
122 (dst_line->flags & line_buf::LFT_INTEGER));
123
124 assert(bit_depth <= 32);
125 uint32_t rounding_mode = _MM_GET_ROUNDING_MODE();
126 _MM_SET_ROUNDING_MODE(_MM_ROUND_NEAREST);
127
128 const float* sp = src_line->f32;
129 si32* dp = dst_line->i32 + dst_line_offset;
130 // There is the possibility that converting to integer will
131 // exceed the dynamic range of 32bit integer; therefore, care must be
132 // exercised.
133 // We look if the floating point number is outside the half-closed
134 // interval [-0.5f, 0.5f). If so, we limit the resulting integer
135 // to the maximum/minimum that number supports.
136 si32 neg_limit = (si32)INT_MIN >> (32 - bit_depth);
137 __m128 mul = _mm_set1_ps((float)(1ull << bit_depth));
138 __m128 fl_up_lim = _mm_set1_ps(-(float)neg_limit); // val < upper
139 __m128 fl_low_lim = _mm_set1_ps((float)neg_limit); // val >= lower
140 __m128i s32_up_lim = _mm_set1_epi32(INT_MAX >> (32 - bit_depth));
141 __m128i s32_low_lim = _mm_set1_epi32(INT_MIN >> (32 - bit_depth));
142
143 if (is_signed)
144 {
145 __m128i zero = _mm_setzero_si128();
146 __m128i bias = _mm_set1_epi32(-(si32)((1ULL << (bit_depth - 1)) + 1));
147 for (int i = (int)width; i > 0; i -= 4, sp += 4, dp += 4) {
148 __m128 t = _mm_loadu_ps(sp);
149 t = _mm_mul_ps(t, mul);
150 __m128i u = _mm_cvtps_epi32(t);
151 u = ojph_mm_max_ge_epi32(u, s32_low_lim, t, fl_low_lim);
152 u = ojph_mm_min_lt_epi32(u, s32_up_lim, t, fl_up_lim);
153 if (NLT_TYPE3)
154 {
155 __m128i c = _mm_cmpgt_epi32(zero, u); //0xFFFFFFFF for -ve value
156 __m128i neg = _mm_sub_epi32(bias, u); //-bias -value
157 neg = _mm_and_si128(c, neg); //keep only - bias - value
158 u = _mm_andnot_si128(c, u); //keep only +ve or 0
159 u = _mm_or_si128(neg, u); //combine
160 }
161 _mm_storeu_si128((__m128i*)dp, u);
162 }
163 }
164 else
165 {
166 __m128i half = _mm_set1_epi32((si32)(1ULL << (bit_depth - 1)));
167 for (int i = (int)width; i > 0; i -= 4, sp += 4, dp += 4) {
168 __m128 t = _mm_loadu_ps(sp);
169 t = _mm_mul_ps(t, mul);
170 __m128i u = _mm_cvtps_epi32(t);
171 u = ojph_mm_max_ge_epi32(u, s32_low_lim, t, fl_low_lim);
172 u = ojph_mm_min_lt_epi32(u, s32_up_lim, t, fl_up_lim);
173 u = _mm_add_epi32(u, half);
174 _mm_storeu_si128((__m128i*)dp, u);
175 }
176 }
177
178 _MM_SET_ROUNDING_MODE(rounding_mode);
179 }
180
182 void sse2_irv_convert_to_integer(const line_buf *src_line,
183 line_buf *dst_line, ui32 dst_line_offset,
184 ui32 bit_depth, bool is_signed, ui32 width)
185 {
186 local_sse2_irv_convert_to_integer<false>(src_line, dst_line,
187 dst_line_offset, bit_depth, is_signed, width);
188 }
189
191 void sse2_irv_convert_to_integer_nlt_type3(const line_buf *src_line,
192 line_buf *dst_line, ui32 dst_line_offset,
193 ui32 bit_depth, bool is_signed, ui32 width)
194 {
195 local_sse2_irv_convert_to_integer<true>(src_line, dst_line,
196 dst_line_offset, bit_depth, is_signed, width);
197 }
198
200 template<int NLT_TYPE>
201 static inline
202 void local_sse2_irv_convert_to_integer_nlt2or4(const line_buf *src_line,
203 line_buf *dst_line, ui32 dst_line_offset,
204 ui32 bit_depth, bool is_signed, ui32 width, const nlt_rec* rec)
205 {
206 ojph_unused(bit_depth);
207 ojph_unused(is_signed);
208
209 assert((src_line->flags & line_buf::LFT_32BIT) &&
210 (src_line->flags & line_buf::LFT_INTEGER) == 0 &&
211 (dst_line->flags & line_buf::LFT_32BIT) &&
212 (dst_line->flags & line_buf::LFT_INTEGER));
213
214 assert(bit_depth <= 32);
215 uint32_t rounding_mode = _MM_GET_ROUNDING_MODE();
216 _MM_SET_ROUNDING_MODE(_MM_ROUND_NEAREST);
217
218 const float* sp = src_line->f32;
219 si32* dp = dst_line->i32 + dst_line_offset;
220
221 __m128 mul = _mm_set1_ps(rec->multiplier);
222 __m128 d_min = _mm_set1_ps(rec->fd_min);
223 __m128 d_max = _mm_set1_ps(rec->fd_max);
224 __m128 delta = _mm_set1_ps(rec->delta);
225 __m128 inv_delta = _mm_set1_ps(rec->inv_delta);
226 const float* lut = rec->dec_points;
227
228 __m128 half_ps = _mm_set1_ps(0.5f);
229
230 if (rec->is_signed())
231 {
232 __m128 half = _mm_set1_ps((float)(1ULL << (rec->get_bit_depth() - 1)));
233 __m128i bias =
234 _mm_set1_epi32(-(si32)((1ULL << (rec->get_bit_depth() - 1)) + 1));
235 __m128i zero = _mm_setzero_si128();
236 for (int i = (int)width; i > 0; i -= 4, sp += 4, dp += 4) {
237 __m128 t = _mm_loadu_ps(sp);
238 t = _mm_add_ps(t, half_ps); // convert to [0, 1]
239 t = _mm_max_ps(t, d_min);
240 t = _mm_min_ps(t, d_max);
241 __m128i k = _mm_cvttps_epi32(
242 _mm_mul_ps(_mm_sub_ps(t, d_min), inv_delta));
243 __m128 d_k = _mm_add_ps(d_min,
244 _mm_mul_ps(_mm_cvtepi32_ps(k), delta));
245 // SSE2 has no gather; perform the LUT lookup 4 times and build
246 // the vector from the individual results
247 si32 kk[4];
248 _mm_storeu_si128((__m128i*)kk, k);
249 __m128 t_k = _mm_set_ps(lut[kk[3]], lut[kk[2]],
250 lut[kk[1]], lut[kk[0]]);
251 __m128 t_kp1 = _mm_set_ps(lut[kk[3] + 1], lut[kk[2] + 1],
252 lut[kk[1] + 1], lut[kk[0] + 1]);
253 __m128 z = _mm_add_ps(t_k,
254 _mm_mul_ps(_mm_mul_ps(_mm_sub_ps(t, d_k), inv_delta),
255 _mm_sub_ps(t_kp1, t_k)));
256 __m128i v =
257 _mm_cvtps_epi32(_mm_sub_ps(_mm_mul_ps(z, mul), half));
258 if (NLT_TYPE == 4)
259 {
260 __m128i c = _mm_cmpgt_epi32(zero, v); // 0xFFFFFFFF for -ve val
261 __m128i neg = _mm_sub_epi32(bias, v); // - bias - value
262 neg = _mm_and_si128(c, neg); // keep only - bias - val
263 v = _mm_andnot_si128(c, v); // keep only +ve or 0
264 v = _mm_or_si128(neg, v); // combine
265 }
266 _mm_storeu_si128((__m128i*)dp, v);
267 }
268 }
269 else
270 {
271 for (int i = (int)width; i > 0; i -= 4, sp += 4, dp += 4) {
272 __m128 t = _mm_loadu_ps(sp);
273 t = _mm_add_ps(t, half_ps); // convert to [0, 1]
274 t = _mm_max_ps(t, d_min);
275 t = _mm_min_ps(t, d_max);
276 __m128i k = _mm_cvttps_epi32(
277 _mm_mul_ps(_mm_sub_ps(t, d_min), inv_delta));
278 __m128 d_k = _mm_add_ps(d_min,
279 _mm_mul_ps(_mm_cvtepi32_ps(k), delta));
280 // SSE2 has no gather; perform the LUT lookup 4 times and build
281 // the vector from the individual results
282 si32 kk[4];
283 _mm_storeu_si128((__m128i*)kk, k);
284 __m128 t_k = _mm_set_ps(lut[kk[3]], lut[kk[2]],
285 lut[kk[1]], lut[kk[0]]);
286 __m128 t_kp1 = _mm_set_ps(lut[kk[3] + 1], lut[kk[2] + 1],
287 lut[kk[1] + 1], lut[kk[0] + 1]);
288 __m128 z = _mm_add_ps(t_k,
289 _mm_mul_ps(_mm_mul_ps(_mm_sub_ps(t, d_k), inv_delta),
290 _mm_sub_ps(t_kp1, t_k)));
291 __m128i v = _mm_cvtps_epi32(_mm_mul_ps(z, mul));
292 _mm_storeu_si128((__m128i*)dp, v);
293 }
294 }
295
296 _MM_SET_ROUNDING_MODE(rounding_mode);
297 }
298
300 void sse2_irv_convert_to_integer_nlt(const line_buf *src_line,
301 line_buf *dst_line, ui32 dst_line_offset,
302 ui32 bit_depth, bool is_signed, ui32 width, const nlt_rec* rec)
303 {
304 using nl = nlt_rec::nonlinearity;
305 if (rec->get_type() == nl::OJPH_NLT_LUT_STYLE_NLT)
306 local_sse2_irv_convert_to_integer_nlt2or4<2>(src_line, dst_line,
307 dst_line_offset, bit_depth, is_signed, width, rec);
308 else if (rec->get_type() == nl::OJPH_NLT_BINARY_COMPLEMENT_PLUS_LUT)
309 local_sse2_irv_convert_to_integer_nlt2or4<4>(src_line, dst_line,
310 dst_line_offset, bit_depth, is_signed, width, rec);
311 else
312 assert(0);
313 }
314
316 // https://github.com/seung-lab/dijkstra3d/blob/master/libdivide.h
317 static inline __m128i sse2_mm_srai_epi64(__m128i a, int amt, __m128i m)
318 {
319 // note than m must be obtained using
320 // __m128i m = _mm_set1_epi64x(1ULL << (63 - amt));
321 __m128i x = _mm_srli_epi64(a, amt);
322 x = _mm_xor_si128(x, m);
323 __m128i result = _mm_sub_epi64(x, m);
324 return result;
325 }
326
328 static inline __m128i sse2_cvtlo_epi32_epi64(__m128i a, __m128i zero)
329 {
330 __m128i t;
331 t = _mm_cmplt_epi32(a, zero); // get -ve
332 t = _mm_unpacklo_epi32(a, t);
333 return t;
334 }
335
337 static inline __m128i sse2_cvthi_epi32_epi64(__m128i a, __m128i zero)
338 {
339 __m128i t;
340 t = _mm_cmplt_epi32(a, zero); // get -ve
341 t = _mm_unpackhi_epi32(a, t);
342 return t;
343 }
344
346 void sse2_rev_convert(const line_buf *src_line,
347 const ui32 src_line_offset,
348 line_buf *dst_line,
349 const ui32 dst_line_offset,
350 si64 shift, ui32 width)
351 {
352 if (src_line->flags & line_buf::LFT_32BIT)
353 {
354 if (dst_line->flags & line_buf::LFT_32BIT)
355 {
356 const si32 *sp = src_line->i32 + src_line_offset;
357 si32 *dp = dst_line->i32 + dst_line_offset;
358 __m128i sh = _mm_set1_epi32((si32)shift);
359 for (int i = (width + 3) >> 2; i > 0; --i, sp+=4, dp+=4)
360 {
361 __m128i s = _mm_loadu_si128((__m128i*)sp);
362 s = _mm_add_epi32(s, sh);
363 _mm_storeu_si128((__m128i*)dp, s);
364 }
365 }
366 else
367 {
368 const si32 *sp = src_line->i32 + src_line_offset;
369 si64 *dp = dst_line->i64 + dst_line_offset;
370 __m128i zero = _mm_setzero_si128();
371 __m128i sh = _mm_set1_epi64x(shift);
372 for (int i = (width + 3) >> 2; i > 0; --i, sp+=4, dp+=4)
373 {
374 __m128i s, t;
375 s = _mm_loadu_si128((__m128i*)sp);
376
377 t = sse2_cvtlo_epi32_epi64(s, zero);
378 t = _mm_add_epi64(t, sh);
379 _mm_storeu_si128((__m128i*)dp, t);
380
381 t = sse2_cvthi_epi32_epi64(s, zero);
382 t = _mm_add_epi64(t, sh);
383 _mm_storeu_si128((__m128i*)dp + 1, t);
384 }
385 }
386 }
387 else
388 {
389 assert(src_line->flags | line_buf::LFT_64BIT);
390 assert(dst_line->flags | line_buf::LFT_32BIT);
391 const si64 *sp = src_line->i64 + src_line_offset;
392 si32 *dp = dst_line->i32 + dst_line_offset;
393 __m128i low_bits = _mm_set_epi64x(0, (si64)ULLONG_MAX);
394 __m128i sh = _mm_set1_epi64x(shift);
395 for (int i = (width + 3) >> 2; i > 0; --i, sp+=4, dp+=4)
396 {
397 __m128i s, t;
398 s = _mm_loadu_si128((__m128i*)sp);
399 s = _mm_add_epi64(s, sh);
400
401 t = _mm_shuffle_epi32(s, _MM_SHUFFLE(0, 0, 2, 0));
402 t = _mm_and_si128(low_bits, t);
403
404 s = _mm_loadu_si128((__m128i*)sp + 1);
405 s = _mm_add_epi64(s, sh);
406
407 s = _mm_shuffle_epi32(s, _MM_SHUFFLE(2, 0, 0, 0));
408 s = _mm_andnot_si128(low_bits, s);
409
410 t = _mm_or_si128(s, t);
411 _mm_storeu_si128((__m128i*)dp, t);
412 }
413 }
414 }
415
417 void sse2_rev_convert_nlt_type3(const line_buf *src_line,
418 const ui32 src_line_offset,
419 line_buf *dst_line,
420 const ui32 dst_line_offset,
421 si64 shift, ui32 width)
422 {
423 if (src_line->flags & line_buf::LFT_32BIT)
424 {
425 if (dst_line->flags & line_buf::LFT_32BIT)
426 {
427 const si32 *sp = src_line->i32 + src_line_offset;
428 si32 *dp = dst_line->i32 + dst_line_offset;
429 __m128i sh = _mm_set1_epi32((si32)(-shift));
430 __m128i zero = _mm_setzero_si128();
431 for (int i = (width + 3) >> 2; i > 0; --i, sp += 4, dp += 4)
432 {
433 __m128i s = _mm_loadu_si128((__m128i*)sp);
434 __m128i c = _mm_cmplt_epi32(s, zero); // 0xFFFFFFFF for -ve value
435 __m128i v_m_sh = _mm_sub_epi32(sh, s); // - shift - value
436 v_m_sh = _mm_and_si128(c, v_m_sh); // keep only - shift - value
437 s = _mm_andnot_si128(c, s); // keep only +ve or 0
438 s = _mm_or_si128(s, v_m_sh); // combine
439 _mm_storeu_si128((__m128i*)dp, s);
440 }
441 }
442 else
443 {
444 const si32 *sp = src_line->i32 + src_line_offset;
445 si64 *dp = dst_line->i64 + dst_line_offset;
446 __m128i sh = _mm_set1_epi64x(-shift);
447 __m128i zero = _mm_setzero_si128();
448 for (int i = (width + 3) >> 2; i > 0; --i, sp += 4, dp += 4)
449 {
450 __m128i s, t, u, c, v_m_sh;
451 s = _mm_loadu_si128((__m128i*)sp);
452
453 t = _mm_cmplt_epi32(s, zero); // find -ve 32bit -1
454 u = _mm_unpacklo_epi32(s, t); // correct 64bit data
455 c = _mm_unpacklo_epi32(t, t); // 64bit -1 for -ve value
456
457 v_m_sh = _mm_sub_epi64(sh, u); // - shift - value
458 v_m_sh = _mm_and_si128(c, v_m_sh); // keep only - shift - value
459 u = _mm_andnot_si128(c, u); // keep only +ve or 0
460 u = _mm_or_si128(u, v_m_sh); // combine
461
462 _mm_storeu_si128((__m128i*)dp, u);
463 u = _mm_unpackhi_epi32(s, t); // correct 64bit data
464 c = _mm_unpackhi_epi32(t, t); // 64bit -1 for -ve value
465
466 v_m_sh = _mm_sub_epi64(sh, u); // - shift - value
467 v_m_sh = _mm_and_si128(c, v_m_sh); // keep only - shift - value
468 u = _mm_andnot_si128(c, u); // keep only +ve or 0
469 u = _mm_or_si128(u, v_m_sh); // combine
470
471 _mm_storeu_si128((__m128i*)dp + 1, u);
472 }
473 }
474 }
475 else
476 {
477 assert(src_line->flags | line_buf::LFT_64BIT);
478 assert(dst_line->flags | line_buf::LFT_32BIT);
479 const si64 *sp = src_line->i64 + src_line_offset;
480 si32 *dp = dst_line->i32 + dst_line_offset;
481 __m128i sh = _mm_set1_epi64x(-shift);
482 __m128i zero = _mm_setzero_si128();
483 __m128i half_mask = _mm_set_epi64x(0, (si64)ULLONG_MAX);
484 for (int i = (width + 3) >> 2; i > 0; --i, sp += 4, dp += 4)
485 {
486 // s for source, t for target, p for positive, n for negative,
487 // m for mask, and tm for temp
488 __m128i s, t, p, n, m, tm;
489 s = _mm_loadu_si128((__m128i*)sp);
490
491 tm = _mm_cmplt_epi32(s, zero); // 32b -1 for -ve value
492 m = _mm_shuffle_epi32(tm, _MM_SHUFFLE(3, 3, 1, 1)); // expand to 64b
493 tm = _mm_sub_epi64(sh, s); // - shift - value
494 n = _mm_and_si128(m, tm); // -ve
495 p = _mm_andnot_si128(m, s); // +ve
496 tm = _mm_or_si128(n, p);
497 tm = _mm_shuffle_epi32(tm, _MM_SHUFFLE(0, 0, 2, 0));
498 t = _mm_and_si128(half_mask, tm);
499
500 s = _mm_loadu_si128((__m128i*)sp + 1);
501 tm = _mm_cmplt_epi32(s, zero); // 32b -1 for -ve value
502 m = _mm_shuffle_epi32(tm, _MM_SHUFFLE(3, 3, 1, 1)); // expand to 64b
503 tm = _mm_sub_epi64(sh, s); // - shift - value
504 n = _mm_and_si128(m, tm); // -ve
505 p = _mm_andnot_si128(m, s); // +ve
506 tm = _mm_or_si128(n, p);
507 tm = _mm_shuffle_epi32(tm, _MM_SHUFFLE(2, 0, 0, 0));
508 tm = _mm_andnot_si128(half_mask, tm);
509
510 t = _mm_or_si128(t, tm);
511 _mm_storeu_si128((__m128i*)dp, t);
512 }
513 }
514 }
515
517 template<bool NLT_TYPE3>
518 static inline
519 void local_sse2_irv_convert_to_float(const line_buf *src_line,
520 ui32 src_line_offset, line_buf *dst_line,
521 ui32 bit_depth, bool is_signed, ui32 width)
522 {
523 assert((src_line->flags & line_buf::LFT_32BIT) &&
524 (src_line->flags & line_buf::LFT_INTEGER) &&
525 (dst_line->flags & line_buf::LFT_32BIT) &&
526 (dst_line->flags & line_buf::LFT_INTEGER) == 0);
527
528 assert(bit_depth <= 32);
529 __m128 mul = _mm_set1_ps((float)(1.0 / (double)(1ULL << bit_depth)));
530
531 const si32* sp = src_line->i32 + src_line_offset;
532 float* dp = dst_line->f32;
533 if (is_signed)
534 {
535 __m128i zero = _mm_setzero_si128();
536 __m128i bias = _mm_set1_epi32(-(si32)((1ULL << (bit_depth - 1)) + 1));
537 for (int i = (int)width; i > 0; i -= 4, sp += 4, dp += 4) {
538 __m128i t = _mm_loadu_si128((__m128i*)sp);
539 if (NLT_TYPE3)
540 {
541 __m128i c = _mm_cmplt_epi32(t, zero); // 0xFFFFFFFF for -ve value
542 __m128i neg = _mm_sub_epi32(bias, t); // - bias - value
543 neg = _mm_and_si128(c, neg); // keep only - bias - value
544 c = _mm_andnot_si128(c, t); // keep only +ve or 0
545 t = _mm_or_si128(neg, c); // combine
546 }
547 __m128 v = _mm_cvtepi32_ps(t);
548 v = _mm_mul_ps(v, mul);
549 _mm_storeu_ps(dp, v);
550 }
551 }
552 else
553 {
554 __m128i half = _mm_set1_epi32((si32)(1ULL << (bit_depth - 1)));
555 for (int i = (int)width; i > 0; i -= 4, sp += 4, dp += 4) {
556 __m128i t = _mm_loadu_si128((__m128i*)sp);
557 t = _mm_sub_epi32(t, half);
558 __m128 v = _mm_cvtepi32_ps(t);
559 v = _mm_mul_ps(v, mul);
560 _mm_storeu_ps(dp, v);
561 }
562 }
563 }
564
566 void sse2_irv_convert_to_float(const line_buf *src_line,
567 ui32 src_line_offset, line_buf *dst_line,
568 ui32 bit_depth, bool is_signed, ui32 width)
569 {
570 local_sse2_irv_convert_to_float<false>(src_line, src_line_offset,
571 dst_line, bit_depth, is_signed, width);
572 }
573
575 void sse2_irv_convert_to_float_nlt_type3(const line_buf *src_line,
576 ui32 src_line_offset, line_buf *dst_line,
577 ui32 bit_depth, bool is_signed, ui32 width)
578 {
579 local_sse2_irv_convert_to_float<true>(src_line, src_line_offset,
580 dst_line, bit_depth, is_signed, width);
581 }
582
584 template<int NLT_TYPE>
585 static inline
586 void local_sse2_irv_convert_to_float_nlt2or4(const line_buf *src_line,
587 ui32 src_line_offset, line_buf *dst_line,
588 ui32 bit_depth, bool is_signed, ui32 width, const nlt_rec* rec)
589 {
590 ojph_unused(is_signed);
591
592 assert((src_line->flags & line_buf::LFT_32BIT) &&
593 (src_line->flags & line_buf::LFT_INTEGER) &&
594 (dst_line->flags & line_buf::LFT_32BIT) &&
595 (dst_line->flags & line_buf::LFT_INTEGER) == 0);
596
597 assert(bit_depth <= 32);
598 __m128 mul = _mm_set1_ps((float)(1.0 / (double)(1ULL << bit_depth)));
599 __m128 d_min = _mm_set1_ps(rec->ft_min);
600 __m128 d_max = _mm_set1_ps(rec->ft_max);
601 __m128 delta = _mm_set1_ps(rec->delta);
602 __m128 inv_delta = _mm_set1_ps(rec->inv_delta);
603 const float* lut = rec->enc_points;
604
605 __m128 half_ps = _mm_set1_ps(0.5f);
606
607 const si32* sp = src_line->i32 + src_line_offset;
608 float* dp = dst_line->f32;
609 if (rec->is_signed())
610 {
611 __m128i bias =
612 _mm_set1_epi32(-(si32)((1ULL << (rec->get_bit_depth() - 1)) + 1));
613 __m128i zero = _mm_setzero_si128();
614 for (int i = (int)width; i > 0; i -= 4, sp += 4, dp += 4) {
615 __m128i v = _mm_loadu_si128((__m128i*)sp);
616 if (NLT_TYPE == 4)
617 {
618 __m128i c = _mm_cmpgt_epi32(zero, v); // 0xFFFFFFFF for -ve val
619 __m128i neg = _mm_sub_epi32(bias, v); // - bias - value
620 neg = _mm_and_si128(c, neg); // keep only - bias - val
621 v = _mm_andnot_si128(c, v); // keep only +ve or 0
622 v = _mm_or_si128(neg, v); // combine
623 }
624 __m128 t = _mm_add_ps( // convert to [0, 1]
625 _mm_mul_ps(_mm_cvtepi32_ps(v), mul), half_ps);
626 t = _mm_max_ps(t, d_min);
627 t = _mm_min_ps(t, d_max);
628 __m128i k = _mm_cvttps_epi32(
629 _mm_mul_ps(_mm_sub_ps(t, d_min), inv_delta));
630 __m128 d_k = _mm_add_ps(d_min,
631 _mm_mul_ps(_mm_cvtepi32_ps(k), delta));
632 // SSE2 has no gather; perform the LUT lookup 4 times and build
633 // the vector from the individual results
634 si32 kk[4];
635 _mm_storeu_si128((__m128i*)kk, k);
636 __m128 t_k = _mm_set_ps(lut[kk[3]], lut[kk[2]],
637 lut[kk[1]], lut[kk[0]]);
638 __m128 t_kp1 = _mm_set_ps(lut[kk[3] + 1], lut[kk[2] + 1],
639 lut[kk[1] + 1], lut[kk[0] + 1]);
640 __m128 y = _mm_add_ps(t_k,
641 _mm_mul_ps(_mm_mul_ps(_mm_sub_ps(t, d_k), inv_delta),
642 _mm_sub_ps(t_kp1, t_k)));
643 _mm_storeu_ps(dp, _mm_sub_ps(y, half_ps));
644 }
645 }
646 else
647 {
648 for (int i = (int)width; i > 0; i -= 4, sp += 4, dp += 4) {
649 __m128i v = _mm_loadu_si128((__m128i*)sp);
650 __m128 t = _mm_mul_ps(_mm_cvtepi32_ps(v), mul); // in [0, 1]
651 t = _mm_max_ps(t, d_min);
652 t = _mm_min_ps(t, d_max);
653 __m128i k = _mm_cvttps_epi32(
654 _mm_mul_ps(_mm_sub_ps(t, d_min), inv_delta));
655 __m128 d_k = _mm_add_ps(d_min,
656 _mm_mul_ps(_mm_cvtepi32_ps(k), delta));
657 // SSE2 has no gather; perform the LUT lookup 4 times and build
658 // the vector from the individual results
659 si32 kk[4];
660 _mm_storeu_si128((__m128i*)kk, k);
661 __m128 t_k = _mm_set_ps(lut[kk[3]], lut[kk[2]],
662 lut[kk[1]], lut[kk[0]]);
663 __m128 t_kp1 = _mm_set_ps(lut[kk[3] + 1], lut[kk[2] + 1],
664 lut[kk[1] + 1], lut[kk[0] + 1]);
665 __m128 y = _mm_add_ps(t_k,
666 _mm_mul_ps(_mm_mul_ps(_mm_sub_ps(t, d_k), inv_delta),
667 _mm_sub_ps(t_kp1, t_k)));
668 _mm_storeu_ps(dp, _mm_sub_ps(y, half_ps));
669 }
670 }
671 }
672
674 void sse2_irv_convert_to_float_nlt(const line_buf *src_line,
675 ui32 src_line_offset, line_buf *dst_line,
676 ui32 bit_depth, bool is_signed, ui32 width, const nlt_rec* rec)
677 {
678 using nl = nlt_rec::nonlinearity;
679 if (rec->get_type() == nl::OJPH_NLT_LUT_STYLE_NLT)
680 local_sse2_irv_convert_to_float_nlt2or4<2>(src_line,
681 src_line_offset, dst_line, bit_depth, is_signed, width, rec);
682 else if (rec->get_type() == nl::OJPH_NLT_BINARY_COMPLEMENT_PLUS_LUT)
683 local_sse2_irv_convert_to_float_nlt2or4<4>(src_line,
684 src_line_offset, dst_line, bit_depth, is_signed, width, rec);
685 else
686 assert(0);
687 }
688
690 void sse2_rct_forward(const line_buf *r,
691 const line_buf *g,
692 const line_buf *b,
693 line_buf *y, line_buf *cb, line_buf *cr,
694 ui32 repeat)
695 {
696 assert((y->flags & line_buf::LFT_INTEGER) &&
697 (cb->flags & line_buf::LFT_INTEGER) &&
698 (cr->flags & line_buf::LFT_INTEGER) &&
699 (r->flags & line_buf::LFT_INTEGER) &&
700 (g->flags & line_buf::LFT_INTEGER) &&
701 (b->flags & line_buf::LFT_INTEGER));
702
703 if (y->flags & line_buf::LFT_32BIT)
704 {
705 assert((y->flags & line_buf::LFT_32BIT) &&
706 (cb->flags & line_buf::LFT_32BIT) &&
707 (cr->flags & line_buf::LFT_32BIT) &&
708 (r->flags & line_buf::LFT_32BIT) &&
709 (g->flags & line_buf::LFT_32BIT) &&
710 (b->flags & line_buf::LFT_32BIT));
711 const si32 *rp = r->i32, * gp = g->i32, * bp = b->i32;
712 si32 *yp = y->i32, * cbp = cb->i32, * crp = cr->i32;
713 for (int i = (repeat + 3) >> 2; i > 0; --i)
714 {
715 __m128i mr = _mm_load_si128((__m128i*)rp);
716 __m128i mg = _mm_load_si128((__m128i*)gp);
717 __m128i mb = _mm_load_si128((__m128i*)bp);
718 __m128i t = _mm_add_epi32(mr, mb);
719 t = _mm_add_epi32(t, _mm_slli_epi32(mg, 1));
720 _mm_store_si128((__m128i*)yp, _mm_srai_epi32(t, 2));
721 t = _mm_sub_epi32(mb, mg);
722 _mm_store_si128((__m128i*)cbp, t);
723 t = _mm_sub_epi32(mr, mg);
724 _mm_store_si128((__m128i*)crp, t);
725
726 rp += 4; gp += 4; bp += 4;
727 yp += 4; cbp += 4; crp += 4;
728 }
729 }
730 else
731 {
732 assert((y->flags & line_buf::LFT_64BIT) &&
733 (cb->flags & line_buf::LFT_64BIT) &&
734 (cr->flags & line_buf::LFT_64BIT) &&
735 (r->flags & line_buf::LFT_32BIT) &&
736 (g->flags & line_buf::LFT_32BIT) &&
737 (b->flags & line_buf::LFT_32BIT));
738 __m128i zero = _mm_setzero_si128();
739 __m128i v2 = _mm_set1_epi64x(1ULL << (63 - 2));
740 const si32 *rp = r->i32, *gp = g->i32, *bp = b->i32;
741 si64 *yp = y->i64, *cbp = cb->i64, *crp = cr->i64;
742 for (int i = (repeat + 3) >> 2; i > 0; --i)
743 {
744 __m128i mr32 = _mm_load_si128((__m128i*)rp);
745 __m128i mg32 = _mm_load_si128((__m128i*)gp);
746 __m128i mb32 = _mm_load_si128((__m128i*)bp);
747 __m128i mr, mg, mb, t;
748 mr = sse2_cvtlo_epi32_epi64(mr32, zero);
749 mg = sse2_cvtlo_epi32_epi64(mg32, zero);
750 mb = sse2_cvtlo_epi32_epi64(mb32, zero);
751
752 t = _mm_add_epi64(mr, mb);
753 t = _mm_add_epi64(t, _mm_slli_epi64(mg, 1));
754 _mm_store_si128((__m128i*)yp, sse2_mm_srai_epi64(t, 2, v2));
755 t = _mm_sub_epi64(mb, mg);
756 _mm_store_si128((__m128i*)cbp, t);
757 t = _mm_sub_epi64(mr, mg);
758 _mm_store_si128((__m128i*)crp, t);
759
760 yp += 2; cbp += 2; crp += 2;
761
762 mr = sse2_cvthi_epi32_epi64(mr32, zero);
763 mg = sse2_cvthi_epi32_epi64(mg32, zero);
764 mb = sse2_cvthi_epi32_epi64(mb32, zero);
765
766 t = _mm_add_epi64(mr, mb);
767 t = _mm_add_epi64(t, _mm_slli_epi64(mg, 1));
768 _mm_store_si128((__m128i*)yp, sse2_mm_srai_epi64(t, 2, v2));
769 t = _mm_sub_epi64(mb, mg);
770 _mm_store_si128((__m128i*)cbp, t);
771 t = _mm_sub_epi64(mr, mg);
772 _mm_store_si128((__m128i*)crp, t);
773
774 rp += 4; gp += 4; bp += 4;
775 yp += 2; cbp += 2; crp += 2;
776 }
777 }
778 }
779
781 void sse2_rct_backward(const line_buf *y,
782 const line_buf *cb,
783 const line_buf *cr,
784 line_buf *r, line_buf *g, line_buf *b,
785 ui32 repeat)
786 {
787 assert((y->flags & line_buf::LFT_INTEGER) &&
788 (cb->flags & line_buf::LFT_INTEGER) &&
789 (cr->flags & line_buf::LFT_INTEGER) &&
790 (r->flags & line_buf::LFT_INTEGER) &&
791 (g->flags & line_buf::LFT_INTEGER) &&
792 (b->flags & line_buf::LFT_INTEGER));
793
794 if (y->flags & line_buf::LFT_32BIT)
795 {
796 assert((y->flags & line_buf::LFT_32BIT) &&
797 (cb->flags & line_buf::LFT_32BIT) &&
798 (cr->flags & line_buf::LFT_32BIT) &&
799 (r->flags & line_buf::LFT_32BIT) &&
800 (g->flags & line_buf::LFT_32BIT) &&
801 (b->flags & line_buf::LFT_32BIT));
802 const si32 *yp = y->i32, *cbp = cb->i32, *crp = cr->i32;
803 si32 *rp = r->i32, *gp = g->i32, *bp = b->i32;
804 for (int i = (repeat + 3) >> 2; i > 0; --i)
805 {
806 __m128i my = _mm_load_si128((__m128i*)yp);
807 __m128i mcb = _mm_load_si128((__m128i*)cbp);
808 __m128i mcr = _mm_load_si128((__m128i*)crp);
809
810 __m128i t = _mm_add_epi32(mcb, mcr);
811 t = _mm_sub_epi32(my, _mm_srai_epi32(t, 2));
812 _mm_store_si128((__m128i*)gp, t);
813 __m128i u = _mm_add_epi32(mcb, t);
814 _mm_store_si128((__m128i*)bp, u);
815 u = _mm_add_epi32(mcr, t);
816 _mm_store_si128((__m128i*)rp, u);
817
818 yp += 4; cbp += 4; crp += 4;
819 rp += 4; gp += 4; bp += 4;
820 }
821 }
822 else
823 {
824 assert((y->flags & line_buf::LFT_64BIT) &&
825 (cb->flags & line_buf::LFT_64BIT) &&
826 (cr->flags & line_buf::LFT_64BIT) &&
827 (r->flags & line_buf::LFT_32BIT) &&
828 (g->flags & line_buf::LFT_32BIT) &&
829 (b->flags & line_buf::LFT_32BIT));
830 __m128i v2 = _mm_set1_epi64x(1ULL << (63 - 2));
831 __m128i low_bits = _mm_set_epi64x(0, (si64)ULLONG_MAX);
832 const si64 *yp = y->i64, *cbp = cb->i64, *crp = cr->i64;
833 si32 *rp = r->i32, *gp = g->i32, *bp = b->i32;
834 for (int i = (repeat + 3) >> 2; i > 0; --i)
835 {
836 __m128i my, mcb, mcr, tr, tg, tb;
837 my = _mm_load_si128((__m128i*)yp);
838 mcb = _mm_load_si128((__m128i*)cbp);
839 mcr = _mm_load_si128((__m128i*)crp);
840
841 tg = _mm_add_epi64(mcb, mcr);
842 tg = _mm_sub_epi64(my, sse2_mm_srai_epi64(tg, 2, v2));
843 tb = _mm_add_epi64(mcb, tg);
844 tr = _mm_add_epi64(mcr, tg);
845
846 __m128i mr, mg, mb;
847 mr = _mm_shuffle_epi32(tr, _MM_SHUFFLE(0, 0, 2, 0));
848 mr = _mm_and_si128(low_bits, mr);
849 mg = _mm_shuffle_epi32(tg, _MM_SHUFFLE(0, 0, 2, 0));
850 mg = _mm_and_si128(low_bits, mg);
851 mb = _mm_shuffle_epi32(tb, _MM_SHUFFLE(0, 0, 2, 0));
852 mb = _mm_and_si128(low_bits, mb);
853
854 yp += 2; cbp += 2; crp += 2;
855
856 my = _mm_load_si128((__m128i*)yp);
857 mcb = _mm_load_si128((__m128i*)cbp);
858 mcr = _mm_load_si128((__m128i*)crp);
859
860 tg = _mm_add_epi64(mcb, mcr);
861 tg = _mm_sub_epi64(my, sse2_mm_srai_epi64(tg, 2, v2));
862 tb = _mm_add_epi64(mcb, tg);
863 tr = _mm_add_epi64(mcr, tg);
864
865 tr = _mm_shuffle_epi32(tr, _MM_SHUFFLE(2, 0, 0, 0));
866 tr = _mm_andnot_si128(low_bits, tr);
867 mr = _mm_or_si128(mr, tr);
868 tg = _mm_shuffle_epi32(tg, _MM_SHUFFLE(2, 0, 0, 0));
869 tg = _mm_andnot_si128(low_bits, tg);
870 mg = _mm_or_si128(mg, tg);
871 tb = _mm_shuffle_epi32(tb, _MM_SHUFFLE(2, 0, 0, 0));
872 tb = _mm_andnot_si128(low_bits, tb);
873 mb = _mm_or_si128(mb, tb);
874
875 _mm_store_si128((__m128i*)rp, mr);
876 _mm_store_si128((__m128i*)gp, mg);
877 _mm_store_si128((__m128i*)bp, mb);
878
879 yp += 2; cbp += 2; crp += 2;
880 rp += 4; gp += 4; bp += 4;
881 }
882 }
883 }
884 }
885}
886
887#endif
ojph::line_buf::LFT_INTEGER
@ LFT_INTEGER
Definition ojph_mem.h:170
ojph::line_buf::LFT_32BIT
@ LFT_32BIT
Definition ojph_mem.h:168
ojph::line_buf::LFT_64BIT
@ LFT_64BIT
Definition ojph_mem.h:169
ojph::local::sse2_rct_backward
void sse2_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)
ojph::local::sse2_irv_convert_to_integer
void sse2_irv_convert_to_integer(const line_buf *src_line, line_buf *dst_line, ui32 dst_line_offset, ui32 bit_depth, bool is_signed, ui32 width)
ojph::local::sse2_irv_convert_to_float_nlt_type3
void sse2_irv_convert_to_float_nlt_type3(const line_buf *src_line, ui32 src_line_offset, line_buf *dst_line, ui32 bit_depth, bool is_signed, ui32 width)
ojph::local::sse2_irv_convert_to_integer_nlt_type3
void sse2_irv_convert_to_integer_nlt_type3(const line_buf *src_line, line_buf *dst_line, ui32 dst_line_offset, ui32 bit_depth, bool is_signed, ui32 width)
ojph::local::sse2_rev_convert
void sse2_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)
ojph::local::sse2_rev_convert_nlt_type3
void sse2_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)
ojph::local::sse2_irv_convert_to_integer_nlt
void sse2_irv_convert_to_integer_nlt(const line_buf *src_line, line_buf *dst_line, ui32 dst_line_offset, ui32 bit_depth, bool is_signed, ui32 width, const nlt_rec *rec)
ojph::local::sse2_rct_forward
void sse2_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)
ojph::local::sse2_irv_convert_to_float_nlt
void sse2_irv_convert_to_float_nlt(const line_buf *src_line, ui32 src_line_offset, line_buf *dst_line, ui32 bit_depth, bool is_signed, ui32 width, const nlt_rec *rec)
ojph::local::sse2_irv_convert_to_float
void sse2_irv_convert_to_float(const line_buf *src_line, ui32 src_line_offset, line_buf *dst_line, ui32 bit_depth, bool is_signed, ui32 width)
ojph::si64
int64_t si64
Definition ojph_defs.h:57
ojph::si32
int32_t si32
Definition ojph_defs.h:55
ojph::ui32
uint32_t ui32
Definition ojph_defs.h:54
ojph_arch.h
ojph_colour.h
ojph_defs.h
ojph_unused
#define ojph_unused(x)
Definition ojph_defs.h:78
ojph_mem.h
ojph_params.h
ojph_params_local.h
ojph::local::nlt_rec
Definition ojph_params_local.h:840
ojph::local::nlt_rec::nonlinearity
ojph::param_nlt::nonlinearity nonlinearity
Definition ojph_params_local.h:841