PBRT阅读：第八章 胶片和图象管线 第8.4-8.5节

http://www.opengpu.org/forum.php?mod=viewthread&tid=5306

8.4 感知上的问题和色调映射

在计算机图形学的早期年代，着色模型总是返回位于0到1之间的颜色值，跟实际的物理量没有什么联系。因此，像素值也是位于这个范围，只要对像素值进行比例变换，就可以将图像直接显示到带有RGB(范围0 ~ 255)帧缓存的CRT显示器上。在真实世界中，经常存在辐射度值等级大约从0.01到1000（即从最亮部分到最暗部分的变化有5个等级）的场景。值得注意的是，人类视觉系统(HVS)却可以很好地处理这种亮度的极端变化，因为人的眼睛对局部的亮度对比要比绝对的亮度更敏感。而计算机显示设备不仅不接受辐射亮度值作为输入，而且还不能够显示非常亮或非常暗的颜色。在理想观看条件下，这些设备通常只能显示两个亮度变化的级别。

因为基于物理的渲染算法所生成的真实感图像同样受制于辐射亮度值和设备显示能力的不匹配，所以，为了使所显示的图像跟实际场景尽可能地接近，解决这个图像显示问题至关重要。这是一个非常活跃的研究领域，即所谓的色调映射(tone mapping)，来寻找将多出来的亮度等级加以压缩以便显示的好方法。该领域深入到对人类视觉系统的研究，以引导对图像显示技术的开发。通过对HVS性质的探讨，人们开发出了能够对显示设备能力进行良好补偿(compensating for display limitation)的色调映射算法。本节将介绍几个这类算法。

8.4.1光亮度和光度学(luminance and Photometry)

因为色调映射算法通常基于人眼对明亮度的感知，所有大多数色调映射算子使用光亮度(luminance)单位，光亮度用来测量一个光谱功率分布对人类观察者而言的明亮程度。例如，光亮度可以用来说明这样的事实：对于人类，有一定绿色波长能量的SPD要比具有相同的蓝色波长能量的SPD看上去要明亮一些。

光亮度跟辐射亮度有紧密的联系：给定一个光谱辐射亮度值，其光亮度值可以用一个简单的转换公式计算出来。实际上，第5章所定义的辐射度学中的所有量在光度学里都有对应。光度学(photometry)是研究可见的电磁辐射以及HVS对其感知的学问。每一个辐射度学中的量都可以转换成相应的光度学中的量，而所依据的就是对光谱分布和光谱反应曲线V(λ)的乘积进行积分，其中光谱反应曲线描述了人眼对不同波长的相对敏感度。

我们用Y来表示光亮度，它跟光谱辐射亮度L(λ)的关系是：

Y = ʃλL(λ)V(λ)dλ

光亮度和光谱反应曲线V(λ)跟颜色的XYZ 表示有紧密的联系(见第5.1.2节)。CIE所选定的三刺激值曲线Y(λ)跟V(λ)成比例关系，使得：

Y = 683 ʃλL(λ)Y(λ)dλ

这样一来，给定了XYZ表示，就可以用一个比例因子而得到每个图像像素的光亮度。光亮度的单位是坎德拉/平方米(cd/m2), 坎德拉(candela)是跟辐射强度相对应的色度学里的量。坎德拉/平方米又常被叫做尼特(nit)。下面是常见的一些光亮度值：

光亮度(nits)

600,000 地平线上的太阳

120,000 60瓦灯泡

8000 晴朗的天空

100-1000 典型的办公室

1-100 典型的计算机显示器

1-10 街道上的照明

0.25 有云遮掩的月光

人眼有两种类型的感光细胞：杆体细胞和锥体细胞。杆体细胞有助于在较暗的环境（即微光水平：10-6 到10 cd/m2)中的视觉活动。杆体细胞对颜色不敏感，并且不擅长分辨微小的细节。而锥体细胞可以处理从10-2到108cd/m2的光照范围（明视水平）。有三类锥体细胞，分别对不同的波长的光敏感。计算机显示器通常可以显示的光亮度范围大约是1到100 cd/m2。

8.4.2 光晕(Bloom)

