Add SIMD implementation of Shift() function

[ihhub]

We have an implementation of Shift() function in src/image_function.h and src/image_function.cpp files. We want to add SIMD implementation of the same function. The code must be located in src/image_function_simd.h and src/image_function_simd.cpp files. Please to any of functions inside these files.

[ihhub]

@0x72D0 this should be interesting and easy to do ;)

[0x72D0]

I will give it a try

[ihhub]

This issue has higher priority than Hough transform and it should be easier to implement so you could do this first :)

[ihhub]

By the way, unit tests for Shift() function are not trully correct. They need some modification. From another side I could describe a way how to teat this function properly.

[0x72D0]

are you sure that it's suppose to be more easy? the solution I found to replace the main for loop for the shift function look like that:

void Shift(uint32_t rowSizeIn, uint32_t rowSizeOut, uint8_t * outY, const uint8_t * outYEnd, const uint8_t * inY,
               uint32_t coeff0, uint32_t coeff1, uint32_t coeff2, uint32_t coeff3, uint32_t emptyLeftArea,
               uint32_t emptyRightArea, uint32_t simdWidth, uint32_t totalSimdWidth, uint32_t nonSimdWidth)
    {
        const simd zero = _mm256_setzero_si256();

        const simd vcoeff0 = _mm256_set1_epi32(coeff0);
        const simd vcoeff1 = _mm256_set1_epi32(coeff1);
        const simd vcoeff2 = _mm256_set1_epi32(coeff2);
        const simd vcoeff3 = _mm256_set1_epi32(coeff3)

        for( ; outY != outYEnd; outY+=rowSizeOut, inY+=rowSizeIn ) {
            memset( outY, 0, emptyLeftArea );
            uint8_t * inX = inY + emptyLeftArea

            const simd * src    = reinterpret_cast <const simd*> (inX);
            simd       * dst    = reinterpret_cast <const simd*> (outY);
            const simd * dstEnd = dst + simdWith;

            for( ; dst != dstEnd; ++dst ) {

                // take the vector that we multiply by coeff0
                const simd * src0    = reinterpret_cast <const simd*> (inX);
                const simd data_coeff0  = _mm256_loadu_si256( src0 );
                const simd data_lo_coeff0  = _mm256_unpacklo_epi8( data_coeff0, zero );
                const simd data1_coeff0  = _mm256_unpacklo_epi16( data_lo_coeff0, zero );
                simd data_1 = _mm256_mullo_epi32(data1_coeff0, vcoeff0);

                const simd data2_coeff0  = _mm256_unpackhi_epi16( data_lo_coeff0, zero );
                simd data_2 = _mm256_mullo_epi32(data2_coeff0, vcoeff0);

                const simd data_hi_coeff0  = _mm256_unpackhi_epi8( data_coeff0, zero );
                const simd data3_coeff0  = _mm256_unpacklo_epi16( data_hi_coeff0, zero );
                simd data_3 = _mm256_mullo_epi32(data3_coeff0, vcoeff0);

                const simd data4_coeff0  = _mm256_unpackhi_epi16( data_hi_coeff0, zero );
                simd data_4 = _mm256_mullo_epi32(data4_coeff0, vcoeff0);

                // take the vector that we multiply by coeff1
                const simd * src1    = reinterpret_cast <const simd*> (inX+rowSizeIn);
                const simd data_coeff1  = _mm256_loadu_si256( src1 );
                const simd data_lo_coeff1  = _mm256_unpacklo_epi8( data_coeff1, zero );
                const simd data1_coeff1  = _mm256_unpacklo_epi16( data_lo_coeff1, zero );
                data_1 = _mm256_add_epi32(data_1, _mm256_mullo_epi32(data1_coeff1, vcoeff1));

                const simd data2_coeff1  = _mm256_unpackhi_epi16( data_lo_coeff1, zero );
                data_2 = _mm256_add_epi32(data_2, _mm256_mullo_epi32(data2_coeff1, vcoeff1));

                const simd data_hi_coeff1  = _mm256_unpackhi_epi8( data_coeff1, zero );
                const simd data3_coeff1  = _mm256_unpacklo_epi16( data_hi_coeff1, zero );
                data_3 = _mm256_add_epi32(data_3, _mm256_mullo_epi32(data3_coeff1, vcoeff1));

                const simd data4_coeff1  = _mm256_unpackhi_epi16( data_hi_coeff1, zero );
                data_4 = _mm256_add_epi32(data_4, _mm256_mullo_epi32(data4_coeff1, vcoeff1));

                // take the vector that we multiply by coeff3
                const simd * src2    = reinterpret_cast <const simd*> (++inX);
                const simd data_coeff3  = _mm256_loadu_si256( src2 );
                const simd data_lo_coeff3  = _mm256_unpacklo_epi8( data_coeff3, zero );
                const simd data1_coeff3  = _mm256_unpacklo_epi16( data_lo_coeff3, zero );
                data_1 = _mm256_add_epi32(data_1, _mm256_mullo_epi32(data1_coeff3, vcoeff3));

                const simd data2_coeff3  = _mm256_unpackhi_epi16( data_lo_coeff3, zero );
                data_2 = _mm256_add_epi32(data_2, _mm256_mullo_epi32(data2_coeff3, vcoeff3));

                const simd data_hi_coeff3  = _mm256_unpackhi_epi8( data_coeff3, zero );
                const simd data3_coeff3  = _mm256_unpacklo_epi16( data_hi_coeff3, zero );
                data_3 = _mm256_add_epi32(data_3, _mm256_mullo_epi32(data3_coeff3, vcoeff3));

                const simd data4_coeff3  = _mm256_unpackhi_epi16( data_hi_coeff3, zero );
                data_4 = _mm256_add_epi32(data_4, _mm256_mullo_epi32(data4_coeff3, vcoeff3));

                // take the vector that we multiply by coeff2
                const simd * src3    = reinterpret_cast <const simd*> (inX+rowSizeIn);
                const simd data_coeff2  = _mm256_loadu_si256( src3 );
                const simd data_lo_coeff2  = _mm256_unpacklo_epi8( data_coeff2, zero );
                const simd data1_coeff2  = _mm256_unpacklo_epi16( data_lo_coeff2, zero );
                data_1 = _mm256_add_epi32(data_1, _mm256_mullo_epi32(data1_coeff2, vcoeff2));

                const simd data2_coeff2  = _mm256_unpackhi_epi16( data_lo_coeff2, zero );
                data_2 = _mm256_add_epi32(data_2, _mm256_mullo_epi32(data2_coeff2, vcoeff2));

                const simd data_hi_coeff2  = _mm256_unpackhi_epi8( data_coeff2, zero );
                const simd data3_coeff2  = _mm256_unpacklo_epi16( data_hi_coeff2, zero );
                data_3 = _mm256_add_epi32(data_3, _mm256_mullo_epi32(data3_coeff2, vcoeff2));

                const simd data4_coeff2  = _mm256_unpackhi_epi16( data_hi_coeff2, zero );
                data_4 = _mm256_add_epi32(data_4, _mm256_mullo_epi32(data4_coeff2, vcoeff2));

                data_1 = _mm256_add_epi32(data_1, 8192);
                data_1 = _mm256_srai_epi32(data_1, 14);

                data_2 = _mm256_add_epi32(data_2, 8192);
                data_2 = _mm256_srai_epi32(data_2, 14);

                data_3 = _mm256_add_epi32(data_3, 8192);
                data_3 = _mm256_srai_epi32(data_3, 14);

                data_4 = _mm256_add_epi32(data_4, 8192);
                data_4 = _mm256_srai_epi32(data_4, 14);

                simd data_1_2 =  _mm256_packus_epi32 (data_1, data_2);
                simd data_3_4 =  _mm256_packus_epi32 (data_3, data_4);

                _mm256_storeu_si256(dst, _mm256_packus_epi32 (data_1_2, data_3_4));

            }

            if( nonSimdWidth > 0 )
            {
                const uint8_t * inX     = inY + totalSimdWidth;
                uint8_t       * outX    = outY + totalSimdWidth;
                const uint8_t * outXEnd = outX + nonSimdWidth;

                for( ; outX != outXEnd; ++outX ) {
                    uint32_t data  = (*(inX))  * coeff0;
                    data += (*(inX+rowSizeIn)) * coeff3;
                    data += (*(++inX))         * coeff1; // here we increment inX
                    data += (*(inX+rowSizeIn)) * coeff2;

                    *outX = static_cast<uint8_t>((data + 8192) >> 14); // 8192 is 0.5 after calculations
                }
            }

            memset( outY, 0, emptyRightArea );

        }

so what was youre idea of the shift function? maybe I over complicated the problem?

[0x72D0]

I didn't test it yet so I don't know if this algorithm is fine

[ihhub]

@0x72D0 you might not overcomplicate things :) From another side I didn't mention (which is my bad) that CPU code is with optimisation avoiding float usage. Original logic behind is:

