opencv cvSVD思考（测试代码和matlab对比）
2018-02-27 17:55
489 查看
matlab对矩阵A进行奇异值分解//matlab:
A =
1 1 1 1
0 1 0 0
0 0 1 1
0 0 0 0
>> [S V D]=svd(A)
S =
0.8460 0.2261 0.4828 0
0.1922 0.7154 -0.6718 0
0.4973 -0.6611 -0.5618 0
0 0 0 1.0000
V =
2.3244 0 0 0
0 1.1472 0 0
0 0 0.5304 0
0 0 0 0
D =
0.3640 0.1971 0.9103 0
0.4467 0.8207 -0.3563 0
0.5779 -0.3792 -0.1490 -0.7071
0.5779 -0.3792 -0.1490 0.7071
*/opencv cvSVD 原理
opencv 奇异值分解（SVD）测试代码：
需要配置opencv1.0和vc6
_cxcore_my.h
#ifndef _CXCORE_INTERNAL_H_
#define _CXCORE_INTERNAL_H_
#if defined _MSC_VER && _MSC_VER >= 1200
/* disable warnings related to inline functions */
#pragma warning( disable: 4711 4710 4514 )
#endif
typedef unsigned long ulong;
#ifdef __BORLANDC__
#define WIN32
#define CV_DLL
#undef _CV_ALWAYS_PROFILE_
#define _CV_ALWAYS_NO_PROFILE_
#endif
#include "cxcore.h"
#include "cxmisc.h"
#include "_cxipp_my.h"
#include <math.h>
#include <assert.h>
#include <string.h>
#include <stdlib.h>
#include <stdio.h>
#include <limits.h>
#include <float.h>
// -128.f ... 255.f
extern const float icv8x32fTab[];
#define CV_8TO32F(x) icv8x32fTab[(x)+128]
extern const ushort icv8x16uSqrTab[];
#define CV_SQR_8U(x) icv8x16uSqrTab[(x)+255]
extern const char* icvHersheyGlyphs[];
extern const signed char icvDepthToType[];
#define icvIplToCvDepth( depth ) \
icvDepthToType[(((depth) & 255) >> 2) + ((depth) < 0)]
extern const uchar icvSaturate8u[];
#define CV_FAST_CAST_8U(t) (assert(-256 <= (t) && (t) <= 512), icvSaturate8u[(t)+256])
#define CV_MIN_8U(a,b) ((a) - CV_FAST_CAST_8U((a) - (b)))
#define CV_MAX_8U(a,b) ((a) + CV_FAST_CAST_8U((b) - (a)))
typedef CvFunc2D_3A1I CvArithmBinMaskFunc2D;
typedef CvFunc2D_2A1P1I CvArithmUniMaskFunc2D;
/****************************************************************************************\
* Complex arithmetics *
\****************************************************************************************/
struct CvComplex32f;
struct CvComplex64f;
struct CvComplex32f
{
float re, im;
CvComplex32f() {}
CvComplex32f( float _re, float _im=0 ) : re(_re), im(_im) {}
explicit CvComplex32f( const CvComplex64f& v );
//CvComplex32f( const CvComplex32f& v ) : re(v.re), im(v.im) {}
//CvComplex32f& operator = (const CvComplex32f& v ) { re = v.re; im = v.im; return *this; }
operator CvComplex64f() const;
};
struct CvComplex64f
{
double re, im;
CvComplex64f() {}
CvComplex64f( double _re, double _im=0 ) : re(_re), im(_im) {}
explicit CvComplex64f( const CvComplex32f& v );
//CvComplex64f( const CvComplex64f& v ) : re(v.re), im(v.im) {}
//CvComplex64f& operator = (const CvComplex64f& v ) { re = v.re; im = v.im; return *this; }
operator CvComplex32f() const;
};
inline CvComplex32f::CvComplex32f( const CvComplex64f& v ) : re((float)v.re), im((float)v.im) {}
inline CvComplex64f::CvComplex64f( const CvComplex32f& v ) : re(v.re), im(v.im) {}
inline CvComplex32f operator + (CvComplex32f a, CvComplex32f b)
{
return CvComplex32f( a.re + b.re, a.im + b.im );
}
inline CvComplex32f& operator += (CvComplex32f& a, CvComplex32f b)
{
a.re += b.re;
a.im += b.im;
return a;
}
inline CvComplex32f operator - (CvComplex32f a, CvComplex32f b)
{
return CvComplex32f( a.re - b.re, a.im - b.im );
}
inline CvComplex32f& operator -= (CvComplex32f& a, CvComplex32f b)
{
a.re -= b.re;
a.im -= b.im;
return a;
}
inline CvComplex32f operator - (CvComplex32f a)
{
return CvComplex32f( -a.re, -a.im );
}
inline CvComplex32f operator * (CvComplex32f a, CvComplex32f b)
{
return CvComplex32f( a.re*b.re - a.im*b.im, a.re*b.im + a.im*b.re );
}
inline double abs(CvComplex32f a)
{
return sqrt( (double)a.re*a.re + (double)a.im*a.im );
}
inline CvComplex32f conj(CvComplex32f a)
{
return CvComplex32f( a.re, -a.im );
}
inline CvComplex32f operator / (CvComplex32f a, CvComplex32f b)
{
double t = 1./((double)b.re*b.re + (double)b.im*b.im);
return CvComplex32f( (float)((a.re*b.re + a.im*b.im)*t),
(float)((-a.re*b.im + a.im*b.re)*t) );
}
inline CvComplex32f operator * (double a, CvComplex32f b)
{
return CvComplex32f( (float)(a*b.re), (float)(a*b.im) );
}
inline CvComplex32f operator * (CvComplex32f a, double b)
{
return CvComplex32f( (float)(a.re*b), (float)(a.im*b) );
}
inline CvComplex32f::operator CvComplex64f() const
{
return CvComplex64f(re,im);
}
inline CvComplex64f operator + (CvComplex64f a, CvComplex64f b)
{
return CvComplex64f( a.re + b.re, a.im + b.im );
}
inline CvComplex64f& operator += (CvComplex64f& a, CvComplex64f b)
{
a.re += b.re;
a.im += b.im;
return a;
}
inline CvComplex64f operator - (CvComplex64f a, CvComplex64f b)
{
return CvComplex64f( a.re - b.re, a.im - b.im );
}
inline CvComplex64f& operator -= (CvComplex64f& a, CvComplex64f b)
{
a.re -= b.re;
a.im -= b.im;
return a;
}
inline CvComplex64f operator - (CvComplex64f a)
{
return CvComplex64f( -a.re, -a.im );
}
inline CvComplex64f operator * (CvComplex64f a, CvComplex64f b)
{
return CvComplex64f( a.re*b.re - a.im*b.im, a.re*b.im + a.im*b.re );
}
inline double abs(CvComplex64f a)
{
return sqrt( (double)a.re*a.re + (double)a.im*a.im );
}
inline CvComplex64f operator / (CvComplex64f a, CvComplex64f b)
{
double t = 1./((double)b.re*b.re + (double)b.im*b.im);
return CvComplex64f( (a.re*b.re + a.im*b.im)*t,
(-a.re*b.im + a.im*b.re)*t );
}
inline CvComplex64f operator * (double a, CvComplex64f b)
{
return CvComplex64f( a*b.re, a*b.im );
}
inline CvComplex64f operator * (CvComplex64f a, double b)
{
return CvComplex64f( a.re*b, a.im*b );
}
inline CvComplex64f::operator CvComplex32f() const
{
return CvComplex32f((float)re,(float)im);
}
inline CvComplex64f conj(CvComplex64f a)
{
return CvComplex64f( a.re, -a.im );
}
inline CvComplex64f operator + (CvComplex64f a, CvComplex32f b)
{
return CvComplex64f( a.re + b.re, a.im + b.im );
}
inline CvComplex64f operator + (CvComplex32f a, CvComplex64f b)
{
return CvComplex64f( a.re + b.re, a.im + b.im );
}
inline CvComplex64f operator - (CvComplex64f a, CvComplex32f b)
{
return CvComplex64f( a.re - b.re, a.im - b.im );
}
inline CvComplex64f operator - (CvComplex32f a, CvComplex64f b)
{
return CvComplex64f( a.re - b.re, a.im - b.im );
}
inline CvComplex64f operator * (CvComplex64f a, CvComplex32f b)
{
return CvComplex64f( a.re*b.re - a.im*b.im, a.re*b.im + a.im*b.re );
}
inline CvComplex64f operator * (CvComplex32f a, CvComplex64f b)
{
return CvComplex64f( a.re*b.re - a.im*b.im, a.re*b.im + a.im*b.re );
}
typedef CvStatus (CV_STDCALL * CvCopyMaskFunc)(const void* src, int src_step,
void* dst, int dst_step, CvSize size,
const void* mask, int mask_step);
CvCopyMaskFunc icvGetCopyMaskFunc( int elem_size );
CvStatus CV_STDCALL icvSetZero_8u_C1R( uchar* dst, int dststep, CvSize size );
CvStatus CV_STDCALL icvScale_32f( const float* src, float* dst, int len, float a, float b );
CvStatus CV_STDCALL icvScale_64f( const double* src, double* dst, int len, double a, double b );
CvStatus CV_STDCALL icvLUT_Transform8u_8u_C1R( const uchar* src, int srcstep, uchar* dst,
int dststep, CvSize size, const uchar* lut );
CvStatus CV_STDCALL icvLUT_Transform8u_16u_C1R( const uchar* src, int srcstep, ushort* dst,
int dststep, CvSize size, const ushort* lut );
CvStatus CV_STDCALL icvLUT_Transform8u_32s_C1R( con
4000
st uchar* src, int srcstep, int* dst,
int dststep, CvSize size, const int* lut );
CvStatus CV_STDCALL icvLUT_Transform8u_64f_C1R( const uchar* src, int srcstep, double* dst,
int dststep, CvSize size, const double* lut );
CvStatus CV_STDCALL icvLUT_Transform8u_8u_C2R( const uchar* src, int srcstep, uchar* dst,
int dststep, CvSize size, const uchar* lut );
CvStatus CV_STDCALL icvLUT_Transform8u_8u_C3R( const uchar* src, int srcstep, uchar* dst,
int dststep, CvSize size, const uchar* lut );
CvStatus CV_STDCALL icvLUT_Transform8u_8u_C4R( const uchar* src, int srcstep, uchar* dst,
int dststep, CvSize size, const uchar* lut );
typedef CvStatus (CV_STDCALL * CvLUT_TransformFunc)( const void* src, int srcstep, void* dst,
int dststep, CvSize size, const void* lut );
CV_INLINE CvStatus
icvLUT_Transform8u_8s_C1R( const uchar* src, int srcstep, char* dst,
int dststep, CvSize size, const char* lut )
{
return icvLUT_Transform8u_8u_C1R( src, srcstep, (uchar*)dst,
dststep, size, (const uchar*)lut );
}
CV_INLINE CvStatus
icvLUT_Transform8u_16s_C1R( const uchar* src, int srcstep, short* dst,
int dststep, CvSize size, const short* lut )
{
return icvLUT_Transform8u_16u_C1R( src, srcstep, (ushort*)dst,
dststep, size, (const ushort*)lut );
}
CV_INLINE CvStatus
icvLUT_Transform8u_32f_C1R( const uchar* src, int srcstep, float* dst,
int dststep, CvSize size, const float* lut )
{
return icvLUT_Transform8u_32s_C1R( src, srcstep, (int*)dst,
dststep, size, (const int*)lut );
}
#endif /*_CXCORE_INTERNAL_H_*/ _cxipp_my.h#ifndef _CXCORE_IPP_H_
#define _CXCORE_IPP_H_
/****************************************************************************************\
* Copy/Set *
\****************************************************************************************/
/* temporary disable ipp zero and copy functions as they affect subsequent functions' performance */
IPCVAPI_EX( CvStatus, icvCopy_8u_C1R, "ippiCopy_8u_C1R", 0/*CV_PLUGINS1(CV_PLUGIN_IPPI)*/,
( const uchar* src, int src_step,
uchar* dst, int dst_step, CvSize size ))
IPCVAPI_EX( CvStatus, icvSetByte_8u_C1R, "ippiSet_8u_C1R", 0/*CV_PLUGINS1(CV_PLUGIN_IPPI)*/,
( uchar value, uchar* dst, int dst_step, CvSize size ))
IPCVAPI_EX( CvStatus, icvCvt_32f64f, "ippsConvert_32f64f",
CV_PLUGINS1(CV_PLUGIN_IPPS), ( const float* src, double* dst, int len ))
IPCVAPI_EX( CvStatus, icvCvt_64f32f, "ippsConvert_64f32f",
CV_PLUGINS1(CV_PLUGIN_IPPS), ( const double* src, float* dst, int len ))
#define IPCV_COPYSET( flavor, arrtype, scalartype ) \
IPCVAPI_EX( CvStatus, icvCopy##flavor, "ippiCopy" #flavor, \
CV_PLUGINS1(CV_PLUGIN_IPPI), \
( const arrtype* src, int srcstep, \
arrtype* dst, int dststep, CvSize size, \
const uchar* mask, int maskstep )) \
IPCVAPI_EX( CvStatus, icvSet##flavor, "ippiSet" #flavor, \
0/*CV_PLUGINS1(CV_PLUGIN_OPTCV)*/, \
( arrtype* dst, int dststep, \
const uchar* mask, int maskstep, \
CvSize size, const arrtype* scalar ))
IPCV_COPYSET( _8u_C1MR, uchar, int )
IPCV_COPYSET( _16s_C1MR, ushort, int )
IPCV_COPYSET( _8u_C3MR, uchar, int )
IPCV_COPYSET( _8u_C4MR, int, int )
IPCV_COPYSET( _16s_C3MR, ushort, int )
IPCV_COPYSET( _16s_C4MR, int64, int64 )
IPCV_COPYSET( _32f_C3MR, int, int )
IPCV_COPYSET( _32f_C4MR, int, int )
IPCV_COPYSET( _64s_C3MR, int64, int64 )
IPCV_COPYSET( _64s_C4MR, int64, int64 )
/****************************************************************************************\
* Arithmetics *
\****************************************************************************************/
#define IPCV_BIN_ARITHM( name ) \
IPCVAPI_EX( CvStatus, icv##name##_8u_C1R, \
"ippi" #name "_8u_C1RSfs", CV_PLUGINS1(CV_PLUGIN_IPPI), \
( const uchar* src1, int srcstep1, const uchar* src2, int srcstep2, \
uchar* dst, int dststep, CvSize size, int scalefactor )) \
IPCVAPI_EX( CvStatus, icv##name##_16u_C1R, \
"ippi" #name "_16u_C1RSfs", CV_PLUGINS1(CV_PLUGIN_IPPI), \
( const ushort* src1, int srcstep1, const ushort* src2, int srcstep2,\
ushort* dst, int dststep, CvSize size, int scalefactor )) \
IPCVAPI_EX( CvStatus, icv##name##_16s_C1R, \
"ippi" #name "_16s_C1RSfs", CV_PLUGINS1(CV_PLUGIN_IPPI), \
( const short* src1, int srcstep1, const short* src2, int srcstep2, \
short* dst, int dststep, CvSize size, int scalefactor )) \
IPCVAPI_EX( CvStatus, icv##name##_32s_C1R, \
"ippi" #name "_32s_C1R", CV_PLUGINS1(CV_PLUGIN_IPPI), \
( const int* src1, int srcstep1, const int* src2, int srcstep2, \
int* dst, int dststep, CvSize size )) \
IPCVAPI_EX( CvStatus, icv##name##_32f_C1R, \
"ippi" #name "_32f_C1R", CV_PLUGINS1(CV_PLUGIN_IPPI), \
( const float* src1, int srcstep1, const float* src2, int srcstep2, \
float* dst, int dststep, CvSize size )) \
IPCVAPI_EX( CvStatus, icv##name##_64f_C1R, \
"ippi" #name "_64f_C1R", CV_PLUGINS1(CV_PLUGIN_IPPI), \
( const double* src1, int srcstep1, const double* src2, int srcstep2,\
double* dst, int dststep, CvSize size ))
IPCV_BIN_ARITHM( Add )
IPCV_BIN_ARITHM( Sub )
#undef IPCV_BIN_ARITHM
/****************************************************************************************\
* Logical operations *
\****************************************************************************************/
#define IPCV_LOGIC( name ) \
IPCVAPI_EX( CvStatus, icv##name##_8u_C1R, \
"ippi" #name "_8u_C1R", 0/*CV_PLUGINS1(CV_PLUGIN_IPPI)*/, \
( const uchar* src1, int srcstep1, const uchar* src2, int srcstep2, \
uchar* dst, int dststep, CvSize size ))
IPCV_LOGIC( And )
IPCV_LOGIC( Or )
IPCV_LOGIC( Xor )
#undef IPCV_LOGIC
IPCVAPI_EX( CvStatus, icvNot_8u_C1R, "ippiNot_8u_C1R", CV_PLUGINS1(CV_PLUGIN_IPPI),
( const uchar* src, int step1, uchar* dst, int step, CvSize size ))
/****************************************************************************************\
* Image Statistics *
\****************************************************************************************/
///////////////////////////////////////// Mean //////////////////////////////////////////
#define IPCV_DEF_MEAN_MASK( flavor, srctype ) \
IPCVAPI_EX( CvStatus, icvMean_##flavor##_C1MR, \
"ippiMean_" #flavor "_C1MR", CV_PLUGINS1(CV_PLUGIN_IPPCV), \
( const srctype* img, int imgstep, const uchar* mask, \
int maskStep, CvSize size, double* mean )) \
IPCVAPI_EX( CvStatus, icvMean_##flavor##_C2MR, \
"ippiMean_" #flavor "_C2MR", 0/*CV_PLUGINS1(CV_PLUGIN_OPTCV)*/, \
( const srctype* img, int imgstep, const uchar* mask, \
int maskStep, CvSize size, double* mean )) \
IPCVAPI_EX( CvStatus, icvMean_##flavor##_C3MR, \
"ippiMean_" #flavor "_C3MR", 0/*CV_PLUGINS1(CV_PLUGIN_OPTCV)*/, \
( const srctype* img, int imgstep, const uchar* mask, \
int maskStep, CvSize size, double* mean )) \
IPCVAPI_EX( CvStatus, icvMean_##flavor##_C4MR, \
"ippiMean_" #flavor "_C4MR", 0/*CV_PLUGINS1(CV_PLUGIN_OPTCV)*/, \
( const srctype* img, int imgstep, const uchar* mask, \
int maskStep, CvSize size, double* mean ))
IPCV_DEF_MEAN_MASK( 8u, uchar )
IPCV_DEF_MEAN_MASK( 16u, ushort )
IPCV_DEF_MEAN_MASK( 16s, short )
IPCV_DEF_MEAN_MASK( 32s, int )
IPCV_DEF_MEAN_MASK( 32f, float )
IPCV_DEF_MEAN_MASK( 64f, double )
#undef IPCV_DEF_MEAN_MASK
//////////////////////////////////// Mean_StdDev ////////////////////////////////////////
#undef IPCV_MEAN_SDV_PLUGIN
#define ICV_MEAN_SDV_PLUGIN 0 /* CV_PLUGINS1(IPPCV) */
#define IPCV_DEF_MEAN_SDV( flavor, srctype ) \
IPCVAPI_EX( CvStatus, icvMean_StdDev_##flavor##_C1R, \
"ippiMean_StdDev_" #flavor "_C1R", ICV_MEAN_SDV_PLUGIN, \
( const srctype* img, int imgstep, CvSize size, double* mean, double* sdv ))\
IPCVAPI_EX( CvStatus, icvMean_StdDev_##flavor##_C2R, \
"ippiMean_StdDev_" #flavor "_C2R", ICV_MEAN_SDV_PLUGIN, \
( const srctype* img, int imgstep, CvSize size, double* mean, double* sdv ))\
IPCVAPI_EX( CvStatus, icvMean_StdDev_##flavor##_C3R, \
"ippiMean_StdDev_" #flavor "_C3R", ICV_MEAN_SDV_PLUGIN, \
( const srctype* img, int imgstep, CvSize size, double* mean, double* sdv ))\
IPCVAPI_EX( CvStatus, icvMean_StdDev_##flavor##_C4R, \
"ippiMean_StdDev_" #flavor "_C4R", ICV_MEAN_SDV_PLUGIN, \
( const srctype* img, int imgstep, CvSize size, double* mean, double* sdv ))\
\
IPCVAPI_EX( CvStatus, icvMean_StdDev_##flavor##_C1MR, \
"ippiMean_StdDev_" #flavor "_C1MR", ICV_MEAN_SDV_PLUGIN, \
( const srctype* img, int imgstep, \
const uchar* mask, int maskStep, \
CvSize size, double* mean, double* sdv )) \
IPCVAPI_EX( CvStatus, icvMean_StdDev_##flavor##_C2MR, \
"ippiMean_StdDev_" #flavor "_C2MR", ICV_MEAN_SDV_PLUGIN, \
( const srctype* img, int imgstep, const uchar* mask, int maskStep, \
CvSize size, double* mean, double* sdv )) \
IPCVAPI_EX( CvStatus, icvMean_StdDev_##flavor##_C3MR, \
"ippiMean_StdDev_" #flavor "_C3MR", ICV_MEAN_SDV_PLUGIN, \
( const srctype* img, int imgstep, \
const uchar* mask, int maskStep, \
CvSize size, double* mean, double* sdv )) \
IPCVAPI_EX( CvStatus, icvMean_StdDev_##flavor##_C4MR, \
"ippiMean_StdDev_" #flavor "_C4MR", ICV_MEAN_SDV_PLUGIN, \
( const srctype* img, int imgstep, \
const uchar* mask, int maskStep, \
CvSize size, double* mean, double* sdv ))
IPCV_DEF_MEAN_SDV( 8u, uchar )
IPCV_DEF_MEAN_SDV( 16u, ushort )
IPCV_DEF_MEAN_SDV( 16s, short )
IPCV_DEF_MEAN_SDV( 32s, int )
IPCV_DEF_MEAN_SDV( 32f, float )
IPCV_DEF_MEAN_SDV( 64f, double )
#undef IPCV_DEF_MEAN_SDV
#undef IPCV_MEAN_SDV_PLUGIN
//////////////////////////////////// MinMaxIndx /////////////////////////////////////////
#define IPCV_DEF_MIN_MAX_LOC( flavor, srctype, extrtype ) \
IPCVAPI_EX( CvStatus, icvMinMaxIndx_##flavor##_C1R, \
"ippiMinMaxIndx_" #flavor "_C1R", CV_PLUGINS1(CV_PLUGIN_IPPCV), \
( const srctype* img, int imgstep, \
CvSize size, extrtype* minVal, extrtype* maxVal, \
CvPoint* minLoc, CvPoint* maxLoc )) \
\
IPCVAPI_EX( CvStatus, icvMinMaxIndx_##flavor##_C1MR, \
"ippiMinMaxIndx_" #flavor "_C1MR", CV_PLUGINS1(CV_PLUGIN_IPPCV),\
( const srctype* img, int imgstep, \
const uchar* mask, int maskStep, \
CvSize size, extrtype* minVal, extrtype* maxVal, \
CvPoint* minLoc, CvPoint* maxLoc ))
IPCV_DEF_MIN_MAX_LOC( 8u, uchar, float )
IPCV_DEF_MIN_MAX_LOC( 16u, ushort, float )
IPCV_DEF_MIN_MAX_LOC( 16s, short, float )
IPCV_DEF_MIN_MAX_LOC( 32s, int, double )
IPCV_DEF_MIN_MAX_LOC( 32f, int, float )
IPCV_DEF_MIN_MAX_LOC( 64f, int64, double )
#undef IPCV_DEF_MIN_MAX_LOC
////////////////////////////////////////// Sum //////////////////////////////////////////
#define IPCV_DEF_SUM_NOHINT( flavor, srctype ) \
IPCVAPI_EX( CvStatus, icvSum_##flavor##_C1R, \
"ippiSum_" #flavor "_C1R", CV_PLUGINS1(CV_PLUGIN_IPPI), \
( const srctype* img, int imgstep, \
CvSize size, double* sum )) \
IPCVAPI_EX( CvStatus, icvSum_##flavor##_C2R, \
"ippiSum_" #flavor "_C2R", CV_PLUGINS1(CV_PLUGIN_IPPI), \
( const srctype* img, int imgstep, \
CvSize size, double* sum )) \
IPCVAPI_EX( CvStatus, icvSum_##flavor##_C3R, \
"ippiSum_" #flavor "_C3R", CV_PLUGINS1(CV_PLUGIN_IPPI), \
( const srctype* img, int imgstep, \
CvSize size, double* sum )) \
IPCVAPI_EX( CvStatus, icvSum_##flavor##_C4R, \
"ippiSum_" #flavor "_C4R", CV_PLUGINS1(CV_PLUGIN_IPPI), \
( const srctype* img, int imgstep, \
CvSize size, double* sum ))
#define IPCV_DEF_SUM_NOHINT_NO_IPP( flavor, srctype ) \
IPCVAPI_EX( CvStatus, icvSum_##flavor##_C1R, \
"ippiSum_" #flavor "_C1R", 0, \
( const srctype* img, int imgstep, \
CvSize size, double* sum )) \
IPCVAPI_EX( CvStatus, icvSum_##flavor##_C2R, \
"ippiSum_" #flavor "_C2R", 0, \
( const srctype* img, int imgstep, \
CvSize size, double* sum )) \
IPCVAPI_EX( CvStatus, icvSum_##flavor##_C3R, \
"ippiSum_" #flavor "_C3R", 0, \
( const srctype* img, int imgstep, \
CvSize size, double* sum )) \
IPCVAPI_EX( CvStatus, icvSum_##flavor##_C4R, \
"ippiSum_" #flavor "_C4R", 0, \
( const srctype* img, int imgstep, \
CvSize size, double* sum ))
#define IPCV_DEF_SUM_HINT( flavor, srctype ) \
IPCVAPI_EX( CvStatus, icvSum_##flavor##_C1R, \
"ippiSum_" #flavor "_C1R", CV_PLUGINS1(CV_PLUGIN_IPPI), \
( const srctype* img, int imgstep, \
CvSize size, double* sum, CvHintAlgorithm )) \
IPCVAPI_EX( CvStatus, icvSum_##flavor##_C2R, \
"ippiSum_" #flavor "_C2R", CV_PLUGINS1(CV_PLUGIN_IPPI)
17a7a
, \
( const srctype* img, int imgstep, \
CvSize size, double* sum, CvHintAlgorithm )) \
IPCVAPI_EX( CvStatus, icvSum_##flavor##_C3R, \
"ippiSum_" #flavor "_C3R", CV_PLUGINS1(CV_PLUGIN_IPPI), \
( const srctype* img, int imgstep, \
CvSize size, double* sum, CvHintAlgorithm )) \
IPCVAPI_EX( CvStatus, icvSum_##flavor##_C4R, \
"ippiSum_" #flavor "_C4R", CV_PLUGINS1(CV_PLUGIN_IPPI), \
( const srctype* img, int imgstep, \
CvSize size, double* sum, CvHintAlgorithm ))
IPCV_DEF_SUM_NOHINT( 8u, uchar )
IPCV_DEF_SUM_NOHINT( 16s, short )
IPCV_DEF_SUM_NOHINT_NO_IPP( 16u, ushort )
IPCV_DEF_SUM_NOHINT( 32s, int )
IPCV_DEF_SUM_HINT( 32f, float )
IPCV_DEF_SUM_NOHINT( 64f, double )
#undef IPCV_DEF_SUM_NOHINT
#undef IPCV_DEF_SUM_HINT
////////////////////////////////////////// CountNonZero /////////////////////////////////
#define IPCV_DEF_NON_ZERO( flavor, srctype ) \
IPCVAPI_EX( CvStatus, icvCountNonZero_##flavor##_C1R, \
"ippiCountNonZero_" #flavor "_C1R", 0/*CV_PLUGINS1(CV_PLUGIN_OPTCV)*/, \
( const srctype* img, int imgstep, CvSize size, int* nonzero ))
IPCV_DEF_NON_ZERO( 8u, uchar )
IPCV_DEF_NON_ZERO( 16s, ushort )
IPCV_DEF_NON_ZERO( 32s, int )
IPCV_DEF_NON_ZERO( 32f, int )
IPCV_DEF_NON_ZERO( 64f, int64 )
#undef IPCV_DEF_NON_ZERO
////////////////////////////////////////// Norms /////////////////////////////////
#define IPCV_DEF_NORM_NOHINT_C1( flavor, srctype ) \
IPCVAPI_EX( CvStatus, icvNorm_Inf_##flavor##_C1R, \
"ippiNorm_Inf_" #flavor "_C1R", CV_PLUGINS1(CV_PLUGIN_IPPI), \
( const srctype* img, int imgstep, \
CvSize size, double* norm )) \
IPCVAPI_EX( CvStatus, icvNorm_L1_##flavor##_C1R, \
"ippiNorm_L1_" #flavor "_C1R", CV_PLUGINS1(CV_PLUGIN_IPPI), \
( const srctype* img, int imgstep, \
CvSize size, double* norm )) \
IPCVAPI_EX( CvStatus, icvNorm_L2_##flavor##_C1R, \
"ippiNorm_L2_" #flavor "_C1R", 0/*CV_PLUGINS1(CV_PLUGIN_IPPI)*/, \
( const srctype* img, int imgstep, \
CvSize size, double* norm )) \
IPCVAPI_EX( CvStatus, icvNormDiff_Inf_##flavor##_C1R, \
"ippiNormDiff_Inf_" #flavor "_C1R", 0/*CV_PLUGINS1(CV_PLUGIN_IPPI)*/, \
( const srctype* img1, int imgstep1, \
const srctype* img2, int imgstep2, \
CvSize size, double* norm )) \
IPCVAPI_EX( CvStatus, icvNormDiff_L1_##flavor##_C1R, \
"ippiNormDiff_L1_" #flavor "_C1R", CV_PLUGINS1(CV_PLUGIN_IPPI), \
( const srctype* img1, int imgstep1, \
const srctype* img2, int imgstep2, \
CvSize size, double* norm )) \
IPCVAPI_EX( CvStatus, icvNormDiff_L2_##flavor##_C1R, \
"ippiNormDiff_L2_" #flavor "_C1R", 0/*CV_PLUGINS1(CV_PLUGIN_IPPI)*/, \
( const srctype* img1, int imgstep1, \
const srctype* img2, int imgstep2, \
CvSize size, double* norm ))
#define IPCV_DEF_NORM_HINT_C1( flavor, srctype ) \
IPCVAPI_EX( CvStatus, icvNorm_Inf_##flavor##_C1R, \
"ippiNorm_Inf_" #flavor "_C1R", CV_PLUGINS1(CV_PLUGIN_IPPI), \
( const srctype* img, int imgstep, \
CvSize size, double* norm )) \
IPCVAPI_EX( CvStatus, icvNorm_L1_##flavor##_C1R, \
"ippiNorm_L1_" #flavor "_C1R", CV_PLUGINS1(CV_PLUGIN_IPPI), \
( const srctype* img, int imgstep, \
CvSize size, double* norm, CvHintAlgorithm )) \
IPCVAPI_EX( CvStatus, icvNorm_L2_##flavor##_C1R, \
"ippiNorm_L2_" #flavor "_C1R", CV_PLUGINS1(CV_PLUGIN_IPPI), \
( const srctype* img, int imgstep, \
CvSize size, double* norm, CvHintAlgorithm )) \
IPCVAPI_EX( CvStatus, icvNormDiff_Inf_##flavor##_C1R, \
"ippiNormDiff_Inf_" #flavor "_C1R", CV_PLUGINS1(CV_PLUGIN_IPPI), \
( const srctype* img1, int imgstep1, \
const srctype* img2, int imgstep2, \
CvSize size, double* norm )) \
IPCVAPI_EX( CvStatus, icvNormDiff_L1_##flavor##_C1R, \
"ippiNormDiff_L1_" #flavor "_C1R", CV_PLUGINS1(CV_PLUGIN_IPPI), \
( const srctype* img1, int imgstep1, \
const srctype* img2, int imgstep2, \
CvSize size, double* norm, CvHintAlgorithm )) \
IPCVAPI_EX( CvStatus, icvNormDiff_L2_##flavor##_C1R, \
"ippiNormDiff_L2_" #flavor "_C1R", CV_PLUGINS1(CV_PLUGIN_IPPI), \
( const srctype* img1, int imgstep1, \
const srctype* img2, int imgstep2, \
CvSize size, double* norm, CvHintAlgorithm ))
#define IPCV_DEF_NORM_MASK_C1( flavor, srctype ) \
IPCVAPI_EX( CvStatus, icvNorm_Inf_##flavor##_C1MR, \
"ippiNorm_Inf_" #flavor "_C1MR", CV_PLUGINS1(CV_PLUGIN_IPPCV), \
( const srctype* img, int imgstep, \
const uchar* mask, int maskstep, \
CvSize size, double* norm )) \
IPCVAPI_EX( CvStatus, icvNorm_L1_##flavor##_C1MR, \
"ippiNorm_L1_" #flavor "_C1MR", CV_PLUGINS1(CV_PLUGIN_IPPCV), \
( const srctype* img, int imgstep, \
const uchar* mask, int maskstep, \
CvSize size, double* norm )) \
IPCVAPI_EX( CvStatus, icvNorm_L2_##flavor##_C1MR, \
"ippiNorm_L2_" #flavor "_C1MR", CV_PLUGINS1(CV_PLUGIN_IPPCV), \
( const srctype* img, int imgstep, \
const uchar* mask, int maskstep, \
CvSize size, double* norm )) \
IPCVAPI_EX( CvStatus, icvNormDiff_Inf_##flavor##_C1MR, \
"ippiNormDiff_Inf_" #flavor "_C1MR", CV_PLUGINS1(CV_PLUGIN_IPPCV), \
( const srctype* img1, int imgstep1, \
const srctype* img2, int imgstep2, \
const uchar* mask, int maskstep, \
CvSize size, double* norm )) \
IPCVAPI_EX( CvStatus, icvNormDiff_L1_##flavor##_C1MR, \
"ippiNormDiff_L1_" #flavor "_C1MR", CV_PLUGINS1(CV_PLUGIN_IPPCV), \
( const srctype* img1, int imgstep1, \
const srctype* img2, int imgstep2, \
const uchar* mask, int maskstep, \
CvSize size, double* norm )) \
IPCVAPI_EX( CvStatus, icvNormDiff_L2_##flavor##_C1MR, \
"ippiNormDiff_L2_" #flavor "_C1MR", CV_PLUGINS1(CV_PLUGIN_IPPCV), \
( const srctype* img1, int imgstep1, \
const srctype* img2, int imgstep2, \
const uchar* mask, int maskstep, \
CvSize size, double* norm ))
IPCV_DEF_NORM_NOHINT_C1( 8u, uchar )
IPCV_DEF_NORM_MASK_C1( 8u, uchar )
IPCV_DEF_NORM_NOHINT_C1( 16u, ushort )
IPCV_DEF_NORM_MASK_C1( 16u, ushort )
IPCV_DEF_NORM_NOHINT_C1( 16s, short )
IPCV_DEF_NORM_MASK_C1( 16s, short )
IPCV_DEF_NORM_NOHINT_C1( 32s, int )
IPCV_DEF_NORM_MASK_C1( 32s, int )
IPCV_DEF_NORM_HINT_C1( 32f, float )
IPCV_DEF_NORM_MASK_C1( 32f, float )
IPCV_DEF_NORM_NOHINT_C1( 64f, double )
IPCV_DEF_NORM_MASK_C1( 64f, double )
#undef IPCV_DEF_NORM_HONINT_C1
#undef IPCV_DEF_NORM_HINT_C1
#undef IPCV_DEF_NORM_MASK_C1
////////////////////////////////////// AbsDiff ///////////////////////////////////////////
#define IPCV_ABS_DIFF( flavor, arrtype ) \
IPCVAPI_EX( CvStatus, icvAbsDiff_##flavor##_C1R, \
"ippiAbsDiff_" #flavor "_C1R", CV_PLUGINS1(CV_PLUGIN_IPPCV),\
( const arrtype* src1, int srcstep1, \
const arrtype* src2, int srcstep2, \
arrtype* dst, int dststep, CvSize size ))
IPCV_ABS_DIFF( 8u, uchar )
IPCV_ABS_DIFF( 16u, ushort )
IPCV_ABS_DIFF( 16s, short )
IPCV_ABS_DIFF( 32s, int )
IPCV_ABS_DIFF( 32f, float )
IPCV_ABS_DIFF( 64f, double )
#undef IPCV_ABS_DIFF
////////////////////////////// Comparisons //////////////////////////////////////////
#define IPCV_CMP( arrtype, flavor ) \
IPCVAPI_EX( CvStatus, icvCompare_##flavor##_C1R, \
"ippiCompare_" #flavor "_C1R", CV_PLUGINS1(CV_PLUGIN_IPPI), \
( const arrtype* src1, int srcstep1, const arrtype* src2, int srcstep2, \
arrtype* dst, int dststep, CvSize size, int cmp_op )) \
IPCVAPI_EX( CvStatus, icvCompareC_##flavor##_C1R, \
"ippiCompareC_" #flavor "_C1R", CV_PLUGINS1(CV_PLUGIN_IPPI), \
( const arrtype* src1, int srcstep1, arrtype scalar, \
arrtype* dst, int dststep, CvSize size, int cmp_op )) \
IPCVAPI_EX( CvStatus, icvThreshold_GT_##flavor##_C1R, \
"ippiThreshold_GT_" #flavor "_C1R", CV_PLUGINS1(CV_PLUGIN_IPPI), \
( const arrtype* pSrc, int srcstep, arrtype* pDst, int dststep, \
CvSize size, arrtype threshold )) \
IPCVAPI_EX( CvStatus, icvThreshold_LT_##flavor##_C1R, \
"ippiThreshold_LT_" #flavor "_C1R", CV_PLUGINS1(CV_PLUGIN_IPPI), \
( const arrtype* pSrc, int srcstep, arrtype* pDst, int dststep, \
CvSize size, arrtype threshold ))
IPCV_CMP( uchar, 8u )
IPCV_CMP( short, 16s )
IPCV_CMP( float, 32f )
#undef IPCV_CMP
/****************************************************************************************\
* Utilities *
\****************************************************************************************/
////////////////////////////// Copy Pixel <-> Plane /////////////////////////////////
#define IPCV_PIX_PLANE( flavor, arrtype ) \
IPCVAPI_EX( CvStatus, icvCopy_##flavor##_C2P2R, \
"ippiCopy_" #flavor "_C2P2R", 0/*CV_PLUGINS1(CV_PLUGIN_IPPI)*/, \
( const arrtype* src, int srcstep, arrtype** dst, int dststep, CvSize size )) \
IPCVAPI_EX( CvStatus, icvCopy_##flavor##_C3P3R, \
"ippiCopy_" #flavor "_C3P3R", CV_PLUGINS1(CV_PLUGIN_IPPI), \
( const arrtype* src, int srcstep, arrtype** dst, int dststep, CvSize size )) \
IPCVAPI_EX( CvStatus, icvCopy_##flavor##_C4P4R, \
"ippiCopy_" #flavor "_C4P4R", CV_PLUGINS1(CV_PLUGIN_IPPI), \
( const arrtype* src, int srcstep, arrtype** dst, int dststep, CvSize size )) \
IPCVAPI_EX( CvStatus, icvCopy_##flavor##_CnC1CR, \
"ippiCopy_" #flavor "_CnC1CR", 0/*CV_PLUGINS1(CV_PLUGIN_OPTCV)*/, \
( const arrtype* src, int srcstep, arrtype* dst, int dststep, \
CvSize size, int cn, int coi )) \
IPCVAPI_EX( CvStatus, icvCopy_##flavor##_C1CnCR, \
"ippiCopy_" #flavor "_CnC1CR", 0/*CV_PLUGINS1(CV_PLUGIN_OPTCV)*/, \
( const arrtype* src, int srcstep, arrtype* dst, int dststep, \
CvSize size, int cn, int coi )) \
IPCVAPI_EX( CvStatus, icvCopy_##flavor##_P2C2R, \
"ippiCopy_" #flavor "_P2C2R", 0/*CV_PLUGINS1(CV_PLUGIN_IPPI)*/, \
( const arrtype** src, int srcstep, arrtype* dst, int dststep, CvSize size )) \
IPCVAPI_EX( CvStatus, icvCopy_##flavor##_P3C3R, \
"ippiCopy_" #flavor "_P3C3R", CV_PLUGINS1(CV_PLUGIN_IPPI), \
( const arrtype** src, int srcstep, arrtype* dst, int dststep, CvSize size )) \
IPCVAPI_EX( CvStatus, icvCopy_##flavor##_P4C4R, \
"ippiCopy_" #flavor "_P4C4R", CV_PLUGINS1(CV_PLUGIN_IPPI), \
( const arrtype** src, int srcstep, arrtype* dst, int dststep, CvSize size ))
IPCV_PIX_PLANE( 8u, uchar )
IPCV_PIX_PLANE( 16s, ushort )
IPCV_PIX_PLANE( 32f, int )
IPCV_PIX_PLANE( 64f, int64 )
#undef IPCV_PIX_PLANE
/****************************************************************************************/
/* Math routines and RNGs */
/****************************************************************************************/
IPCVAPI_EX( CvStatus, icvInvSqrt_32f, "ippsInvSqrt_32f_A21",
CV_PLUGINS1(CV_PLUGIN_IPPVM),
( const float* src, float* dst, int len ))
IPCVAPI_EX( CvStatus, icvSqrt_32f, "ippsSqrt_32f_A21, ippsSqrt_32f",
CV_PLUGINS2(CV_PLUGIN_IPPVM,CV_PLUGIN_IPPS),
( const float* src, float* dst, int len ))
IPCVAPI_EX( CvStatus, icvInvSqrt_64f, "ippsInvSqrt_64f_A50",
CV_PLUGINS1(CV_PLUGIN_IPPVM),
( const double* src, double* dst, int len ))
IPCVAPI_EX( CvStatus, icvSqrt_64f, "ippsSqrt_64f_A50, ippsSqrt_64f",
CV_PLUGINS2(CV_PLUGIN_IPPVM,CV_PLUGIN_IPPS),
( const double* src, double* dst, int len ))
IPCVAPI_EX( CvStatus, icvLog_32f, "ippsLn_32f_A21, ippsLn_32f",
CV_PLUGINS2(CV_PLUGIN_IPPVM,CV_PLUGIN_IPPS),
( const float *x, float *y, int n ) )
IPCVAPI_EX( CvStatus, icvLog_64f, "ippsLn_64f_A50, ippsLn_64f",
CV_PLUGINS2(CV_PLUGIN_IPPVM,CV_PLUGIN_IPPS),
( const double *x, double *y, int n ) )
IPCVAPI_EX( CvStatus, icvExp_32f, "ippsExp_32f_A21, ippsExp_32f",
CV_PLUGINS2(CV_PLUGIN_IPPVM,CV_PLUGIN_IPPS),
( const float *x, float *y, int n ) )
IPCVAPI_EX( CvStatus, icvExp_64f, "ippsExp_64f_A50, ippsExp_64f",
CV_PLUGINS2(CV_PLUGIN_IPPVM,CV_PLUGIN_IPPS),
( const double *x, double *y, int n ) )
IPCVAPI_EX( CvStatus, icvFastArctan_32f, "ippibFastArctan_32f",
CV_PLUGINS1(CV_PLUGIN_IPPCV),
( const float* y, const float* x, float* angle, int len ))
/****************************************************************************************/
/* Error handling functions */
/****************************************************************************************/
IPCVAPI_EX( CvStatus, icvCheckArray_32f_C1R,
"ippiCheckArray_32f_C1R", 0/*CV_PLUGINS1(CV_PLUGIN_OPTCV)*/,
( const float* src, int srcstep,
CvSize size, int flags,
double min_val, double max_val ))
IPCVAPI_EX( CvStatus, icvCheckArray_64f_C1R,
"ippiCheckArray_64f_C1R", 0/*CV_PLUGINS1(CV_PLUGIN_OPTCV)*/,
( const double* src, int srcstep,
CvSize size, int flags,
double min_val, double max_val ))
/****************************************************************************************/
/* Affine transformations of matrix/image elements */
/****************************************************************************************/
#define IPCV_TRANSFORM( suffix, ipp_suffix, cn ) \
IPCVAPI_EX( CvStatus, icvColorTwist##suffix##_C##cn##R, \
"ippiColorTwist" #ipp_suffix "_C" #cn \
"R,ippiColorTwist" #ipp_suffix "_C" #cn "R", \
CV_PLUGINS2(CV_PLUGIN_IPPI, CV_PLUGIN_IPPCC), \
( const void* src, int srcstep, void* dst, int dststep, \
CvSize roisize, const float* twist_matrix ))
IPCV_TRANSFORM( _8u, 32f_8u, 3 )
IPCV_TRANSFORM( _16u, 32f_16u, 3 )
IPCV_TRANSFORM( _16s, 32f_16s, 3 )
IPCV_TRANSFORM( _32f, _32f, 3 )
IPCV_TRANSFORM( _32f, _32f, 4 )
#undef IPCV_TRANSFORM
#define IPCV_TRANSFORM_N1( suffix ) \
IPCVAPI_EX( CvStatus, icvColorToGray##suffix, \
"ippiColorToGray" #suffix ",ippiColorToGray" #suffix, \
CV_PLUGINS2(CV_PLUGIN_IPPI,CV_PLUGIN_IPPCC), \
( const void* src, int srcstep, void* dst, int dststep, \
CvSize roisize, const float* coeffs ))
IPCV_TRANSFORM_N1( _8u_C3C1R )
IPCV_TRANSFORM_N1( _16u_C3C1R )
IPCV_TRANSFORM_N1( _16s_C3C1R )
IPCV_TRANSFORM_N1( _32f_C3C1R )
IPCV_TRANSFORM_N1( _8u_AC4C1R )
IPCV_TRANSFORM_N1( _16u_AC4C1R )
IPCV_TRANSFORM_N1( _16s_AC4C1R )
IPCV_TRANSFORM_N1( _32f_AC4C1R )
#undef IPCV_TRANSFORM_N1
/****************************************************************************************/
/* Matrix routines from BLAS/LAPACK compatible libraries */
/****************************************************************************************/
IPCVAPI_C_EX( void, icvBLAS_GEMM_32f, "sgemm, mkl_sgemm", CV_PLUGINS2(CV_PLUGIN_MKL,CV_PLUGIN_MKL),
(const char *transa, const char *transb, int *n, int *m, int *k,
const void *alpha, const void *a, int *lda, const void *b, int *ldb,
const void *beta, void *c, int *ldc ))
IPCVAPI_C_EX( void, icvBLAS_GEMM_64f, "dgemm, mkl_dgemm", CV_PLUGINS2(CV_PLUGIN_MKL,CV_PLUGIN_MKL),
(const char *transa, const char *transb, int *n, int *m, int *k,
const void *alpha, const void *a, int *lda, const void *b, int *ldb,
const void *beta, void *c, int *ldc ))
IPCVAPI_C_EX( void, icvBLAS_GEMM_32fc, "cgemm, mkl_cgemm", CV_PLUGINS2(CV_PLUGIN_MKL,CV_PLUGIN_MKL),
(const char *transa, const char *transb, int *n, int *m, int *k,
const void *alpha, const void *a, int *lda, const void *b, int *ldb,
const void *beta, void *c, int *ldc ))
IPCVAPI_C_EX( void, icvBLAS_GEMM_64fc, "zgemm, mkl_zgemm", CV_PLUGINS2(CV_PLUGIN_MKL,CV_PLUGIN_MKL),
(const char *transa, const char *transb, int *n, int *m, int *k,
const void *alpha, const void *a, int *lda, const void *b, int *ldb,
const void *beta, void *c, int *ldc ))
#define IPCV_DFT( init_flavor, fwd_flavor, inv_flavor ) \
IPCVAPI_EX( CvStatus, icvDFTInitAlloc_##init_flavor, "ippsDFTInitAlloc_" #init_flavor, \
CV_PLUGINS1(CV_PLUGIN_IPPS), ( void**, int, int, CvHintAlgorithm )) \
\
IPCVAPI_EX( CvStatus, icvDFTFree_##init_flavor, "ippsDFTFree_" #init_flavor, \
CV_PLUGINS1(CV_PLUGIN_IPPS), ( void* )) \
\
IPCVAPI_EX( CvStatus, icvDFTGetBufSize_##init_flavor, "ippsDFTGetBufSize_" #init_flavor,\
CV_PLUGINS1(CV_PLUGIN_IPPS), ( const void* spec, int* buf_size )) \
\
IPCVAPI_EX( CvStatus, icvDFTFwd_##fwd_flavor, "ippsDFTFwd_" #fwd_flavor, \
CV_PLUGINS1(CV_PLUGIN_IPPS), ( const void* src, void* dst, \
const void* spec, void* buffer )) \
\
IPCVAPI_EX( CvStatus, icvDFTInv_##inv_flavor, "ippsDFTInv_" #inv_flavor, \
CV_PLUGINS1(CV_PLUGIN_IPPS), ( const void* src, void* dst, \
const void* spec, void* buffer ))
IPCV_DFT( C_32fc, CToC_32fc, CToC_32fc )
IPCV_DFT( R_32f, RToPack_32f, PackToR_32f )
IPCV_DFT( C_64fc, CToC_64fc, CToC_64fc )
IPCV_DFT( R_64f, RToPack_64f, PackToR_64f )
#undef IPCV_DFT
#endif /*_CXCORE_IPP_H_*/
测试代码#include <STDIO.H>
#include <iostream>
#include "_cxcore_my.h"
#pragma comment(lib, "cxcore.lib")
/////////////////////////////////////////////////////////////////////////////////////////
#define icvGivens_64f( n, x, y, c, s ) \
{ \
int _i; \
double* _x = (x); \
double* _y = (y); \
\
for( _i = 0; _i < n; _i++ ) \
{ \
double t0 = _x[_i]; \
double t1 = _y[_i]; \
_x[_i] = t0*c + t1*s; \
_y[_i] = -t0*s + t1*c; \
} \
}
/* y[0:m,0:n] += diag(a[0:1,0:m]) * x[0:m,0:n] */
static void
icvMatrAXPY_64f( int m, int n, const double* x, int dx,
const double* a, double* y, int dy )
{
int i, j;
for( i = 0; i < m; i++, x += dx, y += dy )
{
double s = a[i];
for( j = 0; j <= n - 4; j += 4 )
{
double t0 = y[j] + s*x[j];
double t1 = y[j+1] + s*x[j+1];
y[j] = t0;
y[j+1] = t1;
t0 = y[j+2] + s*x[j+2];
t1 = y[j+3] + s*x[j+3];
y[j+2] = t0;
y[j+3] = t1;
}
for( ; j < n; j++ ) y[j] += s*x[j];
}
}
/* y[1:m,-1] = h*y[1:m,0:n]*x[0:1,0:n]'*x[-1] (this is used for U&V reconstruction)
y[1:m,0:n] += h*y[1:m,0:n]*x[0:1,0:n]'*x[0:1,0:n] */
static void
icvMatrAXPY3_64f( int m, int n, const double* x, int l, double* y, double h )
{
int i, j;
for( i = 1; i < m; i++ )
{
double s = 0;
y += l;
for( j = 0; j <= n - 4; j += 4 )
s += x[j]*y[j] + x[j+1]*y[j+1] + x[j+2]*y[j+2] + x[j+3]*y[j+3];
for( ; j < n; j++ ) s += x[j]*y[j];
s *= h;
y[-1] = s*x[-1];
for( j = 0; j <= n - 4; j += 4 )
{
double t0 = y[j] + s*x[j];
double t1 = y[j+1] + s*x[j+1];
y[j] = t0;
y[j+1] = t1;
t0 = y[j+2] + s*x[j+2];
t1 = y[j+3] + s*x[j+3];
y[j+2] = t0;
y[j+3] = t1;
}
for( ; j < n; j++ ) y[j] += s*x[j];
}
}
#define icvGivens_32f( n, x, y, c, s ) \
{ \
int _i; \
float* _x = (x); \
float* _y = (y); \
\
for( _i = 0; _i < n; _i++ ) \
{ \
double t0 = _x[_i]; \
double t1 = _y[_i]; \
_x[_i] = (float)(t0*c + t1*s); \
_y[_i] = (float)(-t0*s + t1*c);\
} \
}
static void
icvMatrAXPY_32f( int m, int n, const float* x, int dx,
const float* a, float* y, int dy )
{
int i, j;
for( i = 0; i < m; i++, x += dx, y += dy )
{
double s = a[i];
for( j = 0; j <= n - 4; j += 4 )
{
double t0 = y[j] + s*x[j];
double t1 = y[j+1] + s*x[j+1];
y[j] = (float)t0;
y[j+1] = (float)t1;
t0 = y[j+2] + s*x[j+2];
t1 = y[j+3] + s*x[j+3];
y[j+2] = (float)t0;
y[j+3] = (float)t1;
}
for( ; j < n; j++ )
y[j] = (float)(y[j] + s*x[j]);
}
}
static void
icvMatrAXPY3_32f( int m, int n, const float* x, int l, float* y, double h )
{
int i, j;
for( i = 1; i < m; i++ )
{
double s = 0;
y += l;
for( j = 0; j <= n - 4; j += 4 )
s += x[j]*y[j] + x[j+1]*y[j+1] + x[j+2]*y[j+2] + x[j+3]*y[j+3];
for( ; j < n; j++ ) s += x[j]*y[j];
s *= h;
y[-1] = (float)(s*x[-1]);
for( j = 0; j <= n - 4; j += 4 )
{
double t0 = y[j] + s*x[j];
double t1 = y[j+1] + s*x[j+1];
y[j] = (float)t0;
y[j+1] = (float)t1;
t0 = y[j+2] + s*x[j+2];
t1 = y[j+3] + s*x[j+3];
y[j+2] = (float)t0;
y[j+3] = (float)t1;
}
for( ; j < n; j++ ) y[j] = (float)(y[j] + s*x[j]);
}
}
/* accurate hypotenuse calculation */
static double
pythag( double a, double b )
{
a = fabs( a );
b = fabs( b );
if( a > b )
{
b /= a;
a *= sqrt( 1. + b * b );
}
else if( b != 0 )
{
a /= b;
a = b * sqrt( 1. + a * a );
}
return a;
}
/****************************************************************************************/
/****************************************************************************************/
#define MAX_ITERS 30
static void icvSVD_32f( float* a, int lda, int m, int n,
float* w,
float* uT, int lduT, int nu,
float* vT, int ldvT,
float* buffer )
{
float* e;
float* temp;
float *w1, *e1;
float *hv;
double ku0 = 0, kv0 = 0;
double anorm = 0;
float *a1, *u0 = uT, *v0 = vT;
double scale, h;
int i, j, k, l;
int nm, m1, n1;
int nv = n;
int iters = 0;
float* hv0 = (float*)cvStackAlloc( (m+2)*sizeof(hv0[0])) + 1;
e = buffer;
w1 = w;
e1 = e + 1;
nm = n;
temp = buffer + nm;
memset( w, 0, nm * sizeof( w[0] ));
memset( e, 0, nm * sizeof( e[0] ));
m1 = m;
n1 = n;
/* transform a to bi-diagonal form */
for( ;; )
{
int update_u;
int update_v;
if( m1 == 0 )
break;
scale = h = 0;
update_u = uT && m1 > m - nu;
hv = update_u ? uT : hv0;
for( j = 0, a1 = a; j < m1; j++, a1 += lda )
{
double t = a1[0];
scale += fabs( hv[j] = (float)t );
}
if( scale != 0 )
{
double f = 1./scale, g, s = 0;
for( j = 0; j < m1; j++ )
{
double t = (hv[j] = (float)(hv[j]*f));
s += t * t;
}
g = sqrt( s );
f = hv[0];
if( f >= 0 )
g = -g;
hv[0] = (float)(f - g);
h = 1. / (f * g - s);
memset( temp, 0, n1 * sizeof( temp[0] ));
/* calc temp[0:n-i] = a[i:m,i:n]'*hv[0:m-i] */
icvMatrAXPY_32f( m1, n1 - 1, a + 1, lda, hv, temp + 1, 0 );
for( k = 1; k < n1; k++ ) temp[k] = (float)(temp[k]*h);
/* modify a: a[i:m,i:n] = a[i:m,i:n] + hv[0:m-i]*temp[0:n-i]' */
icvMatrAXPY_32f( m1, n1 - 1, temp + 1, 0, hv, a + 1, lda );
*w1 = (float)(g*scale);
}
w1++;
/* store -2/(hv'*hv) */
if( update_u )
{
if( m1 == m )
ku0 = h;
else
hv[-1] = (float)h;
}
a++;
n1--;
if( vT )
vT += ldvT + 1;
if( n1 == 0 )
break;
scale = h = 0;
update_v = vT && n1 > n - nv;
hv = update_v ? vT : hv0;
for( j = 0; j < n1; j++ )
{
double t = a[j];
scale += fabs( hv[j] = (float)t );
}
if( scale != 0 )
{
double f = 1./scale, g, s = 0;
for( j = 0; j < n1; j++ )
{
double t = (hv[j] = (float)(hv[j]*f));
s += t * t;
}
g = sqrt( s );
f = hv[0];
if( f >= 0 )
g = -g;
hv[0] = (float)(f - g);
h = 1. / (f * g - s);
hv[-1] = 0.f;
/* update a[i:m:i+1:n] = a[i:m,i+1:n] + (a[i:m,i+1:n]*hv[0:m-i])*... */
icvMatrAXPY3_32f( m1, n1, hv, lda, a, h );
*e1 = (float)(g*scale);
}
e1++;
/* store -2/(hv'*hv) */
if( update_v )
{
if( n1 == n )
kv0 = h;
else
hv[-1] = (float)h;
}
a += lda;
m1--;
if( uT )
uT += lduT + 1;
}
m1 -= m1 != 0;
n1 -= n1 != 0;
/* accumulate left transformations */
if( uT )
{
m1 = m - m1;
uT = u0 + m1 * lduT;
for( i = m1; i < nu; i++, uT += lduT )
{
memset( uT + m1, 0, (m - m1) * sizeof( uT[0] ));
uT[i] = 1.;
}
for( i = m1 - 1; i >= 0; i-- )
{
double s;
int lh = nu - i;
l = m - i;
hv = u0 + (lduT + 1) * i;
h = i == 0 ? ku0 : hv[-1];
assert( h <= 0 );
if( h != 0 )
{
uT = hv;
icvMatrAXPY3_32f( lh, l-1, hv+1, lduT, uT+1, h );
s = hv[0] * h;
for( k = 0; k < l; k++ ) hv[k] = (float)(hv[k]*s);
hv[0] += 1;
}
else
{
for( j = 1; j < l; j++ )
hv[j] = 0;
for( j = 1; j < lh; j++ )
hv[j * lduT] = 0;
hv[0] = 1;
}
}
uT = u0;
}
/* accumulate right transformations */
if( vT )
{
n1 = n - n1;
vT = v0 + n1 * ldvT;
for( i = n1; i < nv; i++, vT += ldvT )
{
memset( vT + n1, 0, (n - n1) * sizeof( vT[0] ));
vT[i] = 1.;
}
for( i = n1 - 1; i >= 0; i-- )
{
double s;
int lh = nv - i;
l = n - i;
hv = v0 + (ldvT + 1) * i;
h = i == 0 ? kv0 : hv[-1];
assert( h <= 0 );
if( h != 0 )
{
vT = hv;
icvMatrAXPY3_32f( lh, l-1, hv+1, ldvT, vT+1, h );
s = hv[0] * h;
for( k = 0; k < l; k++ ) hv[k] = (float)(hv[k]*s);
hv[0] += 1;
}
else
{
for( j = 1; j < l; j++ )
hv[j] = 0;
for( j = 1; j < lh; j++ )
hv[j * ldvT] = 0;
hv[0] = 1;
}
}
vT = v0;
}
for( i = 0; i < nm; i++ )
{
double tnorm = fabs( w[i] );
tnorm += fabs( e[i] );
if( anorm < tnorm )
anorm = tnorm;
}
anorm *= FLT_EPSILON;
/* diagonalization of the bidiagonal form */
for( k = nm - 1; k >= 0; k-- )
{
double z = 0;
iters = 0;
for( ;; ) /* do iterations */
{
double c, s, f, g, x, y;
int flag = 0;
/* test for splitting */
for( l = k; l >= 0; l-- )
{
if( fabs( e[l] ) <= anorm )
{
flag = 1;
break;
}
assert( l > 0 );
if( fabs( w[l - 1] ) <= anorm )
break;
}
if( !flag )
{
c = 0;
s = 1;
for( i = l; i <= k; i++ )
{
f = s * e[i];
e[i] = (float)(e[i]*c);
if( anorm + fabs( f ) == anorm )
break;
g = w[i];
h = pythag( f, g );
w[i] = (float)h;
c = g / h;
s = -f / h;
if( uT )
icvGivens_32f( m, uT + lduT * (l - 1), uT + lduT * i, c, s );
}
}
z = w[k];
if( l == k || iters++ == MAX_ITERS )
break;
/* shift from bottom 2x2 minor */
x = w[l];
y = w[k - 1];
g = e[k - 1];
h = e[k];
f = 0.5 * (((g + z) / h) * ((g - z) / y) + y / h - h / y);
g = pythag( f, 1 );
if( f < 0 )
g = -g;
f = x - (z / x) * z + (h / x) * (y / (f + g) - h);
/* next QR transformation */
c = s = 1;
for( i = l + 1; i <= k; i++ )
{
g = e[i];
y = w[i];
h = s * g;
g *= c;
z = pythag( f, h );
e[i - 1] = (float)z;
c = f / z;
s = h / z;
f = x * c + g * s;
g = -x * s + g * c;
h = y * s;
y *= c;
if( vT )
icvGivens_32f( n, vT + ldvT * (i - 1), vT + ldvT * i, c, s );
z = pythag( f, h );
w[i - 1] = (float)z;
/* rotation can be arbitrary if z == 0 */
if( z != 0 )
{
c = f / z;
s = h / z;
}
f = c * g + s * y;
x = -s * g + c * y;
if( uT )
icvGivens_32f( m, uT + lduT * (i - 1), uT + lduT * i, c, s );
}
e[l] = 0;
e[k] = (float)f;
w[k] = (float)x;
} /* end of iteration loop */
if( iters > MAX_ITERS )
break;
if( z < 0 )
{
w[k] = (float)(-z);
if( vT )
{
for( j = 0; j < n; j++ )
vT[j + k * ldvT] = -vT[j + k * ldvT];
}
}
} /* end of diagonalization loop */
/* sort singular values and corresponding vectors */
for( i = 0; i < nm; i++ )
{
k = i;
for( j = i + 1; j < nm; j++ )
if( w[k] < w[j] )
k = j;
if( k != i )
{
float t;
CV_SWAP( w[i], w[k], t );
if( vT )
for( j = 0; j < n; j++ )
CV_SWAP( vT[j + ldvT*k], vT[j + ldvT*i], t );
if( uT )
for( j = 0; j < m; j++ )
CV_SWAP( uT[j + lduT*k], uT[j + lduT*i], t );
}
}
}
static void icvSVD_64f( double* a, int lda, int m, int n,
double* w,
double* uT, int lduT, int nu,
double* vT, int ldvT,
double* buffer )
{
double* e;
double* temp;
double *w1, *e1;
double *hv;
double ku0 = 0, kv0 = 0;
double anorm = 0;
double *a1, *u0 = uT, *v0 = vT;
double scale, h;
int i, j, k, l;
int nm, m1, n1;
int nv = n;
int iters = 0;
double* hv0 = (double*)cvStackAlloc( (m+2)*sizeof(hv0[0])) + 1;
e = buffer;
w1 = w;
e1 = e + 1;
nm = n;
temp = buffer + nm;
memset( w, 0, nm * sizeof( w[0] ));
memset( e, 0, nm * sizeof( e[0] ));
m1 = m;
n1 = n;
/* transform a to bi-diagonal form */
for( ;; )
{
int update_u;
int update_v;
if( m1 == 0 )
break;
scale = h = 0;
update_u = uT && m1 > m - nu;
hv = update_u ? uT : hv0;
for( j = 0, a1 = a; j < m1; j++, a1 += lda )
{
double t = a1[0];
scale += fabs( hv[j] = t );
}
if( scale != 0 )
{
double f = 1./scale, g, s = 0;
for( j = 0; j < m1; j++ )
{
double t = (hv[j] *= f);
s += t * t;
}
g = sqrt( s );
f = hv[0];
if( f >= 0 )
g = -g;
hv[0] = f - g;
h = 1. / (f * g - s);
memset( temp, 0, n1 * sizeof( temp[0] ));
/* calc temp[0:n-i] = a[i:m,i:n]'*hv[0:m-i] */
icvMatrAXPY_64f( m1, n1 - 1, a + 1, lda, hv, temp + 1, 0 );
for( k = 1; k < n1; k++ ) temp[k] *= h;
/* modify a: a[i:m,i:n] = a[i:m,i:n] + hv[0:m-i]*temp[0:n-i]' */
icvMatrAXPY_64f( m1, n1 - 1, temp + 1, 0, hv, a + 1, lda );
*w1 = g*scale;
}
w1++;
/* store -2/(hv'*hv) */
if( update_u )
{
if( m1 == m )
ku0 = h;
else
hv[-1] = h;
}
a++;
n1--;
if( vT )
vT += ldvT + 1;
if( n1 == 0 )
break;
scale = h = 0;
update_v = vT && n1 > n - nv;
hv = update_v ? vT : hv0;
for( j = 0; j < n1; j++ )
{
double t = a[j];
scale += fabs( hv[j] = t );
}
if( scale != 0 )
{
double f = 1./scale, g, s = 0;
for( j = 0; j < n1; j++ )
{
double t = (hv[j] *= f);
s += t * t;
}
g = sqrt( s );
f = hv[0];
if( f >= 0 )
g = -g;
hv[0] = f - g;
h = 1. / (f * g - s);
hv[-1] = 0.;
/* update a[i:m:i+1:n] = a[i:m,i+1:n] + (a[i:m,i+1:n]*hv[0:m-i])*... */
icvMatrAXPY3_64f( m1, n1, hv, lda, a, h );
*e1 = g*scale;
}
e1++;
/* store -2/(hv'*hv) */
if( update_v )
{
if( n1 == n )
kv0 = h;
else
hv[-1] = h;
}
a += lda;
m1--;
if( uT )
uT += lduT + 1;
}
m1 -= m1 != 0;
n1 -= n1 != 0;
/* accumulate left transformations */
if( uT )
{
m1 = m - m1;
uT = u0 + m1 * lduT;
for( i = m1; i < nu; i++, uT += lduT )
{
memset( uT + m1, 0, (m - m1) * sizeof( uT[0] ));
uT[i] = 1.;
}
for( i = m1 - 1; i >= 0; i-- )
{
double s;
int lh = nu - i;
l = m - i;
hv = u0 + (lduT + 1) * i;
h = i == 0 ? ku0 : hv[-1];
assert( h <= 0 );
if( h != 0 )
{
uT = hv;
icvMatrAXPY3_64f( lh, l-1, hv+1, lduT, uT+1, h );
s = hv[0] * h;
for( k = 0; k < l; k++ ) hv[k] *= s;
hv[0] += 1;
}
else
{
for( j = 1; j < l; j++ )
hv[j] = 0;
for( j = 1; j < lh; j++ )
hv[j * lduT] = 0;
hv[0] = 1;
}
}
uT = u0;
}
/* accumulate right transformations */
if( vT )
{
n1 = n - n1;
vT = v0 + n1 * ldvT;
for( i = n1; i < nv; i++, vT += ldvT )
{
memset( vT + n1, 0, (n - n1) * sizeof( vT[0] ));
vT[i] = 1.;
}
for( i = n1 - 1; i >= 0; i-- )
{
double s;
int lh = nv - i;
l = n - i;
hv = v0 + (ldvT + 1) * i;
h = i == 0 ? kv0 : hv[-1];
assert( h <= 0 );
if( h != 0 )
{
vT = hv;
icvMatrAXPY3_64f( lh, l-1, hv+1, ldvT, vT+1, h );
s = hv[0] * h;
for( k = 0; k < l; k++ ) hv[k] *= s;
hv[0] += 1;
}
else
{
for( j = 1; j < l; j++ )
hv[j] = 0;
for( j = 1; j < lh; j++ )
hv[j * ldvT] = 0;
hv[0] = 1;
}
}
vT = v0;
}
for( i = 0; i < nm; i++ )
{
double tnorm = fabs( w[i] );
tnorm += fabs( e[i] );
if( anorm < tnorm )
anorm = tnorm;
}
anorm *= DBL_EPSILON;
/* diagonalization of the bidiagonal form */
for( k = nm - 1; k >= 0; k-- )
{
double z = 0;
iters = 0;
for( ;; ) /* do iterations */
{
double c, s, f, g, x, y;
int flag = 0;
/* test for splitting */
for( l = k; l >= 0; l-- )
{
if( fabs(e[l]) <= anorm )
{
flag = 1;
break;
}
assert( l > 0 );
if( fabs(w[l - 1]) <= anorm )
break;
}
if( !flag )
{
c = 0;
s = 1;
for( i = l; i <= k; i++ )
{
f = s * e[i];
e[i] *= c;
if( anorm + fabs( f ) == anorm )
break;
g = w[i];
h = pythag( f, g );
w[i] = h;
c = g / h;
s = -f / h;
if( uT )
icvGivens_64f( m, uT + lduT * (l - 1), uT + lduT * i, c, s );
}
}
z = w[k];
if( l == k || iters++ == MAX_ITERS )
break;
/* shift from bottom 2x2 minor */
x = w[l];
y = w[k - 1];
g = e[k - 1];
h = e[k];
f = 0.5 * (((g + z) / h) * ((g - z) / y) + y / h - h / y);
g = pythag( f, 1 );
if( f < 0 )
g = -g;
f = x - (z / x) * z + (h / x) * (y / (f + g) - h);
/* next QR transformation */
c = s = 1;
for( i = l + 1; i <= k; i++ )
{
g = e[i];
y = w[i];
h = s * g;
g *= c;
z = pythag( f, h );
e[i - 1] = z;
c = f / z;
s = h / z;
f = x * c + g * s;
g = -x * s + g * c;
h = y * s;
y *= c;
if( vT )
icvGivens_64f( n, vT + ldvT * (i - 1), vT + ldvT * i, c, s );
z = pythag( f, h );
w[i - 1] = z;
/* rotation can be arbitrary if z == 0 */
if( z != 0 )
{
c = f / z;
s = h / z;
}
f = c * g + s * y;
x = -s * g + c * y;
if( uT )
icvGivens_64f( m, uT + lduT * (i - 1), uT + lduT * i, c, s );
}
e[l] = 0;
e[k] = f;
w[k] = x;
} /* end of iteration loop */
if( iters > MAX_ITERS )
break;
if( z < 0 )
{
w[k] = -z;
if( vT )
{
for( j = 0; j < n; j++ )
vT[j + k * ldvT] = -vT[j + k * ldvT];
}
}
} /* end of diagonalization loop */
/* sort singular values and corresponding values */
for( i = 0; i < nm; i++ )
{
k = i;
for( j = i + 1; j < nm; j++ )
if( w[k] < w[j] )
k = j;
if( k != i )
{
double t;
CV_SWAP( w[i], w[k], t );
if( vT )
for( j = 0; j < n; j++ )
CV_SWAP( vT[j + ldvT*k], vT[j + ldvT*i], t );
if( uT )
for( j = 0; j < m; j++ )
CV_SWAP( uT[j + lduT*k], uT[j + lduT*i], t );
}
}
}
void cvSVD_my( CvArr* aarr, CvArr* warr, CvArr* uarr, CvArr* varr, int flags )
{
uchar* buffer = 0;
int local_alloc = 0;
CV_FUNCNAME( "cvSVD" );
__BEGIN__;
CvMat astub, *a = (CvMat*)aarr;
CvMat wstub, *w = (CvMat*)warr;
CvMat ustub, *u;
CvMat vstub, *v;
CvMat tmat;
uchar* tw = 0;
int type;
int a_buf_offset = 0, u_buf_offset = 0, buf_size, pix_size;
int temp_u = 0, /* temporary storage for U is needed */
t_svd; /* special case: a->rows < a->cols */
int m, n;
int w_rows, w_cols;
int u_rows = 0, u_cols = 0;
int w_is_mat = 0;
if( !CV_IS_MAT( a ))
CV_CALL( a = cvGetMat( a, &astub ));
if( !CV_IS_MAT( w ))
CV_CALL( w = cvGetMat( w, &wstub ));
if( !CV_ARE_TYPES_EQ( a, w ))
CV_ERROR( CV_StsUnmatchedFormats, "" );
if( a->rows >= a->cols )
{
m = a->rows;
n = a->cols;
w_rows = w->rows;
w_cols = w->cols;
t_svd = 0;
}
else
{
CvArr* t;
CV_SWAP( uarr, varr, t );
flags = (flags & CV_SVD_U_T ? CV_SVD_V_T : 0)|
(flags & CV_SVD_V_T ? CV_SVD_U_T : 0);
m = a->cols;
n = a->rows;
w_rows = w->cols;
w_cols = w->rows;
t_svd = 1;
}
u = (CvMat*)uarr;
v = (CvMat*)varr;
w_is_mat = w_cols > 1 && w_rows > 1;
if( !w_is_mat && CV_IS_MAT_CONT(w->type) && w_cols + w_rows - 1 == n )
tw = w->data.ptr;
if( u )
{
if( !CV_IS_MAT( u ))
CV_CALL( u = cvGetMat( u, &ustub ));
if( !(flags & CV_SVD_U_T) )
{
u_rows = u->rows;
u_cols = u->cols;
}
else
{
u_rows = u->cols;
u_cols = u->rows;
}
if( !CV_ARE_TYPES_EQ( a, u ))
CV_ERROR( CV_StsUnmatchedFormats, "" );
if( u_rows != m || (u_cols != m && u_cols != n))
CV_ERROR( CV_StsUnmatchedSizes, !t_svd ? "U matrix has unappropriate size" :
"V matrix has unappropriate size" );
temp_u = (u_rows != u_cols && !(flags & CV_SVD_U_T)) || u->data.ptr==a->data.ptr;
if( w_is_mat && u_cols != w_rows )
CV_ERROR( CV_StsUnmatchedSizes, !t_svd ? "U and W have incompatible sizes" :
"V and W have incompatible sizes" );
}
else
{
u = &ustub;
u->data.ptr = 0;
u->step = 0;
}
if( v )
{
int v_rows, v_cols;
if( !CV_IS_MAT( v ))
CV_CALL( v = cvGetMat( v, &vstub ));
if( !(flags & CV_SVD_V_T) )
{
v_rows = v->rows;
v_cols = v->cols;
}
else
{
v_rows = v->cols;
v_cols = v->rows;
}
if( !CV_ARE_TYPES_EQ( a, v ))
CV_ERROR( CV_StsUnmatchedFormats, "" );
if( v_rows != n || v_cols != n )
CV_ERROR( CV_StsUnmatchedSizes, t_svd ? "U matrix has unappropriate size" :
"V matrix has unappropriate size" );
if( w_is_mat && w_cols != v_cols )
CV_ERROR( CV_StsUnmatchedSizes, t_svd ? "U and W have incompatible sizes" :
"V and W have incompatible sizes" );
}
else
{
v = &vstub;
v->data.ptr = 0;
v->step = 0;
}
type = CV_MAT_TYPE( a->type );
pix_size = CV_ELEM_SIZE(type);
buf_size = n*2 + m;
if( !(flags & CV_SVD_MODIFY_A) )
{
a_buf_offset = buf_size;
buf_size += a->rows*a->cols;
}
if( temp_u )
{
u_buf_offset = buf_size;
buf_size += u->rows*u->cols;
}
buf_size *= pix_size;
if( buf_size <= CV_MAX_LOCAL_SIZE )
{
buffer = (uchar*)cvStackAlloc( buf_size );
local_alloc = 1;
}
else
{
CV_CALL( buffer = (uchar*)cvAlloc( buf_size ));
}
if( !(flags & CV_SVD_MODIFY_A) )
{
cvInitMatHeader( &tmat, m, n, type,
buffer + a_buf_offset*pix_size );
if( !t_svd )
cvCopy( a, &tmat );
else
cvT( a, &tmat );
a = &tmat;
}
if( temp_u )
{
cvInitMatHeader( &ustub, u_cols, u_rows, type, buffer + u_buf_offset*pix_size );
u = &ustub;
}
if( !tw )
tw = buffer + (n + m)*pix_size;
if( type == CV_32FC1 )
{
icvSVD_32f( a->data.fl, a->step/sizeof(float), a->rows, a->cols,
(float*)tw, u->data.fl, u->step/sizeof(float), u_cols,
v->data.fl, v->step/sizeof(float), (float*)buffer );
}
else if( type == CV_64FC1 )
{
icvSVD_64f( a->data.db, a->step/sizeof(double), a->rows, a->cols,
(double*)tw, u->data.db, u->step/sizeof(double), u_cols,
v->data.db, v->step/sizeof(double), (double*)buffer );
}
else
{
CV_ERROR( CV_StsUnsupportedFormat, "" );
}
if( tw != w->data.ptr )
{
int shift = w->cols != 1;
cvSetZero( w );
if( type == CV_32FC1 )
for( int i = 0; i < n; i++ )
((float*)(w->data.ptr + i*w->step))[i*shift] = ((float*)tw)[i];
else
for( int i = 0; i < n; i++ )
((double*)(w->data.ptr + i*w->step))[i*shift] = ((double*)tw)[i];
}
if( uarr )
{
if( !(flags & CV_SVD_U_T))
cvT( u, uarr );
else if( temp_u )
cvCopy( u, uarr );
/*CV_CHECK_NANS( uarr );*/
}
if( varr )
{
if( !(flags & CV_SVD_V_T))
cvT( v, varr );
/*CV_CHECK_NANS( varr );*/
}
CV_CHECK_NANS( w );
__END__;
if( buffer && !local_alloc )
cvFree( &buffer );
}
inline void cvDoubleMatPrint_my( const CvMat* mat )
{
int i, j;
for( i = 0; i < mat->rows; i++ )
{
for( j = 0; j < mat->cols; j++ )
{
printf( "%lf ",cvmGet( mat, i, j ) );
}
printf( "\n" );
}
}
参考代码：#include <iostream>
#include <algorithm>
#include <string>
#include <io.h>
#include "highgui.h"
#include "cv.h"
#pragma comment(lib, "cxcore.lib")
using namespace std;
inline void cvDoubleMatPrint( const CvMat* mat );
void SVD(double *SA, double *SU, double *SS, double *SV);
int main()
{
double A[16]={1,1,1,1,0,1,0,0,0,0,1,1,0,0,0,0}, U[16], S[16], V[16];
SVD(A, U, S, V);
printf("V:============\n");
for(int j=0,i=0;i<16;i++)
{
cout<<V[i]<<',';
j++;
if(j%4==0)
cout<<endl;
}
return 0;
}
inline void cvDoubleMatPrint( const CvMat* mat )
{
int i, j;
for( i = 0; i < mat->rows; i++ )
{
for( j = 0; j < mat->cols; j++ )
{
printf( "%lf ",cvmGet( mat, i, j ) );
}
printf( "\n" );
}
}
void SVD(double *SA, double *SU, double *SS, double *SV)
{
CvMat A, U, S, V;
cvInitMatHeader(&A,4,4,CV_64FC1,SA,CV_AUTOSTEP);
cvInitMatHeader(&U,4,4,CV_64FC1,SU,CV_AUTOSTEP);
cvInitMatHeader(&S,4,4,CV_64FC1,SS,CV_AUTOSTEP);
cvInitMatHeader(&V,4,4,CV_64FC1,SV,CV_AUTOSTEP);
cvSVD(&A, &U, &S, &V, CV_SVD_U_T);
printf("A:============\n");
cvDoubleMatPrint( &A);
printf("U:============\n");
cvDoubleMatPrint( &U);
printf("S:============\n");
cvDoubleMatPrint( &S);
printf("V:============\n");
cvDoubleMatPrint( &V);
}
A =
1 1 1 1
0 1 0 0
0 0 1 1
0 0 0 0
>> [S V D]=svd(A)
S =
0.8460 0.2261 0.4828 0
0.1922 0.7154 -0.6718 0
0.4973 -0.6611 -0.5618 0
0 0 0 1.0000
V =
2.3244 0 0 0
0 1.1472 0 0
0 0 0.5304 0
0 0 0 0
D =
0.3640 0.1971 0.9103 0
0.4467 0.8207 -0.3563 0
0.5779 -0.3792 -0.1490 -0.7071
0.5779 -0.3792 -0.1490 0.7071
*/opencv cvSVD 原理
opencv 奇异值分解（SVD）测试代码：
需要配置opencv1.0和vc6
_cxcore_my.h
#ifndef _CXCORE_INTERNAL_H_
#define _CXCORE_INTERNAL_H_
#if defined _MSC_VER && _MSC_VER >= 1200
/* disable warnings related to inline functions */
#pragma warning( disable: 4711 4710 4514 )
#endif
typedef unsigned long ulong;
#ifdef __BORLANDC__
#define WIN32
#define CV_DLL
#undef _CV_ALWAYS_PROFILE_
#define _CV_ALWAYS_NO_PROFILE_
#endif
#include "cxcore.h"
#include "cxmisc.h"
#include "_cxipp_my.h"
#include <math.h>
#include <assert.h>
#include <string.h>
#include <stdlib.h>
#include <stdio.h>
#include <limits.h>
#include <float.h>
// -128.f ... 255.f
extern const float icv8x32fTab[];
#define CV_8TO32F(x) icv8x32fTab[(x)+128]
extern const ushort icv8x16uSqrTab[];
#define CV_SQR_8U(x) icv8x16uSqrTab[(x)+255]
extern const char* icvHersheyGlyphs[];
extern const signed char icvDepthToType[];
#define icvIplToCvDepth( depth ) \
icvDepthToType[(((depth) & 255) >> 2) + ((depth) < 0)]
extern const uchar icvSaturate8u[];
#define CV_FAST_CAST_8U(t) (assert(-256 <= (t) && (t) <= 512), icvSaturate8u[(t)+256])
#define CV_MIN_8U(a,b) ((a) - CV_FAST_CAST_8U((a) - (b)))
#define CV_MAX_8U(a,b) ((a) + CV_FAST_CAST_8U((b) - (a)))
typedef CvFunc2D_3A1I CvArithmBinMaskFunc2D;
typedef CvFunc2D_2A1P1I CvArithmUniMaskFunc2D;
/****************************************************************************************\
* Complex arithmetics *
\****************************************************************************************/
struct CvComplex32f;
struct CvComplex64f;
struct CvComplex32f
{
float re, im;
CvComplex32f() {}
CvComplex32f( float _re, float _im=0 ) : re(_re), im(_im) {}
explicit CvComplex32f( const CvComplex64f& v );
//CvComplex32f( const CvComplex32f& v ) : re(v.re), im(v.im) {}
//CvComplex32f& operator = (const CvComplex32f& v ) { re = v.re; im = v.im; return *this; }
operator CvComplex64f() const;
};
struct CvComplex64f
{
double re, im;
CvComplex64f() {}
CvComplex64f( double _re, double _im=0 ) : re(_re), im(_im) {}
explicit CvComplex64f( const CvComplex32f& v );
//CvComplex64f( const CvComplex64f& v ) : re(v.re), im(v.im) {}
//CvComplex64f& operator = (const CvComplex64f& v ) { re = v.re; im = v.im; return *this; }
operator CvComplex32f() const;
};
inline CvComplex32f::CvComplex32f( const CvComplex64f& v ) : re((float)v.re), im((float)v.im) {}
inline CvComplex64f::CvComplex64f( const CvComplex32f& v ) : re(v.re), im(v.im) {}
inline CvComplex32f operator + (CvComplex32f a, CvComplex32f b)
{
return CvComplex32f( a.re + b.re, a.im + b.im );
}
inline CvComplex32f& operator += (CvComplex32f& a, CvComplex32f b)
{
a.re += b.re;
a.im += b.im;
return a;
}
inline CvComplex32f operator - (CvComplex32f a, CvComplex32f b)
{
return CvComplex32f( a.re - b.re, a.im - b.im );
}
inline CvComplex32f& operator -= (CvComplex32f& a, CvComplex32f b)
{
a.re -= b.re;
a.im -= b.im;
return a;
}
inline CvComplex32f operator - (CvComplex32f a)
{
return CvComplex32f( -a.re, -a.im );
}
inline CvComplex32f operator * (CvComplex32f a, CvComplex32f b)
{
return CvComplex32f( a.re*b.re - a.im*b.im, a.re*b.im + a.im*b.re );
}
inline double abs(CvComplex32f a)
{
return sqrt( (double)a.re*a.re + (double)a.im*a.im );
}
inline CvComplex32f conj(CvComplex32f a)
{
return CvComplex32f( a.re, -a.im );
}
inline CvComplex32f operator / (CvComplex32f a, CvComplex32f b)
{
double t = 1./((double)b.re*b.re + (double)b.im*b.im);
return CvComplex32f( (float)((a.re*b.re + a.im*b.im)*t),
(float)((-a.re*b.im + a.im*b.re)*t) );
}
inline CvComplex32f operator * (double a, CvComplex32f b)
{
return CvComplex32f( (float)(a*b.re), (float)(a*b.im) );
}
inline CvComplex32f operator * (CvComplex32f a, double b)
{
return CvComplex32f( (float)(a.re*b), (float)(a.im*b) );
}
inline CvComplex32f::operator CvComplex64f() const
{
return CvComplex64f(re,im);
}
inline CvComplex64f operator + (CvComplex64f a, CvComplex64f b)
{
return CvComplex64f( a.re + b.re, a.im + b.im );
}
inline CvComplex64f& operator += (CvComplex64f& a, CvComplex64f b)
{
a.re += b.re;
a.im += b.im;
return a;
}
inline CvComplex64f operator - (CvComplex64f a, CvComplex64f b)
{
return CvComplex64f( a.re - b.re, a.im - b.im );
}
inline CvComplex64f& operator -= (CvComplex64f& a, CvComplex64f b)
{
a.re -= b.re;
a.im -= b.im;
return a;
}
inline CvComplex64f operator - (CvComplex64f a)
{
return CvComplex64f( -a.re, -a.im );
}
inline CvComplex64f operator * (CvComplex64f a, CvComplex64f b)
{
return CvComplex64f( a.re*b.re - a.im*b.im, a.re*b.im + a.im*b.re );
}
inline double abs(CvComplex64f a)
{
return sqrt( (double)a.re*a.re + (double)a.im*a.im );
}
inline CvComplex64f operator / (CvComplex64f a, CvComplex64f b)
{
double t = 1./((double)b.re*b.re + (double)b.im*b.im);
return CvComplex64f( (a.re*b.re + a.im*b.im)*t,
(-a.re*b.im + a.im*b.re)*t );
}
inline CvComplex64f operator * (double a, CvComplex64f b)
{
return CvComplex64f( a*b.re, a*b.im );
}
inline CvComplex64f operator * (CvComplex64f a, double b)
{
return CvComplex64f( a.re*b, a.im*b );
}
inline CvComplex64f::operator CvComplex32f() const
{
return CvComplex32f((float)re,(float)im);
}
inline CvComplex64f conj(CvComplex64f a)
{
return CvComplex64f( a.re, -a.im );
}
inline CvComplex64f operator + (CvComplex64f a, CvComplex32f b)
{
return CvComplex64f( a.re + b.re, a.im + b.im );
}
inline CvComplex64f operator + (CvComplex32f a, CvComplex64f b)
{
return CvComplex64f( a.re + b.re, a.im + b.im );
}
inline CvComplex64f operator - (CvComplex64f a, CvComplex32f b)
{
return CvComplex64f( a.re - b.re, a.im - b.im );
}
inline CvComplex64f operator - (CvComplex32f a, CvComplex64f b)
{
return CvComplex64f( a.re - b.re, a.im - b.im );
}
inline CvComplex64f operator * (CvComplex64f a, CvComplex32f b)
{
return CvComplex64f( a.re*b.re - a.im*b.im, a.re*b.im + a.im*b.re );
}
inline CvComplex64f operator * (CvComplex32f a, CvComplex64f b)
{
return CvComplex64f( a.re*b.re - a.im*b.im, a.re*b.im + a.im*b.re );
}
typedef CvStatus (CV_STDCALL * CvCopyMaskFunc)(const void* src, int src_step,
void* dst, int dst_step, CvSize size,
const void* mask, int mask_step);
CvCopyMaskFunc icvGetCopyMaskFunc( int elem_size );
CvStatus CV_STDCALL icvSetZero_8u_C1R( uchar* dst, int dststep, CvSize size );
CvStatus CV_STDCALL icvScale_32f( const float* src, float* dst, int len, float a, float b );
CvStatus CV_STDCALL icvScale_64f( const double* src, double* dst, int len, double a, double b );
CvStatus CV_STDCALL icvLUT_Transform8u_8u_C1R( const uchar* src, int srcstep, uchar* dst,
int dststep, CvSize size, const uchar* lut );
CvStatus CV_STDCALL icvLUT_Transform8u_16u_C1R( const uchar* src, int srcstep, ushort* dst,
int dststep, CvSize size, const ushort* lut );
CvStatus CV_STDCALL icvLUT_Transform8u_32s_C1R( con
4000
st uchar* src, int srcstep, int* dst,
int dststep, CvSize size, const int* lut );
CvStatus CV_STDCALL icvLUT_Transform8u_64f_C1R( const uchar* src, int srcstep, double* dst,
int dststep, CvSize size, const double* lut );
CvStatus CV_STDCALL icvLUT_Transform8u_8u_C2R( const uchar* src, int srcstep, uchar* dst,
int dststep, CvSize size, const uchar* lut );
CvStatus CV_STDCALL icvLUT_Transform8u_8u_C3R( const uchar* src, int srcstep, uchar* dst,
int dststep, CvSize size, const uchar* lut );
CvStatus CV_STDCALL icvLUT_Transform8u_8u_C4R( const uchar* src, int srcstep, uchar* dst,
int dststep, CvSize size, const uchar* lut );
typedef CvStatus (CV_STDCALL * CvLUT_TransformFunc)( const void* src, int srcstep, void* dst,
int dststep, CvSize size, const void* lut );
CV_INLINE CvStatus
icvLUT_Transform8u_8s_C1R( const uchar* src, int srcstep, char* dst,
int dststep, CvSize size, const char* lut )
{
return icvLUT_Transform8u_8u_C1R( src, srcstep, (uchar*)dst,
dststep, size, (const uchar*)lut );
}
CV_INLINE CvStatus
icvLUT_Transform8u_16s_C1R( const uchar* src, int srcstep, short* dst,
int dststep, CvSize size, const short* lut )
{
return icvLUT_Transform8u_16u_C1R( src, srcstep, (ushort*)dst,
dststep, size, (const ushort*)lut );
}
CV_INLINE CvStatus
icvLUT_Transform8u_32f_C1R( const uchar* src, int srcstep, float* dst,
int dststep, CvSize size, const float* lut )
{
return icvLUT_Transform8u_32s_C1R( src, srcstep, (int*)dst,
dststep, size, (const int*)lut );
}
#endif /*_CXCORE_INTERNAL_H_*/ _cxipp_my.h#ifndef _CXCORE_IPP_H_
#define _CXCORE_IPP_H_
/****************************************************************************************\
* Copy/Set *
\****************************************************************************************/
/* temporary disable ipp zero and copy functions as they affect subsequent functions' performance */
IPCVAPI_EX( CvStatus, icvCopy_8u_C1R, "ippiCopy_8u_C1R", 0/*CV_PLUGINS1(CV_PLUGIN_IPPI)*/,
( const uchar* src, int src_step,
uchar* dst, int dst_step, CvSize size ))
IPCVAPI_EX( CvStatus, icvSetByte_8u_C1R, "ippiSet_8u_C1R", 0/*CV_PLUGINS1(CV_PLUGIN_IPPI)*/,
( uchar value, uchar* dst, int dst_step, CvSize size ))
IPCVAPI_EX( CvStatus, icvCvt_32f64f, "ippsConvert_32f64f",
CV_PLUGINS1(CV_PLUGIN_IPPS), ( const float* src, double* dst, int len ))
IPCVAPI_EX( CvStatus, icvCvt_64f32f, "ippsConvert_64f32f",
CV_PLUGINS1(CV_PLUGIN_IPPS), ( const double* src, float* dst, int len ))
#define IPCV_COPYSET( flavor, arrtype, scalartype ) \
IPCVAPI_EX( CvStatus, icvCopy##flavor, "ippiCopy" #flavor, \
CV_PLUGINS1(CV_PLUGIN_IPPI), \
( const arrtype* src, int srcstep, \
arrtype* dst, int dststep, CvSize size, \
const uchar* mask, int maskstep )) \
IPCVAPI_EX( CvStatus, icvSet##flavor, "ippiSet" #flavor, \
0/*CV_PLUGINS1(CV_PLUGIN_OPTCV)*/, \
( arrtype* dst, int dststep, \
const uchar* mask, int maskstep, \
CvSize size, const arrtype* scalar ))
IPCV_COPYSET( _8u_C1MR, uchar, int )
IPCV_COPYSET( _16s_C1MR, ushort, int )
IPCV_COPYSET( _8u_C3MR, uchar, int )
IPCV_COPYSET( _8u_C4MR, int, int )
IPCV_COPYSET( _16s_C3MR, ushort, int )
IPCV_COPYSET( _16s_C4MR, int64, int64 )
IPCV_COPYSET( _32f_C3MR, int, int )
IPCV_COPYSET( _32f_C4MR, int, int )
IPCV_COPYSET( _64s_C3MR, int64, int64 )
IPCV_COPYSET( _64s_C4MR, int64, int64 )
/****************************************************************************************\
* Arithmetics *
\****************************************************************************************/
#define IPCV_BIN_ARITHM( name ) \
IPCVAPI_EX( CvStatus, icv##name##_8u_C1R, \
"ippi" #name "_8u_C1RSfs", CV_PLUGINS1(CV_PLUGIN_IPPI), \
( const uchar* src1, int srcstep1, const uchar* src2, int srcstep2, \
uchar* dst, int dststep, CvSize size, int scalefactor )) \
IPCVAPI_EX( CvStatus, icv##name##_16u_C1R, \
"ippi" #name "_16u_C1RSfs", CV_PLUGINS1(CV_PLUGIN_IPPI), \
( const ushort* src1, int srcstep1, const ushort* src2, int srcstep2,\
ushort* dst, int dststep, CvSize size, int scalefactor )) \
IPCVAPI_EX( CvStatus, icv##name##_16s_C1R, \
"ippi" #name "_16s_C1RSfs", CV_PLUGINS1(CV_PLUGIN_IPPI), \
( const short* src1, int srcstep1, const short* src2, int srcstep2, \
short* dst, int dststep, CvSize size, int scalefactor )) \
IPCVAPI_EX( CvStatus, icv##name##_32s_C1R, \
"ippi" #name "_32s_C1R", CV_PLUGINS1(CV_PLUGIN_IPPI), \
( const int* src1, int srcstep1, const int* src2, int srcstep2, \
int* dst, int dststep, CvSize size )) \
IPCVAPI_EX( CvStatus, icv##name##_32f_C1R, \
"ippi" #name "_32f_C1R", CV_PLUGINS1(CV_PLUGIN_IPPI), \
( const float* src1, int srcstep1, const float* src2, int srcstep2, \
float* dst, int dststep, CvSize size )) \
IPCVAPI_EX( CvStatus, icv##name##_64f_C1R, \
"ippi" #name "_64f_C1R", CV_PLUGINS1(CV_PLUGIN_IPPI), \
( const double* src1, int srcstep1, const double* src2, int srcstep2,\
double* dst, int dststep, CvSize size ))
IPCV_BIN_ARITHM( Add )
IPCV_BIN_ARITHM( Sub )
#undef IPCV_BIN_ARITHM
/****************************************************************************************\
* Logical operations *
\****************************************************************************************/
#define IPCV_LOGIC( name ) \
IPCVAPI_EX( CvStatus, icv##name##_8u_C1R, \
"ippi" #name "_8u_C1R", 0/*CV_PLUGINS1(CV_PLUGIN_IPPI)*/, \
( const uchar* src1, int srcstep1, const uchar* src2, int srcstep2, \
uchar* dst, int dststep, CvSize size ))
IPCV_LOGIC( And )
IPCV_LOGIC( Or )
IPCV_LOGIC( Xor )
#undef IPCV_LOGIC
IPCVAPI_EX( CvStatus, icvNot_8u_C1R, "ippiNot_8u_C1R", CV_PLUGINS1(CV_PLUGIN_IPPI),
( const uchar* src, int step1, uchar* dst, int step, CvSize size ))
/****************************************************************************************\
* Image Statistics *
\****************************************************************************************/
///////////////////////////////////////// Mean //////////////////////////////////////////
#define IPCV_DEF_MEAN_MASK( flavor, srctype ) \
IPCVAPI_EX( CvStatus, icvMean_##flavor##_C1MR, \
"ippiMean_" #flavor "_C1MR", CV_PLUGINS1(CV_PLUGIN_IPPCV), \
( const srctype* img, int imgstep, const uchar* mask, \
int maskStep, CvSize size, double* mean )) \
IPCVAPI_EX( CvStatus, icvMean_##flavor##_C2MR, \
"ippiMean_" #flavor "_C2MR", 0/*CV_PLUGINS1(CV_PLUGIN_OPTCV)*/, \
( const srctype* img, int imgstep, const uchar* mask, \
int maskStep, CvSize size, double* mean )) \
IPCVAPI_EX( CvStatus, icvMean_##flavor##_C3MR, \
"ippiMean_" #flavor "_C3MR", 0/*CV_PLUGINS1(CV_PLUGIN_OPTCV)*/, \
( const srctype* img, int imgstep, const uchar* mask, \
int maskStep, CvSize size, double* mean )) \
IPCVAPI_EX( CvStatus, icvMean_##flavor##_C4MR, \
"ippiMean_" #flavor "_C4MR", 0/*CV_PLUGINS1(CV_PLUGIN_OPTCV)*/, \
( const srctype* img, int imgstep, const uchar* mask, \
int maskStep, CvSize size, double* mean ))
IPCV_DEF_MEAN_MASK( 8u, uchar )
IPCV_DEF_MEAN_MASK( 16u, ushort )
IPCV_DEF_MEAN_MASK( 16s, short )
IPCV_DEF_MEAN_MASK( 32s, int )
IPCV_DEF_MEAN_MASK( 32f, float )
IPCV_DEF_MEAN_MASK( 64f, double )
#undef IPCV_DEF_MEAN_MASK
//////////////////////////////////// Mean_StdDev ////////////////////////////////////////
#undef IPCV_MEAN_SDV_PLUGIN
#define ICV_MEAN_SDV_PLUGIN 0 /* CV_PLUGINS1(IPPCV) */
#define IPCV_DEF_MEAN_SDV( flavor, srctype ) \
IPCVAPI_EX( CvStatus, icvMean_StdDev_##flavor##_C1R, \
"ippiMean_StdDev_" #flavor "_C1R", ICV_MEAN_SDV_PLUGIN, \
( const srctype* img, int imgstep, CvSize size, double* mean, double* sdv ))\
IPCVAPI_EX( CvStatus, icvMean_StdDev_##flavor##_C2R, \
"ippiMean_StdDev_" #flavor "_C2R", ICV_MEAN_SDV_PLUGIN, \
( const srctype* img, int imgstep, CvSize size, double* mean, double* sdv ))\
IPCVAPI_EX( CvStatus, icvMean_StdDev_##flavor##_C3R, \
"ippiMean_StdDev_" #flavor "_C3R", ICV_MEAN_SDV_PLUGIN, \
( const srctype* img, int imgstep, CvSize size, double* mean, double* sdv ))\
IPCVAPI_EX( CvStatus, icvMean_StdDev_##flavor##_C4R, \
"ippiMean_StdDev_" #flavor "_C4R", ICV_MEAN_SDV_PLUGIN, \
( const srctype* img, int imgstep, CvSize size, double* mean, double* sdv ))\
\
IPCVAPI_EX( CvStatus, icvMean_StdDev_##flavor##_C1MR, \
"ippiMean_StdDev_" #flavor "_C1MR", ICV_MEAN_SDV_PLUGIN, \
( const srctype* img, int imgstep, \
const uchar* mask, int maskStep, \
CvSize size, double* mean, double* sdv )) \
IPCVAPI_EX( CvStatus, icvMean_StdDev_##flavor##_C2MR, \
"ippiMean_StdDev_" #flavor "_C2MR", ICV_MEAN_SDV_PLUGIN, \
( const srctype* img, int imgstep, const uchar* mask, int maskStep, \
CvSize size, double* mean, double* sdv )) \
IPCVAPI_EX( CvStatus, icvMean_StdDev_##flavor##_C3MR, \
"ippiMean_StdDev_" #flavor "_C3MR", ICV_MEAN_SDV_PLUGIN, \
( const srctype* img, int imgstep, \
const uchar* mask, int maskStep, \
CvSize size, double* mean, double* sdv )) \
IPCVAPI_EX( CvStatus, icvMean_StdDev_##flavor##_C4MR, \
"ippiMean_StdDev_" #flavor "_C4MR", ICV_MEAN_SDV_PLUGIN, \
( const srctype* img, int imgstep, \
const uchar* mask, int maskStep, \
CvSize size, double* mean, double* sdv ))
IPCV_DEF_MEAN_SDV( 8u, uchar )
IPCV_DEF_MEAN_SDV( 16u, ushort )
IPCV_DEF_MEAN_SDV( 16s, short )
IPCV_DEF_MEAN_SDV( 32s, int )
IPCV_DEF_MEAN_SDV( 32f, float )
IPCV_DEF_MEAN_SDV( 64f, double )
#undef IPCV_DEF_MEAN_SDV
#undef IPCV_MEAN_SDV_PLUGIN
//////////////////////////////////// MinMaxIndx /////////////////////////////////////////
#define IPCV_DEF_MIN_MAX_LOC( flavor, srctype, extrtype ) \
IPCVAPI_EX( CvStatus, icvMinMaxIndx_##flavor##_C1R, \
"ippiMinMaxIndx_" #flavor "_C1R", CV_PLUGINS1(CV_PLUGIN_IPPCV), \
( const srctype* img, int imgstep, \
CvSize size, extrtype* minVal, extrtype* maxVal, \
CvPoint* minLoc, CvPoint* maxLoc )) \
\
IPCVAPI_EX( CvStatus, icvMinMaxIndx_##flavor##_C1MR, \
"ippiMinMaxIndx_" #flavor "_C1MR", CV_PLUGINS1(CV_PLUGIN_IPPCV),\
( const srctype* img, int imgstep, \
const uchar* mask, int maskStep, \
CvSize size, extrtype* minVal, extrtype* maxVal, \
CvPoint* minLoc, CvPoint* maxLoc ))
IPCV_DEF_MIN_MAX_LOC( 8u, uchar, float )
IPCV_DEF_MIN_MAX_LOC( 16u, ushort, float )
IPCV_DEF_MIN_MAX_LOC( 16s, short, float )
IPCV_DEF_MIN_MAX_LOC( 32s, int, double )
IPCV_DEF_MIN_MAX_LOC( 32f, int, float )
IPCV_DEF_MIN_MAX_LOC( 64f, int64, double )
#undef IPCV_DEF_MIN_MAX_LOC
////////////////////////////////////////// Sum //////////////////////////////////////////
#define IPCV_DEF_SUM_NOHINT( flavor, srctype ) \
IPCVAPI_EX( CvStatus, icvSum_##flavor##_C1R, \
"ippiSum_" #flavor "_C1R", CV_PLUGINS1(CV_PLUGIN_IPPI), \
( const srctype* img, int imgstep, \
CvSize size, double* sum )) \
IPCVAPI_EX( CvStatus, icvSum_##flavor##_C2R, \
"ippiSum_" #flavor "_C2R", CV_PLUGINS1(CV_PLUGIN_IPPI), \
( const srctype* img, int imgstep, \
CvSize size, double* sum )) \
IPCVAPI_EX( CvStatus, icvSum_##flavor##_C3R, \
"ippiSum_" #flavor "_C3R", CV_PLUGINS1(CV_PLUGIN_IPPI), \
( const srctype* img, int imgstep, \
CvSize size, double* sum )) \
IPCVAPI_EX( CvStatus, icvSum_##flavor##_C4R, \
"ippiSum_" #flavor "_C4R", CV_PLUGINS1(CV_PLUGIN_IPPI), \
( const srctype* img, int imgstep, \
CvSize size, double* sum ))
#define IPCV_DEF_SUM_NOHINT_NO_IPP( flavor, srctype ) \
IPCVAPI_EX( CvStatus, icvSum_##flavor##_C1R, \
"ippiSum_" #flavor "_C1R", 0, \
( const srctype* img, int imgstep, \
CvSize size, double* sum )) \
IPCVAPI_EX( CvStatus, icvSum_##flavor##_C2R, \
"ippiSum_" #flavor "_C2R", 0, \
( const srctype* img, int imgstep, \
CvSize size, double* sum )) \
IPCVAPI_EX( CvStatus, icvSum_##flavor##_C3R, \
"ippiSum_" #flavor "_C3R", 0, \
( const srctype* img, int imgstep, \
CvSize size, double* sum )) \
IPCVAPI_EX( CvStatus, icvSum_##flavor##_C4R, \
"ippiSum_" #flavor "_C4R", 0, \
( const srctype* img, int imgstep, \
CvSize size, double* sum ))
#define IPCV_DEF_SUM_HINT( flavor, srctype ) \
IPCVAPI_EX( CvStatus, icvSum_##flavor##_C1R, \
"ippiSum_" #flavor "_C1R", CV_PLUGINS1(CV_PLUGIN_IPPI), \
( const srctype* img, int imgstep, \
CvSize size, double* sum, CvHintAlgorithm )) \
IPCVAPI_EX( CvStatus, icvSum_##flavor##_C2R, \
"ippiSum_" #flavor "_C2R", CV_PLUGINS1(CV_PLUGIN_IPPI)
17a7a
, \
( const srctype* img, int imgstep, \
CvSize size, double* sum, CvHintAlgorithm )) \
IPCVAPI_EX( CvStatus, icvSum_##flavor##_C3R, \
"ippiSum_" #flavor "_C3R", CV_PLUGINS1(CV_PLUGIN_IPPI), \
( const srctype* img, int imgstep, \
CvSize size, double* sum, CvHintAlgorithm )) \
IPCVAPI_EX( CvStatus, icvSum_##flavor##_C4R, \
"ippiSum_" #flavor "_C4R", CV_PLUGINS1(CV_PLUGIN_IPPI), \
( const srctype* img, int imgstep, \
CvSize size, double* sum, CvHintAlgorithm ))
IPCV_DEF_SUM_NOHINT( 8u, uchar )
IPCV_DEF_SUM_NOHINT( 16s, short )
IPCV_DEF_SUM_NOHINT_NO_IPP( 16u, ushort )
IPCV_DEF_SUM_NOHINT( 32s, int )
IPCV_DEF_SUM_HINT( 32f, float )
IPCV_DEF_SUM_NOHINT( 64f, double )
#undef IPCV_DEF_SUM_NOHINT
#undef IPCV_DEF_SUM_HINT
////////////////////////////////////////// CountNonZero /////////////////////////////////
#define IPCV_DEF_NON_ZERO( flavor, srctype ) \
IPCVAPI_EX( CvStatus, icvCountNonZero_##flavor##_C1R, \
"ippiCountNonZero_" #flavor "_C1R", 0/*CV_PLUGINS1(CV_PLUGIN_OPTCV)*/, \
( const srctype* img, int imgstep, CvSize size, int* nonzero ))
IPCV_DEF_NON_ZERO( 8u, uchar )
IPCV_DEF_NON_ZERO( 16s, ushort )
IPCV_DEF_NON_ZERO( 32s, int )
IPCV_DEF_NON_ZERO( 32f, int )
IPCV_DEF_NON_ZERO( 64f, int64 )
#undef IPCV_DEF_NON_ZERO
////////////////////////////////////////// Norms /////////////////////////////////
#define IPCV_DEF_NORM_NOHINT_C1( flavor, srctype ) \
IPCVAPI_EX( CvStatus, icvNorm_Inf_##flavor##_C1R, \
"ippiNorm_Inf_" #flavor "_C1R", CV_PLUGINS1(CV_PLUGIN_IPPI), \
( const srctype* img, int imgstep, \
CvSize size, double* norm )) \
IPCVAPI_EX( CvStatus, icvNorm_L1_##flavor##_C1R, \
"ippiNorm_L1_" #flavor "_C1R", CV_PLUGINS1(CV_PLUGIN_IPPI), \
( const srctype* img, int imgstep, \
CvSize size, double* norm )) \
IPCVAPI_EX( CvStatus, icvNorm_L2_##flavor##_C1R, \
"ippiNorm_L2_" #flavor "_C1R", 0/*CV_PLUGINS1(CV_PLUGIN_IPPI)*/, \
( const srctype* img, int imgstep, \
CvSize size, double* norm )) \
IPCVAPI_EX( CvStatus, icvNormDiff_Inf_##flavor##_C1R, \
"ippiNormDiff_Inf_" #flavor "_C1R", 0/*CV_PLUGINS1(CV_PLUGIN_IPPI)*/, \
( const srctype* img1, int imgstep1, \
const srctype* img2, int imgstep2, \
CvSize size, double* norm )) \
IPCVAPI_EX( CvStatus, icvNormDiff_L1_##flavor##_C1R, \
"ippiNormDiff_L1_" #flavor "_C1R", CV_PLUGINS1(CV_PLUGIN_IPPI), \
( const srctype* img1, int imgstep1, \
const srctype* img2, int imgstep2, \
CvSize size, double* norm )) \
IPCVAPI_EX( CvStatus, icvNormDiff_L2_##flavor##_C1R, \
"ippiNormDiff_L2_" #flavor "_C1R", 0/*CV_PLUGINS1(CV_PLUGIN_IPPI)*/, \
( const srctype* img1, int imgstep1, \
const srctype* img2, int imgstep2, \
CvSize size, double* norm ))
#define IPCV_DEF_NORM_HINT_C1( flavor, srctype ) \
IPCVAPI_EX( CvStatus, icvNorm_Inf_##flavor##_C1R, \
"ippiNorm_Inf_" #flavor "_C1R", CV_PLUGINS1(CV_PLUGIN_IPPI), \
( const srctype* img, int imgstep, \
CvSize size, double* norm )) \
IPCVAPI_EX( CvStatus, icvNorm_L1_##flavor##_C1R, \
"ippiNorm_L1_" #flavor "_C1R", CV_PLUGINS1(CV_PLUGIN_IPPI), \
( const srctype* img, int imgstep, \
CvSize size, double* norm, CvHintAlgorithm )) \
IPCVAPI_EX( CvStatus, icvNorm_L2_##flavor##_C1R, \
"ippiNorm_L2_" #flavor "_C1R", CV_PLUGINS1(CV_PLUGIN_IPPI), \
( const srctype* img, int imgstep, \
CvSize size, double* norm, CvHintAlgorithm )) \
IPCVAPI_EX( CvStatus, icvNormDiff_Inf_##flavor##_C1R, \
"ippiNormDiff_Inf_" #flavor "_C1R", CV_PLUGINS1(CV_PLUGIN_IPPI), \
( const srctype* img1, int imgstep1, \
const srctype* img2, int imgstep2, \
CvSize size, double* norm )) \
IPCVAPI_EX( CvStatus, icvNormDiff_L1_##flavor##_C1R, \
"ippiNormDiff_L1_" #flavor "_C1R", CV_PLUGINS1(CV_PLUGIN_IPPI), \
( const srctype* img1, int imgstep1, \
const srctype* img2, int imgstep2, \
CvSize size, double* norm, CvHintAlgorithm )) \
IPCVAPI_EX( CvStatus, icvNormDiff_L2_##flavor##_C1R, \
"ippiNormDiff_L2_" #flavor "_C1R", CV_PLUGINS1(CV_PLUGIN_IPPI), \
( const srctype* img1, int imgstep1, \
const srctype* img2, int imgstep2, \
CvSize size, double* norm, CvHintAlgorithm ))
#define IPCV_DEF_NORM_MASK_C1( flavor, srctype ) \
IPCVAPI_EX( CvStatus, icvNorm_Inf_##flavor##_C1MR, \
"ippiNorm_Inf_" #flavor "_C1MR", CV_PLUGINS1(CV_PLUGIN_IPPCV), \
( const srctype* img, int imgstep, \
const uchar* mask, int maskstep, \
CvSize size, double* norm )) \
IPCVAPI_EX( CvStatus, icvNorm_L1_##flavor##_C1MR, \
"ippiNorm_L1_" #flavor "_C1MR", CV_PLUGINS1(CV_PLUGIN_IPPCV), \
( const srctype* img, int imgstep, \
const uchar* mask, int maskstep, \
CvSize size, double* norm )) \
IPCVAPI_EX( CvStatus, icvNorm_L2_##flavor##_C1MR, \
"ippiNorm_L2_" #flavor "_C1MR", CV_PLUGINS1(CV_PLUGIN_IPPCV), \
( const srctype* img, int imgstep, \
const uchar* mask, int maskstep, \
CvSize size, double* norm )) \
IPCVAPI_EX( CvStatus, icvNormDiff_Inf_##flavor##_C1MR, \
"ippiNormDiff_Inf_" #flavor "_C1MR", CV_PLUGINS1(CV_PLUGIN_IPPCV), \
( const srctype* img1, int imgstep1, \
const srctype* img2, int imgstep2, \
const uchar* mask, int maskstep, \
CvSize size, double* norm )) \
IPCVAPI_EX( CvStatus, icvNormDiff_L1_##flavor##_C1MR, \
"ippiNormDiff_L1_" #flavor "_C1MR", CV_PLUGINS1(CV_PLUGIN_IPPCV), \
( const srctype* img1, int imgstep1, \
const srctype* img2, int imgstep2, \
const uchar* mask, int maskstep, \
CvSize size, double* norm )) \
IPCVAPI_EX( CvStatus, icvNormDiff_L2_##flavor##_C1MR, \
"ippiNormDiff_L2_" #flavor "_C1MR", CV_PLUGINS1(CV_PLUGIN_IPPCV), \
( const srctype* img1, int imgstep1, \
const srctype* img2, int imgstep2, \
const uchar* mask, int maskstep, \
CvSize size, double* norm ))
IPCV_DEF_NORM_NOHINT_C1( 8u, uchar )
IPCV_DEF_NORM_MASK_C1( 8u, uchar )
IPCV_DEF_NORM_NOHINT_C1( 16u, ushort )
IPCV_DEF_NORM_MASK_C1( 16u, ushort )
IPCV_DEF_NORM_NOHINT_C1( 16s, short )
IPCV_DEF_NORM_MASK_C1( 16s, short )
IPCV_DEF_NORM_NOHINT_C1( 32s, int )
IPCV_DEF_NORM_MASK_C1( 32s, int )
IPCV_DEF_NORM_HINT_C1( 32f, float )
IPCV_DEF_NORM_MASK_C1( 32f, float )
IPCV_DEF_NORM_NOHINT_C1( 64f, double )
IPCV_DEF_NORM_MASK_C1( 64f, double )
#undef IPCV_DEF_NORM_HONINT_C1
#undef IPCV_DEF_NORM_HINT_C1
#undef IPCV_DEF_NORM_MASK_C1
////////////////////////////////////// AbsDiff ///////////////////////////////////////////
#define IPCV_ABS_DIFF( flavor, arrtype ) \
IPCVAPI_EX( CvStatus, icvAbsDiff_##flavor##_C1R, \
"ippiAbsDiff_" #flavor "_C1R", CV_PLUGINS1(CV_PLUGIN_IPPCV),\
( const arrtype* src1, int srcstep1, \
const arrtype* src2, int srcstep2, \
arrtype* dst, int dststep, CvSize size ))
IPCV_ABS_DIFF( 8u, uchar )
IPCV_ABS_DIFF( 16u, ushort )
IPCV_ABS_DIFF( 16s, short )
IPCV_ABS_DIFF( 32s, int )
IPCV_ABS_DIFF( 32f, float )
IPCV_ABS_DIFF( 64f, double )
#undef IPCV_ABS_DIFF
////////////////////////////// Comparisons //////////////////////////////////////////
#define IPCV_CMP( arrtype, flavor ) \
IPCVAPI_EX( CvStatus, icvCompare_##flavor##_C1R, \
"ippiCompare_" #flavor "_C1R", CV_PLUGINS1(CV_PLUGIN_IPPI), \
( const arrtype* src1, int srcstep1, const arrtype* src2, int srcstep2, \
arrtype* dst, int dststep, CvSize size, int cmp_op )) \
IPCVAPI_EX( CvStatus, icvCompareC_##flavor##_C1R, \
"ippiCompareC_" #flavor "_C1R", CV_PLUGINS1(CV_PLUGIN_IPPI), \
( const arrtype* src1, int srcstep1, arrtype scalar, \
arrtype* dst, int dststep, CvSize size, int cmp_op )) \
IPCVAPI_EX( CvStatus, icvThreshold_GT_##flavor##_C1R, \
"ippiThreshold_GT_" #flavor "_C1R", CV_PLUGINS1(CV_PLUGIN_IPPI), \
( const arrtype* pSrc, int srcstep, arrtype* pDst, int dststep, \
CvSize size, arrtype threshold )) \
IPCVAPI_EX( CvStatus, icvThreshold_LT_##flavor##_C1R, \
"ippiThreshold_LT_" #flavor "_C1R", CV_PLUGINS1(CV_PLUGIN_IPPI), \
( const arrtype* pSrc, int srcstep, arrtype* pDst, int dststep, \
CvSize size, arrtype threshold ))
IPCV_CMP( uchar, 8u )
IPCV_CMP( short, 16s )
IPCV_CMP( float, 32f )
#undef IPCV_CMP
/****************************************************************************************\
* Utilities *
\****************************************************************************************/
////////////////////////////// Copy Pixel <-> Plane /////////////////////////////////
#define IPCV_PIX_PLANE( flavor, arrtype ) \
IPCVAPI_EX( CvStatus, icvCopy_##flavor##_C2P2R, \
"ippiCopy_" #flavor "_C2P2R", 0/*CV_PLUGINS1(CV_PLUGIN_IPPI)*/, \
( const arrtype* src, int srcstep, arrtype** dst, int dststep, CvSize size )) \
IPCVAPI_EX( CvStatus, icvCopy_##flavor##_C3P3R, \
"ippiCopy_" #flavor "_C3P3R", CV_PLUGINS1(CV_PLUGIN_IPPI), \
( const arrtype* src, int srcstep, arrtype** dst, int dststep, CvSize size )) \
IPCVAPI_EX( CvStatus, icvCopy_##flavor##_C4P4R, \
"ippiCopy_" #flavor "_C4P4R", CV_PLUGINS1(CV_PLUGIN_IPPI), \
( const arrtype* src, int srcstep, arrtype** dst, int dststep, CvSize size )) \
IPCVAPI_EX( CvStatus, icvCopy_##flavor##_CnC1CR, \
"ippiCopy_" #flavor "_CnC1CR", 0/*CV_PLUGINS1(CV_PLUGIN_OPTCV)*/, \
( const arrtype* src, int srcstep, arrtype* dst, int dststep, \
CvSize size, int cn, int coi )) \
IPCVAPI_EX( CvStatus, icvCopy_##flavor##_C1CnCR, \
"ippiCopy_" #flavor "_CnC1CR", 0/*CV_PLUGINS1(CV_PLUGIN_OPTCV)*/, \
( const arrtype* src, int srcstep, arrtype* dst, int dststep, \
CvSize size, int cn, int coi )) \
IPCVAPI_EX( CvStatus, icvCopy_##flavor##_P2C2R, \
"ippiCopy_" #flavor "_P2C2R", 0/*CV_PLUGINS1(CV_PLUGIN_IPPI)*/, \
( const arrtype** src, int srcstep, arrtype* dst, int dststep, CvSize size )) \
IPCVAPI_EX( CvStatus, icvCopy_##flavor##_P3C3R, \
"ippiCopy_" #flavor "_P3C3R", CV_PLUGINS1(CV_PLUGIN_IPPI), \
( const arrtype** src, int srcstep, arrtype* dst, int dststep, CvSize size )) \
IPCVAPI_EX( CvStatus, icvCopy_##flavor##_P4C4R, \
"ippiCopy_" #flavor "_P4C4R", CV_PLUGINS1(CV_PLUGIN_IPPI), \
( const arrtype** src, int srcstep, arrtype* dst, int dststep, CvSize size ))
IPCV_PIX_PLANE( 8u, uchar )
IPCV_PIX_PLANE( 16s, ushort )
IPCV_PIX_PLANE( 32f, int )
IPCV_PIX_PLANE( 64f, int64 )
#undef IPCV_PIX_PLANE
/****************************************************************************************/
/* Math routines and RNGs */
/****************************************************************************************/
IPCVAPI_EX( CvStatus, icvInvSqrt_32f, "ippsInvSqrt_32f_A21",
CV_PLUGINS1(CV_PLUGIN_IPPVM),
( const float* src, float* dst, int len ))
IPCVAPI_EX( CvStatus, icvSqrt_32f, "ippsSqrt_32f_A21, ippsSqrt_32f",
CV_PLUGINS2(CV_PLUGIN_IPPVM,CV_PLUGIN_IPPS),
( const float* src, float* dst, int len ))
IPCVAPI_EX( CvStatus, icvInvSqrt_64f, "ippsInvSqrt_64f_A50",
CV_PLUGINS1(CV_PLUGIN_IPPVM),
( const double* src, double* dst, int len ))
IPCVAPI_EX( CvStatus, icvSqrt_64f, "ippsSqrt_64f_A50, ippsSqrt_64f",
CV_PLUGINS2(CV_PLUGIN_IPPVM,CV_PLUGIN_IPPS),
( const double* src, double* dst, int len ))
IPCVAPI_EX( CvStatus, icvLog_32f, "ippsLn_32f_A21, ippsLn_32f",
CV_PLUGINS2(CV_PLUGIN_IPPVM,CV_PLUGIN_IPPS),
( const float *x, float *y, int n ) )
IPCVAPI_EX( CvStatus, icvLog_64f, "ippsLn_64f_A50, ippsLn_64f",
CV_PLUGINS2(CV_PLUGIN_IPPVM,CV_PLUGIN_IPPS),
( const double *x, double *y, int n ) )
IPCVAPI_EX( CvStatus, icvExp_32f, "ippsExp_32f_A21, ippsExp_32f",
CV_PLUGINS2(CV_PLUGIN_IPPVM,CV_PLUGIN_IPPS),
( const float *x, float *y, int n ) )
IPCVAPI_EX( CvStatus, icvExp_64f, "ippsExp_64f_A50, ippsExp_64f",
CV_PLUGINS2(CV_PLUGIN_IPPVM,CV_PLUGIN_IPPS),
( const double *x, double *y, int n ) )
IPCVAPI_EX( CvStatus, icvFastArctan_32f, "ippibFastArctan_32f",
CV_PLUGINS1(CV_PLUGIN_IPPCV),
( const float* y, const float* x, float* angle, int len ))
/****************************************************************************************/
/* Error handling functions */
/****************************************************************************************/
IPCVAPI_EX( CvStatus, icvCheckArray_32f_C1R,
"ippiCheckArray_32f_C1R", 0/*CV_PLUGINS1(CV_PLUGIN_OPTCV)*/,
( const float* src, int srcstep,
CvSize size, int flags,
double min_val, double max_val ))
IPCVAPI_EX( CvStatus, icvCheckArray_64f_C1R,
"ippiCheckArray_64f_C1R", 0/*CV_PLUGINS1(CV_PLUGIN_OPTCV)*/,
( const double* src, int srcstep,
CvSize size, int flags,
double min_val, double max_val ))
/****************************************************************************************/
/* Affine transformations of matrix/image elements */
/****************************************************************************************/
#define IPCV_TRANSFORM( suffix, ipp_suffix, cn ) \
IPCVAPI_EX( CvStatus, icvColorTwist##suffix##_C##cn##R, \
"ippiColorTwist" #ipp_suffix "_C" #cn \
"R,ippiColorTwist" #ipp_suffix "_C" #cn "R", \
CV_PLUGINS2(CV_PLUGIN_IPPI, CV_PLUGIN_IPPCC), \
( const void* src, int srcstep, void* dst, int dststep, \
CvSize roisize, const float* twist_matrix ))
IPCV_TRANSFORM( _8u, 32f_8u, 3 )
IPCV_TRANSFORM( _16u, 32f_16u, 3 )
IPCV_TRANSFORM( _16s, 32f_16s, 3 )
IPCV_TRANSFORM( _32f, _32f, 3 )
IPCV_TRANSFORM( _32f, _32f, 4 )
#undef IPCV_TRANSFORM
#define IPCV_TRANSFORM_N1( suffix ) \
IPCVAPI_EX( CvStatus, icvColorToGray##suffix, \
"ippiColorToGray" #suffix ",ippiColorToGray" #suffix, \
CV_PLUGINS2(CV_PLUGIN_IPPI,CV_PLUGIN_IPPCC), \
( const void* src, int srcstep, void* dst, int dststep, \
CvSize roisize, const float* coeffs ))
IPCV_TRANSFORM_N1( _8u_C3C1R )
IPCV_TRANSFORM_N1( _16u_C3C1R )
IPCV_TRANSFORM_N1( _16s_C3C1R )
IPCV_TRANSFORM_N1( _32f_C3C1R )
IPCV_TRANSFORM_N1( _8u_AC4C1R )
IPCV_TRANSFORM_N1( _16u_AC4C1R )
IPCV_TRANSFORM_N1( _16s_AC4C1R )
IPCV_TRANSFORM_N1( _32f_AC4C1R )
#undef IPCV_TRANSFORM_N1
/****************************************************************************************/
/* Matrix routines from BLAS/LAPACK compatible libraries */
/****************************************************************************************/
IPCVAPI_C_EX( void, icvBLAS_GEMM_32f, "sgemm, mkl_sgemm", CV_PLUGINS2(CV_PLUGIN_MKL,CV_PLUGIN_MKL),
(const char *transa, const char *transb, int *n, int *m, int *k,
const void *alpha, const void *a, int *lda, const void *b, int *ldb,
const void *beta, void *c, int *ldc ))
IPCVAPI_C_EX( void, icvBLAS_GEMM_64f, "dgemm, mkl_dgemm", CV_PLUGINS2(CV_PLUGIN_MKL,CV_PLUGIN_MKL),
(const char *transa, const char *transb, int *n, int *m, int *k,
const void *alpha, const void *a, int *lda, const void *b, int *ldb,
const void *beta, void *c, int *ldc ))
IPCVAPI_C_EX( void, icvBLAS_GEMM_32fc, "cgemm, mkl_cgemm", CV_PLUGINS2(CV_PLUGIN_MKL,CV_PLUGIN_MKL),
(const char *transa, const char *transb, int *n, int *m, int *k,
const void *alpha, const void *a, int *lda, const void *b, int *ldb,
const void *beta, void *c, int *ldc ))
IPCVAPI_C_EX( void, icvBLAS_GEMM_64fc, "zgemm, mkl_zgemm", CV_PLUGINS2(CV_PLUGIN_MKL,CV_PLUGIN_MKL),
(const char *transa, const char *transb, int *n, int *m, int *k,
const void *alpha, const void *a, int *lda, const void *b, int *ldb,
const void *beta, void *c, int *ldc ))
#define IPCV_DFT( init_flavor, fwd_flavor, inv_flavor ) \
IPCVAPI_EX( CvStatus, icvDFTInitAlloc_##init_flavor, "ippsDFTInitAlloc_" #init_flavor, \
CV_PLUGINS1(CV_PLUGIN_IPPS), ( void**, int, int, CvHintAlgorithm )) \
\
IPCVAPI_EX( CvStatus, icvDFTFree_##init_flavor, "ippsDFTFree_" #init_flavor, \
CV_PLUGINS1(CV_PLUGIN_IPPS), ( void* )) \
\
IPCVAPI_EX( CvStatus, icvDFTGetBufSize_##init_flavor, "ippsDFTGetBufSize_" #init_flavor,\
CV_PLUGINS1(CV_PLUGIN_IPPS), ( const void* spec, int* buf_size )) \
\
IPCVAPI_EX( CvStatus, icvDFTFwd_##fwd_flavor, "ippsDFTFwd_" #fwd_flavor, \
CV_PLUGINS1(CV_PLUGIN_IPPS), ( const void* src, void* dst, \
const void* spec, void* buffer )) \
\
IPCVAPI_EX( CvStatus, icvDFTInv_##inv_flavor, "ippsDFTInv_" #inv_flavor, \
CV_PLUGINS1(CV_PLUGIN_IPPS), ( const void* src, void* dst, \
const void* spec, void* buffer ))
IPCV_DFT( C_32fc, CToC_32fc, CToC_32fc )
IPCV_DFT( R_32f, RToPack_32f, PackToR_32f )
IPCV_DFT( C_64fc, CToC_64fc, CToC_64fc )
IPCV_DFT( R_64f, RToPack_64f, PackToR_64f )
#undef IPCV_DFT
#endif /*_CXCORE_IPP_H_*/
测试代码#include <STDIO.H>
#include <iostream>
#include "_cxcore_my.h"
#pragma comment(lib, "cxcore.lib")
/////////////////////////////////////////////////////////////////////////////////////////
#define icvGivens_64f( n, x, y, c, s ) \
{ \
int _i; \
double* _x = (x); \
double* _y = (y); \
\
for( _i = 0; _i < n; _i++ ) \
{ \
double t0 = _x[_i]; \
double t1 = _y[_i]; \
_x[_i] = t0*c + t1*s; \
_y[_i] = -t0*s + t1*c; \
} \
}
/* y[0:m,0:n] += diag(a[0:1,0:m]) * x[0:m,0:n] */
static void
icvMatrAXPY_64f( int m, int n, const double* x, int dx,
const double* a, double* y, int dy )
{
int i, j;
for( i = 0; i < m; i++, x += dx, y += dy )
{
double s = a[i];
for( j = 0; j <= n - 4; j += 4 )
{
double t0 = y[j] + s*x[j];
double t1 = y[j+1] + s*x[j+1];
y[j] = t0;
y[j+1] = t1;
t0 = y[j+2] + s*x[j+2];
t1 = y[j+3] + s*x[j+3];
y[j+2] = t0;
y[j+3] = t1;
}
for( ; j < n; j++ ) y[j] += s*x[j];
}
}
/* y[1:m,-1] = h*y[1:m,0:n]*x[0:1,0:n]'*x[-1] (this is used for U&V reconstruction)
y[1:m,0:n] += h*y[1:m,0:n]*x[0:1,0:n]'*x[0:1,0:n] */
static void
icvMatrAXPY3_64f( int m, int n, const double* x, int l, double* y, double h )
{
int i, j;
for( i = 1; i < m; i++ )
{
double s = 0;
y += l;
for( j = 0; j <= n - 4; j += 4 )
s += x[j]*y[j] + x[j+1]*y[j+1] + x[j+2]*y[j+2] + x[j+3]*y[j+3];
for( ; j < n; j++ ) s += x[j]*y[j];
s *= h;
y[-1] = s*x[-1];
for( j = 0; j <= n - 4; j += 4 )
{
double t0 = y[j] + s*x[j];
double t1 = y[j+1] + s*x[j+1];
y[j] = t0;
y[j+1] = t1;
t0 = y[j+2] + s*x[j+2];
t1 = y[j+3] + s*x[j+3];
y[j+2] = t0;
y[j+3] = t1;
}
for( ; j < n; j++ ) y[j] += s*x[j];
}
}
#define icvGivens_32f( n, x, y, c, s ) \
{ \
int _i; \
float* _x = (x); \
float* _y = (y); \
\
for( _i = 0; _i < n; _i++ ) \
{ \
double t0 = _x[_i]; \
double t1 = _y[_i]; \
_x[_i] = (float)(t0*c + t1*s); \
_y[_i] = (float)(-t0*s + t1*c);\
} \
}
static void
icvMatrAXPY_32f( int m, int n, const float* x, int dx,
const float* a, float* y, int dy )
{
int i, j;
for( i = 0; i < m; i++, x += dx, y += dy )
{
double s = a[i];
for( j = 0; j <= n - 4; j += 4 )
{
double t0 = y[j] + s*x[j];
double t1 = y[j+1] + s*x[j+1];
y[j] = (float)t0;
y[j+1] = (float)t1;
t0 = y[j+2] + s*x[j+2];
t1 = y[j+3] + s*x[j+3];
y[j+2] = (float)t0;
y[j+3] = (float)t1;
}
for( ; j < n; j++ )
y[j] = (float)(y[j] + s*x[j]);
}
}
static void
icvMatrAXPY3_32f( int m, int n, const float* x, int l, float* y, double h )
{
int i, j;
for( i = 1; i < m; i++ )
{
double s = 0;
y += l;
for( j = 0; j <= n - 4; j += 4 )
s += x[j]*y[j] + x[j+1]*y[j+1] + x[j+2]*y[j+2] + x[j+3]*y[j+3];
for( ; j < n; j++ ) s += x[j]*y[j];
s *= h;
y[-1] = (float)(s*x[-1]);
for( j = 0; j <= n - 4; j += 4 )
{
double t0 = y[j] + s*x[j];
double t1 = y[j+1] + s*x[j+1];
y[j] = (float)t0;
y[j+1] = (float)t1;
t0 = y[j+2] + s*x[j+2];
t1 = y[j+3] + s*x[j+3];
y[j+2] = (float)t0;
y[j+3] = (float)t1;
}
for( ; j < n; j++ ) y[j] = (float)(y[j] + s*x[j]);
}
}
/* accurate hypotenuse calculation */
static double
pythag( double a, double b )
{
a = fabs( a );
b = fabs( b );
if( a > b )
{
b /= a;
a *= sqrt( 1. + b * b );
}
else if( b != 0 )
{
a /= b;
a = b * sqrt( 1. + a * a );
}
return a;
}
/****************************************************************************************/
/****************************************************************************************/
#define MAX_ITERS 30
static void icvSVD_32f( float* a, int lda, int m, int n,
float* w,
float* uT, int lduT, int nu,
float* vT, int ldvT,
float* buffer )
{
float* e;
float* temp;
float *w1, *e1;
float *hv;
double ku0 = 0, kv0 = 0;
double anorm = 0;
float *a1, *u0 = uT, *v0 = vT;
double scale, h;
int i, j, k, l;
int nm, m1, n1;
int nv = n;
int iters = 0;
float* hv0 = (float*)cvStackAlloc( (m+2)*sizeof(hv0[0])) + 1;
e = buffer;
w1 = w;
e1 = e + 1;
nm = n;
temp = buffer + nm;
memset( w, 0, nm * sizeof( w[0] ));
memset( e, 0, nm * sizeof( e[0] ));
m1 = m;
n1 = n;
/* transform a to bi-diagonal form */
for( ;; )
{
int update_u;
int update_v;
if( m1 == 0 )
break;
scale = h = 0;
update_u = uT && m1 > m - nu;
hv = update_u ? uT : hv0;
for( j = 0, a1 = a; j < m1; j++, a1 += lda )
{
double t = a1[0];
scale += fabs( hv[j] = (float)t );
}
if( scale != 0 )
{
double f = 1./scale, g, s = 0;
for( j = 0; j < m1; j++ )
{
double t = (hv[j] = (float)(hv[j]*f));
s += t * t;
}
g = sqrt( s );
f = hv[0];
if( f >= 0 )
g = -g;
hv[0] = (float)(f - g);
h = 1. / (f * g - s);
memset( temp, 0, n1 * sizeof( temp[0] ));
/* calc temp[0:n-i] = a[i:m,i:n]'*hv[0:m-i] */
icvMatrAXPY_32f( m1, n1 - 1, a + 1, lda, hv, temp + 1, 0 );
for( k = 1; k < n1; k++ ) temp[k] = (float)(temp[k]*h);
/* modify a: a[i:m,i:n] = a[i:m,i:n] + hv[0:m-i]*temp[0:n-i]' */
icvMatrAXPY_32f( m1, n1 - 1, temp + 1, 0, hv, a + 1, lda );
*w1 = (float)(g*scale);
}
w1++;
/* store -2/(hv'*hv) */
if( update_u )
{
if( m1 == m )
ku0 = h;
else
hv[-1] = (float)h;
}
a++;
n1--;
if( vT )
vT += ldvT + 1;
if( n1 == 0 )
break;
scale = h = 0;
update_v = vT && n1 > n - nv;
hv = update_v ? vT : hv0;
for( j = 0; j < n1; j++ )
{
double t = a[j];
scale += fabs( hv[j] = (float)t );
}
if( scale != 0 )
{
double f = 1./scale, g, s = 0;
for( j = 0; j < n1; j++ )
{
double t = (hv[j] = (float)(hv[j]*f));
s += t * t;
}
g = sqrt( s );
f = hv[0];
if( f >= 0 )
g = -g;
hv[0] = (float)(f - g);
h = 1. / (f * g - s);
hv[-1] = 0.f;
/* update a[i:m:i+1:n] = a[i:m,i+1:n] + (a[i:m,i+1:n]*hv[0:m-i])*... */
icvMatrAXPY3_32f( m1, n1, hv, lda, a, h );
*e1 = (float)(g*scale);
}
e1++;
/* store -2/(hv'*hv) */
if( update_v )
{
if( n1 == n )
kv0 = h;
else
hv[-1] = (float)h;
}
a += lda;
m1--;
if( uT )
uT += lduT + 1;
}
m1 -= m1 != 0;
n1 -= n1 != 0;
/* accumulate left transformations */
if( uT )
{
m1 = m - m1;
uT = u0 + m1 * lduT;
for( i = m1; i < nu; i++, uT += lduT )
{
memset( uT + m1, 0, (m - m1) * sizeof( uT[0] ));
uT[i] = 1.;
}
for( i = m1 - 1; i >= 0; i-- )
{
double s;
int lh = nu - i;
l = m - i;
hv = u0 + (lduT + 1) * i;
h = i == 0 ? ku0 : hv[-1];
assert( h <= 0 );
if( h != 0 )
{
uT = hv;
icvMatrAXPY3_32f( lh, l-1, hv+1, lduT, uT+1, h );
s = hv[0] * h;
for( k = 0; k < l; k++ ) hv[k] = (float)(hv[k]*s);
hv[0] += 1;
}
else
{
for( j = 1; j < l; j++ )
hv[j] = 0;
for( j = 1; j < lh; j++ )
hv[j * lduT] = 0;
hv[0] = 1;
}
}
uT = u0;
}
/* accumulate right transformations */
if( vT )
{
n1 = n - n1;
vT = v0 + n1 * ldvT;
for( i = n1; i < nv; i++, vT += ldvT )
{
memset( vT + n1, 0, (n - n1) * sizeof( vT[0] ));
vT[i] = 1.;
}
for( i = n1 - 1; i >= 0; i-- )
{
double s;
int lh = nv - i;
l = n - i;
hv = v0 + (ldvT + 1) * i;
h = i == 0 ? kv0 : hv[-1];
assert( h <= 0 );
if( h != 0 )
{
vT = hv;
icvMatrAXPY3_32f( lh, l-1, hv+1, ldvT, vT+1, h );
s = hv[0] * h;
for( k = 0; k < l; k++ ) hv[k] = (float)(hv[k]*s);
hv[0] += 1;
}
else
{
for( j = 1; j < l; j++ )
hv[j] = 0;
for( j = 1; j < lh; j++ )
hv[j * ldvT] = 0;
hv[0] = 1;
}
}
vT = v0;
}
for( i = 0; i < nm; i++ )
{
double tnorm = fabs( w[i] );
tnorm += fabs( e[i] );
if( anorm < tnorm )
anorm = tnorm;
}
anorm *= FLT_EPSILON;
/* diagonalization of the bidiagonal form */
for( k = nm - 1; k >= 0; k-- )
{
double z = 0;
iters = 0;
for( ;; ) /* do iterations */
{
double c, s, f, g, x, y;
int flag = 0;
/* test for splitting */
for( l = k; l >= 0; l-- )
{
if( fabs( e[l] ) <= anorm )
{
flag = 1;
break;
}
assert( l > 0 );
if( fabs( w[l - 1] ) <= anorm )
break;
}
if( !flag )
{
c = 0;
s = 1;
for( i = l; i <= k; i++ )
{
f = s * e[i];
e[i] = (float)(e[i]*c);
if( anorm + fabs( f ) == anorm )
break;
g = w[i];
h = pythag( f, g );
w[i] = (float)h;
c = g / h;
s = -f / h;
if( uT )
icvGivens_32f( m, uT + lduT * (l - 1), uT + lduT * i, c, s );
}
}
z = w[k];
if( l == k || iters++ == MAX_ITERS )
break;
/* shift from bottom 2x2 minor */
x = w[l];
y = w[k - 1];
g = e[k - 1];
h = e[k];
f = 0.5 * (((g + z) / h) * ((g - z) / y) + y / h - h / y);
g = pythag( f, 1 );
if( f < 0 )
g = -g;
f = x - (z / x) * z + (h / x) * (y / (f + g) - h);
/* next QR transformation */
c = s = 1;
for( i = l + 1; i <= k; i++ )
{
g = e[i];
y = w[i];
h = s * g;
g *= c;
z = pythag( f, h );
e[i - 1] = (float)z;
c = f / z;
s = h / z;
f = x * c + g * s;
g = -x * s + g * c;
h = y * s;
y *= c;
if( vT )
icvGivens_32f( n, vT + ldvT * (i - 1), vT + ldvT * i, c, s );
z = pythag( f, h );
w[i - 1] = (float)z;
/* rotation can be arbitrary if z == 0 */
if( z != 0 )
{
c = f / z;
s = h / z;
}
f = c * g + s * y;
x = -s * g + c * y;
if( uT )
icvGivens_32f( m, uT + lduT * (i - 1), uT + lduT * i, c, s );
}
e[l] = 0;
e[k] = (float)f;
w[k] = (float)x;
} /* end of iteration loop */
if( iters > MAX_ITERS )
break;
if( z < 0 )
{
w[k] = (float)(-z);
if( vT )
{
for( j = 0; j < n; j++ )
vT[j + k * ldvT] = -vT[j + k * ldvT];
}
}
} /* end of diagonalization loop */
/* sort singular values and corresponding vectors */
for( i = 0; i < nm; i++ )
{
k = i;
for( j = i + 1; j < nm; j++ )
if( w[k] < w[j] )
k = j;
if( k != i )
{
float t;
CV_SWAP( w[i], w[k], t );
if( vT )
for( j = 0; j < n; j++ )
CV_SWAP( vT[j + ldvT*k], vT[j + ldvT*i], t );
if( uT )
for( j = 0; j < m; j++ )
CV_SWAP( uT[j + lduT*k], uT[j + lduT*i], t );
}
}
}
static void icvSVD_64f( double* a, int lda, int m, int n,
double* w,
double* uT, int lduT, int nu,
double* vT, int ldvT,
double* buffer )
{
double* e;
double* temp;
double *w1, *e1;
double *hv;
double ku0 = 0, kv0 = 0;
double anorm = 0;
double *a1, *u0 = uT, *v0 = vT;
double scale, h;
int i, j, k, l;
int nm, m1, n1;
int nv = n;
int iters = 0;
double* hv0 = (double*)cvStackAlloc( (m+2)*sizeof(hv0[0])) + 1;
e = buffer;
w1 = w;
e1 = e + 1;
nm = n;
temp = buffer + nm;
memset( w, 0, nm * sizeof( w[0] ));
memset( e, 0, nm * sizeof( e[0] ));
m1 = m;
n1 = n;
/* transform a to bi-diagonal form */
for( ;; )
{
int update_u;
int update_v;
if( m1 == 0 )
break;
scale = h = 0;
update_u = uT && m1 > m - nu;
hv = update_u ? uT : hv0;
for( j = 0, a1 = a; j < m1; j++, a1 += lda )
{
double t = a1[0];
scale += fabs( hv[j] = t );
}
if( scale != 0 )
{
double f = 1./scale, g, s = 0;
for( j = 0; j < m1; j++ )
{
double t = (hv[j] *= f);
s += t * t;
}
g = sqrt( s );
f = hv[0];
if( f >= 0 )
g = -g;
hv[0] = f - g;
h = 1. / (f * g - s);
memset( temp, 0, n1 * sizeof( temp[0] ));
/* calc temp[0:n-i] = a[i:m,i:n]'*hv[0:m-i] */
icvMatrAXPY_64f( m1, n1 - 1, a + 1, lda, hv, temp + 1, 0 );
for( k = 1; k < n1; k++ ) temp[k] *= h;
/* modify a: a[i:m,i:n] = a[i:m,i:n] + hv[0:m-i]*temp[0:n-i]' */
icvMatrAXPY_64f( m1, n1 - 1, temp + 1, 0, hv, a + 1, lda );
*w1 = g*scale;
}
w1++;
/* store -2/(hv'*hv) */
if( update_u )
{
if( m1 == m )
ku0 = h;
else
hv[-1] = h;
}
a++;
n1--;
if( vT )
vT += ldvT + 1;
if( n1 == 0 )
break;
scale = h = 0;
update_v = vT && n1 > n - nv;
hv = update_v ? vT : hv0;
for( j = 0; j < n1; j++ )
{
double t = a[j];
scale += fabs( hv[j] = t );
}
if( scale != 0 )
{
double f = 1./scale, g, s = 0;
for( j = 0; j < n1; j++ )
{
double t = (hv[j] *= f);
s += t * t;
}
g = sqrt( s );
f = hv[0];
if( f >= 0 )
g = -g;
hv[0] = f - g;
h = 1. / (f * g - s);
hv[-1] = 0.;
/* update a[i:m:i+1:n] = a[i:m,i+1:n] + (a[i:m,i+1:n]*hv[0:m-i])*... */
icvMatrAXPY3_64f( m1, n1, hv, lda, a, h );
*e1 = g*scale;
}
e1++;
/* store -2/(hv'*hv) */
if( update_v )
{
if( n1 == n )
kv0 = h;
else
hv[-1] = h;
}
a += lda;
m1--;
if( uT )
uT += lduT + 1;
}
m1 -= m1 != 0;
n1 -= n1 != 0;
/* accumulate left transformations */
if( uT )
{
m1 = m - m1;
uT = u0 + m1 * lduT;
for( i = m1; i < nu; i++, uT += lduT )
{
memset( uT + m1, 0, (m - m1) * sizeof( uT[0] ));
uT[i] = 1.;
}
for( i = m1 - 1; i >= 0; i-- )
{
double s;
int lh = nu - i;
l = m - i;
hv = u0 + (lduT + 1) * i;
h = i == 0 ? ku0 : hv[-1];
assert( h <= 0 );
if( h != 0 )
{
uT = hv;
icvMatrAXPY3_64f( lh, l-1, hv+1, lduT, uT+1, h );
s = hv[0] * h;
for( k = 0; k < l; k++ ) hv[k] *= s;
hv[0] += 1;
}
else
{
for( j = 1; j < l; j++ )
hv[j] = 0;
for( j = 1; j < lh; j++ )
hv[j * lduT] = 0;
hv[0] = 1;
}
}
uT = u0;
}
/* accumulate right transformations */
if( vT )
{
n1 = n - n1;
vT = v0 + n1 * ldvT;
for( i = n1; i < nv; i++, vT += ldvT )
{
memset( vT + n1, 0, (n - n1) * sizeof( vT[0] ));
vT[i] = 1.;
}
for( i = n1 - 1; i >= 0; i-- )
{
double s;
int lh = nv - i;
l = n - i;
hv = v0 + (ldvT + 1) * i;
h = i == 0 ? kv0 : hv[-1];
assert( h <= 0 );
if( h != 0 )
{
vT = hv;
icvMatrAXPY3_64f( lh, l-1, hv+1, ldvT, vT+1, h );
s = hv[0] * h;
for( k = 0; k < l; k++ ) hv[k] *= s;
hv[0] += 1;
}
else
{
for( j = 1; j < l; j++ )
hv[j] = 0;
for( j = 1; j < lh; j++ )
hv[j * ldvT] = 0;
hv[0] = 1;
}
}
vT = v0;
}
for( i = 0; i < nm; i++ )
{
double tnorm = fabs( w[i] );
tnorm += fabs( e[i] );
if( anorm < tnorm )
anorm = tnorm;
}
anorm *= DBL_EPSILON;
/* diagonalization of the bidiagonal form */
for( k = nm - 1; k >= 0; k-- )
{
double z = 0;
iters = 0;
for( ;; ) /* do iterations */
{
double c, s, f, g, x, y;
int flag = 0;
/* test for splitting */
for( l = k; l >= 0; l-- )
{
if( fabs(e[l]) <= anorm )
{
flag = 1;
break;
}
assert( l > 0 );
if( fabs(w[l - 1]) <= anorm )
break;
}
if( !flag )
{
c = 0;
s = 1;
for( i = l; i <= k; i++ )
{
f = s * e[i];
e[i] *= c;
if( anorm + fabs( f ) == anorm )
break;
g = w[i];
h = pythag( f, g );
w[i] = h;
c = g / h;
s = -f / h;
if( uT )
icvGivens_64f( m, uT + lduT * (l - 1), uT + lduT * i, c, s );
}
}
z = w[k];
if( l == k || iters++ == MAX_ITERS )
break;
/* shift from bottom 2x2 minor */
x = w[l];
y = w[k - 1];
g = e[k - 1];
h = e[k];
f = 0.5 * (((g + z) / h) * ((g - z) / y) + y / h - h / y);
g = pythag( f, 1 );
if( f < 0 )
g = -g;
f = x - (z / x) * z + (h / x) * (y / (f + g) - h);
/* next QR transformation */
c = s = 1;
for( i = l + 1; i <= k; i++ )
{
g = e[i];
y = w[i];
h = s * g;
g *= c;
z = pythag( f, h );
e[i - 1] = z;
c = f / z;
s = h / z;
f = x * c + g * s;
g = -x * s + g * c;
h = y * s;
y *= c;
if( vT )
icvGivens_64f( n, vT + ldvT * (i - 1), vT + ldvT * i, c, s );
z = pythag( f, h );
w[i - 1] = z;
/* rotation can be arbitrary if z == 0 */
if( z != 0 )
{
c = f / z;
s = h / z;
}
f = c * g + s * y;
x = -s * g + c * y;
if( uT )
icvGivens_64f( m, uT + lduT * (i - 1), uT + lduT * i, c, s );
}
e[l] = 0;
e[k] = f;
w[k] = x;
} /* end of iteration loop */
if( iters > MAX_ITERS )
break;
if( z < 0 )
{
w[k] = -z;
if( vT )
{
for( j = 0; j < n; j++ )
vT[j + k * ldvT] = -vT[j + k * ldvT];
}
}
} /* end of diagonalization loop */
/* sort singular values and corresponding values */
for( i = 0; i < nm; i++ )
{
k = i;
for( j = i + 1; j < nm; j++ )
if( w[k] < w[j] )
k = j;
if( k != i )
{
double t;
CV_SWAP( w[i], w[k], t );
if( vT )
for( j = 0; j < n; j++ )
CV_SWAP( vT[j + ldvT*k], vT[j + ldvT*i], t );
if( uT )
for( j = 0; j < m; j++ )
CV_SWAP( uT[j + lduT*k], uT[j + lduT*i], t );
}
}
}
void cvSVD_my( CvArr* aarr, CvArr* warr, CvArr* uarr, CvArr* varr, int flags )
{
uchar* buffer = 0;
int local_alloc = 0;
CV_FUNCNAME( "cvSVD" );
__BEGIN__;
CvMat astub, *a = (CvMat*)aarr;
CvMat wstub, *w = (CvMat*)warr;
CvMat ustub, *u;
CvMat vstub, *v;
CvMat tmat;
uchar* tw = 0;
int type;
int a_buf_offset = 0, u_buf_offset = 0, buf_size, pix_size;
int temp_u = 0, /* temporary storage for U is needed */
t_svd; /* special case: a->rows < a->cols */
int m, n;
int w_rows, w_cols;
int u_rows = 0, u_cols = 0;
int w_is_mat = 0;
if( !CV_IS_MAT( a ))
CV_CALL( a = cvGetMat( a, &astub ));
if( !CV_IS_MAT( w ))
CV_CALL( w = cvGetMat( w, &wstub ));
if( !CV_ARE_TYPES_EQ( a, w ))
CV_ERROR( CV_StsUnmatchedFormats, "" );
if( a->rows >= a->cols )
{
m = a->rows;
n = a->cols;
w_rows = w->rows;
w_cols = w->cols;
t_svd = 0;
}
else
{
CvArr* t;
CV_SWAP( uarr, varr, t );
flags = (flags & CV_SVD_U_T ? CV_SVD_V_T : 0)|
(flags & CV_SVD_V_T ? CV_SVD_U_T : 0);
m = a->cols;
n = a->rows;
w_rows = w->cols;
w_cols = w->rows;
t_svd = 1;
}
u = (CvMat*)uarr;
v = (CvMat*)varr;
w_is_mat = w_cols > 1 && w_rows > 1;
if( !w_is_mat && CV_IS_MAT_CONT(w->type) && w_cols + w_rows - 1 == n )
tw = w->data.ptr;
if( u )
{
if( !CV_IS_MAT( u ))
CV_CALL( u = cvGetMat( u, &ustub ));
if( !(flags & CV_SVD_U_T) )
{
u_rows = u->rows;
u_cols = u->cols;
}
else
{
u_rows = u->cols;
u_cols = u->rows;
}
if( !CV_ARE_TYPES_EQ( a, u ))
CV_ERROR( CV_StsUnmatchedFormats, "" );
if( u_rows != m || (u_cols != m && u_cols != n))
CV_ERROR( CV_StsUnmatchedSizes, !t_svd ? "U matrix has unappropriate size" :
"V matrix has unappropriate size" );
temp_u = (u_rows != u_cols && !(flags & CV_SVD_U_T)) || u->data.ptr==a->data.ptr;
if( w_is_mat && u_cols != w_rows )
CV_ERROR( CV_StsUnmatchedSizes, !t_svd ? "U and W have incompatible sizes" :
"V and W have incompatible sizes" );
}
else
{
u = &ustub;
u->data.ptr = 0;
u->step = 0;
}
if( v )
{
int v_rows, v_cols;
if( !CV_IS_MAT( v ))
CV_CALL( v = cvGetMat( v, &vstub ));
if( !(flags & CV_SVD_V_T) )
{
v_rows = v->rows;
v_cols = v->cols;
}
else
{
v_rows = v->cols;
v_cols = v->rows;
}
if( !CV_ARE_TYPES_EQ( a, v ))
CV_ERROR( CV_StsUnmatchedFormats, "" );
if( v_rows != n || v_cols != n )
CV_ERROR( CV_StsUnmatchedSizes, t_svd ? "U matrix has unappropriate size" :
"V matrix has unappropriate size" );
if( w_is_mat && w_cols != v_cols )
CV_ERROR( CV_StsUnmatchedSizes, t_svd ? "U and W have incompatible sizes" :
"V and W have incompatible sizes" );
}
else
{
v = &vstub;
v->data.ptr = 0;
v->step = 0;
}
type = CV_MAT_TYPE( a->type );
pix_size = CV_ELEM_SIZE(type);
buf_size = n*2 + m;
if( !(flags & CV_SVD_MODIFY_A) )
{
a_buf_offset = buf_size;
buf_size += a->rows*a->cols;
}
if( temp_u )
{
u_buf_offset = buf_size;
buf_size += u->rows*u->cols;
}
buf_size *= pix_size;
if( buf_size <= CV_MAX_LOCAL_SIZE )
{
buffer = (uchar*)cvStackAlloc( buf_size );
local_alloc = 1;
}
else
{
CV_CALL( buffer = (uchar*)cvAlloc( buf_size ));
}
if( !(flags & CV_SVD_MODIFY_A) )
{
cvInitMatHeader( &tmat, m, n, type,
buffer + a_buf_offset*pix_size );
if( !t_svd )
cvCopy( a, &tmat );
else
cvT( a, &tmat );
a = &tmat;
}
if( temp_u )
{
cvInitMatHeader( &ustub, u_cols, u_rows, type, buffer + u_buf_offset*pix_size );
u = &ustub;
}
if( !tw )
tw = buffer + (n + m)*pix_size;
if( type == CV_32FC1 )
{
icvSVD_32f( a->data.fl, a->step/sizeof(float), a->rows, a->cols,
(float*)tw, u->data.fl, u->step/sizeof(float), u_cols,
v->data.fl, v->step/sizeof(float), (float*)buffer );
}
else if( type == CV_64FC1 )
{
icvSVD_64f( a->data.db, a->step/sizeof(double), a->rows, a->cols,
(double*)tw, u->data.db, u->step/sizeof(double), u_cols,
v->data.db, v->step/sizeof(double), (double*)buffer );
}
else
{
CV_ERROR( CV_StsUnsupportedFormat, "" );
}
if( tw != w->data.ptr )
{
int shift = w->cols != 1;
cvSetZero( w );
if( type == CV_32FC1 )
for( int i = 0; i < n; i++ )
((float*)(w->data.ptr + i*w->step))[i*shift] = ((float*)tw)[i];
else
for( int i = 0; i < n; i++ )
((double*)(w->data.ptr + i*w->step))[i*shift] = ((double*)tw)[i];
}
if( uarr )
{
if( !(flags & CV_SVD_U_T))
cvT( u, uarr );
else if( temp_u )
cvCopy( u, uarr );
/*CV_CHECK_NANS( uarr );*/
}
if( varr )
{
if( !(flags & CV_SVD_V_T))
cvT( v, varr );
/*CV_CHECK_NANS( varr );*/
}
CV_CHECK_NANS( w );
__END__;
if( buffer && !local_alloc )
cvFree( &buffer );
}
inline void cvDoubleMatPrint_my( const CvMat* mat )
{
int i, j;
for( i = 0; i < mat->rows; i++ )
{
for( j = 0; j < mat->cols; j++ )
{
printf( "%lf ",cvmGet( mat, i, j ) );
}
printf( "\n" );
}
}
void SVD_my(double *SA, double *SW, double *SU, double *SV)
{
CvMat A, W, U, V;
cvInitMatHeader(&A,4,4,CV_64FC1,SA,CV_AUTOSTEP);
cvInitMatHeader(&W,4,4,CV_64FC1,SW,CV_AUTOSTEP);
cvInitMatHeader(&U,4,4,CV_64FC1,SU,CV_AUTOSTEP);
cvInitMatHeader(&V,4,4,CV_64FC1,SV,CV_AUTOSTEP);
cvSVD_my(&A, &W, &U, &V, CV_SVD_U_T);
printf("A:============\n");
cvDoubleMatPrint_my( &A);
printf("W:============\n");
cvDoubleMatPrint_my( &W);
printf("U:============\n");
cvDoubleMatPrint_my( &U);
printf("V:============\n");
cvDoubleMatPrint_my( &V);
}
void main()
{
printf("start\n");
double A[16]={1,1,1,1,
0,1,0,0,
0,0,1,1,
0,0,0,0}, W[16], U[16], V[16];
int i,j;
SVD_my(A, W, U, V);
printf("main:V:============\n");
for(i=0;i<16;i++)
{
std::cout<<V[i]<<',';
j++;
if(j%4==0)
std::cout<<std::endl;
}
printf("stop\n");
}参考代码：#include <iostream>
#include <algorithm>
#include <string>
#include <io.h>
#include "highgui.h"
#include "cv.h"
#pragma comment(lib, "cxcore.lib")
using namespace std;
inline void cvDoubleMatPrint( const CvMat* mat );
void SVD(double *SA, double *SU, double *SS, double *SV);
int main()
{
double A[16]={1,1,1,1,0,1,0,0,0,0,1,1,0,0,0,0}, U[16], S[16], V[16];
SVD(A, U, S, V);
printf("V:============\n");
for(int j=0,i=0;i<16;i++)
{
cout<<V[i]<<',';
j++;
if(j%4==0)
cout<<endl;
}
return 0;
}
inline void cvDoubleMatPrint( const CvMat* mat )
{
int i, j;
for( i = 0; i < mat->rows; i++ )
{
for( j = 0; j < mat->cols; j++ )
{
printf( "%lf ",cvmGet( mat, i, j ) );
}
printf( "\n" );
}
}
void SVD(double *SA, double *SU, double *SS, double *SV)
{
CvMat A, U, S, V;
cvInitMatHeader(&A,4,4,CV_64FC1,SA,CV_AUTOSTEP);
cvInitMatHeader(&U,4,4,CV_64FC1,SU,CV_AUTOSTEP);
cvInitMatHeader(&S,4,4,CV_64FC1,SS,CV_AUTOSTEP);
cvInitMatHeader(&V,4,4,CV_64FC1,SV,CV_AUTOSTEP);
cvSVD(&A, &U, &S, &V, CV_SVD_U_T);
printf("A:============\n");
cvDoubleMatPrint( &A);
printf("U:============\n");
cvDoubleMatPrint( &U);
printf("S:============\n");
cvDoubleMatPrint( &S);
printf("V:============\n");
cvDoubleMatPrint( &V);
}
个别小知识介绍：对A进行SVD分解后，矩阵V的最后一列其实是齐次方程Ax=0的一组解 //输出： A:============ 1.000000 1.000000 1.000000 1.000000 0.000000 1.000000 0.000000 0.000000 0.000000 0.000000 1.000000 1.000000 0.000000 0.000000 0.000000 0.000000 U:============ 2.324366 0.000000 0.000000 0.000000 0.000000 1.147184 0.000000 0.000000 0.000000 0.000000 0.530368 0.000000 0.000000 0.000000 0.000000 0.000000 S:============ -0.846041 -0.192165 -0.497279 0.000000 -0.226091 -0.715409 0.661115 0.000000 -0.482801 0.671761 0.561818 0.000000 0.000000 0.000000 0.000000 1.000000 V:============ -0.363988 -0.197084 -0.910314 0.000000 -0.446662 -0.820705 0.356281 0.000000 -0.577930 0.379210 0.148985 0.707107 -0.577930 0.379210 0.148985 -0.707107 V:============ -0.363988,-0.197084,-0.910314,0, -0.446662,-0.820705,0.356281,0, -0.57793,0.37921,0.148985,0.707107, -0.57793,0.37921,0.148985,-0.707107, Press any key to continue //matlab: A = 1 1 1 1 0 1 0 0 0 0 1 1 0 0 0 0 >> [S V D]=svd(A) S = 0.8460 0.2261 0.4828 0 0.1922 0.7154 -0.6718 0 0.4973 -0.6611 -0.5618 0 0 0 0 1.0000 V = 2.3244 0 0 0 0 1.1472 0 0 0 0 0.5304 0 0 0 0 0 D = 0.3640 0.1971 0.9103 0 0.4467 0.8207 -0.3563 0 0.5779 -0.3792 -0.1490 -0.7071 0.5779 -0.3792 -0.1490 0.7071 */
内容来自用户分享和网络整理，不保证内容的准确性，如有侵权内容，可联系管理员处理 
相关文章推荐
- [转]尺度不变特征变换（SIFT算法）Matlab程序代码测试例子的说明(Lowe的代码）
- 模式识别 - 特征归一化 及 测试 代码(Matlab)
- Eclipse+Java+OpenCV246环境搭建和代码测试
- 测试数据精准对比的思考
- OpenCV：金字塔Lucas-Kanade光流测试代码
- Matlab 和 C 动态开辟、存储和相互文件读取的测试代码
- PHP4和PHP5性能测试和对比 测试代码与环境
- 尺度不变特征变换（SIFT算法）Matlab程序代码测试例子的说明(Lowe的代码）
- 思考--测试应该从代码中读懂什么
- matlab转化为opencv代码过程中按行扩大矩阵
- PHP4和PHP5性能测试和对比 测试代码与环境
- Eclipse+Java+OpenCV246环境搭建和代码测试
- 形态学图像处理（matlab代码测试）
- SIFT算法:Matlab程序代码测试例子的说明(Lowe的代码）
- 【转】 尺度不变特征变换（SIFT算法）Matlab程序代码测试例子说明(Lowe的代码）
- 人脸识别---光照预处理（Opencv代码和Matlab代码）
- MATLAB代码转化成C+opencv时的注意事项
- OpenCV stereo matching 代码 matlab实现视差显示
- lua：部分常用操作的效率对比及代码优化建议（附测试代码）
- Eclipse+Java+OpenCV246环境搭建和代码测试