OpenGL

开发平台：
Visual C++

globular.cpp：源码内容
							// qlobular.cpp
//
// Copyright (C) 2008, Celestia Development Team
// Initial code by Dr. Fridger Schrempp <fridger.schrempp@desy.de>
//
// Simulation of globular clusters, theoretical framework by  
// Ivan King, Astron. J. 67 (1962) 471; ibid. 71 (1966) 64
// 
// This program is free software; you can redistribute it and/or
// modify it under the terms of the GNU General Public License
// as published by the Free Software Foundation; either version 2
// of the License, or (at your option) any later version.
#include <fstream>
#include <algorithm>
#include <cstdio>
#include <cassert>
#include "celestia.h"
#include <celmath/mathlib.h>
#include <celmath/perlin.h>
#include <celmath/intersect.h>
#include "astro.h"
#include "globular.h"
#include <celutil/util.h>
#include <celutil/debug.h>
#include "gl.h"
#include "vecgl.h"
#include "render.h"
#include "texture.h"
#include <math.h>
using namespace std;
static int cntrTexWidth = 512, cntrTexHeight = 512; 
static int starTexWidth = 128, starTexHeight = 128;
static Color colorTable[256];
static const unsigned int GLOBULAR_POINTS  = 8192;
static const float LumiShape = 3.0f, Lumi0 = exp(-LumiShape);
// Reference values ( = data base averages) of core radius, King concentration 
// and mu25 isophote radius:  
 
