bb.donnay-software.com

Code: Select all

Public Function TriangulationDelaunay(point As Variant) As Variant
Attribute TriangulationDelaunay.VB_Description = "Построение системы треугольников по набору точек point=mmatrix(XY)."
Attribute TriangulationDelaunay.VB_ProcData.VB_Invoke_Func = "\n28"
Dim n As Long
Dim nt As Long
Dim h As Long
Dim i As Long
Dim j As Long
Dim k As Long
Dim l As Long
Dim m As Long
Dim hn As Long
Dim hf As Boolean
Dim kt As Long
Dim kr As Double
Dim counta As Long
Dim counttr As Long
Dim minX As Double
Dim maxK As Double
Dim minlen As Double
Dim x1 As Double
Dim x2 As Double
Dim x3 As Double
Dim y1 As Double
Dim y2 As Double
Dim y3 As Double
Dim ax As Double
Dim ay As Double
Dim bx As Double
Dim by As Double
Dim xa As Double
Dim ya As Double
Dim sa As Double
Dim sb As Double
Dim sc As Double
Dim x As Double
Dim y As Double
Dim x0 As Double
Dim y0 As Double
Dim k1 As Double
Dim k2 As Double
Dim xx As Double
Dim yy As Double
Dim len1 As Double
Dim t1 As Double
Dim t2 As Double
Dim t3 As Double
Dim xc As Double
Dim yc As Double
Dim r2c As Double
Dim r2 As Double
Dim alive() As Long
Dim tri() As Long
Dim res() As Long
  n = UBound(point)
  nt = n * 10
  counta = 0
ReDim alive(1 To nt, 1 To 4)
ReDim tri(1 To nt, 1 To 3)
  i = 1
  minX = point(1, 1)
  For h = 2 To n
    If point(h, 1) < minX Then
      minX = point(h, 1)
      i = h
    End If
  Next h
  alive(1, 1) = i
  maxK = 0
  For h = 1 To n
    If h <> i Then
      If point(h, 1) = point(i, 1) Then
          j = h
          h = n
      Else
        kr = Abs((point(h, 2) - point(i, 2)) / (point(h, 1) - point(i, 1)))
        If kr > maxK Then
          maxK = kr
          j = h
        End If
      End If
    End If
  Next h
  alive(1, 2) = j
  alive(1, 3) = -1
  counta = 1
  counttr = 0
  h = 1
  kt = counta
  Do While (counta > 0 And kt < nt - 2)
    m = 0
    hf = False
    hn = 0
    For h = 1 To kt
      If (alive(h, 3) <> 0) And alive(h, 4) = 0 Then
        m = h
        h = kt
      End If
    Next h
    If m > 0 Then
      counta = counta - 1
      i = alive(m, 1)
      j = alive(m, 2)
      k = alive(m, 3)
      x1 = point(i, 1)
      y1 = point(i, 2)
      x2 = point(j, 1)
      y2 = point(j, 2)
      hn = 0
      For h = 1 To n
        hf = False
        If (h <> i) And (h <> j) And (h <> k) Then
          x3 = point(h, 1)
          y3 = point(h, 2)
          sc = x1 * (y2 - y3) + x2 * (y3 - y1) + x3 * (y1 - y2)
          If sc <> 0 Then
            t1 = x1 ^ 2 + y1 ^ 2
            t2 = x2 ^ 2 + y2 ^ 2
            t3 = x3 ^ 2 + y3 ^ 2
            sa = t1 * (y2 - y3) + t2 * (y3 - y1) + t3 * (y1 - y2)
            sb = t1 * (x2 - x3) + t2 * (x3 - x1) + t3 * (x1 - x2)
            xc = 0.5 * sa / sc
            yc = -0.5 * sb / sc
            r2c = (x1 - xc) ^ 2 + (y1 - yc) ^ 2
            For l = 1 To n
              If (l <> i) And (l <> j) And (l <> h) Then
                hf = True
                x = point(l, 1)
                y = point(l, 2)
                r2 = (x - xc) ^ 2 + (y - yc) ^ 2
                If r2 < r2c Then
                  hf = False
                  hn = 0
                  l = n
                Else
                  hf = True
                End If
              End If
            Next l
          End If
        End If
        If hf Then
          hn = h
          h = n
        End If
      Next h
      If hf Then
        alive(m, 4) = hn
        k = 0
        For h = 1 To kt
          If (alive(h, 1) = i And alive(h, 2) = hn) Or (alive(h, 1) = hn And alive(h, 2) = i) Then
            If alive(h, 3) <> 0 Then k = h
            h = kt
          End If
        Next h
        If k = 0 Then
          kt = kt + 1
          alive(kt, 1) = i
          alive(kt, 2) = hn
          alive(kt, 3) = j
          counta = counta + 1
        ElseIf k > 0 Then
          alive(k, 4) = j
          counta = counta - 1
        End If
        k = 0
        For h = 1 To kt
          If (alive(h, 1) = j And alive(h, 2) = hn) Or (alive(h, 1) = hn And alive(h, 2) = j) Then
            If alive(h, 3) <> 0 Then k = h
            h = kt
          End If
        Next h
        If k = 0 Then
          kt = kt + 1
          alive(kt, 1) = j
          alive(kt, 2) = hn
          alive(kt, 3) = i
          counta = counta + 1
        ElseIf k > 0 Then
          alive(k, 4) = i
          counta = counta - 1
        End If
        counttr = counttr + 1
        tri(counttr, 1) = i
        tri(counttr, 2) = j
        tri(counttr, 3) = hn
      Else
        alive(m, 4) = -1
      End If
    End If
  Loop
ReDim res(1 To counttr, 1 To 3)
  For h = 1 To counttr
    For m = 1 To 3
      res(h, m) = tri(h, m)
    Next m
  Next h
  TriangulationDelaunay = res
End Function

Code: Select all

unit DeloneTriangulation;

interface

type
TTPoint=record
   x,y:integer;
   end;
TRib=record
   p1,p2:integer;
   end;
TTriangle=record
   p1,p2,p3:integer;
   end;
TCircle=record
   xo,yo,R:single;
   end;

