|
|
楼主 |
发表于 2011-4-28 13:57:03
|
显示全部楼层
以下是完整的代码：
#include "WavesSpectrumPass.h"
#include "ShaderManager.h"
#include <osg/Camera>
#include <osg/FrameBufferObject>
#include <osg/Geode>
#include <osg/geometry>
#include <math.h>
#ifndef M_PI
#define M_PI 3.14159265
#endif
static void drawQuad()
{
glBegin(GL_TRIANGLE_STRIP);
glVertex4f(-1.0, -1.0, 0.0, 0.0);
glVertex4f(+1.0, -1.0, 1.0, 0.0);
glVertex4f(-1.0, +1.0, 0.0, 1.0);
glVertex4f(+1.0, +1.0, 1.0, 1.0);
glEnd();
}
class Quad : public osg: rawable
{
public:
Quad() {}
/** Copy constructor using CopyOp to manage deep vs shallow copy.*/
Quad(const Quad& quad,const osg::CopyOp& copyop=osg::CopyOp::SHALLOW_COPY):
osg::Drawable(quad,copyop) {}
META_Object(zchxOcean,Quad)
virtual void drawImplementation(osg::RenderInfo&) const
{
drawQuad();
}
protected:
virtual ~Quad() {}
};
WavesSpectrumPass::WavesSpectrumPass(int textureSize,float windSpeed,float inverseWaveAge,float amplitude, float GRID1_SIZE,float GRID2_SIZE,float GRID3_SIZE,float GRID4_SIZE):
theTextureSize (textureSize), //纹理尺寸大小
theDirtyFlag (false),
theWindSpeed (windSpeed),//风速 wind speed in meters per second (at 10m above surface)
theInverseWaveAge (inverseWaveAge),//逆波龄 sea state (inverse wave age)
theAmplitude (amplitude),//波浪振幅因子 wave amplitude factor (should be one)
theSpectrum12Tex (NULL),
theSpectrum34Tex (NULL),
theSpectrum12 (NULL),
theSpectrum34 (NULL),
N_SLOPE_VARIANCE (10),
theGRID1_SIZE (GRID1_SIZE),
theGRID2_SIZE (GRID2_SIZE),
theGRID3_SIZE (GRID3_SIZE),
theGRID4_SIZE (GRID4_SIZE),
cm (0.23f),
km (370.0f)
{
thePasses = 8;
theRootGroup = new osg::Group;
createTextures();
createShader();
build();
}
WavesSpectrumPass::~WavesSpectrumPass()
{
}
osg::ref_ptr<osg::Geode> WavesSpectrumPass::createQuad()
{
osg::ref_ptr<osg::Geode> quadGeode = new osg::Geode;
quadGeode->addDrawable(new Quad);
return quadGeode;
}
void WavesSpectrumPass::createTextures()
{
//createTextures
theSpectrum12Tex = new osg::Texture2D;
theSpectrum12Tex->setDataVariance(osg::Object::DYNAMIC); // protect from being optimized away as static state.
theSpectrum12Tex->setTextureSize(theTextureSize, theTextureSize);
theSpectrum12Tex->setInternalFormat(GL_RGBA16F_ARB);
theSpectrum12Tex->setSourceFormat(GL_RGB);
theSpectrum12Tex->setSourceType(GL_FLOAT);
theSpectrum12Tex->setFilter(osg::Texture2D::MIN_FILTER,osg::Texture2D::NEAREST);
theSpectrum12Tex->setFilter(osg::Texture2D::MAG_FILTER,osg::Texture2D::NEAREST);
theSpectrum12Tex->setWrap(osg::Texture2D::WRAP_S,osg::Texture2D::REPEAT);
theSpectrum12Tex->setWrap(osg::Texture2D::WRAP_T,osg::Texture2D::REPEAT);
theSpectrum34Tex = new osg::Texture2D;
theSpectrum34Tex->setDataVariance(osg::Object::DYNAMIC); // protect from being optimized away as static state.
theSpectrum34Tex->setTextureSize(theTextureSize, theTextureSize);
theSpectrum34Tex->setInternalFormat(GL_RGBA16F_ARB);
theSpectrum34Tex->setSourceFormat(GL_RGB);
theSpectrum34Tex->setSourceType(GL_FLOAT);
theSpectrum34Tex->setFilter(osg::Texture2D::MIN_FILTER,osg::Texture2D::NEAREST);
theSpectrum34Tex->setFilter(osg::Texture2D::MAG_FILTER,osg::Texture2D::NEAREST);
theSpectrum34Tex->setWrap(osg::Texture2D::WRAP_S,osg::Texture2D::REPEAT);
theSpectrum34Tex->setWrap(osg::Texture2D::WRAP_T,osg::Texture2D::REPEAT);
//createOutputTextures
theSlopeVarianceTex = new osg::Texture3D;
theSlopeVarianceTex->setTextureSize(N_SLOPE_VARIANCE, N_SLOPE_VARIANCE,N_SLOPE_VARIANCE);
theSlopeVarianceTex->setInternalFormat(GL_LUMINANCE16F_ARB);
theSlopeVarianceTex->setSourceFormat(GL_LUMINANCE_ALPHA);
theSlopeVarianceTex->setSourceType(GL_FLOAT);
theSlopeVarianceTex->setFilter(osg::Texture2D::MIN_FILTER,osg::Texture2D: INEAR);
theSlopeVarianceTex->setFilter(osg::Texture2D::MAG_FILTER,osg::Texture2D::LINEAR);
theSlopeVarianceTex->setWrap(osg::Texture3D::WRAP_S,osg::Texture3D::CLAMP_TO_EDGE);
theSlopeVarianceTex->setWrap(osg::Texture3D::WRAP_T,osg::Texture3D::CLAMP_TO_EDGE);
theSlopeVarianceTex->setWrap(osg::Texture3D::WRAP_R,osg::Texture3D::CLAMP_TO_EDGE);
}
void WavesSpectrumPass::build()
{
generateWavesSpectrumTex();
generateSlopeVarianceTex();
}
/************************************************************************/
/* 生成波浪谱纹理 WAVES SPECTRUM GENERATION */
/************************************************************************/
float WavesSpectrumPass::sqr(float x)
{
return x * x;
}
float WavesSpectrumPass: mega(float k)
{
return sqrt(9.81 * k * (1.0 + sqr(k / km))); // Eq 24
}
//波浪谱算法
// 1/kx and 1/ky in meters
float WavesSpectrumPass::spectrum(float kx, float ky, bool omnispectrum)
{
float U10 = theWindSpeed;
float Omega = theInverseWaveAge;
// phase speed
float k = sqrt(kx * kx + ky * ky);
float c = omega(k) / k;
// spectral peak
float kp = 9.81 * sqr(Omega / U10); // after Eq 3
float cp = omega(kp) / kp;
// friction velocity
float z0 = 3.7e-5 * sqr(U10) / 9.81 * pow(U10 / cp, 0.9f); // Eq 66
float u_star = 0.41 * U10 / log(10.0 / z0); // Eq 60
float Lpm = exp(- 5.0 / 4.0 * sqr(kp / k)); // after Eq 3
float gamma = Omega < 1.0 ? 1.7 : 1.7 + 6.0 * log(Omega); // after Eq 3 // log10 or log??
float sigma = 0.08 * (1.0 + 4.0 / pow(Omega, 3.0f)); // after Eq 3
float Gamma = exp(-1.0 / (2.0 * sqr(sigma)) * sqr(sqrt(k / kp) - 1.0));
float Jp = pow(gamma, Gamma); // Eq 3
float Fp = Lpm * Jp * exp(- Omega / sqrt(10.0) * (sqrt(k / kp) - 1.0)); // Eq 32
float alphap = 0.006 * sqrt(Omega); // Eq 34
float Bl = 0.5 * alphap * cp / c * Fp; // Eq 31
float alpham = 0.01 * (u_star < cm ? 1.0 + log(u_star / cm) : 1.0 + 3.0 * log(u_star / cm)); // Eq 44
float Fm = exp(-0.25 * sqr(k / km - 1.0)); // Eq 41
float Bh = 0.5 * alpham * cm / c * Fm * Lpm; // Eq 40 (fixed)
if (omnispectrum) {
return theAmplitude * (Bl + Bh) / (k * sqr(k)); // Eq 30
}
float a0 = log(2.0) / 4.0; float ap = 4.0; float am = 0.13 * u_star / cm; // Eq 59
float Delta = tanh(a0 + ap * pow(c / cp, 2.5f) + am * pow(cm / c, 2.5f)); // Eq 57
float phi = atan2(ky, kx);
if (kx < 0.0) {
return 0.0;
} else {
Bl *= 2.0;
Bh *= 2.0;
}
return theAmplitude * (Bl + Bh) * (1.0 + Delta * cos(2.0 * phi)) / (2.0 * M_PI * sqr(sqr(k))); // Eq 67
}
long WavesSpectrumPass::lrandom(long *seed)
{
*seed = (*seed * 1103515245 + 12345) & 0x7FFFFFFF;
return *seed;
}
float WavesSpectrumPass::frandom(long *seed)
{
long r = lrandom(seed) >> (31 - 24);
return r / (float)(1 << 24);
}
inline float WavesSpectrumPass::grandom(float mean, float stdDeviation, long *seed)
{
float x1, x2, w, y1;
static float y2;
static int use_last = 0;
if (use_last) {
y1 = y2;
use_last = 0;
} else {
do {
x1 = 2.0f * frandom(seed) - 1.0f;
x2 = 2.0f * frandom(seed) - 1.0f;
w = x1 * x1 + x2 * x2;
} while (w >= 1.0f);
w = sqrt((-2.0f * log(w)) / w);
y1 = x1 * w;
y2 = x2 * w;
use_last = 1;
}
return mean + y1 * stdDeviation;
}
//取得波浪谱纹理数据
void WavesSpectrumPass::getSpectrumSample(int i, int j, float lengthScale, float kMin, float *result)
{
static long seed = 1234;
float dk = 2.0 * M_PI / lengthScale;
float kx = i * dk;
float ky = j * dk;
if (abs(kx) < kMin && abs(ky) < kMin) {
result[0] = 0.0;
result[1] = 0.0;
} else {
float S = spectrum(kx, ky);
float h = sqrt(S / 2.0) * dk;
float phi = frandom(&seed) * 2.0 * M_PI;
result[0] = h * cos(phi);
result[1] = h * sin(phi);
}
}
// generates the waves spectrum 产生波浪谱纹理12和波浪谱纹理34
void WavesSpectrumPass::generateWavesSpectrumTex()
{
if (theSpectrum12 != NULL) {
delete[] theSpectrum12;
delete[] theSpectrum34;
}
theSpectrum12 = new float[theTextureSize * theTextureSize * 4];
theSpectrum34 = new float[theTextureSize * theTextureSize * 4];
for (int y = 0; y < theTextureSize; ++y) {
for (int x = 0; x < theTextureSize; ++x) {
int offset = 4 * (x + y * theTextureSize);
int i = x >= theTextureSize / 2 ? x - theTextureSize : x;
int j = y >= theTextureSize / 2 ? y - theTextureSize : y;
getSpectrumSample(i, j, theGRID1_SIZE, M_PI / theGRID1_SIZE, theSpectrum12 + offset);
getSpectrumSample(i, j, theGRID2_SIZE, M_PI * theTextureSize / theGRID1_SIZE, theSpectrum12 + offset + 2);
getSpectrumSample(i, j, theGRID3_SIZE, M_PI * theTextureSize / theGRID2_SIZE, theSpectrum34 + offset);
getSpectrumSample(i, j, theGRID4_SIZE, M_PI * theTextureSize / theGRID3_SIZE, theSpectrum34 + offset + 2);
}
}
osg::ref_ptr<osg::Image> spectrum12Image = new osg::Image;
spectrum12Image->setImage(theTextureSize, theTextureSize, 1, GL_RGBA16F_ARB, GL_RGBA, GL_FLOAT, (unsigned char*)theSpectrum12, osg::Image::USE_NEW_DELETE);
theSpectrum12Tex->setImage(spectrum12Image);
osg::ref_ptr<osg::Image> spectrum34Image = new osg::Image;
spectrum34Image->setImage(theTextureSize, theTextureSize, 1, GL_RGBA16F_ARB, GL_RGBA, GL_FLOAT, (unsigned char*)theSpectrum34, osg::Image::USE_NEW_DELETE);
theSpectrum34Tex->setImage(spectrum34Image);
}
/************************************************************************/
/* 生成海面斜率方差纹理 */
/************************************************************************/
float WavesSpectrumPass::getSlopeVariance(float kx, float ky, float *spectrumSample)
{
float kSquare = kx * kx + ky * ky;
float real = spectrumSample[0];
float img = spectrumSample[1];
float hSquare = real * real + img * img;
return kSquare * hSquare * 2.0;
}
// precomputes filtered slope variances in a 3d texture, based on the wave spectrum 基于波浪谱计算生成海面斜率方差纹理
void WavesSpectrumPass::generateSlopeVarianceTex()
{
// slope variance due to all waves, by integrating over the full spectrum
float theoreticSlopeVariance = 0.0f;
float k = 5e-3;
while (k < 1e3) {
float nextK = k * 1.001;
theoreticSlopeVariance += k * k * spectrum(k, 0, true) * (nextK - k);
k = nextK;
}
// slope variance due to waves, by integrating over the spectrum part
// that is covered by the four nested grids. This can give a smaller result
// than the theoretic total slope variance, because the higher frequencies
// may not be covered by the four nested grid. Hence the difference between
// the two is added as a "delta" slope variance in the "variances" shader,
// to be sure not to lose the variance due to missing wave frequencies in
// the four nested grids
float totalSlopeVariance = 0.0;
for (int y = 0; y < theTextureSize; ++y) {
for (int x = 0; x < theTextureSize; ++x) {
int offset = 4 * (x + y * theTextureSize);
float i = 2.0 * M_PI * (x >= theTextureSize / 2 ? x - theTextureSize : x);
float j = 2.0 * M_PI * (y >= theTextureSize / 2 ? y - theTextureSize : y);
totalSlopeVariance += getSlopeVariance(i / theTextureSize, j / theTextureSize, theSpectrum12 + offset);
totalSlopeVariance += getSlopeVariance(i / theTextureSize, j / theTextureSize, theSpectrum12 + offset + 2);
totalSlopeVariance += getSlopeVariance(i / theTextureSize, j / theTextureSize, theSpectrum34 + offset);
totalSlopeVariance += getSlopeVariance(i / theTextureSize, j / theTextureSize, theSpectrum34 + offset + 2);
}
}
float slopeVarianceDelta=0.5 * (theoreticSlopeVariance - totalSlopeVariance);
//三维纹理按Z轴逐个摄像头独自渲染
for (int layer = 0; layer < N_SLOPE_VARIANCE; layer++)
{
//创建FBO摄像头
osg::ref_ptr<osg::Camera> rttCamera = new osg::Camera;
// clearing
rttCamera->setClearColor(osg::Vec4(0.0f,0.0f,0.0f,0.0f));
rttCamera->setClearMask(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
// projection and view
rttCamera->setProjectionMatrix(osg::Matrix::ortho2D(0,N_SLOPE_VARIANCE,0,N_SLOPE_VARIANCE));
rttCamera->setReferenceFrame(osg::Transform::ABSOLUTE_RF);
rttCamera->setViewMatrix(osg::Matrix::identity());
// viewport
rttCamera->setViewport(0, 0, N_SLOPE_VARIANCE, N_SLOPE_VARIANCE);
rttCamera->setRenderOrder(osg::Camera: RE_RENDER);
rttCamera->setRenderTargetImplementation(osg::Camera::FRAME_BUFFER_OBJECT);
//相当于//glFramebufferTexture3DEXT(target,attachment,textarget,texture,level,zoffset)
rttCamera->attach(osg::Camera::BufferComponent(osg::Camera::COLOR_BUFFER0), theSlopeVarianceTex.get(),0,layer,false);
//创建几何体
osg::ref_ptr<osg::Geode> box=createQuad();
osg::ref_ptr<osg::StateSet> stateset = box->getOrCreateStateSet();
stateset->setAttributeAndModes(theProgram.get(), osg::StateAttribute::ON | osg::StateAttribute::OVERRIDE );
stateset->setMode(GL_LIGHTING,osg::StateAttribute::OFF);
stateset->setTextureAttributeAndModes(0, theSpectrum12Tex.get(), osg::StateAttribute::ON);
stateset->setTextureAttributeAndModes(1, theSpectrum34Tex.get(), osg::StateAttribute::ON);
float nSlopeVariance = N_SLOPE_VARIANCE;
stateset->addUniform(new osg::Uniform("N_SLOPE_VARIANCE", nSlopeVariance));
stateset->addUniform(new osg::Uniform("spectrum_1_2_Sampler", 0));
stateset->addUniform(new osg::Uniform("spectrum_3_4_Sampler", 1));
stateset->addUniform(new osg::Uniform("FFT_SIZE", theTextureSize));
stateset->addUniform(new osg::Uniform("GRID_SIZES", osg::Vec4(theGRID1_SIZE,theGRID2_SIZE,theGRID3_SIZE,theGRID4_SIZE)));
stateset->addUniform(new osg::Uniform("slopeVarianceDelta", slopeVarianceDelta));
float fC = layer;
stateset->addUniform(new osg::Uniform("c", fC));
rttCamera->addChild(box.get());
theRootGroup->addChild(rttCamera.get());
}
}
void WavesSpectrumPass::createShader()
{
static const char* variances_vert = {
"varying vec2 uv; \n"
" \n"
"void main() { \n"
" uv = gl_Vertex.zw; \n"
" gl_Position = vec4(gl_Vertex.xy, 0.0, 1.0); \n"
"} \n"
};
static const char* variances_frag = {
"#define M_PI 3.14159265 \n"
"uniform float N_SLOPE_VARIANCE; \n"
" \n"
"uniform sampler2D spectrum_1_2_Sampler; \n"
"uniform sampler2D spectrum_3_4_Sampler; \n"
"uniform int FFT_SIZE; \n"
" \n"
"uniform vec4 GRID_SIZES; \n"
"uniform float slopeVarianceDelta; \n"
" \n"
"uniform float c; \n"
" \n"
" \n"
"varying vec2 uv; \n"
" \n"
"vec2 getSlopeVariances(vec2 k, float A, float B, float C, vec2 spectrumSample) { \n"
" float w = 1.0 - exp(A * k.x * k.x + B * k.x * k.y + C * k.y * k.y); \n"
" vec2 kw = k * w; \n"
" return kw * kw * dot(spectrumSample, spectrumSample) * 2.0; \n"
"} \n"
" \n"
"void main() { \n"
" const float SCALE = 10.0; \n"
" float a = floor(uv.x * N_SLOPE_VARIANCE); \n"
" float b = floor(uv.y * N_SLOPE_VARIANCE); \n"
" float A = pow(a / (N_SLOPE_VARIANCE - 1.0), 4.0) * SCALE; \n"
" float C = pow(c / (N_SLOPE_VARIANCE - 1.0), 4.0) * SCALE; \n"
" float B = (2.0 * b / (N_SLOPE_VARIANCE - 1.0) - 1.0) * sqrt(A * C); \n"
" A = -0.5 * A; \n"
" B = - B; \n"
" C = -0.5 * C; \n"
" \n"
" vec2 slopeVariances = vec2(slopeVarianceDelta); \n"
" for (int y = 0; y < FFT_SIZE; ++y) { \n"
" for (int x = 0; x < FFT_SIZE; ++x) { \n"
" int i = x >= FFT_SIZE / 2 ? x - FFT_SIZE : x; \n"
" int j = y >= FFT_SIZE / 2 ? y - FFT_SIZE : y; \n"
" vec2 k = 2.0 * M_PI * vec2(i, j); \n"
" \n"
" vec4 spectrum12 = texture2D(spectrum_1_2_Sampler, vec2(float(x) + 0.5, float(y) + 0.5) / float(FFT_SIZE)); \n"
" vec4 spectrum34 = texture2D(spectrum_3_4_Sampler, vec2(float(x) + 0.5, float(y) + 0.5) / float(FFT_SIZE)); \n"
" \n"
" slopeVariances += getSlopeVariances(k / GRID_SIZES.x, A, B, C, spectrum12.xy); \n"
" slopeVariances += getSlopeVariances(k / GRID_SIZES.y, A, B, C, spectrum12.zw); \n"
" slopeVariances += getSlopeVariances(k / GRID_SIZES.z, A, B, C, spectrum34.xy); \n"
" slopeVariances += getSlopeVariances(k / GRID_SIZES.w, A, B, C, spectrum34.zw); \n"
" } \n"
" } \n"
" gl_FragColor = slopeVariances.xxxy; \n"
"} \n"
};
theProgram = ShaderManager::instance().createProgram("variances", variances_vert, variances_frag,"", false );
} |
|
IP卡
狗仔卡
发表于 2011-4-27 09:44:26
提升卡
置顶卡
沉默卡
喧嚣卡
变色卡
显身卡
发表于 2011-4-28 08:17:31
