cesium后处理的中浓度积分高度雾(下)-优化后

原创已于 2026-04-07 21:37:46 修改 · 1.8k 阅读

42 ·

本内容遵循CC 4.0 BY-SA版权协议

GEO检测

于 2024-11-03 23:10:57 首次发布

该文章已生成可运行项目，

一、回顾

上一篇我们将了如何在后处理中获取世界坐标还有高程，还了解了两个简单的线性雾，从上一篇的线性高度雾我们可以发现，一旦相机不是大角度俯视，那么雾的表现就会失真，如果相机进入雾的范围则更假，因为雾气只是简单的贴在物体和地形上，并不会对整个视野造成影响。

点击去往视频教程

二、浓度积分高度雾公式推导

1.为什么用积分

下图中，雾气浓度随着高程增加不断变低(先不用管怎么变)，对于同一个点，我们从三个不同的角度观看，离目标点的距离都是相等的，那么最终看到的目标颜色会相同吗？或者说这个点最终的雾值会一样吗？其实我们都知道它们应该是不一样的，并且 F1 < F2 < F3 。情况3下视线被雾气遮挡的最多，情况2则最少。

那么怎么模拟这种情况呢？这就需要用到积分了。

现在我们假设雾气随高度变化的函数为 $f = 1 - \frac{1}{100} \cdot h$ ，这里的f指的是雾气浓度fog density。从这个函数我们可以知道在高程为0时雾气浓度为1，在高程为100时，雾气浓度为0，图像大概下面所示，随着高程的增加雾浓度不断降低。

有了浓度随高程变化的公式，我们可以很轻松的算出不同高程的雾气浓度，例如高程50米雾气浓度为0.5。那么假设，现在我们从100米垂直看向地面0米的地方，那么这段距离的雾总量是多少呢，我们知道其实就是这个函数的线与坐标轴围成的三角形的面积，我们当然可以用 0.5 * 100 * 1.0 = 50 来得到最终的结果。

但是我们还有更科学的计算方式，那就是积分，雾气浓度随着高度变化的函数我们有了，积分上限为100，下限为0积分，公式为：

$\int_{0}^{100}f(h)d_{h} = \int_{0}^{100}(1-\frac{1}{100}\cdot h)d_{h} = h-\frac{1}{100\times 2}\cdot h^{2}|^{100}_0$

最终结果当然也是50，但是积分的存在可以帮我们计算曲线函数围成的面积。

2.公式推导

我们知道看向一个物体，最终物体被雾气遮挡的程度取决于我们视线上雾气的总和，如果我们能知道随着视线距离变化的雾气浓度就好了，那么直接对视线整个距离积分就可以得到雾气总和。我们知道高度雾的雾气浓度是随着高度变化而变化的，如果我们知道高度与视线距离的关系我们就其实也知道雾气浓度随距离变换的公式了。

例如我们知道雾气浓度与高度的公式如下，（H是常量，举个线性关系的例子）：

$f(h)=1-\frac{h}{H}$

如果再有高度与经过的雾气距离的关系，就能得到雾气浓度与经过的雾气距离的关系：

$f(s)=1-\frac{h(s)}{H}$

而高度与雾气距离的关系是可以固定的一个函数，如下图

h0为积分起始高度，h1为积分结束高度，则有高度与积分距离s的关系h(s)。

$h(s) = h_{0} + \frac{h_{1} - h_{0}}{L}\cdot s$

则有雾气浓度与距离的关系：

$f(s)=1-\frac{h(s)}{H} = 1-\frac{h_0}{H} + \frac{h_1-h_0}{HL}\cdot s$

而有了雾气浓度与距离的函数f(s)，接下来对这个函数求导不就可以得到这段视线上雾气的总量F：
如果雾气浓度与高度的关系发生了变化，我们只需要改变f(h)函数即可，其他的步骤完全一致。

$F = \int_{0}^{L}f(h)d_h = \int_{0}^{L}1-\frac{h}{H}d_h =\int_0^L\left(1-\frac{h(s)}{H}\right)ds$

$= \int_0^L\left(1-\frac{h(s)}{H}\right)ds =\int_0^L(1 - \frac{h_0}{H} + \frac{h_1-h_0}{LH}\cdot s)ds$

