月份：2023年2月

MFC实现三维图像绘制

绘制物体的面时考虑两个问题：

1. 如何绘制？

2. 如何消隐不可见面？

针对这两个问题，分别采取有效边表算法和Z-Buffer缓冲器算法。

有效边表：

#include “Pi3.h”
#include “T2.h”

class CAET
{
public:
CAET();
virtual ~CAET();
public:
double x;//当前扫描线与有效边交点的x坐标
int yMax;//边的最大y值
double k;//斜率的倒数(x的增量)
CPi3 ps;//边的起点
CPi3 pe;//边的终点
CAET *pNext;
CT2 ts;//起点纹理坐标
CT2 te;//终点纹理坐标
};
桶表：

#include “AET.h”
class CBucket
{
public:
CBucket();
virtual ~CBucket();
public:
int ScanLine; //扫描线
CAET *pET; //桶上的边表指针
CBucket *pNext;
};
具体绘制实现（这里涉及到物体的光照和纹理模型）：

#include “Bucket.h”
#include “Vector3.h”
#include “Lighting.h”
#include “Pi3.h”
#include “Texture.h”
#include “T2.h”

class CFill
{
public:
CFill(void);
virtual ~CFill(void);
void SetPoint(CPi3 *p, int, CT2*);//类的初始化
void SetPoint(CPi3* p, int);//类的初始化
void CreateBucket();//创建桶
void CreateEdge();//边表
void AddEt(CAET *);//合并ET表
void EtOrder();//ET表排序
void Gouraud(CDC *, int,CLighting, CMaterial*,CP3, CTexture*,BOOL);//填充多边形
void Gouraud(CDC*, int, CLighting, CMaterial*, CP3);//填充多边形
CRGB Interpolation(double, double, double, CRGB, CRGB);//颜色线性插值
double Interpolation(double, double, double, double, double);//深度线性插值
CVector3 Interpolation(double, double, double, CVector3, CVector3);//法向量线性插值
CT2 Interpolation(double m, double m0, double m1, CT2 T0, CT2 T1);//纹理地址线性插值
void ClearMemory();//清理内存
void DeleteAETChain(CAET* pAET);//删除边表
CRGB GetTextureColor(double u, double v, CTexture* pTexture);//获得纹理颜色
void InitialDepthBuffer(int nWidth, int nHeight, double zDepth);//初始化深度缓冲器
public:
int PNum;//顶点个数
CT2* T;//纹理动态数组
CPi3 *P;//顶点坐标动态数组
CAET *pHeadE, *pCurrentE, *pEdge; //有效边表结点指针
CBucket *pHeadB, *pCurrentB; //桶表结点指针
double** zBuffer;//深度缓冲器
int nWidth, nHeight;//缓冲区宽度与高度
};

#include “pch.h”
#include “Fill.h”
#include “math.h”
#define ROUND(f) int(floor(f+0.5))//四舍五入宏定义

CFill::CFill(void)
{
PNum = 0;
P = NULL;
pEdge = NULL;
pHeadB = NULL;
pHeadE = NULL;
}

CFill::~CFill(void)
{
if (P != NULL)
{
delete[] P;
P = NULL;
}
for (int i = 0;i < nWidth;i++)
{
delete[] zBuffer[i];
zBuffer[i] = NULL;
}
if (zBuffer != NULL)
{
delete zBuffer;
zBuffer = NULL;
}
ClearMemory();
}

void CFill::SetPoint(CPi3 *p, int m, CT2* t)
{
if (P != NULL)
{
delete[] P;
P = NULL;
}
P = new CPi3[m];//创建一维动态数组
T = new CT2[m];
for (int i = 0; i < m; i++)
{
P[i] = p[i];
P[i].y = ROUND(P[i].y);
T[i] = t[i];
}
PNum = m;
}

void CFill::SetPoint(CPi3* p, int m)
{
if (P != NULL)
{
delete[] P;
P = NULL;
}
P = new CPi3[m];//创建一维动态数组
for (int i = 0; i < m; i++)
{
P[i] = p[i];
P[i].y = ROUND(P[i].y);
}
PNum = m;
}