在介绍色调映射算法之前，我们介绍一种可以欺骗HVS的技术，使HVS看到显示器上的图像要比实际的图像更明亮一些。如果人眼所看到的环境中的某一部分比其它部分明亮得太多，就会出现一种称为"光晕"(Bloom)的效果，即在明亮物体的周围区域有模糊的辉光。人们还不明白这其中的原因，但一般认为这是光在人眼内发生散射而产生的现象。计算机图形学研究者们发现在图像渲染中模拟这种效果会极大地提高图像的真实感。当图像的某一部分出现这种辉光时，HVS就很自然地感知到这部分图像要比其它部分明亮许多。

作为一种选项，我们所介绍的图像管线可以对图像施加bloom效果。这个效果所使用的滤波器是以实验为依据的，而非基于人类视觉系统的模型，但是在实际应用中效果良好。其基本的理念是使用一个非常宽的滤波器，可以作用到图像中所有的像素上，并且快速衰减。因为这个滤波器有很宽的支撑，那些很亮的像素对周围的许多像素都有能量贡献。因为它衰减很快，就不会改变那些有相似亮度值的图像区域，但非常亮的像素可以使低滤波权值失效并影响到其它像素。然后我们再用一个用户提供的权值，将这个bloom图像跟原图像合成。

该滤波器有两个参数：bloomRadius，指定滤波器覆盖图像的比例值；bloomWeight，指定用来混合bloom图像和原图像的权值。如果bloomRadius为0，滤波器就被视为关闭。在实际应用中，使用大约为0.1或0.2的值是使用该滤波器很好的初始选择。

<Possibly apply bloom effect to image> =

if (bloomRadius > 0.f && bloomWeight > 0.f) {

<Compute image space extent of bloom effect>

<Initialize bloom filter table>

<Apply bloom filter to image pixels>

<Mix bloom effect into each pixel>

<Free memory allocated for bloom effect>

}

首先，要确定滤波器的宽度（以像素个数计）。用bloomRadius乘以图像x,y方向的分辨率的最大值，就得到这个值。

< Compute image space extent of bloom effect> =

int bloomSupport = Float2Int(bloomRadius * max (xResolution, yResolution));

int bloomWidth = bloomSupport / 2;

因为要对滤波器函数做很多次求值，有必要预先计算好一个表。这里的实现所用的是一个放射式对称的滤波函数：

f (x, y) = ( 1 - (x2+y2)1/2 / d) 4

其中d是滤波器的宽度。该滤波器是不可分的(即不可分离成两个1D滤波函数的乘积），所以每个要输出的像素需要滤波的像素个数等于滤波器宽度的平方，这样一来，预先计算好滤波值表就更有必要了。

<Initialize bloom filter table> =

float *bloomFilter = new float[bloomWidth * bloomWidth];

for (int i = 0; i < bloomWidth * bloomWidth; ++i) {

float dist = sqrtf(float(i)) / float(bloomWidth);

bloomFilter = powf( max(0.f, 1.f - dist), 4.f);

}

这里的实现先计算一个暂时的图像来存放bloom贡献值。我们不能在计算bloom值的时候更新原图像，这一点很重要：因为这样做会使用已经累加过bloom值的像素值来计算附近像素的bloom值，从而得到错误的结果。

<Apply bloom filter to image pixels> =

float *bloomImage = new float[3 * nPix];

for (int y = 0; y < yResolution; ++y)

for (int x = 0; x < xResolution; ++x) {

<Compute bloom for pixel(x,y)>

}

}

为了计算bloom图像中的像素值，先要找到所有的可能对其有贡献值的那些像素。然后，对这些像素进行滤波计算：

<Compute bloom for pixel(x,y)> =

<Compute extent of pixels contributing bloom>

int offset = y * xResolution + x;

float sumWt = 0.;

for (int by = y0; by <= y1; ++by)

for (int bx = x0; bx <= x1; ++bx) {

<Accumulate bloom from pixel(bx, by)>

}

bloomImage[3*offset ] /= sumWt;

bloomImage[3*offset+1] /= sumWt;

bloomImage[3*offset+2] /= sumWt;

我们以当前像素位置为中心，在每个方向上偏移，偏移量为滤波器宽度，就可以找到有贡献值的像素范围：

