Красный матовый шар, освещённый солнцем... [ПАСКАЛЬ]

[ManHack]

Суть задачи такова:
нужно построить красный матовый шарик, освещённый вектором света.
Проблема в том, что светлое место на шарике должно переходить в полутень, а потом в тень плавно.
Но линии начала полутени и тени изогнуты по-разному (надо сделать так, чтобы они были изогнуты одинаково).
И, кроме того, при некоторых координатах, линия тени, идущая по экватору шара, становится не круглой (почти треугольной).
Это тоже надо исправить.

Вот уже написанный текст программы на Паскале:

Program graph2; {310 190}
uses Graph;
const
cir: record
x, y, r: longint;
kp, kh, kz :real;
end =(
x: 160;
y: 180;
r: 60;
kp: 1;
kh: 1;
kz: 0.5

);
sun: record
x, y, z: integer;
Ip, Ih: real;
end =(
x: -1; {70}
y: 1; {120}
z: 1; {80}
Ip: 8;
Ih: 56
);

{nabl: record
x, y, z :integer;
end =(
x: 20;
y: -40;
z: 20
);}

n = 200;
nah = 10;

long = 100;
dlin = 200;


var
xc, yc, zc, xs, ys, zs: integer;
{--------------}
function fX(x, y: real): integer;
begin
fX := round(x + y/(2*sqrt(2)) + nah);
end;



function fY(y, z: real): integer;
begin
fY := GetMaxY - nah - round(z + y/(2*sqrt(2)));
end;



function Color(I: real): integer;
var
j: integer;
begin
j := round(I);
if j>95 then
j := 95;
Color := j + 127;
end;




{---------------}

procedure SetPalette;
var
j : integer;
cR, cG, cB : integer;
begin
for j := 0 to 95 do begin
if j < 64 then begin
cR := j; cG := j - j; cB := j - j;
end;
{else if (j-64) < 32 then begin
cR := 255; cG := 255; cB := 255;
end;}

SetRGBPalette(j + 127, cR, cG, cB );
end;
end;

procedure Ris(x, y, z: real; var first: Boolean; col:integer);
begin
{
if first then begin
MoveTo(fX(x, y), fY(y, z));
first := false;
end
else
SetColor(col);
LineTo(fX(x, y), fY(y,z));
}
{if (cir.x*cir.x - x*x) + (cir.y*cir.y - y*y) + (cir.z*cir.z - z*z) <=
cir.r then }
putpixel(fX(x, y), fY(y, z),col);

end;

procedure RisNew(x, y:integer; col:integer);
begin
putpixel(x, y, col);

end;


procedure Display;
var
x,y:integer;
begin
{OX}
line(nah, GetMaxY - nah, nah + long, GetMaxY - nah);
line(nah + long, GetMaxY - nah, long, GetMaxY - 2*nah);
line(nah + long, GetMaxY - nah, long, GetMaxY);
OutTextXY(Long + nah, GetMaxY - nah, 'x');
{OY}
line(nah, GetMaxY - nah, fX(0, nah + long), fY(nah + long, 0));
line(fX(0, nah + long), fY(nah + long, 0), fX(0, nah + long) -
nah, fY(nah + long, 0));
line(fX(0, nah + long), fY(nah + long, 0), fX(0, nah + long),
fY(nah + long, 0) + nah);
OutTextXY(fX(0, nah + long), fY(nah + long, 0), 'y');
{OZ}
line(nah, GetMaxY - nah, nah, GetMaxY - nah - long);
line(nah, GetMaxY - nah - long, 0, GetMaxY - long);
line(nah, GetMaxY - nah - long, 2*nah, GetMaxY - long);
OutTextXY(0, GetMaxY - nah - long, 'z');

{0 }
OutTextXY(10, GetMaxY-7, '0');

SetColor(Yellow);
For x := 0 to dlin do
Line(fX(nah + x, nah), fY(nah, 0), fX(nah + x, nah + dlin),
fY(15 + dlin, 0 ));
SetColor(White);
end;



procedure Sphere;
var
z, I, Cosinus, a1, b1, a2, b2, a3, b3: real;
col: integer;
x, y: longint;
begin

{ sun.x := round(sun.x+sun.y/sqrt(2)); // искривление координат
sun.y := round(sun.y/sqrt(2));} // для задания вида в перспективе
for x := 0 to GetMaxX do begin
for y := 0 to GetMaxY do begin
if sqr(cir.x - x) + sqr(cir.y - y) < sqr(cir.r) then begin

z := cir.r - sqrt(sqr(cir.x - x) + sqr(cir.y - y)); // координата Z шарика вводится искуственно,
// пользователь задаёт положение шарика только координатами X и Y