void CFill::CreateBucket()//创建桶表
{
int yMin, yMax;
yMin = yMax = P[0].y;
for (int i = 0; i < PNum; i++)//查找多边形所覆盖的最小和最大扫描线
{
if (P[i].y < yMin)
{
yMin = P[i].y;//扫描线的最小值
}
if (P[i].y > yMax)
{
yMax = P[i].y;//扫描线的最大值
}
}
for (int y = yMin; y <= yMax; y++)
{
if (yMin == y)//如果是扫描线的最小值
{
pHeadB = new CBucket;//建立桶的头结点
pCurrentB = pHeadB;//pCurrentB为CBucket当前结点指针
pCurrentB->ScanLine = yMin;
pCurrentB->pET = NULL;//没有链接边表
pCurrentB->pNext = NULL;
}
else//其他扫描线
{
pCurrentB->pNext = new CBucket;//建立桶的其他结点
pCurrentB = pCurrentB->pNext;
pCurrentB->ScanLine = y;
pCurrentB->pET = NULL;
pCurrentB->pNext = NULL;
}
}
}

void CFill::CreateEdge()//创建边表
{
for (int i = 0; i < PNum; i++)
{
pCurrentB = pHeadB;
int j = (i + 1) % PNum;//边的第2个顶点，P[i]和P[j]点对构成边
if (P[i].y < P[j].y)//边的终点比起点高
{
pEdge = new CAET;
pEdge->x = P[i].x;//计算ET表的值
pEdge->yMax = P[j].y;
pEdge->k = (P[j].x – P[i].x) / (P[j].y – P[i].y);//代表1/k
pEdge->ps = P[i];//绑定顶点和颜色
pEdge->pe = P[j];
pEdge->ts= T[i];//绑定纹理
pEdge->te = T[j];
pEdge->pNext = NULL;
while (pCurrentB->ScanLine != P[i].y)//在桶内寻找当前边的yMin
{
pCurrentB = pCurrentB->pNext;//移到yMin所在的桶结点
}
}
if (P[j].y < P[i].y)//边的终点比起点低
{
pEdge = new CAET;
pEdge->x = P[j].x;
pEdge->yMax = P[i].y;
pEdge->k = (P[i].x – P[j].x) / (P[i].y – P[j].y);
pEdge->ps = P[i];
pEdge->pe = P[j];
pEdge->ts = T[i];//绑定纹理
pEdge->te = T[j];
pEdge->pNext = NULL;
while (pCurrentB->ScanLine != P[j].y)
{
pCurrentB = pCurrentB->pNext;
}
}
if (P[i].y != P[j].y)
{
pCurrentE = pCurrentB->pET;
if (pCurrentE == NULL)
{
pCurrentE = pEdge;
pCurrentB->pET = pCurrentE;
}
else
{
while (pCurrentE->pNext != NULL)
{
pCurrentE = pCurrentE->pNext;
}
pCurrentE->pNext = pEdge;
}
}
}
}