<Compute extent of pixel contributing bloom> =

int x0 = max(0, x - bloomWidth);

int x1 = min(x+bloomWidth, xResolution-1);

int y0 = max(0, y - bloomWidth);

int y1 = min(y+bloomWidth, yResolution-1);

当前像素并不参与bloom计算，因为我们要计算的是其它像素对该像素的贡献值。

<Accumlate bloom from pixel(bx, by)> =

int dx = x - bx, dy = y - by;

if (dx == 0 && dy == 0) continue;

int dist2 = dx * dx + dy * dy;

if (dist2 < bloomWidth * bloomWidth) {

int bloomOffset = bx + by * xResolution;

float wt = bloomFilter[dist2];

sumWt += wt;

for (int j = 0; j < 3; ++j)

bloomImage[3 * offset + j] += wt * rgb[3*bloomOffset+j];

}

计算出bloom图像之后，就可以用bloomWeight值跟原图像进行混合：

<Mix bloom effect into each pixel> =

for (int i = 0; i < 3 * nPix; ++i)

rgb = Lerp(bloomWeight, rgb, bloomImage);

<Free memory allocated for bloom effect> =

delete[] bloomFilter;

delete[] bloomImage;

8.4.3 色调映射接口

色调重现的基本方法是求出一个将像素值映射到可显示的范围的比例变换函数。对于简单的色调映射算子而言，常常使用一个简单的作用于图像中所有像素的函数。这类算子被称为空间均匀（或全局）算子。它们给出了从图像光亮度到显示光亮度的单调映射。更复杂的方法所使用的函数会根据每个像素的亮度和附近像素的亮度而发生变化，称之为有空间变化的（或局部）算子，它们并不一定是单调映射。

有空间变化的算子会比空间均匀算子更具成效，这一点很有趣。这些算子之所以效果良好，是因为人眼对局部的对比要比绝对的光亮度更敏感。由于这个特性，就可以对图像中不同部分赋予不同的像素值，虽然这些图像不同的部分原有相同的绝对光亮度值（却被算子赋给了不同的值），但对于人类观察者而言，应用算子后的结果并没有什么不妥。

因此，许多色调映射算子的基本目标就是要尽可能地在被显示的图像中保留局部对比，而不是保留绝对的亮度。在所有图像（亮的，或暗的）区域中确保使用有足够差别的颜色至关重要，观察者这样才能看到不同的颜色，并且保证大范围的图像强度不会被映射到相同的像素值上。所以说，一个在场景中比另一个物体明亮两倍的物体在显示器上并不一定也有两倍的亮度。再重复一下，有对比的局部变化对人类视觉系统来说是最为重要的。

HVS的对光亮度变化的总体上的敏感性是根据适应亮度(adaption luminance)来变化的，我们用记号Ya来表示适应亮度。适应亮度在图像不同的部分可以有不同的值。在下面的讨论中，我们用显示器适应亮度(display adaption luminance, Yad)来表示人类观察者观看计算机显示器时的适应亮度，用世界适应亮度（world adaption luminance, yaw)来表示人观察实际场景时的适应亮度。

由于在昏暗环境下杆体细胞占主导地位，HVS在黑暗中的特性就全然不同。例如，对颜色的感知减弱了，所有的东西看上去只是不同程度的深灰色。还有，空间视敏度（spatial acuity）也被减弱了：在1000 nits的适应亮度下，HVS可以分辨50周/视角的空间细节，而在0.001 nits的情况下，只能分辨大约2.2周/视角。用于微光水平的色调重现算子常常用一些模糊技术来模拟这个效果。

所有的色调映射算子都是从ToneMap基类继承下来的，该基类提供了TopMap::Map()接口函数：

<ToneMap Declarations> =

class ToneMap {

public:

<ToneMap Interface>

};

TopMap::Map()函数的参数包括一个指向图像像素光亮度值的数组指针，图像分辨率，显示器可显示的最大光亮度值。该函数负责为每个像素计算出一个比例因子，并把它存放在scale数组中。所计算出的比例因子要使被因子乘过的像素光亮度值位于[0,maxDisplayY]范围之内。

<ToneMap Interface> =

