Mesh.pas

Top  Previous  Next

unit Mesh;

 

interface

 

uses

Windows, Messages, Classes, Graphics, Forms, ExtCtrls, Menus, Math,

Controls, Dialogs, SysUtils, OpenGL,TextureGL;

 

Type

   TGeo3DMesh = class;

 

       TGeoRect = record

          dLeft, dTop, dRight, dBottom : double;

   end;

 

   PGeo3DPt = ^TGeo3DPt;

   TGeo3DPt = record

          dX, dY,dZ : glfloat;

   end;

 

   TGeo3DPtArray = array of TGeo3DPt;

   //四叉树方式的点树结构

   TGeoMeshNode = class

   public

          Mesh : TGeo3DMesh;

          nLevel : integer;                                //第几级

 

       nXWidth, nYWidth : integer;

          Pts : TList;

       dOffset : single;                                //高程值的判断阈值

 

       SubNodes : TList;

       constructor Create(_Mesh : TGeo3DMesh);

       destructor Destroy;override;

       procedure Draw(MODE: glenum);

       procedure CalcSubRoot;

   end;

 

   TGLNrm=array[0..2] of glfloat;

 

   TGLvrt=record

     crd:TGLNrm;

     nrm:TGLNrm;

     clr:TGLNrm;

     txt:array[0..1] of glfloat;

   end;

 

   TGLfcs=record

     vrt:array[0..2]of word;

     nrm:TGLNrm;

   end;

 

       TMeshConfig = record

       smt:boolean;

       clr:boolean;

       txt:boolean;

       mode:glenum;                //图形类型,线、点,多边形等

       fa : Boolean;                //是否显示法向量

               nCellCount : integer;

       dBaseHeight : single;

     end;

//三维模型核心

   TGeo3DMesh = class

   public

     Base,TexBmp : TBitmap;

     //所有网格上的点

        Pts : Array of TGeo3DPt;

     nXWidth, nYWidth : integer;                                //水平方向的网格数,缺省为125*125

     Root : TGeoMeshNode;

     procedure GetDataFromVtrs;

     procedure ClearRootPts;

       public

     Texture : TTextureGl;

 

     MeshConfig : TMeshConfig;

     cVrt:longint;                                             //暑. 怵.

     cFcs:longint;                                            //暑. 沭.

     Vrts:array of TGLVrt;                                         //锑耨桠 忮瘌桧

     Fcss:array of TGLFcs;                                         //锑耨桠 沭囗彘

 

     MeshRect : TGeoRect;

     MinH, MaxH : glfloat;                                           // 记录当前范围的最大值和最小值

     ShowMinRange, ShowMaxRange : glfloat;     // 记录显示范围的[0, 1]

     bShowRange : Boolean;                     // 是否只显示设定范围内的信息

 

     constructor Create(BB, TB : string);

     destructor Destroy;override;

        procedure CalcRect;

     procedure calc_normals_fr;

     procedure calc_normals_sm;

 

     procedure filter_rs(x,y,z:glfloat);

     procedure filter_sd(x,y,z:glfloat);

 

     procedure load_lst_obj(var f:textfile);

     procedure load_gms_obj(var f:textfile);

     procedure Save_Lst_obj(var f:textfile);

     procedure LoadTexture(BMP : TBitmap);

     procedure loadfromfile_lst(filename:string);

     procedure loadfromfile_gms(filename:string);

     procedure Save_To_file_lst(filename:string);

       //从位图得到Mesh的高程和颜色值

     procedure CreateMeshFromBitmap(nWidth,nheight,mash:glfloat;cpw,cph:byte;var b:tbitmap);overload;

     procedure CreateMeshFromBitmap(b:tbitmap);overload;

 

     function gt_ln(x1,y1,z1,x2,y2,z2:real;var code:boolean):tGlnrm;

 

     procedure draw;

     Procedure CalRange;

     //建立四叉树

     procedure BuildTree;

   end;

 

   function get_Normal_fl(var p1,p2,p3:TGLNrm):TGLNrm;

   function get_j(r,g,b:glfloat):glfloat;

   function PixelInOtr(x1,x2,x:glfloat):boolean;

   function get_dl_line(x1,y1,z1,x2,y2,z2:glfloat):glfloat;

   function g_ang(x,y:glfloat):real;

   function pt_in_tr(x1,y1,x2,y2,x3,y3,x,y:glfloat):boolean;

   function getpoint(p1,p2,pt1:TGLnrm;nrm:TGLnrm):TGLnrm;

   procedure rotate_point(angle:glfloat;var x,y:glfloat);

   Procedure GetCoords(x1,y1,x2,y2,x3,y3,x4,y4:real; Var x,y : glfloat; code:boolean);

 

implementation

 

function getNormal(x1,y1,z1,x2,y2,z2,x3,y3,z3 : glfloat):TGLArrayf3;

var

wrki, vx1, vy1, vz1, vx2, vy2, vz2 : GLfloat;

nx, ny, nz : GLfloat;

wrkVector : tpoint;

f:textfile;

begin

    vx1 := x1 - x2;

    vy1 := y1 - y2;

    vz1 := z1 - z2;

    vx2 := x2 - x3;

    vy2 := y2 - y3;

    vz2 := z2 - z3;

 

    nx := vy1 * vz2 - vz1 * vy2;

    ny := vz1 * vx2 - vx1 * vz2;

    nz := vx1 * vy2 - vy1 * vx2;

    wrki := sqrt (nx * nx + ny * ny + nz * nz);

 

    If wrki = 0 then wrki := 1;

 

       result[0] := nx/ wrki;

       result[1] := ny/ wrki;

       result[2] := nz/ wrki;

