2021年4月20日火曜日

OpenCV warpPerspective を魔改造


ラスター変換で、図面の画像を縮小変換した場合に、Cubic や Lanczos4 といった周りの画素値から平均を求める手法だと、線が薄くなって消えてしまいます。縮小変換しても線を消さないようにするために、周りの画素値の最小値を取ってやれば良いと思い魔改造する事に・・・。

魔改造を施したのは OpenCV 4.5.2 のバージョンです。

/modules/imgproc/src/imgwarp.cpp
を改変していきます。尚、OpenCL や OpenVX など、コンパイル条件で多様なルーチンを選択できますが、最小限のコードのみを対象としました。
warpPerspectiveでは cv::INTER_AREA というフラグを指定しても cv::INTER_LINEAR というフラグに変換されて処理されていたので、cv::INTER_AREA が指定された場合に、線画用の縮小処理が実行されるという方針で実装しました

手っ取り早くBilinear のコードをベースに改造します。
係数配列を追加します。
static float AreasupTab_f[INTER_TAB_SIZE2][4][4];
static short AreasupTab_i[INTER_TAB_SIZE2][4][4];

初期化関数を追加します。不要なんですが、体裁だけ整えます。
static inline void interpolateAreasup( float x, float* coeffs )
{
    coeffs[0] = coeffs[1] = coeffs[2] = coeffs[3] = 0;
}

