CN101458822A - 一种3d模型的计算全息图快速生成方法 - Google Patents

Abstract

本发明公开了一种3D模型的计算全息图快速生成方法，在世界坐标系中的三角面片经过平移和一系列的旋转后，变换到“标准位置”，同时，对三角面片实施的变换，同样在全息面上进行，使全息面与物面中的三角面片的相对位置在三角面片坐标变换后保持不变。通过菲涅耳衍射积分公式推导出全息面上任意一点光场与三角面片顶点光场间的关系，然后将特定三角面片顶点坐标和光场代入公式求出全息面上各采样点光场。本发明极大地节省了计算量、避免了物面的采样误差且可生成任意分辨率的计算全息图。

Description

一种3D模型的计算全息图快速生成方法

技术领域

本发明涉及利用全息技术取得光学图像的方法，具体是指一种3D模型的计算全息图的快速生成方法。

背景技术

目前，大部分三维物体的显示是通过投影到二维平面上实现的，丢失了景深信息，观察者只能通过遮挡关系等线索感知三维物体。三维立体显示技术的出现使真实再现三维场景成为可能。

通过给左、右眼呈现两幅不同的视差图像，可以提供视觉感知的深度信息，达到三维显示的目的。双色立体图、立体眼镜、偏光镜等都是采用“视差图像对”原理的立体显示技术，但都需要佩戴眼镜，使用不便；同时视差图像对生成时，对焦距等参数的选择也有一些限制，某些场景的三维表现不理想。

根据双眼的视角不同，多视图立体显示器通过双眼不同的观察角度获取具有视差的左右两幅图像而获得立体视觉感受。但它对观测区域有特定要求，头部的转动可能出现重影等问题。因此，该技术需要头部跟踪系统，相关技术复杂；并且该系统只适合单人收看。

全息图能同时记录物光波的强度和相位信息，可以逼真地再现物体的三维像，其应用范围覆盖了商品虚拟展示、广告、娱乐、教育及艺术等领域，显示了巨大的市场潜力。计算全息图的产生不需要复杂和精密的相干光设备，并且可以产生现实世界中并不存在的物体的全息图，且能非常容易地进行参数的调整。通过相干光照射，可以在空间再现原场景的光场分布，因此，其逼真度很高，观察者不需佩戴附属观察设备，可供多人观察。

在二维的物平面上，肉眼可观察尺寸物体的计算全息图需要接近光波波长量级采样，计算全息图要分别在物平面上和全息面上采样后计算产生。由于光波的波长相对肉眼可观察的尺寸要小很多，因此采样点数量很大，从而导致计算量非常巨大，需要耗费很长的时间。

很多软件工具可以创作3D模型，它可以仿真现实场景，在广告、游戏等领域已经广泛应用。计算全息图的优势之一是不需要实际场景，而3D模型作为一种虚拟场景，其计算全息图的生成技术目前较少且发展不成熟。通常的方法是将3D模型分层，按二维的方法分别计算全息图，然后叠加在一起获得三维物体的全息图。这样将进一步增加计算量。

综上所述，3D模型的计算全息图具有很大的实用价值，但现有的技术没有利用3D模型本身的构造特点来设计快速算法。其技术瓶颈在于计算量非常巨大，需要耗费的时间很长；同时，计算量巨大也相应地对硬件提出了更高的要求，导致计算设备的投资增加。

发明内容

本发明在于克服现有技术的缺点和不足，提供一种3D模型的计算全息图的快速生成方法。3D模型由三角面片组成，从其顶点控制其外观表现出发，通过推导建立了顶点参数与计算全息图之间的联系，得到一组生成计算全息图的公式。只要将3D模型各顶点的数据代入公式，就可以求出计算全息图，从而节省了计算量、适用性强、精度高、设备要求低，且可生成任意分辨率的计算全息图。

本发明的目的是通过以下方案实施的：一种3D模型的计算全息图快速生成方法，包括以下步骤，如图9所示：

(1)根据需要确定全息面的位置与大小，选择全息面上的各采样点的位置，即得到各采样点的坐标，则各采样点及其坐标表示为S(x_s，y_s，z_s)，置此时各采样点S(x_s，y_s，z_s)的光场初值U_s(x_s，y_s，z_s)为0；

(2)由于3D模型由三角面片组成，其顶点控制其外观表现，标记3D模型中通过“背面消隐”后保留的三角面片为0，表示其没被处理过。

使用“背面消隐”技术，可以确定3D模型中哪些三角面片的光场可以投射到全息面上。在3D模型中，根据各三角面片的顶点的排列顺序和3D模型的规定，可以确定三角面片的正面。一般3D模型中规定，通过三角面片的任意两边对应的矢量叉乘后得到其法线。该法线与全息面法线的点乘除以它们的模可得两法线夹角的余弦，如果该余弦非正，则背对全息面，标记为1，其它标记为0；

(3)选择任意标记为0的三角面片，标记为1，在世界坐标系中，设标记为1的三角面片的三个顶点坐标分别是A_w(x_wa，y_wa，z_wa)、B_w(x_wb，y_wb，z_wb)和C_w(x_wc，y_wc，z_wc)，如图1，将该标记为1的三角面片△W中最长的边变换到x轴上，将最长边定为A_wC_w，这样将保持在第一象限和第二象限分别有一个子区域。计算空间任意两点(x₁，y₁，z₁)和(x₂，y₂，z₂)距离的公式为：

$d=\sqrt{{(x_1-x_2)}^2+{(y_1-y_2)}^2+{(z_1-z_2)}^2}---\left(1\right)$

对任意的三角形面片，根据其三个顶点的坐标，可以分别计算三条边的长度，这样可以确定边A_wC_w。

接下来要确定点B_W至边A_wC_w的垂足V(x_v，y_v，z_v)，作为移动到原点的基点。将图1中各点看作从原点发出的向量终点，则直线A_wC_w可以表示为：

$\stackrel{\→}{\mathrm{OV}}=\stackrel{\→}{{\mathrm{OC}}_w}+t\·\stackrel{\→}{n}---\left(2\right)$

其中

表示直线A_wC_w方向向量

$\stackrel{\→}{n}={\stackrel{\→}{\mathrm{OA}}}_w-\stackrel{\→}{{\mathrm{OC}}_w}---\left(3\right)$

t表示一个(-∞，+∞)范围内的实数。由垂直关系，可知：

$\stackrel{\→}{{\mathrm{VB}}_w}\·\stackrel{\→}{n}=0$

即：

$(\stackrel{\→}{\mathrm{OV}}-{\stackrel{\→}{\mathrm{OB}}}_w)\·\stackrel{\→}{n}=0---\left(4\right)$

将式(2)与式(3)代入上式，有：

$[\stackrel{\→}{{\mathrm{OC}}_w}+t\·(\stackrel{\→}{{\mathrm{OA}}_w}-\stackrel{\→}{{\mathrm{OC}}_w})-\stackrel{\→}{{\mathrm{OB}}_w}]\·({\stackrel{\→}{\mathrm{OA}}}_w-\stackrel{\→}{{\mathrm{OC}}_w})=0$

根据点A_w、B_w、C_w的坐标，得到各向量分量，代入上式后：

$[x_{\mathrm{wc}}+t\·(x_{\mathrm{wa}}-x_{\mathrm{wc}})-x_{\mathrm{wb}},y_{\mathrm{wc}}+t\·(y_{\mathrm{wa}}-y_{\mathrm{wc}})-y_{\mathrm{wb}},z_{\mathrm{wc}}+t\·(z_{\mathrm{wa}}-z_{\mathrm{wc}})-z_{\mathrm{wb}}]\·\left[\begin{array}{c}{x}_{\mathrm{wa}}-x_{\mathrm{wc}}\\ y_{\mathrm{wa}}-y_{\mathrm{wc}}\\ z_{\mathrm{wa}}-z_{\mathrm{wc}}\end{array}\right]=0$

即：

(x_wc-x_wb)(x_wa-x_wc)+t·(x_wa-x_wc)²+(y_wc-y_wb)(y_wa-y_wc)+t·(y_wa-y_wc)²

+(z_wc-z_wb)(z_wa-z_wc)+t·(z_wa-z_wc)²＝0

整理后得：

$t=\frac{(x_{\mathrm{wb}}-x_{\mathrm{wc}})(x_{\mathrm{wa}}-x_{\mathrm{wc}})+(y_{\mathrm{wb}}-y_{\mathrm{wc}})(y_{\mathrm{wa}}-y_{\mathrm{wc}})+(z_{\mathrm{wb}}-z_{\mathrm{wc}})(z_{\mathrm{wa}}-z_{\mathrm{wc}})}{{(x_{\mathrm{wa}}-x_{\mathrm{wc}})}^2+{(y_{\mathrm{wa}}-y_{\mathrm{wc}})}^2+{(z_{\mathrm{wa}}-z_{\mathrm{wc}})}^2}---\left(5\right)$

这样就求出了t，将其代入式(2)可确定垂足V的坐标(x_v，y_v，z_v)，即：