$=\int_{0}^{L}1-\frac{h_0}{H}ds - \int_{0}^{L}\frac{h_1-h_0}{LH}ds$

$=(1 - \frac{h_0}{H})\cdot L - \frac{(h_1-h_0)}{2H} \cdot L$

$=(1 - \frac{(h_1+h_0)}{2H}) \cdot L$

三、shader代码

现在我们知道了求取视线上雾总量的公式，接下来就是将公式转为shader代码了。

观察最终的公式：

$F=(1 - \frac{(h_1+h_0)}{2H}) \cdot L$

我们需要知道h0，h1，L，H四个的值，其中H为雾气浓度与高度关系函数中的常量，算已知量。

而h1雾气积分的起始位置，h0，L则需要根据片元位置、相机位置、雾气高度等参数请求。

下面代码中fogStartHeight为h0，fogEndHeight为h1，getFogLayerDistance函数的返回值fogDistance为L.


float getFogLayerDistance(float startHeight, float endHeight, float totalDistance, float fogMaxHeight, out float fogStartHeight, out float fogEndHeight) {
	fogStartHeight = 0.0;
	fogEndHeight = 0.0;

	if(totalDistance <= 0.0) {
		return 0.0;
	}

	// 判断有效视线路径的起点和终点是否落在雾层内
	bool startInside = startHeight <= fogMaxHeight;
	bool endInside = endHeight <= fogMaxHeight;

	// 整段路径都在雾层内，直接返回整段长度
	if(startInside && endInside) {
		fogStartHeight = startHeight;
		fogEndHeight = endHeight;
		return totalDistance;
	}

	// 整段路径都在雾层外，不产生高度雾
	if(!startInside && !endInside) {
		return 0.0;
	}

	float heightDelta = endHeight - startHeight;
	if(abs(heightDelta) <= 0.0001) {
		return 0.0;
	}

	// 只有一部分穿过雾层时，求出与雾层上边界的交点比例
	float t = clamp((fogMaxHeight - startHeight) / heightDelta, 0.0, 1.0);
	if(startInside) {
		fogStartHeight = startHeight;
		fogEndHeight = fogMaxHeight;
		return totalDistance * t;
	}

	fogStartHeight = fogMaxHeight;
	fogEndHeight = endHeight;
	return totalDistance * (1.0 - t);
}

float getFogDistanceFromCamera(vec3 positionToCamera, float cameraHeight, float pixelHeight, float fogMaxHeight, out float fogStartHeight, out float fogEndHeight) {
	fogStartHeight = 0.0;
	fogEndHeight = 0.0;

	// 先得到相机到当前像素的总视线长度
	float totalDistance = length(positionToCamera);
	if(totalDistance <= u_fogStartDis) {
		return 0.0;
	}

	// 再裁剪出这段有效路径中真正位于雾层内的部分
	return getFogLayerDistance(startHeight, pixelHeight, totalDistance, fogMaxHeight, fogStartHeight, fogEndHeight);
}

获取了h0，h1，L这三个参数，那么根据公式带入雾浓度与高度的函数f(h)，就可以了。

下面是f(h) = 1 - h / H 雾浓度关系的高度雾。

// 线性浓度积分高度雾
// f(h) = 1.0 - h / fogMaxHeight
float linearHeightFog(vec3 positionToCamera, float cameraHeight, float pixelHeight, float fogMaxHeight) {
	float globalDensity = u_globalDensity / 10.0;
	float safeFogMaxHeight = max(fogMaxHeight, 0.0001);
	float fogStartHeight;
	float fogEndHeight;
	float fogDistance = getFogDistanceFromCamera(positionToCamera, cameraHeight, pixelHeight, safeFogMaxHeight, fogStartHeight, fogEndHeight);
	if(fogDistance <= 0.0) {
		return 0.0;
	}

	float averageDensity = max(0.0, 1.0 - 0.5 * (fogStartHeight + fogEndHeight) / safeFogMaxHeight);
	float fog = globalDensity * fogDistance * averageDensity;

	fog = mix(0.0, 1.0, fog / (fog + 1.0));
	return fog;
}

初始化viewer，复制上面的代码，加载地形。没问题的话以及可以显示出下面的结果了。有问题的话记得检查一下是否开启了深度测试。

四、添加雾气产生距离

1.错误做法