end;

{var a,b:TGLArrayf3;

begin

  //make two vectors

  a[0]:=x2 - x1; a[1]:=y2-y1; a[2]:= z2-z1;

  b[0]:=x3 - x1; b[1]:=y3-y1; b[2] := z3 - z1;

  //calculate cross-product

  result[0]:=a[1]*b[2]-a[2]*b[1];

  result[1]:=a[2]*b[0]-a[0]*b[2];

  result[2]:=a[0]*b[1]-a[1]*b[0];

end;

}

 

 

procedure TGeo3DMesh.filter_rs(x,y,z:glfloat);

var i:integer;

begin

for i:=0 to cVrt do

begin

   Vrts[i].crd[0]:=Vrts[i].crd[0]/x;

   Vrts[i].crd[1]:=Vrts[i].crd[1]/y;

   Vrts[i].crd[2]:=Vrts[i].crd[2]/z;

end;

calc_normals_fr;

calc_normals_sm;

end;

 

procedure TGeo3DMesh.filter_sd(x,y,z:glfloat);

var i:integer;

begin

       for i:=0 to cVrt do begin

           Vrts[i].crd[0]:=Vrts[i].crd[0]+x;

           Vrts[i].crd[1]:=Vrts[i].crd[1]+y;

           Vrts[i].crd[2]:=Vrts[i].crd[2]+z;

       end;

end;

 

procedure TGeo3DMesh.Save_To_file_lst(filename:string);

var f:textfile;

begin

assignfile(f,filename);

rewrite(f);

Save_Lst_obj(f);

closefile(f);

end;

 

procedure TGeo3DMesh.load_lst_obj(var f:textfile);

var i:integer;

begin

 

readln(f,cVrt);

GetMem(Vrts,cVrt*SizeOf(TGLVrt));

for i:=0 to cVrt-1 do

   readln(f,Vrts[i].crd[0],Vrts[i].crd[1],Vrts[i].crd[2],

            Vrts[i].clr[0],Vrts[i].clr[1],Vrts[i].clr[2],

            Vrts[i].txt[0],Vrts[i].txt[1]);

 

readln(f,cFcs);

GetMem(Fcss,cFcs*SizeOf(TGLFcs));

for i:=0 to cFcs-1 do

   readln(f,Fcss[i].vrt[0],Fcss[i].vrt[1],Fcss[i].vrt[2]);

 

calc_normals_fr;

calc_normals_sm;

 

end;

 

procedure TGeo3DMesh.Save_Lst_obj(var f:textfile);

var i:integer;

   j:byte;

begin

writeln(f,'{');

writeln(f,'{');

writeln(f,0,' ',1,' ',4);

for i:=0 to cVrt-1 do

   writeln(f,Vrts[i].crd[0]:15:4,Vrts[i].crd[2]:15:4,Vrts[i].crd[1]:15:4);

writeln(f,cVrt:10);

for i:=0 to cVrt-1 do

   writeln(f,Vrts[i].txt[0]:15:4,' ',Vrts[i].txt[1]:15:4);

writeln(f,cVrt:10);

for i:=0 to cVrt-1 do

   writeln(f,Vrts[i].clr[0]:15:4,' ',Vrts[i].clr[1]:15:4,' ',Vrts[i].clr[2]:15:4,' ',1);

 

writeln(f,cFcs:10);

for i:=0 to cFcs-1 do

   writeln(f,Fcss[i].vrt[0]+1:6,' ',Fcss[i].vrt[2]+1:6,' ',Fcss[i].vrt[1]+1:6,' ',

             Fcss[i].vrt[0]+1:6,' ',Fcss[i].vrt[2]+1:6,' ',Fcss[i].vrt[1]+1:6,' ',

             Fcss[i].vrt[0]+1:6,' ',Fcss[i].vrt[2]+1:6,' ',Fcss[i].vrt[1]+1:6);

writeln(f,'}');

end;

 

procedure TGeo3DMesh.load_gms_obj(var f:textfile);

var s:string;

   i,j:integer;

   c:glfloat;

begin

readln(f,s);

 

if s='TriMesh()' then

begin

   readln(f,s);

 

   readln(f,i,j); cVrt:=i; cFcs:=j;

 

   GetMem(Vrts,cVrt*SizeOf(TGLVrt));

   GetMem(Fcss,cFcs*SizeOf(TGLFcs));

 

   readln(f,s);  //忮瘌桧

   for i:=0 to cVrt-1 do

   begin

     read(f,c);   vrts[i].crd[0]:=c/100;

     read(f,c);   vrts[i].crd[1]:=c/100;

     readln(f,c); vrts[i].crd[2]:=c/100;

   end;

   readln(f,s);

 

   readln(f,s);  //沭囗

   for i:=0 to cFcs-1 do

   begin

     read(f,j);   fcss[i].vrt[0]:=j-1;

     read(f,j);   fcss[i].vrt[1]:=j-1;

     readln(f,j); fcss[i].vrt[2]:=j-1;

   end;

   readln(f,s);

 

   readln(f,s);  //沭. 眇.

   for i:=0 to cFcs-1 do

   begin

     read(f,c);   //fcss[i].nrm[0]:=c;

     read(f,c);   //fcss[i].nrm[1]:=c;

     readln(f,c); //fcss[i].nrm[2]:=c;

   end;

   readln(f,s);

 

   readln(f,s);  //胥. 眇.

   for i:=0 to cVrt-1 do

   begin

     read(f,c);   //vrts[i].nrm[0]:=c;

     read(f,c);   //vrts[i].nrm[1]:=c;

     readln(f,c); //vrts[i].nrm[2]:=c;

   end;

   readln(f,s);

 

   readln(f,s);

 

   calc_normals_fr;

   calc_normals_sm;

 

