MATLAB_视频处理学习记录_03

#include <stdio.h>
#include <stdlib.h>

static double fact( int n ) {  // The factorial function!
double f = 1;
while(n>1)
f *= n--;
return f;
}

static inline unsigned char f2b( double f ) {  // float to byte
if(f<0) f = 0; if(f>1) f = 1;
int i = f*256; if(i>255) i = 255;
return i;
}

void HSV2RGB(double h, double s, double v, unsigned char rgb[3] ) {  // the classic hue scale
if (s == 0) {
rgb[0] = rgb[1] = rgb[2] = f2b(v);
} else {
double v_h = h * 6;
double v_i = floor(v_h);
double v_1 = v * (1 - s);
double v_2 = v * (1 - s * (v_h - v_i));
double v_3 = v * (1 - s * (1 - (v_h - v_i)));
double v_r,v_g,v_b;
if (v_i == 0) {v_r = v;   v_g = v_3; v_b = v_1;}
else if (v_i == 1) {v_r = v_2; v_g = v;   v_b = v_1;}
else if (v_i == 2) {v_r = v_1; v_g = v;   v_b = v_3;}
else if (v_i == 3) {v_r = v_1; v_g = v_2; v_b = v  ;}
else if (v_i == 4) {v_r = v_3; v_g = v_1; v_b = v  ;}
else               {v_r = v;   v_g = v_1; v_b = v_2;};

rgb[0] = f2b(v_r);
rgb[1] = f2b(v_g);
rgb[2] = f2b(v_b);
}
}

double zernike_polynomials( int m, int n, double ro, double th ) {  // the polynomial hitself
double Rmnro = 0;
bool even = m>=0;
if(m<0) m = -m;
if( (n-m)%2 ) return 0;

for(int k=0;k<=(n-m)/2;++k) {
Rmnro += pow(ro,n-2*k)*
( pow(-1,k) * fact(n-k) ) /
( fact(        k) *
fact((n+m)/2-k) *
fact((n-m)/2-k)
);
}

if(even) return Rmnro * cos(m*th);
else     return Rmnro * sin(m*th);
}

void main() {
const double NO_VALUE = 42; const int B  = 64;
const int SX = 1024;
const int SY = 1024;

const int LX = 3; const int LY = 7;
int pairs[21][2] =
{
{ 0,0},
{-1,1},{ 1,1},
{-2,2},{ 0,2},{ 2,2},
{-3,3},{-1,3},{ 1,3},{ 3,3},
{-4,4},{-2,4},{ 0,4},{ 2,4},{ 4,4},
{-5,5},{-3,5},{-1,5},{ 1,5},{ 3,5},{ 5,5},
};

unsigned char * img = new unsigned char[LX*SX*LY*SY*3];

int ix,iy,q = 0;
for(iy=0;iy<LY;++iy) for(ix=0;ix<LX;++ix) {
int m = pairs[q][0];
int n = pairs[q][1];
for(int j=0;j<SY;++j) {
double y = double(SY/2-j)/(SY/2-B);
for(int i=0;i<SX;++i) {
double x = double(SX/2-i)/(SX/2-B);
double ro = sqrt(x*x+y*y);
double th = atan2(y,x);
double z = NO_VALUE;
if(ro<=1)
z = zernike_polynomials(m,n,ro,th);
unsigned char * rgb = img+3*( i+SX*(ix+LX*(j+SY*iy)) );
if(z == NO_VALUE )
rgb[0] = rgb[1] = rgb[2] = 255;
else {
if(z>=0) z =  pow( z,0.8);      // little enance
else     z = -pow(-z,0.8);
z = (z+1)/2;                    // normalize
if(z<0) z = 0;
if(z>1) z = 1;

HSV2RGB( 2.0*z/3.0, 0.8, 0.9, rgb );
}
}
}
printf("%03d%%\n",q*100/(LX*LY)); ++q;
}

FILE * fp = fopen("c:\\temp\\zernike.ppm","wb");
fprintf(fp,"P6\n%d %d\n255\n",SX*LX,SY*LY);
fwrite(img,1,LX*SX*LY*SY*3,fp);
fclose(fp);
delete[] img;
}