平面和球面非线性折射和反射的实时光线跟踪方法属于计算机图形技术领域。
1996年，M.Levoy和P.Hanrahan.等人在SIGGRAPH上发表论文，提出了用4维光场来实现图形、图像的绘制的方法。在预处理中用多维光场计算图像的映像，然后当实时显示时，在每一帧中重新使用映像结果。但光场方法中仍有图像信息存储、走样和快速查找等问题有待进一步解决。
视点无关的环境映射的方法是W.Heidrich等人在1998年提出来的。利用环境映射，通过硬件的支持来实现真实感图形的显示，能达到交互的速率，但是却不能达到实时。它一般多用于处理反射现象，并且常常存在明显的误差。
1998年Ofek和Rappoport提出了利用反射虚物体和扩张映射(explosion map)来实现曲面反射的光线跟踪效果，它针对反射体将所有可能处在反射区域内的顶点进行变换，生成虚物体，然后同真实物体一样进行统一处理；但是他们不能处理折射问题，且扩张映射的误差较大，在递归过程中误差的累积和传播十分严重，生成的结果图像质量不高。虽然他们采用插值方法加以改进，但效果有限，却增加了不少计算量。
目前有许多交互式图形系统利用硬件，每秒钟能处理大量的三维面片，但基本上都只能提供深度缓存(Z-Buffer)算法的绘制效果，不能处理反射和折射场景的实时真实感图形问题。
本发明的特征在于：对于被折射或反射的多面体，首先计算出该多面体各个顶点的虚顶点，然后根据该多面体原有的点和面之间的拓扑关系依次把所有虚顶点连接成由许多虚面构成且最终形成反射或折射的光学映射虚物体，简称为虚物体；在显示图形时，把该虚物体用递归算法递归地投射到反射或折射表面上，从而得到与光线跟踪类似的图像效果；其中，虚顶点P′可表达为：P′=V+lB-V|B-V|,]]>其中V为视点，l为光线L经过折射或反射到达与该虚顶点对应的一个场景中的点P时所经过的路径的长度总和；B为投影平面上的一个采样点。
对于场景在折射体外的平面折射体而言，它依次含有如下步骤：(1)计算各虚顶点P′：P′=V+lI1-V|I1-V|;]]>l＝|V-I1|+|I1-I2|+|I2-P|，I1=Q+dU-Q|U-Q|,]]>Q＝V-w1N，U=V-C3+CTVCT(P-V)(P-V),]]>其中，l为光线经过折射到达场景中的点P时所经过的路径的长度总和；w1为视点V与折射体的水平距离；d为点I1与视点V的垂直距离；Q为折射体入射面下端与视点V处于同一个水平位置上的一点；U为线VP与折射体入射面的交点，折射体入射面的方程设为C0x+C1y+C2z+C3＝0，并记C＝(C0，C1，C2)T；I1、I2分别为光线L与平面折射体投影平面即入、出两个折射面的交点，N为点I1处的单位法向矢量。
对于场景在折射体内的平面折射体而言，它依次含有如下步骤：(1)计算各虚顶点P′：P′=V+lI-V|I-V|;]]>l＝|V-I|+|I-P|，I=Q+dU-Q|U-Q|,]]>Q＝V-w1N，U=V-C3+CTVCT(P-V)(P-V),]]>其中，I为光线L与平面折射体的入射面的交点，该入射面的方程设为C0x+C1y+C2z+C3＝0，并记C＝(C0，C1，C2)T；N为点I处的单位法向矢量；d为点I1与视点V的垂直距离；Q为折射体入射面下端与视点V处于同一个水平位置上的一点；U为线VP与折射体入射面的交点；对于折射面相互垂直的三维平面折射而言，它依次含有如下步骤：(1)计算各虚顶点P′：P′=V+dI1-V|I1-V|;]]>d＝|V-I1|+|I1-I2|+|I2-P|，I1＝V+dx+dy+dz，dx＝xex，dy＝yey，dz＝zez，ey＝N2，ez＝-N1，ey叉乘ez再单位化得到ex，其中，P是场景中即物体上的点，V是视点，N1是入射面外法向量，N2是出射面外法向量，入射面与出射面垂直，N1、N2是单位向量，ex，ey，ez是局部坐标系下的基，x，y，z是局部坐标系下的坐标，dx，dy，dz是I1点相对于视点的偏移量；I2=I1-N1C1-x2+y2η2(x2+y2+z2)-e^Cx2+y2ηx2+y2+z2,]]>其中，e^=V-I1-N1[(V-I1)·N1]|V-I1-N1[(V-I1)·N1]|,C=η(a-y)yx2+y2+z2,]]>η是三维折射体的折射率，α为视点到折射体出射面的垂直距离。对于球面反射体而言，它依次含有如下步骤：(1)计算各虚顶点P′：P′=V+(x+y)I-V|I-V|,]]>I为光线L照到球面时的反射点，它可表达为：I＝T·LI，其中，LI＝[LIx，LIy，Pz，1]T，它是球面上点I在局部坐标系下的坐标；；T为从局部坐标系到世界坐标系的变换矩阵：T=uxvxwx0uyvywy0uzvzwz00001,]]>
其中wx＝Nx，wy＝Ny，wz＝Nz；N＝(Nx，Ny，Nz)T为视点V、球面的球心O和场景中的点P组成的平面VPO的法向量；以O为原点，以N为z轴，以OP为y轴，以OP×N为x轴，建立局部坐标系；ux＝vywz-wyvz，uy＝wxvz-vxwz，uz＝vxwy-wxvy；vx、vy和vz分别是单位矢量(P-O)/|P-O|的三个坐标分量；x，y分别为球面上的点I到视点V和场景中的点P之间的距离：x＝|V-I|，y＝|I-P|；应用公式I＝T·LI便可把LI局部坐标系变换为世界坐标系。
对于球面折射体而言，它依次含有如下步骤：(1)计算各虚顶点P′：P′=V+(x+|I1-I2|+y)I1-V|I1-V|,]]>其中，I1、I2分别是光线射到球面折射体的入射点和出射点，用公式I＝T·LI，把I1、I2的局部坐标变换成世界坐标系的坐标；x，y分别为视点V到入射点I1，以及出射点I2到场景中的一点P的距离，即x＝|V-I1|，y＝|I2-P|；T为从局部坐标系到世界坐标系的变换矩阵：T=uxvxwx0uyvywy0uzvzwz00001,]]>其中，ux＝vywz-wyvz，uy＝wxvz-vxwz，uz＝vxwy-wxvy；vx、vy和vz分别是单位矢量(P-O)/|P-O|的三个坐标分量；wx＝Nx，wy＝Ny，wz＝Nz；N＝(Nx，Ny，N2)T为视点V、球面的球心O和场景中的点P组成的平面VPO的法向量；以O为坐标原点，以N为z轴，以OP为y轴，以OP×N为x轴，建立局部坐标系；所述的用递归算法来显示图形的过程，含有如下步骤：算法1预处理第1步，对场景中的物体构造二叉空间剖分(BSP)树，如果场景中含有球形反射或折射体，则用紧密包容该球的正方体的6个面代替该球参与构造BSP树；第2步，收集场景中各种反射和折射体，并形成反射和折射体对象表；第3步，根据当前的视点，确定反射、折射体表中的各反射、折射体的先后顺序；由于BSP树与视点无关，因而能很快的得到反射、折射体的顺序；第4步，计算虚顶点并生成所有的虚物体，产生虚场景；算法2为折射体或反射体产生虚场景：i)是否到达最大的递归层次？是则返回；
ii)如果为反射体，计算反射虚场景；如果为折射体，计算折射虚场景；iii)根据视点对虚场景中的反射和折射体进行排序；iv)如果虚场景中还存在反射或折射体，则对虚场景中的反射、折射体递归调用本算法，转i)；否则转v)；v)返回；算法2有以下特点：(1)上面的算法是一个递归算法，递归主要是针对场景中的反射、折射体进行的；(2)对于反射体和折射体是统一处理的；(3)对于反射体，生成虚场景使用的是反射公式，对应的面片的顶点顺序要改变，即反射虚物体面片的顶点顺序与原来的顺序相反，以便确保平面的外法矢指向正确的方向；(4)在任一虚场景中，都要根据视点和BSP树重新确定其中的反射和折射体的顺序，以便后面能够正确地实时显示；(5)如果场景中的反射或折射体的个数大于1算法的终止条件是递归深度不大于所给定的最大递归深度值；如果视点发生变化，对于纯反射虚场景，是不用重算的；但对于任意一个折射体及其以后的递归过程中出现的虚场景必须重新计算；下面是视点发生变化时重新计算虚场景的算法：算法3视点位置改变时重新计算虚场景i)是否是折射体？ii)是，调用算法2；iii)否则，则递归调用本算法；要实现实时绘制场景，可以采用支持OpenGL编程(OpenGL-compatible)的三维图形加速卡进行硬件加速的办法；这时，我们要用到图形硬件的软件接口——OpenGL图形系统；前面的虚物体的生成，为我们实时绘制打下了基础；下面是我们使用的实时绘制算法；算法4实时绘制和显示(本算法是针对真实场景中的某个反射体或折射体的)x)设置OpenGL的绘制参数(禁止颜色缓存写，禁止深度缓存写，禁止深度比较，禁止使用Alpha混合，使所绘制的模版缓冲区的值+1)；xi)绘制反射、折射体本身；xii)如果其虚场景中含有反射或折射体，则转i)，递归地绘制虚场景中的反射、折射体；xiii)设置OpenGL的绘制参数(使所绘制的模版缓冲区的值不变，颜色缓存可写，深度缓存可写，深度比较生效，其他的绘制参数不变)；xiv)在当前模板缓冲区中绘制当前虚场景中的其它虚物体；xv)设置OpenGL的绘制参数(颜色缓存可写，深度缓存可写，深度比较生效，Alpha混合生效，使所绘制的模版缓冲区的值-1)；xvi)绘制反射、折射体本身；
xvii)设置OpenGL的绘制参数(恢复深度缓存可写，颜色缓存可写，禁止Alpha混合，禁止摸板操作)；xviii)返回。
实验证明：它达到了预期目的。
图2任意平面反射的示意图。
图3折射第一种情况--场景在折射体外的示意图。
图4折射第二种情况--场景在折射体内的示意图。
图5折射面相互垂直的三维折射光路图。
图6球面反射的光路示意图。
图7球面折射的光路示意图。
图8确定球面的入射点和出射点示意图。
图9方盒子(a)及其虚物体(b)示意图。
图10为算法1预处理框图。
图11为算法2平面折射体或反射体产生虚场景的程序流程框图。
图12为算法3视点位置改变时重新计算虚场景的程序流程框图。
图13为算法4实时绘制和显示的程序流程框图。
当一条L经过若干次反射或折射，与某一物体的平面M相交于V点时，与传统光线跟踪算法不同，我们不是利用光照模型计算对应象素点的颜色值，而是计算点V的光学映射虚顶点V’。假设光线L经过若干次反射或折射到达V时，经过的路径的长度总和为：d＝d1+d2+Λ+d5则虚顶点V’的计算公式为：V′=A+dB-A|B-A|,]]>其中A是视点，B是投影平面上一个采样点。这就是我们提出的“光学映射虚对象的概念”。
对于被反射或折射的多面体，我们首先计算出每个顶点的虚顶点，然后根据拓扑关系依次将虚顶点连接成虚物体的边，最终形成反射或折射的光学映射虚物体(简称虚物体)。例如一个方盒子，见图9(a)，其顶点P1，P2，...，P8的虚顶点依次为P′1，p′2，...，P′8，则我们可以得到该方盒子对应的虚物体，如图9(b)所示。在显示图形时，将这些虚物体投影到反射或折射表面上。这个过程可以递归地进行，从而得到与光线跟踪类似的图像效果。为了得到高质量的真实感图像，被反射或折射的平面的尺寸不能太大，否则应当将平面细分成更小的矩形或三角形。对于场景中被反射或折射的曲面物体，则应当先将曲面离散成若干小平面，用多面体近似表示曲面物体。
2、虚顶点的计算公式对于平面反射问题，我们可以采用类似光束跟踪算法的思想计算各个顶点的虚顶点：假设有一个过点R＝[Vx，Vy，Vz，1]T的平面，如图2所示，其单位法矢量为N＝[a，b，c，0]T，NTN＝1，其平面方程为ax+by+cz+d＝0其中，d＝-NTR，V为视点。
三维场景中的任意一个点P＝[Px，Py，Pz，1]T相对于该平面作反射变换后到达P′＝[P′x，P′y，P′z，1]T。根据几何关系，可以得到：P′＝TmP其中，Tm=1-2a2-2ab-2ac-2ad-2ab1-2b2-2bc-2bd-2ac-2bc1-2c2-2cd0001,]]>P′是P点的反射虚顶点。
折射问题比平面反射问题要复杂得多。平面折射体的非线性折射有两种情况：物体在折射体外(如图3所示)和物体在折射体内(见图4)两种。下面我们分别予以讨论。
1)二维平面折射(a)物体在折射体外如图3所示，V是视点，P是场景中任意选取的一点，它在折射体的外面，P′是与P对应的虚顶点。求其虚物体的公式推导如下：已知：w1，w2，w3，h，令η＝Sinα/Sinβ＝η2/η1，ξ＝1+w3/w1，根据几何关系，可得Ad4+Bd3+Cd2+Dd+E＝0                                       (1)其中，A＝(1-η2)ξ2，B＝2η2hξ-2hξ，C＝h2+w22-η2h2-w12η2ξ2，D＝2η2hξw12，
E＝-w12η2h2。
从(1)式中可以解出d。根据d的值，然后可以计算出折射光线和折射平面的交点I1和I2，从而求出虚顶点P′。
虚顶点P′求解过程如下：(1)连接V和P，求出VP与平面的交点U：U=V-C3+CTVCT(P-V)(P-V),]]>其中，该平面的方程为C0x+C1y+C2z+C3＝0，另外，C＝(C0，C1，C2)T；(2)Q＝V-w1N；(3)---I1=Q+dU-Q|U-Q|.]]>同理，可以求出I2。
令l＝|V-I1|+|I1-I2|+|I2-P|，则P点的虚顶点P′=V+lI1-V|I1-V|.]]>我们用光线跟踪的方法产生了两组点，然后用上面的公式，反算I1，显然，如果I1相同，则I2也相同，具体的结果如下表所示。
(b)物体在折射体内如鱼缸中的鱼，物体(鱼)是在折射体(鱼缸)内面。这种情况的计算公式应该与第一钟情况是不同的。如图4所示，图中V是视点，P是场景中任意选取的点，它在折射体内部，P′是P对应的虚顶点。类似地，我们可以根据几何关系得到如下一元四次方程：Ad4+Bd3+Cd2+Dd+E＝0                                       (2)其中，A＝η2-1，B＝-2hA，C＝h2A+η2w12-w22，D＝-2η2hw12，
E＝w12η2h2。
从(2)式中可以解出d，然后便可以计算出折射光线和折射平面的交点I，从而求出虚顶点P′。虚顶点P′求解过程如下：(1)连接V和P，VP与镜平面相交于U，U可照搬前面的公式得到；(2)Q＝V-w1N；(3)---I=Q+dU-Q|U-Q|;]]>(4)---p′=V+lI-V|I-V|,]]>其中，l＝|V-I|+|I-P|。
我们用光线跟踪的方法产生了两组点，然后用上面的公式，反算出交点I。具体的结果如下表所示。
(c)解的选取显然，(1)式和(2)式都是多解的方程，它们都至多有4个实根。因此，解方程后，将要对解的结果进行选择。根据光学的原理，光在介质中从一点到达另一点时，将选择最短路径传播。这样，我们连接视点V和场景中的点P，得到一条直线，离该直线最近的解，就是我们所求的点。
2)三维平面折射如图5，已知：视点V，场景中任意选取的点P，折射率η，AD＝a，VA＝z，PF＝b，e=|(VP→)x|,f=|(VP→)y|,g=|(VP→)z|.]]>在图5中，光线从视点V出发，经过介质上表面I1点折射进入介质，然后到达介质右侧面，从I2点折射进入空气，并最终达到场景中的点P。
由Snell折射定律和几何关系，我们就可以得到下面的两个方程：
根据上面的两个方程可以求出x和y，然后由x、y以及z和折射率η的值就可以依次算出I1和I2的三维坐标，从而得到P的虚顶点P’的位置。
(3)和(4)式是一个二元4次非线性方程组，为了提高计算速度，我们采用下面的方法进行近似处理：首先，用下式计算出d：f0+f1d+f2d2+f3d3+f4d4＝0，忽略d的高阶(＞2)项，得f0+f1d+f2d2＝0，其中，f0=(14C7-14C1-C2C1)2(C4-C7)2-2b2C1(C7-C1)(C4-C7)(14C7-14C1+C2C1)+]]>b4C12(C7-C1)2,]]>f1=-2C9b2C1(C4-C7)(C2C5+14C6-14C5+C7-C1)+2b4C12C9(C6-C5)-]]>{2C9b2C1(2C5-C6)+[2C9b2C5+(2C6-2C5)b2C1]C8}(14C7-14C1+C2C1)+]]>2b4C1C5(C7-C1)2+2(C4-C7)(2C5-C6)(14C7-14C1-C2C1)2+]]>2(14C7-14C1-C2C1)(C4-C7)2(C7-C1-C2C5+14C6-14C5),]]>f2=(14C7-14C1-C2C1)2[2(C4-C7)(2-η2)+(2C5-C6)2]+4(14C7-14C1-C2C1)(14C6-]]>14C5-C2C5+C7-C1)(C4-C7)(2C5-C6)+[2(C7-C14-C2C1)(C7-C1-C2+14η2+]]>C6-C5-14)+(C7-C1-C2C5+C6-C54)2](C4-C7)2-2(C7-C1)b2C1(C4-C7)(-14+]]>C2+14η2+C6-C5+C7-C1)-{(2C7-2C1)b2C1(2C5-C6)+(C4-C7)[(2C7-2C1)b2C5+]]>2(C6-C5)b2C1]}(C2C5+C6-C54+C7-C1)-{2(C7-C1)b2C1(2-η2)+[2(C7-C1)b2C5+]]>2(C6-C5)b2C1](2C5-C6)+[2(C7-C1)b2+2(C6-C5)b2C5+(2η2-2)b2C1](C4-C7)}(14C7+]]>C2C1-14C1)+b4C12[2(C7-C1)(η2-1)+(C6-C5)2]+4b4C1C5(C7-C1)(C6-C5)+]]>b4(2C1+C52)(C7-C1)2.]]>在上面的3个式子当中，C1=(a-12)2,C2=(g-z)2,C3=C1+z2,C4=2C1+z2,C5=-2a+1,]]>C6=η2C5,C7=η2C3,C8=C4-C7,C9=C7-C1.]]>
然后将d代入下式，便可以计算出当x＝0时y的值，我们把它记为y0，即y0＝a-(0.5+d)。
通常d是一个很小的量，|d|≤0.5；在上面的计算中，我们忽略了d的高阶(≥3)项，如果保留d的更高阶(3、4阶)项，可以提高计算精度，但是这也会使计算量略有增加。将(4)式在x＝0点展开并作线性近似：y＝y0+y1x，其中，将y＝y0+y1x代入(3)式，得从上式中求出x，并将x的值代入式(4)并求解一个关于y的四次方程，该方程经过整理为：H4y4+H3y3+H2y2+H1y+H0＝0，其中，H4＝e2(2-η2)，H3＝-2aex(2-η2)，H2＝a2x2(2-η2)+e2(1-η2)(x2+z2)-b2x2，H1＝-2aex(1-η2)(x2+z2)，H0＝a2x2(1-η2)(x2+z2)。
再由x、y和z以及η可以计算出点I1和I2的三维坐标，从而求出虚顶点P′。
由式(3)和式(4)组成的方程组来求解x和y的时候，还有另一种方案：在前面得到y0后，将(3)式在x＝0点展开并作线性近似，得到：y＝y0+y1x，其中，y1=(η2-1)y02+η2z2e(η2-1)y0;]]>然后，将y＝y0+y1x代入(4)式，整理得到T4x4+T3x3+T2x2+T1x+T0＝0，其中e2y02[(1-η2)+(2-η2)y12]-b2y02,]]>从上式中求出x，并将x的值代入式(3)，可直接求解出y：y=±x2+(g-z)2x2(e-x)2-η2(x2+z2)η2-1.]]>在后面的计算实验中采用的是前一种方案。
具体计算I1、I2及虚顶点P’的过程如下：：a).求I1设物体上的点是P，视点是V，入射面外法向量n1，出射面外法向量n2，入射面与出射面垂直，n1，n2是单位向量。局部坐标系的基是ex，ey，ez，其中ey＝nz，ex＝-n1，ey叉乘ez，再单位化可得ex。
计算出的x，y，z是局部坐标系下的坐标。设dx，dy，dz是I1点相对于视点的偏移量，有：dx＝xex，dy＝yey，dz＝zez，则I1＝V+dx+dy+dz。
b).求I2已知入射光线方向矢量IN＝I1-V，折射面法矢量n1，折射率η，可计算折射光线方向矢量OUT。由折射光线上一点I1以及其方向，和出射面法向量n2及面上一点，可以求出折射光线与出射平面的交点I2。我们可以给出I2的公式如下：I2=I1-n1C1-x2+y2η2(x2+y2+z2)-e^Cx2+y2ηx2+y2+z2,]]>其中e^=V-I1-n1[(V-I1)·n1]|V-I1-n1[(V-I1)·n1]|,C=η(a-y)yx2+y2+z2.]]>c).求P’从视点到物点，光线走过的距离是d＝‖V-I1‖I1-I2‖+‖I2-P‖，于可以得到虚顶点P’的计算公式：P,=V+d(I1-V)||I1-V||.]]>我们用光线跟踪的方法产生了两组点，然后用上面的公式，反算I1，显然，如果I1相近，则I2也相近。具体的结果如下表所示，其中ε为相对误差，即：用近似公式反算的虚顶点P’到精确的虚顶点p’的距离与视点到虚顶点的距离之比ε＝|P′-p′|/|V-p′|。
3)球面反射和折射(1)球面反射已知：视点V，场景中任意选取的点P，设反射球的半径为r，球心在0点，如图6所示。令a＝|V-P|，b＝|V-O|，c＝|O-P|。
根据几何光学原理可以得到：2bcxyzcos(α+β)＝bcxysin2(α+β)-b2y2z-c2x2z                (5)x(r2+y2-c2)＝y(x2+r2-b2)                                     (6)其中，z＝(a2-(x-y)2)/(4abc)。
取视点V为坐标原点，视点到球心0的距离为单位长度。从(6)式可以解出y，它可以看成是x的函数。将y在x＝b-r处作泰勒展开，并取线性近似y＝A+Bx，                                                    (7)其中，A=(c-b)(c-r)/c,B=b(c-r)c(b-r).]]>将(7)代入(5)式，可以得到如下四次方程：x(A+Bx)[a2-(x-A-Bx)2]cos(α+β)＝2bcx(A+Bx)sin2(α+β)-[b2(A+Bx)2+c2x2][(a2-(x-A-Bx)2]/(2bc)，(8)解出x，然后代入(7)可以求出y，从而计算出虚顶点P’。
下面介绍如何计算虚顶点P’：
I.由于V，0，P三点在一个平面上，我们以0点为原点，平面VP0的法向量N＝[Nx，Ny，Nz]T为z轴，以0P为y轴，0P×N为x轴，建立局部坐标系。这样我们可以得到从世界坐标系到局部坐标系的转换矩阵：T=uxvxwx0uyvywy0uzvzwz00001,]]>其中，wx＝Nx，wy＝Ny，wz＝Nz；ux＝vywz-wyvz，uy＝wxvz-vxwz，uz＝vxwy-wxvy；vx、vy和vz分别是单位矢量(P-O)/|P-O|的三个坐标分量；x，y分别为球面上的点I到视点V和场景中的点P之间的距离：x＝|V-I|，y＝|I-P|；应用公式I＝T.LI便可把LI从局部坐标系变换为世界坐标系下的坐标。
II.在局部坐标系中，分别以0，V，P为圆心，R(球的半径)，x，y为半径做三个圆，分别记为c1，c2，c3。
III.求出c1和c2的两个交点p1，p3，c1和c3的两个交点p3，p4。分别计算p1p3，p1p4，p2p3和p2P4的距离，取距离最小的一对点坐标的平均值，即为球面上点I在局部坐标系上的坐标，记为LI[LIx，LIy，Pz，1]T。应用公式I＝T·LI将LI的局部坐标变换成世界坐标系下的坐标，结果记为I。IV.计算P′=V+(x+y)I-V|I-V|.]]>我们用光线跟踪的方法产生了两组点，然后用上面的公式，反算I1和I2。具体的结果如下表所示。
(2)球面折射如图7所示，球的半径为r，根据几何光学可以得到：{(a2+r2-x2)(b2+r2-y2)x2cosθ4abr2-[x2-(x2+r2-a2)24r2](xycosθab+]]>(9)2η2)+x2}2=x2sin2θ[r2x+y](a2-x2)]2a2b2r2[x2-(x2+r2-a2)24r2],]]>y(x2+r2-a2)＝x(y2+r2-b)。                          (10)不失一般性，取视点V为坐标原点，视点到球心0的距离为单位长度。令x＝a-r+t，                                         (11)对(10)式在t＝0处作泰勒展开，并作线性近似，可得：y=b-r+a(b-r)b(a-r)t,]]>将(11)和(12)式代入(9)式，并忽略t的(5次以上的)高阶项，得y＝G0+G1t+G2t2+G3t3+G4t4，                         (13)其中，G0＝T02-C1S0R0；G1＝2T0T1-C1(S0R1+S1R0)；G2＝2T0T2+T1T1-C1(S0R2+S1R1+S2R0)；G3＝2T0T3+2T1T2-C1(S0R3+S1R2+S2R1+S3R0)；G4＝2T0T4+2T1T3+T22-C1(S0R4+S1R3+S2R2+S3R1+S4R0)；C0＝a-r；C1＝sin2θ/(a2b2r2)；C2＝-1/(4r2)；C3＝cosθ/(4abr2)；C4＝cosθ/(ab)；C5＝2/η2+C4AC0；C6＝C4(A+BC0)；C7＝C4B；S0＝C02+C2(C02-a2+r2)2；S1＝2C0+4C2C0(C02-a2+r2)；S2＝1+2C2(3C02-a2+r2)；S3＝4C2C0；S4＝C2；Q0＝r2C0+A(a2-C02)；Q1＝r2-2AC0+B(a2-C02)；
Q2＝-A-2BC0；Q3＝-B；U0＝Q0C02；U1＝2Q0C0+Q1C02；U2＝Q0+2Q1C0+Q2C02；U3＝Q1+2Q2C0+Q3C02；U4＝Q2+2Q3C0；R0＝Q0U0；R1＝Q0U1+Q1U0；R2＝Q0U2+Q1U1+Q2U0；R3＝Q0U3+Q1U2+Q2U1+Q3U0；R4＝Q0U4+Q1U3+Q2U2+Q3U1；P0＝(a2+r2-C02)(b2+r2-A2)；P1＝-2AB(a2+r2-C02)-2C0(b2+r2-A2)；P2＝-B2(a2+r2-C02)+4C0AB-(b2+r2-A2)；P3＝2C0B2+2AB；P4＝BB；T0＝C3C02P0-C5S0+C02；T1＝C3(C02P1+2C0P0)-C5S1-C6S0+2C0；T2＝C3(C02P2+2C0P1+P0)-C5S2-C6S1-C7S0+1；T3＝C3(C02P3+2C0P2+P1)-C5S3-C6S2-C7S1；T4＝C3(C02P4+2C0P3+P2)-C5S4-C6S3-C7S2；从(13)式可以计算出t，将t代入(11)和(12)式即可分别求出x和y，这样就可以求出虚顶点P’。
确定I1和I2的方法：得到x和y后，分别以V和P为圆心，x和y为半径画圆，在折射平面里和球面分别交于I1′，I1″和I2′，I2″，见图8。球面的入射点I1取自I1′和I1″中的某一个，出射点I2取自I2′和I2″中的某一个。也就是说，I1和I2有四种可能的取法。此时，我们将应用折射定律，选取出我们所要的那一种可能情况。注意到，如果I1和I2是精确的入射点和出射点，则sinβsinα=1η,]]>或者说，1-cos2β1-cos2α=1η2,]]>同时，cosβ=||I1I2||2r,cosα=-x2+r2-a22xr=-y2+r2-b22yr,]]>因而1-||I1I2||24r21-(x2+r2-a22xr)2=1η2,1-||I1I2||24r21-(y2+r2-b22yr)2=1η2.]]>于是，我们的选取准则是，取使得|1-||I1I2||24r21-(x2+r2-a22xr)2-1η2|+|1-||I1I2||24r21-(y2+r2-b22yr)2-1η2|]]>最小的那一组I1I2。
下面介绍如何计算虚顶点P’：I.由于V，O，P，I1，I2五点在一个平面上，我们以O点为原点，平面VPO的法向量N＝[Nx，Ny，Nz]T为z轴，以OP为y轴，N×OP为x轴，建立局部坐标系。这样我们可以得到从局部坐标系到世界坐标系的变换矩阵：T=uxvxwx0uyvywy0uzvzwz00001,]]>其中，wx＝Nx，wy＝Ny，wz＝Nz；uz＝vywz-wyvz，uy＝wxvz-vxwz，uz＝vxwy-wxvy；vx、vy和vz分别是单位矢量 的三个坐标分量。
II.应用公式I＝T·LI将I1和I2的局部坐标转换成世界坐标系的坐标，结果仍然用I1和I2表示。
III.计算P′=V+(x+|I2-I1|+y)I1-V|I1-V|.]]>我们用光线跟踪的方法产生了两组点，然后用上面的公式，反算I1和I2。具体的结果如下表所示。
3、算法实现算法1：预处理，如图10所示。
第1步，对场景中的物体构造二叉空间剖分(BSP)树，如果场景中含有球形反射或折射体，则用紧密包容该球的正方体的6个面代替该球参与构造BSP树；第2步，收集场景中所有的反射和折射体，并形成反射和折射体对象表；第3步，根据当前的视点，确定反射、折射体表中的各反射、折射体的先后顺序；由于BSP树与视点无关，因而能很快的得到反射、折射体的顺序；第4步，计算虚顶点并生成所有的虚物体，产生虚场景。
算法2为平面折射体或反射体产生虚场景，如图11所示。有以下特点：(1)上面的算法是一个递归算法，递归主要是针对场景中的反射、折射体进行的；(2)对于反射体和折射体是统一处理的；(3)对于反射体，生成虚场景使用的是反射公式，对应的面片的顶点顺序要改变，即反射虚物体面片的顶点顺序与原来的顺序相反，以便确保平面的外法矢指向正确的方向；(4)在任一虚场景中，都要根据视点和BSP树重新确定其中的反射和折射体的顺序，以便后面能够正确地实时显示；(5)如果场景中的反射或折射体的个数大于1算法的终止条件是递归深度不大于所给定的最大递归深度值；如果视点发生变化，对于纯反射虚场景，是不用重算的；但对于任意一个折射体及其以后的递归过程中出现的虚场景必须重新计算；下面是视点发生变化时重新计算虚场景的算法：算法3为视点位置改变时重新计算虚场景，如图12所示。
要实现实时绘制场景，就必须采用硬件加速的办法。这时，我们要用到图形硬件的公用的软件接口OpenGL。前面的虚物体的生成，为我们实时绘制打下了基础。下面是我们使用的实时绘制算法。
算法4为实时绘制和显示，如图13所示。