logplus/logPlus/DepPairs.cpp

1236 lines
33 KiB
C++
Raw Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

//////////////////////////////////////////////////////////////////
///
/// 测录井信息综合平台
///
/// 编制:王卫
///
/// 版权:中国石化集团石油工程技术研究院
///
/// 2010年10月-
///
// Depth.cpp : implementation file
//
#include "DepPairs.h"
#include "MemRdWt.h"
#include "math.h"
#include "BaseFun.h"
#include "geometryutils.h"
#include <qdir.h>
#include "qapplication.h"
#include <QMessageBox.h>
int GetDepPairs(float* curve1,float stdep1,float endep1,float rlev1,
float* curve2,float stdep2,float endep2,float rlev2,float cmins,
float winLength,float stepLength, float searchLength,
float *FirstDep,
float *SecondDep,
float *corc)
{
DepPairs depPairs(curve1,stdep1,endep1,rlev1,
curve2,stdep2,endep2,rlev2,cmins,
winLength,stepLength,searchLength);
depPairs.FirstDep=FirstDep;
depPairs.SecondDep=SecondDep;
depPairs.corc=corc;
float stdep=FirstDep[0];
float endep=SecondDep[0];
if(stdep==endep){
stdep=stdep1;
endep=endep1;
stdep=stdep1;
if(stdep>stdep2) stdep=stdep2;
endep=endep1;
if(endep<endep2) endep=endep2;
}
depPairs.StartDepth=stdep;
depPairs.EndDepth=endep;
return depPairs.AdmRun(curve1,curve2,NULL);
}
int getDepPairs(float* curve1,float stdep1,float endep1,float rlev1,
float* curve2,float stdep2,float endep2,float rlev2,
float *FirstDep,
float *SecondDep)
{
DepPairs depPairs(curve1,stdep1,endep1,rlev1,
curve2,stdep2,endep2,rlev2,NULL,
0,0,0);
depPairs.FirstDep=FirstDep;
depPairs.SecondDep=SecondDep;
depPairs.corc=NULL;
float stdep=FirstDep[0];
float endep=SecondDep[0];
if(stdep==endep){
stdep=stdep1;
endep=endep1;
stdep=stdep1;
if(stdep>stdep2) stdep=stdep2;
endep=endep1;
if(endep<endep2) endep=endep2;
}
depPairs.StartDepth=stdep;
depPairs.EndDepth=endep;
depPairs.ipar=1;
return depPairs.AdmRun(curve1,curve2,NULL);
}
int AutoCompZone(char *Filename,char *Filename1,char *CurveName1,char *CurveName2,
float *FirstDep,
float *SecondDep)
{
CMemRdWt well;
if(!well.Open(Filename,CMemRdWt::modeRead)) {
return 0;
};
CMemRdWt well1;
if(!well1.Open(Filename1,CMemRdWt::modeRead)) {
return 0;
}
int index=well.OpenCurve(CurveName1);
int index1=well1.OpenCurve(CurveName2);
if(index>-1&&index1>-1) {
Slf_CURVE curve,curve1;
well.GetCurveInfo(index,&curve);
int count=(curve.EndDepth-curve.StartDepth)/curve.DepLevel+1.5;
float *buffer=new float[count+1];
well.ReadCurveToFloatBuf(index,curve.StartDepth,count,buffer);
well1.GetCurveInfo(index1,&curve1);
int count1=(curve1.EndDepth-curve1.StartDepth)/curve1.DepLevel+1.5;
float *buffer1=new float[count1+1];
well1.ReadCurveToFloatBuf(index1,curve1.StartDepth,count1,buffer1);
int num=getDepPairs(buffer,curve.StartDepth,curve.EndDepth,curve.DepLevel,
buffer1,curve1.StartDepth,curve1.EndDepth,curve1.DepLevel,
FirstDep,
SecondDep
);
delete buffer;
delete buffer1;
return num;
}
return 0;
}
int AutoComp(char *Filename,char *Filename1,char *CurveName1,char *CurveName2,
float cmins,
float winLength,float stepLength, float searchLength,
float *FirstDep,
float *SecondDep,
float *corc)
{
CMemRdWt well;
if(!well.Open(Filename,CMemRdWt::modeRead)) {
return 0;
};
CMemRdWt well1;
if(!well1.Open(Filename1,CMemRdWt::modeRead)) {
return 0;
}
int index=well.OpenCurve(CurveName1);
int index1=well1.OpenCurve(CurveName2);
if(index>-1&&index1>-1) {
Slf_CURVE curve,curve1;
well.GetCurveInfo(index,&curve);
int count=(curve.EndDepth-curve.StartDepth)/curve.DepLevel+1.5;
float *buffer=new float[count+1];
well.ReadCurveToFloatBuf(index,curve.StartDepth,count,buffer);
well1.GetCurveInfo(index1,&curve1);
int count1=(curve1.EndDepth-curve1.StartDepth)/curve1.DepLevel+1.5;
float *buffer1=new float[count1+1];
well1.ReadCurveToFloatBuf(index1,curve1.StartDepth,count1,buffer1);
int num=GetDepPairs(buffer,curve.StartDepth,curve.EndDepth,curve.DepLevel,
buffer1,curve1.StartDepth,curve1.EndDepth,curve1.DepLevel,
cmins,
winLength,
stepLength,
searchLength,
FirstDep,
SecondDep,
corc
);
delete buffer;
delete buffer1;
return num;
}
return 0;
}
bool AutoCompute(char *Filename,char *Filename1,char *CurveName1,char *CurveName2)
{
CMemRdWt well;
if(!well.Open(Filename,CMemRdWt::modeRead)) {
return false;
};
CMemRdWt well1;
if(!well1.Open(Filename1)) {
return false;
}
int iIndex=well.OpenTable("OG_RESULT");
if(iIndex<0) {
return false;
}
int index=well.OpenCurve(CurveName1);
int index1=well1.OpenCurve(CurveName2);
if(index1>-1&&index1>-1) {
Slf_CURVE curve,curve1;
well.GetCurveInfo(index,&curve);
well1.GetCurveInfo(index1,&curve1);
int count=(curve.EndDepth-curve.StartDepth)/curve.DepLevel+1.5;
float *buffer=new float[count+1];
int count1=(curve1.EndDepth-curve1.StartDepth)/curve1.DepLevel+1.5;
float *buffer1=new float[count1+1];
well.ReadCurveToFloatBuf(index,curve.StartDepth,count,buffer);
well1.ReadCurveToFloatBuf(index1,curve1.StartDepth,count1,buffer1);
float cmins=0.9;
float winLength=4;
float stepLength=0.5;
float searchLength=8;
float* FirstDep=new float[5000];
float* SecondDep=new float[5000];
float* corc=new float[5000];
int num=GetDepPairs(buffer,curve.StartDepth,curve.EndDepth,curve.DepLevel,
buffer1,curve1.StartDepth,curve1.EndDepth,curve1.DepLevel,
cmins,
winLength,
stepLength,
searchLength,
FirstDep,
SecondDep,
corc
);
delete buffer;
delete buffer1;
char buf[80];
LAYER_DATA m_Result;
int iIndex1=well1.OpenTable("RESULT");
int len=well.GetTableRecordLength(iIndex);
if (iIndex1 < 0)
{
iIndex1=well1.Open_Set_Table("RESULT",0,35,
"NO,SDEP,EDEP,ZONE,RESULTNO,THICK,TT,MDEPTH1,MDEPTH2,MDEPTH3,MDEPTH4,MDEPTH5,MDEPTH6,MDEPTH7,MDEPTH8,MDEPTH9,MDEPTH10,OIL,POOROIL,OILWATER,WATEROIL,GAS,POORGAS,GASWATER,WATERGAS,DEST1,DEST2,DEST3,DEST4,DEST5,DEST6,DEST7,DEST8,DEST9,DEST10,SDEP1,EDEP1,SDEP2,EDEP2,SDEP3,EDEP3,SDEP4,EDEP4,SDEP5,EDEP5,SDEP6,EDEP6,SDEP7,EDEP7,SDEP8,EDEP8,SDEP9,EDEP9,SDEP10,EDEP10",
"4,4,4,8,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,48,48,48,48,48,48,48,48,48,48,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4",//字段长度
"1,4,4,6,1,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,6,6,6,6,6,6,6,6,6,6,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4",//字段类型
"0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0");//字段备注
}
int iIndex2=well.OpenTable("RESULT");
if (iIndex2 < 0)
{
iIndex2=well.Open_Set_Table("RESULT",0,35,
"NO,SDEP,EDEP,ZONE,RESULTNO,THICK,TT,MDEPTH1,MDEPTH2,MDEPTH3,MDEPTH4,MDEPTH5,MDEPTH6,MDEPTH7,MDEPTH8,MDEPTH9,MDEPTH10,OIL,POOROIL,OILWATER,WATEROIL,GAS,POORGAS,GASWATER,WATERGAS,DEST1,DEST2,DEST3,DEST4,DEST5,DEST6,DEST7,DEST8,DEST9,DEST10,SDEP1,EDEP1,SDEP2,EDEP2,SDEP3,EDEP3,SDEP4,EDEP4,SDEP5,EDEP5,SDEP6,EDEP6,SDEP7,EDEP7,SDEP8,EDEP8,SDEP9,EDEP9,SDEP10,EDEP10",
"4,4,4,8,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,48,48,48,48,48,48,48,48,48,48,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4",//字段长度
"1,4,4,6,1,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,6,6,6,6,6,6,6,6,6,6,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4",//字段类型
"0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0");//字段备注
}
count=well.GetTableRecordCount(iIndex);
well.SetTableRecordCount(iIndex2,count);
well1.SetTableRecordCount(iIndex1,count);
int j=0;
for(int k=0;k<count;k++) {
well.ReadTable(iIndex,k+1,&buf);
well.WriteTable(iIndex2,k+1,&buf);
memmove(&m_Result,buf,sizeof(LAYER_DATA));
if(m_Result.StartDepth<curve.StartDepth||m_Result.StartDepth>curve.EndDepth) {
continue;
}
if(m_Result.EndDepth<curve.StartDepth||m_Result.EndDepth>curve.EndDepth) {
continue;
}
well1.WriteTable(iIndex1,j+1,&buf);
j++;
}
well1.SetTableRecordCount(iIndex1,j);
well.CloseTable(iIndex);
well1.CloseTable(iIndex1);
well.CloseTable(iIndex2);
char strName[256];
memset(strName, 0, sizeof (strName));
strcat(strName,"OG_RESULT1");
well1.EShiftTableDepth(strName, num,FirstDep,SecondDep);
delete FirstDep;
delete SecondDep;
delete corc;
return true;
}
return false;
}
////////////////////////////////////////////////////////////////////////////
// DepPairs对象
// 显示边框
// 构造函数
DepPairs::DepPairs(float* curve1,float stdep1,float endep1,float rlev1,
float* curve2,float stdep2,float endep2,float rlev2,float cmins,
float winLength,float stepLength, float searchLength)
{
Curve1.stdep=stdep1;
Curve1.endep=endep1;
Curve1.rlev=rlev1;
Curve1.val=curve1;
Curve2.stdep=stdep2;
Curve2.endep=endep2;
Curve2.rlev=rlev2;
Curve2.val=curve2;
// num=i;
DepLevel=rlev1;
if(DepLevel>rlev2) DepLevel=rlev2;
StartDepth=stdep1;
if(StartDepth>stdep2) StartDepth=stdep2;
EndDepth=endep1;
if(EndDepth<endep2) EndDepth=endep2;
// int num=(EndDepth-StartDepth)/DepLevel+1.5;
m_nDepPairsNum=0;
wl=winLength;
sl=searchLength;
step=stepLength;
dr1=3.;
dr2=3.;
cmin=cmins;
elast=1.;
emp=0;
size=1;
type=0;
way=0;
ipar=0;
SetRef=0;
mCurrent=-1;
ShiftNum=0;
}
// 析构函数
DepPairs::~DepPairs()
{
}
// read depth pair from file
void DepPairs::qs(float *datx,float *daty,float *datc,int left,int right)
{
int i,j;
float x,y;
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;
y=daty[i];daty[i]=daty[j];daty[j]=y;
y=datc[i];datc[i]=datc[j];datc[j]=y;
i++;j--;
}
}while(i<=j);
if(left<j) qs(datx,daty,datc,left,j);
if(i<right) qs(datx,daty,datc,i,right);
}
void DepPairs::detscl(float sdeps,float edeps)
{
int k;
float tempdep;
tempdep = sdeps;
float xmin1=99999.;
float xmin2=99999.;
float xmax1=-99999.;
float xmax2=-99999.;
k=(int)((sdeps-sdeps)/ DepLevel+0.5);
float val=0;
while(tempdep<edeps)
{
val=Curve1.GetValue(tempdep);
if(val<xmin1) xmin1=val;
if(val>xmax1) xmax1=val;
val=Curve2.GetValue(tempdep);
if(val<xmin2) xmin2=val;
if(val>xmax2) xmax2=val;
k++;
tempdep=sdeps + k * DepLevel;
}
if(xmin1==xmax1)
{
xmin1=0.;
xmax1=1.;
}
if(xmin2==xmax2)
{
xmin2=0.;
xmax2=1.;
}
return;
}
float DepPairs::ampl(int n,float *x)
{
int i;
float xmin,xmax;
xmin=9999.;
xmax=-9999.;
for(i=0;i<n;i++) {
if(x[i]<xmin)xmin=x[i];
if(x[i]>xmax)xmax=x[i];
}
return xmax-xmin;
}
float DepPairs::corr(int n,float *x,float *y,int size,int emp)
{
float sx,sy,xave,yave,sqx,sqy,sxy,dp,dx,dy,c;
float r1,r2,amp1,amp2,ratio;
int i;
sx=0.;
sy=0.;
for(i=0;i<n;i++) {
sx=sx+x[i];
sy=sy+y[i];
}
xave=sx/n;
yave=sy/n;
sqx=0.;
sqy=0.;
for(i=0;i<n;i++) {
sqx=sqx+(x[i]-xave)*(x[i]-xave);
sqy=sqy+(y[i]-yave)*(y[i]-yave);
}
dx=sqrt(sqx/n);
dy=sqrt(sqy/n);
sxy=0.;
for(i=0;i<n;i++) {
sxy=sxy+(x[i]-xave)*(y[i]-yave);
}
dp=sxy/n;
c=0.;
if(dx*dy!=0.)c=dp/(dx*dy);
if(size==0.) {
r1=1.;
r2=1.;
if(dx!=0.)r1=dy/dx;
if(dy!=0.)r2=dx/dy;
if(r1>=r2)ratio=r1;
if(r1<r2)ratio=r2;
c=c/ratio;
}
if(emp!=0) {
amp1=ampl(n,x);
amp2=ampl(n,y);
ratio=amp1;
if(ratio>amp2)ratio=amp2;
c=c*ratio;
}
return c;
}
void DepPairs::cor(int iwl,int isl,float *x,float *y,float *c,int size,int emp)
{
int i,ix=isl,iy;
for(i=0;i<2*isl;i++) {
ix=isl;
iy=i;
c[i]=corr(iwl,&x[ix],&y[iy],size,emp);
}
}
double DepPairs::fun(DPoint &p1, DPoint &p2, DPoint &p3)
{
double s;
s= (p2.x-p1.x)*(p3.y-p1.y)+(p1.y-p2.y)*(p3.x-p1.x);
return s;
}
float DepPairs::act(float *x,int n,int way)
{
int i,i1,i2;
float sum,sq,ave,z,sqd;
sum=0.;
for(i=0;i<n;i++) sum=sum+x[i];
ave=sum/n;
sq=0.;
for(i=0;i<n;i++) sq=sq+(x[i]-ave)*(x[i]-ave);
sqd=sqrt(sq/n);
if(way==0) return sqd;
i1=0;
i2=0;
for(i=0;i<n;i++) {
if(i>=n/2&&x[i]>=ave)i1++;
if(i>=n/2&&x[i]<ave)i2++;
if(i<n/2&&x[i]>=ave)i2++;
if(i<n/2&&x[i]<ave)i1++;
}
z=sqd;
if(way==1&&i2>i1)z=-z;
if(way==2&&i2<=i1)z=-z;
return z;
}
void DepPairs::actvty(int iwl,int isl,float *dr,int way,float *x,float *c)
{
int n,n2,i;
n= (int)(*dr / DepLevel);
n2=n/2;
n=2*n2+1;
for(i=0;i<iwl+2*isl;i++) c[i]=act(&x[i],n,way);
for(i=0;i<iwl+2*isl;i++) x[i]=c[i];
}
void DepPairs::adjust(int n,int way,float *c)
{
int i;
for(i=0;i<n;i++) {
if(way==0)c[i]=abs(c[i]);
if(way==2)c[i]=-c[i];
}
}
int DepPairs::MMax(float depth,int isl,float *x,float *y,float *c,float cmins)
{
int i,iup,idown,n,idir1=0,idir2=0;
float xdep,ydep;
n=2*isl+1;
iup=1;
idown=-1;
if( c[0] < c[1] ) idir1 = iup;
if( c[0] >= c[1] ) idir1 = idown;
for(i=1;i<n-1;i++)
{
if( c[i] < c[i+1] )idir2=iup;
if( c[i] >= c[i+1] )idir2=idown;
if( idir1 == idown && idir2 == iup ) idir1=idir2;
if( idir1 == iup && idir2 == idown )
{
xdep=depth;
ydep=depth+(i-isl)*DepLevel;
if ( c[i] >= 0.5) {
fprintf(fp_temp,"%8.3f %8.3f %5.3f\n",ydep,xdep,c[i]);
// break;
}
idir1=idir2;
}
}
return 1;
}
int DepPairs::MMax(float ydep,float dep0,int isl,float *c,float cmins)
{
float xdep;
int ss=-1;
int n=1.5*isl+1;
float maxcor=-99999;
for(int i=0;i<n;i++)
{
if( c[i] >maxcor) {
maxcor=c[i];
ss=i;
}
}
if(ss>-1&&c[ss]>cmins)
{
xdep=dep0+ss*DepLevel;
fprintf(fp_temp,"%8.3f %8.3f %5.3f\n",xdep,ydep,c[ss]);
}
return 1;
}
int DepPairs::signum(float x)
{
int y;
y=0;
if(x<0.)y=-1;
if(x>0.)y=1;
return y;
}
int DepPairs::optim(int n,float *x1,float *y1,float *c1,float cmins,float elast)
{
float xmax=999999.;
int xnull=0;
float tfine=0,check,diff=0,fcrif=0,fine=0,fmin=0;
int i,iw,j,ip1,ip2,imin,nd,ids,nd1;
n=1;
float a,b,d;
int aa=0;
rewind(fp_temp);
while(!feof(fp_temp))
{
fscanf(fp_temp,"%f%f%f\n",&a,&b,&d);
if(ipar) aa++;
else if(d>=cmins) aa++;
}
if(aa<=0) return -1;
float *x=new float[aa+1000];
float *y=new float[aa+1000];
float *c=new float[aa+1000];
float *f=new float[aa+1000];
int *ncs=new int[aa+1000];
rewind(fp_temp);
while(!feof(fp_temp))
{
if(fscanf(fp_temp,"%f%f%f\n",&x[n],&y[n],&c[n])!=3)
{
break;
};
if(ipar) n++;
else if(c[n]>=cmins) n++;
}
fclose(fp_temp);
QString dir=GetLogdataPath();
dir+="\\temp";
CreateDir((char *)dir.toStdString().c_str());
dir+="\\admtemp";
fp_temp=fopen(dir.toStdString().c_str(),"w+");
if(ipar)
{
for(ip1=0;ip1<n;ip1++)
{
ncs[ip1]=0;
}
for(int i=1;i<n;i++)
{
for(j=i-1;j>0;j--)
{
if(x[i]<=x[j]){
ncs[i]=1;
break;
}
}
if(x[i]-x[i-1]>0&&x[i]-x[i-1]<0.5)ncs[i]=1;
if(!ncs[i])
{
for(j=i+1;j<n-1;j++)
{
if(x[i]>=x[j]){
ncs[i]=1;
break;
}
}
}
}
rewind(fp_temp);
j=0;
for(int i=1;i<n;i++){
if(ncs[i]) continue;
fprintf(fp_temp,"%f %f %f\n",x[i],y[i],c[i]);
j++;
}
}
else {
x[0]=xmax; y[0]=xmax; c[0]=0.;
x[n]=xmax; y[n]=xmax; c[n++]=0.;
f[0]=0; ncs[0]=xnull;
fine=elast;
for(ip1=1;ip1<n;ip1++)
{
f[ip1]=0;
ncs[ip1]=0;
}
for(ip1=1;ip1<n;ip1++)
{
if(x[ip1]>=xmax&&ip1!=n-1) goto l2100;
fmin=xmax;
tfine=0.;
for(iw=0;iw<ip1;iw++)
{
ip2=ip1-iw-1;
if(iw!=0&&x[ip2+1]<xmax) tfine=tfine+fine*c[ip2+1];
if(x[ip2]>=xmax&&ip2!=0) goto l1900;
if(tfine>fmin) goto l2000;
if(x[ip1]<xmax&&x[ip2]<xmax) goto l1500;
fcrif=f[ip2]+tfine;
goto l1800;
l1500:;
check=signum(x[ip1]-x[ip2])*signum(y[ip1]-y[ip2]);
if(check<=0.0) goto l1900;
diff=abs(( x[ip1]-x[ip2]) - (y[ip1]-y[ip2]));
fcrif=f[ip2]+tfine+diff;
l1800:;
if(fcrif>fmin) goto l1900;
fmin=fcrif;
imin=ip2;
l1900:;
}
l2000:;
f[ip1]=fmin;
ncs[ip1]=imin;
l2100:;
}
nd=0;
ids=n-1;
while(1)
{
ids=ncs[ids];
if(ids==0) break;
f[nd++]=ids;
}
for(i=0;i<nd;i++)
{
j=nd-i-1;
ncs[i]=f[j];
}
for(nd1=0,i=0;i<nd;i++)
{
if(x[ncs[i]]<xmax)
{
x[nd1]=x[ncs[i]];
y[nd1]=y[ncs[i]];
c[nd1++]=c[ncs[i]];
}
}
rewind(fp_temp);
fprintf(fp_temp,"%f %f %f\n",x[0],y[0],c[0]);
j=0;
float d1=0;
j++;
for(int i=1;i<nd1;i++) {
d=x[i]-x[i-1];
d1=y[j]-y[i-1];
if(d!=d1||i==nd1-1)
{
fprintf(fp_temp,"%f %f %f\n",x[i],y[i],c[i]);
j++;
}
}
}
delete x;
delete y;
delete c;
delete f;
delete ncs;
return j;
}
int DepPairs::GetDots(int n1,float sdeps,float edeps,float DepLevel,float *datx1,int n)
{
int count=(edeps-sdeps)/DepLevel+1.5;
DPoint*p1=new DPoint[count+1];
for(int i=0;i<count;i++)
{
if(datx1)
{
p1[i]=DPoint(sdeps+i*DepLevel,datx1[i]);
}
else p1[i]=DPoint(sdeps+i*DepLevel,Curve1.GetValue(sdeps+i*DepLevel));
}
int datn=GetGDots(count,p1,n);
delete p1;
return datn;
}
int DepPairs::GetGDots(int n1,DPoint *x,int n)
{
int k=0;
double s1, s2;
int *spos=new int[n1];
for(int i=1; i<n1-1; i++)
{
s1 = x[i].y-x[i-1].y;//fun(x[i-2], x[i-1], x[i]);
s2 = x[i].y-x[i+1].y;//fun(x[i-1], x[i], x[i+1]);
if(s1<0&&s2<0){
spos[k]=i;
k++;
}
}
int kk=0;
for(int j=0;j<n;j++)
{
if(kk)k=kk;
kk=0;
for(int i=1; i<k-1; i++)
{
s1 = x[spos[i]].y-x[spos[i-1]].y;//fun(x[i-2], x[i-1], x[i]);
s2 = x[spos[i]].y-x[spos[i+1]].y;//fun(x[i-1], x[i], x[i+1]);
if(s1>0&&s2<0||s1>0&&s2>0||s1<0&&s2>0){
continue;
}
spos[kk]=spos[i];
kk++;
}
}
if(n<=0) kk=k;
for(int i=0;i<kk-1;i++)
{
int ppos=-1;
float d=99999;
float dd;
int len=spos[i+1]-1-spos[i];
int isl=sl/DepLevel;
if(len>isl) len=isl;
for(int m=spos[i];m<spos[i]+len;m++)
{
dd=(x[m].y-x[m+1].y);
if(dd<d)
{
ppos=m;
d=dd;
}
}
if(ppos>0){
spos[i]=ppos-1;
}
}
for(int i=0; i<kk; i++)
{
fprintf(fp_temp,"%f %f %f\n",x[spos[i]].x,x[spos[i]].x,1);
}
delete spos;
return kk;
}
int DepPairs::paral(int n1,float *x1,float *y1,float *c1,int retype)
{
float a,b,d;
float *x=new float[n1+1];
float *y=new float[n1+1];
float *c=new float[n1+1];
float *spos=new float[n1+1];
memset(spos,0,n1*sizeof(float));
rewind(fp_temp);
int n=0;
while(!feof(fp_temp))
{
fscanf(fp_temp,"%f%f%f",&x[n],&y[n],&c[n]);
n++;
if(n>=n1) break;
}
float xmax=999999.,disp1,disp2,disp;
int i,k,j;
if(n<3) {
delete x;
delete y;
delete c;
delete spos;
return -1;
}
i=-1;
float v1,v2,v3;
fclose(fp_temp);
QString dir=GetLogdataPath();
dir+="\\temp";
CreateDir((char *)dir.toStdString().c_str());
dir+="\\admtemp";
fp_temp=fopen(dir.toStdString().c_str(),"w+");
while(1) {
i=i+1;
if(i>=n) break;
disp1=y[i]-x[i];
disp2=y[i+1]-x[i+1];
if(disp1==disp2) {
i+=2;
for(j=i;j<n;j++) {
disp=y[j]-x[j];
if(disp==disp1)x[j-1]=xmax;
if(disp!=disp1) break;
}
if(disp==disp1) break;
i=j-1;
}
}
for(k=0,i=0;i<n;i++) {
if(x[i]<xmax) {
fprintf(fp_temp,"%f %f %f\n",x[i],y[i],c[i]);
k++;
}
}
delete x;
delete y;
delete c;
delete spos;
return k;
}
void DepPairs::DepthSort(float *xxx,int number)
{
int i5,j5;
float tempf5;
for (i5=0; i5<number; i5++)
for (j5=i5; j5<number; j5++)
if (xxx[i5] > xxx[j5])
{
tempf5 = xxx[i5];
xxx[i5] = xxx[j5];
xxx[j5] = tempf5;
}
}
int DepPairs::GetCompare(int datn,float sdeps,float edeps,float* datx1,float* daty1,float* datc1,float cmins,float xelast)
{
float *x=new float[datn+10];
float *y=new float[datn+10];
float *c=new float[datn+10];
rewind(fp_temp);
float wl2,xdr2;
int num1,number,numt,i,iwl,isl,iwl2,istep,jsize,ixdr2,idr,idr2;
int n1,kno;
int cccc;
if(wl==0)wl=10.;
if(sl==0)sl=4.;
if(step==0)step=1;
if(dr1==0)dr1=3.;
if(dr2==0)dr2=3.;
if(type!=0 &&sl<dr2/2)sl=dr2/2;
iwl=(int)(wl/DepLevel+0.5);
iwl2=iwl/2;
iwl=2*iwl2+1;
isl=(int)(sl/DepLevel+0.5);
istep=(int)(step/DepLevel+0.5);
idr=(int)(dr2/DepLevel+0.5);
idr2=idr/2;
idr=2*idr2+1;
xdr2=idr2*DepLevel;
if( type == 0 ) xdr2=0;
ixdr2=idr2;
if( type == 0 ) ixdr2=0;
jsize=iwl+2*isl;
kno=1;
if(emp!=0||size==0) detscl(sdeps,edeps);
num1 = 0;
number = iwl + 2 * isl+2*ixdr2;
numt = number - istep;
wl2=iwl2*DepLevel;
int n=0;
int maxnum=-99999;
iwl=(int)(wl/DepLevel+0.5);
iwl2=iwl/2;
iwl=2*iwl2+1;
while(!feof(fp_temp))
{
fscanf(fp_temp,"%f%f%f",&x[n],&y[n],&c[n]);
if(n)iwl=(x[n]-x[n-1])/DepLevel;
if(maxnum<iwl) maxnum=iwl;
n++;
if(n>=datn) break;
}
fclose(fp_temp);
QString dir=GetLogdataPath();
dir+="\\temp";
CreateDir((char *)dir.toStdString().c_str());
dir+="\\admtemp";
fp_temp=fopen(dir.toStdString().c_str(),"w+");
number = maxnum + 2 * isl+2*ixdr2;
float *datx=new float[number+1];
float *daty=new float[number+1];
float *datc=new float[number+1];
cccc = 0;
for(int j=0;j<n-1;j++){
int k=isl;
for (i=0; i<number; i++) {
datx[i]=0.;
daty[i]=0.;
datc[i]=0.;
}
for(float dep=x[j];dep<x[j+1];dep+=DepLevel)
{
if(datx1)
{
int pos=(dep-sdeps)/DepLevel;
if(dep<sdeps||dep>edeps) datx[k++]=-9999.0;
else datx[k++]=datx1[pos];
}
else datx[k++]=Curve1.GetValue(dep);
}
k=0;
for(float dep=x[j]-isl*DepLevel/2;dep<x[j+1]+isl*DepLevel;dep+=DepLevel)
{
if(daty1)
{
int pos=(dep-sdeps)/DepLevel;
if(dep<sdeps||dep>edeps) daty[k++]=-9999.0;
else daty[k++]=daty1[pos];
}
else daty[k++]=Curve2.GetValue(dep);
}
iwl=(x[j+1]-x[j])/DepLevel;
cor(iwl,isl,datx,daty,datc,size,emp);
MMax(x[j],x[j]-isl*DepLevel/2,isl,datc,cmins);
}
datn=optim(0,datx1,daty1,datc1,cmins,xelast);
delete x;
delete y;
delete c;
return datn;
}
int DepPairs::joincurve(float sdeps,float edeps,float *datx1,float *daty1,float *datc1,
float wl,float sl,float step,float dr1,float dr2,
int type,int size,int way,int emp,int ipar,
float cmins,float xelast)
{
float depth,depth1,dept,deps;
float wl2,xdr2;
int num1,number,numt,i,iwl,isl,iwl2,istep,jsize,ixdr2,idr,idr2;
int n1,kno;
int cccc;
if(wl==0)wl=10.;
if(sl==0)sl=4.;
if(step==0)step=1;
if(dr1==0)dr1=3.;
if(dr2==0)dr2=3.;
if(type!=0 &&sl<dr2/2)sl=dr2/2;
iwl=(int)(wl/DepLevel+0.5);
iwl2=iwl/2;
iwl=2*iwl2+1;
isl=(int)(sl/DepLevel+0.5);
istep=(int)(step/DepLevel+0.5);
idr=(int)(dr2/DepLevel+0.5);
idr2=idr/2;
idr=2*idr2+1;
xdr2=idr2*DepLevel;
if( type == 0 ) xdr2=0;
ixdr2=idr2;
if( type == 0 ) ixdr2=0;
jsize=iwl+2*isl;
kno=1;
if(emp!=0||size==0) detscl(sdeps,edeps);
num1 = 0;
number = iwl + 2 * isl+2*ixdr2;
numt = number - istep;
wl2=iwl2*DepLevel;
depth=sdeps-step+wl2+sl+xdr2;
depth1 = depth + step + step/2;
dept = iwl2 + isl;
int ret=TRUE;
QString dir=GetLogdataPath();
dir+="\\temp";
CreateDir((char *)dir.toStdString().c_str());
dir+="\\admtemp";
fp_temp=fopen(dir.toStdString().c_str(),"w+");
if(fp_temp==NULL)
{
QMessageBox::warning(nullptr, "提示", "Not create file :admtemp !\n");
return -1;
}
// ::GetProgressBar()->setVisible(true);
// ::GetProgressBar()->setRange(sdeps,edeps);
// ::GetProgressBarHint()->setVisible(true);
char showname[200];
sprintf(showname,"自动产生深度对比线...");
int datn=0;
if(ipar){
// ::GetProgressBarHint()->setText("查找层界面....");
QApplication::processEvents();
datn=GetDots(n1,sdeps,edeps,DepLevel,datx1,1);
// ::GetProgressBar()->setValue(sdeps+(edeps-sdeps)*0.1);
QApplication::processEvents();
// ::GetProgressBarHint()->setText("开始地层对比...");
// ::GetProgressBar()->setValue(sdeps+(edeps-sdeps)*0.3);
QApplication::processEvents();
datn=GetCompare(datn,sdeps,edeps,datx1,daty1,datc1,cmins,xelast);
// ::GetProgressBar()->setValue(edeps);
QApplication::processEvents();
}
else{
float *datx=new float[number+1];
float *daty=new float[number+1];
float *datc=new float[number+1];
cccc = 0;
for (i=0; i<number; i++) {
datx[i]=0.;
daty[i]=0.;
datc[i]=0.;
}
while(1)
{
depth=depth+step;
// ::GetProgressBar()->setValue(depth);
// ::GetProgressBarHint()->setText(showname);
QApplication::processEvents();
if(depth>edeps-wl2-sl-xdr2) {
break;
}
if ( depth > depth1 )
{
if ( numt > 0 )
{
for(i=0;i<numt;i++)
{
datx[i] = datx[istep+i];
daty[i] = daty[istep+i];
}
num1 = numt;
}
}
for(i=num1;i<number;i++)
{
deps=depth+(i-dept)*DepLevel;
if(emp!=0||size==0)
{
datx[i]=(Curve1.GetValue(deps)-xmin1)/(xmax1-xmin1);
daty[i]=(Curve2.GetValue(deps)-xmin2)/(xmax2-xmin2);
}
else
{
datx[i]=Curve1.GetValue(deps);
daty[i]=Curve2.GetValue(deps);
}
}
if(type!=0)
{
actvty(iwl,isl,&dr1,way,datx,datc);
actvty(iwl,isl,&dr2,way,daty,datc);
}
cor(iwl,isl,datx,daty,datc,size,emp);
n1=iwl+isl*2;
if(type==0) adjust(n1,way,datc);
cccc ++;
MMax(depth,isl,datx,daty,datc,cmins);
}
fseek(fp_temp,0,0);
delete datx;
delete daty;
delete datc;
datn=optim(0,datx1,daty1,datc1,cmins,xelast);
if(datn>0)datn=paral(datn,datx1,daty1,datc1,ipar);
}
fclose(fp_temp);
// ::GetProgressBar()->setVisible(false);
// ::GetProgressBarHint()->setVisible(false);
QApplication::processEvents();
return datn;
}
int DepPairs::adm(float stem,float etem ,float cmins,
float winLength,float stepLength, float searchLength,
float *datx1,float *daty1,float *datc1
)
{
char str[120];
int i;
int datn;
cmin = cmins;
datn = joincurve(stem,etem,datx1,daty1,datc1,
winLength,searchLength,stepLength,dr1,dr2,type,size,
way,emp,ipar,cmins,elast);
if ( datn < 0 )
{
sprintf(str,"%f--%f层段相关性太差相关系数均小于截止值%f ",stem,etem,cmins);
QMessageBox::warning(nullptr, "提示", str);
}
else {
float *datx,*daty,*datc;
m_nDepPairsNum=0;
datx = new float[datn+100];
if ( datx == NULL )
{
QMessageBox::warning(nullptr, "提示", "自动校深:没有足够的内存.");
return m_nDepPairsNum;
}
daty = new float[datn+100];
if ( daty == NULL )
{
QMessageBox::warning(nullptr, "提示", "自动校深:没有足够的内存.");
delete datx;
return m_nDepPairsNum;
}
datc = new float[datn+100];
if ( datc == NULL )
{
QMessageBox::warning(nullptr, "提示", "自动校深:没有足够的内存.");
delete datx;
delete daty;
return m_nDepPairsNum;
}
QString dir=GetLogdataPath();
dir+="\\temp";
CreateDir((char *)dir.toStdString().c_str());
dir+="\\admtemp";
fp_temp=fopen(dir.toStdString().c_str(),"rt");
if(fp_temp==NULL)
{
QMessageBox::warning(nullptr, "提示", "Not create file :admtemp !\n");
delete [] datx;
delete [] daty;
delete [] datc;
QDir r;
r.remove(dir.toStdString().c_str());
return -1;
}
i=0;
while(!feof(fp_temp))
{
fscanf(fp_temp,"%f%f%f",&datx[i],&daty[i],&datc[i]);
i++;
if(i>=datn) break;
}
if(!i) {
QMessageBox::warning(nullptr, "提示", "自动对比文件读写错误!");
fclose(fp_temp);
delete [] datx;
delete [] daty;
delete [] datc;
QDir r;
r.remove(dir.toStdString().c_str());
return -1;
}
datn=i;
fclose(fp_temp);
qs(datx,daty,datc,0,datn-1);
if(last_edep!=stem)
{
float ddep=0;
if(datn) {
ddep=datx[0]-stem;
FirstDep[m_nDepPairsNum]=daty[0]-ddep;
}
else FirstDep[m_nDepPairsNum]=stem;
SecondDep[m_nDepPairsNum]=stem;
corc[m_nDepPairsNum]=1.;
m_nDepPairsNum++;
}
for(i=0;i<datn;i++)
{
FirstDep[m_nDepPairsNum]=daty[i];
SecondDep[m_nDepPairsNum]=datx[i];
if(corc) corc[m_nDepPairsNum]=datc[i];
m_nDepPairsNum++;
}
float ddep=0;
if(datn) {
ddep=etem-datx[datn-1];
FirstDep[m_nDepPairsNum]=daty[datn-1]+ddep;
}
else FirstDep[m_nDepPairsNum]=etem;
SecondDep[m_nDepPairsNum]=etem;
if(corc) corc[m_nDepPairsNum]=1;
m_nDepPairsNum++;
delete [] datx;
delete [] daty;
delete [] datc;
QDir r;
r.remove(dir.toStdString().c_str());
}
return m_nDepPairsNum;
}
int DepPairs::AdmRun(float *curve1,float *curve2,float *cc)
{
return adm(StartDepth,EndDepth ,cmin,wl,step,sl,curve1,curve2,cc);
}