var
  Points:array [0..16383] of TTPoint;
  PointsCount:integer;
  Ribs:array [0..16383] of TRib;
  RibsCount:integer;
  Triangles:array [0..16383] of TTriangle;
  TrianglesCount:integer;

procedure Triangulation;  

implementation

function Side(i,j,k:integer):integer;
var x1,y1,x2,y2,xo,yo,dx,dy,a,b,v:single;
begin
x1:=Points[i].x;
y1:=Points[i].y;
x2:=Points[j].x;
y2:=Points[j].y;
xo:=Points[k].x;
yo:=Points[k].y;
dx:=x2-x1;
dy:=y2-y1;
if abs(dx)>abs(dy) then
	 begin
   a:=dy/dx;
	 b:=y1-a*x1;
   v:=a*xo+b;
   if yo>v then result:=0 else result:=1;
	 end
	 else
	 begin
	 a:=dx/dy;
	 b:=x1-a*y1;
   v:=a*yo+b;
   if xo>v then result:=0 else result:=1;
	 end;
end;

function TriangleExists(p1,p2,p3:integer):boolean;
var i:integer;
begin
result:=true;
for i:=TrianglesCount-1 downto 0 do if ((p1=triangles[i].p1) or (p1=triangles[i].p2) or (p1=triangles[i].p3)) and
                                       ((p2=triangles[i].p1) or (p2=triangles[i].p2) or (p2=triangles[i].p3)) and
                                       ((p3=triangles[i].p1) or (p3=triangles[i].p2) or (p3=triangles[i].p3)) then exit;
result:=false;   
end;

function  SolveCircle(x1,y1,x2,y2,x3,y3:single):TCircle;
var ma,mb,dx1,dy1,dx2,dy2,dm:single;
begin
dx1:=x2-x1;dy1:=y2-y1;
dx2:=x3-x2;dy2:=y3-y2;
if abs(dx1)<0.01 then begin x1:=x1-0.1;dx1:=x2-x1;end;
if abs(dx2)<0.01 then begin x3:=x3+0.1;dx2:=x3-x2;end;
if abs(dy1)<0.01 then begin y1:=y1-0.1;dy1:=y2-y1;end;
if abs(dy2)<0.01 then begin y3:=y3+0.1;dy2:=y3-y2;end;
ma:=dy1/dx1;
mb:=dy2/dx2;
dm:=mb-ma;
if abs(dm)<0.0000001 then begin y3:=y3+0.1;dy2:=y3-y2;mb:=dy2/dx2;dm:=mb-ma;end;
result.xo:=(ma*mb*(y1-y3)+mb*(x1+x2)-ma*(x2+x3))*0.5/dm;
result.yo:=-1/mb*(result.xo-(x2+x3)*0.5)+(y2+y3)*0.5;
dx1:=x1-result.xo;
dy1:=y1-result.yo;
result.R:=sqrt(dx1*dx1+dy1*dy1);
end;

function FindPoint(r1,r2:integer):integer;
var i,j:integer;
cr:TCircle;
b:boolean;
x2,y2,v:single;
begin
result:=-1;
for i:=0 to pointsCount-1 do
 if (i<>r1) and (i<>r2) and (not TriangleExists(r1,r2,i)) then
   begin
   cr:=SolveCircle(points[r1].x,points[r1].y,points[r2].x,points[r2].y,points[i].x,points[i].y);
   b:=true;
   for j:=0 to pointsCount-1 do if (j<>r1) and (j<>r2) and (j<>i) then
      begin
      x2:=points[j].x-cr.xo;
      y2:=points[j].y-cr.yo;
      v:=sqrt(x2*x2+y2*y2);
      if v<cr.R then begin b:=false;break;end;
      end;
   if b then
      begin
      result:=i;
      exit;
      end;
   end;
end;

procedure FindFirstRib;
var i,j,k,n:integer;
st_1,st_0:boolean;
begin
for i:=0 to pointsCount-2 do
	 begin
	 for j:=i+1 to pointsCount-1 do
		  begin
	   	st_1:=false;
	  	st_0:=false;
	  	for k:=0 to pointsCount-1 do
		  	 begin
		  	 if (k<>i) and (k<>j) then
			  	  begin
			  	  n:=Side(i,j,k);
			  	  if n=1 then st_1:=true;
			  	  if n=0 then st_0:=true;
			  	  end;
			   end;
	  	if not (st_1 and st_0) then
			   begin
         Ribs[0].p1:=i;
		     Ribs[0].p2:=j;
		     RibsCount:=1;
         TrianglesCount:=0;
         exit;
			   end;
		  end;
	 end;
end;

procedure Triangulation;
var i,p1,p2,n:integer;
begin
FindFirstRib;
i:=0;
while (i<RibsCount) do
   begin
   p1:=Ribs[i].p1;
   p2:=Ribs[i].p2;
   n:=FindPoint(p1,p2);
   if n>=0 then
       begin
		   Ribs[RibsCount].p1:=p1;Ribs[RibsCount].p2:=n;
		   Ribs[RibsCount+1].p1:=p2;Ribs[RibsCount+1].p2:=n;
		   RibsCount:=RibsCount+2;
       triangles[TrianglesCount].p1:=p1;
       triangles[TrianglesCount].p2:=p2;
       triangles[TrianglesCount].p3:=n;
       TrianglesCount:=TrianglesCount+1;
		   end;
	 i:=i+1;
   end;
end;

end.

Code: Select all

unit DeloneTriangulation;

interface

type
TTPoint=record
   x,y:integer;
   end;
TRib=record
   p1,p2:integer;
   end;
TTriangle=record
   p1,p2,p3:integer;
   end;
TCircle=record
   xo,yo,R:single;
   end;

var
  Points:array [1..16384] of TTPoint;
  PointsCount:integer;
  Ribs:array [1..16384] of TRib;
  RibsCount:integer;
  Triangles:array [1..16384] of TTriangle;
  TrianglesCount:integer;

procedure Triangulation;  

implementation