end;

 

end;

 

 

function TGeo3DMesh.gt_ln(x1,y1,z1,x2,y2,z2:real;var code:boolean):tGlnrm;

var i:integer;

   p1,p2,p:TGLnrm;

begin

code:=false;

p1[0]:=x1; p1[1]:=y1; p1[2]:=z1;

p2[0]:=x2; p2[1]:=y2; p2[2]:=z2;

for i:=0 to cfcs-1 do

begin

    p:=getpoint(p1,p2,vrts[fcss[i].vrt[0]].crd,fcss[i].nrm);

    if PixelInOtr(p1[1],p2[1],p[1]) then

    if pt_in_tr(vrts[fcss[i].vrt[0]].crd[0],vrts[fcss[i].vrt[0]].crd[2],

                vrts[fcss[i].vrt[1]].crd[0],vrts[fcss[i].vrt[1]].crd[2],

                vrts[fcss[i].vrt[2]].crd[0],vrts[fcss[i].vrt[2]].crd[2],

                p[0],p[2]) then

    begin

       gt_ln:=p;

       code:=true;

       break;

    end;

end;

end;

 

 

procedure TGeo3DMesh.loadfromfile_gms(filename:string);

var f:textfile;

   s:string;

label m1;

begin

assignfile(f,filename);

reset(f);

repeat

   ReadLn(f, S);

   if s='New object' then

   begin

    load_gms_obj(f);

    goto m1;

   end;

until (eof(f));

m1:

closefile(f);

end;

 

 

procedure TGeo3DMesh.loadfromfile_lst(filename:string);

var f:textfile;

   s:string;

label m1;

begin

assignfile(f,filename);

reset(f);

load_lst_obj(f);

closefile(f);

end;

 

 

procedure TGeo3DMesh.calc_normals_fr;

var

       i:integer;

begin

   with MeshRect do begin

          dLeft := MaxDouble;

       dTop := MaxDouble;

       dRight := -MaxDouble;

       dBottom := -MaxDouble;

   end;

 

         for i:=0 to cFcs-1 do begin

            fcss[i].nrm := get_Normal_fl(vrts[fcss[i].vrt[0]].crd,  vrts[fcss[i].vrt[1]].crd,vrts[fcss[i].vrt[2]].crd);

       if fcss[i].nrm[0] < MeshRect.dLeft then MeshRect.dLeft := fcss[i].nrm[0];

       if fcss[i].nrm[0] > MeshRect.dRight then MeshRect.dRight := fcss[i].nrm[0];

       if fcss[i].nrm[2] < MeshRect.dTop then MeshRect.dTop := fcss[i].nrm[2];

       if fcss[i].nrm[2] > MeshRect.dBottom then MeshRect.dBottom := fcss[i].nrm[2];

       end;

end;

 

procedure TGeo3DMesh.calc_normals_sm;

var i,j:integer;

   n:array of word;

begin

       getmem(n,sizeof(word)*cVrt);

 

       for i:=0 to cVrt-1 do begin

          Vrts[i].nrm[0]:=0;

          Vrts[i].nrm[1]:=0;

          Vrts[i].nrm[2]:=0;

          n[i]:=0;

       end;

 

       for i:=0 to cFcs-1 do begin

                 for j:=0 to 2 do begin

                  Vrts[Fcss[i].vrt[j]].nrm[0]:=(Vrts[Fcss[i].vrt[j]].nrm[0]+Fcss[i].nrm[0]);

                  Vrts[Fcss[i].vrt[j]].nrm[1]:=(Vrts[Fcss[i].vrt[j]].nrm[1]+Fcss[i].nrm[1]);

                  Vrts[Fcss[i].vrt[j]].nrm[2]:=(Vrts[Fcss[i].vrt[j]].nrm[2]+Fcss[i].nrm[2]);

                  n[Fcss[i].vrt[j]]:=n[Fcss[i].vrt[j]]+1;

                 end;

       end;

 

       for i:=0 to cVrt-1 do begin

           Vrts[i].nrm[0]:=Vrts[i].nrm[0]/n[i];

           Vrts[i].nrm[1]:=Vrts[i].nrm[1]/n[i];

           Vrts[i].nrm[2]:=Vrts[i].nrm[2]/n[i];

       end;

 

       freemem(n);

end;

 

procedure TGeo3DMesh.draw;

var

       i : integer;

   dY, dRange : glfloat;

   dR : Double;

       n1, n2: TGLArrayf3;   //normals

begin

       //先在比DEM低一些的地方把网格画出来,上面铺地图位图,再把DEM画在空中

       if Assigned(Texture) then Texture.Enable;

 

   //