简单的将一定距离内的雾总量变为0是行不通的，这样会导致分界线特别明显割裂。比如下面的代码，如果点到相机的距离小于雾气产生距离，那么我们就让雾总量为0。

if(length(positionToCamera) < u_fogStartDis){
	fog = 0.0;
}

结果就是下图，非常割裂。

2.正确做法

只需要修改之前封装的getFogDistanceFromCamera函数，扣除起雾距离后再计算h0，h1，L即可

float getFogDistanceFromCamera(vec3 positionToCamera, float cameraHeight, float pixelHeight, float fogMaxHeight, out float fogStartHeight, out float fogEndHeight) {
	fogStartHeight = 0.0;
	fogEndHeight = 0.0;

	// 先得到相机到当前像素的总视线长度
	float totalDistance = length(positionToCamera);
	if(totalDistance <= u_fogStartDis) {
		return 0.0;
	}

	// 扣除起雾距离后，计算雾真正开始作用位置的高度
	float startRatio = clamp(u_fogStartDis / totalDistance, 0.0, 1.0);
	float startHeight = mix(cameraHeight, pixelHeight, startRatio);
	float endHeight = pixelHeight;
	float effectiveDistance = totalDistance - u_fogStartDis;

	// 再裁剪出这段有效路径中真正位于雾层内的部分
	return getFogLayerDistance(startHeight, endHeight, effectiveDistance, fogMaxHeight, fogStartHeight, fogEndHeight);
}

结果，可以看到现在就有过渡效果了，不再那么割裂。

五、引申其他浓度关系的雾

其实浓度积分高度雾的公式非常的简单，就如下面所示。

$F = \int_{0}^{L}f(h)d_h$

如果想实现指数浓度积分高度雾也非常的简单，只需要将雾浓度随高度变化的函数换成指数函数就可以了。

例如换成 $f = a\cdot e^{-bh}$ 。a为全局浓度，b为浓度衰减系数。然后最终代码如下所示：

由于heightDelta在指数关系中会成为分母，所以需要分为2中情况，视线相对水平的详情下我们任务视线上的雾浓度不变，可以回退为f(h) = X，其中X是一个常量，为公式f(h) = 1.0 - h / fogMaxHeight带入h求得，高度雾函数如下：

// 指数浓度积分高度雾
// f(h) = exp(-u_heightFalloff * h)
float exponentialHeightFog(vec3 positionToCamera, float cameraHeight, float pixelHeight, float fogMaxHeight) {
	float globalDensity = u_globalDensity / 10.0;
	float safeHeightFalloff = max(u_heightFalloff, 0.0001);
	float fogStartHeight;
	float fogEndHeight;
	float fogDistance = getFogDistanceFromCamera(positionToCamera, cameraHeight, pixelHeight, fogMaxHeight, fogStartHeight, fogEndHeight);
	if(fogDistance <= 0.0) {
		return 0.0;
	}

	float heightDelta = fogEndHeight - fogStartHeight;
	float fog;
	if(abs(heightDelta) <= 0.001) {
		fog = globalDensity * exp(-safeHeightFalloff * fogStartHeight) * fogDistance;
	} else {
		fog = globalDensity * fogDistance * (exp(-safeHeightFalloff * fogStartHeight) - exp(-safeHeightFalloff * fogEndHeight)) / (safeHeightFalloff * heightDelta);
	}

	fog = mix(0.0, 1.0, fog / (fog + 1.0));
	return fog;
}

六、总结

通过对视线上的雾气浓度进行积分得到视线上的雾气总量，这比较真实的反映了高度雾的效果，既有高度上的变化也有距离上的变化。

原本文章的代码如果片元到相机的矢量与地面接近水平的情况下，该片元上的雾浓度计算会出问题，当前文章已对该问题进行修复，并且仿照ceisum的函数使用了getEllipsoidHeight来计算高度。

总体vue代码如下。initCesiumViewer为本人封装的代码，需要替换为初始viewer的代码。