function Side(i,j,k:integer):integer;
var x1,y1,x2,y2,xo,yo,dx,dy,a,b,v:single;
begin
x1:=Points[i].x;
y1:=Points[i].y;
x2:=Points[j].x;
y2:=Points[j].y;
xo:=Points[k].x;
yo:=Points[k].y;
dx:=x2-x1;
dy:=y2-y1;
if abs(dx)>abs(dy) then
	 begin
   a:=dy/dx;
	 b:=y1-a*x1;
   v:=a*xo+b;
   if yo>v then result:=0 else result:=1;
	 end
	 else
	 begin
	 a:=dx/dy;
	 b:=x1-a*y1;
   v:=a*yo+b;
   if xo>v then result:=0 else result:=1;
	 end;
end;

function TriangleExists(p1,p2,p3:integer):boolean;
var i:integer;
begin
result:=true;
for i:=TrianglesCount downto 1 do if ((p1=triangles[i].p1) or (p1=triangles[i].p2) or (p1=triangles[i].p3)) and
                                     ((p2=triangles[i].p1) or (p2=triangles[i].p2) or (p2=triangles[i].p3)) and
                                     ((p3=triangles[i].p1) or (p3=triangles[i].p2) or (p3=triangles[i].p3)) then exit;
result:=false;   
end;

function  SolveCircle(x1,y1,x2,y2,x3,y3:single):TCircle;
var ma,mb,dx1,dy1,dx2,dy2,dm:single;
begin
dx1:=x2-x1;dy1:=y2-y1;
dx2:=x3-x2;dy2:=y3-y2;
if abs(dx1)<0.01 then begin x1:=x1-0.1;dx1:=x2-x1;end;
if abs(dx2)<0.01 then begin x3:=x3+0.1;dx2:=x3-x2;end;
if abs(dy1)<0.01 then begin y1:=y1-0.1;dy1:=y2-y1;end;
if abs(dy2)<0.01 then begin y3:=y3+0.1;dy2:=y3-y2;end;
ma:=dy1/dx1;
mb:=dy2/dx2;
dm:=mb-ma;
if abs(dm)<0.0000001 then begin y3:=y3+0.1;dy2:=y3-y2;mb:=dy2/dx2;dm:=mb-ma;end;
result.xo:=(ma*mb*(y1-y3)+mb*(x1+x2)-ma*(x2+x3))*0.5/dm;
result.yo:=-1/mb*(result.xo-(x2+x3)*0.5)+(y2+y3)*0.5;
dx1:=x1-result.xo;
dy1:=y1-result.yo;
result.R:=sqrt(dx1*dx1+dy1*dy1);
end;

function FindPoint(r1,r2:integer):integer;
var i,j:integer;
cr:TCircle;
b:boolean;
x2,y2,v:single;
begin
result:=-1;
for i:=1 to pointsCount do
 if (i<>r1) and (i<>r2) and (not TriangleExists(r1,r2,i)) then
   begin
   cr:=SolveCircle(points[r1].x,points[r1].y,points[r2].x,points[r2].y,points[i].x,points[i].y);
   b:=true;
   for j:=1 to pointsCount do if (j<>r1) and (j<>r2) and (j<>i) then
      begin
      x2:=points[j].x-cr.xo;
      y2:=points[j].y-cr.yo;
      v:=sqrt(x2*x2+y2*y2);
      if v<cr.R then begin b:=false;break;end;
      end;
   if b then
      begin
      result:=i;
      exit;
      end;
   end;
end;

procedure FindFirstRib;
var i,j,k,n:integer;
st_1,st_0:boolean;
begin
for i:=1 to pointsCount-1 do
	 begin
	 for j:=i+1 to pointsCount do
		  begin
	   	st_1:=false;
	  	st_0:=false;
	  	for k:=1 to pointsCount do
		  	 begin
		  	 if (k<>i) and (k<>j) then
			  	  begin
			  	  n:=Side(i,j,k);
			  	  if n=1 then st_1:=true;
			  	  if n=0 then st_0:=true;
			  	  end;
			   end;
	  	if not (st_1 and st_0) then
			   begin
         Ribs[1].p1:=i;
		     Ribs[1].p2:=j;
		     RibsCount:=1;
         exit;
			   end;
		  end;
	 end;
end;

procedure Triangulation;
var i,p1,p2,n:integer;
begin
TrianglesCount:=0;
FindFirstRib;
i:=1;
while (i<RibsCount+1) do
   begin
   p1:=Ribs[i].p1;
   p2:=Ribs[i].p2;
   n:=FindPoint(p1,p2);
   if n>=0 then
       begin
		   Ribs[RibsCount].p1:=p1;Ribs[RibsCount].p2:=n;
		   Ribs[RibsCount+1].p1:=p2;Ribs[RibsCount+1].p2:=n;
		   RibsCount:=RibsCount+2;
       triangles[TrianglesCount].p1:=p1;
       triangles[TrianglesCount].p2:=p2;
       triangles[TrianglesCount].p3:=n;
       TrianglesCount:=TrianglesCount+1;
		   end;
	 i:=i+1;
   end;
end;

end.

Code: Select all

/*
*Обращаться очень просто - там есть массив Points в который записывать точки, переменная PointsCount в 
*которую записывать количество точек. Вам нужно заполнить точками этот массив перед вызовом функции и 
*присвоить соответствующее значение переменной PointsCount. Затем нужно вызвать одну единственную 
*функцию - Triangulation. В функцию не передаются никакие параметры. И после того как функция отработает - 
*получаем заполненный массив рёбер Ribs и заполненный массив треугольников Triangles - которые вы 
*можете прочитать после выполнения функции Triangulation. Количества рёбер и треугольников так же можно 
*прочитать в соответствующих переменных RibsCount и TrianglesCount. В массиве Triangles треугольники 
*описаны тремя целыми числами - это номера точек в массиве Points. То есть допустим если треугольник 
*описан числами 1,2,3 это значит что координаты точек нужно брать из ячеек points[1], points[2] и points[3] 
*соответственно. 

Ribs - это массив ребер. Он заполняется в процессе работы функции триангуляции, нужен для работы самой процедуры 
но может затем использовать и для любых других нужд. В этом массиве просто записан список ребер которые обнаруживаются 
в ходе триангуляции. Каждое ребро описано двумя целыми числами - это номера 2-х точек задающих ребро. Номера точек 
- это их индексы в массиве points.

Triangles - это массив треугольников. Он заполняется в процессе работы функции триангуляции, используется для работы 
самой функции а так же является результатом работы функции. В этом массиве просто записан список треугольников в порядке 
как они обнаруживались в ходе работы функции. Каждый треугольник описан тремя целыми числами. Каждое из этих чисел имеет 
тот же смысл что и в массиве Ribs - это номера точек задающих вершины треугольника. Номера точек - это их индексы в массиве points.

Кстати, из глобальных переменных у меня там только 6 - три массива и три целочисленных переменных. 
Это массив Triangles, массив Ribs, массив Points, переменная TrianglesCount, перменная RibsCount, переменная PointsCount. 
Всё - на этом список исчерпан. ВСЕ остальные переменные у меня там локальные.
*/