$\left[\begin{array}{c}{x}_v\\ y_v\\ z_v\end{array}\right]=\left[\begin{array}{c}{x}_{\mathrm{wc}}+t\·(x_{\mathrm{wa}}-x_{\mathrm{wc}})\\ y_{\mathrm{wc}}+t\·(y_{\mathrm{wa}}-y_{\mathrm{wc}})\\ z_{\mathrm{wc}}+t\·(z_{\mathrm{wa}}-z_{\mathrm{wc}})\end{array}\right]---\left(6\right)$

然后把垂足V移动到世界坐标系的原点0，同时对应地将整个三角面片△W进行移动，△W变换成△W₁，对应的坐标为：

$\left\{\begin{array}{c}{x}_{1a}=x_{\mathrm{wa}}-x_v,y_{1a}=y_{\mathrm{wa}}-y_v,z_{1a}=z_{\mathrm{wa}}-z_v\\ x_{1b}=x_{\mathrm{wb}}-x_v,y_{1b}=y_{\mathrm{wb}}-y_v,z_{1b}=z_{\mathrm{wb}}-z_v\\ x_{1c}=x_{\mathrm{wc}}-x_v,y_{1c}=y_{\mathrm{wc}}-y_v,z_{1c}=z_{\mathrm{wc}}-z_v\end{array}\right.---\left(7\right)$

则该标记为1的三角面片平移后三个顶点的坐标分别是A₁(x_1a，y_1a，z_1a)、B₁(x_1b，y_1b，z_1b)和C₁(x_1c，y_1c，z_1c)，如图2；

(4)对步骤(3)中平移后标记为1的三角面片△W₁进行旋转：

(4.1)将△W₁绕z轴旋转α，使边A₁C₁旋转至x-z平面内，△W₁变换成△W₂，绕z轴旋转的矩阵为：

$M_1=\left[\begin{array}{ccc}\cos \mathrm{\α}& \sin \mathrm{\α}& 0\\ -\sin \mathrm{\α}& \cos \mathrm{\α}& 0\\ 0& 0& 1\end{array}\right]---\left(8\right)$

其中 $\cos \mathrm{\α}=\frac{x_{1a}}{\sqrt{x_{1a}^2+y_{1a}^2}},$   $\sin \mathrm{\α}=\frac{y_{1a}}{\sqrt{x_{1a}^2+y_{1a}^2}}$

旋转后三角面片△W₂如图3所示，△W₂中各顶点坐标为：

$\left[\begin{array}{c}{x}_{2a}\\ y_{2a}\\ z_{2a}\end{array}\right]=M_1\·\left[\begin{array}{c}{x}_{1a}\\ y_{1a}\\ z_{1a}\end{array}\right],$   $\left[\begin{array}{c}{x}_{2b}\\ y_{2b}\\ z_{2b}\end{array}\right]=M_1\·\left[\begin{array}{c}{x}_{1b}\\ y_{1b}\\ z_{1b}\end{array}\right],$   $\left[\begin{array}{c}{x}_{2c}\\ y_{2c}\\ z_{2c}\end{array}\right]=M_1\·\left[\begin{array}{c}{x}_{1c}\\ y_{1c}\\ z_{1c}\end{array}\right]---\left(9\right)$

(4.2)然后三角面片△W₂绕y轴旋转β，使边A₂C₂与x轴重合，△W₂变换成△W₃，绕y轴旋转的矩阵为：

$M_2=\left[\begin{array}{ccc}\cos \mathrm{\β}& 0& \sin \mathrm{\β}\\ 0& 1& 0\\ -\sin \mathrm{\β}& 0& \cos \mathrm{\β}\end{array}\right]---\left(10\right)$

其中 $\cos \mathrm{\β}=\frac{x_{2a}}{\sqrt{x_{2a}^2+z_{2a}^2}},$   $\sin \mathrm{\β}=\frac{z_{2a}}{\sqrt{x_{2a}^2+z_{2a}^2}}$

旋转后三角面片△W₃如图4所示，旋转后△W₃中各顶点坐标为：

$\left[\begin{array}{c}{x}_{3a}\\ y_{3a}\\ z_{3a}\end{array}\right]=M_2\·\left[\begin{array}{c}{x}_{2a}\\ y_{2a}\\ z_{2a}\end{array}\right],$   $\left[\begin{array}{c}{x}_{3b}\\ y_{3b}\\ z_{3b}\end{array}\right]=M_2\·\left[\begin{array}{c}{x}_{2b}\\ y_{2b}\\ z_{2b}\end{array}\right],$   $\left[\begin{array}{c}{x}_{3c}\\ y_{3c}\\ z_{3c}\end{array}\right]=M_2\·\left[\begin{array}{c}{x}_{2c}\\ y_{2c}\\ z_{2c}\end{array}\right]---\left(11\right)$

(4.3)然后三角面片△W₃绕x轴旋转γ至x-y平面内，△W₃变换成△W₄，绕x轴旋转的矩阵为：

$M_3=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos \mathrm{\γ}& \sin \mathrm{\γ}\\ 0& -\sin \mathrm{\γ}& \cos \mathrm{\γ}\end{array}\right]---\left(12\right)$

其中 $\cos \mathrm{\γ}=\frac{y_{3b}}{\sqrt{y_{3b}^2+z_{3b}^2}},$   $\sin \mathrm{\γ}=\frac{z_{3b}}{\sqrt{y_{3b}^2+z_{3b}^2}}$

旋转后三角面片△ABC如图5所示，旋转后其各顶点坐标为：

$\left[\begin{array}{c}{x}_a\\ y_a\\ z_a\end{array}\right]=M_3\·\left[\begin{array}{c}{x}_{3a}\\ y_{3a}\\ z_{3a}\end{array}\right],$   $\left[\begin{array}{c}{x}_b\\ y_b\\ z_b\end{array}\right]=M_3\·\left[\begin{array}{c}{x}_{3b}\\ y_{3b}\\ z_{3b}\end{array}\right],$   $\left[\begin{array}{c}{x}_c\\ y_c\\ z_c\end{array}\right]=M_3\·\left[\begin{array}{c}{x}_{3c}\\ y_{3c}\\ z_{3c}\end{array}\right]---\left(13\right)$

则该标记为1的三角面片旋转后三个顶点的坐标分别是A(x_a，y_a，z_a)、B(x_b，y_b，z_b)和C(x_c，y_c，z_c)。这样，最后三角形面片便旋转到了x-y平面内，并以y轴为分界线分割成两个子区域，这里把图5所示三角面片在坐标系中的位置，称为“标准位置”，在该位置，整个三角面片在x-y平面内，因此没有指出z轴。

式(9)、式(11)、式(13)可合并为简洁的形式：

$\left[\begin{array}{c}{x}_a\\ y_a\\ z_a\end{array}\right]{=M}_3\·M_2\·M_1\·\left[\begin{array}{c}{x}_{1a}\\ y_{1a}\\ z_{1a}\end{array}\right],$   $\left[\begin{array}{c}{x}_b\\ y_b\\ z_b\end{array}\right]{=M}_3\·M_2\·M_1\·\left[\begin{array}{c}{x}_{1b}\\ y_{1b}\\ z_{1b}\end{array}\right],$   $\left[\begin{array}{c}{x}_c\\ y_c\\ z_c\end{array}\right]{=M}_3\·M_2\·M_1\·\left[\begin{array}{c}{x}_{1c}\\ y_{1c}\\ z_{1c}\end{array}\right]---\left(14\right)$

世界坐标系中的任意三角面片经步骤(3)和步骤(4)的调整最终都可以达到“标准位置”。

(5)根据步骤(3)的V(x_v，y_v，z_v)和步骤(4)中的M₁、M₂、M₃，将全息面进行平移和旋转，使全息面与该标记为1的三角形面片的相对位置，在三角面片坐标变换后保持不变，具体的方法是：

要计算的是全息面上各采样点的光场，因此，对全息面位置的变换，体现在对各采样点S(x_s，y_s，z_s)位置的变换，该变换与该标记为1的三角形面片的变换完全相同，以保持全息面与该标记为1的三角形面片的相对位置不变。设变换后的采样点及其坐标为S′(u，v，η)，则通过下式求变换后采样点的位置：

$\left[\begin{array}{c}u\\ v\\ \mathrm{\η}\end{array}\right]{=M}_3\·M_2\·M_1\·\left[\begin{array}{c}{x}_s{-x}_v\\ y_s-y_v\\ z_s-z_v\end{array}\right]---\left(15\right)$

(6)计算全息面上各采样点S′(u，v，η)的光场U(u，v，η)，然后累加到全息面上各采样点原来的光场值U_s(x_s，y_s，z_s)上：

ΔABC可看作一个衍射孔径，全息面可看作衍射屏，如图6。根据菲涅耳衍射积分公式：

$U(u,v,\mathrm{\η})=\frac{e^{\mathrm{jk\η}}}{\mathrm{j\λ\η}}{e}^{j\frac{k}{2\mathrm{\η}}(u^2+v^2)}\underset{\mathrm{\Δ}}{\∫\∫}\left\{U(x,y)e^{j\frac{k}{2\mathrm{\η}}(x^2+y^2)}\right\}{e}^{-j\frac{2\mathrm{\π}}{\mathrm{\λ\η}}(\mathrm{ux}+\mathrm{vy})}\mathrm{dxdy}---\left(16\right)$