static const float R_c_ref = 0.83f, C_ref = 2.1f, R_mu25 = 40.32f;  
// min/max c-values of globular cluster data
static const float MinC = 0.50f, MaxC = 2.58f, BinWidth = (MaxC - MinC) / 8.0f + 0.02f;
// P1 determines the zoom level, where individual cluster stars start to appear.
// The smaller P2 (< 1), the faster stars show up when resolution increases.
static const float P1 = 65.0f, P2 = 0.75f;
static const float RRatio_min = pow(10.0f, 1.7f);
static float CBin, RRatio, XI, DiskSizeInPixels, Rr = 1.0f, Gg = 1.0f, Bb = 1.0f;
static GlobularForm** globularForms = NULL;
static Texture* globularTex = NULL;
static Texture* centerTex[8] = {NULL};
static void InitializeForms();
static GlobularForm* buildGlobularForms(float);
static bool formsInitialized = false;
static bool decreasing (const GBlob& b1, const GBlob& b2)
{
	return (b1.radius_2d > b2.radius_2d);
}
static void GlobularTextureEval(float u, float v, float /*w*/, unsigned char *pixel)
{
	// use an exponential luminosity shape for the individual stars
    // giving sort of a halo for the brighter (i.e.bigger) stars.
    float lumi = exp(- LumiShape * sqrt(u * u + v * v)) - Lumi0;
	if (lumi <= 0.0f)
		lumi = 0.0f;
    int pixVal = (int) (lumi * 255.99f);
#ifdef HDR_COMPRESS
    pixel[0] = 127;
    pixel[1] = 127;
    pixel[2] = 127;
#else
    pixel[0] = 255;
    pixel[1] = 255;
    pixel[2] = 255;
#endif
    pixel[3] = pixVal;
}
float relStarDensity(float eta)
{
	/*! As alpha blending weight (relStarDensity) I take the theoretical
	 *  number of globular stars in 2d projection at a distance 
	 *  rho = r / r_c = eta * r_t from the center (cf. King_1962's Eq.(18)), 
	 *  divided by the area = PI * rho * rho . This number density of stars 
	 *  I normalized to 1 at rho=0.   
	 *  The resulting blending weight increases strongly -> 1 if the 
	 *  2d number density of stars rises, i.e for rho -> 0. 
	 */	 	 
     // Since the central "cloud" is due to lack of visual resolution,
	 // rather than cluster morphology, we limit it's size by
	 // taking max(C_ref, CBin). Smaller c gives a shallower distribution!
     
     float rRatio = max(RRatio_min, RRatio);
     float Xi = 1.0f / sqrt(1.0f + rRatio * rRatio);           	 
     float XI2 = Xi * Xi; 
	 float rho2 = 1.0001f + eta * eta * rRatio * rRatio; //add 1e-4 as regulator near rho=0
	 return ((log(rho2) + 4.0f * (1.0f - sqrt(rho2)) * Xi) / (rho2 - 1.0f) + XI2) / (1.0f - 2.0f * Xi + XI2);
}
static void CenterCloudTexEval(float u, float v, float /*w*/, unsigned char *pixel)
{	
	/*! For reasons of speed, calculate central "cloud" texture only for  
	 *  8 bins of King_1962 concentration, c = CBin, XI(CBin), RRatio(CBin). 
	 */
	 
	// Skyplane projected King_1962 profile at center (rho = eta = 0):	
	float c2d   = 1.0f - XI;	
		
	float eta = sqrt(u * u + v * v);  // u,v = (-1..1)
	
	// eta^2 = u * u  + v * v = 1 is the biggest circle fitting into the quadratic
	// procedural texture. Hence clipping 
		
	if (eta >= 1.0f)
		eta  = 1.0f;
	// eta = 1 corresponds to tidalRadius:
	float rho   = eta  * RRatio;
	float rho2  = 1.0f + rho * rho;		
    // Skyplane projected King_1962 profile (Eq.(14)), vanishes for eta = 1:
	// i.e. absolutely no globular stars for r > tidalRadius:	
		
	float profile_2d = (1.0f / sqrt(rho2) - 1.0f)/c2d + 1.0f ;	
	profile_2d = profile_2d * profile_2d;	
#ifdef HDR_COMPRESS
    pixel[0] = 127;
    pixel[1] = 127; 
    pixel[2] = 127;
#else    
    pixel[0] = 255;
    pixel[1] = 255;
    pixel[2] = 255;
#endif
	pixel[3] = (int) (relStarDensity(eta) * profile_2d * 255.99f);
}
Globular::Globular() :
    detail (1.0f),
    customTmpName (NULL),
    form (NULL),
	r_c (R_c_ref),
	c (C_ref),
	tidalRadius(0.0f)
{
	recomputeTidalRadius();
}
unsigned int Globular::cSlot(float conc) const
{
	// map the physical range of c, minC <= c <= maxC, 
	// to 8 integers (bin numbers), 0 < cSlot <= 7:
	
	if (conc <= MinC)
		conc  = MinC;
	if (conc >= MaxC)
		conc  = MaxC;
  	 
    return (unsigned int) floor((conc - MinC) / BinWidth);
}
		
const char* Globular::getType() const
{
    return "Globular";
}
void Globular::setType(const std::string& /*typeStr*/)
{
}
float Globular::getDetail() const
{
    return detail;
}
void Globular::setDetail(float d)
{
    detail = d;
}
string Globular::getCustomTmpName() const
{
    if (customTmpName == NULL)
        return "";
    else
        return *customTmpName;
}
void Globular::setCustomTmpName(const string& tmpNameStr)
{
    if (customTmpName == NULL)
        customTmpName = new string(tmpNameStr);
    else
        *customTmpName = tmpNameStr;
}
float Globular::getCoreRadius() const
{
    return r_c;
}
void Globular::setCoreRadius(const float coreRadius)
{
	r_c = coreRadius;
	recomputeTidalRadius();
}
	
