Програмирование изображения поверхности z=f(x,y) с удалением невидимых линий используя алгоритм плавающего горизонта

Точки (x,y,z) вычерчиваемой части поверхности принадлежат некоторому прямоугольному параллелепипеду x min  <= x <= x max , y min  <= y <= y max , z min  <= z <= z max . Точка наблюдения (x V  ,y V  ,z V ) удовлетворяет условиям x V   ?  x max , y V  ?  y max , z V  ?  z max . Функции ex и ey переводят мировые координаты (x,y,z) в координаты проекции (x ? ,y ? ) . Плоскость проекции имеет единичный вектор нормали u = (cos a , sin a , 0) и находится на расстоянии d от точки (x V , y V  ,z V ) (см. рассмотренный ранее пример 3.3 на центральную проекцию). Координаты центральной проекции вычисляются по формулам:

x '=- d (( x - x V )(- sin a )+( y - y V ) cos a )/(( x - x V ) cos a +( y - y V ) sin a ) ;

y'=-d(z-z V )/((x-x V )cos a +(y-y V )sin a ) .

Вычислим наибольшие и наименьшие координаты отображаемого параллелепипеда

ex min  = ex(x max  ,y min  ,z min ), ex max  = ex(x min , y max , z max ),

ey min  = ey(x max , y max , z min ), ey max  = ey(x min , y min , z max ).

Прямоугольник ex min  <= ex <= ex max , ey min  <= ey <= ey max будет окном, отображаемым на область вывода, состоящую из точек, имеющих целые координаты

X = (ex-exmin)*getmaxx()/(exmax-exmin),Y = (ey-eymin)*getmaxy()/(eymax-eymin).

Рассматривается прямоугольная сетка в плоскости Oxy с шагом hx = (xmax-xmin)/nx по x , и hy = (ymax -ymin )/ny - по y .

Каждому прямоугольнику (( x j , y i ), ( x j , y i +1 ), ( x j +1 , y i +1 ), ( x j +1 , y i )), 0 <= j <= nx -1, 0 <= j <= n - y -1 на плоскости Oxy соответствует четырехугольник ( f ( x j , y i ), f ( x j , y i +1 ), f ( x j +1 , y i +1 ), f ( x j +1 , y i )) в пространстве Oxyz . Будем приближать поверхность z = f ( x , y ) многогранником, состоящим из таких четырехугольников.

Элементом изображения является отрезок. Отрезки выводятся от ближнего к дальним. Для вывода видимой части отрезка используется алгоритм плавающего горизонта (функция sline и vectsec ).

Код на C++
#include <graphics.h>
#include <stdio.h>
#include <math.h>
#include <conio.h>
#define nx 60
#define ny 60
 
float xmax = 3, xmin = -3, ymax = 3, ymin = -3,  zmax = 3, zmin = -3;
float xv = 5, yv = 3.5 , zv = 4 ;	// положение наблюдателя
float d = 10 ;				// расстояние до плоскости проекции
float cosa, sina ;			// меридиана точки наблюдения
float exmax, exmin, eymax, eymin ;
int gmex, gmey ;
int fimin[640], fimax[640] ;
float cvals[ny+1], cval ;
 
// Генерация точек отрезка по алгоритму простого ЦДА
void sline( int X0, int Y0, int X1, int Y1 )
{
int i, n;
float x, y, dx=0, dy=0;
n = abs(X1-X0); if (n< abs(Y1-Y0) ) n = abs(Y1-Y0);
if (n>0) {dx=(0.+X1-X0)/n; dy=(0.+Y1-Y0)/n;}
x=X0; y=Y0;
if (dx>0) x+=0.5; else if (dx<0) x-=0.5;
if (dy>0) y+=0.5; else if (dy<0) y-=0.5;
for(i=0;i<=n;i++)
{
	if (y>fimax[x])	fimax[x]=y;	// Изменяем верхний горизонт
	if (y<fimin[x])	fimin[x]=y;	// Изменяем нижний горизонт
	x+=dx; y+=dy;
}
}
 
// Вывод видимой части отрезка
//с удалением невидимых линий используя алгоритм плавающего горизонта
void vectsec( int X0, int Y0, int X1, int Y1 )
{
int begx, begy, endx, endy, i;
sline(X0, Y0, X1, Y1);
begx=X0; begy=Y0;
endx=X1; endy=Y1;
if (begx>X1) {begx=X1; begy=Y1;	endx=X0; endy=Y0;}
moveto( begx, gmey-fimax[begx]);
for (i=begx;i<endx;i++) lineto(i, gmey-fimax[i]);
moveto(begx, gmey-fimin[begx]);
for (i=begx;i<endx;i++) lineto(i, gmey-fimin[i]);
}
 
