Mathematica图像处理

MD Document：2/25/2016 4:05:56 AM by Jimbowhy
CSDN PuBlISheD: http://blog.csdn.net/WinsenJiansbomber/article/details/50769511

当我听着大佑的歌《將進酒》，写着《Make自动化编译工具》，想着自己的存在，感觉很好。突发想法觉得音量是不太小了？就要去调大点音量，然而立刻有个声音告诉我，好吧，其实音量刚刚好，因为刚刚感觉好就是证明啊！—— by Jimbowhy 2/29/2016 8:46:50 AM

Mathematica 和 Matlab 一样，是世界上数一数二的科学数值编程软件，最先接触到 Mathematica 是在2000的时候，那时计算机还讲486的年代，Mathematica 的版本好像还是 3.0 的时候，但是它的程序界面先出来一个像分形一样带有尖刺的球体留下了深刻的印象。它就像是一个童话里的世界，永远都忘不了，童话世界？对我就是这么说的。可惜当年不会用它，因为还没有什么数学常识，后来正式使用它的时候已经是10年后了，直到上了大学参加数学建模的时候需要。《Mathematica全书》就是在那时看完的，它就是 Mathematica 的帮助文档的中文版，有英文基础可以直接上文档，现在使用的 Mathematica 9.0 中文版已经包含中文帮助文档。

基本功能

MM的功能十分强大，可以上图，也可以玩音乐，比如说这个 fur elise:

BPM = 70; W = 60/BPM; O == W/8; Q = W/4; H = W/2;
Sound[{SoundNote["E6", Q], SoundNote["Eb6", Q],SoundNote["E6", Q, "Violin"], SoundNote["Eb6", Q], SoundNote["E6", Q],SoundNote["B5", Q], SoundNote["D6", Q], SoundNote["C6", Q], SoundNote["A5", H],SoundNote["A3", Q], SoundNote["E4", Q], SoundNote["A4", Q],SoundNote["C5", Q], SoundNote["E5", Q], SoundNote["A5", Q], SoundNote["B5", Q],SoundNote["E3", Q], SoundNote["E4", Q], SoundNote["Ab4", Q],SoundNote["E5", Q], SoundNote["C6", Q], SoundNote["B5", Q],SoundNote[{"E5", "A3"}, Q], SoundNote["A5", Q],SoundNote["B5", Q], SoundNote["C6", Q], SoundNote["D6", Q],SoundNote[{"C4", "E6"}, H], SoundNote["C5", Q],SoundNote["G5", Q], SoundNote["F6", Q], SoundNote["E6", Q], SoundNote["D6", H],SoundNote["F5", Q], SoundNote["E6", Q], SoundNote["D6", Q], SoundNote["C6", H],SoundNote["E5", Q], SoundNote["D6", Q], SoundNote["C6", Q], SoundNote["B5", H],SoundNote["E5", Q], SoundNote["E5", Q]}]
Export["C:\\writing\\Mathica\\fur_elise.mid", %, "Sound"]

毕竟不是传做音乐的货啊，在时值不同的音符组合时就没有办法处理，还是使用 javax.sound.midi 包实现MDI播放的，开发人员有点懒。

这里就主要介绍它的图片处理功能。通过 Import/Export 命令可以导入导出各种图片格式文件，当然这两个命令还有其它更多的用途。新建一个作业本，输入以下命令，通过 Shift+Enter 执行它，就可以加载自带的 Lena 人物图片，同时将文件导出到我正在写作的文件夹。

i = Import["ExampleData/lena.tif"]
Export["C:\\writing\\Mathica\\lena.png", i]

然后就可以用能MD写作工具将图片文件放到MD文档中，像这样，这就是去Office时代的便利：

![Girl Lena][100]
[100]: lena.png "Girl Lena"

Import还支持HTTP协议，可以从网络上加载图片数据、音频、视频等等，在导入时，可以通过 Elements 来查询支持的导入元素，例如查询AVI文件支持的导入元素就是有 ImageList，这个在后面制作 GIF 图片时有用处。也可以导入视频的指定帧：

In[59]:= Import[ "C:\writing\Mathica\a.avi", "Elements"]Out[59]= {"Animation", "BitDepth", "ColorSpace", "Data", "Duration", "FrameCount", "FrameRate","Frames", "GraphicsList", "ImageList", "ImageSize", "VideoEncoding"}In[65]:= Import[ "C:\writing\Mathica\a.avi", {"Frames", {3, 5, 9}}]

Mathematica 每次执行命令时都会产生一个输出，有一个数字号码即 In2 方括号提示的数值，通过 % 可以引用给其它需要输入数据的函数。

Mathematica 可以生成动画序列，表达式中的 /. 表示正则替换，用它来更新画圆的位置,设置旋转的角度：

ani = Animate[Graphics[{EdgeForm[Directive[Dashed, Red]],Black, Disk[{0, 0}, 2, {0, Pi}],EdgeForm[Directive[Dashed, Red]],White, Disk[{0, 0}, 2, {Pi, 2 Pi}],EdgeForm[None],Black, Disk[{-1, 0}, 1],White, Disk[{1, 0}, 1],Red, Disk[{-1, 0}, 0.3],Red, Disk[{1, 0}, 0.3],},PlotRange -> 2.1, ImageSize -> 150] /. Disk[x__] :> Rotate[Disk[x], d Degree, {0, 0}], {d, 0, 360}]

也可以生成的序列图片保存为GIF动画文件，Animate 生成 gif 需要点技巧，先保存为动画视频再导入视频的帧，最后将序列帧导出为gif。哦，很抱歉，写到这里得到一条消息，外公去世了，老人家卧床多年，活动不便，可说得上是久经折磨，我想这算是个好消息吧。虽说病痛不能分担，但可以爱惜健康减少病痛，而对于生老却是无能为力的事，如果人生得意就尽欢吧，幸福的事就动画一样可以重复播放。

Ani2Gif[ani_, name_String, step_Integer] := Export[name, Import[ Export[name <> If[$OperatingSystem, "MacOSX", ".mov", ".avi"], ani], "ImageList"][[1 ;; -1 ;; step]]]Ani2Gif[ani, "C:\\writing\\Mathica\\the_eight_diagrams.gif",1]

表达式中的 [[1 ;; -1 ;; step]]] 是将图片序列元素重新按 step 间隔提出来组成新的序列，它的语法结构如下：

expr[[m;;n;;s]] 给出从第 m 个元素到第 n 个元素且步长为 s 组成的子列表

通常导出的图片带有控件，如果不想要控件部分可以通过 ImageCrop 函数进行修剪，小边10像素，大边52像素就差不多了：

$gif = "C:\\writing\\Mathica\\the_eight_diagrams.gif"
size = Import[$gif, "ImageSize"]
img = Import[$gif, "ImageList"][[1 ;; 30 ;; 1]]
ImageCrop[ ImageCrop[ img, size + {-10, -10}, {Right, Bottom}], size + {-20, -62}, {Left, Top}]

太极两仪GIF动画

另外 Animate 动画导出时是带有正反向双重图片的，所以一为了节省内存，去掉一半就可以了。

gif = "C:\\writing\\Mathica\\the_eight_diagrams.gif";
size = Import[gif, "ImageSize"];
frames = Import[gif, "ImageCount"]
img = Import[gif, "ImageList"][[1 ;; frames/2 ;; 1]];
F[x_] := ImageCrop[ImageCrop[x, size + {-10, -10}, {Right, Bottom}], size + {-20, -62}, {Left, Top}]
croped = Map[F, img]
Export[StringReplace[gif, ".gif" -> "-crop.gif"], %]

有些时候，在编写语句的时候需要将语句中的某些内容进行预处理，如果因为这样而定义一个函数显得有点繁琐，可以使用 Function 来定义内联的函数。MM 还提供了一个参数传递占位功能 Slot，可以在语句中设置参数点位符，在后面再接收参数传递。如下例子中，StringJoin 函数放置了一个参数占位 #，在函数结束方括号后使用了 & 来接收参数列表，这样在执行时参数就会流入 # 号的位置。#还可以带序号，用来引用参数列表中指定序号的参数。

Map[StringJoin[#, "o"] &, {"x", "y", "z"}]

分号在MM系统中是列表求值的意思，它还可以用来消隐不希望输出的内容。另外在遇到需要分段定义的函数时，可以使用条件测试 Condition /; 来设置分段条件，例如以下定义的一个分段函数：

F[x_] := x*10 /; x < 10
F[x_] := x^1.2 /; x < 30
F[x_] := x
Plot[F[x], {x, -10, 40}, PlotRange -> Automatic]

当变量较多比较混乱时，为了避免变量的相互影响，可以使用 Module 来限制变量的作用域，例如定义函数时就可以这样用：

f[x0_] :=
Module[{x = x0},While[x > 0, x = Log[x]];x
]

色彩模型

不同的图片具有不同的色彩空间，计算机上典型的色彩空间就有 RGB，这种色彩模式下，每个像素都由 Red/Green/Blue 三个纯色分量构成。在常用的24-bit颜色显示系统中每个分量用8-bit表示，取值范围是 0-255，可以表示160万种色彩，各个分量通过显示设备混合后形成人眼可以观察的颜色感觉信号。而人眼接收到的颜色信号其实是色调 Hue、饱和度 Saturation 和亮度 Brightness。从色觉上看，可以接收的光波从 300nm 到 900nm 波长范围，颜色分为类有红橙黄绿青蓝紫 red, orange, yellow, green, cyan, blue, violet。

为了更直观地反映人眼的色觉系统，美国色彩学家孟塞尔 H.A.Munseu 于1915年提出HSI模型是，以色调、饱和度和强度 Intensity 三种基本特征量来反映了人的视觉系统感知彩色的方式。由此衍生出了几个类似的模型，有 HSL，HSB，HSV，BLV 分别表示 brightness/lightness/value，HSV和HSB是同一个模式的不同叫法。这里了解一下和色彩光学相关的几个词汇：

Hue           色调，视觉感官产生的红绿蓝等等颜色的感觉；
Intensity，
radiance      光强，光线照射在一个地方的能量大小；
Luminance (Y) 流明，光密度的度量，单位是坎德拉/平米，candela per square meter (cd/m2)；
Luma (Y′)     卢马，R′, G′, B′ 伽马射线较正使用的度量；
Brightness    亮度，视觉感官对光照能量大小产生的感觉；
Colorfulness  色彩，视觉感官对图像各种色彩混合产生的感觉；
Chroma        色度，色彩和光强度；
Saturation    饱和度，对每种色彩来说，饱和度越高，这种色就越鲜艳明显。

对于HSL、HSV模型可以用一段圆柱体来表达，也可以用两个拼接在一起的圆锥来表达。而RGB则可以用一个立方体表达，以一角为起点O，RGB三分量都为0，颜色为黑色，辐射开的三条边所连接的角为RGB三分量的最大值，O的对角O’就是RGB三分量同时为最大值的点即白色。可以想像，将O-O’的对角线立起来，那么就可以在下面投影出一个六边形色彩平面，这条对角线就是灰度线：

而圆柱模型竖直方向是强度分量，取值范围 0-1往上表示颜色越亮越白，往下表示颜色越暗越黑；半径方向是饱和度分量，越近中心线表示色彩越淡，中心线就表示灰度色；旋转角度表示色调分量，0°度表示红色调，120°度表示绿色调，240°表示蓝色调。各分量取值范围 s∈[0,1]，l∈[0,1]，h∈[0, 360°)。采用HSI模型的一个好处就是，通过改变饱和度或色调分量就可以改变图像的整体颜色，可以很方便地调制出具有特定情感色彩的图像，因此软件可以利用它来做一些颜色渐变的动画，取色器工具等等。

HSI颜色转换为RGB颜色时，因为HSI色彩空间的精度是无限的取决于软件，而RGB是显示设备上的模型，设备制造成型它的色值范围就固定了。因此需要定义RGB的最大分量值和最小分量值，并计算色度值。

M = Max(R,G,B)
m = Min(R,G,B)
C = M-m

当 C=0，即意味是单色系统，不存在转换的可能，当M取值来自RGB的不同分量时，HSI的计算公式就可以按分段函数 piecewise 来表示：

注意公式H’的分段函数还没有完全罗列出来，如果R为最大分量，|G-B|<C，这里有两种情况 G>B 和 B>G，在软件的算法处理中就需要分别对待。如 G>B，那么+0，否则算式就要+1。极端情况下G=B，即是纯红色，色调为0°=0*60°，同理其它两条算式也是需要根据两个分量的大小进行处理，这里约定减号后面的分量较前面的大，算式就分别+1，+3，+5。算式中 mod6 结合下一步的60°就意味着色调值保持在 [0,360°) 这个范围。

在前面提到的几个类似的模型中，差别主要体现在亮度和饱和度度分量的计算方法，C=0时，饱和度均为0：

V = M
S = C/V  ;HSVL = (M+m)/2
S = C/(1-|2L-1|) ;HSL

HSI模型转换到RGB时，就是通过已经的HSI三个分量来找到RGB的对应分量。以HSL模型的算法为例。已知 S、H、L，可以求出M/m和色度C，还有两个小分量的差值X；通过H的度数信息可以得到H’，由前面的模型公式可知，H的度数还与 RGB 分量的大小相关，可以来用确定 X/M/m与RGB分量的对应关系。

C = S * (1-|2L-1|)
H'= H / 60°
m = L - C/2
M = L + C/2
X = C (1-|H' mod2-1|)H'=[0,2) >M=R，[2,4) >M=G，[4,6) >M=B
H'=[0,1) >m=G，[5,6) >m=B，[1,3) >m=B，[3,5) >m=R，[6,6) >m=G，

通过以上条件可以确定RGB中的其中2个分量，最一个分量就是 X+m，这个过程用公式表达比较费解，所以就用文字描述了。

YUV彩色模型

YUV这种彩色模型较前面几种特殊，单独一个小节进行介绍。一般的视频采集芯片输出的码流一般都是YUV数据流的形式，因此做视频采集与处理，必需要理解YUV模型和YUV数据流。

在早期的彩色电视系统中，还没有采用彩色像素数字显示系统，而是靠电子显像管的电子扫描系统成像的，因此不能直接使用RGB模型，需要采用模拟电子电路通用的彩色模型 YUV，配套的编码系统就是逐行倒相 Phase Alternating Line (PAL)制。与之相关有 Y’UV, YUV, YCbCr, YPbPr等等，好多人傻傻分不清。Wikipedia 上解析算比较完整的，Y’UV 和 YUV 两种叫法是过去只有模拟的显示系统时使用色彩模型，YPbPr 则是用于模拟分量信号视频系统的色彩模型。而到了现代，数字计算机普及，YUV 在各种图象软件上有广泛的应用，特别适用于静态图片的处理，于是就产生了 YPbPr 的数字版 YCbCr。

YUV模型定义了一个亮度分量 luma (Y’) 和两个色度分量 UV，分别称为蓝差信号和红差信号，这些分量并不某个词的缩写，只是它们就被这样用而已。Cb/Pb 和 Cr/Pr 也都是蓝差信号和红差信号，只是它们的算法有差异，在不歧义的情况下可以用UV代替它们。 Y’ 上的分号表示伽马修正量，不带伽马修正的亮度分量称为 luminance (Y)。在早期，电视等等图像系统是黑白色彩，只有一个亮度信号，YUV的出现主要是为了在彩色信号系统上兼容旧的黑白电视系统。黑白电视在接收到YUV信号时，只使用亮度分量，色度分量丢弃。分离亮度分量是YUV色彩模型的一个重要特点，因为人眼对黑白色调敏感，对色度分量反而不是那么敏感，在YUV采样过程就是利用这个原理对色度分量进行缩减采样以实现数据压缩。

和RGB模型类似地，YUV模型也可以按 12-bit, 16-bit 或 24-bit 来编码每个像素。常用的模式有以下四种，420 这样的数字表示各色度分量与亮度分量的采样比：

YUV444   12 bytes per 4 pixels
YUV422    8 bytes per 4 pixels
YUV411    6 bytes per 4 pixels
YUV420p   6 bytes per 4 pixels, reordered

关于YUV的采样缩减，YUV444 作为正常采样的情况，就是对每一个像素进行进行三个分量的采样。YUV422 表示在扫描线方向隔一像素采样色度分量，即2:1(4:2)的缩减色度采样，对于一行四个像素的采样过程来计，亮度采样有4个，而色度分量各为2个，亦即422数字的含义。YUV420这种也是扫描线方向2:1缩减采样，只是对于色度分量采用事隔行交换扫描。YUV411 则表示在扫描线方向每四个像素中采样4个亮度分量，色度分量各采样一个，实际使用中这样的缩减采样手段比RGB的采样更有效地节省存储空间，也有利节省带宽，对图像质量的影响亦可以接受。可以想像一块待扫描的图像区域，对于不同的YUV模式，扫描状态可以用下图表示，●实心圆表示亮度色度都采样，○空心圆表示只采样亮度，□表示只采样色度分量：

●●●●●●●●    ●○●○●○●○    ○○○○○○○○    ○○○●○○○●
●●●●●●●●    ●○●○●○●○    □  □  □  □    ○○○●○○○●
●●●●●●●●    ●○●○●○●○    ○○○○○○○○    ○○○●○○○●
●●●●●●●●    ●○●○●○●○    □  □  □  □   ○○○●○○○●YUV4:4:4       YUV4:2:2         YUM4:2:0        YUV4:1:1

由于YUV采样过程的特殊性，在编译像素分量值是，YUV三分量是分开处理的，并没有像RGB那样逐个像素存储。对一块图像区域来讲，前后两行的像素可能共享了相同的色度分量值，典型的YUM420就是这样的情况，编码流先送入亮度分量，再送入色度分量。

YUV和RGB的转换公式比较混乱，因为不同的系统差异，运算速度要求等等需求不同，公式有所调整因此产生各种版本。但最重要的是亮度分量这条公式，它是整个YUV的基础，有了亮度分量，那么用RGB的蓝色分量减去亮度就是蓝差色度分量，用红色分量减去亮度就得到红差色度分量。在电路系统中，RGB信号分量会被量化为[0,1]的范围，因此就会有以下一组标准公式：

    Y = 0.299 R + 0.587 G + 0.114 B          Y∈[0,1]
B - Y =-0.299 R - 0.587 G + 0.886 B = U      U∈[-0.886,0.886]
R - Y = 0.701 R - 0.587 G - 0.114 B = V      V∈[-0.701,0.701]

对于像我一样懒数学的人来讲，Mathematica 就是福音，像上面这样的方程，输入MM后只需要 Solve 一下就可以给出解的表达式：

ClearAll[eqs, Y, U, V, R, G, B]
eqs = {Y == 0.299 R + 0.587 G + 0.114 B,U == B - Y,V == R - Y}
Solve[eqs, {Y, U, V}]Out[50]= {{Y ->  0.114 B + 0.587 G + 0.299 R, U ->  0.886 B - 0.587 G - 0.299 R, V -> -0.114 B - 0.587 G + 0.701 R
}}

常用的代数操作有提公因式 Factor，展开式 Expand，消元 Eliminate，约分 Cancel，化简 Simplify，求实数解 NSolve 等等。

在1982年国际电信联盟-无线电通信部门ITU-R(International Telecommunication Union – Radiocommunication sector)发布的 Rec.601(Recommendation BT.601)规范中，使用 YCbCr 4:2:2 采样方案，规定了以下几个常数：

Wr=0.299    Wg=0.587    Wb=0.114    Umax=0.436  Vmax=0.615

在计算RGB信号的UV分量时，就使用的是缩减值域的公式，留余量以避免色值在转换时饱和失真：

U = (B-Y)*UMax/(1-Wb)≈0.492(B-Y)
V = (R-Y)*Vmax/(1-Wr)≈0.877(R-Y)

那么各分量的取值范围就是 Y/R/G/B∈[0,1]、U∈[-0.436,0.436]、V∈[-0.615,0.615]。那么相应的转换公式就成了另一个版本，这就是模拟系统 PAL 或 NTSC 用的标准的转换公式，符号上的分号’表示这是一个伽马修正量，通常取值[0，1]：

Y' =  0.114 B' + 0.587 G' + 0.299 R'
U' =  0.436 B' - 0.289 G' - 0.147 R'
V' = -0.100 B' - 0.515 G' + 0.615 R'
R' =  Y' + 1.140 V'
G' =  Y' - 0.395 U' - 0.581 V'
B' =  Y' + 2.033 U'

用行列式表达就是：

[Y']   [ 0.299  0.587  0.114] [R']
[U'] = [-0.147 -0.289  0.436]·[G']
[V']   [ 0.615 -0.515 -0.100] [B'][R']   [ 1.000      0  1.140] [Y']
[G'] = [ 1.000 -0.395 -0.581]·[U']
[B']   [ 1.000  2.033      0] [V']

在高清电视 HDTV High Definition TV 系统上，BT.709 标准则采用了不同的转换系数，规定 Wr=0.2126、Wb=0.0722，Wg=1-Wr-Wb=0.7152：

[Y']   [ 0.2126  0.7152  0.0722] [R']
[U'] = [-0.0999 -0.3361  0.4360]·[G']
[V']   [ 0.6150 -0.5586 -0.0564] [B'][R']   [ 1      0   1.2803] [Y']
[G'] = [ 1 -0.2148 -0.3806]·[U']
[B']   [ 1  2.1280       0] [V']

这些公式可以通过Mathematica进行计算得到：

ClearAll[eqs, Y, U, V, R, G, B, Wr, Wg, Wb]
Wr = 0.2126; Wb = 0.0722; Wg = 1 - Wr - Wb; Umax = 0.436; Vmax = 0.615;
eqs = {Y == Wr R + Wg G + Wb B,U == (B - Y) Umax/(1 - Wb),V == (R - Y) Vmax/(1 - Wr)};
Simplify[eqs]
Solve[eqs, {R, G, B}]

在数字系统中公式则需要数字化版本的 YCbCr，转换公式则如下，分量值域为R/G/B∈[0,255]，Y∈[16，235]，Cb/Cr∈[16,240]，值域留用余量兼顾到模拟系统，避免在转换过程出现色值饱和而产生图像失真：

[ Y]   [ 0.257  0.504  0.098] [R]   [ 16]
[Cb] = [-0.148 -0.291  0.439]·[G] + [128]
[Cr]   [ 0.439 -0.368 -0.071] [B]   [128][ R]   [ 1.164      0  1.596] [(Y -   16)]
[ G] = [ 1.164 -0.391 -0.813]·[(Cb - 128)]
[ B]   [ 1.164  2.018      0] [(Cr - 128)]

在 Keith Jack 的《Video Demystified》(ISBN 1-878707-09-4)书中则给RGB转换的另一种形式，只是重排了矩阵：

[ B]   [ 1.164      0  2.018] [( Y -  16)]
[ G] = [ 1.164 -0.813 -0.391]·[(Cr - 128)]
[ R]   [ 1.164  1.596      0] [(Cb - 128)]

在高清电视 HDTV High Definition TV 系统上，则采用了不同的转换系数：

[ Y]   [ 0.183  0.614  0.062] [R]   [ 16]
[Cb] = [-0.101 -0.339  0.439]·[G] + [128]
[Cr]   [ 0.439 -0.399 -0.040] [B]   [128][ R]   [ 1.164      0  1.793] [(Y -   16)]
[ G] = [ 1.164 -0.213 -0.533]·[(Cb - 128)]
[ B]   [ 1.164  2.112      0] [(Cr - 128)]

对于是纯数字图像处理系统，如JPG图片应用则可以不留余量，YCbCr和RGB全作用域转换，所有六个分量取全作用域[0,255]，因此重新定义了系数以避免各分量超出限值：

[ Y]   [ 0.299  0.587  0.114] [R]   [  0]
[Cb] = [-0.169 -0.331  0.500]·[G] + [128]
[Cr]   [ 0.500 -0.419 -0.081] [B]   [128][ R]   [ 1.000      0  1.400] [Y]
[ G] = [ 1.000 -0.343 -0.711]·[(Cb - 128)]
[ B]   [ 1.000  1.765      0] [(Cr - 128)]

还有 YPbPr 模型的转换，由于计算机使用较少，就不作介绍了，有兴趣可以参考后面的链接信息。

色彩分离

一般情况下图片处理会先将不同的色彩进行分类处理。这时就可以使用色彩分离函数 ColorSeparate，它可以任意指定以下的各种色彩空间进行色彩分离，默认模式为 RGB：

"Grayscale","Gray","Intensity"  灰度级强度
"RGB"                   红色、绿色、蓝色
"RGBA"                  红色、绿色、蓝色、alpha
"CMYK"                  青色、品红、黄色、黑色，工业印刷色彩
"Cyan","Magenta","Yellow","Black"   青色、品红、黄色、黑色
"HSB"                   色调、饱和度、亮度
"Red","Green","Blue"    红色、绿色或者蓝色通道
"Hue"                   色调
"Saturation"            饱和度
"Brightness"            亮度
"Mean"                  所有通道的均值，即黑白图
"Alpha"                 alpha 通道
"Luminance"             亮度通道
"XYZ"                   CIE XYZ 色彩空间的通道
"LAB"                   CIE LAB 色彩空间的通道
"LUV"                   CIE LUV 色彩空间的通道

不同的模式使用不同的计算公式，如常用的黑白灰度图可以用这条转换公式：

Y = 0.299*R + 0.587*G + 0.114*B

最近在着迷任天堂发迹的那个时代，为此就来一个超级玛丽做DEMO了：

{i, ColorSeparate[i, "RGB"]}

通过 ColorSeparate 分色后产生三个色调图，即灰度图，分量图片可以通过 ColorCombine 函数来合并。可以看到最后一幅图最明显，因为图片背景色是蓝色调，RGB中的蓝色成分的值 B=255，是最大值的分量值，形成分量图片时就显示为白色了，而蓝色分量值较小的区域就显示为较深的灰度色，或黑色如果 B=0。注意白云，在RGB图片中，白色即是三分量为 255 的颜色，生成分量图时就显示为白色。类似地还有黑色，因为黑色意味RGB三个分量均为 0 ，所以分量图上同样显示为黑色。

以上图片是经过组合的，手动编排为一个 2x2 的图片矩阵，再使用 ImageAssemble 合成为一张图，这样方便写作编辑。

ImageAssemble[%3]
Export["C:\writing\Mathica\mario.png", %4]

以下是各种模型的分量图，可以更直观感觉不同色彩模型的差别：

img = Import["ExampleData/peacock.tif"]
{Flatten[{"RGB ", img, ColorSeparate[img, "RGB"]}], Flatten[{"HSB ", img, ColorSeparate[img, "HSB"]}], Flatten[{"LAB ", img, ColorSeparate[img, "LAB"]}], Flatten[{"XYZ ", img, ColorSeparate[img, "XYZ"]}], Flatten[{"LUV ", img, ColorSeparate[img, "LUV"]}], Flatten[{"CMYK", img, ColorSeparate[img, "CMYK"]}]}
Export["C:\\writing\\Mathica\\peacock.jpg", %]

[ColorSeparation.jpg]

图片颜色的直方图统计通过 ImageHistogram 函数来计算，对于FC的图形来说可用的颜色不多，也就几十种，所以直方图显示出来的就是一条条分明的色块，这些色块代表了不同的色值在图片中出现的频度。右半部分是负片效果，通过 ColorNegate 函数即可处理得到。

ColorConvert 函数可以将图片转换到指定的色彩空间，而 ImageColorSpace 函数则可以读取图片所在的色彩空间。

一些工具与方法

在处理图片前学习一些有用的绘图函数，Graphics 可以用来绘制需要的图形，通过Graphics函数，可以将任意的图形对象绘到。比如说一个实心圆，扇形，椭圆，自定义的边线，Bézier 曲线等等：

{Graphics[{Hue[1], Disk[]}],Graphics[{Orange, Disk[{0, 0}, 1, {Pi/4, 3 Pi/4}]}],Graphics[{Orange, Disk[{0, 0}, {3, 4}]}],Graphics[{EdgeForm[{Thick, Black, Dashed}], Pink, Disk[]}]}
f[x_] := ImageResize[x, {196, 196}]
list = Map[f, %%]
ImageAssemble[%]
Export["C:/writing/Mathica/graphics.jpg", %]

在图形处理中常常会需要一些控件来修改参数，以方便查看参数对图形的影响。带控件的操作函数，Manipulate 就是为这个准备的。还可以使用定位器 Locator 来定位图片的像素，定点后可以用来在指定点进行绘图，或通过 ImageValue 函数读取像素的色值。比如，给正弦函数添加两个控件来控制其周期可以幅度：

Manipulate[Plot[Sin[a x + b], {x, 0, 6}], {a, 1, 4}, {b, 0, 10}]
Manipulate[Graphics[Line[{{0, 0}, p}], PlotRange -> 2], {{p, {1, 1}}, Locator}]

又或者使用取值范围来产生控件，并且控件的标题设置为“Period”，然后这个控件的值会反馈到 t 变量上。因为这个变量会被传入正弦函数，所以会影响函数输出的图形，下面的例子就是通过拖动控制点来改变曲线的同期：

Manipulate[Plot[{Line[{{0, 0}, p}], 2^3 Sin[x 2 π/p[[1]]]}, {x,0,π}], {{p, {0.1, 1}}, Locator}]

也可以产生拾色器控件，例中变量 c 和拾色器关系，当通过控件改变色值时变量 c 传入 Graphics 函数进行图形绘制：

Manipulate[Graphics[{c, Disk[]}], {{c, Blue, "Change to"}, Blue}]

以上都是 Manipulate 函数常用的功能，它可以看作是一个交互界面，通过它可以产生其它类型的控件，如二维滑尺，输入栏，选择框，甚至是弹出菜单也可以。

Manipulate 这种动态更新的能力其实是通过动态功能 Dynamic，DynamicModule 实现的，后者是前者的模块化版本，增加了局域变量的管理。动态的基本原理就是为需要动态的对象指定一个监视程序，当对象被改变时监视程序就更新输出，也可以通过监视程序来修改对象的内容。如例中两个同步的滑块，因为它们都绑定了同一个动态监视程序 Dynamic[x]，任何时候修改 x 值时，它们都会即时更新：

{Slider[Dynamic[x]], Slider[Dynamic[x]], Dynamic[x]}

但是，如果 Dynamic 的参数不是可以赋值的对象时，问题就来了，如 Dynamic[1 - x] 是不能和 Slider 这样的控制绑定的，因为动态对象不能赋值。这就要使用 Dynamic 的第二个参数了，它接收一个函数定义，通过它来改变动态执行的计算。

Slider[Dynamic[1 - y, (y = #) &]]

表达式中出现了纯函数的另一种定义，(x=1-#)&，Dynamic 调用这个函数时，参数就会被赋予 # 所表示的变量。这个函数实际上完成了，对变量 y 的赋值，然后根据 y 的值计算得到第一个变量表达式的值 1 -y，并返回给 Slider。因为动态具有这种实时更新的功能，所以下面的例句可以实现无限循环的自我更新，或有限的循环：

Dynamic[x = If[x == 0, 1, Max[0, x - 0.01]]]
Dynamic[x = Max[0, x - 0.01]]

在这个例子中，有一个很重要的事实被忽略了，那就是 Dynamic 的 HoldFirst 特性。简单来说，HoldFirst 就是先放着待会再处理的意思，Dynamic 是一个在前端运行的格式化函数，因此不能把它用在需要表达式的地方。如布下的例子就不会像期待一样正常运行，Plot 这个命令需要有具体的 x 的数值才能画出一个图形，但是在被画图的函数里的 Dynamic[x]。在内核中不会被计算成任何内容，它只保持为 Dynamic[x] 不变，使得 Plot 命令不能作出任何有意义的事情：

DynamicModule[{x}, {Slider[Dynamic[x]], Plot[Sin[Dynamic[x] i], {i, 0, π}]}]

由于 HoldFirst 特性，在正确的位置使用 Dynamic 就显得非常重要了，比较下面的两条相似的语句，那个才是正确的呢：

DynamicModule[{x = 0.5}, {Slider[Dynamic[x]], Dynamic[x]}]
DynamicModule[{x = 0.5}, {Dynamic[Slider[x]], Dynamic[x]}]

第二条是错误的用法，那是因为当用 Dynamic 封装 Slider[x] 时，在计算它的内容时，x 的值被代入，结果只是 Slider 的第一变量，一个具体数字，没有跟踪剩下的变量名，这种情况下的滑块是一个静态滑块的动态显示。对于控件，有一个简单的规则来决定在什么地方放 Dynamic。任何控件函数，例如 Slider（滑块）、Checkbox （复选框）或 PopupMenu（弹出菜单）的第一个自变量都几乎总是 Dynamic[x]。理解以上的区别，就可以使用 Dynamic 封装函数的方法动态地选择值域来进行绘图：

Manipulate[Dynamic[Plot[3 E ^(x^2), {x, -y, y}]], {y}, {{y, 1.6}, .1^10, 10}]

如有一组数据，需要通过动态来调整取值：

data = {.1, .5, .3, .9, .2};
Dynamic[data]
Table[Slider[Dynamic[data[[i]]]], {i, 5}]

，执行后你能看见在滑块周围有一个错误的指示，它们不能被移动，并且上方的动态输出从来不变化。因为 Dynamic 它有 HoldFirst 这个特性，不计算它的第一个自变量，所以出乎意料的是，这行不通！问题出现在，变量 i 被 Table 命令赋予了一个临时的数值，但是那个数值从来没有被使用，因为 Dynamic 是 HoldFirst 的。看一下那个滑块列表的 InputForm 就会明白这个问题：

In[59]:= Table[Slider[Dynamic[data[[i]]]], {i, 5}] // InputFormOut[59]//InputForm=
{Slider[Dynamic[data[[i]]]], Slider[Dynamic[data[[i]]]], Slider[Dynamic[data[[i]]]], Slider[Dynamic[data[[i]]]], Slider[Dynamic[data[[i]]]]}

所以根据这个输入内容可以采取其它方法来修正，比如用配合 With 就可以正常运行：

Table[With[{i = i}, Slider[Dynamic[data[[i]]]]], {i, 5}]

需要注意在等待 Dynamic 的输出计算的过程中，前端是冻结的，连打字动作都是不可能的。所以重点是把放在 Dynamic 里面的表达式限制为那些能相对迅速完成的计算，为了避免锁住前端，动态计算在内部以时间约束函数 TimeConstrained 封装，使用默认 5 秒钟超时数值，能用 DynamicEvaluationTimeout 选项来改变。在某种极端情况下，TimeConstrained 会失败以放弃计算，在这种情况下，几秒钟之后，前端会显示出一个对话框允许你终止动态更新直到引起问题的那个输出被删除为止。

如果确实需要在 Dynamic 里有计算慢的内容，那么可以在这种情况下禁止同步更新选项 SynchronousUpdating->False 就可以允许动态以一种不会锁住前端的方式来计算，当然这会使得计算相对没那么快。下面这个计算这个例子占时 10 秒钟，但由于禁止了同步更新，所以在 Dynmaic 运算的过程中还可以做其它事。

Dynamic[{DateList[], Pause[10]; DateList[]}, SynchronousUpdating -> False]

有了这些基础，就可以来模拟 Animate 的效果了，如下两段代码产生相似的动画：

Animate[Plot[Sin[x + y], {x, y, y + 10}, Axes -> {True, False}], {y, 0, 2 π}]Module[{t}, t = 0; Dynamic[t = Mod[t + 0.05, 2 π];Plot[Sin[x], {x, t, t + 8}, Axes -> {True, False}, ImageSize -> {480, 220}]
]]

为了导出图片的序列以生成 gif 动画，可以通过另一个动态来跟踪并导出图片帖：

SetDirectory[NotebookDirectory[]];
Module[{i}, i = 0;Dynamic[i = i + 1; Export[ToString[i] <> "_Trig.gif", %]]

为了将图片帖生组合生成gif动画，可以使用下面的代码片段，根据图片序列的数据调整 x y 参数：

SetDirectory[NotebookDirectory[]];
Table[ToString[4 x + y] <> "_Trig.gif", {x, 0, 18}, {y, 1, 4}]
Flatten[%]
Import[#] & /@ %
Partition[%, 4]
ImageAssemble[#] & /@ %
Export["trigonometry.gif", %]

上面这个gif动画中的动感线条是不很像鲤鱼跃龙门？确实有几分相似，这种相似并不是偶然的，其实三角函数是自然界普遍的存在。现在越来发越喜欢三角函数，因为它够讨人喜欢，更因为它真的好玩，生活的每个角落都是它的影子，就连走路时腿部的关节运动、手部的摆动都有它的存在。下面这段代码则是三角函数原理的动画展示:

Module[{x, y, t, s}, t = 0; s = 1.2;Column[{Slider[Dynamic[s], {1, 2}];Dynamic[Graphics[{{LightBlue, Rectangle[{-Dynamic[s], -Dynamic[s]}, {Dynamic[s], Dynamic[s]}]},Point[{1.2, y}], Text[Style["by Jimbowhy", 18], { 2.6, -1}],Text[Style["(x,y) Sin(θ)=y/r", 18], {1.1 x, 1.1 y}], Text[Style["r", 18], { 2 x/5, y/2}],Text[Style["θ", 18], {Cos[t/2]/4, Sin[t/2]/4}],Circle[{0, 0}, 0.3, {0, t}],{Thin, Dashed, Line[{{x, y}, {1.2, y}}]},Plot[Sin[-2 (z - 1.2) + t], {z, 1.2, 2 π + 1.2}][[1]],{Thin, Arrowheads[0.05], Arrow[{{0, 0}, {x, y}}]},{Thin, Dashed, Line[{{0, y}, {x, y}}]}, {Thin, Dashed, Line[{{x, 0}, {x, y}}]},{Thin, Line[{{-1, 0}, {3.8, 0}}]}, {Thin, Line[{{0, 1}, {0, -1}}]},Circle[{x = Cos[t], y = Sin[t = If[t >= 2 π, 0, t + 0.02]]},0.0]}, PlotRange -> {{-1.2, 4}, {-1.2, 1.2}}, ImageSize -> {480, 220}]]}]]

通过前面的代码可以将其处理后保存为gif动画，这张图就比较有意思了，三角函数的的可视化演示，早几年前刚开始玩编程的时候就想玩这个的了，专为数学渣准备的神器：

那么以上这么多功能如何以完成的操作界面呈现呢，那就要通过布局控件，如 Panel（面板）、TabView（选项卡）、Row（行）、Column（列）、Grid（格子）等等，以及其它的格式函数就可以组织好操作界面，Manipulate 就是这样做的。

另外还可以通过 SystemDialogInput 来调用系统提供的常用对话框，只需要参会字符串参数就可以打开如拾色器 Color、目录选择 Direcotry,文件选择 FileOpen/FileSave，甚至还可以是录音对话框 RecordSound。

接下来是关于图片处理的函数。前面提到的色彩分离，是将整个图片的色彩的重新归类处理，而对于每个像素，可以使用 ColorConvert 函数，将其转换到不同的色彩空间，得到的返回值就是对应色彩模型的各个分量值，如下例显示了同一个色在 HSB、LAB 色彩空间的分量值：

Manipulate[{ColorConvert[c,Hue],ColorConvert[c,"LAB"]}, {{c, Blue, "Pick"}, Blue}]

使用 ImageApply 函数来应用图片的像素计算，例如下面将图片的各像素各通道分量应用 Max/Min 函数进行计算，并将像素所有通道分量值替换为最大/最小的分量值，这个值由 Max/Min 函数计算后得出：

ImageApply[Max, Import["ExampleData/lena.mgf"]]
ImageApply[Min, Import["ExampleData/lena.mgf"]]Y[R_, G_, B_] := (0.299*R + 0.587*G + 0.114*B)
ImageApply[Y, ColorSeparate[img]]

ImageApply 函数接收的参数个与按通道数一样，而且通道值按顺序传入处理函数的参数列表中。

在图片处理的技术中，有一种叫伪色增强 Pseudo Coloring 的概念。医院使用的电子扫描成像设备如磁共振成像 MRI(Nucler Magnetic Resonance Imaging)，电子计算机断层扫描 CT(Computed Tomography)，公共场所安检使用的X射线扫描设备，还有电子显微镜等等，它们是利用X射线、γ射线、超声波对物体进行探测，形成只有强弱而无色彩的信号反馈，为了将这些结果呈现出来，就需要将信号的强弱关系映射了眼睛可以看见的色彩空间上，这个过程就是色彩化，也就是伪色增强的概念。MM 中有一个函数叫 Colorize，就是用来色彩化的函数，系统定义了多套色彩化参数方案，可以通过 ColorData 来获取，它可以返回色彩方案的 ColorDataFunction 函数，表示色值的映射关系，也可以返回色彩方案的名称列表 Name，也可以返回色彩方案的图像等等。

比如下例使用的红外线测温度的使用的热力映射方案 TemperatrueMap，ThermometerColors：

fiber = Import["http://www.forth.com/images/sidebar/image_x1.jpg"];
Colorize[fiber, ColorFunction -> ColorData["TemperatureMap"]]
Colorize[fiber, ColorFunction -> ColorData["ThermometerColors"]]

各个色彩方案也可以以列表的形式列出来：

Manipulate[map, {{map, "SunsetColors", "Color Map"},ColorData["Gradients"]}]

以下代码则示范如何查询系统定义的色彩方案的名称，以及四个色彩方案的色卡图片。Gradients、Physical、Indexed 和 Named 是四个系统定义方案集合，它们包含各种色彩方案的定义：

ColorData["Gradients", "Name"]
ColorData["Physical", "Name"]
ColorData["Named", "Name"]ColorData["WebSafe", "Image"]
ColorData["HTML", "Image"]
ColorData["TemperatureMap", "Image"]
ColorData["CMYKColors", "Image"]

图片处理范例

现在来综合一下前面讲的这些功能，将它们组合起实现一个可以调整图片色调的程序。程序需要一个取色器，用来选取目标色调。还需要 Locator 来定位像素，通过 ImageValue 取得像素的色值后将其转换到 HSB 色彩空间，在 ImageApply 运算时用来比较源图的像素是否属与定位点上的像素同色调。运算过程还可以加入一个误差值，增加误差可以扩大匹配的图像范围，减小误差则缩小匹配的图像范围。这里为何选择 HSB 色彩模型作为运算的方式呢？正如前面提到的，HSB 这个模型更能反映人眼的成像模式，通过比较图像的色调就更能模拟眼睛对色彩的感官效应。源图就选用前面用到的超级玛丽：

img = Import["C:\\writing\\VirtualNes\\MARIO.jpg"];
{hue, sat, brt} = ColorSeparate[img, "HSB"];
Htr[from_, to_, tol_] := {ImageApply[Function[{hraw}, If[Abs[from - hraw] > tol, from, to]], hue], sat,brt};
Manipulate[ColorCombine[Htr[ImageValue[hue, p], ColorConvert[cto, Hue][[1]], tol], "HSB"], {{cto, Blue, "Color"}, Red}, {{p, {125, 45}}, Locator}, {{tol, 0.1, "Tolerance"}, 0, 0.5}]
Export["C:\\writing\\Mathica\\colorized_mario e.gif", %]

需要注意的是，HSB 这个色彩模型并不能完全反映眼睛的成像模型，它只是在一定程度上模拟而已。所以上面的程序不见得会很完美，但是如果图片的辨识度够高，效果还是很不错的，换个肤色什么的那是不在话下。就效果上来看，和Photoshop的 Hue/Saturation 图像调整功相同，因为算法都是同一套。不够完美的地方就是，明明看着图片上是着色分明的一块，针对它调整着色时，却发现其它区域，通常是浅色区域的变化更大。而Photoshop的 Selective Color 选色调整功能使用的是 CMYK 色彩模型，效果显然是专业级的，比 HSL 模型的调整要精准，但调整的强度要弱小一点，这就是三分量与四分量的差别。

就 HSL 模型下，如何做到更精准的色调修改呢？前面讲的色调分离就可以更好地利用起来，可以查看 HSB 分享通道中，B 通道是对比最明显的，也可以试试其它色彩模型的通道分离，如 CMYK 的 Y 即黄色通道也是，RGB 模型中的 B 通道，XYZ 的 Z 通道的对比度都是最高的。这些通道用起来，精准度是也会是最高的。因此可以使用其它色彩模型高对比度的通道来作为蒙版，当图像的像素在蒙版区域时，就当作是目标区域进行处理，否则忽略它，这样就有针对性了。

I HAVE THE POWER & TO BE CONTINUED…….

资源参考

HSI Color Space
Color vision
YUV Color Space
Showcasing Manipulate via .GIF animations
Recommended 8-Bit YUV Formats for Video Rendering (Windows)
Using SIMD Instructions For Image Processing
YUV Color Conversion - equasys GmbH
Linear Algebra and Its Applications Fourth Edition by Gilbert Strang
Mathematica Second Edition by Eugene Don, Ph.D. Queens College, CUNY