Алгоритм удаления невидимых линий и алгоритм закраски Гуро

#include <cstdlib>
#include <iostream>
#include <math.h>
#include <conio.h>
#include <graphics.h>
 
using namespace std;
 
int const n=40;   // Количество линий сетки по фи и пси.
int const m=20;
 
double const a=5;    // Параметры эллипсоида
double const b=6;
double const c=3;
 
double phiz=0*M_PI/180,  phiy=30*M_PI/180, phix=30*M_PI/180;  //  Поворот Объктной системы координат относительно мировой (в радианах)
double S[4]={0,0,1,0};//Вектор направленный на наблюдателя
double V[4]={-1,-1,-1,0};//Вектор направления падающего луча 
double K=0.8;// Коэфициент отражения
double Jf=0.3;// Интенсивность фонoвого освещения
double J=0.7; // Интенсивность источника
//======================================================
const int xemin=0, xemax=800, yemin=0, yemax=800;
const double xmin=-10, xmax=10, ymin=-10, ymax=10;
double KZ, xc, yc, xec, yec;
//======================================================
double Co_m[4][4]; //Матрица перехода из объектных координат в мировые
double Rz[4][4]={{cos(phiz),-sin(phiz),0,0},{sin(phiz),cos(phiz),0,0},{0,0,1,0},{0,0,0,1}};// Матрица поворота по оси z
double Ry[4][4]={{cos(phiy),0,sin(phiy),0},{0,1,0,0},{-sin(phiy),0,cos(phiy),0},{0,0,0,1}};// Матрица поворота по оси x
double Rx[4][4]={{1,0,0,0},{0,cos(phix),-sin(phix),0},{0,sin(phix),cos(phix),0},{0,0,0,1}};
//======================================================
void umn_matr(double a[4][4], double b[4][4], double c[4][4])// Процедура умножения матриц
{ 
    double s;
    for(int i=0; i<4; i++)
        for(int j=0; j<4; j++)
        {
            c[i][j]=0;
            s=0;
            for(int k=0; k<4; k++)
                s+=a[i][k]*b[k][j];
            c[i][j]=s;
        }                
}
//=========================================================
void str_na_matr(double d[4], double a[4][4], double c[4]){// Процедура умножения строки на матрицу
     double s;
      for(int i=0; i<4; i++){
                 s=0;
                 for(int k=0; k<4; k++)
                 s+=d[k]*a[k][i];
         c[i]=s;
         }
}
//==========================================================
double vision(double N[4])//коэффициент освещённости плоскости
{
    double dN=sqrt(N[0]*N[0]+N[0]*N[0]+N[0]*N[0]);
    double dV=sqrt(V[0]*V[0]+V[0]*V[0]+V[0]*V[0]);
    double dS=sqrt(S[0]*S[0]+S[0]*S[0]+S[0]*S[0]);
    double vis=(Jf+J*(V[0]*N[0]+V[1]*N[1]+V[2]*N[2])/dN/dV)*K*(S[0]*N[0]+S[1]*N[1]+N[2]*S[2])/dN/dS;
    if(vis>0)
        return vis;
    else
        return 0;
}
//======================================================
bool visible(double N[4],double S[4])
{
     return ((S[0]*N[0]+S[1]*N[1]+S[2]*N[2])>0);
}
//======================================================
double fx(double phi, double psi) // Задание уравнения эллипоида
{
       return (a*cos(phi))*cos(psi);
}
 
double fy(double phi, double psi)
{
       return (b*sin(phi))*cos(psi);
}
 
double fz(double psi)
{
       return c*sin(psi);
}
//====================================================
 
int main(int argc, char *argv[])
{
    double q00[4],q10[4],q01[4],q11[4], N[4], u[4], v[4];
    double qm00[4],qm10[4],qm01[4],qm11[4], Nm[4], Sm[4];
    double phi,psi; double KV;
    double x,y;
    q00[3]=1; q10[3]=1; q01[3]=1; q11[3]=1; N[3]=0; Sm[3]=0;
    umn_matr(Rx,Ry,Co_m);
    double kx=(xemax-xemin)/(xmax-xmin);
    double ky=(yemax-yemin)/(ymax-ymin);
    if (kx<ky) KZ=kx;
    else KZ=ky;
    xc=(xmin+xmax)/2;
    yc=(ymin+ymax)/2;
    xec=(xemin+xemax)/2;
    yec=(yemin+yemax)/2;
    
    initwindow(xemax-xemin, yemax-yemin);
    
    for(int i=0; i<n; i++)
    {
        phi=(i*2*M_PI)/n;
        for(int j=0; j<m; j++)
        {
            psi=(j*2*M_PI)/m;
            q00[0]=fx(phi,psi); q00[1]=fy(phi,psi); q00[2]=fz(psi);
                    
            phi=phi+(2*M_PI)/n;
            q10[0]=fx(phi,psi); q10[1]=fy(phi,psi); q10[2]=fz(psi);
                    
            psi=psi+(2*M_PI)/m;
            q11[0]=fx(phi,psi); q11[1]=fy(phi,psi); q11[2]=fz(psi);
                    
            phi=phi-(2*M_PI)/n;
            q01[0]=fx(phi,psi); q01[1]=fy(phi,psi); q01[2]=fz(psi);
                    
            for(int k=0; k<4; k++)
            {
                u[k]=q10[k]-q00[k];
                v[k]=q01[k]-q00[k];
            }
            N[0]=u[1]*v[2]-v[1]*u[2];
            N[1]=-u[0]*v[2]+v[0]*v[2];
            N[2]=u[0]*v[1]-v[0]*u[1];
            N[3]=0;
            
            str_na_matr(q00,Co_m,qm00);
            str_na_matr(q10,Co_m,qm10);
            str_na_matr(q01,Co_m,qm01);
            str_na_matr(q11,Co_m,qm11);
            str_na_matr(N,Co_m,Nm);
                    //cout<<"S={"<<Sm[0]<<";"<<Sm[1]<<";"<<Sm[2]<<"}\n";
                    //cout<<"N={"<<Nm[0]<<";"<<Nm[1]<<";"<<Nm[2]<<"}\n";
            if(visible(Nm,S))
            { 
                qm00[0]=qm00[0]*KZ+xec; qm00[1]=qm00[1]*KZ+yec;
                qm01[0]=qm01[0]*KZ+xec; qm01[1]=qm01[1]*KZ+yec;
                qm10[0]=qm10[0]*KZ+xec; qm10[1]=qm10[1]*KZ+yec;
                qm11[0]=qm11[0]*KZ+xec; qm11[1]=qm11[1]*KZ+yec;
                
                KV=int(vision(Nm)*255);
                
                x=(qm00[0]+qm01[0]+qm10[0]+qm11[0])/4;
                y=(qm00[1]+qm01[1]+qm10[1]+qm11[1])/4;
                
                setcolor(COLOR(255,255,255));
                
                setfillstyle(1,COLOR(KV,KV,KV));
                floodfill(x,y,COLOR(255,255,255));
 
                line(qm00[0],qm00[1],qm01[0],qm01[1]);
                line(qm01[0],qm01[1],qm11[0],qm11[1]);
                line(qm11[0],qm11[1],qm10[0],qm10[1]);
                line(qm10[0],qm10[1],qm00[0],qm00[1]);
                
                
            }
        }
    }
    system("PAUSE");
    return EXIT_SUCCESS;
}