其中j²＝-1，λ为波长，k为波数，定义为k＝2π/λ，e为自然常数，U(x，y)是x-y平面内三角面片上的光场的分布。

三角面片上的相位分布不影响强度分布，而人眼对相位不敏感，因此想从计算全息图重建的是强度分布。可以让物光场U(x，y)的相位为

这样式(16)简化为：

$U(u,v,\mathrm{\η})=\frac{e^{\mathrm{jk\η}}}{\mathrm{j\λ\η}}{e}^{j\frac{k}{2\mathrm{\η}}(u^2+v^2)}\underset{\mathrm{\Δ}}{\∫\∫}\left|U(x,y)\right|e^{-j\frac{2\mathrm{\π}}{\mathrm{\λ\η}}(\mathrm{ux}+\mathrm{vy})}\mathrm{dxdy}---\left(17\right)$

|U(x，y)|表示一个0相位的光场，可以通过三角面片的顶点光场插值获得。由于接下来的分析不涉及U(x，y)的相位，因此以下不加区别地使用符号|U|与U。

任意一个三角型面片都可以旋转到图5的状态，它实际上由两个三角形组成，所以式(17)可以变形为：

$U(u,v,\mathrm{\η})=\frac{e^{\mathrm{jk\η}}}{\mathrm{j\λ\η}}{e}^{j\frac{k}{2\mathrm{\η}}(u^2+v^2)}[\underset{\mathrm{\Δ}1}{\∫\∫}{U}_1(x,y)e^{-j\frac{2\mathrm{\π}}{\mathrm{\λ\η}}(\mathrm{ux}+\mathrm{vy})}\mathrm{dxdy}+\underset{\mathrm{\Δ}2}{\∫\∫}{U}_2(x,y)e^{-j\frac{2\mathrm{\π}}{\mathrm{\λ\η}}(\mathrm{ux}+\mathrm{vy})}\mathrm{dxdy}---\left(18\right)$

式(18)中U₁(x，y)和U₂(x，y)分别代表三角形△1和△2上的光场。

进一步地，设：

$U_1(u,v,\mathrm{\η})=\underset{\mathrm{\Δ}1}{\∫\∫}{U}_1(x,y)e^{-j\frac{2\mathrm{\π}}{\mathrm{\λ\η}}(\mathrm{ux}+\mathrm{vy})}\mathrm{dxdy}---\left(19\right)$

$U_2(u,v,\mathrm{\η})=\underset{\mathrm{\Δ}2}{\∫\∫}{U}_2(x,y)e^{-j\frac{2\mathrm{\π}}{\mathrm{\λ\η}}(\mathrm{ux}+\mathrm{vy})}\mathrm{dxdy}---\left(20\right)$

则式(18)可表示为：

$U(u,v,\mathrm{\η})=\frac{e^{\mathrm{jk\η}}}{\mathrm{j\λ\η}}{e}^{j\frac{k}{2\mathrm{\η}}(u^2+v^2)}[U_1(u,v,\mathrm{\η})+U_2(u,v,\mathrm{\η})]---\left(21\right)$

如图7所示，已知三角形各顶点的坐标及光场，通过公式可推导出三角型中任一点的坐标与光场U₁(x，y)、U₂(x，y)的关系。这里点A、B、C的坐标分别是(x_a，0)、(0，y_b)、(x_c，0)，对应光场为U_a、U_b、U_c，U₀是原点处光场，在对△1和△2的光场插值时使用：

$U_o=\frac{1}{x_a-x_c}(U_c{x}_a-U_a{x}_c)---\left(22\right)$

按y方向插值得到：

$U_r=\frac{1}{y_b}[U_{b}y+U_o(y_b-y)]---\left(23\right)$

$U_p=\frac{1}{y_b}[U_{b}y+U_c(y_b-y)]---\left(24\right)$

$U_q=\frac{1}{y_b}[U_{b}y+U_a(y_b-y)]---\left(25\right)$

光场U_r是为插值引入的中间变量，与其对应的点在y轴上，其y坐标与插值点的y坐标相等，U_p、U_q对应三角形边上的点P、Q处光场。

根据三角形的比例关系，点P、Q的x坐标分别为：

$x_q=x_a(1-\frac{y}{y_b})---\left(26\right)$

$x_p=x_c(1-\frac{y}{y_b})---\left(27\right)$

按x方向插值得到两个三角形的光场：

$U_1=\frac{1}{x_q}[U_{q}x+U_r(x_q-x)]---\left(28\right)$

$U_2=\frac{1}{x_p}[U_{p}x+U_r(x_p-x)]---\left(29\right)$

为了计算式(19)和式(20)，将二重积分转化成两个一重积分，先按x方向积分，再按y方向积分：

根据式(28)可将式(19)转化为：

$U_1(u,v,\mathrm{\η})=\underset{0}{\overset{y_b}{\∫}}\left\{\underset{0}{\overset{x_a(1-y/y_b)}{\∫}}\frac{1}{x_q}[U_{q}x+U_r(x_q-x)]e^{-j\frac{2\mathrm{\π}}{\mathrm{\λ\η}}\mathrm{ux}}\mathrm{dx}\right\}{e}^{-j\frac{2\mathrm{\π}}{\mathrm{\λ\η}}\mathrm{vy}}\mathrm{dy}---\left(30\right)$

根据式(29)可将式(20)转化为：

$U_2(u,v,\mathrm{\η})=\underset{0}{\overset{y_b}{\∫}}\left\{\underset{x_c(1-y/y_b)}{\overset{0}{\∫}}\frac{1}{x_p}[U_{p}x+U_r(x_p-x)]e^{-j\frac{2\mathrm{\π}}{\mathrm{\λ\η}}\mathrm{ux}}\mathrm{dx}\right\}{e}^{-j\frac{2\mathrm{\π}}{\mathrm{\λ\η}}\mathrm{vy}}\mathrm{dy}---\left(31\right)$

(a)当u≠0且v≠0时，根据公式(30)、(31)求定积分后，将结果带入式(21)可得：

$U(u,v,\mathrm{\η})=\frac{1}{8}\·j\·\mathrm{\λ}\·y_b\·\mathrm{\η}\·e^{\frac{\mathrm{jk}}{2\mathrm{\η}}(u^2+v^2+2{\mathrm{\η}}^2)}\·\frac{f_1}{f_2}---(32-1)$

将f₁表示成u、v与η的三元多项式的形式，多项式的项分别为u⁴、u³·v、u³·η、u²·v²、u²·v·η、u·v³、u·v²·η与v³·η：

$f_1=(f_{11}\·r_1^{u/z}+f_{12}\·r_2^{v/z}+f_{13}\·r_3^{u/z})\·u^4+(f_{14}\·r_1^{u/z}+f_{15}\·r_2^{v/z}+f_{16}\·r_3^{u/z})\·u^3\·v+$

$(f_{17}\·r_1^{u/z}+f_{18}\·r_2^{v/z}+f_{19}\·r_3^{u/z})\·u^3\·\mathrm{\η}+(f_{31}\·r_1^{u/z}+f_{32}\·r_2^{v/z}+f_{33}\·r_3^{u/z})\·u^2\·v^2$

$+(f_{34}\·r_1^{u/z}+f_{35}\·r_2^{v/z}+f_{36}\·r_3^{u/z})\·u^2\·v\·\mathrm{\η}+(f_{37}\·r_1^{u/z}+f_{38}\·r_3^{v/z})\·u\·v^3+$

$(f_{41}\·r_1^{u/z}+f_{42}\·r_3^{u/z})\·u\·v^2\·\mathrm{\η}+(f_{43}\·r_1^{u/z}+f_{44}\·r_3^{u/z})\·v^3\·\mathrm{\η}$

f₁的多项式系数计算时，需要先计算其系数r₁、r₂、r₃，其计算公式为：

$r_1=e^{-2\mathrm{j\π}{x}_c/\mathrm{\λ}},$    $r_2=e^{-2\mathrm{j\π}{y}_b/\mathrm{\λ}},$    $r_3=e^{-2\mathrm{j\π}{x}_a/\mathrm{\λ}}$

计算多项式项u⁴的系数时，需要先计算其系数f₁₁、f₁₂、f₁₃，其计算公式为：

$f_{11}=2\mathrm{\π}\·(x_a^2\·x_c^2-x_a^3\·x_c)\·U_c,$    $f_{12}=2\mathrm{\π}\·(x_a^2\·x_c-2x_a^2\·x_c^2{+x}_c^3\·x_a)\·U_b$

$f_{13}=2\mathrm{\π}\·(x_a^2\·x_c^2-x_c^3\·x_a)\·U_a$

计算多项式项u³·v的系数时，需要先计算其系数f₁₄、f₁₅、f₁₆，其计算公式为：

$f_{14}=2\mathrm{\π}\·y_b\·(x_a^3-2x_a\·x_c^2{+x}_a^2\·x_c)\·U_c,$    $f_{15}=2\mathrm{\π}\·y_b\·(x_c^2\·x_a-x_c^3-x_a^3+x_a^2\·x_c)\·U_b,$

$f_{16}=2\mathrm{\π}\·y_b\·(x_c^2\·x_a-2x_a^2+x_a^2)\·U_a$

计算多项式项u³·η的系数时，需要先计算其系数f₁₇、f₁₈、f₁₉，其计算公式为：

$f_{17}=j\·\mathrm{\λ}\·[(x_a^2\·x_c-x_a^3)\·(U_b-U_c)+x_a^2\·x_c\·U_a]$

$f_{18}=j\·\mathrm{\λ}\·[(2x_c^2\·x_a-x_a^2\·x_c-x_c^3)\·U_a+(x_a^3-x_a^2\·x_c-x_c^2\·x_a+x_c^3)\·U_b$

     $+(2x_a^2\·x_c-x_c^2\·x_a-x_a^3)\·U_c]$

$f_{19}=j\·\mathrm{\λ}\·\left[(x_c^3-2x_c^2\·x_a)\·U_a+(x_c^2\·x_a-x_c^3)\·U_b+x_a\·x_c^2\·U_c\right]$

计算多项式项u²·v²的系数时，需要先计算其系数f₃₁、f₃₂、f₃₃，其计算公式为：

$f_{31}=2\mathrm{\π}\·y_b^2\·(x_c^2-2x_a^2{+x}_a\·x_c)\·U_c,$    $f_{32}=2\mathrm{\π}\·y_b^2\·{(x_a-x_c)}^2\·U_b$

$f_{33}=2\mathrm{\π}\·y_b^2\·(x_a^2-2x_c^2{+x}_a\·x_c)\·U_a$

计算多项式项u²·v·η的系数时，需要先计算其系数f₃₄、f₃₅、f₃₆，其计算公式为：

$f_{34}=j\·\mathrm{\λ}\·y_b\·[2(x_a^2-x_a\·x_c)\·U_b-(x_a^2+2x_a{x}_c)\·U_a+({4x}_a\·x_c-x_a^2)\·U_c]$

f₃₅＝j·λ·y_b·(x_a-x_c)²·(U_a-2U_b+U_c)

$f_{36}=j\·\mathrm{\λ}\·y_b\·[({4x}_a\·x_c-x_c^2)\·U_a+2(x_c^2-x_a\·x_c)\·U_b-(x_c^2+2x_a\·x_c)\·U_c]$

计算多项式项u·v³的系数时，需要先计算其系数f₃₇、f₃₈，其计算公式为：

$f_{37}=2\mathrm{\π}\·y_b^3\·(x_a-x_c)\·U_c,$    $f_{38}=2\mathrm{\π}\·y_b^3\·(x_c-x_a)\·U_a$

计算多项式项u·v²·η的系数时，需要先计算其系数f₄₁、f₄₂，其计算公式为：

$f_{41}=j\·\mathrm{\λ}\·y_b^2\·[({2x}_a+x_c)\·U_a+(x_c-x_a)\·U_b-(x_a+2x_c)\·U_c],f_{42}=-f_{41}$

计算多项式项v³·η的系数时，需要先计算其系数f₄₃、f₄₄，其计算公式为：

$f_{43}=j\·\mathrm{\λ}\·y_b^3\·(U_c-U_a),$     f₄₄＝-f₄₃

将f₂表示成u与v的二元多项式的形式：

f₂＝f₂₁·u⁶+f₂₂·u⁵·v+f₂₃·u⁴·v²+f₂₄·u³·v³+f₂₅·u²·v⁴

多项式项u⁶的系数f₂₁用下式计算：

$f_{21}={\mathrm{\π}}^3\·x_c^2\·(x_a^3-x_a^2\·x_c)$

多项式项u⁵·v的系数f₂₂用下式计算：

$f_{22}={2\mathrm{\π}}^3\·y_b\·(x_c^3\·x_a-x_c\·x_a^3)$

多项式项u⁴·v²的系数f₂₃用下式计算：

$f_{23}={\mathrm{\π}}^3\·y_b^2\·(x_a-x_c)\·(x_c^2+4x_a\·x_c+x_a^2)$

多项式项u³·v³的系数f₂₄用下式计算：

$f_{24}={2\mathrm{\π}}^3\·y_b^3\·(x_c^2-x_a^2)$

多项式项u²·v⁴的系数f₂₅用下式计算：

$f_{25}={\mathrm{\π}}^3\·y_b^4\·(x_a-x_c)$

作为公式(32-1)中的分子f₂要求不为0，这里限定u≠0且v≠0来满足该条件。

公式(32-1)在使用时，需要引用上面计算系数r₁、r₂、r₃、f₁₁、f₁₂、f₁₃、f₁₄、f₁₅、f₁₆、f₁₇、f₁₈、f₁₉、f₂₁、f₂₂、f₂₃、f₂₄、f₂₅、f₃₁、f₃₂、f₃₃、f₃₄、f₃₅、f₃₆、f₃₇、f₃₈、f₄₁、f₄₂、f₄₃、f₄₄的一组公式，称这一组公式为与公式(32-1)相关联的公式。对于某个三角面片，在计算所有采样点光场时，系数r₁、r₂、r₃、f₁₁、f₁₂、f₁₃、f₁₄、f₁₅、f₁₆、f₁₇、f₁₈、f₁₉、f₂₁、f₂₂、f₂₃、f₂₄、f₂₅、f₃₁、f₃₂、f₃₃、f₃₄、f₃₅、f₃₆、f₃₇、f₃₈、f₄₁、f₄₂、f₄₃、f₄₄这些值只需要计算一次，因为它们的计算式中并不包含采样点的位置信息，它们只与三角面片本身相关，而与采样点的坐标无关。可见，引入这些系数符号，可以节约计算开销。

(b)当u＝0且v≠0时，根据公式(30)、(31)求定积分后可得下列更简洁的形式：

$U(u,v,\mathrm{\η})=\frac{1}{8}\·e^{\frac{\mathrm{jk}}{2\mathrm{\η}}(u^2+v^2+2{\mathrm{\η}}^2)}\·[f_3\·(r_4-1)\·z^2\·v^{-3}+(f_4-f_5\·r_4)\·z\·v^{-2}+f_6\·v^{-1}]---(32-2)$

其中，r₄、f₃、f₄、f₅、f₆及下面要使用的r5这些系数是为了方便计算引入的符号，它们的计算公式如下：

$r_4=e^{-2\mathrm{j\π}{\mathrm{vy}}_b/\mathrm{\λ}/\mathrm{\η}},$    $r_5=2j\·\mathrm{\λ}\·y_b^{-1}\·{\mathrm{\π}}^{-2}\·(x_c-x_a)$

$f_3=(U_a+U_c-2U_b)\·(x_c-x_a)\·{\mathrm{\λ}}^2\·y_b^{-2}\·{\mathrm{\π}}^{-3}$

f₄＝r₅·(U_a+U_c-U_b)，f₅＝r₅·U_b，f₆＝2·(U_a+U_c)·(x_c-x_a)·π^-1

公式(32-2)在使用时，需要引用上面的一组公式，称这一组公式为与公式(32-2)相关联的公式。对于某个三角面片，在计算所有采样点光场时，r₅、f₃、f₄、f₅、f₆这些值只需要计算一次，因为它们也只与三角面片本身相关。

(c)当u＝0且v＝0时，根据公式(30)、(31)求定积分后可得下列更简洁的形式：

$U(u,v,\mathrm{\η})=\frac{y_b\·(x_a-x_c)\·(U_a+U_b+U_c)}{6j\·\mathrm{\λ}\·\mathrm{\η}}\·e^{\frac{\mathrm{jk}}{2\mathrm{\η}}(u^2+v^2+{2\mathrm{\η}}^2)}---(32-3)$

在上述(a)、(b)、(c)各情形中，j²＝-1，λ为波长，k为波数，定义为k＝2π/λ，e为自然常数，U_a、U_b、U_c是x-y平面内三角面片三个顶点A、B、C的光场；上述公式(32-1)、公式(32-2)、公式(32-3)是在不同条件下用不同的形式计算同一个量，统称为公式(32)，与公式(32-1)、公式(32-2)相关联的公式，统称为与公式(32)相关联的公式。若要求取衍射到空间任一点的光场，只要给出该点的坐标(u，v，η)，以及三角面片顶点的坐标与顶点处的光场，然后利用以上的一组公式，就可以求出其光场U(u，v，η)。最后将U(u，v，η)对应累加到全息面上各采样点的原光场值U_s(x_s，y_s，z_s)上。

(7)判断是否有标记为0的三角面片，若有，转步骤(3)；若无，则求各采样点S(x_s，y_s，z_s)处参考光的光场，叠加到U_s(x_s，y_s，z_s)上，产生计算全息图。

为更好的实现本发明，步骤(7)中所述参考光为球面光时，设参考点光源所在位置为(x_r，y_r，z_r)，发出的参考光在(x_s，y_s，z_s)处的光场为：

$R(x_s,y_s,z_s)=R_0(1/\sqrt{{(z_s-z_r)}^2+{(x_s-x_r)}^2+{(y_s-y_r)}^2})e^{\mathrm{jk}\sqrt{{(z_s-z_r)}^2+{(x_s-x_r)}^2+{(y_s-y_r)}^2}}---\left(33\right)$

其中，R(x_s，y_s，z_s)表示坐标为(x_s，y_s，z_s)的采样点的参考光光场，R₀表示与所选光源相关的复常数，其模为点光源发出的初始球面光波振幅，相位为点光源发出的初始球面光波的相位。

步骤(7)中所述参考光为平面光时，也需根据平面光波的特点计算该参考光在(x_s，y_s，z_s)处的光场。从光源发出平面光波，到达距离τ处的分布R_τ为：

$R_{\mathrm{\τ}}=R_p{e}^{\mathrm{jk\τ}}---\left(34\right)$

其中R_p表示与所选平面光源相关的复常数，其模为光源位置的平面光波的振幅，相位为光源位置的平面光波的初始相位。设平面光波的单位方向向量为

朝向全息面，从坐标原点0到采样点S的向量为

则从S(x_s，y_s，z_s)到光源所在平面的垂直距离为

与

的点积，将该距离作为τ代入公式(34)，就得到了在点(x_s，y_s，z_s)处的光场R(x_s，y_s，z_s)：

$R(x_s,y_s,z_s)=R_p{e}^{\mathrm{jk}(x_s{x}_p+y_s{y}_p+z_s{z}_p)}---\left(35\right)$

本发明的作用原理是：相比较于将3D模型分层，然后分别在物面和全息面上采样的传统方法，本发明的特征是不对3D模型分层与采样，而是逐面片分别计算其全息面光场，然后将它们叠加在一起而得到计算全息图。由于3D模型由三角面片组成，其顶点控制其外观表现。每个三角面片可看成一个衍射孔径，然后可以利用菲涅耳衍射积分公式计算采样点的光场。由于三角面片内部的点光场由高洛德着色模式，通过双线性插值得到，因此，在该孔径上的光场分布可以表示为三个顶点的解析式，这样，全息面光场的分布仅取决于三角面片的三个顶点，本发明通过菲涅耳衍射积分公式推导出了全息面上任意一点光场与三角面片顶点光场间的关系，当计算全息面上采样点光场时，只要将特定三角面片顶点坐标和光场代入式(32)即可求出结果。但由于式(32)的推导是假设三角面片在“标准位置”的前提上进行的，因此，在世界坐标系中的三角面片需要变换到“标准位置”才可直接适用，同时，对三角面片实施的变换，同样在全息面上进行，使全息面与3D模型中三角形面片的相对位置，在坐标变换后保持不变。

本发明与现有技术相比，具有如下优点和有益效果：

(1)极大地节省了计算量，用本发明生成3D模型的计算全息图，可以以0(1)的时间复杂度计算空间任意点的光场，从而极大地提高了计算速度。设3D模型由M片三角面片组成，用传统方法要分别在3D模型的物平面上和全息面上采样后计算生成，设物平面上和全息面上的采样点数量都为N，全息图上每一个采样点和3D模型上的所有点相联系，则要计算N*N次光场的叠加，时间为0(N*N)；而用本发明，只要将特定三角面片顶点坐标和光场代入公式，就可以求出该面片在采样点的光场，时间为0(M)，各面片的光场叠加后就得到了采样点的光场，时间为0(M*N)。由于M<<N，因此，新方法极大地节省了计算量。

(2)避免了物面的采样误差，从而有利于3D模型的精确重建。现有技术需要将3D模型分层，然后在各层次上按照二维物面全息图的方法计算。在二维的物面上，需要根据物面的实际大小与光波波长对其采样，然后用离散化的方法计算。采样点的光场表示受到表示方法的限制，必然有量化误差；在离散化的运算过程中，中间的计算步骤也必然伴随着舍入误差。而新方法不需要在物面采样，避免了量化误差；根据新方法的一组公式，可从3D模型的顶点信息直接计算全息面上各采样点的光强，因此计算过程是基于“连续”的方法，而不是“离散”的方法，避免了中间计算步骤的舍入误差。

(3)可渐进地生成任意分辨率的计算全息图。在现有技术中，当需要保持3D场景中物体的细节时，需要较高的物面采样率；反之可以用较低采样率。如果当前生成了低采样率的计算全息图，当需要获取高采样率的计算全息图时，物面的重采样导致低采样率下的计算结果难于利用。由于新方法不需要物面采样，可以先生成低分辨率的计算全息图，然后再计算更多的全息面采样点的光场，加入到低分辨率的计算全息图中，从而生成高分辨率的计算全息图。

附图说明

图1是世界坐标中的三角面片示意图；

图2是三角面片移动到坐标系原点的示意图；

图3是三角面片绕Z轴旋转至X-Z的示意图；

图4是三角面片绕Y轴旋转至A₂C₂与X轴重合的示意图；

图5是三角面片绕X轴旋转至X-Y平面的示意图；

图6是三角面片上光场的衍射示意图；

图7是在三角面片上插值的示意图；

图8是3D模型通过本发明生成计算全息图的示意图；

图9是本发明3D模型的计算全息图产生的流程图。

具体实施方式

下面结合实施例及附图，对本发明作进一步地详细说明，但本发明的实施方式不限于此。

实施例：

如图8所示，3D模型由4个面构成，分别是△ABC、△ABD、△ACD、△BCD。点A的坐标为(-0.85，0，-0.68)，点B的坐标为(0，0，-0.34)，点C的坐标为(-0.22，0.36，0.46)，点D的坐标为(-0.63，0.95，0)，坐标系的刻度单位为毫米。相应各顶点的光场振幅：A点为32，B点为85，C点为172，D点为40，这里光场振幅的单位为伏特/米，以下在没有指明的情况下都取该单位。这里选取的参考光为点光源，发出球面光波，R₀＝0.001，波长为625nm，位置坐标为(0，1，300)。全息面到世界坐标系中的原点距离为500毫米，全息面上一个采样点的S(0.6，0.3，500)的光场U_s(0.6，0.3，500)初值为0。以该点的光场计算为例，其它各采样点的计算过程相同。

使用“背面消隐”技术，确定3D模型中哪些三角面片的光场可以投射到全息面上。在3D模型中，根据各三角面片的顶点的排列顺序，可以确定三角面片的正面，通过三角面片的任意两边对应的矢量叉乘后得到其法线。该法线与全息面法线的点乘除以它们的模可得两法线夹角的余弦，如果该余弦非正，则背对全息面，标记为1。其它标记为0。本3D模型全息面法线向量为(0，0，1)，求得△ABC、△ACD、△ABD、△BCD的法线向量分别为(0.122，0.759，-0.306)、(0.838，0.178，-0.519)、(-0.323，-0.503，0.808)、(-0.638，-0.429，0.018)，它们与全息面法线夹角的余弦分别为-0.372、-0.518、0.804、0.023。因此，只有△ABD、△BCD需要进行计算，标记为0，而△ABC、△ACD标记为1。

先选择三角面片△BCD进行处理，根据步骤(3)对三角面片△BCD进行平移。将其标记为1。根据式(1)计算三条边的长度，分别为BD＝1.190、CB＝0.904、CD＝0.853，单位为毫米，其中长边为BD，其上垂足根据式(6)确定，为(-0.335，0.505，-0.159)，这是进行移动的基点。根据式(7)平移后，点B、C、D的坐标分别为：(0.335，-0.505，-0.181)、(0.115，-0.145，0.619)、(-0.295，0.445，0.159)。

然后根据步骤(4)旋转到标准位置，这需要三次旋转，其旋转矩阵为：

$M_1=\left[\begin{array}{ccc}-0.553& 0.833& 0\\ -0.833& -0.553& 0\\ 0& 0& 1\end{array}\right],$   $M_2=\left[\begin{array}{ccc}0.958& 0& 0.286\\ 0& 1& 0\\ -0.286& 0& 0.958\end{array}\right],$   $M_3=\left[\begin{array}{ccc}1& 0& 0\\ 0& -0.024& 0.999\\ 0& -0.999& -0.024\end{array}\right]$

根据式(14)旋转△BCD的三个顶点后的坐标依次分别为：(-0.633，0，0)、(0，0.646，0)、(0.557，0，0)。至此，△BCD位于标准位置。

根据步骤(5)，全息面各采样点S(x_s，y_s，z_s)作同样的平移和旋转，变换后的全息面上各采样点坐标为S′(u，v，η)，使变换后的各采样点与标准位置的△BCD的相对位置同原采样点与初始3D模型中△BCD的相对位置保持不变。全息面采样点S(0.6，0.3，500)作同样的平移和旋转，将上面求出的△BCD的平移基点坐标与旋转矩阵代入式(15)，可得变换后的采样点S′的坐标(142.299，479.366，-10.919)。

根据步骤(6)，将处于标准位置的△BCD及采样点S′(u，v，η)的相关参数代入公式(32-1)及相关联公式，可以计算出三角面片△BCD作用在全息面上任意采样点S′(u，v，η)的光场U(u，v，η)。这里，点D、C、B分别对应公式中的A、B、C。先计算出下列与△BCD参数相关的数值，见表一：

表一：

 

f11＝1.248×10-10	f12＝-5.396×10-10	f13＝6.675×10-11	f14＝-1.621×10-10
				f15＝7.514×10-11	f16＝-1.008×10-7	f17＝-2.598×10-14j	f18＝-3.231×10-14j
f19＝4.748×10-14j	f21＝4.587×10-15	f22＝-1.277×10-15	f23＝-1.077×10-14
				f24＝1.512×1015	f25＝6.426×10-15	f31＝-1.276×10-10	f32＝6.387×10-10
f33＝-8.849×10-11	f34＝4.097×10-14j	f35＝-1.302×10-13j	f36＝8.926×10-14j
				f37＝1.713×10-10	f38＝-8.063×10-11	f41＝-3.395×10-14j	f42＝3.395×10-14j
f43＝7.885×10-15j	f44＝-7.885×10-15j	r1＝0.568-0.823j	r2＝0.568+0.823j
				r3＝0.886+0.465j

将上述值及采样点S′的坐标(142.299，479.366，-10.919)代入公式(32)及其相关联公式，可求出U(142.299，479.366，-10.919)＝-1.001×10^-5+7.581×10^-6j。

最后把全息面上各采样点S′(u，v，η)的光场U(u，v，η)累加到全息面旋转前各采样点的光场值U_s(x_s，y_s，z_s)上，此刻完成了对△BCD的处理。将U(142.299，479.366，-10.919)累加到全息面上采样点S(0.6，0.3，500)的原光场值U_s(0.6，0.3，500)上，即U_s(0.6，0.3，500)＝-1.001×10^-5+7.581×10^-6j。

根据步骤(7)判断，还需再处理三角面片△ABD，将其标记为1。根据步骤(3)对三角面片△ABD进行平移。根据式(1)计算三条边的长度，分别为：BD＝1.120、AB＝0.916、AD＝1.189，单位为毫米，其中长边为BD，其上垂足根据式(6)确定，为(-0.187，0.282，-0.239)，这是△ABD及采样点进行移动的基点。根据式(7)平移后，点A、B、D的坐标分别为：(-0.663，-0.282，-0.441)、(0.187，-0.282，-0.101)、(-0.443，0.668，0.239)。

然后根据步骤(4)旋转到标准位置，这需要三次旋转，其旋转矩阵为：

$M_1=\left[\begin{array}{ccc}-0.553& 0.833& 0\\ -0.833& -0.553& 0\\ 0& 0& 1\end{array}\right],$ $M_2=\left[\begin{array}{ccc}0.958& 0& 0.286\\ 0& 1& 0\\ -0.286& 0& 0.958\end{array}\right],$ $M_3=\left[\begin{array}{ccc}1& 0& 0\\ 0& 0.839& -0.545\\ 0& 0.545& 0.839\end{array}\right]$

根据式(14)旋转△ABD的三个顶点后的坐标依次为：(0，0.845，0)、(-0.353，0，0)、(0.837，0，0)。至此，△ABD位于标准位置。

根据步骤(5)，全息面采样点的S(0.6，0.3，500)与△ABD作同样的平移和旋转，以保持相对位置不变。将上面求出的△ABD的平移基点与旋转矩阵代入式(15)，可得采样点S′的坐标为(142.579，-261.730，401.757)。

根据步骤(6)，将处于标准位置的△ABD及采样点S′(u，v，η)的相关参数代入公式(32-1)及相关联公式，先计算出下列与△ABD本身相关的数值，见表二：

表二：

 

f11＝1.572×10-10	f12＝-8.413×10-11	f13＝3.119×10-11	f14＝5.889×10-11
				f15＝-1.165×10-10	f16＝-2.976×10-7	f17＝2.229×10-14j	f18＝-3.823×10-14j
f19＝2.280×10-15j	f21＝3.221×10-15	f22＝8.917×10-15	f23＝-9.397×10-15
				f24＝-2.155×10-14	f25＝1.881×10-14	f31＝-5.995×10-10	f32＝2.033×10-10
f33＝2.798×10-11	f34＝-5.528×10-14j	f35＝4.745×10-14j	f36＝7.834×10-15j
				f37＝3.835×10-10	f38＝-1.805×10-10	f41＝1.682×10-15j	f42＝-1.682×10-15j
f43＝1.765×10-14j	f44＝-1.765×10-14j	r1＝0.886+0.465j	r2＝1.000
				r3＝-0.355+0.935j

计算出全息面上各采样点S′(u，v，η)的光场U(u，v，η)，然后累加到全息面上各采样点旋转前的光场值U_s(x_s，y_s，z_s)上，此刻完成了对△ABD的处理。将上述值及采样点S′(142.579，-261.730，401.757)的坐标代入公式(32)及相关联公式，可求出U(142.579，-261.730，401.757)＝-0.0024-0.0002j。最后将U(142.579，-261.730，401.757)累加到全息面上采样点S(0.6，0.3，500)的原光场值U_s(0.6，0.3，500)上，即U_s(0.6，0.3，500)＝-0.00241-0.000193j。

根据步骤(7)判断，所有三角面片都标记为1，没有需要处理的三角面片了。接下来叠加参考光。因为这里选取的参考光为点光源，发出球面光波，R₀＝0.001，所以将点光源R的坐标(0，1，300)与采样点S的坐标(0.6，0.3，500)代入公式(33)，计算出在点S(0.6，0.3，500)处的参考光场为-0.00404+0.00294j。将其叠加到U_s(0.6，0.3，500)上，最后获得计算全息图上的一个采样点S(0.6，0.3，500)的光场-0.00645+0.00281j。

在全息面上设置若干采样点，都按上述步骤计算其在计算全息图上的光场，就得到了计算全息图。

分别利用本发明与现有方法生成计算全息图的时间对比，见表三，新方法与传统方法均用C语言实现，并运行于Athlon xp 2500+CPU上，内存为256M。

从表三可以看出，新方法的速度比传统方法快2个数量级，并且随着现有技术中3D模型采样间距的减小，其优势更加明显。其原因是新方法计算过程与3D模型采样间距并无直接联系，它通过“连续”函数的积分，直接得到了全息图上采样点的数值。

表三：

 

3D模型采样间距(um)	10	8	5	2	1
						全息图大小(采样点数)	1200	2000	4000	8000	10000
传统方法耗时(秒)	1.42	11.62	27.65	215.32	1610.34
						新方法耗时(秒)	0.014	0.102	0.21	1.29	7.18
加速比	101	114	132	167	224

上述实施例为本发明较佳的实施方式，但本发明的实施方式并不受所述实施例的限制，其他的任何未背离本发明的精神实质与原理下所作的改变、修饰、替代、组合、简化，均应为等效的置换方式，都包含在本发明的保护范围之内。

Claims (3)

1、一种3D模型的计算全息图快速生成方法，包括以下步骤：

(1)根据需要确定全息面的位置与大小，选择全息面上的各采样点的位置，得到各采样点的坐标S(x_s，y_s，z_s)，置此时各采样点S(x_s，y_s，z_s)的光场初值U_s(x_s，y_s，z_s)为0；

(2)标记3D模型中通过“背面消隐”后保留的三角面片为0，表示其没被处理过，其具体方法是：在3D模型中，根据各三角面片的顶点的排列顺序和3D模型的规定，确定三角面片的正面，通过三角面片的任意两边对应的矢量叉乘后得到其法线；该法线与全息面法线的点乘，除以它们的模可得两法线夹角的余弦，如果该余弦非正，则三角面片背对全息面，标记为1，其它标记为0；

(3)选择任意标记为0的三角面片，标记为1，在世界坐标系中，设标记为1的三角面片的三个顶点坐标分别是A_w(x_wa，y_wa，z_wa)、B_w(x_wb，y_wb，z_wb)和C_w(x_wc，y_wc，z_wc)，其中最长边为A_wC_w，点B_W到边A_wC_w的垂足为V(x_v，y_v，z_v)，则根据下列公式对该标记为1的三角面片进行平移，使垂足V移动到世界坐标系的原点0：

$t=\frac{(x_{\mathrm{wb}}-x_{\mathrm{wc}})(x_{\mathrm{wa}}-x_{\mathrm{wc}})+(y_{\mathrm{wb}}-y_{\mathrm{wc}})(y_{\mathrm{wa}}-y_{\mathrm{wc}})+(z_{\mathrm{wb}}-z_{\mathrm{wc}})(z_{\mathrm{wa}}-z_{\mathrm{wc}})}{{(x_{\mathrm{wa}}-x_{\mathrm{wc}})}^2+{(y_{\mathrm{wa}}-y_{\mathrm{wc}})}^2+{(z_{\mathrm{wa}}-z_{\mathrm{wc}})}^2}$

$\left[\begin{array}{c}{x}_v\\ y_v\\ z_v\end{array}\right]=\left[\begin{array}{c}{x}_{\mathrm{wc}}+t\&CenterDot;(x_{\mathrm{wa}}-x_{\mathrm{wc}})\\ y_{\mathrm{wc}}+t\&CenterDot;(y_{\mathrm{wa}}-y_{\mathrm{wc}})\\ z_{\mathrm{wc}}+t\&CenterDot;(z_{\mathrm{wa}}-z_{\mathrm{wc}})\end{array}\right]$

$\left\{\begin{array}{c}{x}_{1a}=x_{\mathrm{wa}}-x_v,y_{1a}=y_{\mathrm{wa}}-y_v,z_{1a}=z_{\mathrm{wa}}-z_v\\ x_{1b}=x_{\mathrm{wb}}-x_v,y_{1b}=y_{\mathrm{wb}}-y_v,z_{1b}=z_{\mathrm{wb}}-z_v\\ x_{1c}=x_{\mathrm{wc}}-x_v,y_{1c}=y_{\mathrm{wc}}-y_v,z_{1c}=z_{\mathrm{wc}}-z_v\end{array}\right.$

则该标记为1的三角面片平移后三个顶点的坐标分别是A₁(x_1a，y_1a，z_1a)、B₁(x_1b，y_1b，z_1b)和C₁(x_1c，y_1c，z_1c)；

(4)根据下列公式对步骤(3)中平移后标记为1的三角面片进行旋转：

$M_1=\left[\begin{array}{ccc}\cos \mathrm{\&alpha;}& \sin \mathrm{\&alpha;}& 0\\ -\sin \mathrm{\&alpha;}& \cos \mathrm{\&alpha;}& 0\\ 0& 0& 1\end{array}\right]$

$M_2=\left[\begin{array}{ccc}\cos \mathrm{\&beta;}& 0& \sin \mathrm{\&beta;}\\ 0& 1& 0\\ -\sin \mathrm{\&beta;}& 0& \cos \mathrm{\&beta;}\end{array}\right]$

$M_3=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos \mathrm{\&gamma;}& \sin \mathrm{\&gamma;}\\ 0& -\sin \mathrm{\&gamma;}& \cos \mathrm{\&gamma;}\end{array}\right]$

M₁为绕z轴旋转的矩阵，M₂为绕y轴旋转的矩阵，M₃为绕x轴旋转的矩阵，其中 $\cos \mathrm{\&alpha;}=\frac{x_{1a}}{\sqrt{x_{1a}^2+y_{1a}^2}},$ $\sin \mathrm{\&alpha;}=\frac{y_{1a}}{\sqrt{x_{1a}^2+y_{1a}^2}},$ $\cos \mathrm{\&beta;}=\frac{x_{2a}}{\sqrt{x_{2a}^2+z_{2a}^2}},$ $\sin \mathrm{\&beta;}=\frac{z_{2a}}{\sqrt{x_{2a}^2+z_{2a}^2}},$ $\cos \mathrm{\&gamma;}=\frac{y_{3b}}{\sqrt{y_{3b}^2+z_{3b}^2}},$ $\sin \mathrm{\&gamma;}=\frac{x_{3b}}{\sqrt{y_{3b}^2+z_{3b}^2}},$ α是该三角面片绕z轴旋转至其最长边到x-z平面内的角度，β是该三角面片绕y轴旋转至其最长边与x轴重合的角度，γ是该三角面片绕x轴旋转至x-y平面内的角度；

$\left[\begin{array}{c}{x}_a\\ y_a\\ z_a\end{array}\right]{=M}_3\&CenterDot;M_2\&CenterDot;M_1\&CenterDot;\left[\begin{array}{c}{x}_{1a}\\ y_{1a}\\ z_{1a}\end{array}\right],\left[\begin{array}{c}{x}_b\\ y_b\\ z_b\end{array}\right]=M_3\&CenterDot;M_2\&CenterDot;M_1\&CenterDot;\left[\begin{array}{c}{x}_{1b}\\ y_{1b}\\ z_{1b}\end{array}\right],\left[\begin{array}{c}{x}_c\\ y_c\\ z_c\end{array}\right]{=M}_3\&CenterDot;M_2\&CenterDot;M_1\&CenterDot;\left[\begin{array}{c}{x}_{1c}\\ y_{1c}\\ z_{1c}\end{array}\right]$

则该标记为1的三角面片旋转后三个顶点的坐标分别是A(x_a，y_a，z_a)、B(x_b，y_b，z_b)和C(x_c，y_c，z_c)；

(5)根据步骤(3)的V(x_v，y_v，z_v)和步骤(4)中的M₁、M₂、M₃，对全息面中各采样点S(x_s，y_s，z_s)进行平移和旋转，以保持全息面与该标记为1的三角形面片的相对位置在三角面片坐标变换后保持不变；

设变换后的采样点及其坐标为S′(u，v，η)，则通过下式求变换后采样点的位置：

$\left[\begin{array}{c}u\\ v\\ \mathrm{\&eta;}\end{array}\right]{=M}_3\&CenterDot;M_2\&CenterDot;M_1\&CenterDot;\left[\begin{array}{c}{x}_s{-x}_v\\ y_s-y_v\\ z_s-z_v\end{array}\right]$

(6)根据下列公式计算该标记为1的三角面片作用到全息面上各采样点S′(u，v，η)的光场U(u，v，η)，然后累加到全息面上各采样点旋转前的光场值U_s(x_s，y_s，z_s)上：

(a)当u≠0且v≠0时：

$U(u,v,\mathrm{\&eta;})=\frac{1}{8}\&CenterDot;j\&CenterDot;\mathrm{\&lambda;}\&CenterDot;y_b\&CenterDot;\mathrm{\&eta;}\&CenterDot;e^{\frac{\mathrm{jk}}{2\mathrm{\&eta;}}(u^2+v^2+{2\mathrm{\&eta;}}^2)}\&CenterDot;\frac{f_1}{f_2}$

将f₁表示成u、v与η的三元多项式的形式，多项式的项分别为u⁴、u³·v、u³·η、u²·v²、u²·v·η、u·v³、u·v²·η与v³·η：

$f=(f_{11}\&CenterDot;r_1^{u/z}+f_{12}\&CenterDot;r_2^{v/z}+f_{13}\&CenterDot;r_3^{u/z})\&CenterDot;u^4+(f_{14}\&CenterDot;r_1^{u/z}+f_{15}\&CenterDot;r_2^{v/z}+f_{16}\&CenterDot;r_3^{u/z})\&CenterDot;u^3\&CenterDot;v+$

$(f_{17}\&CenterDot;r_1^{u/z}+f_{18}\&CenterDot;r_2^{v/z}+f_{19}\&CenterDot;r_3^{u/z})\&CenterDot;u^3\&CenterDot;\mathrm{\&eta;}+(f_{31}\&CenterDot;r_1^{u/z}+f_{32}\&CenterDot;r_2^{v/z}+f_{33}\&CenterDot;r_3^{u/z})\&CenterDot;u^2\&CenterDot;v^2$

$+(f_{34}\&CenterDot;r_1^{u/z}+f_{35}\&CenterDot;r_2^{v/z}+f_{36}\&CenterDot;r_3^{u/z})\&CenterDot;u^2\&CenterDot;v\&CenterDot;\mathrm{\&eta;}+(f_{37}\&CenterDot;r_1^{u/z}+f_{38}\&CenterDot;r_3^{u/z})\&CenterDot;u\&CenterDot;v^3+$

$+(f_{41}\&CenterDot;r_1^{u/z}+f_{42}\&CenterDot;r_3^{u/z})\&CenterDot;u\&CenterDot;v^2\&CenterDot;\mathrm{\&eta;}+(f_{43}\&CenterDot;r_1^{u/z}+f_{44}\&CenterDot;r_3^{u/z})\&CenterDot;v^3\&CenterDot;\mathrm{\&eta;}$

f₁的多项式系数计算时，需要先计算其系数r₁、r₂、r₃，其计算公式为：

$r_1=e^{-2\mathrm{j\&pi;}{x}_c/\mathrm{\&lambda;}},$ $r_2=e^{-2\mathrm{j\&pi;}{y}_b/\mathrm{\&lambda;}},$ $r_3=e^{-2\mathrm{j\&pi;}{x}_a/\mathrm{\&lambda;}}$

计算多项式项u⁴的系数时，需要先计算其系数f₁₁、f₁₂、f₁₃，其计算公式为：

$f_{11}=2\mathrm{\&pi;}\&CenterDot;(x_a^2\&CenterDot;x_c^2-x_a^3\&CenterDot;x_c)\&CenterDot;U_c,$                $f_{12}=2\mathrm{\&pi;}\&CenterDot;(x_a^3\&CenterDot;x_c-2x_a^2\&CenterDot;x_c^2{+x}_c^3\&CenterDot;x_a)\&CenterDot;U_b$

$f_{13}=2\mathrm{\&pi;}\&CenterDot;(x_a^2\&CenterDot;x_c^2-x_c^3\&CenterDot;x_a)\&CenterDot;U_a$

计算多项式项u³·v的系数时，需要先计算其系数f₁₄、f₁₅、f₁₆，其计算公式为：

$f_{14}=2\mathrm{\&pi;}\&CenterDot;y_b\&CenterDot;(x_a^3-2x_a\&CenterDot;x_c^2{+x}_a^2\&CenterDot;x_c)\&CenterDot;U_c,$

$f_{15}=2\mathrm{\&pi;}\&CenterDot;y_b\&CenterDot;(x_c^2\&CenterDot;x_a-x_c^3-x_a^3+x_a^2\&CenterDot;x_c)\&CenterDot;U_b,$

$f_{16}=2\mathrm{\&pi;}\&CenterDot;y_b\&CenterDot;(x_c^2\&CenterDot;x_a-2x_a^2+x_a^3)\&CenterDot;U_a$

计算多项式项u³·η的系数时，需要先计算其系数f₁₇、f₁₈、f₁₉，其计算公式为：

$f_{17}=j\&CenterDot;\mathrm{\&lambda;}\&CenterDot;[(x_a^2\&CenterDot;x_c-x_a^3)\&CenterDot;(U_b-U_c)+x_a^2\&CenterDot;x_c\&CenterDot;U_a]$

$f_{18}=j\&CenterDot;\mathrm{\&lambda;}\&CenterDot;[({2x}_c^2\&CenterDot;x_a-x_a^2\&CenterDot;x_c-x_c^3)\&CenterDot;U_a+(x_a^3-x_a^2\&CenterDot;x_c-x_c^2\&CenterDot;x_a+x_c^3){\&CenterDot;U}_b]$

$+({2x}_a^2\&CenterDot;x_c-x_c^2\&CenterDot;x_a-x_a^3)\&CenterDot;U_c]$