#include "appevent.ch"
#include "axcdxcmx.ch"
#include "collat.ch"
#include "common.ch"
#include "dbedit.ch"
#include "dbfdbe.ch"
#include "dcapp.ch"
#include "dcbitmap.ch"
#include "dccargo.ch"
#include "dccursor.ch"
#include "dcdialog.ch"
#include "dcdir.ch"
#include "dcfiles.ch"
#include "dcgra.ch"
#include "dcgraph.ch"        // графика
#include "BdColors.Ch"       // графика
#include "dccolors.ch"       // графика
#include "dcprint.ch"        // графика
*#INCLUDE "rmchart.CH"       // графика
#include "dcicon.ch"
#include "dcmsg.ch"
#include "dcpick.ch"
#include "deldbe.ch"
#include "directry.ch"
#include "dmlb.ch"
#include "express.ch"
#include "fileio.ch"
#include "font.ch"
#include "gra.ch"
#include "inkey.ch"
#include "memvar.ch"
#include "natmsg.ch"
#include "prompt.ch"
#include '_dcdbfil.ch'
*#INCLUDE "dcads.CH"
#include "set.ch"
#include "std.ch"
#include "xbp.ch"
#include '_dcappe.ch'
#include 'dcscope.ch'
#include '_dcstru.ch'
#include 'dcfields.ch'
#include 'dccolor.ch'
#INCLUDE "dll.CH"

#pragma library( "ascom10.lib" )
#pragma library( "dclip1.lib" )
#pragma library( "dclip2.lib" )
#pragma library( "dclipx.lib" )
#pragma library( "xbtbase1.lib" )
#pragma library( "xbtbase2.lib" )
#pragma library( "xppui2.lib" )

#INCLUDE "dll.CH"
#INCLUDE "dcdialog.CH"
#DEFINE SRCCOPY  0xCC0020

STATIC snHdll

PROCEDURE AppSys
// Рабочий стол остается окном приложения
RETURN

*****************************************************************************
FUNCTION Main()

PUBLIC GetList[0], GetOptions, oSay, hDC1, hDC2, oStatic1, oStatic2, aPixel

PUBLIC nColorB := AutomationTranslateColor(GraMakeRGBColor({0,0,0}),.f.)          // Черные пиксели
PUBLIC nColorG := AutomationTranslateColor(GraMakeRGBColor({200,200,200}),.f.)    // Серые  пиксели
PUBLIC nColorR := AutomationTranslateColor(GraMakeRGBColor({255,050,039}),.f.)    // Ярко-красный пиксель

DC_IconDefault(1000)

*********** Формирование массива точек

PUBLIC X := {}, Y := {}, Z := {}, mFlagCircle := .T.                              // Координаты X,Y,Z точек облак
PUBLIC TrianglesP1:= {}, TrianglesP2:= {}, TrianglesP3:= {}                       // Массивы номеров точек вершин треугольников
PUBLIC RibsP1:= {}, RibsP2:= {}, Points:= {}                                      // Массивы номеров точек ребер
PUBLIC TrianglesCount:=0, RibsCount:=0, PointsCount:=20                           // Кол-во треугольников, ребер, точек

PUBLIC X_MaxW := 1300, Y_MaxW := 700                                              // Размер графического окна для самого графика в пикселях

PUBLIC nXSize := X_MaxW
PUBLIC nYSize := Y_MaxW

** Максимальные значения x,y,z

maxX = nXSize-10
maxY = nYSize-10
maxZ = 1000

** Имя графического файла для рисования

mFileName = 'Delone.jpg'

H = 20  // Высота кнопки
W =  8  // Ширина кнопки
D =  5  // Расстояние между кнопками

@ 0,0 DCSTATIC TYPE XBPSTATIC_TYPE_BITMAP ;
      CAPTION mFileName ;
      OBJECT oStatic1 ;
      PREEVAL {|o|o:autoSize := .t.} ;
      EVAL {|o|hDC1 := GetWindowDC(o:getHWnd()), ;
               o:motion := {|a,b,o|ShowColor( hDC1, a, oSay, o )}, ;
               aPixel := Array(o:caption:xSize,o:caption:ySize)}

@ 0,0 DCSTATIC TYPE XBPSTATIC_TYPE_BITMAP;
      CAPTION mFileName ;
      PREEVAL {|o|o:autoSize := .t.} ;
      OBJECT oStatic2 ;
      EVAL {|o|hDC2 := GetWindowDC(o:getHWnd())}

@  40,2 DCSAY '' SAYSIZE 350,20 FONT '10.Lucida Console' OBJECT oSay
@  40,2 DCPUSHBUTTON                    CAPTION 'Очистка'                SIZE 100, H ACTION {||ClearImage(hDC2,aPixel)}
@ DCGUI_ROW, DCGUI_COL + D DCPUSHBUTTON CAPTION 'Генерация облака точек' SIZE 200, H ACTION {||GenPoints(hDC2,PointsCount,.T.)}
@ DCGUI_ROW, DCGUI_COL + D DCPUSHBUTTON CAPTION 'Триангуляция (сетка)'   SIZE 150, H ACTION {||Triangulation(hDC2)}
@ DCGUI_ROW, DCGUI_COL + D DCPUSHBUTTON CAPTION 'Триангуляция (цвет)'    SIZE 150, H ACTION {||Shading(hDC2)}
@ DCGUI_ROW, DCGUI_COL + D DCPUSHBUTTON CAPTION 'Поиск 1-го ребра'       SIZE 100, H ACTION {||FindFirstRib(hDC2)}
@ DCGUI_ROW, DCGUI_COL + 80*D DCCHECKBOX mFlagCircle PROMPT  'Рисовать окружности?'