void CFill::Gouraud(CDC *pDC, int light_type,CLighting light, CMaterial* material, CP3 EyePoint, CTexture* pTexture,BOOL IfBump)//填充多边形
{
CAET *pT1 = NULL, *pT2 = NULL;
pHeadE = NULL;
for (pCurrentB = pHeadB; pCurrentB != NULL; pCurrentB = pCurrentB->pNext)
{
for (pCurrentE = pCurrentB->pET; pCurrentE != NULL; pCurrentE = pCurrentE->pNext)
{
pEdge = new CAET;
pEdge->x = pCurrentE->x;
pEdge->yMax = pCurrentE->yMax;
pEdge->k = pCurrentE->k;
pEdge->ps = pCurrentE->ps;
pEdge->pe = pCurrentE->pe;
pEdge->pNext = NULL;
pEdge->ts = pCurrentE->ts;
pEdge->te = pCurrentE->te;
pEdge->pNext = NULL;
AddEt(pEdge);
}
EtOrder();
pT1 = pHeadE;
if (pT1 == NULL)
{
return;
}
while (pCurrentB->ScanLine >= pT1->yMax)//下闭上开
{
CAET * pAETTEmp = pT1;
pT1 = pT1->pNext;
delete pAETTEmp;
pHeadE = pT1;
if (pHeadE == NULL)
return;
}
if (pT1->pNext != NULL)
{
pT2 = pT1;
pT1 = pT2->pNext;
}
while (pT1 != NULL)
{
if (pCurrentB->ScanLine >= pT1->yMax)//下闭上开
{
CAET* pAETTemp = pT1;
pT2->pNext = pT1->pNext;
pT1 = pT2->pNext;
delete pAETTemp;
}
else
{
pT2 = pT1;
pT1 = pT2->pNext;
}
}
CRGB Ca, Cb, Cf;//Ca、Cb代表边上任意点的颜色，Cf代表面上任意点的颜色
double za, zb, zf;
CVector3 va, vb, vf;
CT2 ta, tb, tf;//ta和tb代表边上任意点的纹理地址，tf代表面上任意点的纹理地址
ta = Interpolation(pCurrentB->ScanLine, pHeadE->ps.y, pHeadE->pe.y, pHeadE->ts, pHeadE->te);
tb = Interpolation(pCurrentB->ScanLine, pHeadE->pNext->ps.y, pHeadE->pNext->pe.y, pHeadE->pNext->ts, pHeadE->pNext->te);
Ca = Interpolation(pCurrentB->ScanLine, pHeadE->ps.y, pHeadE->pe.y, pHeadE->ps.c, pHeadE->pe.c);
Cb = Interpolation(pCurrentB->ScanLine, pHeadE->pNext->ps.y, pHeadE->pNext->pe.y, pHeadE->pNext->ps.c, pHeadE->pNext->pe.c);
za = Interpolation(pCurrentB->ScanLine, pHeadE->ps.y, pHeadE->pe.y, pHeadE->ps.z, pHeadE->pe.z);
zb = Interpolation(pCurrentB->ScanLine, pHeadE->pNext->ps.y, pHeadE->pNext->pe.y, pHeadE->pNext->ps.z, pHeadE->pNext->pe.z);
va = Interpolation(pCurrentB->ScanLine, pHeadE->ps.y, pHeadE->pe.y, pHeadE->ps.v, pHeadE->pe.v);
vb = Interpolation(pCurrentB->ScanLine, pHeadE->pNext->ps.y, pHeadE->pNext->pe.y, pHeadE->pNext->ps.v, pHeadE->pNext->pe.v);
BOOL Flag = FALSE;
double xb, xe;//扫描线和有效边相交区间的起点和终点坐标
for (pT1 = pHeadE; pT1 != NULL; pT1 = pT1->pNext)
{
if (Flag == FALSE)
{
xb = pT1->x;
Flag = TRUE;
}
else
{
xe = pT1->x;
for (double x = xb; x < xe; x++)//左闭右开
{
if (ROUND(x) + nWidth / 2 < nWidth && ROUND(x) + nWidth / 2 > 0 && pCurrentB->ScanLine + nHeight / 2 < nHeight && pCurrentB->ScanLine + nHeight / 2 > 0) {
zf = Interpolation(x, xb, xe, za, zb);
if (zf >= zBuffer[ROUND(x) + nWidth / 2][pCurrentB->ScanLine + nHeight / 2])//如果当前采样点的深度小于帧缓冲器中原采样点的深度)
{
tf = Interpolation(x, xb, xe, ta, tb);
tf.c = GetTextureColor(ROUND(tf.u), ROUND(tf.v), pTexture);
material->SetDiffuse(tf.c);//用纹理颜色作为材质的漫反射光反射率
material->SetAmbient(tf.c);//用纹理颜色作为材质的环境光反射率
if (light_type == 1) {
Cf = Interpolation(x, xb, xe, Ca, Cb);
}
else if (light_type == 2) {
vf = Interpolation(x, xb, xe, va, vb);
// 凹凸纹理
if (IfBump) {
CRGB frontU = GetTextureColor(ROUND(tf.u – 1), ROUND(tf.v), pTexture);
CRGB backU = GetTextureColor(ROUND(tf.u + 1), ROUND(tf.v), pTexture);
double Bu = (frontU.blue – backU.blue) / 2;
CRGB frontV = GetTextureColor(ROUND(tf.u), ROUND(tf.v – 1), pTexture);
CRGB backV = GetTextureColor(ROUND(tf.u), ROUND(tf.v + 1), pTexture);
double Bv = (frontV.blue – backV.blue) / 2;
CVector3 DNormal = CVector3(Bu, Bv, 0);
double BumpScale = 100.0;
vf = vf + BumpScale * DNormal;
vf = vf.Normalize();
}
Cf = light.Illuminate(EyePoint, CP3(ROUND(x), pCurrentB->ScanLine, zf), vf, material);
}
zBuffer[ROUND(x) + nWidth / 2][pCurrentB->ScanLine + nHeight / 2] = zf;//使用当前采样点的深度更新深度缓冲器
pDC->SetPixelV(ROUND(x), pCurrentB->ScanLine, RGB(Cf.red * 255, Cf.green * 255, Cf.blue * 255));
}
}
}
Flag = FALSE;
}
}
for (pT1 = pHeadE; pT1 != NULL; pT1 = pT1->pNext)//边的连续性
{
pT1->x = pT1->x + pT1->k;
}
}
}

