logplus/DataMgr/src/DataHelper.cpp

1051 lines
25 KiB
C++
Raw Normal View History

2025-10-30 09:50:15 +08:00
#include "DataHelper.h"
#include <QWidget>
#include "cdialog.h"
#include "GeometryUtils.h"
// #include "DataImport.h"
// #include "ObjProject.h"
#include "BaseFun.h"
//DEFAULTTABLE DefauleTable[DefTabNum];
DEFAULTTABLE *DefauleTable=NULL;
int DefTabNum=0;
bool DataHelper::m_IsGc=true;//false;
DataHelper::DataHelper()
{
}
DataHelper::~DataHelper()
{
}
//whp change 2020.3.9 for 统一离散数据入口
//QStringList DataHelper::GetVlidTable(QStringList listFiles)
QStringList DataHelper::GetVlidTable(int curitemType,QStringList listFiles,QStringList &BadlistFiles)
{
int BadFileNum=0;
QStringList myListFile,BadListFile;
QString mes="";
char *pLine=new char[50000];
foreach(QString strFileName, listFiles)
{
QFile file(strFileName);
if (!file.open(QIODevice::ReadOnly )){mes+=QString::fromLocal8Bit("无法打开文件");mes+=strFileName;mes+="\r\n";continue;}
file.readLine(pLine,50000);
QByteArray line=pLine;
line=line.simplified();
while(IsNullLine(line)){
file.readLine(pLine,50000);
line=pLine;
}
QStringList StrList,StrListZD=GetStringList(line);//line.split(",");
int NumLine=0;
int ColNum=StrListZD.size();//列数
int InvalidLineNum=0;//无效行,即数据列数和首行不一致
// QString mes1="\r\n"+strFileName+QString::fromLocal8Bit("数据表中各行字段数不一致\r\n首行字段数=")+QString::number(ColNum);
QString mes1="\r\n"+QString::number(BadFileNum+1)+""+strFileName+"\r\n"+QString::fromLocal8Bit("数据行字段数大于字段名个数\r\n首行字段数=")+QString::number(ColNum);
//if(BadFileNum==10)mes1="............";
if(ColNum<=1)mes1="\r\n"+QString::number(BadFileNum+1)+""+strFileName+"\r\n"+QString::fromLocal8Bit("首行字段数=")+QString::number(ColNum);
if(BadFileNum>5)mes1="";
//whp add 3.13 for 多井时检测井名有效性
bool IsDwell=0;
if(StrListZD.at(0).toUpper().indexOf("WELLNAME")>=0)IsDwell=1;
if(curitemType&&IsDwell)//whp 2020.3.23 for 选择单井导入时不导入多井数据
{
BadFileNum++;
BadListFile.append(strFileName);
file.close();
continue;
}
while(!file.atEnd())
{
file.readLine(pLine,50000);
QByteArray line=pLine;
line=line.simplified();
if(line=="")continue;
if(IsNullLine(line))continue;
StrList=GetStringList(line);
if(InvalidLineNum<10)
{
if(StrList.size()>ColNum)//if(StrList.size()!=ColNum)允许缺失后面数据,不允许数据个数大于字段名个数
{
InvalidLineNum++;
if(BadFileNum<=5)mes1=mes1+QString::fromLocal8Bit("\r\n")+QString::number(NumLine+1)+QString::fromLocal8Bit("行字段数=")+QString::number(StrList.size());
}
}
//whp add 3.13 for 多井时检测井名有效性
if(!InvalidLineNum)
{
if(IsDwell)
{
if(IsValidWellName(StrList.at(0))==0)
{
QString ss=StrList.at(0);
InvalidLineNum++;
}
}
}
NumLine++;
if(InvalidLineNum)break;
}
if(InvalidLineNum>0||ColNum<=1)
{
mes+=mes1;//无效的表数据,不再显示
BadFileNum++;
BadListFile.append(strFileName);
}
else myListFile.append(strFileName);
file.close();
}
delete pLine;
BadlistFiles=BadListFile;//whp add 2020.3.9 for 统一离散数据入口
return myListFile;
}
int DataHelper::StrType(QString str)
{
QByteArray ba = str.toLatin1();//QString 转换为 char*
const char *buf = ba.data();
int len=strlen(str.toStdString().c_str());//str.length();//whp change 2020.3.16
int dig[30];
int DotNum=0;//小数点个数
BOOL IsDigital=1;
for(int i=0;i<len;i++)
{
if(i==0&&buf[i]=='-')continue;//可能是负数
if(buf[i]=='.')DotNum++;
else if(!isdigit(buf[i])){IsDigital=0;break;}
if(DotNum>1){IsDigital=0;break;}
}
if(!IsDigital)return 6;//字符串
else if(DotNum)return 4;//浮点数
else return 1;//整形数
}
DEFAULTTABLE DataHelper::SetTableInf(QString TableName,QString TableHzName,int num,QStringList ZdName,QStringList HzName,QStringList ZdUnit,QStringList ZdType,QStringList ZdLen,QStringList ZdRes)
{
DEFAULTTABLE DefTable;
DefTable.TableName=TableName;
DefTable.TableAliasName=TableHzName;
DefTable.ZdNum=num;
DefTable.tinfo=new Slf_OBJECT_FIELD[num];
for(int i=0;i<num;i++)
{
strcpy(DefTable.tinfo[i].Name,ZdName[i].toStdString().c_str());
strcpy(DefTable.tinfo[i].HZName,HzName[i].toStdString().c_str());
strcpy(DefTable.tinfo[i].Unit,ZdUnit[i].toStdString().c_str());
strcpy(DefTable.tinfo[i].HZUnit,"");
DefTable.tinfo[i].RepCode=ZdType[i].toInt();
DefTable.tinfo[i].CodeLength=ZdLen[i].toInt();
DefTable.tinfo[i].Sequence=0;//2是否连续控制
DefTable.tinfo[i].Start=1;//4字段起始值
DefTable.tinfo[i].End=99999;//4字段起始值
DefTable.tinfo[i].Rlev=1;//4字段采样间隔
DefTable.tinfo[i].SamplePoint=1;//4一个阵列的横向采样样本点数
DefTable.tinfo[i].ArrayNum=1;//4阵列数
DefTable.tinfo[i].Vmin=-99999;//4字段最大值
DefTable.tinfo[i].Vmax=99999;//4字段最小值
DefTable.tinfo[i].DefVal=9999;//4缺省值
//*(DWORD *)&DefTable.tinfo[i].Reserved=ZdRes[i].toInt();
//20201126 GZL change
strcpy(DefTable.tinfo[i].Reserved, ZdRes[i].toStdString().c_str());
}
return DefTable;
}
void DataHelper::InitDefauleTable()
{
if(!m_IsGc)return;
m_IsGc=false;
QString ConfigName = ::GetConfPath()+"DefauleTable.ini";
QFile InFile(ConfigName);
if( !InFile.open(QIODevice::ReadOnly ) )
{
QMessageBox::warning(NULL,"提示","打开缺省表配置文件"+ConfigName+"错误");
return ;
}
QTextStream textstream(&InFile);
DefTabNum=0;
QString line;
while(!textstream.atEnd())
{
line=textstream.readLine(5000);
if(line.indexOf("//")>=0)continue;
if(line.indexOf("#")>=0)
{
DefTabNum++;
}
}
if(!DefTabNum)
{
QMessageBox::warning(NULL,"提示","缺省表配置文件"+ConfigName+"中定义表数为0");
return ;
}
DefauleTable=new DEFAULTTABLE[DefTabNum+1];
InFile.seek(0);
QString TableName,TableAliasName;
//int ZdNum;
QStringList ZdName,HzName,ZdUnit,ZdLen,ZdType,ZdRes,DataList;
DefTabNum=0;
line=textstream.readLine(5000);
// line = line.toUpper();
while(!textstream.atEnd())
{
line=textstream.readLine(5000);
while(line.indexOf("#")>=0)
{
line=line.replace("#","");
DataList=GetStringList(line,1,1,1,1,1);
if(DataList.count()<2)continue;
TableName=DataList.at(0);
TableAliasName=DataList.at(1);
while(!textstream.atEnd())
{
line=textstream.readLine(5000);
if(line.indexOf(",")>=0)
{
continue;//枚举行
}
// line=line.toUpper();
if(line.indexOf("#")>=0)
{
break;
}
DataList=GetStringList(line,1,1,1,1,1);
if(DataList.count()<5)continue;
ZdName.append(DataList.at(0).toUpper());
HzName.append(DataList.at(1).toUpper());
ZdUnit.append(DataList.at(2));
bool IsZdType=0;
for(int j=0;j<11;j++)
{
if(DataList.at(3).toUpper()==QString(Rep_STR[j]))
{
ZdType.append(QString::number(j+1));
IsZdType=1;
break;
}
}
if(!IsZdType)QMessageBox::warning(NULL,"提示","defauletable.ini文件中表"+TableName+"的字段名"+DataList.at(0)+"类型不正确");
ZdLen.append(DataList.at(4));
if(DataList.count()>5){
ZdRes.append(DataList.at(5).toUpper());
}
else
ZdRes.append(DataList.at(3).toUpper());
}
DefauleTable[DefTabNum++]=SetTableInf(TableName,TableAliasName,ZdName.count(),ZdName,HzName,ZdUnit,ZdType,ZdLen,ZdRes);
ZdName.clear();
ZdLen.clear();
ZdType.clear();
HzName.clear();
HzName.clear();
ZdUnit.clear();
ZdRes.clear();
break;
}
}
InFile.close();
}
void DataHelper::DelDefauleTable()
{
if(!m_IsGc)
{
for(int i=0;i<DefTabNum;i++)
{
delete []DefauleTable[i].tinfo;
}
delete []DefauleTable;
DefTabNum=0;
}
m_IsGc=true;
}
int DataHelper::GetSimilarTable(int ZdNum,int *FieldType)
{
int DefaultNo=-1;
//先找字段类型全部一致的进行匹配
for(int i=0;i<DefTabNum;i++)
{
if(ZdNum!=DefauleTable[i].ZdNum)continue;//字段数不同的不匹配
int SameTypeNum=0;
for(int j=0;j<ZdNum;j++)
{
if(DefauleTable[i].tinfo[j].RepCode==FieldType[j])SameTypeNum++;
else break;
}
if(SameTypeNum==ZdNum)return i;
}
//whp del 2020.5.25 会把其他四字段的数据匹配成RESULT不合适
/* //如果字段类型不一致只要可以转换成缺省表的字段类型也可以比如int、float都可转换成stringint可以转换成float
for(int i=0;i<DefTabNum;i++)
{
if(ZdNum!=DefauleTable[i].ZdNum)continue;//字段数不同的不匹配
int SameTypeNum=0;
for(int j=0;j<ZdNum;j++)
{
if(DefauleTable[i].tinfo[j].RepCode>=FieldType[j])SameTypeNum++;
}
if(SameTypeNum==ZdNum)return i;//先找字段类型全部一致,匹配
}*/
return DefaultNo;
}
//剔除空行
bool DataHelper::IsNullLine(QByteArray line)
{
if(line=="")return 1;
QStringList List=GetStringList(line);
QString sss;
bool IsNullLine=1;
for(int kk=0;kk<List.size();kk++)
{
sss=List.at(kk);
if(sss!="")
{
IsNullLine=0;break;
}
}
return IsNullLine;
}
QStringList DataHelper::GetAllWellName()
{
QStringList WellNameList;
// WellNameList.clear();
// PaiObject *pPrj=GetProject();
// if(!pPrj) return WellNameList;
// PaiObject *pObj=pPrj->GetObjectByType(GetClassID_WellsFolder());
// QList<PaiObject *> childs;
// pObj->GetChildren(childs);
// int wellnum=childs.count();
// foreach(PaiObject *pChild,childs)
// {
// CObjWell *pWell=dynamic_cast<CObjWell*>(pChild);
// if(pWell) WellNameList.append(pWell->GetName());
// }
return WellNameList;
}
QStringList DataHelper::GetAllWellRoundName(QString WellName,QStringList& WellRoundFileNameList)
{
QStringList WellRoundNameList;
WellRoundNameList.clear();
WellRoundFileNameList.clear();
// CObjWell * pWell = CDataImport::GetWellByName(WellName);
// QList<PaiObject*> wellroundchildren;
// int count=pWell->GetAllWellRound(wellroundchildren);
// for(int i=0;i<count;i++) {
// CObjWelllogRound *pWR=dynamic_cast<CObjWelllogRound *>(wellroundchildren[i]);
// QString strfileWellName=pWR->GetSlfFileName();
// WellRoundFileNameList.append(pWR->GetSlfFileName());
// QFileInfo fileInfo(pWR->GetSlfFileName());
// QString strWellName = fileInfo.completeBaseName();
// WellRoundNameList.append(strWellName);
// }
return WellRoundNameList;
}
int Function::get_exp(float *result)
{
*result=0.0;
get_token();
if (!*token)
{
//serror(2);
return 2+1;
}
int re=level1(result);
if(re>0)return re;
else return 0;
}
int Function::level1(float *result)
{
int slot,ttok_type;
char temp_token[80];
int re;
if (tok_type==VARIABLE) /* save old token */
{
strcpy(temp_token,token);
ttok_type=tok_type;
slot=toupper(*token)-'A';
get_token();
if (*token !='=')
{
putback(); /* return current token */
/* restore old token */
strcpy(token,temp_token);
tok_type=ttok_type;
}
else
{
get_token(); /* get next part of exp */
re=level2(result);
if(re>0)return re;
vars[slot]=*result;
return 0;
}
}
re=level2(result);
if(re>0)return re;
else return 0;
}
/* add or subtract two terms */
int Function::level2(float *result)
{
register char op;
float hold;
int re=level3(result);
if(re>0)return re;
while ((op=*token)=='+' || op=='-')
{
get_token();
re=level3(&hold);
if(re>0)return re;
arith(op,result,&hold);
}
return 0;
}
/* multiply or divide two factors */
int Function::level3(float *result)
{
register char op;
float hold;
int re=level4(result);
if(re>0)return re;
while ((op=*token)=='*' || op=='/' || op=='%')
{
get_token();
re=level4(&hold);
if(re>0)return re;
arith(op,result,&hold);
}
return 0;
}
/* process integer exponet */
int Function::level4(float *result)
{
register char op;
float hold;
int re=level5(result);
if(re>0)return re;
while ((op=*token)=='!' || op=='^')
{
get_token();
re=level4(&hold);
if(re>0)return re;
arith(op,result,&hold);
}
return 0;
}
/* unary + or - */
int Function::level5(float *result)
{
register char op;
op=0;
if ((tok_type==DELIMITER) && *token=='+' || *token=='-')
{
op=*token;
get_token();
}
int re=level6(result);
if(re>0)return re;
if (op)
unary(op,result);
return 0;
}
/* process parenthesized expression */
int Function::level6(float *result)
{int re;
if ((*token=='(') && (tok_type==DELIMITER))
{
get_token();
re=level1(result);
if(re>0)return re;
if (*token !=')')
{
return 1+1;
//serror(1);
}
get_token();
}
else
re=primitive(result);
return re;
}
/* find value of number or variable */
int Function::primitive(float *result)
{
int re;
switch (tok_type)
{
case VARIABLE:
*result=find_var(token,&re);
if(re>0)return re;
get_token();
return 0;
case MNUMBER:
*result=atof(token);
get_token();
return 0;
default:
//serror(0);
return 1;
}
}
/* perform the specified arithmetic */
void Function::arith(char o,float *r,float *h)
{
register int t,ex;
switch (o)
{
case '-':
*r=*r-*h;
break;
case '+':
*r=*r+*h;
break;
case '*':
*r=(*r) * (*h);
break;
case '/':
if(*h==0) *r=0;
else *r=(*r)/(*h);
break;
case '%':
if(*h==0) *r=0;
else {
t=(*r)/(*h);
*r=*r-t*(*h);
}
break;
case '^':
ex=*r;
if (*h==0)
{
*r=1;
break;
}
//for (t=*h-1;t>0;--t) *r=(*r) *ex;
*r=pow((*r),(*h));
break;
case '!':
if (*r==10.)
*r=log10(*h);
else
*r=log(*h);
/*
if (*h==10.)
*r=log10(*r);
else
*r=log(*r);
*/
break;
}
}
void Function::unary(char o,float *r)
{
if (o=='-') *r=-(*r);
}
/* return a token to its resting place */
void Function::putback()
{
char *t;
t=token;
for (;*t;t++) prog--;
}
/* find the value of a variable */
float Function::find_var(char *s,int *ret)
{
*ret=0;
if (!isalpha(*s))
{
//serror(1);
*ret=2;
return 0.;
}
float a;
// return d_pDoc->INP[0].yy[toupper(*token)-'A'];
return val[toupper(*token)-'A'];
}
/* display an error message */
void Function::serror(int error)
{
const static char *e[]={
"syntax error.",
"unbalanced parentheses",
"no expression present"
};
AfxMessageBox(e[error]);
}
/* get a token */
void Function::get_token()
{
register char *temp;
tok_type=0;
temp=token;
while (iswhite(*prog)) ++prog; /* skip over white space */
char strTmp[11] = "+-*/%!^=()";
if (is_in(*prog,strTmp))
{
tok_type=DELIMITER;
*temp++=*prog++;
/* advance to next position */
}
else
if (isalpha(*prog))
{
while (!isdelim(*prog))
{
*temp++=*prog++;
}
tok_type=VARIABLE;
}
else
if (isdigit(*prog))
{
while (!isdelim(*prog))
{
*temp++=*prog++;
}
tok_type=MNUMBER;
}
*temp='\0';
}
int Function::iswhite(char c)
{
/* look for space and tabs */
if (c==' ' || c==9) return 1;
return 0;
}
int Function::isdelim(char c)
{
char strTmp[11] = "+-*/%!^=()";
if (is_in(c,strTmp) || c==9 || c=='\r' || c==0)
return 1;
return 0;
}
int Function::is_in(char ch,char *s)
{
while (*s) if (*s++==ch) return 1;
return 0;
}
int Function::GetExpress(QString csExpress)
{
int len,dPos[MaxArg],dCurvePos[MaxArg],dCurveNum,dOutPos;
int i,j,flag,num,noi;
int x1,y1,y2;
QString csMessage,cs;
char name[curve_name_len+1];
int iProgramOutCurvePos;
Slf_OBJECT_ENTRY entryy;
Slf_CHANNEL channel;
// notice :
// 进行曲线计算时,作为变量的曲线采样间距应一致
//A-Z 26 user's variable
for (i=0; i<MaxArg; i++)
{
vars[i] = 0;
csName[i]="";
}
csExpress.replace(" ","");//清除tab键
csExpress.replace(" ","");//清除空格
csExpress.toUpper();
len = csExpress.length();
num = csExpress.indexOf("=");
if ( num <= 0)//检查“=”
{
QMessageBox::warning(NULL,"提示","对不起,您输入的表达式\n"+csExpress+"不正确\n请重新输入");
return -1;
}
QString Model=csExpress;
csOutName = csExpress.left(num);
csExpress = csExpress.right(len-num-1);
int Lkh=0,Rkh=0;//左右括号个数
char cExpress[200];
strcpy(cExpress,csExpress.toStdString().c_str());
for(int i=0;i<csExpress.length();i++)//检查不识别的字符
{
if(cExpress[i]<0)//检查汉字
{
QString ss="对不起,该系统不支持汉字表达式计算\r\n所以无法处理您所输入的公式:"+Model+"\r\n请重新输入";
QMessageBox::warning(NULL,"提示",ss);
return -1;
}
char strTmp[11] = "+-*/%!^=()";
if ( !(isupper(cExpress[i]) || isdigit(cExpress[i]) || cExpress[i]=='.'|| cExpress[i]=='_' ||is_in(cExpress[i],strTmp)))
{
QString ss="对不起,您所输入的表达式\r\n"+Model+"\r\n中含有该系统无法识别的字符";//+QString(QLatin1String(cExpress[i]))+"\r\n请重新输入";
QMessageBox::warning(NULL,"提示",ss);
return -1;
}
if(cExpress[i]=='(')Lkh++;
else if(cExpress[i]==')')Rkh++;
}
if(Lkh!=Rkh)
{
QMessageBox::warning(NULL,"提示","计算公式:"+Model+" 中含有不对称的括号\r\n无法进行计算,请重新输入");
return -1;
}
len = csExpress.length();
if(csOutName.length()<0||len<=0)
{
QMessageBox::warning(NULL,"提示","对不起,您输入的表达式\n"+csExpress+"不正确\n请重新输入,如AC1=AC+20\r\n");
return -1;
}
strcpy(name,csOutName.toStdString().c_str());
x1 = 0; y1 = 0;
len = csExpress.length();
while ( x1 < len )
{
flag = 0;
for (i=x1; i<len; i++)
{
if ( isupper(cExpress[i])) //csExpress[i]<0汉字
{
dPos[y1] = i;
flag=1;
break;
}
}
x1=i;
if ( flag )
{
y2 = 0;
for (i=dPos[y1]; i<len; i++)
{
if ( isupper(cExpress[i]) || isdigit(cExpress[i])
|| cExpress[i]=='.'|| cExpress[i]=='_' )
{
csName[y1]+=cExpress[i];
}
else
break;
}
x1 = i;
y1 ++;
}
}
noi=y1;
for (i=0; i<noi; i++)
{
cExpress[dPos[i]]='A'+ i;
len = csName[i].length();
for (j=1; j<len; j++)cExpress[dPos[i]+j]=' ';
/*csExpress.setAt(dPos[i],'A'+ i);
len = csName[i].GetLength();
for (j=1; j<len; j++)
csExpress.SetAt(dPos[i]+j,(TCHAR)' ');*/
}
strcpy(TempProg,cExpress);//csExpress);
return noi;
}
//Resample
float Resample::GetXstep(float x2,float x1)
{
float xstep;
double eps=1.0e-5;
xstep=x2-x1;
if(fabs(xstep)<eps)
{
if(xstep>=0.0) xstep=eps; else xstep=-eps;
}
return(xstep);
}
//---------------------------------------------------------------------
// Input: (x1,y1),(x2,y2),...,(xk,yk),(xk+1,yk+1),...,(xn,yn),
// Region: {x[k],x[k+1]}
// Function: Sk(x)=a[k]+b[k]*(x-x[k])
// Output: a[k],b[k]
//---------------------------------------------------------------------
void Resample::Linear(float *x,float *y,int n,float *a,float *b)
{
for(int i=0;i<n-1;i++)
{
b[i]=(y[i+1]-y[i])/GetXstep(x[i+1],x[i]);
a[i]=y[i];
}
a[n-1]=a[n-2]; b[n-1]=0.0;
}
//---------------------------------------------------------------------
// Input: (x1,y1),(x2,y2),...,(xk,yk),(xk+1,yk+1),...,(xn,yn),
// Region: {x[k],x[k+1]}
// Function: Sk(x)=a[k]+b[k]*(x-x[k])+c[k]*(x-x[k])^2
// Output: a[k],b[k],c[k]
//---------------------------------------------------------------------
void Resample::Parabola(float *x,float *y,int n,float *a,float *b,float *c)
{
float *dx1,*dx2;
dx1=new float [n]; dx2=new float [n];
if(dx1==NULL||dx2==NULL) return;
for(int i=0;i<n-1;i++) dx1[i]=GetXstep(x[i+1],x[i]);
for(int i=0;i<n-2;i++) dx2[i]=GetXstep(x[i+2],x[i]);
dx2[n-2]=GetXstep(x[n-1],x[n-3]);
for(int i=0;i<n-2;i++)
{
a[i]=y[i];
b[i]=-y[i]*(dx1[i]+dx2[i])/dx1[i]/dx2[i]+y[i+1]*dx2[i]/dx1[i]/dx1[i+1]
-y[i+2]*dx1[i]/dx1[i+1]/dx2[i];
c[i]=y[i]/dx1[i]/dx2[i]-y[i+1]/dx1[i]/dx1[i+1]+y[i+2]/dx1[i+1]/dx2[i];
}
a[n-2]=y[n-2];
b[n-2]=-y[n-3]*dx1[n-2]/dx1[n-3]/dx2[n-2]+y[n-2]*(dx1[n-2]-dx1[n-3])/dx1[n-3]/dx1[n-2]
+y[n-1]*dx1[n-3]/dx1[n-2]/dx2[n-2];
c[n-2]=y[n-3]/dx1[n-3]/dx2[n-2]-y[n-2]/dx1[n-3]/dx1[n-2]+y[n-1]/dx1[n-2]/dx2[n-2];
a[n-1]=y[n-1]; b[n-1]=0.0; c[n-1]=0.0;
delete [] dx1;
delete [] dx2;
}
//---------------------------------------------------------------------
// Input: (x1,y1),(x2,y2),...,(xk,yk),(xk+1,yk+1),...,(xn,yn),
// Region: {x[k],x[k+1]}
// Function: Sk(x)=a[k]+b[k]*(x-x[k])+c[k]*(x-x[k])^2+d[k]*(x-x[k])^3
// Output: a[k],b[k],c[k],d[k]
//---------------------------------------------------------------------
void Resample::Akima(float *x,float *y,int n,float *a,float *b,float *c,float *d)
{
float *u,*v;
double eps=1.0e-6;
int dtemp=2;
u=new float [n+dtemp+1]; v=new float [n+dtemp+1];
if(u==NULL||v==NULL) return;
n=n-1;
for(int i=0;i<n;i++) u[i+dtemp]=(y[i+1]-y[i])/GetXstep(x[i+1],x[i]);
u[-1+dtemp]=2*u[0+dtemp]-u[1+dtemp];
u[-2+dtemp]=2*u[-1+dtemp]-u[0+dtemp];
u[n+dtemp]=2*u[n-1+dtemp]-u[n-2+dtemp];
u[n+1+dtemp]=2*u[n+dtemp]-u[n-1+dtemp];
for(int i=0;i<n+1;i++)
{ if(fabs(u[i+1+dtemp]-u[i+dtemp])+fabs(u[i-1+dtemp]-u[i-2+dtemp])<eps)
v[i]=(u[i-1+dtemp]+u[i+dtemp])/2.0;
else
v[i]=(u[i-1+dtemp]*fabs(u[i+1+dtemp]-u[i+dtemp])
+u[i+dtemp]*fabs(u[i-1+dtemp]-u[i-2+dtemp]))
/(fabs(u[i+1+dtemp]-u[i+dtemp])+fabs(u[i-1+dtemp]-u[i-2+dtemp])+1.0e-5);
}
float xstep;
for(int i=0;i<n;i++)
{ a[i]=y[i];
b[i]=v[i];
xstep=GetXstep(x[i+1],x[i]);
c[i]=(3*u[i+dtemp]-2*v[i]-v[i+1])/xstep;
d[i]=(v[i+1]+v[i]-2*u[i+dtemp])/xstep/xstep;
}
a[n]=y[n]; b[n]=0.0; c[n]=0.0; d[n]=0.0;
delete [] u;
delete [] v;
}
//--------------------------------------------------------------------
//
bool Resample::LReSampling(float *aDep,float *aVal,int anp,float *tDep,float *curve,int tnp,int method)
{
int id;
float av;
float xk,x,*a,*b,*c,*d;
a=new float [anp+1]; if(a==NULL) return(false);
b=new float [anp+1]; if(b==NULL) { delete [] a; return(false); }
c=new float [anp+1]; if(c==NULL) { delete [] a; delete [] b; return(false); }
d=new float [anp+1]; if(d==NULL) { delete [] a; delete [] b; delete [] c; return(false); }
av=0.0; for(int k=0; k<anp; k++) av+=(aVal[k]/anp);
// Calculate coefficents
switch(method)
{
case 0:
Linear(aDep,aVal,anp,a,b);
break;
case 1:
Parabola(aDep,aVal,anp,a,b,c);
break;
case 2:
Akima(aDep,aVal,anp,a,b,c,d);
break;
}
// Resampling
for(int kk=0; kk<tnp; kk++)
{
x=tDep[kk];
if(x<aDep[0]) { id=0; }
else if(x>aDep[anp-1]){ id=anp-1; }
else
{
id=0; for(int i=0;i<anp-1;i++) if((aDep[i]<=x)&&(x<aDep[i+1])){ id=i; break; }
}
xk=aDep[id];
switch(method)
{
case 0:
curve[kk]=a[id]+b[id]*(x-xk);
break;
case 1:
curve[kk]=a[id]+b[id]*(x-xk)+c[id]*(x-xk)*(x-xk);
break;
case 2:
curve[kk]=a[id]+b[id]*(x-xk)+c[id]*(x-xk)*(x-xk)+d[id]*(x-xk)*(x-xk)*(x-xk);
break;
}
if(curve[kk]<0.0) curve[kk]=av;
}
curve[tnp-1]=curve[tnp-2];
delete [] a;
delete [] b;
delete [] c;
delete [] d;
return(true);
}
void Resample::qs(float*datx,float*mr,int left,int right,int Col)
{
int i,j;
float x,y,va[1024];
i=left;j=right;
x=datx[(left+right)/2];
do {
while(datx[i]<x&&i<right) i++;
while(datx[j]>x&&j>left ) j--;
if(i<=j) {
y=datx[i];
datx[i]=datx[j];
datx[j]=y;
memmove(va,&mr[i*Col],Col*sizeof(float));
memmove(&mr[i*Col],&mr[j*Col],Col*sizeof(float));
memmove(&mr[j*Col],va,Col*sizeof(float));
i++;j--;
}
} while(i<=j);
if(left<j) qs(datx,mr,left,j,Col);
if(i<right) qs(datx,mr,i,right,Col);
}
// ww chanage
void Resample::qs(float*datx,short*mr,int left,int right,int Col)
{
int i,j;
float x,y,va[1024];
i=left;j=right;
x=datx[(left+right)/2];
do {
while(datx[i]<x&&i<right) i++;
while(datx[j]>x&&j>left ) j--;
if(i<=j) {
y=datx[i];
datx[i]=datx[j];
datx[j]=y;
memmove(va,&mr[i*Col],Col*sizeof(short));
memmove(&mr[i*Col],&mr[j*Col],Col*sizeof(short));
memmove(&mr[j*Col],va,Col*sizeof(short));
i++;j--;
}
} while(i<=j);
if(left<j) qs(datx,mr,left,j,Col);
if(i<right) qs(datx,mr,i,right,Col);
}
void Resample::ReOrder(float *mr,int Row,int Col,float *tDep)
{
// ww change
qs(tDep,mr,0,Row-1,Col);
/*
float vd;
float va[1024];
for(int i=0; i<Row-1; i++)
{
for(int k=i+1; k<Row; k++)
{
if(tDep[k]<tDep[i])
{
for(int j=0; j<Col; j++) va[j]=mr[i*Col+j];
for(int j=0; j<Col; j++) mr[i*Col+j]=mr[k*Col+j];
for(int j=0; j<Col; j++) mr[k*Col+j]=va[j];
vd=tDep[i]; tDep[i]=tDep[k]; tDep[k]=vd;
}
}
}
*/
}
void Resample::ReOrder(short *mr,int Row,int Col,float *tDep)
{
// ww chanage
qs(tDep,mr,0,Row-1,Col);
/*
float vd;
short va[1024];
for(int i=0; i<Row-1; i++)
{
for(int k=i+1; k<Row; k++)
{
if(tDep[k]<tDep[i])
{
for(int j=0; j<Col; j++) va[j]=mr[i*Col+j];
for(int j=0; j<Col; j++) mr[i*Col+j]=mr[k*Col+j];
for(int j=0; j<Col; j++) mr[k*Col+j]=va[j];
vd=tDep[i];
tDep[i]=tDep[k];
tDep[k]=vd;
}
}
}
*/
}
void Resample::ReSampling(float *mr,int Row,int Col,float *aDep,float *tDep,int tnp)
{
int nnn;
float *lbuf,*aVal;
nnn=tnp; if(nnn<Row) nnn=Row;
lbuf=new float [nnn+10]; if(lbuf<0) return;
aVal=new float [nnn+10]; if(aVal<0){ delete [] lbuf; return; }
for(int k=0; k<Row; k++) aVal[k]=tDep[k];
ReOrder(mr,Row,Col,aDep);
//ww change
for(int j=0; j<Col; j++)
{
for(int k=0; k<Row; k++) lbuf[k]=mr[k*Col+j];
LReSampling(aDep,lbuf,Row,tDep,aVal,tnp,0);
for(int k=0; k<tnp; k++) mr[k*Col+j]=aVal[k];
}
delete [] lbuf; delete [] aVal;
}
void Resample::ReSampling(short *mr,int Row,int Col,float *aDep,float *tDep,int tnp)
{
int nnn;
float *lbuf,*aVal;
nnn=tnp; if(nnn<Row) nnn=Row;
lbuf=new float [nnn+10]; if(lbuf<0) return;
aVal=new float [nnn+10]; if(aVal<0){ delete [] lbuf; return; }
ReOrder(mr,Row,Col,tDep);
//ww chanage
for(int j=0; j<Col; j++)
{
for(int k=0; k<Row; k++) lbuf[k]=mr[k*Col+j];
LReSampling(aDep,lbuf,Row,tDep,aVal,tnp,0);
for(int k=0; k<tnp; k++) mr[k*Col+j]=aVal[k];
}
delete [] lbuf; delete [] aVal;
}