DCGETOPTIONS PIXEL

DCREAD GUI FIT TITLE 'Триангуляция Делоне' OPTIONS GetOptions ;
   EVAL {||GenPoints(hDC2,PointsCount,.F.)} SETAPPWINDOW

CLOSE ALL

RETURN NIL
*****************************************************************************

*--------------------

FUNCTION FRND(mMax)
RETURN(1 + INT(RANDOM() / 65535 * mMax))

*--------------------

FUNCTION LB_Warning( message, ctitle )

  LOCAL aMsg := {}
  DEFAULT cTitle TO ''
  IF VALTYPE(message) # 'A'
    aadd(aMsg,message)
  ELSE
    aMsg := message
  ENDIF
  IF LEN(ALLTRIM(cTitle)) > 0
     DC_MsgBox( ,,aMsg,cTitle)
  ELSE
     DC_MsgBox( ,,aMsg,'Универсальная когнитивная аналитическая система "Эйдос-Х++"')
  ENDIF

RETURN NIL

*--------------------

******** Задание количества точек
FUNCTION NPoints(oStatic)

LOCAL GetList[0], GetOptions, oSay

@10,10 DCGROUP oGroup1 CAPTION 'Задайте количество точек:' SIZE 23.0, 2.5
@ 1, 1 DCSAY "" GET PointsCount PICTURE "##########" PARENT oGroup1

   DCGETOPTIONS TABSTOP
   DCREAD GUI ;
      TO lExit ;
      FIT ;
      OPTIONS GetOptions ;
      ADDBUTTONS;
      MODAL ;
      TITLE 'Триангуляция Делоне'

      ********************************************************************
      IF lExit
         ** Button Ok
      ELSE
         QUIT
      ENDIF
      ********************************************************************

RETURN(PointsCount)

* ---------

******** Генерация и отображение облака точек
FUNCTION GenPoints(hDC,PointsCount,mClear)

LOCAL GetList[0], GetOptions, oSay, oDevice

LOCAL hMemoryDC := hDC      // CreateMemoryDC( hDC, nXSize, nYSize )

PUBLIC X := {}, Y := {}, Z := {}                                                  // Координаты X,Y,Z точек облак
PUBLIC TrianglesP1:= {}, TrianglesP2:= {}, TrianglesP3:= {}                       // Массивы номеров точек вершин треугольников
PUBLIC RibsP1:= {}, RibsP2:= {}, Points:= {}                                      // Массивы номеров точек ребер
*PUBLIC TrianglesCount:=0, RibsCount:=0, PointsCount:=20                          // Кол-во треугольников, ребер, точек

PointsCount = NPoints()     // Задание количества точек

IF mClear
   ClearImage(hDC,aPixel)
ENDIF

**** Создать БД для облака точек X,Y,Z

aStructure := { { "Num", "N",  15, 0 }, ;
                { "pX" , "N",  15, 0 }, ;
                { "pY" , "N",  15, 0 }, ;
                { "pZ" , "N",  15, 0 }  }
DbCreate( 'Points_XYZ', aStructure )

CLOSE ALL
USE Points_XYZ EXCLUSIVE NEW

SELECT Points_XYZ

** Максимальные значения x,y,z

maxX = nXSize-10
maxY = nYSize-10
maxZ = 1000

FOR p=1 TO PointsCount
    
    mX = FRND(maxX)
    mY = FRND(maxY)
    mZ = FRND(maxZ)

    AADD(X, mX)
    AADD(Y, mY)
    AADD(Z, mZ)

    APPEND BLANK
    REPLACE Num WITH p
    REPLACE pX  WITH X[p]
    REPLACE pY  WITH Y[p]
    REPLACE pZ  WITH Z[p]

    Circle(hDC,mX,mY,1,nColorR)          // Маленький кружочек (r=1)
    Circle(hDC,mX,mY,2,nColorR)          // Маленький кружочек (r=2)
    Circle(hDC,mX,mY,3,nColorB)          // Маленький кружочек (r=3)
    Circle(hDC,mX,mY,4,nColorG)          // Маленький кружочек (r=4)

NEXT

LB_Warning( 'Построение точек завершено','Триангуляция Делоне' )

CLOSE ALL

RETURN nil

* ---------

******** Градиентная заливка
FUNCTION Shading(oStatic)


RETURN nil

* ---------

*************************************************

function Side(hDC,i,j,k)
LOCAL x1,y1,x2,y2,xo,yo,dx,dy,a,b,v

x1:=X[i]
y1:=Y[i]
x2:=X[j]
y2:=Y[j]
xo:=X[k]
yo:=Y[k]
dx:=x2-x1
dy:=y2-y1

if abs(dx)>abs(dy)
   a:=dy/dx
   b:=y1-a*x1
   v:=a*xo+b
   result = if(yo>v,0,1)
else
   a:=dx/dy
   b:=x1-a*y1
   v:=a*yo+b
   result = if(xo>v,0,1)
endif

*Circle(hDC,xo,yo,5,IF(result=1,nColorB,nColorR))          // Сделать отображение окружности, если это задано
*Line(hDC,x1,y1,x2,y2,nColorR)

RETURN(result)

*--------------------

function TriangleExists(p1,p2,p3)

LOCAL i

IF TrianglesCount = 0
   RETURN(.F.)
ELSE
   for i:=TrianglesCount to 1 STEP -1
       f1=.F.;if p1=trianglesP1[i] .or. p1=trianglesP2[i] .or. p1=trianglesP3[i];f1=.T.;endif
       f2=.F.;if p2=trianglesP1[i] .or. p2=trianglesP2[i] .or. p2=trianglesP3[i];f2=.T.;endif
       f3=.F.;if p3=trianglesP1[i] .or. p3=trianglesP2[i] .or. p3=trianglesP3[i];f3=.T.;endif
       IF f1 .and. f2 .and. f3
          RETURN(.T.)
       ENDIF
   NEXT