$f_{19}=j\&CenterDot;\mathrm{\&lambda;}\&CenterDot;\left[(x_c^3-2x_c^2\&CenterDot;x_a)\&CenterDot;U_a+(x_c^2\&CenterDot;x_a-x_c^3)\&CenterDot;U_b+x_a\&CenterDot;x_c^2\&CenterDot;U_c\right]$

计算多项式项u²·v²的系数时，需要先计算其系数f₃₁、f₃₂、f₃₃，其计算公式为：

$f_{31}=2\mathrm{\&pi;}\&CenterDot;y_b^2\&CenterDot;(x_c^2-2x_a^2{+x}_a\&CenterDot;x_c)\&CenterDot;U_c,$     $f_{32}=2\mathrm{\&pi;}\&CenterDot;y_b^2\&CenterDot;{(x_a-x_c)}^2\&CenterDot;U_b$

$f_{33}=2\mathrm{\&pi;}\&CenterDot;y_b^2\&CenterDot;(x_a^2-2x_c^2{+x}_a\&CenterDot;x_c)\&CenterDot;U_a$

计算多项式项u²·v·η的系数时，需要先计算其系数f₃₄、f₃₅、f₃₆，其计算公式为：

$f_{34}=j\&CenterDot;\mathrm{\&lambda;}\&CenterDot;y_b\&CenterDot;[2(x_a^2-x_a\&CenterDot;x_c)\&CenterDot;U_b-(x_a^2+2x_a{x}_c)\&CenterDot;U_a+({4x}_a\&CenterDot;x_c-x_a^2)\&CenterDot;U_c]$

