// PlanetView.cpp: implementation of the CPlanetView class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include <math.h>

#include "PlanetView.h"


#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

CPlanetView::CPlanetView()
{
	m_viewerDistance = 0;
	m_viewerAngle = 0;
	m_lightAngle = 0;
}

CPlanetView::~CPlanetView()
{
}

void CPlanetView::SetViewerPosition(double angle, double distance)
{
	m_viewerDistance = distance;
	m_viewerAngle = angle;
}

void CPlanetView::SetLightPosition(double angle)
{
	m_lightAngle = angle;
}

void CPlanetView::Render(HBITMAP color, HBITMAP topography, double heightRange, HBITMAP canvas)
{
	BITMAP canvasInfo;
	GetObject(canvas, sizeof(canvasInfo), &canvasInfo);

	// Locate the projection plane so that the image fills the
	// canvas. Projection plane distance is measured from the
	// viewer.
	LONG imageDiameter;
	if(canvasInfo.bmWidth < canvasInfo.bmHeight)
	{
		imageDiameter = canvasInfo.bmWidth;
	}
	else
	{
		imageDiameter = canvasInfo.bmHeight;
	}
	double projectionDistance =
		(imageDiameter / 2) * m_viewerDistance / sqrt(1 - 1 / m_viewerDistance / m_viewerDistance);

	LONG x, y;
	for(y = 0; y < canvasInfo.bmHeight; y++)
	{
		double y0 = y - 0.5 * canvasInfo.bmHeight;
		for(x = 0; x < canvasInfo.bmWidth; x++)
		{
			double x0 = x - 0.5 * canvasInfo.bmWidth;
			COLORREF pixelColor;

			// Solve a quadratic equation to determine where
			// the line of sight intersects the sphere.
			double a = x0 * x0 + y0 * y0 + projectionDistance * projectionDistance;
			double b = projectionDistance * m_viewerDistance;
			double c = m_viewerDistance * m_viewerDistance - 1;

			if(b * b >= a * c)
			{
				double t = (b - sqrt(b * b - a * c)) / a;
				
				// Calculate the latitude and longitude of the intersection.
				double latitude = asin(y0 * t);
				double longitude = m_viewerAngle + asin(x0 * t / cos(latitude));
				pixelColor = GetIlluminatedColor(
					color,
					topography,
					heightRange,
					latitude,
					longitude);
			}
			else
			{
				// This pixel is not part of the projection.
				pixelColor = 0;
			}

			// Set the pixel's color in the canvas.
			DWORD canvasPitch = (canvasInfo.bmWidth * 3 + 3) / 4 * 4;
			BYTE *pixelAddress = (BYTE *)canvasInfo.bmBits + y * canvasPitch + x * 3;
			pixelAddress[0] = (BYTE) (pixelColor & 0xff);
			pixelAddress[1] = (BYTE) ((pixelColor >> 8) & 0xff);
			pixelAddress[2] = (BYTE) ((pixelColor >> 16) & 0xff);
		}
	}
}