ENDIF
RETURN(.F.)

*--------------------

function SolveCircle(hDC,x1,y1,x2,y2,x3,y3)

LOCAL ma,mb,dx1,dy1,dx2,dy2,dm

*MsgBox("SolveCircle: x1,y1=("+str(x1)+","+str(y1)+"),  x2,y2=("+str(x2)+","+str(y2)+") , x3,y3=("+str(x3)+","+str(y3)+")")

dx1:=x2-x1;dy1:=y2-y1
dx2:=x3-x2;dy2:=y3-y2

if abs(dx1)<0.01;x1:=x1-0.1;dx1:=x2-x1;endif
if abs(dx2)<0.01;x3:=x3+0.1;dx2:=x3-x2;endif
if abs(dy1)<0.01;y1:=y1-0.1;dy1:=y2-y1;endif
if abs(dy2)<0.01;y3:=y3+0.1;dy2:=y3-y2;endif

ma:=dy1/dx1
mb:=dy2/dx2
dm:=mb-ma

if abs(dm)<0.0000001;y3:=y3+0.1;dy2:=y3-y2;mb:=dy2/dx2;dm:=mb-ma;endif

xo:=(ma*mb*(y1-y3)+mb*(x1+x2)-ma*(x2+x3))*0.5/dm
yo:=-1/mb*(xo-(x2+x3)*0.5)+(y2+y3)*0.5
dx1:=x1-xo
dy1:=y1-yo
R:=sqrt(dx1*dx1+dy1*dy1)

IF mFlagCircle
   Circle(hDC,xo,yo,R,nColorB)          // Сделать отображение окружности, если это задано
ENDIF

cr := {}
AADD(cr,xo)
AADD(cr,yo)
AADD(cr,R )

RETURN(cr)

*--------------------

function FindPoint(hDC,r1,r2)

LOCAL i,j,cr,b,x2,y2,v,xo,yo,R

for i:=1 to pointsCount

    if .not. TriangleExists(r1,r2,i)

       if i<>r1 .and. i<>r2

*         MsgBox("FindPoint: x1,y1=("+str(X[r1])+","+str(Y[r1])+"),  x2,y2=("+str(X[r2])+","+str(Y[r2])+") , x3,y3=("+str(X[i])+","+str(Y[i])+")")

          cr = SolveCircle(hDC,X[r1],Y[r1],X[r2],Y[r2],X[i],Y[i])

          xo=cr[1]
          yo=cr[2]
          R =cr[3]

          b:=.T.
          for j:=1 to pointsCount

              if j<>r1 .and. j<>r2 .and. j<>i
                 x2:=X[j]-xo
                 y2:=Y[j]-yo
                 v:=sqrt(x2*x2+y2*y2)

*                 Circle(hDC,xo,yo,R,nColorB)          // Сделать отображение окружности, если это задано
*                 Circle(hDC,X[j],Y[j],10,nColorR)     // Сделать отображение окружности, если это задано
*                 MsgBox('STOP')

                 if v<R
                    b:=.F.
                    EXIT
                 endif
              endif
          NEXT
          if b
             RETURN(i)
          endif
       endif
    endif
NEXT
RETURN(-1)

*--------------------

function FindFirstRib(hDC)

LOCAL i,j,k,n,st_1,st_0

for i:=1 to pointsCount-1
    for j:=i+1 to pointsCount
   	st_1:=.F.
  	st_0:=.F.
  	for k:=1 to pointsCount
   	    if k<>i .and. k<>j
  	       n:=Side(hDC,i,j,k)
	       if n=1;st_1:=.T.;endif
  	       if n=0;st_0:=.T.;endif
	    endif
        NEXT
	if st_1 <> st_0

           AADD(RibsP1, i)
           AADD(RibsP2, j)
           RibsCount:=1

           Line(hDC,X[i],Y[i],X[j],Y[j],nColorR)

           CLOSE ALL
           **** Создать БД для координат концов ребер
           aStructure := { { "Num" , "N",  15, 0 }, ;
                           { "pX1" , "N",  15, 0 }, ;
                           { "pY1" , "N",  15, 0 }, ;
                           { "pX2" , "N",  15, 0 }, ;
                           { "pY2" , "N",  15, 0 }, ;
                           { "pID" , "C",  20, 0 }  }
           DbCreate( 'Ribs_XY', aStructure )

           ar := {}
           AADD(ar, i)
           AADD(ar, j)
           ASORT(ar)
           mRibsID = STRTRAN(STR(ar[1])+STR(ar[2]),' ','_')

           CLOSE ALL
           USE Ribs_XY EXCLUSIVE NEW
           SELECT Ribs_XY

           APPEND BLANK
           REPLACE Num WITH i
           REPLACE pX1 WITH X[i]
           REPLACE pY1 WITH Y[i]
           REPLACE pX2 WITH X[j]
           REPLACE pY2 WITH Y[j]
           REPLACE pID WITH mRibsID

           CLOSE ALL
           RETURN NIL
	endif
*       MsgBox('STOP '+STR(j))
    NEXT
NEXT

RETURN NIL

*--------------------

FUNCTION Triangulation(hDC)

LOCAL i,p1,p2,n

FindFirstRib(hDC)

**** Создать БД для координат вершин треугольников

aStructure := { { "Num" , "N",  15, 0 }, ;
                { "pX1" , "N",  15, 0 }, ;
                { "pY1" , "N",  15, 0 }, ;
                { "pZ1" , "N",  15, 0 }, ;
                { "pX2" , "N",  15, 0 }, ;
                { "pY2" , "N",  15, 0 }, ;
                { "pZ2" , "N",  15, 0 }, ;
                { "pX3" , "N",  15, 0 }, ;
                { "pY3" , "N",  15, 0 }, ;
                { "pZ3" , "N",  15, 0 }, ;
                { "pID" , "C",  30, 0 }  }
DbCreate( 'Triang_XYZ', aStructure )

CLOSE ALL
USE Ribs_XY    EXCLUSIVE NEW
USE Triang_XYZ EXCLUSIVE NEW

TrianglesCount:=0

*MsgBox("Кол-во ребер: "+STR(RibsCount))