$f_{35}=j\&CenterDot;\mathrm{\&lambda;}\&CenterDot;y_b\&CenterDot;{(x_a-x_c)}^2\&CenterDot;(U_a-2U_b+U_c)$

$f_{36}=j\&CenterDot;\mathrm{\&lambda;}\&CenterDot;y_b\&CenterDot;[({4x}_a\&CenterDot;x_c-x_c^2)\&CenterDot;U_a+2(x_c^2-x_a\&CenterDot;x_c)\&CenterDot;U_b-(x_c^2+2x_a\&CenterDot;x_c)\&CenterDot;U_c]$

计算多项式项u·v³的系数时，需要先计算其系数f₃₇、f₃₈，其计算公式为：

$f_{37}=2\mathrm{\&pi;}\&CenterDot;y_b^3\&CenterDot;(x_a-x_c)\&CenterDot;U_c,$

$f_{38}=2\mathrm{\&pi;}\&CenterDot;y_b^3\&CenterDot;(x_c-x_a)\&CenterDot;U_a$

计算多项式项u·v²·η的系数时，需要先计算其系数f₄₁、f₄₂，其计算公式为：

$f_{41}=j\&CenterDot;\mathrm{\&lambda;}\&CenterDot;y_b^2\&CenterDot;[({2x}_a+x_c)\&CenterDot;U_a+(x_c-x_a)\&CenterDot;U_b-(x_a+2x_c)\&CenterDot;U_c],f_{42}=-f_{41}$