virtual void Map(const float *y, int xRes, int yRes,

float maxDisplayY, float *scale) const = 0;

如果用户为ApplyImaingPipeLine()提供了色调映射算子，我们就用下列片段来实施色调重现操作。首先，要为每个像素计算光亮度值，然后为每个像素计算比例因子，然后对图像进行比例变换。显示器所能显示的最大光亮度值maxDisplayY被设为固定值100.f。

<Apply tone reproduction to image> =

Tone *toneMap = NULL;

if (toneMapName)

toneMap = MakeToneMap(toneMapName,

toneMapParams ? *toneMapParams : ParamSet());

if (toneMap) {

float maxDisplayY = 100.f;

float *scale = new float[nPix], *lum = new float[nPix];

<Compute pixel luminance values>

toneMap->Map(lum, xResolution, yResolution,

maxDisplayY, scale);

<Apply scale to pixels for tone mapping and map to [0.1] >

delete[] scale;

delete[] lum;

}

ApplyImagingPipeLine()的yWeight变量给出了计算RGB像素值的权值。如果它为NULL，就用标准值。

<Compute pixel luminance values> =

float stdYWeight[3] = { 0.212671f, 0.715160f, 0.072169f};

if (!yWeight) yWeight = stdYWeight;

for (int i = 0; i < nPix; ++i)

