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 fd_min = _mm_set1_ps(rec->fd_min);
223 __m128 fd_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, fd_min);
240 t = _mm_min_ps(t, fd_max);
241 __m128i k = _mm_cvttps_epi32(
242 _mm_mul_ps(_mm_sub_ps(t, fd_min), inv_delta));
243 __m128 d_k = _mm_add_ps(fd_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, fd_min);
275 t = _mm_min_ps(t, fd_max);
276 __m128i k = _mm_cvttps_epi32(
277 _mm_mul_ps(_mm_sub_ps(t, fd_min), inv_delta));
278 __m128 d_k = _mm_add_ps(fd_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 ojph_unused(bit_depth);
592
593 assert((src_line->flags & line_buf::LFT_32BIT) &&
594 (src_line->flags & line_buf::LFT_INTEGER) &&
595 (dst_line->flags & line_buf::LFT_32BIT) &&
596 (dst_line->flags & line_buf::LFT_INTEGER) == 0);
597
598 assert(bit_depth <= 32);
599 __m128 mul = _mm_set1_ps(rec->multiplier);
600 __m128 ft_min = _mm_set1_ps(rec->ft_min);
601 __m128 ft_max = _mm_set1_ps(rec->ft_max);
602 __m128 delta = _mm_set1_ps(rec->delta);
603 __m128 inv_delta = _mm_set1_ps(rec->inv_delta);
604 const float* lut = rec->approx_enc_points;
605
606 __m128 half_ps = _mm_set1_ps(0.5f);
607
608 const si32* sp = src_line->i32 + src_line_offset;
609 float* dp = dst_line->f32;
610 if (rec->is_signed())
611 {
612 __m128i bias =
613 _mm_set1_epi32(-(si32)((1ULL << (rec->get_bit_depth() - 1)) + 1));
614 __m128i zero = _mm_setzero_si128();
615 for (int i = (int)width; i > 0; i -= 4, sp += 4, dp += 4) {
616 __m128i v = _mm_loadu_si128((__m128i*)sp);
617 if (NLT_TYPE == 4)
618 {
619 __m128i c = _mm_cmpgt_epi32(zero, v); // 0xFFFFFFFF for -ve val
620 __m128i neg = _mm_sub_epi32(bias, v); // - bias - value
621 neg = _mm_and_si128(c, neg); // keep only - bias - val
622 v = _mm_andnot_si128(c, v); // keep only +ve or 0
623 v = _mm_or_si128(neg, v); // combine
624 }
625 __m128 t = _mm_add_ps( // convert to [0, 1]
626 _mm_mul_ps(_mm_cvtepi32_ps(v), mul), half_ps);
627 t = _mm_max_ps(t, ft_min);
628 t = _mm_min_ps(t, ft_max);
629 __m128i k = _mm_cvttps_epi32(
630 _mm_mul_ps(_mm_sub_ps(t, ft_min), inv_delta));
631 __m128 d_k = _mm_add_ps(ft_min,
632 _mm_mul_ps(_mm_cvtepi32_ps(k), delta));
633 // SSE2 has no gather; perform the LUT lookup 4 times and build
634 // the vector from the individual results
635 si32 kk[4];
636 _mm_storeu_si128((__m128i*)kk, k);
637 __m128 t_k = _mm_set_ps(lut[kk[3]], lut[kk[2]],
638 lut[kk[1]], lut[kk[0]]);
639 __m128 t_kp1 = _mm_set_ps(lut[kk[3] + 1], lut[kk[2] + 1],
640 lut[kk[1] + 1], lut[kk[0] + 1]);
641 __m128 y = _mm_add_ps(t_k,
642 _mm_mul_ps(_mm_mul_ps(_mm_sub_ps(t, d_k), inv_delta),
643 _mm_sub_ps(t_kp1, t_k)));
644 _mm_storeu_ps(dp, _mm_sub_ps(y, half_ps));
645 }
646 }
647 else
648 {
649 for (int i = (int)width; i > 0; i -= 4, sp += 4, dp += 4) {
650 __m128i v = _mm_loadu_si128((__m128i*)sp);
651 __m128 t = _mm_mul_ps(_mm_cvtepi32_ps(v), mul); // in [0, 1]
652 t = _mm_max_ps(t, ft_min);
653 t = _mm_min_ps(t, ft_max);
654 __m128i k = _mm_cvttps_epi32(
655 _mm_mul_ps(_mm_sub_ps(t, ft_min), inv_delta));
656 __m128 d_k = _mm_add_ps(ft_min,
657 _mm_mul_ps(_mm_cvtepi32_ps(k), delta));
658 // SSE2 has no gather; perform the LUT lookup 4 times and build
659 // the vector from the individual results
660 si32 kk[4];
661 _mm_storeu_si128((__m128i*)kk, k);
662 __m128 t_k = _mm_set_ps(lut[kk[3]], lut[kk[2]],
663 lut[kk[1]], lut[kk[0]]);
664 __m128 t_kp1 = _mm_set_ps(lut[kk[3] + 1], lut[kk[2] + 1],
665 lut[kk[1] + 1], lut[kk[0] + 1]);
666 __m128 y = _mm_add_ps(t_k,
667 _mm_mul_ps(_mm_mul_ps(_mm_sub_ps(t, d_k), inv_delta),
668 _mm_sub_ps(t_kp1, t_k)));
669 _mm_storeu_ps(dp, _mm_sub_ps(y, half_ps));
670 }
671 }
672 }
673
675 // The precise inverse of the look-up table, for a table whose derived
676 // version is fine enough that one step of it crosses at most one entry of
677 // the table that was supplied. The entry that a value falls in is found
678 // with the same arithmetic that the approximation above uses, and the
679 // entry after that is looked at as well, to catch the values that fall
680 // past the entry that was found.
682 template<int NLT_TYPE>
683 static inline
684 void local_sse2_irv_convert_to_float_nlt2or4_precise_1step(
685 const line_buf *src_line, ui32 src_line_offset, line_buf *dst_line,
686 ui32 bit_depth, bool is_signed, ui32 width, const nlt_rec* rec)
687 {
688 ojph_unused(is_signed);
689 ojph_unused(bit_depth);
690
691 assert((src_line->flags & line_buf::LFT_32BIT) &&
692 (src_line->flags & line_buf::LFT_INTEGER) &&
693 (dst_line->flags & line_buf::LFT_32BIT) &&
694 (dst_line->flags & line_buf::LFT_INTEGER) == 0);
695
696 assert(bit_depth <= 32);
697 __m128 mul = _mm_set1_ps(rec->multiplier);
698 __m128 ft_min = _mm_set1_ps(rec->ft_min);
699 __m128 ft_max = _mm_set1_ps(rec->ft_max);
700 __m128 fd_min = _mm_set1_ps(rec->fd_min);
701 __m128 delta = _mm_set1_ps(rec->delta);
702 __m128 inv_delta = _mm_set1_ps(rec->inv_delta);
703 const float* lut = rec->precise_dec_lut;
704 const ui16* indices = rec->precise_enc_dec_indices;
705
706 __m128 half_ps = _mm_set1_ps(0.5f);
707 __m128 one_ps = _mm_set1_ps(1.0f);
708 __m128 zero_ps = _mm_setzero_ps();
709 __m128i one = _mm_set1_epi32(1);
710
711 const si32* sp = src_line->i32 + src_line_offset;
712 float* dp = dst_line->f32;
713 if (rec->is_signed())
714 {
715 __m128i bias =
716 _mm_set1_epi32(-(si32)((1ULL << (rec->get_bit_depth() - 1)) + 1));
717 __m128i zero = _mm_setzero_si128();
718 for (int i = (int)width; i > 0; i -= 4, sp += 4, dp += 4) {
719 __m128i v = _mm_loadu_si128((__m128i*)sp);
720 if (NLT_TYPE == 4)
721 {
722 __m128i c = _mm_cmpgt_epi32(zero, v); // 0xFFFFFFFF for -ve val
723 __m128i neg = _mm_sub_epi32(bias, v); // - bias - value
724 neg = _mm_and_si128(c, neg); // keep only - bias - val
725 v = _mm_andnot_si128(c, v); // keep only +ve or 0
726 v = _mm_or_si128(neg, v); // combine
727 }
728 __m128 t = _mm_add_ps( // convert to [0, 1]
729 _mm_mul_ps(_mm_cvtepi32_ps(v), mul), half_ps);
730 t = _mm_max_ps(t, ft_min);
731 t = _mm_min_ps(t, ft_max);
732 __m128i kk = _mm_cvttps_epi32(
733 _mm_mul_ps(_mm_sub_ps(t, ft_min), inv_delta));
734 // SSE2 has no gather; look the entries up one at a time and build
735 // the vectors from the individual results. The two entries that are
736 // read past the end of the indices, and the two that are read past
737 // the end of the table, are part of their allocations
738 si32 e[4], k[4];
739 _mm_storeu_si128((__m128i*)e, kk);
740 for (int j = 0; j < 4; ++j)
741 k[j] = (si32)indices[e[j]];
742 __m128i kv = _mm_setr_epi32(k[0], k[1], k[2], k[3]);
743 __m128 y0 = _mm_set_ps(lut[k[3]], lut[k[2]], lut[k[1]], lut[k[0]]);
744 __m128 y1 = _mm_set_ps(lut[k[3] + 1], lut[k[2] + 1],
745 lut[k[1] + 1], lut[k[0] + 1]);
746 __m128 y2 = _mm_set_ps(lut[k[3] + 2], lut[k[2] + 2],
747 lut[k[1] + 2], lut[k[0] + 2]);
748 // one step of search: a value that falls past the entry that was
749 // found belongs to the entry after it; SSE2 has no blend, so select
750 // with the mask
751 __m128 adv = _mm_cmpgt_ps(t, y1);
752 y0 = _mm_or_ps(_mm_and_ps(adv, y1), _mm_andnot_ps(adv, y0));
753 y1 = _mm_or_ps(_mm_and_ps(adv, y2), _mm_andnot_ps(adv, y1));
754 kv = _mm_add_epi32(kv,
755 _mm_and_si128(_mm_castps_si128(adv), one));
756 // interpolate, but only over a segment that has a width; the last
757 // entry of the table is repeated, and a segment without a width
758 // contributes nothing
759 __m128 den = _mm_sub_ps(y1, y0);
760 __m128 has = _mm_cmpgt_ps(den, zero_ps);
761 den = _mm_or_ps(_mm_and_ps(has, den), _mm_andnot_ps(has, one_ps));
762 __m128 y = _mm_add_ps(
763 _mm_add_ps(fd_min, _mm_mul_ps(_mm_cvtepi32_ps(kv), delta)),
764 _mm_and_ps(has, _mm_div_ps(
765 _mm_mul_ps(_mm_sub_ps(t, y0), delta), den)));
766 _mm_storeu_ps(dp, _mm_sub_ps(y, half_ps));
767 }
768 }
769 else
770 {
771 for (int i = (int)width; i > 0; i -= 4, sp += 4, dp += 4) {
772 __m128i v = _mm_loadu_si128((__m128i*)sp);
773 __m128 t = _mm_mul_ps(_mm_cvtepi32_ps(v), mul); // in [0, 1]
774 t = _mm_max_ps(t, ft_min);
775 t = _mm_min_ps(t, ft_max);
776 __m128i kk = _mm_cvttps_epi32(
777 _mm_mul_ps(_mm_sub_ps(t, ft_min), inv_delta));
778 // SSE2 has no gather; look the entries up one at a time and build
779 // the vectors from the individual results. The two entries that are
780 // read past the end of the indices, and the two that are read past
781 // the end of the table, are part of their allocations
782 si32 e[4], k[4];
783 _mm_storeu_si128((__m128i*)e, kk);
784 for (int j = 0; j < 4; ++j)
785 k[j] = (si32)indices[e[j]];
786 __m128i kv = _mm_setr_epi32(k[0], k[1], k[2], k[3]);
787 __m128 y0 = _mm_set_ps(lut[k[3]], lut[k[2]], lut[k[1]], lut[k[0]]);
788 __m128 y1 = _mm_set_ps(lut[k[3] + 1], lut[k[2] + 1],
789 lut[k[1] + 1], lut[k[0] + 1]);
790 __m128 y2 = _mm_set_ps(lut[k[3] + 2], lut[k[2] + 2],
791 lut[k[1] + 2], lut[k[0] + 2]);
792 // one step of search: a value that falls past the entry that was
793 // found belongs to the entry after it; SSE2 has no blend, so select
794 // with the mask
795 __m128 adv = _mm_cmpgt_ps(t, y1);
796 y0 = _mm_or_ps(_mm_and_ps(adv, y1), _mm_andnot_ps(adv, y0));
797 y1 = _mm_or_ps(_mm_and_ps(adv, y2), _mm_andnot_ps(adv, y1));
798 kv = _mm_add_epi32(kv,
799 _mm_and_si128(_mm_castps_si128(adv), one));
800 // interpolate, but only over a segment that has a width; the last
801 // entry of the table is repeated, and a segment without a width
802 // contributes nothing
803 __m128 den = _mm_sub_ps(y1, y0);
804 __m128 has = _mm_cmpgt_ps(den, zero_ps);
805 den = _mm_or_ps(_mm_and_ps(has, den), _mm_andnot_ps(has, one_ps));
806 __m128 y = _mm_add_ps(
807 _mm_add_ps(fd_min, _mm_mul_ps(_mm_cvtepi32_ps(kv), delta)),
808 _mm_and_ps(has, _mm_div_ps(
809 _mm_mul_ps(_mm_sub_ps(t, y0), delta), den)));
810 _mm_storeu_ps(dp, _mm_sub_ps(y, half_ps));
811 }
812 }
813 }
814
816 // The precise inverse of the look-up table, for a table whose derived
817 // version is coarse enough that one step of it can cross more than one
818 // entry of the table that was supplied. The entry that a value falls in
819 // is found with a fixed-depth binary search over the entries the step can
820 // reach; the depth is the worst case, ceil(log2(precise_max_steps + 1)),
821 // so every lane performs the same number of steps. SSE2 has no gather,
822 // so the table entries are looked up one at a time and assembled into
823 // vectors. This mirrors local_gen_irv_convert_to_float_nlt2or4_precise_search
824 // in ojph_colour.cpp, so keep the two in step when that one changes.
826 template<int NLT_TYPE>
827 static inline
828 void local_sse2_irv_convert_to_float_nlt2or4_precise_search(
829 const line_buf *src_line, ui32 src_line_offset, line_buf *dst_line,
830 ui32 bit_depth, bool is_signed, ui32 width, const nlt_rec* rec)
831 {
832 ojph_unused(is_signed);
833 ojph_unused(bit_depth);
834
835 assert((src_line->flags & line_buf::LFT_32BIT) &&
836 (src_line->flags & line_buf::LFT_INTEGER) &&
837 (dst_line->flags & line_buf::LFT_32BIT) &&
838 (dst_line->flags & line_buf::LFT_INTEGER) == 0);
839
840 assert(bit_depth <= 32);
841 assert(rec->precise_max_steps > 1);
842 __m128 mul = _mm_set1_ps(rec->multiplier);
843 __m128 ft_min = _mm_set1_ps(rec->ft_min);
844 __m128 ft_max = _mm_set1_ps(rec->ft_max);
845 __m128 fd_min = _mm_set1_ps(rec->fd_min);
846 __m128 delta = _mm_set1_ps(rec->delta);
847 __m128 inv_delta = _mm_set1_ps(rec->inv_delta);
848 const float* lut = rec->precise_dec_lut;
849 const ui16* indices = rec->precise_enc_dec_indices;
850 ui32 num_points = rec->num_points;
851
852 // the search interval is at most precise_max_steps + 1 entries wide;
853 // run the worst-case number of binary-search steps on every lane,
854 // which is ceil(log2(span))
855 ui32 span = rec->precise_max_steps + 1;
856 int steps = 32 - (int)count_leading_zeros(rec->precise_max_steps);
857
858 __m128 half_ps = _mm_set1_ps(0.5f);
859 __m128 one_ps = _mm_set1_ps(1.0f);
860 __m128 zero_ps = _mm_setzero_ps();
861 __m128i span_epi32 = _mm_set1_epi32((si32)span);
862 __m128i num_epi32 = _mm_set1_epi32((si32)num_points);
863
864 const si32* sp = src_line->i32 + src_line_offset;
865 float* dp = dst_line->f32;
866 if (rec->is_signed())
867 {
868 __m128i bias =
869 _mm_set1_epi32(-(si32)((1ULL << (rec->get_bit_depth() - 1)) + 1));
870 __m128i zero = _mm_setzero_si128();
871 for (int i = (int)width; i > 0; i -= 4, sp += 4, dp += 4) {
872 __m128i v = _mm_loadu_si128((__m128i*)sp);
873 if (NLT_TYPE == 4)
874 {
875 __m128i c = _mm_cmpgt_epi32(zero, v); // 0xFFFFFFFF for -ve val
876 __m128i neg = _mm_sub_epi32(bias, v); // - bias - value
877 neg = _mm_and_si128(c, neg); // keep only - bias - val
878 v = _mm_andnot_si128(c, v); // keep only +ve or 0
879 v = _mm_or_si128(neg, v); // combine
880 }
881 __m128 t = _mm_add_ps( // convert to [0, 1]
882 _mm_mul_ps(_mm_cvtepi32_ps(v), mul), half_ps);
883 t = _mm_max_ps(t, ft_min);
884 t = _mm_min_ps(t, ft_max);
885 __m128i kk = _mm_cvttps_epi32(
886 _mm_mul_ps(_mm_sub_ps(t, ft_min), inv_delta));
887 // SSE2 has no gather; look the indices up one at a time
888 si32 e[4], k[4];
889 _mm_storeu_si128((__m128i*)e, kk);
890 for (int j = 0; j < 4; ++j)
891 k[j] = (si32)indices[e[j]];
892 __m128i lo = _mm_setr_epi32(k[0], k[1], k[2], k[3]);
893 // hi = min(k + span, num_points); SSE2 has no 32-bit min, so
894 // select with a signed compare
895 __m128i kspan = _mm_add_epi32(lo, span_epi32);
896 __m128i m = _mm_cmpgt_epi32(kspan, num_epi32);
897 __m128i hi = _mm_or_si128(_mm_and_si128(m, num_epi32),
898 _mm_andnot_si128(m, kspan));
899
900 for (int step = 0; step < steps; ++step)
901 {
902 __m128i mid = _mm_srli_epi32(_mm_add_epi32(lo, hi), 1);
903 si32 m4[4];
904 _mm_storeu_si128((__m128i*)m4, mid);
905 __m128 lmid = _mm_set_ps(lut[m4[3]], lut[m4[2]],
906 lut[m4[1]], lut[m4[0]]);
907 __m128 ge = _mm_cmpge_ps(t, lmid);
908 __m128i g = _mm_castps_si128(ge);
909 lo = _mm_or_si128(_mm_and_si128(g, mid),
910 _mm_andnot_si128(g, lo));
911 hi = _mm_or_si128(_mm_and_si128(g, hi),
912 _mm_andnot_si128(g, mid));
913 }
914
915 si32 l4[4], h4[4];
916 _mm_storeu_si128((__m128i*)l4, lo);
917 _mm_storeu_si128((__m128i*)h4, hi);
918 __m128 y0 = _mm_set_ps(lut[l4[3]], lut[l4[2]], lut[l4[1]], lut[l4[0]]);
919 __m128 y1 = _mm_set_ps(lut[h4[3]], lut[h4[2]], lut[h4[1]], lut[h4[0]]);
920 // interpolate, but only over a segment that has a width; the last
921 // entry of the table is repeated, and a segment without a width
922 // contributes nothing
923 __m128 den = _mm_sub_ps(y1, y0);
924 __m128 has = _mm_cmpgt_ps(den, zero_ps);
925 den = _mm_or_ps(_mm_and_ps(has, den), _mm_andnot_ps(has, one_ps));
926 __m128 y = _mm_add_ps(
927 _mm_add_ps(fd_min, _mm_mul_ps(_mm_cvtepi32_ps(lo), delta)),
928 _mm_and_ps(has, _mm_div_ps(
929 _mm_mul_ps(_mm_sub_ps(t, y0), delta), den)));
930 _mm_storeu_ps(dp, _mm_sub_ps(y, half_ps));
931 }
932 }
933 else
934 {
935 for (int i = (int)width; i > 0; i -= 4, sp += 4, dp += 4) {
936 __m128i v = _mm_loadu_si128((__m128i*)sp);
937 __m128 t = _mm_mul_ps(_mm_cvtepi32_ps(v), mul); // in [0, 1]
938 t = _mm_max_ps(t, ft_min);
939 t = _mm_min_ps(t, ft_max);
940 __m128i kk = _mm_cvttps_epi32(
941 _mm_mul_ps(_mm_sub_ps(t, ft_min), inv_delta));
942 // SSE2 has no gather; look the indices up one at a time
943 si32 e[4], k[4];
944 _mm_storeu_si128((__m128i*)e, kk);
945 for (int j = 0; j < 4; ++j)
946 k[j] = (si32)indices[e[j]];
947 __m128i lo = _mm_setr_epi32(k[0], k[1], k[2], k[3]);
948 __m128i kspan = _mm_add_epi32(lo, span_epi32);
949 __m128i m = _mm_cmpgt_epi32(kspan, num_epi32);
950 __m128i hi = _mm_or_si128(_mm_and_si128(m, num_epi32),
951 _mm_andnot_si128(m, kspan));
952
953 for (int step = 0; step < steps; ++step)
954 {
955 __m128i mid = _mm_srli_epi32(_mm_add_epi32(lo, hi), 1);
956 si32 m4[4];
957 _mm_storeu_si128((__m128i*)m4, mid);
958 __m128 lmid = _mm_set_ps(lut[m4[3]], lut[m4[2]],
959 lut[m4[1]], lut[m4[0]]);
960 __m128 ge = _mm_cmpge_ps(t, lmid);
961 __m128i g = _mm_castps_si128(ge);
962 lo = _mm_or_si128(_mm_and_si128(g, mid),
963 _mm_andnot_si128(g, lo));
964 hi = _mm_or_si128(_mm_and_si128(g, hi),
965 _mm_andnot_si128(g, mid));
966 }
967
968 si32 l4[4], h4[4];
969 _mm_storeu_si128((__m128i*)l4, lo);
970 _mm_storeu_si128((__m128i*)h4, hi);
971 __m128 y0 = _mm_set_ps(lut[l4[3]], lut[l4[2]], lut[l4[1]], lut[l4[0]]);
972 __m128 y1 = _mm_set_ps(lut[h4[3]], lut[h4[2]], lut[h4[1]], lut[h4[0]]);
973 __m128 den = _mm_sub_ps(y1, y0);
974 __m128 has = _mm_cmpgt_ps(den, zero_ps);
975 den = _mm_or_ps(_mm_and_ps(has, den), _mm_andnot_ps(has, one_ps));
976 __m128 y = _mm_add_ps(
977 _mm_add_ps(fd_min, _mm_mul_ps(_mm_cvtepi32_ps(lo), delta)),
978 _mm_and_ps(has, _mm_div_ps(
979 _mm_mul_ps(_mm_sub_ps(t, y0), delta), den)));
980 _mm_storeu_ps(dp, _mm_sub_ps(y, half_ps));
981 }
982 }
983 }
984
986 void sse2_irv_convert_to_float_nlt(const line_buf *src_line,
987 ui32 src_line_offset, line_buf *dst_line,
988 ui32 bit_depth, bool is_signed, ui32 width, const nlt_rec* rec)
989 {
990 using nl = nlt_rec::nonlinearity;
991 if (rec->get_type() == nl::OJPH_NLT_LUT_STYLE_NLT)
992 {
993 if (rec->approx_encoding_nlt)
994 local_sse2_irv_convert_to_float_nlt2or4<2>(src_line,
995 src_line_offset, dst_line, bit_depth, is_signed, width, rec);
996 else if (rec->precise_max_steps == 1)
997 local_sse2_irv_convert_to_float_nlt2or4_precise_1step<2>(src_line,
998 src_line_offset, dst_line, bit_depth, is_signed, width, rec);
999 else
1000 local_sse2_irv_convert_to_float_nlt2or4_precise_search<2>(src_line,
1001 src_line_offset, dst_line, bit_depth, is_signed, width, rec);
1002 }
1003 else if (rec->get_type() == nl::OJPH_NLT_BINARY_COMPLEMENT_PLUS_LUT)
1004 {
1005 if (rec->approx_encoding_nlt)
1006 local_sse2_irv_convert_to_float_nlt2or4<4>(src_line,
1007 src_line_offset, dst_line, bit_depth, is_signed, width, rec);
1008 else if (rec->precise_max_steps == 1)
1009 local_sse2_irv_convert_to_float_nlt2or4_precise_1step<4>(src_line,
1010 src_line_offset, dst_line, bit_depth, is_signed, width, rec);
1011 else
1012 local_sse2_irv_convert_to_float_nlt2or4_precise_search<4>(src_line,
1013 src_line_offset, dst_line, bit_depth, is_signed, width, rec);
1014 }
1015 else
1016 assert(0);
1017 }
1018
1020 void sse2_rct_forward(const line_buf *r,
1021 const line_buf *g,
1022 const line_buf *b,
1023 line_buf *y, line_buf *cb, line_buf *cr,
1024 ui32 repeat)
1025 {
1026 assert((y->flags & line_buf::LFT_INTEGER) &&
1027 (cb->flags & line_buf::LFT_INTEGER) &&
1028 (cr->flags & line_buf::LFT_INTEGER) &&
1029 (r->flags & line_buf::LFT_INTEGER) &&
1030 (g->flags & line_buf::LFT_INTEGER) &&
1031 (b->flags & line_buf::LFT_INTEGER));
1032
1033 if (y->flags & line_buf::LFT_32BIT)
1034 {
1035 assert((y->flags & line_buf::LFT_32BIT) &&
1036 (cb->flags & line_buf::LFT_32BIT) &&
1037 (cr->flags & line_buf::LFT_32BIT) &&
1038 (r->flags & line_buf::LFT_32BIT) &&
1039 (g->flags & line_buf::LFT_32BIT) &&
1040 (b->flags & line_buf::LFT_32BIT));
1041 const si32 *rp = r->i32, * gp = g->i32, * bp = b->i32;
1042 si32 *yp = y->i32, * cbp = cb->i32, * crp = cr->i32;
1043 for (int i = (repeat + 3) >> 2; i > 0; --i)
1044 {
1045 __m128i mr = _mm_load_si128((__m128i*)rp);
1046 __m128i mg = _mm_load_si128((__m128i*)gp);
1047 __m128i mb = _mm_load_si128((__m128i*)bp);
1048 __m128i t = _mm_add_epi32(mr, mb);
1049 t = _mm_add_epi32(t, _mm_slli_epi32(mg, 1));
1050 _mm_store_si128((__m128i*)yp, _mm_srai_epi32(t, 2));
1051 t = _mm_sub_epi32(mb, mg);
1052 _mm_store_si128((__m128i*)cbp, t);
1053 t = _mm_sub_epi32(mr, mg);
1054 _mm_store_si128((__m128i*)crp, t);
1055
1056 rp += 4; gp += 4; bp += 4;
1057 yp += 4; cbp += 4; crp += 4;
1058 }
1059 }
1060 else
1061 {
1062 assert((y->flags & line_buf::LFT_64BIT) &&
1063 (cb->flags & line_buf::LFT_64BIT) &&
1064 (cr->flags & line_buf::LFT_64BIT) &&
1065 (r->flags & line_buf::LFT_32BIT) &&
1066 (g->flags & line_buf::LFT_32BIT) &&
1067 (b->flags & line_buf::LFT_32BIT));
1068 __m128i zero = _mm_setzero_si128();
1069 __m128i v2 = _mm_set1_epi64x(1ULL << (63 - 2));
1070 const si32 *rp = r->i32, *gp = g->i32, *bp = b->i32;
1071 si64 *yp = y->i64, *cbp = cb->i64, *crp = cr->i64;
1072 for (int i = (repeat + 3) >> 2; i > 0; --i)
1073 {
1074 __m128i mr32 = _mm_load_si128((__m128i*)rp);
1075 __m128i mg32 = _mm_load_si128((__m128i*)gp);
1076 __m128i mb32 = _mm_load_si128((__m128i*)bp);
1077 __m128i mr, mg, mb, t;
1078 mr = sse2_cvtlo_epi32_epi64(mr32, zero);
1079 mg = sse2_cvtlo_epi32_epi64(mg32, zero);
1080 mb = sse2_cvtlo_epi32_epi64(mb32, zero);
1081
1082 t = _mm_add_epi64(mr, mb);
1083 t = _mm_add_epi64(t, _mm_slli_epi64(mg, 1));
1084 _mm_store_si128((__m128i*)yp, sse2_mm_srai_epi64(t, 2, v2));
1085 t = _mm_sub_epi64(mb, mg);
1086 _mm_store_si128((__m128i*)cbp, t);
1087 t = _mm_sub_epi64(mr, mg);
1088 _mm_store_si128((__m128i*)crp, t);
1089
1090 yp += 2; cbp += 2; crp += 2;
1091
1092 mr = sse2_cvthi_epi32_epi64(mr32, zero);
1093 mg = sse2_cvthi_epi32_epi64(mg32, zero);
1094 mb = sse2_cvthi_epi32_epi64(mb32, zero);
1095
1096 t = _mm_add_epi64(mr, mb);
1097 t = _mm_add_epi64(t, _mm_slli_epi64(mg, 1));
1098 _mm_store_si128((__m128i*)yp, sse2_mm_srai_epi64(t, 2, v2));
1099 t = _mm_sub_epi64(mb, mg);
1100 _mm_store_si128((__m128i*)cbp, t);
1101 t = _mm_sub_epi64(mr, mg);
1102 _mm_store_si128((__m128i*)crp, t);
1103
1104 rp += 4; gp += 4; bp += 4;
1105 yp += 2; cbp += 2; crp += 2;
1106 }
1107 }
1108 }
1109
1111 void sse2_rct_backward(const line_buf *y,
1112 const line_buf *cb,
1113 const line_buf *cr,
1114 line_buf *r, line_buf *g, line_buf *b,
1115 ui32 repeat)
1116 {
1117 assert((y->flags & line_buf::LFT_INTEGER) &&
1118 (cb->flags & line_buf::LFT_INTEGER) &&
1119 (cr->flags & line_buf::LFT_INTEGER) &&
1120 (r->flags & line_buf::LFT_INTEGER) &&
1121 (g->flags & line_buf::LFT_INTEGER) &&
1122 (b->flags & line_buf::LFT_INTEGER));
1123
1124 if (y->flags & line_buf::LFT_32BIT)
1125 {
1126 assert((y->flags & line_buf::LFT_32BIT) &&
1127 (cb->flags & line_buf::LFT_32BIT) &&
1128 (cr->flags & line_buf::LFT_32BIT) &&
1129 (r->flags & line_buf::LFT_32BIT) &&
1130 (g->flags & line_buf::LFT_32BIT) &&
1131 (b->flags & line_buf::LFT_32BIT));
1132 const si32 *yp = y->i32, *cbp = cb->i32, *crp = cr->i32;
1133 si32 *rp = r->i32, *gp = g->i32, *bp = b->i32;
1134 for (int i = (repeat + 3) >> 2; i > 0; --i)
1135 {
1136 __m128i my = _mm_load_si128((__m128i*)yp);
1137 __m128i mcb = _mm_load_si128((__m128i*)cbp);
1138 __m128i mcr = _mm_load_si128((__m128i*)crp);
1139
1140 __m128i t = _mm_add_epi32(mcb, mcr);
1141 t = _mm_sub_epi32(my, _mm_srai_epi32(t, 2));
1142 _mm_store_si128((__m128i*)gp, t);
1143 __m128i u = _mm_add_epi32(mcb, t);
1144 _mm_store_si128((__m128i*)bp, u);
1145 u = _mm_add_epi32(mcr, t);
1146 _mm_store_si128((__m128i*)rp, u);
1147
1148 yp += 4; cbp += 4; crp += 4;
1149 rp += 4; gp += 4; bp += 4;
1150 }
1151 }
1152 else
1153 {
1154 assert((y->flags & line_buf::LFT_64BIT) &&
1155 (cb->flags & line_buf::LFT_64BIT) &&
1156 (cr->flags & line_buf::LFT_64BIT) &&
1157 (r->flags & line_buf::LFT_32BIT) &&
1158 (g->flags & line_buf::LFT_32BIT) &&
1159 (b->flags & line_buf::LFT_32BIT));
1160 __m128i v2 = _mm_set1_epi64x(1ULL << (63 - 2));
1161 __m128i low_bits = _mm_set_epi64x(0, (si64)ULLONG_MAX);
1162 const si64 *yp = y->i64, *cbp = cb->i64, *crp = cr->i64;
1163 si32 *rp = r->i32, *gp = g->i32, *bp = b->i32;
1164 for (int i = (repeat + 3) >> 2; i > 0; --i)
1165 {
1166 __m128i my, mcb, mcr, tr, tg, tb;
1167 my = _mm_load_si128((__m128i*)yp);
1168 mcb = _mm_load_si128((__m128i*)cbp);
1169 mcr = _mm_load_si128((__m128i*)crp);
1170
1171 tg = _mm_add_epi64(mcb, mcr);
1172 tg = _mm_sub_epi64(my, sse2_mm_srai_epi64(tg, 2, v2));
1173 tb = _mm_add_epi64(mcb, tg);
1174 tr = _mm_add_epi64(mcr, tg);
1175
1176 __m128i mr, mg, mb;
1177 mr = _mm_shuffle_epi32(tr, _MM_SHUFFLE(0, 0, 2, 0));
1178 mr = _mm_and_si128(low_bits, mr);
1179 mg = _mm_shuffle_epi32(tg, _MM_SHUFFLE(0, 0, 2, 0));
1180 mg = _mm_and_si128(low_bits, mg);
1181 mb = _mm_shuffle_epi32(tb, _MM_SHUFFLE(0, 0, 2, 0));
1182 mb = _mm_and_si128(low_bits, mb);
1183
1184 yp += 2; cbp += 2; crp += 2;
1185
1186 my = _mm_load_si128((__m128i*)yp);
1187 mcb = _mm_load_si128((__m128i*)cbp);
1188 mcr = _mm_load_si128((__m128i*)crp);
1189
1190 tg = _mm_add_epi64(mcb, mcr);
1191 tg = _mm_sub_epi64(my, sse2_mm_srai_epi64(tg, 2, v2));
1192 tb = _mm_add_epi64(mcb, tg);
1193 tr = _mm_add_epi64(mcr, tg);
1194
1195 tr = _mm_shuffle_epi32(tr, _MM_SHUFFLE(2, 0, 0, 0));
1196 tr = _mm_andnot_si128(low_bits, tr);
1197 mr = _mm_or_si128(mr, tr);
1198 tg = _mm_shuffle_epi32(tg, _MM_SHUFFLE(2, 0, 0, 0));
1199 tg = _mm_andnot_si128(low_bits, tg);
1200 mg = _mm_or_si128(mg, tg);
1201 tb = _mm_shuffle_epi32(tb, _MM_SHUFFLE(2, 0, 0, 0));
1202 tb = _mm_andnot_si128(low_bits, tb);
1203 mb = _mm_or_si128(mb, tb);
1204
1205 _mm_store_si128((__m128i*)rp, mr);
1206 _mm_store_si128((__m128i*)gp, mg);
1207 _mm_store_si128((__m128i*)bp, mb);
1208
1209 yp += 2; cbp += 2; crp += 2;
1210 rp += 4; gp += 4; bp += 4;
1211 }
1212 }
1213 }
1214 }
1215}
1216
1217#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::ui16
uint16_t ui16
Definition ojph_defs.h:52
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:852
ojph::local::nlt_rec::nonlinearity
ojph::param_nlt::nonlinearity nonlinearity
Definition ojph_params_local.h:853