pyramid.hpp

Go to the documentation of this file.

#pragma once

#include <armadillo>

namespace auxiliary
{
#ifdef USE_CXX11
    enum border_type : arma::uword
#else
    enum border_type
#endif
    {
        constant,   
        reflect,    
        replicate,  
        warp,       
        reflect101, 
        transparent 
    };

    inline arma::uword borderInterpolate(int p, int n, border_type type = reflect101)
    {
        //return (p < 0) ? -p : ((p < n) ? p : (n + n - p - 1));
        if ((unsigned)p < (unsigned)n) return p;

        switch (type) {
        case reflect:
        case reflect101:
            {
                int delta = (type == reflect101);

                if (n == 1) p = 0;

                do {
                    if (p < 0)
                        p = -p - 1 + delta;
                    else
                        p = n - 1 - (p - n) - delta;
                } while ((unsigned)p >= (unsigned)n);
            }
            break;
        case replicate:
            p = (p < 0) ? 0 : n - 1;
            break;
        case warp:
            if (p < 0)
                p -= ((p - n + 1) / n) * n;
            if (p > n)
                p %= n;
            break;
        case constant:
            p = -1;
                break;
        default:
            break;
        };

        return p;
    }

#define castOp(x) ((x + 128) >> 8)

    template <typename T1, typename T2>
    void pyrDown(const T1& in, T2& out)
    {
        const uword KERNEL_SIZE = 5;

        //uword width = std::min((src.n_cols - SZ / 2 - 1) / 2;
        
        circular_buffer<arma::ivec> cols(KERNEL_SIZE);

#ifdef __VXWORKS__
        ivec dummy(out.n_rows); dummy.zeros();
#endif
        for (arma::uword i = 0 ; i < KERNEL_SIZE ; i++)
#ifdef __VXWORKS__
            cols.push_back(dummy);
#else
            cols.push_back(zeros<ivec>(out.n_rows));
#endif

        int sx0 = -(int)KERNEL_SIZE / 2, sx = sx0;

        arma::umat tab(KERNEL_SIZE + 2, 2);
        uword* lptr = tab.colptr(0),
             * rptr = tab.colptr(1);
        for (uword y = 0 ; y <= KERNEL_SIZE + 1 ; y++) {
            lptr[y] = borderInterpolate((int)y + sx0, (int)in.n_rows);
            rptr[y] = borderInterpolate((int)(y + (out.n_rows - 1) * 2) + sx0, (int)in.n_rows);
        }

        // gaussian convolution with 
        for (arma::uword x = 0 ; x < out.n_cols ; x++) {
            typename T2::elem_type* dst = out.colptr(x);

            // vertical convolution and decimation
            for ( ; sx <= (int)x * 2 + 2 ; sx++) {
                ivec& col = cols.next();
                int* colptr = col.memptr();

                // interpolate border
                const typename T2::elem_type* src = in.colptr(borderInterpolate(sx, (int)in.n_cols));

                colptr[0] = src[lptr[2]] * 6 + (src[lptr[1]] + src[lptr[3]]) * 4 + (src[lptr[0]] + src[lptr[4]]);

                for (arma::uword y = 1 ; y < out.n_rows - 1; y++)
                //concurrency::parallel_for(uword(1), out.n_rows - 1, [&](uword y) {
                    colptr[y] = src[y * 2] * 6 + 
                             (src[y * 2 - 1] + src[y * 2 + 1]) * 4 + 
                             (src[y * 2 - 2] + src[y * 2 + 2]);
                //});

                colptr[out.n_rows - 1] = src[rptr[2]] * 6 + 
                                      (src[rptr[1]] + src[rptr[3]]) * 4 + 
                                      (src[rptr[0]] + src[rptr[4]]);
            }

            const int* col0 = cols[0].memptr();
            const int* col1 = cols[1].memptr();
            const int* col2 = cols[2].memptr();
            const int* col3 = cols[3].memptr();
            const int* col4 = cols[4].memptr();

            // horizontal convolution and decimation
#if ENABLE_SSE2
            //__m128i d = _mm_set1_epi16(128);
            //uword y = 0;
            //for ( ; y <= out.n_rows - 16 ; y += 16) {
            //  __m128i c0, c1, c2, c3, c4, t0, t1;
            //  c0 = _mm_packs_epi32(_mm_load_si128((const __m128i*)(col0 + y)),
            //                       _mm_load_si128((const __m128i*)(col0 + y + 4)));
            //  c1 = _mm_packs_epi32(_mm_load_si128((const __m128i*)(col1 + y)),
            //                       _mm_load_si128((const __m128i*)(col1 + y + 4)));
            //  c2 = _mm_packs_epi32(_mm_load_si128((const __m128i*)(col2 + y)),
            //                       _mm_load_si128((const __m128i*)(col2 + y + 4)));
            //  c3 = _mm_packs_epi32(_mm_load_si128((const __m128i*)(col3 + y)),
            //                       _mm_load_si128((const __m128i*)(col3 + y + 4)));
            //  c4 = _mm_packs_epi32(_mm_load_si128((const __m128i*)(col4 + y)),
            //                       _mm_load_si128((const __m128i*)(col4 + y + 4)));

            //  c0 = _mm_add_epi16(r0, r4);
            //  c1 = _mm_add_epi16(_mm_add_epi16(c1, c3), c2);
            //}
#else
            for (arma::uword y = 0 ; y < out.n_rows ; y++)
            //concurrency::parallel_for(uword(0), out.n_rows, [&](uword y) {
                dst[y] = (typename T2::elem_type)castOp(col2[y] * 6 + (col1[y] + col3[y]) * 4 + col0[y] + col4[y]);
            //});
#endif
        }
    }
}