yuvdenoise performance patch
Franz Brauße <[email protected]> Fri, 1 Oct 2010 15:47:46 +0200
| Newsgroups | gmane.comp.video.mjpeg.devel |
|---|---|
| Message-ID | <[email protected]> |
This is a multi-part message in MIME format. --Multipart=_Fri__1_Oct_2010_15_47_46_+0200_3M22tzehqubz742L Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: quoted-printable Hello there, I've sent a previous mail already with an earlier version of the patch. Maybe someone still monitors the list... I got some of my movies denoised pretty nicely using yuvdenoise. The quality is ok, my only concern was the performance. Attached is a patch which contains SSE2-accelerated versions of the (non MC-) functions for temporal and spatial filtering (which I mainly use). Additionally I've reenabled the shortcircuiting of temporal_filter_planes_MC, otherwise it fails with divide by zero when using level 0 (e.g. to only filter the luma plane). The temporal filter function now processes a block of 14 pixels and the spatial function does 4 pixels at a time, using the effect that adjacent pixels share many of their neighbours that now only need to be examined once. Other than that, the computations are duplicated from the original functions. With this revised patch, temporal filtering runs about 8 times as fast and spatial filtering 4 to 5 times with this patch, at least on 64-bit machines; tested on a recent Xeon and Opteron processor. Now more than realtime-filtering is possible on my machine. A 32-bit Pentium M doesn't seem to like my SSE2-version of the spatial filtering, I may however find the time to optimize the function a bit further. For now it is not enabled on non-x86-64 processors. Since there are some rounding issues to be taken into account here, I've added a define "OLD_ROUNDING", which produces the exact same output as the old filter (I've checked the denoise-results for errors using md5sum). If it's not set, the filtering should be a bit faster and more accurate, since floating point rounding towards nearest is used then. Thank you all for your work on mjpegtools, I hope you may find the patch worth incorporating. Regards, Franz Brau=C3=9Fe --Multipart=_Fri__1_Oct_2010_15_47_46_+0200_3M22tzehqubz742L Content-Type: text/x-patch; name="yuvdenoise-sse2-2.patch" Content-Disposition: attachment; filename="yuvdenoise-sse2-2.patch" Content-Transfer-Encoding: quoted-printable diff -urw a/yuvdenoise/main.c b/yuvdenoise/main.c --- a/yuvdenoise/main.c 2007-04-02 17:43:35.000000000 +0200 +++ b/yuvdenoise/main.c 2010-10-01 07:03:16.211611711 +0200 @@ -12,6 +12,11 @@ * * ***********************************************************/ =20 +/* 2010-09-22, Franz Brau=C3=9Fe <[email protected]> + * - added SSE2-accelerated versions of filter_plane_median() + * and temporal_filter_planes() + */ + #include <stdio.h> #include <stdlib.h> #include <string.h> @@ -23,6 +28,12 @@ #include "cpu_accel.h" #include "motionsearch.h" =20 +#if defined(__SSE3__) +# include <pmmintrin.h> +#elif defined(__SSE2__) +# include <emmintrin.h> +#endif + int verbose =3D 1; int width =3D 0; int height =3D 0; @@ -76,6 +87,9 @@ * helper-functions * ***********************************************************/ =20 +static void (*filter_plane_median)(uint8_t *, int, int, int); +static void (*temporal_filter_planes)(int, int, int, int); + =20 void gauss_filter_plane (uint8_t * frame, int w, int h, int t) @@ -129,154 +143,6 @@ } =20 void -temporal_filter_planes (int idx, int w, int h, int t) -{ - uint32_t r, c, m; - int32_t d; - int x; - - uint8_t *f1 =3D frame1[idx]; - uint8_t *f2 =3D frame2[idx]; - uint8_t *f3 =3D frame3[idx]; - uint8_t *f4 =3D frame4[idx]; - uint8_t *f5 =3D frame5[idx]; - uint8_t *f6 =3D frame6[idx]; - uint8_t *f7 =3D frame7[idx]; - uint8_t *of =3D outframe[idx]; - - if (t =3D=3D 0) // shortcircuit filter if t =3D 0... - { - memcpy (of, f4, w * h); - return; - } - - for (x =3D 0; x < (w * h); x++) - { - - r =3D *(f4-1-w); - r +=3D *(f4 -w)*2; - r +=3D *(f4+1-w); - r +=3D *(f4-1 )*2; - r +=3D *(f4 )*4; - r +=3D *(f4+1 )*2; - r +=3D *(f4-1+w); - r +=3D *(f4 +w)*2; - r +=3D *(f4+1+w); - r /=3D 16; - - m =3D *(f4)*(t+1)*2; - c =3D t+1; - - d =3D *(f3-1-w); - d +=3D *(f3 -w)*2; - d +=3D *(f3+1-w); - d +=3D *(f3-1 )*2; - d +=3D *(f3 )*4; - d +=3D *(f3+1 )*2; - d +=3D *(f3-1+w); - d +=3D *(f3 +w)*2; - d +=3D *(f3+1+w); - d /=3D 16; - - d =3D t - abs (r-d); - d =3D d<0? 0:d; - c +=3D d; - m +=3D *(f3)*d*2; - - d =3D *(f2-1-w); - d +=3D *(f2 -w)*2; - d +=3D *(f2+1-w); - d +=3D *(f2-1 )*2; - d +=3D *(f2 )*4; - d +=3D *(f2+1 )*2; - d +=3D *(f2-1+w); - d +=3D *(f2 +w)*2; - d +=3D *(f2+1+w); - d /=3D 16; - - d =3D t - abs (r-d); - d =3D d<0? 0:d; - c +=3D d; - m +=3D *(f2)*d*2; - - d =3D *(f1-1-w); - d +=3D *(f1 -w)*2; - d +=3D *(f1+1-w); - d +=3D *(f1-1 )*2; - d +=3D *(f1 )*4; - d +=3D *(f1+1 )*2; - d +=3D *(f1-1+w); - d +=3D *(f1 +w)*2; - d +=3D *(f1+1+w); - d /=3D 16; - - d =3D t - abs (r-d); - d =3D d<0? 0:d; - c +=3D d; - m +=3D *(f1)*d*2; - - d =3D *(f5-1-w); - d +=3D *(f5 -w)*2; - d +=3D *(f5+1-w); - d +=3D *(f5-1 )*2; - d +=3D *(f5 )*4; - d +=3D *(f5+1 )*2; - d +=3D *(f5-1+w); - d +=3D *(f5 +w)*2; - d +=3D *(f5+1+w); - d /=3D 16; - - d =3D t - abs (r-d); - d =3D d<0? 0:d; - c +=3D d; - m +=3D *(f5)*d*2; - - d =3D *(f6-1-w); - d +=3D *(f6 -w)*2; - d +=3D *(f6+1-w); - d +=3D *(f6-1 )*2; - d +=3D *(f6 )*4; - d +=3D *(f6+1 )*2; - d +=3D *(f6-1+w); - d +=3D *(f6 +w)*2; - d +=3D *(f6+1+w); - d /=3D 16; - - d =3D t - abs (r-d); - d =3D d<0? 0:d; - c +=3D d; - m +=3D *(f6)*d*2; - - d =3D *(f7-1-w); - d +=3D *(f7 -w)*2; - d +=3D *(f7+1-w); - d +=3D *(f7-1 )*2; - d +=3D *(f7 )*4; - d +=3D *(f7+1 )*2; - d +=3D *(f7-1+w); - d +=3D *(f7 +w)*2; - d +=3D *(f7+1+w); - d /=3D 16; - - d =3D t - abs (r-d); - d =3D d<0? 0:d; - c +=3D d; - m +=3D *(f7)*d*2; - - *(of) =3D ((m/c)+1)/2; - - f1++; - f2++; - f3++; - f4++; - f5++; - f6++; - f7++; - of++; - } -} - -void temporal_filter_planes_MC (int idx, int w, int h, int t) { uint32_t sad,min; @@ -301,7 +167,7 @@ uint8_t *f7 =3D frame7[idx]; uint8_t *of =3D outframe[idx]; =20 -#if 0 +#if 1 =20 if (t =3D=3D 0) // shortcircuit filter if t =3D 0... { @@ -581,8 +447,640 @@ *(frame+i)=3D(*(frame+i)*(255-level)+random[i&8191]*level)/255; } =20 +#if defined(__SSE2__) + +static inline __m128i tf0(const __m128i mask, const __m128i l0, const __m1= 28i vt, const __m128i vc, const __m128i vb) { + __m128i k0, k1, k2, k3, d0; /* temp storage, pixel surroundings, 16-bit w= ords */ +=09 + /* even pixels */ + /* extract and add the pixels above and below the current one */ + k0 =3D _mm_and_si128(_mm_srli_si128(vt, 1), mask); + k1 =3D _mm_and_si128(_mm_srli_si128(vb, 1), mask); + k0 =3D _mm_add_epi16(k0, k1); +=09 + /* add together the 4 corner pixels diagonal of the current one */ + k2 =3D _mm_add_epi16(_mm_and_si128(vt, mask), _mm_and_si128(vb, mask)); + k2 =3D _mm_add_epi16(k2, _mm_srli_si128(k2, 2)); +=09 + /* add pixels left and right of the current */ + k3 =3D _mm_and_si128(vc, mask); + k3 =3D _mm_add_epi16(k3, _mm_srli_si128(k3, 2)); +=09 + /* add weighted current pixel and the above results */ + k1 =3D _mm_and_si128(_mm_srli_si128(vc, 1), mask); + d0 =3D _mm_slli_epi16(k1, 1); /* center * 4 */ + d0 =3D _mm_add_epi16(d0, k0); + d0 =3D _mm_add_epi16(d0, k3); + d0 =3D _mm_slli_epi16(d0, 1); /* + above,below,left,right * 2 */ + d0 =3D _mm_add_epi16(d0, k2); /* + diagonal * 1 */ + d0 =3D _mm_srli_epi16(d0, 4); /* all / 16 */ + return d0; +} + +static inline __m128i tf1(const __m128i mask, const __m128i l0, const __m1= 28i vt, const __m128i vc, const __m128i vb) { + __m128i k0, k1, k2, k3, d1; +=09 + k0 =3D _mm_srli_si128(_mm_add_epi16(_mm_and_si128(vt, mask), _mm_and_si12= 8(vb, mask)), 2); +=09 + k1 =3D _mm_and_si128(_mm_srli_si128(vt, 1), mask); + k2 =3D _mm_and_si128(_mm_srli_si128(vb, 1), mask); + k2 =3D _mm_add_epi16(k1, k2); + k2 =3D _mm_add_epi16(k2, _mm_srli_si128(k2, 2)); +=09 + k3 =3D _mm_and_si128(_mm_srli_si128(vc, 1), mask); + k3 =3D _mm_add_epi16(k3, _mm_srli_si128(k3, 2)); +=09 + k1 =3D _mm_and_si128(_mm_srli_si128(vc, 2), mask); + d1 =3D _mm_slli_epi16(k1, 1); + d1 =3D _mm_add_epi16(d1, k0); + d1 =3D _mm_add_epi16(d1, k3); + d1 =3D _mm_slli_epi16(d1, 1); + d1 =3D _mm_add_epi16(d1, k2); + d1 =3D _mm_srli_epi16(d1, 4); + return d1; +} + +/* 8 times as fast on x86_64, 2.2 times as fast on i686 */ +void temporal_filter_planes_sse2(int idx, int w, int h, int t) +{ + int x, k; +=09 + uint8_t *f4 =3D frame4[idx]; + uint8_t *of =3D outframe[idx]; +=09 + uint8_t *f[6] =3D { + frame3[idx], frame2[idx], frame1[idx], frame5[idx], frame6[idx], frame7[= idx] + }; +=09 + if (t =3D=3D 0) // shortcircuit filter if t =3D 0... + { + memcpy (of, f4, w * h); + return; + } +=09 + /* vt: x x x x x x x x x x x x x x x x + * vc: x 0 1 2 3 4 5 6 7 8 9 a b c d x + * vb: x x x x x x x x x x x x x x x x + * + * c0, d0, m0 and m1 store the respective values for even pixels, + * m0 for 0, 2, 4 and 6, m1 for 8, a and c in its lower dwords; + * c1, d1, m2 and m3 store the values for the odd pixels, + * m2 for 1, 3, 5 and 7, m3 for 9, b and d in its lower dwords; + * whereas computation for each pixel is as follows (equivalent to the + * non-SSE-accelerated variant of this function): + * + * for each pixel k1 of the 14 pixels per block: + * vt: k2 k0 k2 + * vc: k3 k1 k3 + * vb: k2 k0 k2 + * + * r =3D *f4++; + * c =3D t + 1; + * m =3D r * (t+1); + * for (i=3D0; i<6; i++) { + * d =3D sum(k2) + 2 * (sum(k0) + sum(k3) + 2 * k1) + * d =3D saturate(t - abs(r-d)); + * c +=3D d; + * m +=3D *f[i]++ * d; + * } + * m *=3D 2; + * *of++ =3D (m / c + 1) / 2; + * + * For each frame, first the 7 even pixels are being processed, then the 7 + * odd ones. + */ +=09 + __m128i vt, vc, vb; /* top-, center- and bottom-line of 3x16 block */ + __m128i c0, c1, m0, m1, m2, m3; /* c: 16-bit words, m: 32-bit dwords */ + __m128i d0, d1, r0, r1; /* 16-bit words, 0: even, 1: odd */ + const __m128i mask =3D _mm_set1_epi16(0x00ff); + const __m128i l0 =3D _mm_set1_epi16(t); +=09 +#ifndef OLD_ROUNDING + _MM_SET_ROUNDING_MODE(_MM_ROUND_NEAREST); +#endif +=09 + for (x =3D 0; x < (w * h); x+=3D14) + { + vt =3D _mm_loadu_si128((__m128i *)(f4 - 1 - w)); + vc =3D _mm_loadu_si128((__m128i *)(f4 - 1 )); + vb =3D _mm_loadu_si128((__m128i *)(f4 - 1 + w)); + f4 +=3D 14; + =09 + r0 =3D tf0(mask, l0, vt, vc, vb); /* even pixels */ + r1 =3D tf1(mask, l0, vt, vc, vb); /* odd pixels */ + =09 + c0 =3D c1 =3D _mm_set1_epi16(t + 1); + =09 + /* m =3D *f4 * (t+1) */ + /* The low 16-bit of the multiplication suffice, because both operands a= re + * only 8-bit-values */ + d0 =3D _mm_mullo_epi16(_mm_and_si128(_mm_srli_si128(vc, 1), mask), c0); + m0 =3D _mm_unpacklo_epi16(d0, _mm_setzero_si128()); + m1 =3D _mm_unpackhi_epi16(d0, _mm_setzero_si128()); + d1 =3D _mm_mullo_epi16(_mm_and_si128(_mm_srli_si128(vc, 2), mask), c0); + m2 =3D _mm_unpacklo_epi16(d1, _mm_setzero_si128()); + m3 =3D _mm_unpackhi_epi16(d1, _mm_setzero_si128()); + =09 + for (k=3D0; k<sizeof(f)/sizeof(*f); k++) { + vt =3D _mm_loadu_si128((__m128i *)(f[k] - 1 - w)); + vc =3D _mm_loadu_si128((__m128i *)(f[k] - 1 )); + vb =3D _mm_loadu_si128((__m128i *)(f[k] - 1 + w)); + f[k] +=3D 14; + =09 + /* even pixels */ + d0 =3D tf0(mask, l0, vt, vc, vb); + /* d =3D l - abs(r-d) */ + d0 =3D _mm_subs_epu16(l0, _mm_sub_epi16(_mm_max_epi16(r0, d0), _mm_min_= epi16(r0, d0))); + c0 =3D _mm_add_epi16(c0, d0); + /* d *=3D *f[k] */ + d0 =3D _mm_mullo_epi16(_mm_and_si128(_mm_srli_si128(vc, 1), mask), d0); + m0 =3D _mm_add_epi32(m0, _mm_unpacklo_epi16(d0, _mm_setzero_si128())); + m1 =3D _mm_add_epi32(m1, _mm_unpackhi_epi16(d0, _mm_setzero_si128())); + =09 + /* odd pixels */ + d1 =3D tf1(mask, l0, vt, vc, vb); + d1 =3D _mm_subs_epu16(l0, _mm_sub_epi16(_mm_max_epi16(r1, d1), _mm_min_= epi16(r1, d1))); + c1 =3D _mm_add_epi16(c1, d1); + d1 =3D _mm_mullo_epi16(_mm_and_si128(_mm_srli_si128(vc, 2), mask), d1); + m2 =3D _mm_add_epi32(m2, _mm_unpacklo_epi16(d1, _mm_setzero_si128())); + m3 =3D _mm_add_epi32(m3, _mm_unpackhi_epi16(d1, _mm_setzero_si128())); + } + =09 + /* extract results from c0, c1 and m0 to m3: + * 14 byte result =3D interleave_bytes((m0,m1) / c0, (m2,m3) / c1) */ +#ifdef OLD_ROUNDING + /* r =3D m*2/c */ + /* multiply each m with 2 */ + m0 =3D _mm_slli_epi32(m0, 1); + m1 =3D _mm_slli_epi32(m1, 1); + m2 =3D _mm_slli_epi32(m2, 1); + m3 =3D _mm_slli_epi32(m3, 1); +#endif + =09 + /* m0-m3 each contain 4 19-bit values ((8-bit)^2 * 8), so a single preci= sion + * float with 23-bit mantissa can hold these without precision loss */ + /* r =3D m/c */ + __m128i k0 =3D _mm_setzero_si128(); + __m128 f0, f1, f2, f3; + f0 =3D _mm_div_ps(_mm_cvtepi32_ps(m0), _mm_cvtepi32_ps(_mm_unpacklo_epi1= 6(c0, k0))); + f1 =3D _mm_div_ps(_mm_cvtepi32_ps(m1), _mm_cvtepi32_ps(_mm_unpackhi_epi1= 6(c0, k0))); + f2 =3D _mm_div_ps(_mm_cvtepi32_ps(m2), _mm_cvtepi32_ps(_mm_unpacklo_epi1= 6(c1, k0))); + f3 =3D _mm_div_ps(_mm_cvtepi32_ps(m3), _mm_cvtepi32_ps(_mm_unpackhi_epi1= 6(c1, k0))); + =09 +#ifdef OLD_ROUNDING + m0 =3D _mm_cvttps_epi32(f0); + m1 =3D _mm_cvttps_epi32(f1); + m2 =3D _mm_cvttps_epi32(f2); + m3 =3D _mm_cvttps_epi32(f3); +#else + m0 =3D _mm_cvtps_epi32(f0); + m1 =3D _mm_cvtps_epi32(f1); + m2 =3D _mm_cvtps_epi32(f2); + m3 =3D _mm_cvtps_epi32(f3); +#endif + =09 + r0 =3D _mm_packs_epi32(m0, m1); /* 7 words f0,f1 */ + r1 =3D _mm_packs_epi32(m2, m3); /* 7 words f2,f3 */ + =09 +#ifdef OLD_ROUNDING + /* (r+1)/2 */ + k0 =3D _mm_set1_epi16(1); + r0 =3D _mm_srli_epi16(_mm_add_epi16(r0, k0), 1); + r1 =3D _mm_srli_epi16(_mm_add_epi16(r1, k0), 1); +#endif + =09 + /* 7 words r0 interleaved with 7 words r1, all converted to bytes */ + r0 =3D _mm_packus_epi16(_mm_unpacklo_epi16(r0, r1), _mm_unpackhi_epi16(r= 0, r1)); + /* write 16, but the 2 bytes overlap will be overwritten by the next pas= s */ + _mm_storeu_si128((__m128i *)of, r0); + of +=3D 14; + } + _mm_empty(); +} +#endif + +void temporal_filter_planes_p (int idx, int w, int h, int t) +{ + uint32_t r, c, m; + int32_t d; + int x; + + uint8_t *f1 =3D frame1[idx]; + uint8_t *f2 =3D frame2[idx]; + uint8_t *f3 =3D frame3[idx]; + uint8_t *f4 =3D frame4[idx]; + uint8_t *f5 =3D frame5[idx]; + uint8_t *f6 =3D frame6[idx]; + uint8_t *f7 =3D frame7[idx]; + uint8_t *of =3D outframe[idx]; + + if (t =3D=3D 0) // shortcircuit filter if t =3D 0... + { + memcpy (of, f4, w * h); + return; + } + + for (x =3D 0; x < (w * h); x++) + { + r =3D *(f4-1-w); + r +=3D *(f4 -w)*2; + r +=3D *(f4+1-w); + r +=3D *(f4-1 )*2; + r +=3D *(f4 )*4; + r +=3D *(f4+1 )*2; + r +=3D *(f4-1+w); + r +=3D *(f4 +w)*2; + r +=3D *(f4+1+w); + r /=3D 16; + + m =3D *(f4)*(t+1)*2; + c =3D t+1; + + d =3D *(f3-1-w); + d +=3D *(f3 -w)*2; + d +=3D *(f3+1-w); + d +=3D *(f3-1 )*2; + d +=3D *(f3 )*4; + d +=3D *(f3+1 )*2; + d +=3D *(f3-1+w); + d +=3D *(f3 +w)*2; + d +=3D *(f3+1+w); + d /=3D 16; + + d =3D t - abs (r-d); + d =3D d<0? 0:d; + c +=3D d; + m +=3D *(f3)*d*2; + + d =3D *(f2-1-w); + d +=3D *(f2 -w)*2; + d +=3D *(f2+1-w); + d +=3D *(f2-1 )*2; + d +=3D *(f2 )*4; + d +=3D *(f2+1 )*2; + d +=3D *(f2-1+w); + d +=3D *(f2 +w)*2; + d +=3D *(f2+1+w); + d /=3D 16; + + d =3D t - abs (r-d); + d =3D d<0? 0:d; + c +=3D d; + m +=3D *(f2)*d*2; + + d =3D *(f1-1-w); + d +=3D *(f1 -w)*2; + d +=3D *(f1+1-w); + d +=3D *(f1-1 )*2; + d +=3D *(f1 )*4; + d +=3D *(f1+1 )*2; + d +=3D *(f1-1+w); + d +=3D *(f1 +w)*2; + d +=3D *(f1+1+w); + d /=3D 16; + + d =3D t - abs (r-d); + d =3D d<0? 0:d; + c +=3D d; + m +=3D *(f1)*d*2; + + d =3D *(f5-1-w); + d +=3D *(f5 -w)*2; + d +=3D *(f5+1-w); + d +=3D *(f5-1 )*2; + d +=3D *(f5 )*4; + d +=3D *(f5+1 )*2; + d +=3D *(f5-1+w); + d +=3D *(f5 +w)*2; + d +=3D *(f5+1+w); + d /=3D 16; + + d =3D t - abs (r-d); + d =3D d<0? 0:d; + c +=3D d; + m +=3D *(f5)*d*2; + + d =3D *(f6-1-w); + d +=3D *(f6 -w)*2; + d +=3D *(f6+1-w); + d +=3D *(f6-1 )*2; + d +=3D *(f6 )*4; + d +=3D *(f6+1 )*2; + d +=3D *(f6-1+w); + d +=3D *(f6 +w)*2; + d +=3D *(f6+1+w); + d /=3D 16; + + d =3D t - abs (r-d); + d =3D d<0? 0:d; + c +=3D d; + m +=3D *(f6)*d*2; + + d =3D *(f7-1-w); + d +=3D *(f7 -w)*2; + d +=3D *(f7+1-w); + d +=3D *(f7-1 )*2; + d +=3D *(f7 )*4; + d +=3D *(f7+1 )*2; + d +=3D *(f7-1+w); + d +=3D *(f7 +w)*2; + d +=3D *(f7+1+w); + d /=3D 16; + + d =3D t - abs (r-d); + d =3D d<0? 0:d; + c +=3D d; + m +=3D *(f7)*d*2; + + *(of) =3D ((m/c)+1)/2; + + f1++; + f2++; + f3++; + f4++; + f5++; + f6++; + f7++; + of++; + } +} + +#if defined(__SSE2__) +/* 4 to 5 times faster */ +void filter_plane_median_sse2(uint8_t *plane, int w, int h, int level) { + int i; + int avg; + int cnt; + uint8_t * p; + uint8_t * d; +=09 + if(level=3D=3D0) return; +=09 + p =3D plane; + d =3D scratchplane1; =20 -void filter_plane_median ( uint8_t * plane, int w, int h, int level) + // remove strong outliers from the image. An outlier is a pixel which lie= s outside + // of max-thres and min+thres of the surrounding pixels. This should not = cause blurring + // and it should leave an evenly spread noise-floor to the image. + for (i=3D0; i<=3D(w*h); i+=3D14) { + __m128i t, c, b, min, max, minmin, maxmax; + =09 + t =3D _mm_loadu_si128((__m128i *)&p[i-1-w]); + c =3D _mm_loadu_si128((__m128i *)&p[i-1 ]); + b =3D _mm_loadu_si128((__m128i *)&p[i-1+w]); + min =3D _mm_min_epu8(t, b); /* k: (0,k), (2,k) */ + max =3D _mm_max_epu8(t, b); + minmin =3D _mm_min_epu8(min, c); /* k: (0,k), (1,k), (2,k) */ + maxmax =3D _mm_max_epu8(max, c); + =09 + /* k: (0,k), (1,k), (2,k). (0,k+2), (1,k+2), (2,k+2) */ + minmin =3D _mm_min_epu8(minmin, _mm_srli_si128(minmin, 2)); + maxmax =3D _mm_max_epu8(maxmax, _mm_srli_si128(maxmax, 2)); + /* k: (0,k), (1,k), (2,k). (0,k+2), (1,k+2), (2,k+2), (0,k+1), (2,k+1) */ + min =3D _mm_min_epu8(minmin, _mm_srli_si128(min, 1)); + max =3D _mm_max_epu8(maxmax, _mm_srli_si128(max, 1)); + =09 + /* limit c to range [min,max] */ + c =3D _mm_max_epu8(min, _mm_min_epu8(max, _mm_srli_si128(c, 1))); + /* write 14 valid pixels, the 2 remaining bytes are overwritten subseque= ntly + * or lie outside the frame area */ + _mm_storeu_si128((__m128i *)&d[i], c); + } +=09 + // in the second stage we try to average similar spatial pixels, only. Th= is, like + // a median, should also not reduce sharpness but flatten the noisefloor.= This + // part is quite similar to what 2dclean/yuvmedianfilter do. But because = of the + // different weights given to the pixels it is less aggressive... + + p =3D scratchplane1; + d =3D scratchplane2; +=09 + // this filter needs values outside of the imageplane, so we just copy th= e first line=20 + // and the last line into the out-of-range area... +=09 + memcpy ( p-w , p, w ); + memcpy ( p-w*2, p, w ); +=09 + memcpy ( p+(w*h) , p+(w*h)-w, w ); + memcpy ( p+(w*h)+w, p+(w*h)-w, w ); +=09 + __m128i lvl =3D _mm_set1_epi16(level); +=09 +#ifndef OLD_ROUNDING + _MM_SET_ROUNDING_MODE(_MM_ROUND_NEAREST); +#endif +=09 + for (i=3D0; i<=3D(w*h); i+=3D4) + { + uint64_t k0, k1, k2, k3, k6; + __m128i c0, c1, v[4], t[4], e[4], a[4]; + =09 + /* p points to pixel a. There are 3 stages, each processing 8 surrounding + * pixels, resulting in the complete neighbourhood of 24 pixels. + * + * 0 1 2 3 4 5 6 7 + * k0: x x x x x x x x + * k1: x x x x x x x x + * k6: x x a b c d x x + * k2: x x x x x x x x + * k3: x x x x x x x x + * + * a,b and c,d each share two 2x4-blocks in their surrounding area, which + * are processed by stages 1 and 2. These blocks are referred to as c0 a= nd + * c1 by the code below. The remaining surrounding pixels are processed = for + * each pixel out of a,b,c,d individually in stage 3. + * Stage 4 assembles the results, adds the weighted averages and weights + * together for each pixel, computes the median and stores it in the + * scratch plane. + * Coordinates (y,x) originate from the top left of the above diagram. */ + =09 + k0 =3D *(uint64_t *)(p-w*2-2); + k1 =3D *(uint64_t *)(p-w*1-2); + k6 =3D *(uint64_t *)(p -2); + k2 =3D *(uint64_t *)(p+w*1-2); + k3 =3D *(uint64_t *)(p+w*2-2); + =09 + v[0] =3D _mm_set1_epi16((k6 >> 16) & 0xff); /* pixel a */ + v[1] =3D _mm_set1_epi16((k6 >> 24) & 0xff); /* pixel b */ + v[2] =3D _mm_set1_epi16((k6 >> 32) & 0xff); /* pixel c */ + v[3] =3D _mm_set1_epi16((k6 >> 40) & 0xff); /* pixel d */ + =09 + // stage 1: c0 for a,b: (0,1) -> (1,4), c1 for c,d: (0,3) -> (1,6) + c0 =3D _mm_set_epi32((k0 >> 8) & 0xff00ff, (k0 >> 16) & 0xff00ff, + (k1 >> 8) & 0xff00ff, (k1 >> 16) & 0xff00ff); + c1 =3D _mm_set_epi32((k0 >> 24) & 0xff00ff, (k0 >> 32) & 0xff00ff, + (k1 >> 24) & 0xff00ff, (k1 >> 32) & 0xff00ff); + t[0] =3D _mm_sub_epi16(_mm_max_epu8(c0, v[0]), _mm_min_epu8(c0, v[0])); + t[1] =3D _mm_sub_epi16(_mm_max_epu8(c0, v[1]), _mm_min_epu8(c0, v[1])); + t[2] =3D _mm_sub_epi16(_mm_max_epu8(c1, v[2]), _mm_min_epu8(c1, v[2])); + t[3] =3D _mm_sub_epi16(_mm_max_epu8(c1, v[3]), _mm_min_epu8(c1, v[3])); + e[0] =3D _mm_subs_epu16(lvl, t[0]); + e[1] =3D _mm_subs_epu16(lvl, t[1]); + e[2] =3D _mm_subs_epu16(lvl, t[2]); + e[3] =3D _mm_subs_epu16(lvl, t[3]); + a[0] =3D _mm_madd_epi16(e[0], c0); + a[1] =3D _mm_madd_epi16(e[1], c0); + a[2] =3D _mm_madd_epi16(e[2], c1); + a[3] =3D _mm_madd_epi16(e[3], c1); + =09 + // stage 2: c0 for a,b: (3,1) -> (4,4), c1 for c,d: (3,3) -> (4,6) + c0 =3D _mm_set_epi32((k2 >> 8) & 0xff00ff, (k2 >> 16) & 0xff00ff, + (k3 >> 8) & 0xff00ff, (k3 >> 16) & 0xff00ff); + c1 =3D _mm_set_epi32((k2 >> 24) & 0xff00ff, (k2 >> 32) & 0xff00ff, + (k3 >> 24) & 0xff00ff, (k3 >> 32) & 0xff00ff); + t[0] =3D _mm_sub_epi16(_mm_max_epu8(c0, v[0]), _mm_min_epu8(c0, v[0])); + t[1] =3D _mm_sub_epi16(_mm_max_epu8(c0, v[1]), _mm_min_epu8(c0, v[1])); + t[2] =3D _mm_sub_epi16(_mm_max_epu8(c1, v[2]), _mm_min_epu8(c1, v[2])); + t[3] =3D _mm_sub_epi16(_mm_max_epu8(c1, v[3]), _mm_min_epu8(c1, v[3])); + t[0] =3D _mm_subs_epu16(lvl, t[0]); + t[1] =3D _mm_subs_epu16(lvl, t[1]); + t[2] =3D _mm_subs_epu16(lvl, t[2]); + t[3] =3D _mm_subs_epu16(lvl, t[3]); + e[0] =3D _mm_add_epi16(t[0], e[0]); + e[1] =3D _mm_add_epi16(t[1], e[1]); + e[2] =3D _mm_add_epi16(t[2], e[2]); + e[3] =3D _mm_add_epi16(t[3], e[3]); + a[0] =3D _mm_add_epi32(_mm_madd_epi16(t[0], c0), a[0]); + a[1] =3D _mm_add_epi32(_mm_madd_epi16(t[1], c0), a[1]); + a[2] =3D _mm_add_epi32(_mm_madd_epi16(t[2], c1), a[2]); + a[3] =3D _mm_add_epi32(_mm_madd_epi16(t[3], c1), a[3]); + =09 + // stage 3: + // pixel a: (0,0) -> (4,0), (2,1), (2,3), (2,4) + c0 =3D _mm_set_epi32(((k0 & 0xff) << 16) | (k1 & 0xff), + ((k2 & 0xff) << 16) | (k3 & 0xff), + (k6 >> 8) & 0xff00ff, + (k6 & 0xff) | ((k6 >> 16) & 0xff0000)); + t[0] =3D _mm_sub_epi16(_mm_max_epu8(c0, v[0]), _mm_min_epu8(c0, v[0])); + t[0] =3D _mm_subs_epu16(lvl, t[0]); + e[0] =3D _mm_add_epi16(t[0], e[0]); + a[0] =3D _mm_add_epi32(_mm_madd_epi16(t[0], c0), a[0]); + =09 + // pixel c: (0,2) -> (4,2), (2,3), (2,5), (2,6) + c1 =3D _mm_set_epi32((k0 & 0xff0000) | ((k1 >> 16) & 0xff), + (k2 & 0xff0000) | ((k3 >> 16) & 0xff), + (k6 >> 24) & 0xff00ff, + ((k6 >> 16) & 0xff) | ((k6 >> 32) & 0xff0000)); + t[2] =3D _mm_sub_epi16(_mm_max_epu8(c1, v[2]), _mm_min_epu8(c1, v[2])); + t[2] =3D _mm_subs_epu16(lvl, t[2]); + e[2] =3D _mm_add_epi16(t[2], e[2]); + a[2] =3D _mm_add_epi32(_mm_madd_epi16(t[2], c1), a[2]); + =09 + // pixel b: (0,5) -> (4,5), (2,1), (2,2), (2,4), (2,5) + c0 =3D _mm_set_epi32(((k0 >> 24) & 0xff0000) | ((k1 >> 40) & 0xff), + ((k2 >> 24) & 0xff0000) | ((k3 >> 40) & 0xff), + (k6 >> 16) & 0xff00ff, + ((k6 >> 8) & 0xff) | ((k6 >> 24) & 0xff0000)); + t[1] =3D _mm_sub_epi16(_mm_max_epu8(c0, v[1]), _mm_min_epu8(c0, v[1])); + t[1] =3D _mm_subs_epu16(lvl, t[1]); + e[1] =3D _mm_add_epi16(t[1], e[1]); + a[1] =3D _mm_add_epi32(_mm_madd_epi16(t[1], c0), a[1]); + =09 + // pixel d: (0,7) -> (4,7), (2,3), (2,4), (2,6), (2,7) + c1 =3D _mm_set_epi32(((k0 >> 40) & 0xff0000) | (k1 >> 56), + ((k2 >> 40) & 0xff0000) | (k3 >> 56), + (k6 >> 32) & 0xff00ff, + ((k6 >> 24) & 0xff) | ((k6 >> 40) & 0xff0000)); + t[3] =3D _mm_sub_epi16(_mm_max_epu8(c1, v[3]), _mm_min_epu8(c1, v[3])); + t[3] =3D _mm_subs_epu16(lvl, t[3]); + e[3] =3D _mm_add_epi16(t[3], e[3]); + a[3] =3D _mm_add_epi32(_mm_madd_epi16(t[3], c1), a[3]); + =09 + /* add them all together (a loop j=3D0 to 4 slows things down with gcc 4= .4.4) */ + uint32_t tmp; + =09 +#if defined(__SSE3__) + int j; + __m128 f0, f1, f2, flvl; + __m128i zero, vv; + flvl =3D _mm_set1_ps((float)level); + zero =3D _mm_setzero_si128(); + vv =3D _mm_set_epi32((k6 >> 40) & 0xff, (k6 >> 32) & 0xff, (k6 >> 24) & = 0xff, (k6 >> 16) & 0xff); + =09 + f0 =3D _mm_hadd_ps(_mm_cvtepi32_ps(a[0]), _mm_cvtepi32_ps(a[1])); + f1 =3D _mm_hadd_ps(_mm_cvtepi32_ps(a[2]), _mm_cvtepi32_ps(a[3])); + f0 =3D _mm_hadd_ps(f0, f1); + f0 =3D _mm_add_ps(f0, _mm_mul_ps(flvl, _mm_cvtepi32_ps(vv))); + =09 +# ifdef OLD_ROUNDING + /* avg *=3D 2 */ + f0 =3D _mm_mul_ps(f0, _mm_set1_ps(2.0f)); +# endif + =09 + for (j=3D0; j<4; j++) + e[j] =3D _mm_unpacklo_epi16(_mm_add_epi16(e[j], _mm_srli_si128(e[j], 8)= ), zero); + f1 =3D _mm_hadd_ps(_mm_cvtepi32_ps(e[0]), _mm_cvtepi32_ps(e[1])); + f2 =3D _mm_hadd_ps(_mm_cvtepi32_ps(e[2]), _mm_cvtepi32_ps(e[3])); + f1 =3D _mm_hadd_ps(f1, f2); + f1 =3D _mm_add_ps(f1, flvl); + =09 + f0 =3D _mm_div_ps(f0, f1); +# ifdef OLD_ROUNDING + /* r =3D (r+1) / 2 */ + vv =3D _mm_cvttps_epi32(f0); + vv =3D _mm_srli_epi32(_mm_add_epi32(vv, _mm_set1_epi32(1)), 1); +# else + vv =3D _mm_cvtps_epi32(f0); +# endif + vv =3D _mm_packus_epi16(_mm_packs_epi32(vv, zero), zero); + tmp =3D _mm_cvtsi128_si32(vv); + =09 +#else + =09 + // p[0] + e[0] =3D _mm_add_epi16(e[0], _mm_srli_si128(e[0], 8)); /* 8 words */ + e[0] =3D _mm_add_epi16(e[0], _mm_srli_si128(e[0], 4)); + e[0] =3D _mm_add_epi16(e[0], _mm_srli_si128(e[0], 2)); + cnt =3D level + (_mm_cvtsi128_si32(e[0]) & 0xffff); + a[0] =3D _mm_add_epi32(a[0], _mm_srli_si128(a[0], 8)); /* 4 dwords */ + a[0] =3D _mm_add_epi32(a[0], _mm_srli_si128(a[0], 4)); + avg =3D p[0] * level * 2 + (_mm_cvtsi128_si32(a[0]) << 1); + tmp =3D ((avg/cnt) + 1) / 2; + =09 + // p[1] + e[1] =3D _mm_add_epi16(e[1], _mm_srli_si128(e[1], 8)); /* 8 words */ + e[1] =3D _mm_add_epi16(e[1], _mm_srli_si128(e[1], 4)); + e[1] =3D _mm_add_epi16(e[1], _mm_srli_si128(e[1], 2)); + cnt =3D level + (_mm_cvtsi128_si32(e[1]) & 0xffff); + a[1] =3D _mm_add_epi32(a[1], _mm_srli_si128(a[1], 8)); /* 4 dwords */ + a[1] =3D _mm_add_epi32(a[1], _mm_srli_si128(a[1], 4)); + avg =3D p[1] * level * 2 + (_mm_cvtsi128_si32(a[1]) << 1); + tmp |=3D (((avg/cnt) + 1) / 2) << 8; + =09 + // p[2] + e[2] =3D _mm_add_epi16(e[2], _mm_srli_si128(e[2], 8)); /* 8 words */ + e[2] =3D _mm_add_epi16(e[2], _mm_srli_si128(e[2], 4)); + e[2] =3D _mm_add_epi16(e[2], _mm_srli_si128(e[2], 2)); + cnt =3D level + (_mm_cvtsi128_si32(e[2]) & 0xffff); + a[2] =3D _mm_add_epi32(a[2], _mm_srli_si128(a[2], 8)); /* 4 dwords */ + a[2] =3D _mm_add_epi32(a[2], _mm_srli_si128(a[2], 4)); + avg =3D p[2] * level * 2 + (_mm_cvtsi128_si32(a[2]) << 1); + tmp |=3D (((avg/cnt) + 1) / 2) << 16; + =09 + // p[3] + e[3] =3D _mm_add_epi16(e[3], _mm_srli_si128(e[3], 8)); /* 8 words */ + e[3] =3D _mm_add_epi16(e[3], _mm_srli_si128(e[3], 4)); + e[3] =3D _mm_add_epi16(e[3], _mm_srli_si128(e[3], 2)); + cnt =3D level + (_mm_cvtsi128_si32(e[3]) & 0xffff); + a[3] =3D _mm_add_epi32(a[3], _mm_srli_si128(a[3], 8)); /* 4 dwords */ + a[3] =3D _mm_add_epi32(a[3], _mm_srli_si128(a[3], 4)); + avg =3D p[3] * level * 2 + (_mm_cvtsi128_si32(a[3]) << 1); + tmp |=3D (((avg/cnt) + 1) / 2) << 24; +#endif + =09 + *(uint32_t *)d =3D tmp; + =09 + d +=3D 4; + p +=3D 4; + } + _mm_empty(); +=09 + memcpy(plane,scratchplane2,w*h); +} +#endif + +void filter_plane_median_p ( uint8_t * plane, int w, int h, int level) { int i; int min; @@ -825,6 +1323,27 @@ * Main Loop * ***********************************************************/ =20 +static void init_accel() { + filter_plane_median =3D filter_plane_median_p; + temporal_filter_planes =3D temporal_filter_planes_p; +=09 +#if defined(__SSE2__) + int d =3D 0; + __asm__ volatile("cpuid" : "=3Dd"(d) : "a"(1) : "ebx", "ecx"); + if ((d & (1 << 26))) { + mjpeg_info("SETTING SSE2 for standard Temporal-Noise-Filter"); + temporal_filter_planes =3D temporal_filter_planes_sse2; + =09 + __asm__ volatile("cpuid" : "=3Dd"(d) : "a"(0x80000001) : "ebx", "ecx"); + if ((d & (1 << 29))) { + /* x86_64 processor */ + mjpeg_info("SETTING SSE2 for Median-Filter"); + filter_plane_median =3D filter_plane_median_sse2; + } + } +#endif +} + int main (int argc, char *argv[]) { @@ -1062,6 +1581,8 @@ /* initialize motion_library */ init_motion_search (); =20 + init_accel(); + /* read every frame until the end of the input stream and process it */ while (Y4M_OK =3D=3D (err =3D y4m_read_frame (fd_in, &istreaminfo, --Multipart=_Fri__1_Oct_2010_15_47_46_+0200_3M22tzehqubz742L Content-Type: text/plain; charset="us-ascii" MIME-Version: 1.0 Content-Transfer-Encoding: 7bit Content-Disposition: inline ------------------------------------------------------------------------------ Start uncovering the many advantages of virtual appliances and start using them to simplify application deployment and accelerate your shift to cloud computing. http://p.sf.net/sfu/novell-sfdev2dev --Multipart=_Fri__1_Oct_2010_15_47_46_+0200_3M22tzehqubz742L Content-Type: text/plain; charset="us-ascii" MIME-Version: 1.0 Content-Transfer-Encoding: 7bit Content-Disposition: inline _______________________________________________ Mjpeg-developer mailing list [email protected] https://lists.sourceforge.net/lists/listinfo/mjpeg-developer --Multipart=_Fri__1_Oct_2010_15_47_46_+0200_3M22tzehqubz742L--