{         glColor3f(0.85, 1, 0.85);

       glEnable(GL_BLEND);

       glEnable(GL_TEXTURE_GEN_S);

       glEnable(GL_TEXTURE_GEN_R);

 

//         glEnable(GL_TEXTURE_GEN_T);

}

   if Assigned(Root) then begin

               Root.Draw(MeshConfig.mode);

   end

   else begin

       for I := 0 to (nXWidth - 1) * nYWidth - 1 do begin

           if I mod nXWidth <> nXWidth - 1 then begin

               glBegin(MeshConfig.mode);

                   n1 := GetNormal(pts[I].dX, pts[I].dY, Pts[I].dZ, pts[I + nXWidth].dX, pts[I + nXWidth].dY, Pts[I+ nXWidth].dZ ,pts[I + 1].dX, pts[I + 1].dY, Pts[I + 1].dZ);

   //                   n1 := GetNormal(pts[I].dX, pts[I].dY, Pts[I].dZ, pts[I + 1].dX, pts[I + 1].dY, Pts[I+ 1].dZ ,pts[I + nXWidth].dX, pts[I + nXWidth].dY, Pts[I + nXWidth].dZ);

                   glNormal3fv(@n1);

                   glTexCoord2f((I mod nYWidth) / nXWidth, (I div nXWidth) / nYWidth);

                   glVertex3f(pts[I].dX, pts[i].dy, pts[i].dZ);

                   glTexCoord2f(((I + nXWidth) mod nXWidth) / nXWidth, ((I + nXWidth) div nYWidth) / nYWidth);

                   glVertex3f(pts[I + nXWidth].dX, pts[I + nXWidth].dy, pts[I + nXWidth].dZ);

                   glTexCoord2f(((I + 1) mod nXWidth) / nXWidth, ((I + 1) div nYWidth) / nYWidth);

                   glVertex3f(pts[I + 1].dX, pts[I + 1].dy, pts[I + 1].dZ);

               glEnd;

 

               glBegin(MeshConfig.mode);

                   n2 := GetNormal(pts[I + nXWidth + 1].dX, pts[I + nXWidth + 1].dY, Pts[I + nXWidth + 1].dZ,pts[I + nXWidth].dX, pts[I + nXWidth].dY, Pts[I + nXWidth].dZ ,pts[I + 1].dX, pts[I + 1].dY, Pts[I + 1].dZ);

   //                   n2 := GetNormal(pts[I + 1].dX, pts[I + 1].dY, Pts[I + 1].dZ,pts[I + nXWidth].dX, pts[I+ nXWidth ].dY, Pts[I+ nXWidth ].dZ,pts[I  + nXWidth +1].dX, pts[I + nXWidth + 1].dY, Pts[I + nXWidth + 1].dZ);

   //                   glNormal3fv(@n2);

                   glTexCoord2f(((I + 1) mod nXWidth) / nXWidth, ((I + 1) div nYWidth) / nYWidth);

                   glVertex3f(pts[I + 1].dX, pts[i + 1].dy, pts[i + 1].dZ);

                   glTexCoord2f(((I + nXWidth) mod nXWidth) / nXWidth, ((I + nXWidth) div nYWidth) / nYWidth);

                   glVertex3f(pts[I + nXWidth].dX, pts[I + nXWidth].dy, pts[I + nXWidth].dZ);

                   glTexCoord2f(((I + nXWidth + 1) mod nXWidth) / nXWidth, ((I + nXWidth + 1) div nYWidth) / nYWidth);

                   glVertex3f(pts[I + nXWidth + 1].dX, pts[I + nXWidth + 1].dy, pts[I + nXWidth + 1].dZ);

               glEnd;

           end;

       end;

   end;

end;

 

procedure TGeo3DMesh.CreateMeshFromBitmap(nWidth,nHeight,mash:glfloat;cpw,cph:byte;var b:tbitmap);

var i,j,x,y:integer;

   wp,hp:glfloat;

   wb,hb:glfloat;

   ceny:glfloat;

begin

       self.MeshConfig.nCellCount := cpw;

 

       cvrt:=cpw*cph;  cFcs:=((cpw-1)*(cph-1))*2;

       getmem(vrts,sizeof(TGlVrt)*cvrt);

       getmem(fcss,sizeof(TGlFcs)*cfcs);

 

       x:=0;

       wp := nWidth / cpw; hp := nHeight/cph;

 

       wb := 1;//b.Width/cpw;

       hb := 1;//b.height/cph;

 

       ceny:=0;

       for i:=0 to cpw-1 do begin

          for j:=0 to cph-1 do begin

 

     vrts[x].clr[0]:=GetRValue(B.Canvas.Pixels[round(i*wb),round(j*hb)])/255;

     vrts[x].clr[1]:=GetgValue(B.Canvas.Pixels[round(i*wb),round(j*hb)])/255;

     vrts[x].clr[2]:=GetbValue(B.Canvas.Pixels[round(i*wb),round(j*hb)])/255;

 

     vrts[x].txt[0]:=(i*wb)/b.Width;

     vrts[x].txt[1]:=(j*hb)/b.height;

 

     vrts[x].crd[0]:=(i*wp)-(nWidth/2);

     vrts[x].crd[1]:=get_j(vrts[x].clr[0]*255

                          ,vrts[x].clr[1]*255

                          ,vrts[x].clr[2]*255)/mash;

     vrts[x].crd[2]:=(j*hp)-(nWidth/2);

 

     ceny:=ceny+vrts[x].crd[1];

     x:=x+1;

   end;

end;

ceny:=ceny/cvrt;

for i:=0 to cvrt-1 do vrts[i].crd[1]:=vrts[i].crd[1]-ceny;

 

x:=0;

for i:=0 to cpw-2 do