计算多项式项v³·η的系数时，需要先计算其系数f₄₃、f₄₄，其计算公式为：

$f_{43}=j\&CenterDot;\mathrm{\&lambda;}\&CenterDot;y_b^3\&CenterDot;(U_c-U_a),$        f₄₄＝-f₄₃

将f₂表示成u与v的二元多项式的形式：

f₂＝f₂₁·u⁶+f₂₂·u⁵·v+f₂₃·u⁴·v²+f₂₄·u³·v³+f₂₅·u²·v⁴

多项式项u⁶的系数f₂₁用下式计算：

$f_{21}={\mathrm{\&pi;}}^3\&CenterDot;x_c^2\&CenterDot;(x_a^3-x_a^2\&CenterDot;x_c)$

多项式项u⁵·v的系数f₂₂用下式计算：

$f_{22}={2\mathrm{\&pi;}}^3\&CenterDot;y_b\&CenterDot;(x_c^3\&CenterDot;x_a-x_c\&CenterDot;x_a^3)$

多项式项u⁴·v²的系数f₂₃用下式计算：

$f_{23}={\mathrm{\&pi;}}^3\&CenterDot;y_b^2\&CenterDot;(x_a-x_c)\&CenterDot;(x_c^2+4x_a\&CenterDot;x_c+x_a^2)$