COLORREF CPlanetView::GetIlluminatedColor(HBITMAP color, HBITMAP topography, double heightRange, double latitude, double longitude)
{
	COLORREF surfaceColor = GetSurfaceColor(color, latitude, longitude);
	COLORREF illuminatedColor;
	
	double normal[3];
	GetSurfaceNormal(
		topography,
		heightRange,
		latitude,
		longitude,
		normal);
	double normalMagnitude =
		sqrt(normal[0] * normal[0] + normal[1] * normal[1] + normal[2] * normal[2]);
	normal[0] /= normalMagnitude;
	normal[1] /= normalMagnitude;
	normal[2] /= normalMagnitude;
	
	double lightSource[3] =
	{
		sin(m_lightAngle),
		0,
		cos(m_lightAngle)
	};
	double incomingCosine =
		(lightSource[0] * normal[0] +
		lightSource[1] * normal[1] +
		lightSource[2] * normal[2]);

	if(incomingCosine > 0)
	{
		double viewer[3] =
		{
			sin(m_viewerAngle),
			0,
			cos(m_viewerAngle)
		};

		double outgoingCosine =
			(viewer[0] * normal[0] +
			viewer[1] * normal[1] +
			viewer[2] * normal[2]);

		double diffuseFactor = 2.5 * pow(incomingCosine, 0.5) * outgoingCosine;
		int red   = (int) (diffuseFactor * (surfaceColor & 0xff));
		if(red > 255)
		{
			red = 255;
		}
		int green = (int) (diffuseFactor * ((surfaceColor >> 8) & 0xff));
		if(green > 255)
		{
			green = 255;
		}
		int blue  = (int) (diffuseFactor * ((surfaceColor >> 16) & 0xff));
		if(blue > 255)
		{
			blue = 255;
		}

		illuminatedColor = (blue << 16) | (green << 8) | red;
	}
	else
	{
		illuminatedColor = 0;
	}

	return illuminatedColor;
}

COLORREF CPlanetView::GetSurfaceColor(HBITMAP color, double latitude, double longitude)
{
	longitude = fmod(longitude, 2 * m_pi);

	BITMAP bmColor;
	int iResult = GetObject(color, sizeof(bmColor), &bmColor);
	long x = (long) (longitude / (2 * m_pi) * bmColor.bmWidth);
	long y = (long) ((0.5 + latitude / m_pi) * bmColor.bmHeight);

	return *(DWORD *)((BYTE *)bmColor.bmBits + y * bmColor.bmWidthBytes + x * 3) & 0x00ffffff;
}

void CPlanetView::GetSurfaceNormal(HBITMAP topography, double heightRange, double latitude, double longitude, double normal[3])
{
	longitude = fmod(longitude, 2 * m_pi);

	BITMAP bmTopo;
	int iResult = GetObject(topography, sizeof(bmTopo), &bmTopo);
	long x = (long) (longitude / (2 * m_pi) * bmTopo.bmWidth);
	long y = (long) ((0.5 + latitude / m_pi) * bmTopo.bmHeight);

	long xWest  = (x + bmTopo.bmWidth - 1) % bmTopo.bmWidth;
	long xEast  = (x + 1) % bmTopo.bmWidth;
	long yNorth = (y < bmTopo.bmHeight)? y + 1 : y;
	long ySouth = (y > 0)? y - 1 : 0;

	BYTE *bmPixels = (BYTE *)bmTopo.bmBits + 256 * sizeof(RGBQUAD);
	int xHeightChange = 
		(*(bmPixels + y * bmTopo.bmWidthBytes + xEast) -
		*(bmPixels + y * bmTopo.bmWidthBytes + xWest)) / 2;
	int yHeightChange = 
		(*(bmPixels + yNorth * bmTopo.bmWidthBytes + x) -
		*(bmPixels + ySouth * bmTopo.bmWidthBytes + x)) / 2;

	double dTheta = 360.0 / bmTopo.bmWidth;
	double dPhi = 180.0 / bmTopo.bmHeight;

	double localNormal[3];
	localNormal[0] = xHeightChange * (heightRange / 255) / cos(latitude) / dTheta;
	localNormal[1] = yHeightChange * (heightRange / 255) / dPhi;
	localNormal[2] = 1;
	double localMagnitude =
		sqrt(localNormal[0] * localNormal[0] +
		localNormal[1] * localNormal[1] +
		localNormal[2] * localNormal[2]);

    double radiusVector[3];
	radiusVector[0] = cos(latitude) * sin(longitude);
	radiusVector[1] = sin(latitude);
	radiusVector[2] = cos(latitude) * cos(longitude);

	normal[0] = radiusVector[0] + (localNormal[0] / localMagnitude);
	normal[1] = radiusVector[1] + (localNormal[1] / localMagnitude);
	normal[2] = radiusVector[2] + (localNormal[2] / localMagnitude - 1);
}

const double CPlanetView::m_pi = 3.14159;