lum = 683.f * (yWeight[0] * rgb[3*i] +

(yWeight[1] * rgb[3*i+1] +

(yWeight[2] * rgb[3*i+2]);

因为色调映射算子返回的比例因子使像素光亮度值lum落在范围[0, maxDisplayY]之间，而显示设备并不用光亮度值为输入值，所以我们还要把结果再缩比为[0,1]之间的值。另外，乘683这个因子很重要，因为色调映射算子假定该因子已经被乘到像素值上了。

<Apply scale to pixels for tone mapping and map to [0,1] > =

float displayTo01 = 683.f / MaxDisplayY;

for (int i = 0; i < xResolution * yResolution; ++i) {

rgb[ 3*i ] *= scale * displayTo01;

rgb[ 3*i+1 ] *= scale * displayTo01;

rgb[ 3*i +2] *= scale * displayTo01;

}

8.4.4 映射最大值为白色

最直截了当的色调重现算子就是“最大值到白色”算子。它对所有像素进行循环，寻找最大的光亮度值，并把所有像素值进行均匀地映射，使得最大光亮度值被映射到显示器的最大光亮度值。

<MaxWhiteOp Declarations> =

class MaxWhiteOp: public ToneMap {

public:

<MaxWhiteOp Public Methods>

};

<MaxWhiteOp Public Methods> =

void Map(const float *y, int xRes, int yRes, float maxDisplayY, float *scale) const {

<Compute maximum luminance of all pixels>

float s = maxDisplayY / maxY;

for (int i = 0; i < xRes * yRes; ++i)

scale = s;

<Compute maximum luminance of all pixels> =

float maxY = 0;

for (int i = 0; i < xRes * yRes; ++i)

maxY = max(maxY, y);

在实际应用中，该算子有两个缺点、首先，它根本没有考虑人类视觉系统：如果将场景中的光加亮100倍再重新渲染，使用该算子所得到的结果全然相同。其次，如果有很少的几个极亮的像素，也会导致图像的其它部分非常暗淡，以至于看不见。但是，它对那些没有太多的动态范围的图像而言还是效果不错的，也可以作为一种基准，来评估更复杂的算子所带来的改进。

8.4.5 基于对比的比例因子

我们所介绍的下一个算子注重保留图像中的亮度对比。它是由Greg Ward(1994a)研发出来的。该算子基于那些研究HVS并开发出仿真模型的研究者们的工作，所使用的模型描述了在给定适应亮度的情况下观察者所能察觉到的最小的光亮度变化--最小可觉差(JND, just noticeable difference)。适应亮度越大，能够被察觉到的光亮度变化也就越大。该算子是这样设置图像光亮度值的：即试图将被显示图像的一个JND和实际环境中的一个JND相对应。

这个均匀比例因子是这样试图保留对比可见性的--给定一个原图像中刚好能被观察者区分开来的区域，它试图对显示像素值进行比例变换，使得观看显示器的人刚好分辨出两个像素值的不同。一个能够增大JND的比例因子是对宝贵的显示动态范围的一种浪费，而减小JND的比例因子又会丢失掉视觉上可以观察出的特征。

研究者们发现，对于给定的一个在明视范围中的一个适应亮度，最小的可见的光亮度变化可以由下面的模型给出：

ΔY(Ya) = 0.0594(1.219 + (Ya)0.4)2.5

该算子的比例因子s要满足：

ΔY(Yad) = s ΔY(Yaw)

其中Yad是显示适应亮度，而Yaw是世界适应亮度。

解出s值，就得到：

s = ((1.219 + (Yad)0.4) / (1.219 + (Yaw)0.4) )2.5

<ContrastOp Declarations> =

class ContrastOp(float day) { displayAdapationY = day;}

void Map(const float *y, int xRes, int yRes, float maxDisplayY, float *scale) const ；

float displayAdaptionY;

}

<ContrastOp Method Definitions> =

void Map(const float *y, int xRes, int yRes, float maxDisplayY, float *scale) const

{

<Compute world adaptation luminance>

<Compute contrast-preserving scale factor, s>

for (int i = 0; i < xRes * yRes; ++i)

scale = s;

}

有一个没有解决的问题就是如何计算世界适应亮度Yaw。最为理想的情况是，该值的计算基于观察者所看到的那一部分场景和观看的时间（HVS要花时间适应亮度的变化）。由于缺乏这类信息，算子只是计算原图像所有光亮度值的log平均值。使用log平均值而非一般平均值的原因是防止很小的明亮区域对图像的其它部分起压倒性的作用。

<Compute world adaptation luminance, Ywa> =

float Ywa = 0;

for (int i = 0; i < xRes * yRes; ++i)

if (y > 0) Ywa += logf(y);

Ywa = expf(Ywa / (xRes * yRes));

利用上面的公式可以算出比例因子s:

<Compute contrast-preserving scale factor, s> =

float s = powf((1.219f + powf(displayAdaptationY, 0.4f)) /

(1.219f + powf(Ywa, 0.4f)), 2.5f);

该算子一般情况下工作良好，但在维持明亮区域的细节方面却麻烦多多。这也难怪，因为任何使用单一比例因子的算子都有这个缺点。该算子适合于室内场景，而且计算效率高。

8.4.6 可变化的适应亮度

前面提到过，我们常常可以使用一个可变的比例因子对显示器的动态范围加以更好地利用。这里我们将实现一个色调重现算子，使之特别适用于具有许多亮度变化等级的高对比的场景。它计算一个在整个图像上面平滑变化的局部适应亮度。然后使用这个局部适应亮度和一个保持亮度对比的色调重现算子来计算一个比例因子，其方法跟前面定义的ContrastOp算子相似。

计算空间上可变的局部适应亮度的困难之处在于，在非常明亮的区域和非常暗淡的区域之间的边界上，很容易产生人为的缺陷。如果色调重现算子使用了受明亮像素影响的适应亮度，那么暗淡的像素值就被映射为黑色，就会在边界上产生光圈效应（halo effect）。

所以，更好的方法是要保证暗淡的像素所使用的适应亮度是仅仅基于附近的暗淡像素的。本节所介绍的算子使用了一种图像处理技术来探测那些适应亮度相差过大的区域之间的边界。在那些适应亮度变化缓慢的局部区域，该算子使用一个局部的可以保持亮度对比的比例因子--即明亮的区域不至于被变为白色，暗淡的区域不至于被映射为黑色。

<HighContrastOp Declarations> =

class HighContrastOp : public ToneMap {

public:

void Map(const float *y, int xRes, int yRes, float maxDisplayY, float *scale) const;

private:

< HighContrastOp Utility Methods>

};

这个色调映射函数基于亮度感知门限(Threshold Versus Intensity, TVI)函数，由它可以得出给定适应水平的最小可觉差,TVI(Ya)。在本质上，它跟ContrastOp所用的JND函数相似，但是它基于更复杂的人类视觉系统模型上的，其中包括在微光水平下的反应模型。

从TVI函数出发，我们可以定义感知容量(perceptual capacity)，即对于给定的适应水平Ya，光亮度范围（Ya,Yb）所覆盖的JND数：

(Ya - Yb) / TVI(Ya)

稍后我们会利用上述关系式来重新映射图像的局部区域，使得被显示的图像仍然保持了感知容量。

为了快速地计算出一对亮度值的感知容量，我们定义下面的辅助容量函数C(Y):

C(Y) = ʃY 0(1 / TVI(Y')) dY'

其中我们做了这样的近似：对于给定光亮值，计算感知容量微分的适应水平值等于该光亮度值。这样，C(Ya) - C(Yb) 就是从Ya到Yb的感知容量。

Ashikhmin对广泛使用的TVI函数做了某些简化，使得更容易计算C(Y)，其结果如下：

<HighContrastOp Utility Methods> =

static float C (float y) {

if (y < 0.0034f)

return y / 0.0014f;

else if (y < 1)

return 2.4483f + log10f(y / 0.0034f) / 0.4027f;

else if (y < 7.2444f)

return 16.563f + (y - 1) / 0.4027f;

else

return 32.0693f + log10f(y / 7.2444f) / 0.0556f;

}

有了C(Y)函数后，我们就可以用给定的光亮度值Y作为参数，确定它与图像中的最小光亮度之间有多少个JND：

C(Y) - C(Ymin)

我们也可以计算出该JND数占整个图像的JND数的比率：

(C(Y) - C(Ymin)) / (C(Ymax) - C(Ymin))

我们就会知道这个世界光亮度要被映射到显示光亮度范围的位置。这样，总体上的色调映射算子如下：

T(Y) = Ymaxd (C(Y) - C(Ymin)) / (C(Ymax) - C(Ymin))

<HighContrastOp Utility Methods> +=

static float T( float y, float CYmin, float CYmax, float maxDisplayY) {

return maxDisplay * (C(y) - CYmin) / (CYmax - CYmin);

}

现在我们可以定义色调映射主函数了。该函数计算图像中所有像素的最大(Ymax)、最小光亮度值(Ymin)。在计算适应亮度时，为了加快搜索速度，该算法还要建立一个金字塔型的图像数据结构，其中原图像要被逐步地被滤波到更低分辨率的自身拷贝中。当算子对所有像素进行循环并计算每个像素的比例因子时，就要用到这个金字塔结构。

<HighContrastOp Method Definitions> =

void Map(const float *y, int xRes, int yRes, float maxDisplayY, float *scale) const {

<Find minimum and maximum image luminances>

<Build luminance image pyramid>

<Apply high-contrast tone mapping operator>

}

<Find minimum and maximum image luminances> =

float minY = y[0], maxY = y[0];

for (int i = 0; i < xRes * yRes; ++i){

minY = min(minY, y);

maxY = max(maxY, y);

}

float CYmin = C(minY), CYmax = C(maxY);

前面所介绍的计算局部适应亮度的大多数方法使用了原图像的一个模糊版本，这样会引起光圈效应。这里所实现方法的背后的洞见是，适应亮度不应该基于围绕像素(x,y)的大小固定的区域，而应该是一个可变的区域：只要光亮度在局部大致保持相等，就可以扩张这个区域，直到遇到明显的光亮度变化为止。这就拥有了两方面的优点：当光亮度变化缓慢时，就可以为一个较大的区域计算适应亮度，就给出了远离具有高对比度的特征的适应亮度的平滑变化。当对比度快速变化时，也能够探测出这种变化，通过计算更局部的适应亮度来避免人为的缺陷。

图像处理的一个标准技术是定义像素的局部对比度lc(x,y)，即像素在两个被模糊的图像版本中的像素值差，其中一个版本的滤波器宽度是另一个的两倍：

lc (s, x, y) = (Bs(x, y) - B2s(x,y)) / Bs(x,y)

这里s是用于模糊图像的滤波器宽度（以像素个数计）, Bs(x, y)是像素(x,y)在模糊图像中的值。我们希望找到围绕每个像素(x,y)的最大局部范围半径s，使得|lc (s, x, y)|小于某个固定值。当它比这个固定值大时，就说明我们超过了可以接受的局部对比度。找到了满足条件的s后，就可以用下面的公式计算适应亮度：

Ya(x,y) = Bs(x,y)

为了快速地找到在模糊图像中的像素值Bs(x,y)，该算子用MIPMap类创建一个图像金字塔，该类在第11.4.2节中被介绍。在本节中，我们只需知道它可以精确而且有效地计算出任意s值所对应的Bs(x,y)。

<Build luminance image pyramid> =

MIPMap<float> pyramid(xRes, yRes, y, false, 4.f, TEXTURE_CLAMP);

下一步计算每个像素的适应亮度和比例因子。注意我们需要把离散像素坐标转换为连续的像素坐标，用来MIPMap查找。

<Apply high-contrast tone mapping operator> =

for (int y = 0; y < yRes; ++y) {

float yc = (float(y) + .5f) / float(yRes);

for (int x =0; x < xRes; ++x) {

float xc = (float(x) + .5f) / float(xRes);

<Compute local adaptation luminance at (x,y)>

<Apply tone mapping based on local adaptation luminance>

}

}

为了计算适应亮度，该函数在模糊图像中查找像素值，并计算局部对比度函数值lc。如果该值大于maxLocalContrast，0.5（一个任意常量，根据经验而定），适应亮度就设置为稍微小一些的区域的平均值，并结束循环。否则，模糊半径就增加一个像素的跨度dwidth，并计算lc的新值。一旦该过程到达了一个很大的模糊半径却没有找到足够大的对比度，就终止循环。

<Compute local adaptation luminance> =

float dwith = 1.f / float(max(xRes, yRes));

float maxWidth = 32.f / float(max(xRes,yRes));

float width = dwidth, prevWidth = 0.f;

float Yadapt;

float prevlc = 0.f;

const float maxLocalContrast = .5f;

while(1) {

<Compute local contrast at (x, y)>

<If maximum contrast is exceeded, compute adaptation luminance>

<Increase search region and prepare to compute contrast again>

}

有了MIPMap图像金字塔，就很容易算出局部对比度了。MIPMap::Lookup()在给定的像素位置上使用一个给定宽度的滤波器。

<Compute local constrast at (x,y)>=

float b0 = pyramid.Lookup(xc, yc, width, 0.f, 0.f, width);

float b1 = pyramid.Lookup(xc, yc, 2.f * width, 0.f, 0.f, 2.f * width);

float lc = fabs((b0 - b1) / b0);

如果局部对比度超出了所允许的最大对比度值，就近似地为宽度s确定一个值，使得lc(s,x,y) = maxLocalConstrast。我们把上次循环得到的局部对比度存放在prevlc, 则有下列关系：

prevlc < maxLocalContrast < lc

我们可以找到插值量t，使之满足在prevlc和lc之间做线性插值恰好等于maxLocalContrast。我们假定在上次循环的宽度和当前宽度之间对比度呈线性变化，就可以容易地算出对比度约束被打破的那个宽度。

<If maximum contrast is exceeded, compute adaptation luminance> =

if (lc > maxLocalContrast) {

float t = (maxLocalContrast - prevlc) / (lc - prevlc);

float w = Lerp (t, prevWidth, width);

Yadapt = pyramid.Lookup(xc, yc, w, 0.f, 0.f, w);

break;

}

<Increase search region and prepare to compute contrast again> =

prevlc = lc;

prevWidth = width;

width += dwidth;

if (width > maxWidth) {

Yadapt = pyramid.Lookup(xc, yc, maxWidth, 0.f, 0.f, maxWidth);

break;

}

给定了色调映射函数T(Ya)，在像素(x,y)的比例因子定义为：

s(x,y) = T(Ya (x,y)) / Ya (x,y)

只要Ya (x,y)变化缓慢，这实际上是一个局部的线性映射。

<Apply tone mapping based on local adaptation luminance> =

scale[x + y*xRes] = T(Yadapt, CYmin, CYmax, maxDisplayY) / Yadapt;

8.4.7 有空间变化的非线性比例变换

最后要实现的色调映射方法并不基于任何关于感知的理论文献，却在实际应用中效果良好。它是基于一个特定公式，是由Reinhard等人发表的。该方法使用空间可变的因子来对每个像素做比例变换：

注意这个算子并不基于光亮度Y，而是XYZ颜色的y分量（也就是说，并没有包括比例因子683）。

这个比例因子把黑色像素映射到0，把最亮的像素映射到1。在两者之间，跟明亮的像素相比，暗淡一些的像素需要相对少一些的亮度变化就能够带来对所输出的像素值的变化。这是跟HVS性质相一致的，因为HVS的反应曲线是对数函数型的，而非线性的。

< NonLinearOp Declarations> =

class NonLinearOp: public ToneMap {

public:

<NonLinearOp Public Methods>

private:

float maxY;

};

构造器可以选择使用适应亮度。如果没有提供该值，则将像素光亮度的log平均值做为适应亮度。

< NonLinearOp Public Methods> =

NonLinearOp(float my) { maxY = my;}

算子的实现是很容易的，唯一的麻烦是需要将因子683从亮度值中去掉。

<NonLinearOp Public Methods> +=

void Map(const float *y, int xRes, int yRes, float maxDisplayY, float *scale) const {

float invY2;

if (maxY <= 0.f) {

<Compute world adaptation luminance, Ywa>

Ywa /= 683.f;

invY2 = 1.f / (Ywa * Ywa);

}

else invY2 = 1.f / (maxY * maxY);

for (int i = 0; i < xRes * yRes; ++i) {

float ys = y / 683.f;

scale = maxDisplayY / 683.f * (1.f + ys * invY2) / (1.f + ys);

}

}

8.5 图像管线的最后阶段

不幸的是，有许多颜色在现代的显示设备中无法再现，例如饱和的橘黄色和紫色，这些颜色被称为“超出色域”(out of gamut)。超出色域的颜色被映射到RGB空间后的值超出了[0,1]范围。关于如何处理超出色域的颜色，并没有令人满意的方法。既然显示设备无法显示这些颜色，需要做的只能是在不同的颜色选择所产生的错误中寻找折中方案。这里的实现方法是对超出色域的颜色重新做比例映射，将三个分量中最大的分量设置为1，其它两个分量做相应的比例变换。

<Handle out-of-gamut RGB values> =

for (int i = 0; i < nPix; ++i) {

float m = max(rgb[3*i], max(rgb[3*i+1], rgb[3*i+2]));

if (m < 1.f)

for (int j = 0; j < 3; ++j)

rgb[3*i+j] / = m;

}

一旦所有的颜色都被映射到可显示范围后，还需为CRT显示器的亮度反应特性调整颜色值。对于这类显示器而言，被显示的亮度并不跟像素值呈线性关系：像素值为100的像素并不比值为50的像素亮两倍。虽然像LCD这类的新型显示设备本质上没有非线性反应特性，但它们一般还是在模拟CRT设备的这个反应特性。

这个非线性反应可以用一个幂函数来表示：

d =v ɣ

其中d是显示亮度，v是加到电子枪上的电压，gamma值ɣ的典型值一般为2.2。因此，给定像素值v，需要把下面的值提供给显示器才可以使被显示的颜色有线性关系：

v' = v 1/ ɣ

<Apply gamma correction to image> =

if (gamma != 1.f) {

float invGamma = 1.f / gamma;

for (int i = 0; i < 3 * nPix; ++i)

rgb = powf(rgb, invGamma);

}

然后，我们将经过gamma校正后的像素值映射到显示器所需要的范围（一般是0到255）。参数maxDisplayValue用来控制这个映射：

<Map image to display range> =

for (int i = 0; i < 3 * nPix; ++i)

rgb *= maxDisplayValue;

最后，在把这些像素值映射到整数值之前，可以对这些值进行抖动(dither)，即对每个像素的颜色分量加上一个随机的微小值。引入这个随机噪声的目的是可以使得同一种颜色的区域到另一个区域的过渡更自然些，从而改进了被显示图像的视觉效果。没有做抖动处理的图像会在那些颜色的过渡区域出现带状的缺陷。

<Dither image> =

if (dither > 0.f)

for (int i = 0; i < 3*nPix; ++i)

rgb += 2.f * dither * (RandomFloat() - .5f);