float Globular::getHalfMassRadius() const
{
	// Aproximation to the half-mass radius r_h [ly]
	// (~ 20% accuracy) 
    
	return std::tan(degToRad(r_c / 60.0f)) * (float) getPosition().distanceFromOrigin() *                                 pow(10.0f, 0.6f * c - 0.4f); 	  
}
float Globular::getConcentration() const
{
	return c;
}
void Globular::setConcentration(const float conc)
{
    c = conc;
    if (!formsInitialized)
		InitializeForms();
	
	// For saving time, account for the c dependence via 8 bins only,
	form = globularForms[cSlot(conc)];
    recomputeTidalRadius();
}
size_t Globular::getDescription(char* buf, size_t bufLength) const
{
	return snprintf(buf, bufLength, _("Globular (core radius: %4.2f', King concentration: %4.2f)"), r_c, c);
}
GlobularForm* Globular::getForm() const
{
    return form;
}
const char* Globular::getObjTypeName() const
{
    return "globular";
}
static const float RADIUS_CORRECTION = 0.025f;
bool Globular::pick(const Ray3d& ray,
                  double& distanceToPicker,
                  double& cosAngleToBoundCenter) const
{
    if (!isVisible())
        return false;
	/*
     * The selection ellipsoid should be slightly larger to compensate for the fact
     * that blobs are considered points when globulars are built, but have size
     * when they are drawn. 
	 */
	Vec3d ellipsoidAxes(getRadius() * (form->scale.x + RADIUS_CORRECTION),
                        getRadius() * (form->scale.y + RADIUS_CORRECTION),
                        getRadius() * (form->scale.z + RADIUS_CORRECTION));
    Quatf qf= getOrientation();
    Quatd qd(qf.w, qf.x, qf.y, qf.z);
    return testIntersection(Ray3d(Point3d() + (ray.origin - getPosition()), ray.direction) * conjugate(qd).toMatrix3(),
                            Ellipsoidd(ellipsoidAxes),
                            distanceToPicker,
                            cosAngleToBoundCenter);
}
bool Globular::load(AssociativeArray* params, const string& resPath)
{
    // Load the basic DSO parameters first
    bool ok = DeepSkyObject::load(params, resPath);
    if (!ok)
		return false;
    if(params->getNumber("Detail", detail))
    	setDetail((float) detail);
			
	string customTmpName;
	if(params->getString("CustomTemplate", customTmpName ))
		setCustomTmpName(customTmpName);
			
	if(params->getNumber("CoreRadius", r_c))
		setCoreRadius(r_c);
			
	if(params->getNumber("KingConcentration", c))
		setConcentration(c);
	
    return true;
}
void Globular::render(const GLContext& context,
                    const Vec3f& offset,
                    const Quatf& viewerOrientation,
                    float brightness,
                    float pixelSize)
{					        
	renderGlobularPointSprites(context, offset, viewerOrientation, brightness, pixelSize);
}
void Globular::renderGlobularPointSprites(const GLContext&,
                                      const Vec3f& offset,
                                      const Quatf& viewerOrientation,
                                      float brightness,
                                      float pixelSize)
{
    if (form == NULL)
        return;
			
	float distanceToDSO = offset.length() - getRadius();
	if (distanceToDSO < 0)
    	distanceToDSO = 0;
			
	float minimumFeatureSize = 0.5f * pixelSize * distanceToDSO;    
	
	DiskSizeInPixels = getRadius() / minimumFeatureSize;
	
	/*
     * Is the globular's apparent size big enough to
     * be noticeable on screen? If it's not, break right here to
     * avoid all the overhead of the matrix transformations and
	 * GL state changes: 
	 */
	
	if (DiskSizeInPixels < 1.0f)	
		return;	    			
    /*
	 * When resolution (zoom) varies, the blended texture opacity is controlled by the 
	 * factor 'pixelWeight'. At low resolution, the latter starts at 1, but tends to 0,  
     * if the resolution increases sufficiently (DiskSizeInPixels >= P1 pixels)! 
	 * The smaller P2 (<1), the faster pixelWeight -> 0, for DiskSizeInPixels >= P1.
	 */	
			
	float pixelWeight = (DiskSizeInPixels >= P1)? 1.0f/(P2 + (1.0f - P2) * DiskSizeInPixels / P1): 1.0f;
							
    // Use same 8 c-bins as in globularForms below!        
    
    unsigned int ic = cSlot(c); 
    CBin = MinC + ((float) ic + 0.5f) * BinWidth; // center value of (ic+1)th c-bin
	RRatio = pow(10.0f, CBin); 
    XI = 1.0f / sqrt(1.0f + RRatio * RRatio);           
    
    if(centerTex[ic] == NULL)
	{
		centerTex[ic] = CreateProceduralTexture( cntrTexWidth, cntrTexHeight, GL_RGBA, CenterCloudTexEval);
	}
	assert(centerTex[ic] != NULL);
		
	if (globularTex == NULL)
    {
        globularTex = CreateProceduralTexture( starTexWidth, starTexHeight, GL_RGBA,
                                               GlobularTextureEval);
    }
    assert(globularTex != NULL);
	glEnable (GL_BLEND);
	glEnable (GL_TEXTURE_2D);
	glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
					    
    Mat3f viewMat = viewerOrientation.toMatrix3();
    Vec3f v0 = Vec3f(-1, -1, 0) * viewMat;
    Vec3f v1 = Vec3f( 1, -1, 0) * viewMat;
    Vec3f v2 = Vec3f( 1,  1, 0) * viewMat;
    Vec3f v3 = Vec3f(-1,  1, 0) * viewMat;
	
	float tidalSize = 2 * tidalRadius;
	Mat3f m =
        Mat3f::scaling(form->scale) * getOrientation().toMatrix3() * 
		Mat3f::scaling(tidalSize);
				     
    vector<GBlob>* points = form->gblobs;
    unsigned int nPoints = 
				 (unsigned int) (points->size() * clamp(getDetail()));
				 				 	
	/* Render central cloud sprite (centerTex). It fades away when
	 * distance from center or resolution increases sufficiently. 	
	 */
	
	centerTex[ic]->bind();	
    
	float br = 2 * brightness;
	glColor4f(Rr, Gg, Bb, min(br * pixelWeight, 1.0f));	
	glBegin(GL_QUADS);		
		
	glTexCoord2f(0, 0);          glVertex(v0 * tidalSize);
    glTexCoord2f(1, 0);          glVertex(v1 * tidalSize);
    glTexCoord2f(1, 1);          glVertex(v2 * tidalSize);
    glTexCoord2f(0, 1);          glVertex(v3 * tidalSize);   					
		
	glEnd();						     
			 				
	/*! Next, render globular cluster via distinct "star" sprites (globularTex) 
	 * for sufficiently large resolution and distance from center of globular.
	 *  
	 * This RGBA texture fades away when resolution decreases (e.g. via automag!), 
	 * or when distance from globular center decreases.
	 */
	
	
	globularTex->bind();
	
	
	int pow2 = 128;			     // Associate "Red Giants" with the 128 biggest star-sprites  
	float starSize =  br * 0.5f; // Maximal size of star sprites -> "Red Giants"
	float clipDistance = 100.0f; // observer distance [ly] from globular, where we 
		                         // start "morphing" the star-sprite sizes towards 
			                     // their physical values
	glBegin(GL_QUADS);
    
	for (unsigned int i = 0; i < nPoints; ++i)
	{  						  
		GBlob         b  = (*points)[i];
		Point3f       p  = b.position * m;		
		float    eta_2d  = b.radius_2d;
		
		/*! Note that the [axis,angle] input in globulars.dsc transforms the 
		 *  2d projected star distance r_2d in the globular frame to refer to the 
		 *  skyplane frame for each globular! That's what I need here.
	     *
	     *  The [axis,angle] input will be needed anyway, when upgrading to 
		 *  account for  ellipticities, with corresponding  inclinations and 
		 *  position angles...
	     */
		if ((i & pow2) != 0)
		{      
			pow2    <<=   1;		   
			starSize /= 1.25f;			
	   			
			if (starSize < minimumFeatureSize)
           		break;
		} 
		float obsDistanceToStarRatio = (p + offset).distanceFromOrigin() / clipDistance;		
		float saveSize = starSize;
		if (obsDistanceToStarRatio < 1.0f)
		{  
			// "Morph" the star-sprite sizes at close observer distance such that
			// the overdense globular core is dissolved upon closing in.
            
			starSize = starSize * min(obsDistanceToStarRatio, 1.0f);     
		}	
         
		/* Colors of normal globular stars are given by color profile.
		 * Associate orange "Red Giant" stars with the largest sprite
		 * sizes (while pow2 = 128).
		 */
		Color col  = (pow2 < 256)? colorTable[255]: colorTable[b.colorIndex];		
		glColor4f(col.red(), col.green(), col.blue(), 
			      min(br * (1.0f - pixelWeight * relStarDensity(eta_2d)), 1.0f));
	   
		glTexCoord2f(0, 0);          glVertex(p + (v0 * starSize));
		glTexCoord2f(1, 0);          glVertex(p + (v1 * starSize));
		glTexCoord2f(1, 1);          glVertex(p + (v2 * starSize));
		glTexCoord2f(0, 1);          glVertex(p + (v3 * starSize));
	   					        
		starSize = saveSize;	
	}
	glEnd();
}
unsigned int Globular::getRenderMask() const
{
    return Renderer::ShowGlobulars;
}
unsigned int Globular::getLabelMask() const
{
    return Renderer::GlobularLabels;
}
		