多项式项u³·v³的系数f₂₄用下式计算：

$f_{24}={2\mathrm{\&pi;}}^3\&CenterDot;y_b^3\&CenterDot;(x_c^2-x_a^2)$

多项式项u²·v⁴的系数f₂₅用下式计算：

$f_{25}={\mathrm{\&pi;}}^3\&CenterDot;y_b^4\&CenterDot;(x_a-x_c)$

(b)当u＝0且v≠0时：

$U(u,v,\mathrm{\&eta;})=\frac{1}{8}\&CenterDot;e^{\frac{\mathrm{jk}}{2\mathrm{\&eta;}}(u^2+v^2+2{\mathrm{\&eta;}}^2)}\&CenterDot;[f_3\&CenterDot;(r_4-1)\&CenterDot;z^2\&CenterDot;v^{-3}+(f_4-f_5\&CenterDot;r_4)\&CenterDot;z\&CenterDot;v^{-2}+f_6\&CenterDot;v^{-1}]$

其中，r₄、f₃、f₄、f₅、f₆及下面要使用的r₅这些系数是为了方便计算引入的符号，它们的计算公式如下：

$r_4=e^{-2\mathrm{j\&pi;}{\mathrm{vy}}_b/\mathrm{\&lambda;}/\mathrm{\&eta;}},$         $r_5=2j\&CenterDot;\mathrm{\&lambda;}\&CenterDot;y_b^{-1}\&CenterDot;{\mathrm{\&pi;}}^{-2}\&CenterDot;(x_c-x_a)$