begin

   for j:=1 to cph-1 do

   begin

     fcss[x].vrt[0]:=(i*cpw)+j-1;

     fcss[x].vrt[1]:=(i*cpw)+j;

     fcss[x].vrt[2]:=( (i+1)*cpw )+j-1;

 

     fcss[x+1].vrt[0]:=( (i+1)*cpw )+j;

     fcss[x+1].vrt[1]:=( (i+1)*cpw )+j-1;

     fcss[x+1].vrt[2]:=(i*cpw)+j;

 

     x:=x+2;

   end;

end;

calc_normals_fr;

//  ShowMessage(FloatToStr(meshrect.dLeft) + ',' + floattostr(meshrect.dTop) + ',' + floattostr(meshrect.dRight) + ',' + floattostr(meshrect.dBottom)); 

calc_normals_sm;

end;

 

function get_Normal_fl(var p1,p2,p3:TGLNrm):TGLNrm;

var

wrki, vx1, vy1, vz1, vx2, vy2, vz2 : GLfloat;

nx, ny, nz : GLfloat;

wrkVector : tpoint;

f:textfile;

begin

    vx1 := p1[0] - p2[0];

    vy1 := p1[1] - p2[1];

    vz1 := p1[2] - p2[2];

    vx2 := p2[0] - p3[0];

    vy2 := p2[1] - p3[1];

    vz2 := p2[2] - p3[2];

 

    nx := vy1 * vz2 - vz1 * vy2;

    ny := vz1 * vx2 - vx1 * vz2;

    nz := vx1 * vy2 - vy1 * vx2;

    wrki := sqrt (nx * nx + ny * ny + nz * nz);

 

    If wrki = 0 then wrki := 1;

 

    get_Normal_fl[0] := nx/ wrki;

    get_Normal_fl[1] := ny/ wrki;

    get_Normal_fl[2] := nz/ wrki;

end;

 

function get_j(r,g,b:glfloat):glfloat;

begin

get_j:=(0.59*g)+(0.3*r)+(0.11*b);

end;

 

function getpoint(p1,p2,pt1:TGLnrm;nrm:TGLnrm):TGLnrm;

var p,l,m,n:glfloat;

begin

p1[0]:=p1[0]-pt1[0];  p1[1]:=p1[1]-pt1[1];  p1[2]:=p1[2]-pt1[2];

 

l:=(p2[0]-pt1[0])-p1[0];

m:=(p2[1]-pt1[1])-p1[1];

n:=(p2[2]-pt1[2])-p1[2];

 

p:=((nrm[0]*p1[0])+(nrm[1]*p1[1])+(nrm[2]*p1[2]))

/((nrm[0]*l)+(nrm[1]*m)+(nrm[2]*n));

 

getpoint[0]:=(p1[0]-(l*p))+pt1[0];

getpoint[1]:=(p1[1]-(m*p))+pt1[1];

getpoint[2]:=(p1[2]-(n*p))+pt1[2];

end;

 

Procedure GetCoords(x1,y1,x2,y2,x3,y3,x4,y4:real; Var

                   x,y : glfloat; code:boolean);

Var k1,k2 : extended;

   b1,b2 : extended;

  xx, yy : extended;

Begin

  If (x1-x2=0) and (x3-x4=0) Then

  Begin

     code:=false;

     Exit;

  End;

  code:=true;

  If (x1-x2=0) and (x3<>x4) Then

  Begin

     k2:=(y3-y4)/(x3-x4);

     b2:=y3-k2*x3;

     x:=x1;

     y:=k2*x+b2;

     Exit;

  End;

  If (x3-x4=0) and (x1<>x2) Then

  Begin

     k1:=(y1-y2)/(x1-x2);

     b1:=y1-k1*x1;

     x:=x3;

     y:=k1*x+b1;

     Exit;

  End;

  If (y1=y2) and ((x1=x2) or (x3=x2)) Then

  Begin

     If x1=x2 then Begin x:=x1; y:=y3 End

     Else Begin x:=x3; y:=y1 End;

     Exit;

  End;

  If (y3=y4) and ((x1=x2) or (x3=x2)) Then

  Begin

     If x3=x4 then Begin x:=x3; y:=y1 End

     Else Begin x:=x1; y:=y3 End;

     Exit;

  End;

  k1:=(y1-y2)/(x1-x2);

  k2:=(y3-y4)/(x3-x4);

  If k1=k2 Then Begin code:=False; Exit End;

  b1:=y1-k1*x1;

  b2:=y3-k2*x3;

  If k1-k2=0 then

  Begin

     code:=False;

     Exit

  End;

  x:=(b2-b1)/(k1-k2);

  y:=k1*x+b1;

  code:=True

End;

 

function PixelInOtr(x1,x2,x:glfloat):boolean;

begin

if x1<x2 then PixelInOtr:=(x1<x)and(x<x2);

if x2<x1 then PixelInOtr:=(x2<x)and(x<x1);

end;

 

function get_dl_line(x1,y1,z1,x2,y2,z2:glfloat):glfloat;

var x,y,z:glfloat;

begin

x:=x2-x1;y:=y2-y1;z:=z2-z1;

get_dl_line:=sqrt((x*x)+(y*y)+(z*z));

end;

 

 

procedure rotate_point(angle:glfloat;var x,y:glfloat);

var a,r:real;

begin

if x=0 then

begin

  r:=abs(y);

  if y>=0 then a:=90;

  if y<0 then a:=270;

end else

begin

  r:=sqrt((x*x)+(y*y));

  a:=abs(arctan(y/x)*180/pi);

end;

 

if (x<0)and(y>0)then a:=180-a; {2- 绁猗ム忪}