aRibsID   := {}
aTriangID := {}

TrianglesCount:=0

i:=1

DO WHILE i < RibsCount+1

   p1:=RibsP1[i]
   p2:=RibsP2[i]
   n:=FindPoint(hDC,p1,p2)          // Не происходит обход цикла и выход из цикла

*  MsgBox('p1='+STR(p1)+',  p2='+STR(p2)+',  n='+STR(n))
*  MsgBox("Номер найденной точки: "+STR(n)+", номер текущего ребра: "+STR(i))

   IF n > 0                         // Не происходит обход цикла и выход из цикла

      ********* Формирование ID ребер и тругольников и обход, если они уже есть

      SELECT Ribs_XY

      AADD(RibsP1, p1);AADD(RibsP2, n);RibsCount++

      ar := {}
      AADD(ar, p1)
      AADD(ar, n )
      ASORT(ar)
      mRibsID = STRTRAN(STR(ar[1])+STR(ar[2]),' ','_')

      APPEND BLANK
      REPLACE Num WITH i
      REPLACE pX1 WITH X[p1]
      REPLACE pY1 WITH Y[p1]
      REPLACE pX2 WITH X[n ]
      REPLACE pY2 WITH Y[n ]
      REPLACE pID WITH mRibsID

      Line(hDC,X[p1],Y[p1],X[n],Y[n],nColorR)

      ar := {}
      AADD(ar, p2)
      AADD(ar, n )
      ASORT(ar)
      mRibsID = STRTRAN(STR(ar[1])+STR(ar[2]),' ','_')

      AADD(RibsP1, p2);AADD(RibsP2, n);RibsCount++

      APPEND BLANK
      REPLACE Num WITH i
      REPLACE pX1 WITH X[p2]
      REPLACE pY1 WITH Y[p2]
      REPLACE pX2 WITH X[n ]
      REPLACE pY2 WITH Y[n ]
      REPLACE pID WITH mRibsID

      Line(hDC,X[p2],Y[p2],X[n],Y[n],nColorR)

      SELECT Triang_XYZ

      ar := {}
      AADD(ar, p1)
      AADD(ar, p2)
      AADD(ar, n )
      ASORT(ar)
      mTriangID = STRTRAN(STR(ar[1])+STR(ar[2])+STR(ar[3]),' ','_')

      AADD (aTriangID, mTriangID)

      AADD(trianglesP1, p1)
      AADD(trianglesP2, p2)
      AADD(trianglesP3, n )
      TrianglesCount++

      APPEND BLANK
      REPLACE Num WITH i
      REPLACE pX1 WITH X[p1]
      REPLACE pY1 WITH Y[p1]
      REPLACE pZ1 WITH Z[p1]
      REPLACE pX2 WITH X[p2]
      REPLACE pY2 WITH Y[p2]
      REPLACE pZ2 WITH Z[p2]
      REPLACE pX3 WITH X[n ]
      REPLACE pY3 WITH Y[n ]
      REPLACE pZ3 WITH Z[n ]
      REPLACE pID WITH mTriangID

      Line(hDC,X[p1],Y[p1],X[n ],Y[n ],nColorB)
      Line(hDC,X[p2],Y[p2],X[n ],Y[n ],nColorB)
      Line(hDC,X[p1],Y[p1],X[p2],Y[p2],nColorB)

   ENDIF

   i++

ENDDO

IF mFlagCircle

*  SELECT Ribs_XY
*  DBGOTOP()
*  DO WHILE .NOT. EOF()
*     Line(hDC,pX1,pY1,pX2,pY2,nColorR)
*     DBSKIP(1)
*  ENDDO

*   SELECT Triang_XYZ
*   DBGOTOP()
*   DO WHILE .NOT. EOF()
*      Line(hDC,pX1,pY1,pX2,pY2,nColorR)
*      Line(hDC,pX1,pY1,pX3,pY3,nColorR)
*      Line(hDC,pX2,pY2,pX3,pY3,nColorR)
*      DBSKIP(1)
*   ENDDO

    FOR j=1 TO LEN(trianglesP1)
        X1 = X[trianglesP1[j]]
        Y1 = Y[trianglesP1[j]]
        X2 = X[trianglesP2[j]]
        Y2 = Y[trianglesP2[j]]
        X3 = X[trianglesP3[j]]
        Y3 = Y[trianglesP3[j]]
        Line(hDC,X1,Y1,X2,Y2,nColorR)
        Line(hDC,X1,Y1,X3,Y3,nColorR)
        Line(hDC,X2,Y2,X3,Y3,nColorR)
    NEXT

ENDIF

CLOSE ALL

LB_Warning( 'Триангуляция завершена','Триангуляция Делоне' )

RETURN NIL

******************************************************************************************
******************************************************************************************
******************************************************************************************

* --------- Графика Роджера ---------------------------

FUNCTION LoadArray( hDC1, aPixel )

LOCAL hMemoryDC
LOCAL i, j, oScrn, nXSize := Len(aPixel), nYSize := Len(aPixel[1])
LOCAL nSeconds := Seconds()

/*
IF !aPixel[1,1] == nil
  DCMSGBOX 'Array is already loaded!'
  RETURN nil
ENDIF
*/

hMemoryDC := CreateMemoryDC( hDC1, nXSize, nYSize)

FOR i := 1 TO nXSize
  FOR j := 1 TO nYSize
    aPixel[i,j] := GetPixel(hMemoryDC,i-1,j-1)
  NEXT
NEXT

MsgBox(Alltrim(Str(Seconds()-nSeconds)) + ' Seconds to load Array')

DC_ClearEvents()

RETURN aPixel

* ---------

FUNCTION ClearImage( hDC2, aPixel )

LOCAL i, j, nXSize := Len(aPixel), nYSize := Len(aPixel[1])
LOCAL nColor := AutomationTranslateColor(GraMakeRGBColor({255,255,255}),.f.)
LOCAL hMemoryDC := hDC2 // CreateMemoryDC( hDC2, nXSize, nYSize )   // Для ускорения работы GetPixel() примерно в 50 раз