// Функция z=f(x,y)
float fz(float x, float y) {return (x*sin(y)-y*sin(x));}
 
// Центральная проекция
// x координата на плоскости проекциии
float ex(float x, float y, float z)
	{return (-d*(-(x-xv)*sina+(y-yv)*cosa))/((x-xv)*cosa+(y-yv)*sina);}
// y координата на плоскости проекции
float ey ( float x, float y, float z )
	{return (-d*(z-zv))/((x-xv)*cosa+(y-yv)*sina);}
// с удалением невидимых линий
// Вычисление экранных координат
void vectfi(float x0, float y0, float z0, float x1, float y1, float z1,int mm)
{
float ex0, ex1, ey0, ey1;
int ix0, ix1, iy0, iy1;
ex0=ex(x0,y0,z0);
ix0=(ex0-exmin)*gmex/(exmax-exmin);
ey0=ey(x0,y0,z0);
iy0=(ey0-eymin)*gmey/(eymax-eymin);
ex1=ex(x1,y1,z1);
ix1=(ex1-exmin)*gmex/(exmax-exmin);
ey1=ey(x1,y1,z1);
iy1=(ey1-eymin)*gmey/(eymax-eymin);
 
// Если mm не равно нулю, то выводим
// изображение с удалением невидимых линий,
// иначе - без удаления невидимых линий
//с удалением невидимых линий используя алгоритм плавающего горизонта
if (mm) vectsec(ix0,iy0,ix1,iy1); else line(ix0,gmey-iy0,ix1,gmey-iy1);
}
 
main()
{
int graphdriver=DETECT, graphmode;
int i, j;
float x,y,hx,hy;
 
// Инициализация графического режима
initgraph(&graphdriver, &graphmode, "");
hy=(ymax-ymin)/ny;
hx=(xmax-xmin)/nx;
gmex=getmaxx(); gmey=getmaxy();
cosa=xv/sqrt(xv*xv+yv*yv);
sina=yv/sqrt(xv*xv+yv*yv );
for(i=0;i<gmex;i++) {fimax[i]=0; fimin[i]=gmey;}
exmax=ex(xmin, ymax, zmax);
exmin=ex(xmax, ymin, zmin);
eymax=ey(xmin, ymin, zmax);
eymin=ey(xmax, ymax, zmin);
 
// Заполнение первой строки
cval=fz(xmax, ymax);
for(i=ny-1;i>=0;i--)
	{
	cvals[i+1]=cval; y=ymin+hy*i;
	cval=fz(xmax, y);
	vectfi(xmax, y+hy, cvals[i+1], xmax, y, cval, 1);
	}
cvals[0] = cval ;
 
// Заполнение остальных строк
for (j=nx-1;j>=0;j--)
	{
	x=xmin+hx*j;
	cval=fz(x, ymax);
	vectfi(x+hx, ymax, cvals[ny], x,ymax,cval, 1);
	for(i=ny-1;i>=0;i--)
		{
		cvals[i+1]=cval; y=ymin+hy*i;
		cval=fz(x, y);
		vectfi(x, y+hy, cvals[i+1], x, y, cval, 1);
		vectfi(x+hx, y, cvals[i], x, y, cval, 1);
		}
	cvals[0] = cval ;
	}
 
// Изображение объемлющего параллелепипеда
vectfi(xmax, ymax, zmax, xmin, ymax, zmax, 0);
vectfi(xmax, ymax, zmax, xmax, ymax, zmin, 0);
vectfi(xmax, ymax, zmax, xmax, ymin, zmax, 0);
vectfi(xmin, ymax, zmin, xmin, ymax, zmax, 0);
vectfi(xmin, ymax, zmin, xmax, ymax, zmin, 0);
vectfi(xmax, ymin, zmin, xmax, ymin, zmax, 0);
vectfi(xmax, ymin, zmin, xmax, ymax, zmin, 0);
vectfi(xmin, ymin, zmin, xmin, ymin, zmax, 1);
vectfi(xmin, ymin, zmin, xmax, ymin, zmin, 1);
vectfi(xmin, ymin, zmin, xmin, ymax, zmin, 1);
vectfi(xmax, ymin, zmax, xmin, ymin, zmax, 1);
vectfi(xmin, ymax, zmax, xmin, ymin, zmax, 1);
getch(); closegraph();
}
При использовании обязательна ссылка на http://DMTSoft.ru