#include <stdio.h>
#include <math.h>

double f(double x, double y, double z);
double fx(double x, double y, double z);
double fy(double x, double y, double z);
double fz(double x, double y, double z);
double compute_y(double x, double z);

int main(int argc, char const *argv[])
{
	int width = 512, height = 512;
	FILE *fp = fopen("heart3d.ppm","w");
	fprintf(fp, "P3 %d %d 255\n", width, height);
	double x, y, z, nx, ny, nz, a;
	// 光源方向（入射光方向）
	double lx = 1, ly = 1, lz = -1;
	double L = sqrt(lx*lx+ly*ly+lz*lz);
	lx = lx / L;
	ly = ly / L;
	lz = lz / L;

	for (int h=0; h<height; h++)
	{
		// z \in [-1.5, 1.5]
		z = - h*3.0/height + 1.5;
		for (int w=0; w<width; w++)
		{
			// z \in [-1.5, 1.5]
			x = w*3.0/width - 1.5;
			int r = 0;

			// 点在心形图内
			if (f(x,0,z) < 0)
			{
				// 计算y
				y = compute_y(x,z);
				// 计算梯度（法向量）
				nx = fx(x,y,z);
				ny = fy(x,y,z);
				nz = fz(x,y,z);
				// light source (1,1,-0.5)
				a = (lx*nx + ly*ny + lz*nz)/sqrt(nx*nx+ny*ny+nz*nz);
				// 根据入射光与法线夹角余弦值计算光强
				r = (int) ((0.5 * a + 0.5)*255);	
			}
			
			fprintf(fp, "%d 0 0 ", r);
		}
		fprintf(fp,"\n");
	}
	return 0;
}

double f(double x, double y, double z){
	double a = x*x + 9.0f/4.0f*y*y + z*z - 1;
	return a*a*a - x*x*z*z*z - 9.0f/80.0f*y*y*z*z*z;
}

double fx(double x, double y, double z){
	double dx = 1e-8;
	return (f(x+dx,y,z) - f(x-dx,y,z)) / (2*dx);
}

double fy(double x, double y, double z){
	double dy = 1e-8;
	return (f(x,y+dy,z) - f(x,y-dy,z)) / (2*dy);
}

double fz(double x, double y, double z){
	double dz = 1e-8;
	return (f(x,y,z+dz) - f(x,y,z-dz)) / (2*dz);
}

double compute_y(double x, double z){
	double z_min = 0.0, z_max = 1.5;
	double eps = 1e-6, r = 1.0, y;
	while(r>eps){
		y = (z_min + z_max) / 2.0;
		if (f(x,y,z) == 0)
			return y;
		if (f(x,y,z) > 0)
			z_max = y;
		else
			z_min = y;
		r = z_max - z_min;
	}

	return 0.5*(z_max + z_min);
}