{a1 := x - cir.x;
a2 := (y - cir.y);
a3 := z;
b1 := -sun.x;
b2 := sun.y;
b3 := sun.z;
if (sqrt(a1*a1 + a2*a2 + a3*a3)*sqrt(b1*b1 + b2*b2 + b3*b3))=0 then
cosinus:=0
else
cosinus:=(a1*b1 + a2*b2 + a3*b3)/
(sqrt(a1*a1 + a2*a2 + a3*a3)*sqrt(b1*b1 + b2*b2 + b3*b3));
if cosinus<0 then
cosinus:=0; }

a1 := x - cir.x;
a2 := y - cir.y;
a3 := z;
b1 := -sun.x;
b2 := sun.z;
b3 := sun.y;
if (sqrt(a1*a1 + a2*a2 + a3*a3)*sqrt(b1*b1 + b2*b2 + b3*b3))=0 then
cosinus:=0
else
cosinus:=(a1*b1 + a2*b2 + a3*b3)/
(sqrt(a1*a1 + a2*a2 + a3*a3)*sqrt(b1*b1 + b2*b2 + b3*b3));
if cosinus<0 then
cosinus:=0;


I:= cir.Kp*sun.Ip +
cir.Kh*sun.Ih*cosinus; // формула задания освещённости, взята из умной книжки

col:=Color(I);

RisNew(x,y,col);
end;
end;
end;
end;


procedure GraphInit;
var
grDriver, grMode : integer;
begin
grDriver := InstallUserDriver( 'SVGA256', nil );
grMode := 2; { 320x200, 256 ”ХЙЃјХ }
InitGraph( grDriver, grMode, '' );
end;


procedure Init_Graph;

{$F+}
function DetectVGA256:Integer;
var
GraphDriver, Graphmode: Integer;
begin {DetectVGA256}
if GraphDriver = VGA then
DetectVga256 := 1
else
DetectVGA256 := grError;
end; {DetectVGA256}
{$F-}

var
VGA256, GD, GM: Integer;

begin {Init_Graph}
VGA256 := InstallUserDriver('VGA256',@DetectVGA256);
Gd := VGA256;
InitGraph(GD,GM,'');

SetBkColor(Black);
SetColor(White);
end; {Init_Graph}

{------------------------}

begin
GraphInit; {}
{Init_Graph; {}
SetPalette;
Display; {}

sphere;

Readln;
CloseGraph;
end.

[pva]

1. Возможно дело в квантизации (неравномерная изогнутость линий). То есть заметен переход от RGB(10,10,10) к RGB(11,11,11), Такой же эффект может возникнуть при округлении (при целочисленной арифметике особенно).
2. Проверь ограничение значения косинуса угла между освещающим вектором и нормалью поверхности. Он может принимать отрицательные значения (как бы освещение изнутри). Влияет на появление артефактов освещения в тени.
3. Проверь ограничения на попадание в шар (при сканировании картинки) - освещение вне шара
4. В аттачменте примерно так оно работает (сделано в математике 4.0)
Понимаю, что текст математики выглядит мягко говоря непонятно... Если есть вопросы - охотно объясню.

[ManHack]

1) Если так, как можно исправить?
2) Для отрицательных косинусов значению косинуса присваивается нуль, это я учёл.
3) Не понял...
4) В идеале линия разделения света и тени должна проходить чётко по экватору шарика, на картинке далеко не так...
А текст там и правда не слишком понятный ><

[pva]

1) Можно использовать диферинг, то есть добавить в цвет шумовую составляющую с разбросом 1-2 значения интенсивности. Глаз "сгладит" соседние пиксели по среднему и получится эффект более плавных переходов.
3) я пробегал каждый пиксель картинки и вычислял в нём 3-мерную нормаль к поверхности. Там, где поверхности нет (вне шара) - за нормаль был принят нулевой вектор.
4) Линия и так чётко по экватору, просто надо помнить, что мы видим проекцию, то есть экватор (круг) вырождается в эллипс или отрезок.

Для пояснения переделал в линейную функцию рассеяния (чтобы лучше было видно границу на экваторе, перекрасил в красный цвет и сделал несколько кадров с разными углами освещения).

[pva]

Без отражения (то есть чисто за счёт матового рассеяния во все стороны) выглядит неестественно, добавляем отражение (будем считать что освещение идёт от бесконечного белого размазанного пятна, расположенного на расстоянии бесконечности)
1. считаем что матовое рассеяние идёт линейно n1.light, где n1 - нормаль к поверхности, light - направление падающей волны, точка - скалярное произведение.
2. считаем что пятно света с учётом всех эффектов даёт экспоненциальную функцию рассеяния. Отражённый луч вычисляем так: reflect1 = light - n1 (2(n1.light)). На экран попадает только его положительно направленная проекция на ось z, то есть reflect = Max[0,reflect1.{0,0,-1}]