void Globular::recomputeTidalRadius()
{
	// Convert the core radius from arcminutes to light years
	// Compute the tidal radius in light years
        
	float coreRadiusLy = std::tan(degToRad(r_c / 60.0f)) * (float) getPosition().distanceFromOrigin();
    tidalRadius = coreRadiusLy * std::pow(10.0f, c);     
}
GlobularForm* buildGlobularForms(float c)
{
	GBlob b;
	vector<GBlob>* globularPoints = new vector<GBlob>;
	float rRatio = pow(10.0f, c); //  = r_t / r_c
	float prob;
	float cc =  1.0f + rRatio * rRatio;
	unsigned int i = 0, k = 0;
	
	// Value of King_1962 luminosity profile at center:	
	
	float prob0 = sqrt(cc) - 1.0f;
	/*! Generate the globular star distribution randomly, according 
	 *  to the  King_1962 surface density profile f(r), eq.(14).
	 *
	 *  rho = r / r_c = eta r_t / r_c, 0 <= eta <= 1,  
  	 *  coreRadius r_c, tidalRadius r_t, King concentration c = log10(r_t/r_c). 
	 */	
	
	while (i < GLOBULAR_POINTS)
	{	
	    /*!  
		 *  Use a combination of the Inverse Transform method and 
		 *  Von Neumann's Acceptance-Rejection method for generating sprite stars 
		 *  with eta distributed according to the exact King luminosity profile.    
		 *
		 * This algorithm leads to almost 100% efficiency for all values of  
		 * parameters and variables!
		 */
	
		float uu = Mathf::frand();
		
		/* First step: eta distributed as inverse power distribution (~1/Z^2) 
		 * that majorizes the exact King profile. Compute eta in terms of uniformly 
		 * distributed variable uu! Normalization to 1 for eta -> 0.
         */ 				
		
		float eta = tan(uu *atan(rRatio))/rRatio;
		float rho = eta * rRatio;		
		float  cH = 1.0f/(1.0f + rho * rho);
        float   Z = sqrt((1.0f + rho * rho)/cc); // scaling variable
		
	    // Express King_1962 profile in terms of the UNIVERSAL variable 0 < Z <= 1,		
		
		prob = (1.0f - 1.0f / Z) / prob0;
		prob = prob * prob;
								
  	    /* Second step: Use Acceptance-Rejection method (Von Neumann) for
		 * correcting the power distribution of eta into the exact, 
		 * desired King form 'prob'!
		 */
        k++;
		if (Mathf::frand() < prob / cH)
		{
		   	/* Generate 3d points of globular cluster stars in polar coordinates:
			 * Distribution in eta (<=> r) according to King's profile. 
		     * Uniform distribution on any spherical surface for given eta. 
			 * Note: u = cos(phi) must be used as a stochastic variable to get uniformity in angle!
			 */
			float u = Mathf::sfrand(); 
			float theta = 2 * (float) PI * Mathf::frand();
			float sthetu2 = sin(theta) * sqrt(1.0f - u * u);
			
			// x,y,z points within -0.5..+0.5, as required for consistency: 
			b.position = 0.5f * Point3f(eta * sqrt(1.0f - u * u) * cos(theta), eta * sthetu2 , eta * u);
			            
			/* 
			 * Note: 2d projection in x-z plane, according to Celestia's
			 * conventions! Hence...
			 */			
			b.radius_2d = eta * sqrt(1.0f - sthetu2 * sthetu2);
			/* For now, implement only a generic spectrum for normal cluster
			 * stars, modelled from Hubble photo of M80. 
			 * Blue Stragglers are qualitatively accounted for...
   			 * assume color index poportional to Z as function of which the King profile
			 * becomes universal!
			 */
			  
			b.colorIndex  = (unsigned int) (Z * 254);
            
		    globularPoints->push_back(b);					
			i++;		
		}
	}
	
	// Check for efficiency of sprite-star generation => close to 100 %!
	//cout << "c =  "<< c <<"  i =  " << i - 1 <<"  k =  " << k - 1 <<                                                       "  Efficiency:  " << 100.0f * i / (float)k<<"%" << endl;	
	
	GlobularForm* globularForm  = new GlobularForm();
	globularForm->gblobs        = globularPoints;
	globularForm->scale         = Vec3f(1.0f, 1.0f, 1.0f);	
		
	return globularForm;
}
void InitializeForms()
{
	// Build RGB color table, using hue, saturation, value as input.
	// Hue in degrees.
	
	// Location of hue transition and saturation peak in color index space: 
	int i0 = 36, i_satmax = 16; 
	// Width of hue transition in color index space: 	
	int i_width = 3;
				
	float sat_l = 0.08f, sat_h = 0.1f, hue_r = 27.0f, hue_b = 220.0f;
	
	// Red Giant star color: i = 255:
	// -------------------------------
	// Convert hue, saturation and value to RGB
	DeepSkyObject::hsv2rgb(&Rr, &Gg, &Bb, 25.0f, 0.65f, 1.0f);
	colorTable[255]  = Color(Rr, Gg, Bb);
		
	// normal stars: i < 255, generic color profile for now, improve later
	// --------------------------------------------------------------------
    // Convert hue, saturation, value to RGB
		
	for (int i = 254; i >=0; i--)
    {
		// simple qualitative saturation profile:	
		// i_satmax is value of i where sat = sat_h + sat_l maximal
	    float x = (float) i / (float) i_satmax, x2 = x ;
		float sat = sat_l + 2 * sat_h /(x2 + 1.0f / x2);
				
		// Fast transition from hue_r to hue_b at i = i0 within a width
		// i_width in color index space:
		float hue = hue_r + 0.5f * (hue_b - hue_r) * (std::tanh((float)(i - i0) / (float) i_width) + 1.0f);
		
		DeepSkyObject::hsv2rgb(&Rr, &Gg, &Bb, hue, sat, 0.85f);
		colorTable[i]  = Color(Rr, Gg, Bb);
	}
	// Define globularForms corresponding to 8 different bins of King concentration c 
   
	globularForms = new GlobularForm*[8];   
	for (unsigned int ic  = 0; ic <= 7; ++ic)
    {
        float CBin = MinC + ((float) ic + 0.5f) * BinWidth;
		globularForms[ic] = buildGlobularForms(CBin); 
	}				
	formsInitialized = true;
}

						