<template>
	<div id="cesium-container">
		<el-card class="panel">
			<div class="title">
				<span>{{ fogType === 'lin' ? '线性' : '指数' }}雾气参数</span>
				<el-switch v-model="fogType" active-value="exp" inactive-value="lin" active-text="指数雾" inactive-text="线性雾" />
			</div>
			<div class="item">
				<span>雾颜色：</span>
				<PickColors v-model:value="u_fogColor" format="rgb" />
			</div>
			<div class="item">
				<span>最大高度：</span><el-slider v-model="u_fogHeight" :min="100" :max="3000" />
			</div>
			<div class="item">
				<span>全局浓度：</span><el-slider v-model="u_globalDensity" :min="0.1" :max="2.0" :step="0.1" />
			</div>
			<div class="item" v-if="fogType == 'exp'">
				<span>浓度衰减系数：</span><el-slider v-model="u_heightFalloff" :min="0.0001" :max="0.01" :step="0.0001" />
			</div>
			<div class="item">
				<span>雾气起始距离：</span><el-slider v-model="u_fogStartDis" :min="100" :max="2000" />
			</div>
		</el-card>
	</div>
</template>

<script setup lang='ts'>
import {
	Cartesian3,
	PostProcessStage,
	Color,
	createWorldTerrainAsync,
} from 'cesium'
import PickColors from 'vue-pick-colors'

const fogType = ref<'lin' | 'exp'>('lin')
const u_fogColor = ref('rgb(204, 209, 214)')
const u_fogHeight = ref(1000)
const u_globalDensity = ref(0.6)
const u_heightFalloff = ref(0.002)
const u_fogStartDis = ref(500)