PUBLIC X := {}, Y := {}, Z := {}, mFlagCircle := .T.                              // Координаты X,Y,Z точек облак
PUBLIC TrianglesP1:= {}, TrianglesP2:= {}, TrianglesP3:= {}                       // Массивы номеров точек вершин треугольников
PUBLIC RibsP1:= {}, RibsP2:= {}, Points:= {}                                      // Массивы номеров точек ребер
PUBLIC TrianglesCount:=0, RibsCount:=0, PointsCount:=20                           // Кол-во треугольников, ребер, точек

FOR i := 0 TO nXSize
    FOR j := 0 TO nYSize
*       SetPixel(hMemoryDC, i, j, nColor)
        SetPixel(hDC2, i, j, nColor)
    NEXT
NEXT

LB_Warning( 'Очистка изображения завершена','Триангуляция Делоне' )

RETURN nil

* ----------

FUNCTION TransferImage( hDC1, hDC2, aPixel )

LOCAL i, j, nColor, lEmptyArray := aPixel[1,1] == nil, ;
      nXSize := Len(aPixel), nYSize := Len(aPixel[1])

LOCAL hMemoryDC := CreateMemoryDC( hDC1, nXSize, nYSize )

FOR i := 0 TO nXSize-1
  FOR j := 0 TO nYSize-1
    IF lEmptyArray
      SetPixel(hDC2,i,j,GetPixel(hMemoryDC,i,j))
    ELSE
      SetPixel(hDC2,i,j,aPixel[i+1,j+1])
    ENDIF
  NEXT
NEXT

RETURN nil

* ----------

FUNCTION FlipImage( hDC1, hDC2, aPixel )

LOCAL i, j, lEmptyArray := aPixel[1,1] == nil, ;
      nXSize := Len(aPixel), nYSize := Len(aPixel[1])

LOCAL hMemoryDC := CreateMemoryDC( hDC1, nXSize, nYSize )

FOR i := 0 TO nXSize-1
  FOR j := 0 TO nYSize-1
    IF lEmptyArray
      SetPixel(hDC2,j,i,GetPixel(hMemoryDC,j,nXSize-i))
    ELSE
      SetPixel(hDC2,j,i,aPixel[i+1,j+1])
    ENDIF
  NEXT
NEXT

RETURN nil

* -----------

FUNCTION RotateImage( hDC1, hDC2, aPixel )

LOCAL i, j, lEmptyArray := aPixel[1,1] == nil, ;
      nXSize := Len(aPixel), nYSize := Len(aPixel[1])

LOCAL hMemoryDC := CreateMemoryDC( hDC1, nXSize, nYSize )

FOR i := 0 TO nXSize-1
  FOR j := 0 TO nYSize-1
    IF lEmptyArray
      SetPixel(hDC2,i,j,GetPixel(hMemoryDC,j,nXSize-i))
    ELSE
      SetPixel(hDC2,i,j,aPixel[j+1,nXSize-i])
    ENDIF
  NEXT
NEXT

RETURN nil

* ---------

FUNCTION CreateMemoryDC( hDC, nXSize, nYSize )

LOCAL hMemoryDC, hBMP

hMemoryDC := CreateCompatibleDC(hDC)    // create compatible memory DC
hBMP      := CreateCompatibleBitmap(hDC,nXSize,nYSize) // create DDB
SelectObject(hMemoryDC,hBMP)                    // put hBMP into memory DC
BitBlt( hMemoryDC,0,0,nXSize,nYSize,hDC,0,0,SRCCOPY ) // copy desktop DC into memory DC

RETURN hMemoryDC

* ---------

STATIC FUNCTION ShowColor( hDC, aCoords, oSay, oStatic )

LOCAL nColor

aCoords[2] := oStatic:currentSize()[2] - aCoords[2]

nColor := GetPixel(hDC,aCoords[1],aCoords[2])

oSay:setCaption('Color: ' + DC_Array2String(GraGetRGBIntensity(AutomationTranslateColor(nColor,.T.))) + ;
   ' Coords: ' + DC_Array2String(aCoords))

RETURN nil

* ----------

#command  GDIFUNCTION <Func>([<x,...>]) ;
       => ;
FUNCTION <Func>([<x>]);;
STATIC scHCall := nil ;;
IF scHCall == nil ;;
  IF snHdll == nil ;;
    snHDll := DllLoad('GDI32.DLL') ;;
  ENDIF ;;
  scHCall := DllPrepareCall(snHDll,DLL_STDCALL,<(Func)>) ;;
ENDIF ;;
RETURN DllExecuteCall(scHCall,<x>)

GDIFUNCTION GetPixel( nHDC, x, y)
GDIFUNCTION SetPixel( nHDC, x, y, n )
DLLFUNCTION GetWindowDC( hwnd ) USING STDCALL FROM USER32.DLL
DLLFUNCTION CreateCompatibleDC( nHDC ) USING STDCALL FROM GDI32.DLL
DLLFUNCTION CreateCompatibleBitmap( nHDC, dw, dh ) USING STDCALL FROM GDI32.DLL
DLLFUNCTION SelectObject(hMemoryDC,hBMP) USING STDCALL FROM GDI32.DLL
DLLFUNCTION BitBlt( hDC,nXDest,nYDest,nXSize,nYSize,hDCSrc,nXSrc,nYSrc,dwROP ) USING STDCALL FROM GDI32.DLL

**********************************************************************************************

FUNCTION Circle(hDC,X0,Y0,R0,nColor)

FOR j = 1 TO 360 STEP 1
    nX = X0 + R0 * COS( j * 3.14159265358979323846 / 180 )
    nY = Y0 - R0 * SIN( j * 3.14159265358979323846 / 180 )
    SetPixel(hDC, nX, nY, nColor)
NEXT

RETURN nil

*--------------

FUNCTION Line(hDC,X1,Y1,X2,Y2,nColor)

FOR nX = X1 TO X2 STEP 0.01
    nY=Y1+(Y2-Y1)/(X2-X1)*(nX-X1)
    SetPixel(hDC, nX, nY, nColor)
NEXT
FOR nY = Y1 TO Y2 STEP 0.01
    nX=X1+(X2-X1)/(Y2-Y1)*(nY-Y1)
    SetPixel(hDC, nX, nY, nColor)
NEXT

RETURN nil