天文望远镜

你可能看过很多炫酷的天体照片，大型射电望远镜拍的那种，五颜六色的，像这样：

但是你肯定从来没有抬头真的在天空上看见过这样的玩意。是它们太暗、太远了的缘故，所以我们看不到吗？ Well，一半是吧。其实就算把它们拉到几百个太阳系半径这么短的距离上来，（假设地球还没因此而消失）你抬头看到的也不可能有这么炫酷。这是因为绝大部分天体的发光频段都在可见光的范围外。

这是什么意思？想想红外线，绝大多数有一定温度的物体，比如人体、地表，都在源源不断地向外辐射着红外线。和可见光一样，红外线也是电磁波，然而人眼却看不到。再想想 X 光，它和红外线处在可见光频段的另外一头，也是人眼看不到的电磁波。没听说过谁去医院拍片子被闪瞎了狗眼。

看看上面这幅图，可见光在全频段上才占了多小的一部分啊！人类怎么可能指望随便一个天体发的光正好是人眼能看到的呢？完全没道理的。那么，那些炫酷的天体图片是怎么来的？难不成是天文机构在纯纯唬人？

其实不然。为了让地球人能够看到天体的样貌，制作这些图片的家伙们往往会把超出可见光频段的电磁波频谱按照一定的规则转换到可见光的频段里。频段里靠近低频的那一端映射到偏红的颜色，高频的那一端则映射到偏蓝的颜色。这有点像是图像编辑软件里的色相调节功能，但不完全一样。色相调节会让高频率的蓝色越过紫色、粉色，最后到达低频率的红色来，而这种转换会保持频率之间的大小关系。为了方便指代，我们将这种保持大小关系的频率转换称为频谱重映射。

频谱重映射

让我们来讨论一下这种映射的数学形式。为了方便起见，我们不妨用区间来表示频段。低频端对应左端点，高频段对应右端点；具体的数值和频率的单位有关，我们暂且用 Hz。

现在，我想将 $[100\text{kHz},200\text{kHz}]$ 映射到 $[400\text{bHz},700\text{bHz}]$（即可见光频段）上去，如何？一个很自然的想法是，可以用一个线性变换来将源区间完美覆盖到目标区间上。这样的线性变换是唯一的。对于这两个区间来说，这变换是 $x\mapsto 3\cdot{10}^6\cdot x+{10}^{11}$。那么，100kHz 会被映射成 400bHz，即红色；200kHz 会被映射成 700bHz，即蓝色；而处于它俩中间的 150kHz 则会被映射成 550bHz，一种大概绿色的颜色。

看上去很合理，right？

Well… actually not. 考虑将 $[10\text{Hz},11\text{Hz}]$ 映射到 $[1\text{Hz},5\text{Hz}]$ 的线性映射，它是 $x\mapsto 5\cdot(x-10)$。那么对于任意小于 10Hz 的频率，这映射将会把它映射到一个负数；甚至对于正好 10Hz，映射后的结果是 0Hz——一个直流信号！这明显不是我们想要的。

为了补救这个缺陷，我们可以把线性变换的对象换成频率的对数。因为所有频率都是非负数，而我们将一个对数值通过指数函数恢复过来时，得到的也总是一个非负数，这就很合理了。新的变换形式是： $$\ln x\mapsto T(\ln x)=k\ln x+b.$$ 更详细地，若欲将 $[a,b]$ 映射到 $[0,1]$，那么映射为： $$x\mapsto\exp\!\left(\frac{\ln x-\ln a}{\ln b-\ln a}\right),$$ 而它的逆变换（即 $[0,1]$ 映射到 $[a,b]$）为： $$x\mapsto e^{(\ln b-\ln a)x+\ln a}=a\ln\!\left(\frac ba\right)^x.$$

特别要注意的一点是，由于对数值是没有单位和零点的，我们可以直接不用考虑实际频率而强行让可见光区段为 $[0,1]$。这就方便了好多。

Unity Shader

前一阵子我们在开发一个关于星体观测的功能性游戏，我就想能不能把这个玩意用在游戏机制里（很遗憾地，组里没人对这个机制感兴趣）。这样就得有一个 shader 来实现这种颜色变化的效果。总之，先来写个 shader 起手式：

Shader "Custom/SpectrumRemapping" {
	Properties {
		_MainTex("Texture", 2D) = "white" {}
	}
	SubShader {
		Pass {
			CGPROGRAM
			#pragma vertex vertex
			#pragma fragment fragment

			#include "UnityCG.cginc"

			struct VertexIn {
				float4 position : POSITION;
				float2 uvCoord : TEXCOORD0;
			};

			struct FragmentIn
			{
				float2 uvCoord : TEXCOORD0;
				float4 position : SV_POSITION;
			};

			FragmentIn vertex(VertexIn v) {
				FragmentIn o;
				o.position = UnityObjectToClipPos(v.position);
				o.uvCoord = v.uvCoord;
				return o;
			}

			sampler2D _MainTex;

			fixed4 fragment(FragmentIn i) : SV_Target{
				fixed4 sample = tex2D(_MainTex, i.uvCoord);
				return sample;
			}
			ENDCG
		}
	}
}