if (x<0)and(y<0)then a:=180+a; {3- 绁猗ム忪}

if (x>0)and(y<0)then a:=360-a; {4- 绁猗ム忪}

 

a:=((a+angle)*pi/180);

x:=r*cos(a);

y:=r*sin(a);

end;

 

 

function g_ang(x,y:glfloat):real;

var a:glfloat;

begin

if x=0 then

begin

  if y>0 then a:=90;

  if y<0 then a:=270;

end else

begin

  a:=abs(arctan(y/x)*180/pi);

  if (x<0)and(y>0)then a:=180-a; {2- 绁猗ム忪}

  if (x<0)and(y<0)then a:=180+a; {3- 绁猗ム忪}

  if (x>0)and(y<0)then a:=360-a; {4- 绁猗ム忪}

  if (y=0)and(x<0)then a:=180;

end;

g_ang:=a;

end;

 

 

function pt_in_tr(x1,y1,x2,y2,x3,y3,x,y:glfloat):boolean;

var angs,b:array[0..2] of real;

   itog:real;

begin

if ((not((x>x1)and(x>x2)and(x>x3)))and(not((x<x1)and(x<x2)and(x<x3))))

and ((not((y>y1)and(y>y2)and(y>y3)))and(not((y<y1)and(y<y2)and(y<y3))))

then begin

   angs[0]:=g_ang(x1-x,y1-y);

   angs[1]:=g_ang(x2-x,y2-y);

   angs[2]:=g_ang(x3-x,y3-y);

 

   b[0]:=abs(angs[0]-angs[1]);

   if b[0]>180 then b[0]:=360-b[0];

   b[1]:=abs(angs[1]-angs[2]);

   if b[1]>180 then b[1]:=360-b[1];

   b[2]:=abs(angs[2]-angs[0]);

   if b[2]>180 then b[2]:=360-b[2];

 

   itog:=b[0]+b[1]+b[2];

   if itog>=360 then pt_in_tr:=true

   else pt_in_tr:=false;

end else pt_in_tr:=false;

end;

 

constructor TGeo3DMesh.Create(BB, TB : string);

begin

       nXWidth := 257; nYWidth := 257;

       MeshConfig.dBaseHeight := - 1.5;

 

   MinH := 0;

   MaxH := 0;

   ShowMinRange := 0;

   ShowMaxRange := 1;

   bShowRange := True;

 

   Base := TBitmap.Create;

   Base.LoadFromFile(BB);

 

       CreateMeshFromBitmap(Base);

 

//    TexBMP := TBitmap.Create;

//    TexBmp.LoadFromFile(TB);

 

   Texture := TTextureGl.Create;

   //by file mapping

       Texture.LoadFrom_bmp_File3(TB);

 

//    LoadTexture(TexBMP);

end;

 

destructor TGeo3DMesh.Destroy;

begin

       if Assigned(Pts) then FreeMem(Pts, SizeOf(TGeo3DPt) * nXWidth * nYWidth);

   if Assigned(Texture) then Texture.Free;

   if Assigned(Base) then Base.Free;

//    if Assigned(TexBMP) then TexBMP.Free;

       if Assigned(Root) then begin

          ClearRootPts;

          Root.Free;

   end;

       inherited;

end;

 

procedure TGeo3DMesh.LoadTexture(BMP : TBitmap);

begin

       Texture.LoadFrom_bmp_File2(BMP);

end;

 

procedure TGeo3DMesh.CalcRect;

var

       i:integer;

begin

   with MeshRect do begin

          dLeft := MaxDouble;

       dTop := MaxDouble;

       dRight := -MaxDouble;

       dBottom := -MaxDouble;

   end;

 

         for i := 0 to nXWidth * nYWidth - 1 do begin

       if TGeo3DPt(Pts[I]).dX < MeshRect.dLeft then MeshRect.dLeft := TGeo3DPt(Pts[I]).dX;

       if TGeo3DPt(Pts[I]).dX > MeshRect.dRight then MeshRect.dRight := TGeo3DPt(Pts[I]).dX;

       if TGeo3DPt(Pts[I]).dZ < MeshRect.dTop then MeshRect.dTop := TGeo3DPt(Pts[I]).dZ;

       if TGeo3DPt(Pts[I]).dZ > MeshRect.dBottom then MeshRect.dBottom := TGeo3DPt(Pts[I]).dZ;

       end;

end;

 

procedure TGeo3DMesh.CalRange;

var

       i : Integer;

   dy : glfloat;

begin

       MinH := vrts[fcss[0].vrt[0]].crd[1];

   MaxH := vrts[fcss[0].vrt[0]].crd[1];

       for i:=1 to cFcs-1 do begin

       dy := vrts[fcss[i].vrt[0]].crd[1];

       if MinH > dy then MinH := dy;

       if MaxH < dy then MaxH := dy;

       end;

end;

 

procedure TGeo3DMesh.GetDataFromVtrs;

var

       I : integer;

begin

       for I := 0 to cFcs - 1 do begin

          if I >= nXWidth * nYWidth then Break;

          Pts[I].dX := vrts[fcss[i].vrt[0]].crd[0];

          Pts[I].dY := vrts[fcss[i].vrt[0]].crd[1];

          Pts[I].dZ := vrts[fcss[i].vrt[0]].crd[2];

   end;

end;

 

procedure TGeo3DMesh.CreateMeshFromBitmap(b: tbitmap);

var

       I : integer;

   fX,fY : Single;                //x方向和y方向的步长

   C : TColor;