$f_3=(U_a+U_c-2U_b)\&CenterDot;(x_c-x_a)\&CenterDot;{\mathrm{\&lambda;}}^2\&CenterDot;y_b^{-2}\&CenterDot;{\mathrm{\&pi;}}^{-3}$

f₄＝r₅·(U_a+U_c-U_b)，f₅＝r₅·U_b，f₆＝2·(U_a+U_c)·(x_c-x_a)·π^-1

(c)当u＝0且v＝0时：

$U(u,v,\mathrm{\&eta;})=\frac{y_b\&CenterDot;(x_a-x_c)\&CenterDot;(U_a+U_b+U_c)}{6j\&CenterDot;\mathrm{\&lambda;}\&CenterDot;\mathrm{\&eta;}}\&CenterDot;e^{\frac{\mathrm{jk}}{2\mathrm{\&eta;}}(u^2+v^2+{2\mathrm{\&eta;}}^2)}$

在上述(a)、(b)、(c)各情形中，j²＝-1，λ为波长，k为波数，定义为k＝2π/λ，e为自然常数，U_a、U_b、U_c是x-y平面内三角面片三个顶点A、B、C的光场；

(7)判断是否有标记为0的三角面片，若有，转步骤(3)；若无，则求各采样点S(x_s，y_s，z_s)处参考光的光场，叠加到U_s(x_s，y_s，z_s)上，产生计算全息图。

2、根据权利要求1所述的一种3D模型的计算全息图快速生成方法，其特征在于：步骤(7)中所述参考光为球面光，设参考点光源所在位置为(x_r，y_r，z_r)，发出的球面参考光在(x_s，y_s，z_s)处的光场为：

$R(x_s,y_s,z_s)=R_0(1/\sqrt{{(z_s-z_r)}^2+{(x_s-x_r)}^2+{(y_s-y_r)}^2})e^{\mathrm{jk}\sqrt{{(z_s-z_r)}^2+{(x_s-x_r)}^2+{(y_s-y_r)}^2}}$

其中，R(x_s，y_s，z_s)表示坐标为(x_s，y_s，z_s)的采样点的参考光光场，R₀表示与所选光源相关的复常数，其模为点光源发出的初始球面光波振幅，相位为点光源发出的初始球面光波的相位。

3、根据权利要求1所述的一种3D模型的计算全息图快速生成方法，其特征在于：步骤(7)中所述参考光为平面光，在点(x_s，y_s，z_s)处的光场R(x_s，y_s，z_s)为：

$R(x_s,y_s,z_s)=R_p{e}^{\mathrm{jk}(x_s{x}_p+y_s{y}_p+z_s{z}_p)}$

其中，设平面光传播的单位方向向量为

朝向全息面，R_p表示与所选平面光源相关的复常数，其模为光源处平面光波振幅，相位为光源处平面光波的初始相位。

Images