这里我们主要关注的就是片元着色器 fixed4 fragment(FragmentIn)，因为它决定了像素的颜色。但在实际开始编写转换函数之前，我们得先搞明白一个事:

相较于 RGB 分量，我们更关注颜色的色相，所以我们不妨先将颜色转换到 HSV 空间里讨论。

如图，这是一个色轮。就像前面所说的那样，红色到蓝色之间这部分粉、紫色的区域是单色光所未有的色相，它们不能代表任何频率。因此，我们希望 $[0,1]$ 所指代的颜色正好对应 RGB 红色到 RGB 蓝色这一段区间；当原始颜色在粉、紫色区间时，我们认为它们对应红外或紫外波段的频率，应当在色相上截取到红或蓝，并在明度上淡出。如图：

现在，首先让我们编写 RGB 和 HSV 互相转换的函数：

static const fixed tau = atan(1) * 8;
static const fixed sqrt3 = sqrt(3);
static const fixed halfSqrt3 = sqrt3 * .5f;

fixed3 RGB2HSV(fixed3 rgb) {
	fixed3 xyz
		= rgb.r * fixed3(1, 0, 1)
		+ rgb.g * fixed3(-.5f, halfSqrt3, 1)
		+ rgb.b * fixed3(-.5f, -halfSqrt3, 1);
	fixed h = frac(atan2(xyz.y, xyz.x) / tau + 1);
	fixed s = length(xyz.xy);
	fixed v = xyz.z / 3;
	return fixed3(h, s, v);
}

fixed3 HSV2RGB(fixed3 hsv) {
	hsv.x *= tau;
	fixed3 xyz = fixed3(
		fixed2(cos(hsv.x), sin(hsv.x)) * hsv.y,
		hsv.z * 3
	);
	fixed3 rgb
		= xyz.x * fixed3(2, -1, -1)
		+ xyz.y * fixed3(0, sqrt3, -sqrt3)
		+ xyz.z * fixed3(1, 1, 1);
	rgb /= 3;
	return rgb;
}

这两段转换函数与标准方法有出入。为了图省事我采用了在原理上比较简单的方法，但是间色的明度会只有基色的 $\sqrt3/2$ 倍。

这里的原理很简单，先把三个基色想象成平面上互呈 $120^\circ$ 夹角、Z 分量为 $1/3$ 的基向量。只消以 RGB 分量作为权重把它们加到一起，所成的向量在 XY 平面上的方向即为色相、长度即为饱和度，而 Z 坐标即为明度。

然后是 HSV 到 XSV（X 为频率）的转换函数：

static const fixed hueLimit = 5.f / 6;

// HSV: Hue - Spectrum - Value
// XSV: Relative Frequency - # - #
fixed3 HSV2XSV(fixed3 hsv) {
	if(hsv.x > hueLimit)
		hsv.x -= 1;
	hsv.x /= hueLimit;
	return hsv;
}

fixed3 XSV2HSV(fixed3 xsv) {
	fixed offset = abs(xsv.x - clamp(xsv.x, 0, 1));
	fixed dim = exp(-offset);
	fixed3 hsv = fixed3(
		clamp(xsv.x, 0, 1) * hueLimit,
		xsv.y,
		xsv.z * dim
	);
	return hsv;
}

这里 hueLimit 就是指的粉、紫色区域在色相环上的方位（正好是 5/6 圈的地方）。在转 XSV 之前，我们要先把红蓝区间归一化。

在从 XSV 转回 HSV 时，由于我们希望可见区段附近的频率渐隐出去，这里用到了一个小技巧： $\left|x-\mathrm{clamp}(x,S)\right|$ 可以计算实数 $x$ 到区间 $S$ 的距离（事实上，clamp 得出的结果是 $S$ 里距 $x$ 最近的点）。然后我们将距离送进 $-\exp$，就得到了一个相对合理的衰减系数。将它乘到明度分量上，就实现了衰减的效果（别忘了截断色相分量）。

最后是线性变换的部分，直接抄公式：

float IntervalTransform(float x, float2 source, float2 spectating) {
	x = (source.y - source.x) * x + source.x;
	x = (x - spectating.x) / (spectating.y - spectating.x);
	return x;
}

把所有步骤拼在一起，加上一个四维向量的属性用来调节源频段和目标频段（XY 分量是源频段对数值，ZW 分量是目标频段对数值）：

Properties {
	_Frequencies("Frequencies (log)", Vector) = (0, 1, 0, 1)
}

// ...

fixed3 RGB2XSV(fixed3 rgb) {
	return HSV2XSV(RGB2HSV(rgb));
}

fixed3 XSV2RGB(fixed3 xsv) {
	return HSV2RGB(XSV2HSV(xsv));
}

float4 _Frequencies;

fixed4 fragment(FragmentIn i) : SV_Target{
	fixed4 sample = tex2D(_MainTex, i.uvCoord);
	fixed3 xsv = RGB2XSV(sample.rgb);
	xsv.x = IntervalTransform(xsv.x, _Frequencies.xy, _Frequencies.zw);
	fixed3 _rgb = XSV2RGB(xsv);
	return fixed4(_rgb, sample.a);
}

最终效果如下：

我把 shader 部分的代码放在了 gist 上，有需自取（国内访问需翻墙）。