begin

       //根据B位图的大小和当前nXWidth,nYWidth的长度来对bitmap进行像素点的采样

   GetMem(Pts, SizeOf(TGeo3DPt) * nXWidth * nYWidth);

   fX := b.Width / nXWidth;

   fY := b.Height / nYWidth;

   for I := 0 to nXWidth * nYWidth - 1 do begin

          C := B.Canvas.Pixels[Round(I mod nXWidth * fx), Round((I div nYWidth) * fy)];

               Pts[I].dX := (I mod nXWidth)/ nXWidth - 0.5;

          Pts[I].dY := (GetRValue(C) + GetGValue(C) + GetBValue(C))/(255 * 18) - 0.5;

       Pts[I].dZ := (I div nYWidth)/ nYWidth - 0.5;       

   end;

       CalcRect;   

end;

 

{ TGeoMeshTree }

 

procedure TGeoMeshNode.CalcSubRoot;

var

       bFound : Boolean;

       I, J : integer;

   Pt : PGeo3DPt;

   N : TGeoMeshNode;

   nPosX, nPosY : integer;                        //分界处

   dYAvg : double;

begin

       //先计算本结点中的Z值起伏情况,如果有变化,则再分割

   bFound := FALSE;

   //先算出所有Y的平均值

   dYAvg := 0.0;

   for I := 0 to Pts.Count - 1 do begin

               dYAvg := dYAvg + PGeo3dPt(Pts.Items[I]).dY;

   end;

   dYAvg := dYAvg / Pts.Count;

       for I := 0 to Pts.Count - 2 do begin

               Pt := PGeo3DPt(Pts.Items[I]);

       if Abs(Pt^.dY - dYAvg) > dOffset then begin

              bFound := TRUE;

           Break;

       end;

   end;

   if bFound then begin

          //划分四个单元

       nPosX := Trunc(nXWidth / 2) + 1;

       nPosY := Trunc(nYWidth / 2) + 1;

       //左上方格

               N := TGeoMeshNode.Create(Mesh);

       N.nXWidth := nPosX;

       N.nYWidth := nPosY;

 

       for I := 0 to nPosX - 1 do begin

              for J := 0 to nPosY - 1 do begin

                  N.Pts.Add(Pts.Items[I * nXWidth + J]);

           end;

       end;

 

       SubNodes.Add(N);

 

       //右上方格

               N := TGeoMeshNode.Create(Mesh);

       N.nXWidth := nPosX;

       N.nYWidth := nPosY;

 

       for I := 0 to nPosX - 1 do begin

              for J := 0 to nPosY - 1 do begin

                  N.Pts.Add(Pts.Items[I * nXWidth + J + nPosX - 1]);

           end;

       end;

 

       SubNodes.Add(N);

 

       //左下方格

               N := TGeoMeshNode.Create(Mesh);

       N.nXWidth := nPosX;

       N.nYWidth := nPosY;

 

       for I := nPosX - 1 to nXWidth - 1 do begin

              for J := 0 to nPosY - 1 do begin

                  N.Pts.Add(Pts.Items[I * nXWidth + J]);

           end;

       end;

 

       SubNodes.Add(N);

 

       //右下方格

               N := TGeoMeshNode.Create(Mesh);

       N.nXWidth := nPosX;

       N.nYWidth := nPosY;

 

       for I := nPosX - 1 to nXWidth - 1 do begin

                       for J := 0 to nPosY - 1 do begin

                  N.Pts.Add(Pts.Items[I * nXWidth + J + nPosX - 1]);

           end;

       end;

 

       SubNodes.Add(N);

 

       for I := 0 to SubNodes.Count - 1 do begin

              if tGeoMeshNode(SubNodes.Items[I]).Pts.Count > 4 then begin

                               TGeoMeshNode(SubNodes.Items[I]).CalcSubRoot;

           end;

       end;

   end

   else begin

          //此时Node中的Y方向变化没有发现,取四个角的数值

               for I := 0 to Pts.Count - 1 do begin

              if (I <> 0) and (I <> nXWidth - 1) and (I <> nXWidth * (nYWidth - 1)) and (I <> nXWidth * nYWidth - 1) then begin

                               //Pts中的Item不能被free,因为所有的Pts都是放在Root中分配和销毁的!

                               Pts.Items[I] := nil;

           end;

       end;

               Pts.Pack;

   end;

end;

 

constructor TGeoMeshNode.Create(_Mesh : TGeo3DMesh);

begin

       Mesh := _Mesh;

       dOffset := 0.0001;

   SubNodes := TList.Create;

   Pts := TList.Create;

end;

 

destructor TGeoMeshNode.Destroy;

var

       I : integer;

begin

{

       不要在此销毁Pts,要在Mesh中显式地调用Node.ClearRootPts

       for I := 0 to Pts.Count - 1 do begin

          FreeMem(PGeo3DPt(Pts.Items[I]),SizeOf(TGeo3DPT));

   end;

}

   Pts.Free;

 

         for I := 0 to SubNodes.Count - 1 do begin

             TGeoMeshNode(SubNodes.Items[I]).Free;

   end;

 

       inherited;

end;

 

procedure TGeo3DMesh.BuildTree;

var

       I : integer;

       Pt : PGeo3DPt;

begin

       Root := TGeoMeshNode.Create(Self);

       for I := 0 to nXWidth * nYWidth - 1 do begin

               GetMem(Pt, SizeOf(TGeo3DPt));

       Pt^ := Pts[I];

               Root.Pts.Add(Pt);

   end;

 

   Root.nXWidth := nXWidth;

   Root.nYWidth := nYWidth;

   Root.CalcSubRoot;