value = pixel1 * coeff1 + pixel2 * coeff2 + pixel3 * coeff3 + pixel3 * coeff3 + 0.5f; // 0.5f is needed for rounding value to nearest

So you could try to implement this first rather going to greater optimisation. I will take deeper look into your code in a few hours.

[ihhub]

With your implementation we read input values 4 times: 2 times for first row and 2 times for second row (when we move by one row below). I'm proposing to use a different approach when we read values only once. I'll try to explain as much as possible (for SSE case).

Step 1: We need 4 simd variables for coefficients which are done by you (int32 or float it doesn't matter now):

1. [coeff1; coeff1; coeff1; coeff1]
2. [coeff2; coeff2; coeff2; coeff2]
3. [coeff3; coeff3; coeff3; coeff3]
4. [coeff4; coeff4; coeff4; coeff4]

These variables are outside any loop.

Step 2 Read data from first row and unpack it into 32bit format (float or int32):

• Initial data is [pixel1; pixel2; pixel3; pixel4; pixel5; pixel6; pixel7; pixel8; pixel9; pixel10; pixel11; pixel12; pixel13; pixel14; pixel15; pixel16]
• We unpack the data into 16 bit format [pixel1; pixel3; pixel5; pixel7; pixel9; pixel11; pixel13; pixel15] and [pixel2; pixel4; pixel6; pixel8; pixel10; pixel12; pixel14; pixel16] (2 variables)
• We unpack again and convert into float or int32: [pixel1; pixel5; pixel9; pixel13] and [pixel3; pixel7; pixel11; pixel15] ; [pixel2; pixel6; pixel10; pixel14] and [pixel4; pixel8; pixel12; pixel16] (4 variables)

Step 3 Find result for uneven numbers (pixel positions aka first pixel, third pixel and etc.). First two variables we multiple by coeff1, second two multiple by coeff2. Then we summarize first and third variable and second and forth variable. This is pretty straightforward.

Step 4 Find result for even numbers. This is tricky but solvable. We need to shift [pixel1; pixel5; pixel9; pixel13] by 8 bit having as a result simd [pixel5; pixel9; pixel13; 0]. Okay, then next we multiple first and second variables by coeff2 and third and forth - by coeff1. Summaries them as on previous step.

Everything looks good except the fact that our calculations for pixel15 are incorrect (by right we should have variable as [pixel5; pixel9; pixel13; pixel17]). Unfortunately it's impossible to do such as we can read only 16 bytes.. In such case we could just read every 15th byte instead of 16.

Step 5 Do the same for second row with coeff3 and coeff4. Step 6. Summarize results from 2 rows and write it into output array. Step 7. Swap already unpacked variables (we can use reference or pointer) from second row to first row. Repeat steps 3-7.

Of course it's just a common explanation to give you an idea :)

[ihhub]

As of now ignore step 7. Let's make working version first even with double data reading.

[0x72D0]

seeing the calcul: value = pixel1 coeff1 + pixel2 coeff2 + pixel3 coeff3 + pixel3 coeff3 + 0.5f; I think this function gonna meet the same fate than the hough transform...

[ihhub]

Yes and no. The implementation of Shift function uses integer values instead of float. I move this issue into Wishlist as of now as it's not urgent :)