Раскраска графического представления полинома Цернике

//---------------------------------------------------------------------------
 
#include <vcl.h>
#pragma hdrstop
 
#include "Unit1.h"
#include <math.h>
#include <windows.h>
//---------------------------------------------------------------------------
#pragma package(smart_init)
#pragma resource "*.dfm"
TForm1 *Form1;
 
float pows(float x, float y)
{
if( y == 0) return 1;
  return pow(x, y);
 
 
 
}
long F(int k)
{
        long f = 1;
        for ( int i = 2; i <= k; i++)
                f *= i;
        return f;
 
 
}
 
 
float R(int m, int n, float rho)
{
        if((n - m) % 2 != 0) return 0.0f;
        float s = 0.0f;
        for (int k = 0; k <= (n - m) / 2; k++)
        {
                long t = (long) pows(-1.0f, (float) k) * F(n - k);
                long d = F(k) *  F((n + m) / 2 - k) * F((n - m) / 2 - k);
                s += (float) (t / d) * (float) pows(rho, (float) (n - 2 * k));
        }
        return s;
}
 
float Z(int m, int n, float rho, float phi)
{
    
if(m < 0 ) return R(-m, n, rho) * sin((float) (m) * phi); 
else return R(m, n, rho) * - cos((float) m * phi);
 
 
 
}
void setPixel(int x, int y, TColor c)
{
 
 Form1->Image1->Canvas->Pen->Color = c;
Form1->Image1->Canvas->Pen->Width = 2;
 
 Form1->Image1->Canvas->MoveTo(x, y);
 Form1->Image1->Canvas->LineTo(x, y);
 
 
 
}
 
 void draw(int m, int n, int cx, int cy, int r)
        {
            int x0 = cx - r;
            int y0 = cy - r;
            int x1 = cx + r;
            int y1 = cy + r;
            int rsq = r * r;
 
            for (int y = y0; y <= y1; y++)
            {
                int dy = y - cy;
                int dysq = dy * dy;
                for (int x = x0; x <= x1; x++)
                {
                    int dx = x - cx;
                    int dsq = dx * dx + dysq;
                    if(dsq > rsq) continue;
 
                    float len = (float) sqrt(dsq);
 
                    float rho = len / (float) r;
 
                    float phi;
                    if( dx == 0 && dy == 0) phi = 0;
                    else
                    {
                        phi = (float) acos((float) dx / len);
                        if( dy < 0) phi = 2.0f * (float) M_PI - phi;
                    }
 
                    float z = Z(m, n, rho, phi);
 
                    int c = (int) ((0.5f + 0.5f * z) * 255.0f);
 
                    if (dx * dx + dysq > rsq) continue;
                    setPixel(x, y, RGB(c, c, c));
 
                }
 
 
 
            }
 
        }
 
//---------------------------------------------------------------------------
__fastcall TForm1::TForm1(TComponent* Owner)
        : TForm(Owner)
{
 
}
//---------------------------------------------------------------------------
 
void __fastcall TForm1::FormCreate(TObject *Sender)
{
        int size = 50;
      int cx = 400;
      for(int y = 0; y <= 5; y++){
       for(int x = -y; x <= y; x+=2) {
 
            draw(x, y, (x * size) + cx, y * size + (2 * size), 20);
       }
 
       }
}
//---------------------------------------------------------------------------