end;

 

procedure TGeoMeshNode.Draw(MODE: glenum);

var

       I : integer;

       n1, n2: TGLArrayf3;   //normals

begin

       if SubNodes.Count > 0 then begin

          for I := 0 to SubNodes.Count - 1 do begin

              TGeoMeshNode(SubNodes.Items[I]).Draw(Mode);

       end;

   end

   else begin

       glBegin(mode);

           n1 := GetNormal(PGeo3DPt(pts.Items[0]).dX, PGeo3DPt(pts.Items[0]).dY, PGeo3DPt(pts.Items[0]).dZ,PGeo3DPt(pts.Items[2]).dX, PGeo3DPt(pts.Items[2]).dY, PGeo3DPt(pts.Items[2]).dZ ,PGeo3DPt(pts.Items[1]).dX, PGeo3DPt(pts.Items[1]).dY, PGeo3DPt(pts.Items[1]).dZ);

//                   n1 := GetNormal(pts[I].dX, pts[I].dY, Pts[I].dZ, pts[I + 1].dX, pts[I + 1].dY, Pts[I+ 1].dZ ,pts[I + nXWidth].dX, pts[I + nXWidth].dY, Pts[I + nXWidth].dZ);

           glNormal3fv(@n1);

           glTexCoord2f(pGeo3DPt(Pts.Items[0]).dX / (Mesh.MeshRect.dRight - Mesh.MeshRect.dLeft) + 0.5, pGeo3DPt(Pts.Items[0]).dZ / (Mesh.MeshRect.dBottom - Mesh.MeshRect.dTop) + 0.5);

           glVertex3f(PGeo3DPt(pts.Items[0]).dX, PGeo3DPt(pts.Items[0]).dy, PGeo3DPt(pts.Items[0]).dZ);

           glTexCoord2f(pGeo3DPt(Pts.Items[2]).dX / (Mesh.MeshRect.dRight - Mesh.MeshRect.dLeft) + 0.5, pGeo3DPt(Pts.Items[2]).dZ / (Mesh.MeshRect.dBottom - Mesh.MeshRect.dTop) + 0.5);

           glVertex3f(PGeo3DPt(pts.Items[2]).dX, PGeo3DPt(pts.Items[2]).dy, PGeo3DPt(pts.Items[2]).dZ);

           glTexCoord2f(pGeo3DPt(Pts.Items[1]).dX / (Mesh.MeshRect.dRight - Mesh.MeshRect.dLeft) + 0.5, pGeo3DPt(Pts.Items[1]).dZ / (Mesh.MeshRect.dBottom - Mesh.MeshRect.dTop) + 0.5);

           glVertex3f(PGeo3DPt(pts.Items[1]).dX, PGeo3DPt(pts.Items[1]).dy, PGeo3DPt(pts.Items[1]).dZ);

       glEnd;

 

       glBegin(mode);

           n1 := GetNormal(PGeo3DPt(pts.Items[3]).dX, PGeo3DPt(pts.Items[3]).dY, PGeo3DPt(pts.Items[3]).dZ,PGeo3DPt(pts.Items[1]).dX, PGeo3DPt(pts.Items[1]).dY, PGeo3DPt(pts.Items[1]).dZ ,PGeo3DPt(pts.Items[2]).dX, PGeo3DPt(pts.Items[2]).dY, PGeo3DPt(pts.Items[2]).dZ);

           glNormal3fv(@n1);

           glTexCoord2f(pGeo3DPt(Pts.Items[1]).dX / (Mesh.MeshRect.dRight - Mesh.MeshRect.dLeft) + 0.5, pGeo3DPt(Pts.Items[1]).dZ / (Mesh.MeshRect.dBottom - Mesh.MeshRect.dTop) + 0.5);

           glVertex3f(PGeo3DPt(pts.Items[1]).dX, PGeo3DPt(pts.Items[1]).dy, PGeo3DPt(pts.Items[1]).dZ);

           glTexCoord2f(pGeo3DPt(Pts.Items[2]).dX / (Mesh.MeshRect.dRight - Mesh.MeshRect.dLeft) + 0.5, pGeo3DPt(Pts.Items[2]).dZ / (Mesh.MeshRect.dBottom - Mesh.MeshRect.dTop) + 0.5);

           glVertex3f(PGeo3DPt(pts.Items[2]).dX, PGeo3DPt(pts.Items[2]).dy, PGeo3DPt(pts.Items[2]).dZ);

           glTexCoord2f(pGeo3DPt(Pts.Items[3]).dX / (Mesh.MeshRect.dRight - Mesh.MeshRect.dLeft) + 0.5, pGeo3DPt(Pts.Items[3]).dZ / (Mesh.MeshRect.dBottom - Mesh.MeshRect.dTop) + 0.5);

           glVertex3f(PGeo3DPt(pts.Items[3]).dX, PGeo3DPt(pts.Items[3]).dy, PGeo3DPt(pts.Items[3]).dZ);

       glEnd;

   end;

end;

 

procedure TGeo3DMesh.ClearRootPts;

var

       I : integer;

begin

       if Assigned(Root) then begin

               for I := 0 to Root.Pts.Count - 1 do begin

              FreeMem(Root.Pts.Items[I], SizeOf(TGeo3DPt));

       end;

       Root.Pts.Clear;

   end;

end;

 

end.