cv::INTER_AREAが指定された時に初期化されるコードを追加します
static const void* initInterTab2D( int method, bool fixpt )
{
    static bool inittab[INTER_MAX+1] = {false};
    float* tab = 0;
    short* itab = 0;
    int ksize = 0;
    if( method == INTER_LINEAR )
        tab = BilinearTab_f[0][0], itab = BilinearTab_i[0][0], ksize=2;
    else if( method == INTER_CUBIC )
        tab = BicubicTab_f[0][0], itab = BicubicTab_i[0][0], ksize=4;
    else if( method == INTER_LANCZOS4 )
        tab = Lanczos4Tab_f[0][0], itab = Lanczos4Tab_i[0][0], ksize=8;
        // --> ここから
    else if( method == INTER_AREA )
        tab = AreasupTab_f[0][0], itab = AreasupTab_i[0][0], ksize=4;
        // --> ここまで追加
    else
        CV_Error( CV_StsBadArg, "Unknown/unsupported interpolation type" );
    //...

係数配列を全部初期化する関数にもcv::INTER_AREA の初期化を追加します
static bool initAllInterTab2D()
{
    return  initInterTab2D( INTER_LINEAR, false ) &&
            initInterTab2D( INTER_LINEAR, true ) &&
            initInterTab2D( INTER_CUBIC, false ) &&
            initInterTab2D( INTER_CUBIC, true ) &&
            initInterTab2D( INTER_LANCZOS4, false ) &&
            initInterTab2D( INTER_LANCZOS4, true ) &&
            initInterTab2D( INTER_AREA, false ) && // この辺
            initInterTab2D( INTER_AREA, true );
}

処理の核心部を追加します
template<class CastOp, typename AT, int ONE>
static void remapAreasup( const Mat& _src, Mat& _dst, const Mat& _xy,
                          const Mat& _fxy, const void* _wtab,
                          int borderType, const Scalar& _borderValue )
{
    typedef typename CastOp::rtype T;
    typedef typename CastOp::type1 WT;
    Size ssize = _src.size(), dsize = _dst.size();
    const int cn = _src.channels();
    const AT* wtab = (const AT*)_wtab;
    const T* S0 = _src.ptr<T>();
    size_t sstep = _src.step/sizeof(S0[0]);
    T cval[CV_CN_MAX];
    CastOp castOp;

    for(int k = 0; k < cn; k++ )
        cval[k] = saturate_cast<T>(_borderValue[k & 3]);

    int borderType1 = borderType != BORDER_TRANSPARENT ? borderType : BORDER_REFLECT_101;

    unsigned width1 = std::max(ssize.width-3, 0), height1 = std::max(ssize.height-3, 0);

    if( _dst.isContinuous() && _xy.isContinuous() && _fxy.isContinuous() )
    {
        dsize.width *= dsize.height;
        dsize.height = 1;
    }

    for(int dy = 0; dy < dsize.height; dy++ )
    {
        T* D = _dst.ptr<T>(dy);
        const short* XY = _xy.ptr<short>(dy);
        const ushort* FXY = _fxy.ptr<ushort>(dy);

        for(int dx = 0; dx < dsize.width; dx++, D += cn )
        {
            int sx = XY[dx*2]-1, sy = XY[dx*2+1]-1;
            //const AT* w = wtab + FXY[dx]*16;
            if( (unsigned)sx < width1 && (unsigned)sy < height1 )
            {
                const T* S = S0 + sy*sstep + sx*cn;
                for(int k = 0; k < cn; k++ )
                {
                    WT sum = std::min<WT>(std::min<WT>(S[0],S[cn]), std::min<WT>(S[cn*2], S[cn*3]));
                    S += sstep;
                    sum = std::min<WT>(sum, std::min<WT>(std::min(S[0],S[cn]), std::min<WT>(S[cn*2], S[cn*3])));
                    S += sstep;
                    sum = std::min<WT>(sum, std::min<WT>(std::min(S[0],S[cn]), std::min<WT>(S[cn*2], S[cn*3])));
                    S += sstep;
                    sum = std::min<WT>(sum, std::min<WT>(std::min(S[0],S[cn]), std::min<WT>(S[cn*2], S[cn*3])));
                    S += 1 - sstep*3;
                    D[k] = static_cast<WT>(sum);
                }
            }
            else
            {
                int x[4], y[4];
                if( borderType == BORDER_TRANSPARENT &&
                    ((unsigned)(sx+1) >= (unsigned)ssize.width ||
                    (unsigned)(sy+1) >= (unsigned)ssize.height) ) {
                      continue;
                }

                if( borderType1 == BORDER_CONSTANT &&
                    (sx >= ssize.width || sx+4 <= 0 ||
                    sy >= ssize.height || sy+4 <= 0))
                {
                    for(int k = 0; k < cn; k++ )
                        D[k] = cval[k];
                    continue;
                }

                for(int i = 0; i < 4; i++ )
                {
                    x[i] = borderInterpolate(sx + i, ssize.width, borderType1)*cn;
                    y[i] = borderInterpolate(sy + i, ssize.height, borderType1);
                }

                for(int k = 0; k < cn; k++, S0++ )
                {
                  WT cv = cval[k], sum = std::numeric_limits<WT>::max();
                    for(int i = 0; i < 4; i++ )
                    {
                        int yi = y[i];
                        const T* S = S0 + yi*sstep;
                        if( yi < 0 )
                            continue;
                        if( x[0] >= 0 )
                            sum = std::min<WT>(sum, S[x[0]]);
                        if( x[1] >= 0 )
                            sum = std::min<WT>(sum, S[x[1]]);
                        if( x[2] >= 0 )
                            sum = std::min<WT>(sum, S[x[2]]);
                        if( x[3] >= 0 )
                            sum = std::min<WT>(sum, S[x[3]]);
                    }
                    D[k] = static_cast<WT>(sum);
                }
                S0 -= cn;
            }
        }
    }
}

remap に cv::INTER_AREA 時の処理を追加します。IPPとかは無視です。
void cv::remap( InputArray _src, OutputArray _dst,
                InputArray _map1, InputArray _map2,
                int interpolation, int borderType, const Scalar& borderValue )
{
    CV_INSTRUMENT_REGION();

    static RemapNNFunc nn_tab[] =
    {
        remapNearest<uchar>, remapNearest<schar>, remapNearest<ushort>, remapNearest<short>,
        remapNearest<int>, remapNearest<float>, remapNearest<double>, 0
    };

    static RemapFunc linear_tab[] =
    {
        remapBilinear<FixedPtCast<int, uchar, INTER_REMAP_COEF_BITS>, RemapVec_8u, short>, 0,
        remapBilinear<Cast<float, ushort>, RemapNoVec, float>,
        remapBilinear<Cast<float, short>, RemapNoVec, float>, 0,
        remapBilinear<Cast<float, float>, RemapNoVec, float>,
        remapBilinear<Cast<double, double>, RemapNoVec, float>, 0
    };

    static RemapFunc cubic_tab[] =
    {
        remapBicubic<FixedPtCast<int, uchar, INTER_REMAP_COEF_BITS>, short, INTER_REMAP_COEF_SCALE>, 0,
        remapBicubic<Cast<float, ushort>, float, 1>,
        remapBicubic<Cast<float, short>, float, 1>, 0,
        remapBicubic<Cast<float, float>, float, 1>,
        remapBicubic<Cast<double, double>, float, 1>, 0
    };

    static RemapFunc lanczos4_tab[] =
    {
        remapLanczos4<FixedPtCast<int, uchar, INTER_REMAP_COEF_BITS>, short, INTER_REMAP_COEF_SCALE>, 0,
        remapLanczos4<Cast<float, ushort>, float, 1>,
        remapLanczos4<Cast<float, short>, float, 1>, 0,
        remapLanczos4<Cast<float, float>, float, 1>,
        remapLanczos4<Cast<double, double>, float, 1>, 0
    };

    // --> ここから
    static RemapFunc areasup_tab[] = 
    {
        remapAreasup<FixedPtCast<int, uchar, INTER_REMAP_COEF_BITS>, short, INTER_REMAP_COEF_SCALE>, 0,
        remapAreasup<Cast<float, ushort>, float, 1>,
        remapAreasup<Cast<float, short>, float, 1>, 0,
        remapAreasup<Cast<float, float>, float, 1>,
        remapAreasup<Cast<double, double>, float, 1>, 0
    };
    // --> ここまで追加
  
    CV_Assert( !_map1.empty() );
    CV_Assert( _map2.empty() || (_map2.size() == _map1.size()));

    CV_OCL_RUN(_src.dims() <= 2 && _dst.isUMat(),
               ocl_remap(_src, _dst, _map1, _map2, interpolation, borderType, borderValue))

    Mat src = _src.getMat(), map1 = _map1.getMat(), map2 = _map2.getMat();
    _dst.create( map1.size(), src.type() );
    Mat dst = _dst.getMat();


    CV_OVX_RUN(
        src.type() == CV_8UC1 && dst.type() == CV_8UC1 &&
        !ovx::skipSmallImages<VX_KERNEL_REMAP>(src.cols, src.rows) &&
        (borderType& ~BORDER_ISOLATED) == BORDER_CONSTANT &&
        ((map1.type() == CV_32FC2 && map2.empty() && map1.size == dst.size) ||
         (map1.type() == CV_32FC1 && map2.type() == CV_32FC1 && map1.size == dst.size && map2.size == dst.size) ||
         (map1.empty() && map2.type() == CV_32FC2 && map2.size == dst.size)) &&
        ((borderType & BORDER_ISOLATED) != 0 || !src.isSubmatrix()),
        openvx_remap(src, dst, map1, map2, interpolation, borderValue));

    CV_Assert( dst.cols < SHRT_MAX && dst.rows < SHRT_MAX && src.cols < SHRT_MAX && src.rows < SHRT_MAX );

    if( dst.data == src.data )
        src = src.clone();

    // --> 使われていない INTER_AREA を利用するので、以下コメントアウトして有効化
    //if( interpolation == INTER_AREA )
    //    interpolation = INTER_LINEAR;

    int type = src.type(), depth = CV_MAT_DEPTH(type);

#if defined HAVE_IPP && !IPP_DISABLE_REMAP
    CV_IPP_CHECK()
    {
        if ((interpolation == INTER_LINEAR || interpolation == INTER_CUBIC || interpolation == INTER_NEAREST) &&
                map1.type() == CV_32FC1 && map2.type() == CV_32FC1 &&
                (borderType == BORDER_CONSTANT || borderType == BORDER_TRANSPARENT))
        {
            int ippInterpolation =
                interpolation == INTER_NEAREST ? IPPI_INTER_NN :
                interpolation == INTER_LINEAR ? IPPI_INTER_LINEAR : IPPI_INTER_CUBIC;

            ippiRemap ippFunc =
                type == CV_8UC1 ? (ippiRemap)ippiRemap_8u_C1R :
                type == CV_8UC3 ? (ippiRemap)ippiRemap_8u_C3R :
                type == CV_8UC4 ? (ippiRemap)ippiRemap_8u_C4R :
                type == CV_16UC1 ? (ippiRemap)ippiRemap_16u_C1R :
                type == CV_16UC3 ? (ippiRemap)ippiRemap_16u_C3R :
                type == CV_16UC4 ? (ippiRemap)ippiRemap_16u_C4R :
                type == CV_32FC1 ? (ippiRemap)ippiRemap_32f_C1R :
                type == CV_32FC3 ? (ippiRemap)ippiRemap_32f_C3R :
                type == CV_32FC4 ? (ippiRemap)ippiRemap_32f_C4R : 0;

            if (ippFunc)
            {
                bool ok;
                IPPRemapInvoker invoker(src, dst, map1, map2, ippFunc, ippInterpolation,
                                        borderType, borderValue, &ok);
                Range range(0, dst.rows);
                parallel_for_(range, invoker, dst.total() / (double)(1 << 16));

                if (ok)
                {
                    CV_IMPL_ADD(CV_IMPL_IPP|CV_IMPL_MT);
                    return;
                }
                setIppErrorStatus();
            }
        }
    }
#endif

    RemapNNFunc nnfunc = 0;
    RemapFunc ifunc = 0;
    const void* ctab = 0;
    bool fixpt = depth == CV_8U;
    bool planar_input = false;

    if( interpolation == INTER_NEAREST )
    {
        nnfunc = nn_tab[depth];
        CV_Assert( nnfunc != 0 );
    }
    else
    {
        if( interpolation == INTER_LINEAR )
            ifunc = linear_tab[depth];
        else if( interpolation == INTER_CUBIC ){
            ifunc = cubic_tab[depth];
            CV_Assert( _src.channels() <= 4 );
        }
        else if( interpolation == INTER_LANCZOS4 ){
            ifunc = lanczos4_tab[depth];
            CV_Assert( _src.channels() <= 4 );
        }
        // --> ここから
        else if( interpolation == INTER_AREA ) {
            ifunc = areasup_tab[depth];
            CV_Assert( _src.channels() <= 4 );
        }
        // --> ここまで追加
        else
            CV_Error( CV_StsBadArg, "Unknown interpolation method" );
        CV_Assert( ifunc != 0 );
        ctab = initInterTab2D( interpolation, fixpt );
    }

    const Mat *m1 = &map1, *m2 = &map2;

    if( (map1.type() == CV_16SC2 && (map2.type() == CV_16UC1 || map2.type() == CV_16SC1 || map2.empty())) ||
        (map2.type() == CV_16SC2 && (map1.type() == CV_16UC1 || map1.type() == CV_16SC1 || map1.empty())) )
    {
        if( map1.type() != CV_16SC2 )
            std::swap(m1, m2);
    }
    else
    {
        CV_Assert( ((map1.type() == CV_32FC2 || map1.type() == CV_16SC2) && map2.empty()) ||
            (map1.type() == CV_32FC1 && map2.type() == CV_32FC1) );
        planar_input = map1.channels() == 1;
    }

    RemapInvoker invoker(src, dst, m1, m2,
                         borderType, borderValue, planar_input, nnfunc, ifunc,
                         ctab);
    parallel_for_(Range(0, dst.rows), invoker, dst.total()/(double)(1<<16));
}

成果
 変換元画像

 cv::INTER_LINEAR で変換

 魔改造 cv::INTER_AREA で変換

0 件のコメント: