半色调印刷品中光传播过程的蒙特卡罗模拟

刘洪豪,王　琪,柴江松,于艺铭

(南京林业大学轻工科学与工程学院,南京 210037)



半色调印刷品中光传播过程的蒙特卡罗模拟

刘洪豪,王琪*,柴江松,于艺铭

(南京林业大学轻工科学与工程学院,南京 210037)

摘要:基于蒙特卡罗方法思想,结合半色调印刷品网点的结构属性,提出一种研究半色调印刷品呈色性质的模拟方法。考虑光的粒子属性,以程序化语言描述光子在半色调印品中的传播历程,涵盖光子从入射到传播终止的全过程,体现承印物、油墨、网点形状以及面积率变化对光散射行为的影响。对175 lpi的不同网点结构模拟结果表明：该算法能合理模拟光与印品相互作用,适用于光子在多种网点结构中的散射模拟。模拟后的光子能量及位置分布均能可视化显示,经统计分析可转化为反射率、光子分布,散射距离等信息,为半色调印品的微观研究及质量控制提供一种新思路。

关键词：半色调印品；光传播；蒙特卡罗模拟

1引言

作为信息载体,印刷品是人们获取信息的重要来源。印刷品是一个复杂的层状结构物质,光在印刷品中的传播形式较为多样,包含有镜面反射、镜面透射、界面内反射、界面外反射、光散射等[1]。现有的对网点反射率相关研究各有侧重点,有重点考虑油墨与纸张结合效应的,更多的是重点描述光传播现象。Murray-Davies模型[2]是最早的反射率模型,其理想化认为反射率只与网点部分反射率及纸张自身反射率相关,忽略了光散射现象的影响,模型精度不高。Yule-Nielsen模型[3]考虑光散射效应的影响,通过引入修正参数n,对Murray-Davies模型进行了修正,一定程度上提高反射率预测精度,但修正参数n是经验值,具有一定随机性。在Yule-Nielsen模型基础上,出现Clapper-Yule模型[4],其考虑了墨层与纸张间散射和界面反射现象,反射率精度明显提高。之后,根据Kubelka-Munk理论[5],研究光线在介质内多重反射产生的影响,Kubelka-Munk 预测模型[6]也被建立。此外,还有将光散射概率[7]、网点扩大函数[8]作为侧重点研究的模型。从上述模型不断发展改进过程中不难发现,影响反射率模型精度主要因素是光散射、反射及多重折射,由于光传播行为的复杂性及不可描述性,上述研究成果对光传播过程进行了适当的简化和假设,未能够将半色调印品呈色的各个因素全面纳入考虑。

蒙特卡罗法是一种以概率统计为基础的随机模拟方法,能够较好实现对光传播过程的定量分析。早期S. A. Prahl和 M. Keijzer等人详细介绍光在介质中传播的蒙特卡罗模拟流程[9],而后,Lihong Wang和S. L. Jacques[10]提出完善的光在多层介质中传播的模拟方法。本文以此为基础,结合半色调印刷品网点的结构属性及光的粒子属性,提出针对半色调印品中光传播过程的蒙特卡罗算法流程,为研究半色调印品中光传播行为提供新思路。

2模拟方法原理

2.1光传播过程模拟

光照射到半色调印品上,通过与油墨、纸张等介质的相互作用,最终分解为携带不同能量的电磁波,因此,模拟光子与介质的整个作用过程可分为光子的生成、介质中传播、光子的终止三个环节。光子的生成是对每个与介质发生相互作用的光子进行初始化；模拟光子在介质中的传播过程主要是对光子的移动步长和散射距离进行随机抽样及分析；光子传播的最终结果有反射、透射、吸收三种形式。

(1)光子初始化

光子初始化是模拟光源生成光子的过程,主要包括光子束能量、入射位置和方向的初始化。在模拟过程中,认为光子入射方向垂直于印品,光子的入射位置为(x,y,z),传播方向使用(ux,uy,uz)表示,其中ux,uy,uz分别为传播方向与x,y,z坐标轴的夹角余弦值。光子束能量以概率密度形式表示,初始化能量人为的规定为w个单位,光子每次在介质中移动都会伴随着吸收与散射,经过一次移动后,光子束能量一部分被吸收,另一部分散射出去进入下一次模拟流程,经过一次散射后部分能量为w′,可表示为

(1)

Fig.1　Illustration of coordinates and angle

入射位置不同会影响光子束在纸张中的能量差异。印品表面分为着墨与空白部分,光子入射位置在着墨部分时,光子需经历油墨层的吸收过程。当油墨层透射率为ti时,经油墨层吸收后,光子束能量为wti。

(2)光子散射移动步长

光子移动过程中主要需考虑传播步长和散射方向的随机抽样。光子在介质中传播步长s的概率密度函数p(s)可表示为[11]

p(s)=σte-σts

(2)

其中σt=σα+σs,σt为作用系数,σα为介质吸收系数,σs为介质散射系数。由此抽样可得

(3)

(3)光子散射方向

α=2πε′

(4)

偏转角β的余弦值符合Henyey-Greenstein函数的概率分布[9]

(5)

式(5)中g为各向异性系数,其范围为[-1,1]。当g为0,表示各向同性散射,g取值越接近-1,表示越接近完全背向散射,g取值越接近1,表示越接近完全前向散射。偏转角β的抽样可通过式(6)实现,其中ε为新生成的随机数,ε∈[0,1]

(6)

当方位角α和偏转角β确定之后,光子散射方向(ux′,uy′,uz′)满足以下关系

(7)

(4)光子传播终止

光子在介质中传播最终结果有吸收、反射、透射三种形式,被称为光子的终止。每个入射的初始化光子,经过在介质中的不同传播,有且仅得到一种结果。

光子束中部分光子在传播过程中会被吸收、出射,表现为光子束能量的损失,经有限次散射传播后光子束能量趋近于零,但当光子束能量极小时,即使光子最终离开介质发生反射或透射,其作用效果也相当有限。为避免计算机无穷次计算,提高模拟效率,提出轮盘筛选规则[10],规定当光子束能量小于一定阈值w″时,该光子束有1/m的机会继续传播,其能量变为mw″,否则该光子束完全被吸收,能量变为零。轮盘筛选规则通过预先设定m的数值,比如10,和生成[0,1]之间的随机数ε实现

(8)

光子在介质中每一次移动步长s都有可能撞击到介质与空气的界面,可能离开介质,出射到空气中,也有可能发生内反射,回到介质中。光子在界面处是否发生内反射取决于光子移动步长以及移动方向与界面的夹角。要判断光子从原位置(x,y,z),沿方向(ux,uy,uz),移动步长s至新位置(x′,y′,z′)时,是否会撞击到界面,可先判断|z′|与pt/2的大小关系,如果|z′|≤pt/2,光子新位置在介质内,未撞击到界面；如果|z′|>pt/2,光子若沿方向(ux,uy,uz)移动会撞击到界面,但不能确定是出射还是内反射。假设ni和nt分别为光子所在介质和出射介质的折射率,当发生全反射时,光子移动方向与界面的撞击夹角αi大于全反射临界角αc,此时光子发生内反射,其中

αi=cos-1(|uz|)

(9)

αc=sin-1(nt/ni)

(10)

其他情况下,光子在界面是否发生内反射由Fresnel反射系数R(αi)确定

(11)

根据Snell准则[12],出射角αt与撞击夹角αi存在以下关系

nisinαi=ntsinαt

(12)

取随机数ε∈[0,1]与R(αi)比较决定光子是出射还是内反射。当ε≤R(αi)时,光子发生内反射；当ε>R(αi)时,光子发生出射,离开介质。若光子发生内反射,可将光子视为移动至界面处,其在界面处的坐标点可通过求经过点(x,y,z)且方向余弦为(ux,uy,uz)的直线与界面的焦点得到。光子对应的移动方向相应地变为(ux,uy,-uz)。若光子出射,离开介质,当光子由上界面出射,则被视为反射光子,需要判断其是否由着墨部分出射,着墨部分出射的光子束能量需经墨层吸收；当光子由下界面出射,则被视为透射光子。

2.2不同网点结构模拟

以同心圆网点为例：同心圆网点类似圆形调幅网点,是一种双圆环结构,其中着墨部分称为“条”,空白部分称为“空”。

Fig.2　Geometry of round and concentric dot

当加网线数为L,网点面积率为a时,圆形网点的半径可由式(13)求得

(13)

同心圆网点重要的特点是具有固定的条空参数比。其着墨部分称为条,网点中空白部分称为空。条空参数比是指加网中构成条空的像素宽度比[8]。设条的宽度为t,空的宽度为k,条空之比为m/n,结合图2理解,构成着墨部分的双圆环宽度都为t,而构成空白部分包括一个圆形和一个圆环,空白圆形半径和空白圆环宽度为k。在加网线数为L,网点面积率为a的情况下求解条空宽度如下

(14)

为方便区别着墨网点和纸张空白部分,在构建坐标系时,可以加网单元网格中心点作为XY坐标系的原点。光子在圆形和圆环着墨网点传播时,都可通过比较入射或者出射点到坐标原点距离和半径大小,从而判断光子是在着墨网点部分或是纸张空白部分。假设光子入射或出射点与原点距离为d,在圆形网点半径为r的情况下：当d≤r时,光子出现在着墨网点区域,油墨层对该光子会有吸收作用；当d>r时,光子出现在纸张空白区域,光子能量不会衰减。

以同样的方法,可分别针对圆形、菱形、椭圆、方形网点进行模拟。

3实验结果及分析

3.1算法实现

根据蒙特卡罗模拟流程,在Matlab中针对不同网点结构进行模拟算法编程。所编写的程序具有较强的适用性、扩展性,可对各项印刷参数进行设定,包括加网线数、网点形状、网点面积率、纸张属性、油墨属性等[13-14]。

过程模拟是算法实现的核心内容,就本文而言就是对光子传播过程的程序化表达,光子传播过程中的各种情况都需要在模拟过程中得到体现。具体而言,半色调印刷品光子传播过程模拟的特殊方面为：光子入射和反射位置点是在承印物还是油墨层,两者对光子的作用是存在差异的。此外,光子在承印介质中的移动步长、散射方向的改变,光子束能量的吸收,阈值和轮盘筛选规则以及光子撞击界面时的判断都是编程的难点。光子每移动一个步长就要对以上所对应的各个情况进行判定,判定结果决定不同的传播历程。算法程序可区别三种传播结果,分别记录光子最终所在的位置和对应的光子束能量,对数据进行分析统计可模拟出反射率、透射率、网点扩大值等信息。

3.2半色调印品网点模拟结果

本文中所编写的模拟算法实用性较高,在不同加网线数下,可实现对不同网点面积率的不同网点结构进行模拟。同样以同心圆网点为例：同心圆网点加网线数设为175 lpi,条空比为1.8/1.3[14],以10%网点面积率间隔,选取10%~80%网点进行模拟,如图3。

由图3可发现：在模拟图中,同心圆网点的形状完整,随着网点面积率的变化,网点形状发生相应变化。总体上,该模拟算法能够直观体现网点微观结构。

此外,不同的网点结构也可模拟再现,以20%网点面积率为例,不同网点结构模拟结果如图4。

3.3算法反射率的验证

在选定了适当的光子数和验证了模拟结果稳定性的基础上,本文设计针对圆形、同心圆网点结构反射率的验证实验。在加网线数为175 lpi的条件下,对圆形网点以10%为间隔,0到100%的网点面积率进行了模拟；考虑到同心圆的特殊结构,仅对同心圆网点为10%~80%部分进行模拟,实验所用同心圆样张由ESKO公司推出。模拟反射率是由模拟结果数据分析得到。测量反射率值由Murray-Davies公式先求得密度值,根据密度与反射率的关系式,转化为对应的反射率。反射率验证结果如图5。

Fig.3　Simulation of concentric dot

Fig.4　Simulation of different dot structures

Fig.5　Verification of round dot

Fig.6　Verification of concentric dot

全阶调范围内,圆形网点模拟网点反射率与测量网点反射率数据较吻合,表明调幅圆形网点的蒙特卡罗模型精确度较高。图6表明：在10%~60%阶调范围,同心圆网点模拟算法较合理,但当网点面积率为60%以上时,模拟反射率与实际测量反射率出现较大偏差,这主要是由于175 lpi同心圆网点在网点面积率为60%时出现网点搭接情况,网点面积率越高,搭接情况越严重,导致数据偏差越大。

3.4光子分布

本文所设计的算法程序为了能够了解光子终止的结果,会将反射、吸收和透射三种类型光子的最终位置分类记录,为分析光子运行情况提供数据基础。实验中以加网线数为175 lpi,面积率为50%的圆形网点作为模拟对象,模拟结果分别以散点形式标注在XZ二维平面和XYZ三维空间内,对应得到圆形网点的二维和三维图为a和b。在图4中“o”代表反射光子位置,“*”代表被吸收光子位置,“△”代表透射光子位置。需要说明的是图3中只体现了1000个光子传播的结果,如果将所有10万个光子都标注出来,图示中将遍布散点,根本无法体现其中的规律,所以选择1000个光子作为代表。

Fig.7　Simulated photons distribution of 50% round dot(a,b)

由于在判定光子离开纸张时,如果光子从上界面出射,会被认为是反射光子,其反射位置点会集中在纸张的上界面,反之,如果光子从下界面出射,则会被认为是透射光子,透射位置点就会集中在纸张的下界面。被吸收的光子散布在整个介质三维空间中。由图中可直观发现,三类光子的大部分都较为集中且有规律性,但也存在少数异点,会分布于更为分散的范围内。同一网点面积率下,三类光子的分布情况存有差异。随着网点面积率的变化,三类光子各自所占的比例也会发生改变。根据对不同网点面积率下圆形网点光子最终结果的统计,得到反映吸收、反射、透射光子的比例关系,见表1。

Tab.1　The proportion of absorption,reflection and transmission

由表1中的数据变化可发现：随着网点面积率的增加,吸收和透射光子比例会增加,反射光子比例会降低。网点面积率越大,入射光子由着墨区域入射的概率越大,进入纸张内部进行传播的光子能量会降低,在与纸张发生作用的过程中,光子能量小于阈值的概率就越大,相应地,吸收光子的比例会有所增加。纸张表面的直接反射是反射光子重要的来源,网点面积率的增加使发生纸张表面直接反射的光子减少,从而导致反射光子比例降低。直接镜面反射光子减少,进入纸张内部的光子数则会增加,发生透射的光子会随之增多。透射光子比例的增加并非意味着透射率的增加,透射率不仅与透射光子的多少有关,还与透射光子能量大小有关。未直接反射的光子会被油墨层吸收,能量会有所损失,再经过多次的传播过程,最终达到纸张下界面的光子能量会偏低。以上数据仅对应实验中特定的条件才成立,不同纸张和油墨组合,三类光子之间的比例关系并非固定。表中数据所体现的是三类光子随着网点面积率变化所发生的对应变化。

4总结

本文提出了基于蒙特卡罗模拟的研究半色调印品中光传播行为的一种方法,该模拟算法能合理体现光与印品的相互作用,有助于印品的微观质量控制。光传播过程模拟包含光与承印物和油墨相互作用的多种形式,基于蒙特卡罗思想,结合网点形状和网点面积率的变化,能够实现对半色调印品中光传播行为进行较为完整的模拟,并从反射率方面对模拟算法进行验证,验证结果表明蒙特卡罗模拟方法精确度较高。虽然本文仅粗略地从反射率方面对175 lpi圆形、同心圆网点的模拟算法进行验证,但至少为印品的微观研究及质量控制提供新思路。将来可以不断加以完善,如可从透射、吸收及位置能量分布方面探究,也可以基于此模拟法探究不同网点结构类型及纸张不同参数对印品的影响。

参考文献：

[1]Jang W,Allebach J P. Characterization of printer MTF[J]. J Imaging Sci Technol,2006,50(3):264-275.

[2]Hersch R D,Collaud F,Cretel F,etal. Spectral Prediction and dot surface estimation models for halftone prints[C]. Proceedings of the Society of Photo-optical Instrumentation Engineers,2004:356-369.

[3]Ruckdeschel F R,Hauser O G. Yule-Nielsen effect in printing:a physical analysis[J]. Applied Optics. 1978,17(21):3376-3383.

[4]Viggiano J A. Ink Penetration,Isomorphic Colorant Mixing and Negative Values of Yule-Nielsen n[C]. Society for Imaging Science and Technology 2010(1): 285-290.

[5]Kubelka P,Munk F. Ein beitrag zur optik der farbanstriche[J].F Tecn Physik,1931,12 :593-601.

[6]施国运,董娜,张逸新. 双面半色调印刷品Kubelka-Munk色彩预测模型[J]. 光散射学报,2009,01:13-17(Shi Guoyun,Dong na,Zhang yixin. The color prediction for halftone recto-verso print based on Kubelka-Munk Model[J]. J Light Scatt,2009,01:13-17)

[7]刘萧豪,王琪. 基于光散射概率函数的半色调反射率模型[J]. 包装工程,2014,01:108-113. (Liu Xiaohao,Wang qi. Halftone printing reflectance model based on light scattering probability function [J]. Packaging Engineering,2014,01:108-113)

[8]葛惊寰,刘振. 基于光学网点扩大的光谱反射率模型研究[J]. 光散射学报,2011,02:109-113.(Ge Jinghuan,Liu zhen. Research on the spectral reflectance model based on optical dot gain[J]. J Light Scatt,2011,02:109-113)

[9]Wang L,Jacques S L,Zheng L. MCML-Monte carlo modeling of light transport in multi-layered tissues[J].Comput. Methods Programs Biomed. 1995,47(2): 131-146.

[10]Damir M,Stanislav B. Monte carlo modeling of light scattering in paper[J]. J Imaging Sci Technol.2009,53(2):020201-020208.

[11]Li T,Gong H,Luo Q M.MCVM-Monte carlo modeling of photon migration in voxelized meda[J]. J Innov Opt Health Sci.2010,3(2):91-102.

[12]Prahl S. A,Keijzer M,Jacques S L,etal.A monte carlo model of light propagation in tissue[J].Dosimetry of Laser Radiation in Medicine and Biology 1989,5:102-111.

[13]Bestem’yanow KP,Gordienko VM. Optical heterodyning study of the propagation dynamics of IR femtosecond laser pulses in a strongly scattering porous medium[J] Quantum Electron. 2004,34(7): 666-668.

[14]王琪. 基于网点结构形态的图像信息印刷复制研究[D],南京林业大学,2013.(Wang Qi.Study on image printing based on dot shape[D]. Nanjing Forestry University,2013)

Monte Carlo Simulation of Light Propagation in Halftone Prints

LIU Hong-hao,WANG Qi*,CHAI Jiang-song,YU Yi-ming

(SchoolofLightIndustryTechnologyandEngineering,NanjingForestryUniversity,Nanjing210037,China)

Abstract:The simulation of light propagation in halftone prints within the Monte Carlo method is provided,considering the structural characteristics of halftone prints. This simulation model describes the whole route of photons propagation in mathematical language,illustrating a real-scattering effect,which is applicable to various dot structures. The 175 lpi simulation results showed that: the model can reasonably simulate the interaction of light with printed materials,which was suit to various dots. Simulated distribution of energy and location was visible. Statistical analysis could be converted to reflectance,scattering distance information,which provides a new idea to research and control for microscopic halftone prints.

Key words:halftone printing; light propagation; Monte Carlo simulation

中图分类号：TS801.1

文献标志码：A

doi:10.13883/j.issn1004-5929.201601001

作者简介：刘洪豪(1992-)，硕士，主要研究：印刷品质量控制模拟模型。Email：mobeikehan@126.com通讯作者：王琪(1971-)，副教授，主要研究：色彩管理与控制技术、印刷防伪技术等。E-mail：wangqi_3639@163.com

收稿日期：2015-05-20; 修改稿日期: 2015-06-18

文章编号：1004-5929(2016)01-0070-07