112 lines
3.8 KiB
C++
112 lines
3.8 KiB
C++
#include <iostream>
|
|
#include "Function.h"
|
|
#include "Function2D.h"
|
|
#include "mkl.h"
|
|
|
|
using namespace Aurora;
|
|
|
|
double Aurora::immse(const Aurora::Matrix &aImageA, const Aurora::Matrix &aImageB) {
|
|
if (aImageA.getDims()!=2|| aImageB.getDims()!=2){
|
|
std::cerr<<"Fail!immse args must all 2d matrix!";
|
|
return 0.0;
|
|
}
|
|
if (!aImageB.compareShape(aImageA)){
|
|
std::cerr<<"Fail!immse args must be same shape!";
|
|
return 0.0;
|
|
}
|
|
if (aImageA.getValueType()!=Normal || aImageB.getValueType() != Normal) {
|
|
std::cerr << "Fail!immse args must be normal value type!";
|
|
return 0.0;
|
|
}
|
|
int size = aImageA.getDataSize();
|
|
auto temp = malloc(size);
|
|
vdSub(size, aImageA.getData(), aImageB.getData(), temp);
|
|
vdSqr(size, temp, temp);
|
|
double result = cblas_dasum(size, temp, 1) / (double) size;
|
|
free(temp);
|
|
return result;
|
|
}
|
|
|
|
Aurora::Matrix Aurora::inv(const Aurora::Matrix &aMatrix) {
|
|
if (aMatrix.getDims() != 2) {
|
|
std::cerr << "Fail!inv args must be 2d matrix!";
|
|
return aMatrix;
|
|
}
|
|
if (aMatrix.getDimSize(0) != aMatrix.getDimSize(1)) {
|
|
std::cerr << "Fail!inv args must be square matrix!";
|
|
return aMatrix;
|
|
}
|
|
if (aMatrix.getValueType() != Normal) {
|
|
std::cerr << "Fail!inv args must be normal value type!";
|
|
return aMatrix;
|
|
}
|
|
int size = aMatrix.getDataSize();
|
|
int *ipiv = new int[aMatrix.getDimSize(0)];
|
|
auto result = malloc(size);
|
|
cblas_dcopy(size,result, 1,aMatrix.getData(), 1);
|
|
LAPACKE_dgetrf(LAPACK_ROW_MAJOR, aMatrix.getDimSize(0), aMatrix.getDimSize(0), result, aMatrix.getDimSize(0), ipiv);
|
|
LAPACKE_dgetri(LAPACK_ROW_MAJOR, aMatrix.getDimSize(0), result, aMatrix.getDimSize(0), ipiv);
|
|
delete[] ipiv;
|
|
return Matrix::New(result,aMatrix);
|
|
}
|
|
|
|
Aurora::Matrix Aurora::inv(Aurora::Matrix&& aMatrix) {
|
|
if (aMatrix.getDims() != 2) {
|
|
std::cerr << "Fail!inv args must be 2d matrix!";
|
|
return aMatrix;
|
|
}
|
|
if (aMatrix.getDimSize(0) != aMatrix.getDimSize(1)) {
|
|
std::cerr << "Fail!inv args must be square matrix!";
|
|
return aMatrix;
|
|
}
|
|
if (aMatrix.getValueType() != Normal) {
|
|
std::cerr << "Fail!inv args must be normal value type!";
|
|
return aMatrix;
|
|
}
|
|
int *ipiv = new int[aMatrix.getDimSize(0)];
|
|
LAPACKE_dgetrf(LAPACK_ROW_MAJOR, aMatrix.getDimSize(0), aMatrix.getDimSize(0), aMatrix.getData(), aMatrix.getDimSize(0), ipiv);
|
|
LAPACKE_dgetri(LAPACK_ROW_MAJOR, aMatrix.getDimSize(0), aMatrix.getData(), aMatrix.getDimSize(0), ipiv);
|
|
delete[] ipiv;
|
|
return aMatrix;
|
|
}
|
|
|
|
Matrix Aurora::interp2(const Matrix& aX, const Matrix& aY, const Matrix& aV, const Matrix& aX1, const Matrix& aY1, InterpnMethod aMethod)
|
|
{
|
|
if (aV.getDims() != 2)
|
|
{
|
|
return Matrix();
|
|
}
|
|
|
|
if (aX1.getDimSize(0) != aY1.getDimSize(0))
|
|
{
|
|
return Matrix();
|
|
}
|
|
|
|
int columnNum = aV.getDimSize(1);
|
|
int rowNum = aV.getDimSize(0);
|
|
|
|
if(aX.getDimSize(0) != columnNum || aY.getDimSize(0) != rowNum)
|
|
{
|
|
return Matrix();
|
|
}
|
|
|
|
int nx1 = aX1.getDimSize(0);
|
|
std::shared_ptr<double> resultData = std::shared_ptr<double>(Aurora::malloc(nx1), Aurora::free);
|
|
for (int i = 0; i < nx1; ++i)
|
|
{
|
|
std::shared_ptr<double> xResultData = std::shared_ptr<double>(Aurora::malloc(columnNum), Aurora::free);
|
|
for(int j =0; j < columnNum; ++j)
|
|
{
|
|
xResultData.get()[j] = interp1(aY,aV($,j).toMatrix(),aY1(i).toMatrix(),aMethod).getData()[0];
|
|
}
|
|
Matrix xResult(xResultData,std::vector<int>{columnNum});
|
|
resultData.get()[i] = interp1(aX,xResult,aX1(i).toMatrix(),aMethod).getData()[0];
|
|
}
|
|
return Matrix(resultData,std::vector<int>{nx1});
|
|
}
|
|
|
|
Matrix Aurora::interpn(const Matrix& aX, const Matrix& aY, const Matrix& aV, const Matrix& aX1, const Matrix& aY1, InterpnMethod aMethod)
|
|
{
|
|
return Aurora::interp2(aY,aX,aV,aY1,aX1,aMethod);
|
|
}
|