void CFill::AddEt(CAET *pNewEdge)//合并ET表
{
CAET *pCE = pHeadE;
if (pCE == NULL)
{
pHeadE = pNewEdge;
pCE = pHeadE;
}
else
{
while (pCE->pNext != NULL)
{
pCE = pCE->pNext;
}
pCE->pNext = pNewEdge;
}
}
void CFill::EtOrder()//边表的冒泡排序算法
{
CAET *pT1 = NULL, *pT2 = NULL;
int Count = 1;
pT1 = pHeadE;
if (NULL == pT1)
{
return;
}
if (NULL == pT1->pNext)//如果该ET表没有再连ET表
{
return;//桶结点只有一条边，不需要排序
}
while (NULL != pT1->pNext)//统计结点的个数
{
Count++;
pT1 = pT1->pNext;
}
for (int i = 1; i < Count; i++)//冒泡排序
{
pT1 = pHeadE;
if (pT1->x > pT1->pNext->x)//按x由小到大排序
{
pT2 = pT1->pNext;
pT1->pNext = pT1->pNext->pNext;
pT2->pNext = pT1;
pHeadE = pT2;
}
else
{
if (pT1->x == pT1->pNext->x)
{
if (pT1->k > pT1->pNext->k)//按斜率由小到大排序
{
pT2 = pT1->pNext;
pT1->pNext = pT1->pNext->pNext;
pT2->pNext = pT1;
pHeadE = pT2;
}
}
}
pT1 = pHeadE;
while (pT1->pNext->pNext != NULL)
{
pT2 = pT1;
pT1 = pT1->pNext;
if (pT1->x > pT1->pNext->x)//按x由小到大排序
{
pT2->pNext = pT1->pNext;
pT1->pNext = pT1->pNext->pNext;
pT2->pNext->pNext = pT1;
pT1 = pT2->pNext;
}
else
{
if (pT1->x == pT1->pNext->x)
{
if (pT1->k > pT1->pNext->k)//按斜率由小到大排序
{
pT2->pNext = pT1->pNext;
pT1->pNext = pT1->pNext->pNext;
pT2->pNext->pNext = pT1;
pT1 = pT2->pNext;
}
}
}
}
}
}