onMounted(async () => {
	const viewer = initCesiumViewer({ baseLayer: useAMapLayer('img') })
	viewer.scene.globe.depthTestAgainstTerrain = true

	viewer.camera.setView({
		destination: Cartesian3.fromDegrees(114.1, 25.1, 4000),
		orientation: {
			pitch: -0.6
		}
	})
	const handler = ElMessage.info({ message: '地形加载中。。。', duration: 0 })
	viewer.scene.globe.terrainProvider = await createWorldTerrainAsync()
	handler.close()
	ElMessage.success({ message: '地形加载完毕' })

	// 后处理可以使用cesium自带的，通过PostProcessStageLibrary中的静态函数拿取
	// const postProcessStages = PostProcessStageLibrary.createBlackAndWhiteStage()
	// 然后通过viewer.postProcessStages.add添加到地图中
	// viewer.postProcessStages.add(postProcessStages)

	// 或者通过直接new PostProcessStage()自定义一个后处理
	const fs = `
    uniform sampler2D colorTexture;  // 颜色纹理
		uniform sampler2D depthTexture;  // 深度纹理
		in vec2 v_textureCoordinates;  // 纹理坐标
		uniform float u_earthRadiusOnCamera;
		uniform float u_cameraHeight;
		uniform float u_fogHeight;
		uniform vec3 u_fogColor;
		uniform float u_globalDensity;
		uniform float u_heightFalloff;
		uniform float u_fogStartDis;
		uniform float isLinearFog;

		// 通过深度纹理与纹理坐标得到世界坐标
		vec4 getWorldCoordinate(sampler2D depthTexture, vec2 texCoords){
      float depthOrLogDepth = czm_unpackDepth(texture(depthTexture, texCoords));
      vec4 eyeCoordinate = czm_windowToEyeCoordinates(gl_FragCoord.xy, depthOrLogDepth);
			eyeCoordinate = eyeCoordinate / eyeCoordinate.w;
      vec4 worldCoordinate = czm_inverseView * eyeCoordinate;
			worldCoordinate = worldCoordinate / worldCoordinate.w;

      return worldCoordinate;
    }

		// 计算粗略的高程，依赖js传递的相机位置处的地球高程u_earthRadiusOnCamera。好处是计算量非常低
		float getRoughHeight(vec4 worldCoordinate){
			float disToCenter = length(vec3(worldCoordinate));
			return disToCenter - u_earthRadiusOnCamera;
		}

		// 计算基于WGS84椭球面的真实高程。相较粗略高程，在远距离和接近地平线时更稳定
		float getEllipsoidHeight(vec4 worldCoordinate){
			vec3 positionWC = vec3(worldCoordinate);
			vec3 oneOverRadii = czm_ellipsoidInverseRadii;
			vec3 oneOverRadiiSquared = oneOverRadii * oneOverRadii;
			// 先将世界坐标缩放到椭球参数空间，得到初始径向投影
			vec3 scaledPositionSquared = positionWC * positionWC * oneOverRadiiSquared;
			float squaredNorm = scaledPositionSquared.x + scaledPositionSquared.y + scaledPositionSquared.z;
			float ratio = sqrt(1.0 / max(squaredNorm, 0.0000001));
			vec3 intersection = positionWC * ratio;

			if(squaredNorm < 0.0000001) {
				return 0.0;
			}

			vec3 gradient = intersection * oneOverRadiiSquared * 2.0;
			float lambda = ((1.0 - ratio) * length(positionWC)) / max(0.5 * length(gradient), 0.0000001);
			float correction = 0.0;
			vec3 multipliers = vec3(1.0);

			// 使用牛顿迭代，把初始投影修正到真正的椭球面上
			for(int i = 0; i < 8; i++) {
				lambda -= correction;
				multipliers = 1.0 / (1.0 + lambda * oneOverRadiiSquared);
				vec3 multipliers2 = multipliers * multipliers;
				vec3 multipliers3 = multipliers2 * multipliers;
				float func = dot(scaledPositionSquared, multipliers2) - 1.0;
				if(abs(func) <= 0.0000001) {
					break;
				}
				float denominator = dot(scaledPositionSquared, multipliers3 * oneOverRadiiSquared);
				float derivative = -2.0 * max(denominator, 0.0000001);
				correction = func / derivative;
			}

			// 回推出椭球面交点，再由交点到原始点的距离得到高程
			vec3 surfacePosition = positionWC * multipliers;
			vec3 heightVector = positionWC - surfacePosition;
			return sign(dot(heightVector, positionWC)) * length(heightVector);
		}

		// 得到a向量在b向量的投影长度，如果同向结果为正，异向结果为复
		float projectVector(vec3 a, vec3 b) {
    	float scale = dot(a, b) / dot(b, b);
			if(abs(scale) <= 0.000001) {
				return 0.0;
			}
			return sign(scale) * length(scale * b);
		}

		float getFogLayerDistance(float startHeight, float endHeight, float totalDistance, float fogMaxHeight, out float fogStartHeight, out float fogEndHeight) {
			fogStartHeight = 0.0;
			fogEndHeight = 0.0;

			if(totalDistance <= 0.0) {
				return 0.0;
			}

			// 判断有效视线路径的起点和终点是否落在雾层内
			bool startInside = startHeight <= fogMaxHeight;
			bool endInside = endHeight <= fogMaxHeight;

			// 整段路径都在雾层内，直接返回整段长度
			if(startInside && endInside) {
				fogStartHeight = startHeight;
				fogEndHeight = endHeight;
				return totalDistance;
			}

			// 整段路径都在雾层外，不产生高度雾
			if(!startInside && !endInside) {
				return 0.0;
			}

			float heightDelta = endHeight - startHeight;
			if(abs(heightDelta) <= 0.0001) {
				return 0.0;
			}

			// 只有一部分穿过雾层时，求出与雾层上边界的交点比例
			float t = clamp((fogMaxHeight - startHeight) / heightDelta, 0.0, 1.0);
			if(startInside) {
				fogStartHeight = startHeight;
				fogEndHeight = fogMaxHeight;
				return totalDistance * t;
			}

			fogStartHeight = fogMaxHeight;
			fogEndHeight = endHeight;
			return totalDistance * (1.0 - t);
		}

		float getFogDistanceFromCamera(vec3 positionToCamera, float cameraHeight, float pixelHeight, float fogMaxHeight, out float fogStartHeight, out float fogEndHeight) {
			fogStartHeight = 0.0;
			fogEndHeight = 0.0;

			// 先得到相机到当前像素的总视线长度
			float totalDistance = length(positionToCamera);
			if(totalDistance <= u_fogStartDis) {
				return 0.0;
			}

			// 扣除起雾距离后，计算雾真正开始作用位置的高度
			float startRatio = clamp(u_fogStartDis / totalDistance, 0.0, 1.0);
			float startHeight = mix(cameraHeight, pixelHeight, startRatio);
			float endHeight = pixelHeight;
			float effectiveDistance = totalDistance - u_fogStartDis;

			// 再裁剪出这段有效路径中真正位于雾层内的部分
			return getFogLayerDistance(startHeight, endHeight, effectiveDistance, fogMaxHeight, fogStartHeight, fogEndHeight);
		}

		// 线性浓度积分高度雾
		// f(h) = 1.0 - h / fogMaxHeight
		float linearHeightFog(vec3 positionToCamera ,float cameraHeight, float pixelHeight, float fogMaxHeight){
			// 密度随高度线性衰减
			float globalDensity = u_globalDensity / 10.0;
			float safeFogMaxHeight = max(fogMaxHeight, 0.0001);
			float fogStartHeight;
			float fogEndHeight;
			float fogDistance = getFogDistanceFromCamera(positionToCamera, cameraHeight, pixelHeight, safeFogMaxHeight, fogStartHeight, fogEndHeight);
			if(fogDistance <= 0.0) {
				return 0.0;
			}

			float averageDensity = max(0.0, 1.0 - 0.5 * (fogStartHeight + fogEndHeight) / safeFogMaxHeight);
			float fog = globalDensity * fogDistance * averageDensity;

			fog = mix(0.0, 1.0, fog / (fog + 1.0));
			return fog;
		}

		// 指数浓度积分高度雾
		// f(h) = exp(-u_heightFalloff * h)
		float exponentialHeightFog(vec3 positionToCamera, float cameraHeight, float pixelHeight, float fogMaxHeight){
			float globalDensity = u_globalDensity / 10.0;
			float safeHeightFalloff = max(u_heightFalloff, 0.0001);
			float fogStartHeight;
			float fogEndHeight;
			float fogDistance = getFogDistanceFromCamera(positionToCamera, cameraHeight, pixelHeight, fogMaxHeight, fogStartHeight, fogEndHeight);
			if(fogDistance <= 0.0) {
				return 0.0;
			}

			float heightDelta = fogEndHeight - fogStartHeight;
			float fog;
			if(abs(heightDelta) <= 0.001) {
				fog = globalDensity * exp(-safeHeightFalloff * fogStartHeight) * fogDistance;
			}
			else {
				fog = globalDensity * fogDistance * (exp(-safeHeightFalloff * fogStartHeight) - exp(-safeHeightFalloff * fogEndHeight)) / (safeHeightFalloff * heightDelta);
			}

			fog = mix(0.0, 1.0, fog / (fog + 1.0));
			return fog;
		}

		void main(void) 
		{
			vec4 color = texture(colorTexture, v_textureCoordinates);
			vec4 positionWC = getWorldCoordinate(depthTexture, v_textureCoordinates);	
			float pixelHeight = getEllipsoidHeight(positionWC);
			vec3 positionToCamera = vec3(vec3(positionWC) - czm_viewerPositionWC);

			float fog;
			if(isLinearFog > 0.5){
				fog = linearHeightFog(positionToCamera, u_cameraHeight, pixelHeight, u_fogHeight);
			}
			else{
				fog = exponentialHeightFog(positionToCamera, u_cameraHeight, pixelHeight, u_fogHeight);
			}
			out_FragColor = mix(color, vec4(u_fogColor, 1.0), fog);
		}`


	const customPostProcessStage = new PostProcessStage({
		fragmentShader: fs,
		uniforms: {
			u_earthRadiusOnCamera: () => Cartesian3.magnitude(viewer.camera.positionWC) - viewer.camera.positionCartographic.height,
			u_cameraHeight: () => viewer.camera.positionCartographic.height,
			u_fogColor: () => Color.fromCssColorString(u_fogColor.value),
			u_fogHeight: () => u_fogHeight.value,
			u_globalDensity: () => u_globalDensity.value,
			u_heightFalloff: () => u_heightFalloff.value,
			u_fogStartDis: () => u_fogStartDis.value,
			isLinearFog: () => fogType.value === 'lin' ? 1 : 0
		}
	})


	viewer.scene.postProcessStages.add(customPostProcessStage)
})


</script>

<style scoped lang='scss'>
.panel {
	position: absolute;
	top: 20px;
	right: 20px;
	width: 340px;
	z-index: 1;

	.title {
		display: flex;
		font-weight: bold;
		margin-bottom: 10px;

		.el-switch {
			margin-left: auto;
		}
	}

	.item {
		display: flex;

		span {
			width: 120px;
		}

		.el-slider {
			width: 180px;
		}
	}
}
</style>

本文章已经生成可运行项目