“好的,既然您想要正确的输出(没有扭曲),那我会选择射线投射器和着色器……但是为了不让您同时被C++和GLSL淹没,我决定在纯软件渲染中完成这个(从中转移到GLSL很容易)……”
“它的工作原理如下:”
循环遍历屏幕上的"所有"像素
您可以使用球体在屏幕上被剪裁后的期望BBOX来提高速度...
从相机焦点向实际屏幕像素投射光线
计算光线与球体之间的交点
您可以使用以下内容:
计算击中表面点和您的PCL之间的SDF(有符号距离函数)
为此,您需要将球坐标转换为笛卡尔坐标。如果您使用的是椭球体而不是球体,请参见:
根据SDF结果对像素进行着色
因此,如果距离为正,则表示像素位于区域外部。如果为零,则处于边缘,如果为负,则位于内部...
这是一个使用WGS84的小而简单的C++/VCL示例:
#include <vcl.h>
#include "GLSL_math.h"
#include "performance.h"
#pragma hdrstop
#include "win_main.h"
#pragma package(smart_init)
#pragma resource "*.dfm"
TMain *Main;
struct _pnt
{
float lat,lon,r;
vec3 p;
};
const int pnts=100;
_pnt pnt[pnts];
void pnt_init()
{
int i;
Randomize();
for (i=0;i<pnts;i++)
{
pnt[i].lon=(2.0*M_PI*Random()) ;
pnt[i].lat=( M_PI*Random())-(0.5*M_PI);
pnt[i].r=1000000.0;
}
}
const float _WGS84_a=6378137.00000;
const float _WGS84_b=6356752.31414;
const float _WGS84_e=8.1819190842622e-2;
const float _WGS84_aa=_WGS84_a*_WGS84_a;
const float _WGS84_bb=_WGS84_b*_WGS84_b;
const float _WGS84_ee=_WGS84_e*_WGS84_e;
const vec3 _WGS84_rr(1.0/_WGS84_aa,1.0/_WGS84_aa,1.0/_WGS84_bb);
bool WGS84_ray_intersection(vec3 p0,vec3 dp,float &l)
{
vec3 r;
float a,b,c,d,l0,l1;
l=-1.0;
r=_WGS84_rr;
a=(dp.x*dp.x*r.x)
+(dp.y*dp.y*r.y)
+(dp.z*dp.z*r.z);
b=(p0.x*dp.x*r.x)
+(p0.y*dp.y*r.y)
+(p0.z*dp.z*r.z); b*=2.0;
c=(p0.x*p0.x*r.x)
+(p0.y*p0.y*r.y)
+(p0.z*p0.z*r.z)-1.0;
d=((b*b)-(4.0*a*c));
if (d<0.0) return false;
d=sqrt(d);
a*=2.0;
l0=(-b+d)/a;
l1=(-b-d)/a;
if (fabs(l0)>fabs(l1)) { a=l0; l0=l1; l1=a; }
if (l0<0.0) { a=l0; l0=l1; l1=a; }
if (l0<0.0) return false;
l=l0;
return true;
}
void XYZtoWGS84(float &lon,float &lat,float &alt,float x,float y,float z)
{
int i;
float a,b,h,l,n,db,s;
a=atan2(y,x);
l=sqrt((x*x)+(y*y));
b=atan2(z,(1.0-_WGS84_ee)*l);
for (i=0;i<100;i++)
{
s=sin(b); db=b;
n=_WGS84_a/sqrt(1.0-(_WGS84_ee*s*s));
h=(l/cos(b))-n;
b=atan2(z,(1.0-divide(_WGS84_ee*n,n+h))*l);
db=fabs(db-b);
if (db<1e-12) break;
}
if (b>0.5*M_PI) b-=2.0*M_PI;
lon=a;
lat=b;
alt=h;
}
void WGS84toXYZ(float &x,float &y,float &z,float lon,float lat,float alt)
{
float a,b,h,l,c,s;
a=lon;
b=lat;
h=alt;
c=cos(b);
s=sin(b);
l=_WGS84_a/sqrt(1.0-(_WGS84_ee*s*s));
x=(l+h)*c*cos(a);
y=(l+h)*c*sin(a);
z=(((1.0-_WGS84_ee)*l)+h)*s;
}
void WGS84_draw(float yaw)
{
int **Pixels=Main->pyx;
int xs=Main->xs;
int ys=Main->ys;
float xs2=0.5*float(xs);
float ys2=0.5*float(ys);
int sx,sy;
vec3 O,X,Y,Z;
float focal_length;
vec3 p0,dp;
vec3 p;
float l,ll;
int i;
DWORD c=0x00204080;
l=3.0*_WGS84_a;
focal_length=(0.5*float(xs))/tan(30.0*M_PI/180.0);
O=vec3(-l*cos(yaw),0.0,-l*sin(yaw));
X=vec3( -sin(yaw),0.0, +cos(yaw));
Y=vec3( 0.0,1.0, 0.0);
Z=vec3( +cos(yaw),0.0, +sin(yaw));
p0=O;
for (i=0;i<pnts;i++)
{
float x,y,z;
WGS84toXYZ(x,y,z,pnt[i].lon,pnt[i].lat,0.0);
pnt[i].p=vec3(x,y,z);
}
for (sx=0;sx<xs;sx++)
for (sy=0;sy<ys;sy++)
{
p=vec3(sx,sy,focal_length);
p.x-=xs2;
p.y-=ys2;
dp.x=dot(p,X);
dp.y=dot(p,Y);
dp.z=dot(p,Z);
dp=normalize(dp);
if (WGS84_ray_intersection(p0,dp,l))
{
p=p0+l*dp;
l=1e100;
for (i=0;i<pnts;i++)
{
ll=length(p-pnt[i].p)-pnt[i].r;
if (l>ll) l=ll;
}
if (l> 0.0) c=0x00204080;
else if (l>-50000.0) c=0x00808000;
else c=0x00802020;
Pixels[sy][sx]=c;
}
}
}
void TMain::draw()
{
bmp->Canvas->Brush->Color=clBlack;
bmp->Canvas->FillRect(TRect(0,0,xs,ys));
tbeg();
static float yaw=0.0; yaw=fmod(yaw+0.1,2.0*M_PI);
WGS84_draw(yaw);
tend();
bmp->Canvas->Font->Color=clYellow;
bmp->Canvas->Brush->Style=bsClear;
bmp->Canvas->TextOutA(5,5,tstr());
bmp->Canvas->Brush->Style=bsSolid;
Main->Canvas->Draw(0,0,bmp);
}
__fastcall TMain::TMain(TComponent* Owner) : TForm(Owner)
{
bmp=new Graphics::TBitmap;
bmp->HandleType=bmDIB;
bmp->PixelFormat=pf32bit;
pyx=NULL;
pnt_init();
}
void __fastcall TMain::FormDestroy(TObject *Sender)
{
if (pyx) delete[] pyx;
delete bmp;
}
void __fastcall TMain::FormResize(TObject *Sender)
{
xs=ClientWidth; xs2=xs>>1;
ys=ClientHeight; ys2=ys>>1;
bmp->Width=xs;
bmp->Height=ys;
if (pyx) delete pyx;
pyx=new int*[ys];
for (int y=0;y<ys;y++) pyx[y]=(int*) bmp->ScanLine[y];
draw();
}
void __fastcall TMain::FormPaint(TObject *Sender)
{
draw();
}
void __fastcall TMain::Timer1Timer(TObject *Sender)
{
draw();
}
只需忽略VCL的东西,将像素渲染转换为您的环境。请注意,像素访问通常非常慢,因此我会将其渲染到位图上,而不是作为32位无符号整数指针进行访问,当整个过程完成时,我再渲染位图...有关更多信息,请参见:
唯一重要的是
WGS84_draw
和它正在使用的子函数。
我使用的唯一非标准库是我的
GLSL_math.h,它模仿了GLSL数学,并可与任何矢量库(如GLM)一起使用...
这里是预览:
正如您所看到的,没有任何扭曲...
如果您不想像这样进行射线投射,则必须将您的点转换为覆盖其半径和球面曲率的三角形多边形,并以此方式呈现它们:
- 使用边缘颜色呈现所有圆
- 使用内部颜色呈现所有圆,但将半径减小边缘宽度
如果您不需要内部,则只需拥有边缘多边形并一次渲染即可...
[编辑1]我设法让此版本的着色器工作:
CPU端代码:
#include <vcl.h>
#include "GLSL_math.h"
#pragma hdrstop
#include "Unit1.h"
#include "gl_simple.h"
#pragma package(smart_init)
#pragma resource "*.dfm"
TForm1 *Form1;
const float _WGS84_a=6378137.00000;
const float _WGS84_b=6356752.31414;
const float _WGS84_e=8.1819190842622e-2;
const float _WGS84_aa=_WGS84_a*_WGS84_a;
const float _WGS84_bb=_WGS84_b*_WGS84_b;
const float _WGS84_ee=_WGS84_e*_WGS84_e;
const vec3 _WGS84_rr(1.0/_WGS84_aa,1.0/_WGS84_aa,1.0/_WGS84_bb);
void XYZtoWGS84(float &lon,float &lat,float &alt,vec3 p)
{
int i;
float a,b,h,l,n,db,s;
a=atan2(p.y,p.x);
l=sqrt((p.x*p.x)+(p.y*p.y));
b=atan2(p.z,(1.0-_WGS84_ee)*l);
for (i=0;i<100;i++)
{
s=sin(b); db=b;
n=_WGS84_a/sqrt(1.0-(_WGS84_ee*s*s));
h=(l/cos(b))-n;
b=atan2(p.z,(1.0-divide(_WGS84_ee*n,n+h))*l);
db=fabs(db-b);
if (db<1e-12) break;
}
if (b>0.5*M_PI) b-=2.0*M_PI;
lon=a;
lat=b;
alt=h;
}
void WGS84toXYZ(float &x,float &y,float &z,float lon,float lat,float alt)
{
float a,b,h,l,c,s;
a=lon;
b=lat;
h=alt;
c=cos(b);
s=sin(b);
l=_WGS84_a/sqrt(1.0-(_WGS84_ee*s*s));
x=(l+h)*c*cos(a);
y=(l+h)*c*sin(a);
z=(((1.0-_WGS84_ee)*l)+h)*s;
}
const int pnts=1000;
struct _pnt{ float x,y,z,r; };
_pnt pnt[pnts];
void pnt_init()
{
int i;
Randomize();
for (i=0;i<pnts;i++)
{
float lat,lon,r,x,y,z;
lon=(2.0*M_PI*Random()) ;
lat=( M_PI*Random())-(0.5*M_PI);
r=250000.0+250000.0*Random();
WGS84toXYZ(x,y,z,lon,lat,0.0);
pnt[i].x=x;
pnt[i].y=y;
pnt[i].z=z;
pnt[i].r=r;
}
}
void gl_draw()
{
int i;
static float yaw=0.0; yaw=fmod(yaw+0.05,2.0*M_PI);
float l=3.0*_WGS84_a;
float focal_length=(1.0)/tan(30.0*M_PI/180.0);
float tm_eye[16];
tm_eye[ 0]=-sin(yaw);
tm_eye[ 1]=0.0;
tm_eye[ 2]=+cos(yaw);
tm_eye[ 3]=0.0;
tm_eye[ 4]=0.0;
tm_eye[ 5]=1.0;
tm_eye[ 6]=0.0;
tm_eye[ 7]=0.0;
tm_eye[ 8]=+cos(yaw);
tm_eye[ 9]=0.0;
tm_eye[10]=+sin(yaw);
tm_eye[11]=0.0;
tm_eye[12]=-l*cos(yaw);
tm_eye[13]=0.0;
tm_eye[14]=-l*sin(yaw);
tm_eye[15]=0.0;
glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
glMatrixMode(GL_PROJECTION);
glLoadIdentity();
glMatrixMode(GL_MODELVIEW);
glLoadIdentity();
glDisable(GL_DEPTH_TEST);
glDisable(GL_CULL_FACE);
glUseProgram(prog_id);
i=glGetUniformLocation(prog_id,"pnts" ); glUniform1i(i,pnts);
i=glGetUniformLocation(prog_id,"pnt" ); glUniform4fv(i,pnts,(float*)(pnt));
i=glGetUniformLocation(prog_id,"_WGS84_rr" ); glUniform3fv(i,1,_WGS84_rr.dat);
i=glGetUniformLocation(prog_id,"width"); glUniform1f(i,50000.0);
i=glGetUniformLocation(prog_id,"focal_length"); glUniform1f(i,focal_length);
i=glGetUniformLocation(prog_id,"tm_eye" ); glUniformMatrix4fv(i,1,false,tm_eye);
glBegin(GL_QUADS);
glVertex2f(-1.0,-1.0);
glVertex2f(-1.0,+1.0);
glVertex2f(+1.0,+1.0);
glVertex2f(+1.0,-1.0);
glEnd();
glUseProgram(0);
glFlush();
SwapBuffers(hdc);
}
__fastcall TForm1::TForm1(TComponent* Owner):TForm(Owner)
{
gl_init(Handle);
int hnd,siz; char vertex[4096],geometry[4096],fragment[4096];
hnd=FileOpen("WGS84_region.glsl_vert",fmOpenRead); siz=FileSeek(hnd,0,2); FileSeek(hnd,0,0); FileRead(hnd,vertex ,siz); vertex [siz]=0; FileClose(hnd);
hnd=FileOpen("WGS84_region.glsl_frag",fmOpenRead); siz=FileSeek(hnd,0,2); FileSeek(hnd,0,0); FileRead(hnd,fragment,siz); fragment[siz]=0; FileClose(hnd);
glsl_init(vertex,NULL,fragment);
mm_log->Text=glsl_log;
pnt_init();
}
void __fastcall TForm1::FormDestroy(TObject *Sender)
{
gl_exit();
glsl_exit();
}
void __fastcall TForm1::FormResize(TObject *Sender)
{
gl_resize(ClientWidth-mm_log->Width,ClientHeight);
}
void __fastcall TForm1::FormPaint(TObject *Sender)
{
gl_draw();
}
void __fastcall TForm1::Timer1Timer(TObject *Sender)
{
gl_draw();
}
Vertex:
#version 400 core
layout(location = 0) in vec2 pos;
uniform mat4 tm_eye;
uniform float focal_length;
out vec3 ray_p0;
out vec3 ray_dp;
void main()
{
vec3 p;
ray_p0=tm_eye[3].xyz;
ray_dp=normalize((tm_eye*vec4(pos,focal_length,0.0)).xyz);
gl_Position=vec4(pos,0.0,1.0);
}
Fragment:
#version 400 core
uniform int pnts;
uniform vec4 pnt[1000];
uniform vec3 _WGS84_rr;
uniform float width;
in vec3 ray_p0;
in vec3 ray_dp;
out vec4 col;
float ray_l=1e30;
bool WGS84_ray_intersection(vec3 p0,vec3 dp)
{
vec3 r;
float a,b,c,d,l0,l1;
ray_l=-1.0;
r=_WGS84_rr;
a=(dp.x*dp.x*r.x)
+(dp.y*dp.y*r.y)
+(dp.z*dp.z*r.z);
b=(p0.x*dp.x*r.x)
+(p0.y*dp.y*r.y)
+(p0.z*dp.z*r.z); b*=2.0;
c=(p0.x*p0.x*r.x)
+(p0.y*p0.y*r.y)
+(p0.z*p0.z*r.z)-1.0;
d=((b*b)-(4.0*a*c));
if (d<0.0) return false;
d=sqrt(d);
a*=2.0;
l0=(-b+d)/a;
l1=(-b-d)/a;
if (abs(l0)>abs(l1)) { a=l0; l0=l1; l1=a; }
if (l0<0.0) { a=l0; l0=l1; l1=a; }
if (l0<0.0) return false;
ray_l=l0;
return true;
}
void main()
{
int i;
vec3 p,c;
float l,ll;
if (!WGS84_ray_intersection(ray_p0,ray_dp)) discard;
p=ray_p0+ray_l*ray_dp;
l=1e30;
for (i=0;i<pnts;i++)
{
ll=length(p-pnt[i].xyz)-pnt[i].w;
if (l>ll) l=ll;
}
if (l> 0.0) c=vec3(0.1,0.3,0.5);
else if (l>-width) c=vec3(0.1,0.8,0.9);
else c=vec3(0.3,0.3,0.3);
col=vec4(c,1.0);
}
预览(1000分):
然而,由于我的图形处理实现限制,我无法通过1021点。因此,如果你需要更多的话,你必须像之前提到的那样使用纹理(不要忘记使用非夹紧像素格式,如
GL_LUMINANCE32F_ARB
...)。
在我的设置中,速度和分辨率似乎都不是问题(nVidia GTX 550 Ti)。
另外,请注意,上面的示例尚未包含任何纵横比校正...
r
对于所有点都是相同的吗?你意识到为什么在 n==2 的情况下,KD 树是相当无用/误导的吗? - Mad Physicist