CRGB CFill::GetTextureColor(double u, double v, CTexture* pTexture)//获得纹理颜色
{
/*检测图片的边界，防止越界*/
if (u < 0) u = 0; if (v < 0) v = 0;
if (u > pTexture->bmp.bmWidth – 1) u = pTexture->bmp.bmWidth – 1;
if (v > pTexture->bmp.bmHeight – 1) v = pTexture->bmp.bmHeight – 1;
/*查找对应纹理空间的颜色值*/
v = pTexture->bmp.bmHeight – 1 – v;
int position = ROUND(v * pTexture->bmp.bmWidthBytes + 4 * u);//颜色分量位置
COLORREF color = RGB(pTexture->image[position + 2], pTexture->image[position + 1], pTexture->image[position]);//获取颜色值
return CRGB(GetRValue(color) / 255.0, GetGValue(color) / 255.0, GetBValue(color) / 255.0);
}

CRGB CFill::Interpolation(double t, double t1, double t2, CRGB c1, CRGB c2)//颜色线性插值
{
CRGB c;
c = (t – t2) / (t1 – t2)*c1 + (t – t1) / (t2 – t1)*c2;
return c;
}

double CFill::Interpolation(double t, double t1, double t2, double z1,double z2)//线性插值
{
double z;
z = (t – t2) / (t1 – t2) * z1 + (t – t1) / (t2 – t1) * z2;
return z;
}

CVector3 CFill::Interpolation(double x, double xa, double xb, CVector3 va, CVector3 vb)//法向量线性插值
{
CVector3 v;
v = (x – xb) / (xa – xb) * va + (x – xa) / (xb – xa) * vb;
return v;
}

CT2 CFill::Interpolation(double m, double m0, double m1, CT2 T0, CT2 T1)//纹理地址线性插值
{
CT2 Texture;
Texture = (m1 – m) / (m1 – m0) * T0 + (m – m0) / (m1 – m0) * T1;
return Texture;
}

void CFill::ClearMemory()//安全删除所有桶和桶上面的边
{
DeleteAETChain(pHeadE);
CBucket *pBucket = pHeadB;
while (pBucket != NULL)// 针对每一个桶
{
CBucket * pBucketTemp = pBucket->pNext;
DeleteAETChain(pBucket->pET);
delete pBucket;
pBucket = pBucketTemp;
}
pHeadB = NULL;
pHeadE = NULL;
}

void CFill::DeleteAETChain(CAET* pAET)
{
while (pAET != NULL)
{
CAET* pAETTemp = pAET->pNext;
delete pAET;
pAET = pAETTemp;
}
}

void CFill::InitialDepthBuffer(int nWidth, int nHeight, double zDepth)//初始化深度缓冲
{
this->nWidth = nWidth, this->nHeight = nHeight;
zBuffer = new double* [nWidth];
for (int i = 0;i < nWidth;i++)
zBuffer[i] = new double[nHeight];
for (int i = 0;i < nWidth;i++)//初始化深度缓冲
for (int j = 0;j < nHeight;j++)
zBuffer[i][j] = zDepth;
}

————————————————
版权声明：本文为CSDN博主「仿生葡萄」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
原文链接：https://blog.csdn.net/uchiha_rin/article/details/122132213

本文为原创文章，转载请注明出处！

小波探测算法

目前针对声波远探测测井异常数据的识别主要方法为波形幅度法，该方法采用波形幅度值和预设的临界值对各个深度点的波形进行坏道识别，认为某深度点的波形幅度值大于临界值，则确定该深度点的波形为坏道。事实上，声波远探测可以简化视为，随着时间增加，其波形幅度出现正负交替特征，叠加后的平均幅度(或者幅度绝对值的平均)变化并不能客观反映坏道特征。为此，本发明采用基于平均幅度的总体标准差方法，进行水平井声波远探测测井资料坏道信号识别，从而实现资料质量控制。

基于平均幅度的总体标准差计算公式为：

式中Xi为一个波列中某个时刻的数据，n为波列中数据个数，QC_Ave为平均幅度，QC_Std为总体标准差。

根据实际情况设置识别阈值σ，当时则视为坏道信号。

如图5在5973.5m、5979m、5994m等深度点平均幅度QC_Ave没有明显变化，依靠QC_Ave不能有效指示坏道信号；而QC_Std在上述井段出现明显的指状尖峰特征，指示存在波形异常，判断为坏道信号，根据波形和频率分析结果（图5中第三道与第四道）已证实上述深度的存在声波远探测测井数据异常。sheng

本文为原创文章，转载请注明出处！

OpendTect基础模块

为用户提供了一个开放式的平台。用户可在平台上进行二次开发，增加新的模块和应用程序。基础模块中，OpendTect提供了以下功能：

1、数据的输入与输出：OpendTect可支持SEGY、LAS、ASCII格式的地震、钻井、测井、层位等常用数据的输入与输出。

2、地震反射层位的追踪与编辑：通过基础模块，可以在OpendTect系统中开展2D、3D地震层位的对比追踪。同时可以对其它解释软件所解释的地震反射层位进行编辑。

3、常规地震属性计算：在基础模块中可以通过等时、沿层等方式对“三瞬”、均方根振幅、平均频率、平均能量等常规地震属性进行计算、分析。

4、特有地震属性计算：OpendTect针对断层、小型褶皱及烃类异常等开发了如相似性、中值曲率、吸收系数等特有的地震属性。这些特有的地震属性较之常规属性不仅更具针对性，应用效果也得到改善。

5、频谱分解：应用频谱分解针对不同频率段提取地震属性已经被广泛应用。OpendTect在基础模块中也提供了频谱分解技术。可将时深域转换为频率域，开展不同频段地震属性提取。

6、三维可视化：在OpendTect基础平台上，可对地震数据体、层位、属性计算结果等进行立体、透视、透明等三维可视化显示。

双平方根算子地震波偏移

  1. void main()
  2. {
  3.   int b1, b2, b3, nz, ny, nx, i,j,k,kk, nshot, ndim, is,order,n;
  4.   float sz, sy, sx, o1, o2, o3, dz, dy, dx, vel0, br1, br2, br3,vel1,vel2,vel3;
  5.   float *v，*inface1,*inface2;
  7.   dz=10;
  8.   dy=20;
  9.   dx=20;
  10.   nz=80;
  11.   ny=100;
  12.   nz=100;
  13.   b1=1;
  14.   b2=1;
  15.   b3=1;
  16.   order=2;
  17.   vel1=800;
  18.   vel2=2000;
  19.   vel3=4000;
  20.   vel0=1.0/vel1;
  21.   n=nx*ny*nz;
  22.   v=new float[n];
  23.   inface1=new float[nx*ny];
  24.   inface2=new float[nx*ny];
  26.   for(i=0;i<n;i++)
  27.   {
  28.   v[i]=1.0/vel1;
  29.   }
  30.   for(i=0;i<nx*ny;i++)
  31.   {
  32.           inface1[i]=20+i*nz;
  33.   }
  34.   for(j=0;j<ny;j++)
  35.   {
  36.   for(i=0;i<nx;i++)
  37.   {
  38.                   int kk=int(700-3*i/20);
  39.                   inface2[j*nx+i]=kk;
  40.   for(k=20;k<kk;k++)
  41.   {
  42.   v[(i*ny+j)*nz+k]=1.0/vel2;
  43.   }
  44.   for(k=kk;k<nz;k++)
  45.   {
  46.   v[(i*ny+j)*nz+k]=1.0/vel3;
  47.   }
  48.   }
  49.   }
  50.   ifstream vel(“D:double_v”);
  51.   for(i=0;i<n;i++)
  52.   {
  53.   vel<<v[i]<<“t”;
  54.   }
  56.   int ngx,ngy,nsx,nsy;
  57.   float dgy,dgx,dsx,dsy;

本文为原创文章，转载请注明出处！