根据一3d对象的2d投影再生和处理该对象的设备和方法

专利名称：根据一3d对象的2d投影再生和处理该对象的设备和方法
技术领域：
本发明是关于处理3D对象的2D投影的设备和方法，特别是关于作2D投影图象的几何分析的设备和方法。
发明的背景描述处理3D对象的2D投影的现存技术方法的状态的发表及其相关技术在后面引用的参考文献中有说明。
本说明中所提到所有发表物和这里引用的发表的揭示内容均在此被采用作为参考。
发明概述本发明是试图提供用于从一3D景象的第一和第二参考视图产生一3D景象的新视图的图象变换设备和方法。
本发明还试图提供为由一3D景象的第一、第二和第三视图产生其3D景象的3D表示的再现方法和设备。
本发明还试图提供为处理3D对象的2D投影的改善的设备和方法。
本发明另外还试图提供用于根据对一3D对象的三个视图定义的三线性张量作此3D对象的图象变换的方法。
本发明还试图提供由一3D景象的第一和第二参考视图产生此3D景象一新视图的图象变换方法，此方法包括提供一3D景象的第一和第二参考视图，利用分别有关第一参考视图、第二参考视图和新视图的几何信息来产生一表明第一、第二和新视图之间的几何关系的三线性张量，及以根据第一和第二参考视图内的第一和第二对应位置以及三线性张量计算多个的各自对应于此不同的第一和第二对应位置的新视图位置来产生此新视图。
此提供步骤可包括例如在第一和第二参考图象中进行扫描。
此利用步骤可包含例如分别在第一参考视图、第二参考视图和新视图中产生第一、第二和第三位置组的步骤。
这样，按照本发明一优选实施例而提供有用于由一3D景象的第一和第二参考视图产生其新视图的图象变换方法，此方法包含提供一3D景象的第一和第二参考视图，利用各自有关第一参考视图、第二参考视图和新视图的几何信息来产生表明第一、第二和新视图之间的几何关系的三线性张量，和以根据第一和第二参考视图中的第一和第二对应位置及三线性张量计算多个各自对应于不同的第一和第二对应位置的新视图位置来产生此新视图。
而按照本发明一优选实施例，此提供步骤包括在第一和第二参考图象中进行扫描。
而且按照本发明一优选实施例，此利用步骤包含分别在第一参考视图、第二参考视图和新视图中生成一组第一、第二和第三对应位置的步骤。
按照本发明一优选实施例还提供从一3D景象的第一、第二和第三视图生成此3D景象的3D表示的3D景象再现方法，此方法包含提供一3D景象的第一、第二和第三视图，利用有关此第一、第二和第三视图的几何信息产生表明第一、第二和第三视图之间的几何关系的三线性张量，和由此三线性张量产生此3D景象的3D表示。
而按照本发明一优选实施例，此产生一3D表示的步骤包含由此三线性张量计算第一和第二视图的表射偏振几何表示，和由此表射偏振几何显示生成此3D表示。
按照本发明的另一优选实施例还提供图象变换设备，用于由一3D景象的第一和第二参考视图产生此3D景象的新视图，此设备包含用于提供一3D景象的第一和第二参考视图的设备，用于利用各自有关第一参考视图、第二参考视图和新视图的几何信息来产生表明第一、第二和新视图之间的几何关系的三线性张量的三线性张量发生器，和通过基于第一和第二参考视图中的第一和第二对应位置和三线性张量计算多个各自对应于此不同的第一和第二对应的新视图位置产生该新视图的新视图发生器。
按照本发明的另一优选实施例还提供用于由一3D景象的第一、第二和第三视图产生此3D景象的3D表示的3D景象再现设备，此设备包含用于提供一3D景象的第一、第二和第三视图的设备，利用有关第一、第二和第三视图的几何信息产生表明第一、第二和第三视图之间的几何关系的三线性张量的三线性张量发生器，和用于由三线性张量产生3D景象的3D表示的3D景象表示发生器。
按照本发明的另一优选实施例还提供由一3D对象的至少一个而最好为三个2D投影来产生关于此3D对象的信息的方法，此方法包含提供一3D对象的至少一个而最好三个2D投影，产生一以aijk＝vi′bjk-vj″aik(i，j，k＝1，2，3)定义的数矩阵，其中aij和bjk分别为矩阵A和B的元，vi′和vj″分别为向量v′和v″的元素，而这些矩阵和向量一起描述此3D对象的三个视图的摄象机参数，并利用此矩阵生成有关此3D对象的信息。
按照本发明一优选实施例还提供一由一3D对象的具有n个对应点Pi和Pi′(i＝1…n)的第一和第二现有视图产生此3D对象的一新视图的方法，此方法包含产生一张量aijk和将点Pi和Pi′(i＝1…n)的值插入到三线性公式，由此来提取表示新视图中的位置的x″、y″值，和根据插入步骤的结果产生新的视图。
按照本发明一优选实施例还提供由一3D对象的至少一个而最好三个2D投影再现此3D对象的方法，此方法包括提供一3D对象的至少一个和最好三个2D投影，产生由αijk＝Vi′bjk-Vj″aik(i，j，k＝1，2，3)表述的数的阵列，其中aij和bik分别为矩阵A和B的元及Vi′和Vi″为向量v′和v″的元素，而这些矩阵和向量一起说明此3D对象的三个视图的摄象机参数，将阵列置换成与3D空间内的三个对应平面相关连的三个单应性矩阵，和利用此三个单应性矩阵再现3D对象。
按照本发明一优选实施例还提供一视觉识别方法，包含提供一3D对象的三个视图间相互间存在三线性关系的三个透视图，和利用视图之间的三线性关系通过对准来进行视觉识别。
按照本发明一优选实施例此方法还包含对此3D对象作再投影。
按照本发明一优选实施例，有关3D对象的信息还包含3D对象的再现。
另外，按照本发明一优选实施例，有关3D对象的信息包含至少一个不再现3D对象而产生3D对象的新视图。
而按照本发明一优选实施例，此至少一个和最好三个的2D投影包含至少一个航摄照片。
按照本发明一优选实施例，此至少一个和最好三个2D投影还包含至少一个卫星照片。
而按照本发明一优选实施例，有关3D对象的信息包括3D对象的至少一个座标。
而按照本发明一优选实施例，3D对象包括一宇宙空间对象。
按照本发明一优选实施例，此3D对象还包含一如轮船的大的对象。
按照本发明一优选实施例，3D对象还包括不实在的对象。
按照本发明一优选实施例，还提供由一3D对象的至少一个和最好三个2D投影产生有关此3D对象的信息的设备，此设备包含用于提供一3D对象的至少一个和最好三个2D投影的设备，一阵列发生器用于生成一以αijk＝Vi′bjk-Vj″aik(i，j，k＝1，2，3)表述的数的阵列，其中aij和bjk分别为矩阵A和B的元及Vi和Vi″分别为向量V′和V″的元素，而这些矩阵和向量一起描述此3D对象的三个视图，和一利用此阵列产生有关3D对象的信息的阵列分析器。
按照本发明的一优选实施例，还提供用于由一3D对象的具有n个对应点Pi和Pi′(i＝1…n)的第一和第二现有视图产生此3D对象的一个新视图的设备，此设备包含用于产生一张量αijk的设备，和用于将点Pi和Pi′(i＝1…n)插入到三线性公式，由此来提取表示新视图中位置的X″，Y″值，和根据此插入步骤的结果产生新视图的设备。
按照本发明一优选实施例还提供由一3D对象的至少一个和最好三个2D投影再现此3D对象的设备，此设备包含提供一3D对象的至少一个和最好三个2D投影的设备，一阵列发生器用于生成一以αijk＝vi′bjk-vj″aik(i，j，k＝1，2，3)，表述的数的阵列，其中aij和bjk分别为矩阵A和B的元及Vi和Vj″分别为向量V′和V″的元素，而这些矩阵和向量一起说明此3D对象三个视图的摄象机参数，一将此阵列置换为与3D空间中三个对应平面相关连的三个单应性矩阵的阵列变换器，和进行利用此三个单应性矩阵来再现此3D对应操作的3D对象再现器。
按照本发明一优选实施例还提供视觉识别设备，包含有用于提供一3D对象的相互间存在三线性关系的三个透视图的设备，和用于利用视图间的三线性关系通过对准来进行视觉识别的设备。
而按照本发明一优选实施例，执行上述方法的至少一个结果被采用来进行下述应用中的至少一个由航空和卫星摄影绘制地图及宇宙空间和船装配场座标测量，生产部件的座标测量(CMM)，生产部件的自动光学检测，无人操纵单元对准，无人操纵轨迹识别，3D自动操纵反馈，景象的3D模拟，对象的3D模拟，反相工程，和3D数字化。
用于识别的代数运算过程下面将参照

图1-5来作说明。
附图的简要说明结合附图，从下面更详细的说明中将理解本发明。
图1表示在有噪音的情况下表射偏振交叉方法的性能与三线性函数方法的性能之间比较的2个示图，其中点线表示表射偏振交叉法，而实线表示三线性函数法；
图2表示2个模型三维景物的视图及第三重投影视图；图3为利用三线性结果的重投影示图；图4为利用表射偏振线的交叉的重投影示图；图5为利用视图线性组合方法的重投影的示意图；图6为根据本发明一优选实施例构造和运行的3D景物再现设备的简化的功能模块图，其中该设备用于从一3D景物的至少三个视图生成该3D景物的3D表示；图7为根据本发明一优选实施例运行的一优选3D景物再现方法的简化流程图，其中该方法用于图6的设备；图8为根据本发明一优选实施例构造和运行的图象变换设备的简化的功能模块图，该设备用于从一3D景物的至少2个参考视图生成该3D视图的一新视图；图9为根据本发明一优选实施例运行的一优选图象变换方法的简化的流程图，该方法用于图8的设备；图10为也示于图2中的鞋的3D再现的示意图；图11为用于确保零件质量的优选方法和设备的说明的简化方块图；及图12为用于产生数字地图的优化方法和设备的简化方块图。
在此还附有助于理解本发明优选实施例的附件，即附件A图6所示的三线性计算单元50及表射偏振几何生成单元60的优选操作的“C”语言的表列；及附件B根据本发明的一优选实施例构造和运行的3D再现设备的优选软件操作的表列。
优选实施例的详细说明一针孔照相机，如35mm静物摄影机或摄象机，产生被观测的三维(3D)世界的二维(2D)投影。所得图象可按几何和光度标准加以分析。几何标准是指在3D中特征细节(点，线)的位置与它们各自在2D图象中的位置之间的几何关系。光度标准是指景象与图象中的亮度(象素的灰度值)之间的辐射度(表面的反射特征，照亮景物的光源的光谱特性等)关系。
3D世界被模拟化成一群点，2D世界也是如此。换言之，人们辨认图象中对应于3D中点特征的点特征。在此二组点之间存在着几何关系。当摄象机在3D中运动时人们获得同一3D世界的(亦即相同的3D点群的)更多的图象，而相关的问题就在于对应的2D图象点组与3D点组之间的几何关系。换句话说，如果Pi表明一组3D点，而pi，pi′，pi″为三组(在三个图象中)被安排成带有同一附标i的点对应于同一3D点Pi的2D点的话，则这组图象即单独地揭示大量此3D组的信息。
这里有二个要关心的问题(i)为恢复3D点的位置(再现)哪些特征必须在2D图象中加以测定(多少点跨越多少图象)；(ii)已知图象对应关系，能否合成新的视图(图象变换问题)而无须明确重现对象。本发明即为处理这两个问题。
3D至2D的总体内容是Computer Vision(计算机视觉)中的一活跃范围(运动结构和主体观测)，并且是摄影测绘工程中的唯一工具(由航空和卫星照片描制地图，宇宙空间和船坞装配场地中的座标测定等)。
3D点与2D点之间的关系可被表示成一3×4变换矩阵p≈MP式中p＝(x，y，1)T表示齐性座标中的2D点(亦即图象平面为2D投影面)，其中x，y为图象中点的位置的水平和垂直分量(相对于某任意的原点，例如图象平面的几何中心)；p由一4-向量表示，这是指将对象空间看成为被包围在三维的投影空间内。矩阵M为3×4，代表摄象机参数(例如，以MQ＝0定义的向量Q为摄象机的投影中心-所有光线的会聚点)。符号≈是指按比例相等。
因为是将图象和对象置入投影空间内(这就是说这些可经受刚性的运动和能被伸展及剪切)所以能将对象空间进行变换来使得与第一视图相关连的摄象机变换成为[I，O]，其中I为识别矩阵。对于三个视图，有
p＝[I，O]Pp′≈[A，υ′]Pp″≈[B，υ″]P式中P＝(x，y，l，k)T，图象点为p＝(x，y，1)T，p′＝(x′，y′，1)T和p″＝(x″，y″，1)T。应指出，在图象中测量的是p、p′、p″的位置，P的位置(即K的值)是未知数。A和B是3×3的矩阵(相互不是独立的)，而v′、v″为共同描述三个视图的摄象机参数的3-向量。应指出，摄象机参数A、B、v′、v″是未知的。换句话说，给出三个带有对应的点pi，pi′，pi″(i＝1…n)的视图，就存在一组数Ki，和参数A、B、v′、v″使得在上述等式对于所有i均成立。
I.三线性张量以下式描述的27数字阵列αijk=vi′bjk-vj″aik]]>i，j，k＝1，2，3式中ars、brs为A，B的元素，而vr′、vr″为v′、v″的元素，该阵列被称之为“三线性张量”。这些数字对3D和2D世界的特定投影表示均是不变的，亦即它们是此三视图所固有的(这是张量的普遍特性之一，它们不取决于为表示所选择的依据)。
应指出，阵列αijk提供作许多应用的基础。下面将说明由图象测量值p、p′、p″产生阵列的优选方法。A.由图象测量值恢复αijk一适当的三数组p，p′，p″满足多个三线性关系。第一标注以固定某些下标而改变其他来区别向量和矩阵。例如，αijk是一组标量，αij是一由9个向量组成的组(K变化而i，j维持固定)，αi是一由3个矩阵组成的组(α1…，α2…和α3…)等等。
对一给定的顺序(即不置换视图)具有9个三线性函数，以4个线性地独立的公式的组出现。这些公式的系数为阵列αijk的元素。例如以下4公式即是三个视图的线性独立的可变函数x″α13T.p-x″x′α33T.p+x′α31T.p-α11T.p=0,]]>y″α13T.p-y″x′α33T.p+x′α32T.p-α12T.p=0,]]>x″α23T.p-x″y′α33T.p+y′α31T.p-α21T.p=0,]]>y″α23T.p-y″y′α33T.p+y′α32T.p-α22T.p=0.]]>其于的函数在这里提出的等式5-11中说明。
由于每一对应的三个值组p，p′，p″提供四个线性独立的等式，所以七个跨越三视图的对应点唯一地确定(按比例)张量αijk。
注1这一点的意思就是，如果要识别7个跨越三个视图的对应点的三组值-例如说对象是一人的面孔而跟踪跨越三视图的鼻棱其在第一视图中的位置表示为(x，y)，在第二视图中的位置表示为(x′，y′)和第三视图中鼻的位置为(x″，y″)-则每一这种三组值即提供所试图获得的27个参数的4个线性独立的等式。这样即得到大于足够求解αijk需要的28个等式。
注2为线性系统求解可有许多方法，特别是在提供的等式多于未知数时(亦即在给出多于7个对应的三组值时)，例如最小二乘方技术，和为在最小二乘方技术中消除外层的增强法。
II.秩4约束跨越多于3个视图的张量连接给出一组m＞3的视图，则整个视组图的各个3视图子集的三线性张量不是不独立的，而事实上是相互作线性相关的。这种关系下面将称之“秩4约束”。
这里将三视图的张量安排为一标以G的9×3矩阵中的张量，其中列j＝1，2，3含有元素α1，j1，α1，j，2，……α3，j，3。
假定将任一组m视图安排成为ψ1，ψ2，ψj，j＝3……m。考虑视图<ψ1，ψ2，ψj>的m-2张量并使对应G矩阵被标以Gj。则Gj矩阵对-9×3(m-2)矩阵的链接不管m如何均有秩4(通常，如果Gj矩阵是独立的，链接本该是较小维次，即当m≥5时为9)。
秩4约束意味着链接矩阵四个最大的主要元素表示视图ψ1和ψ2之间的几何结构在某种程度上为静态最佳状态。A.应用1图象变换例如二个视图带有对应点pi，pi′，i＝1，…n，希望生成此同一对象的新的视图。已知张量aijk，将pi，pi′值插入三线性公式以提取x″，y″的值，而不管获得αijk的方式。
例如只要知道7点的pi″(如i＝1…7)的位置就可得到张量αijk。然后对于各个第8点就能由p，p′对和此张量来求得p″。B.应用2再现三视图的张量或m＞3视图的链接矩阵的四个最大的基本元素(见段II)可被用来恢复二视图之间的表射偏振几何关系，并由此被用来再现景象的3D模型。
假定Ej＝αij为3个3×3矩阵，亦即Ej＝α1，j，1，α1，j，2，…α3，j，3。这一将27数字阵列变成为三个矩阵的特定置换将产生下列结果此三矩阵E1、E2、E3为与3D空间中某三个平面相关连的投影变换(单应性矩阵)。
这就是说，存在有空间中这样的三个平面πj，即如果p与π1(例如)为共面的而p，p′为第一2视图中的图象位置，则E1p≈p′。对j＝2，3也同样。
第一二视图间的“基本矩阵”F(或摄影测量文献中所采用的“共面性约束”)满足EjTF+FTEj＝0 j＝1，2，3一个用于求解F而提供线性系统等式的方程。
同样，m＞3的链接矩阵的四个最大基本元素(第II部分)可被用于如下面所述求解F。假定Ak(k＝1…4)为此链接矩阵的四个最大基本元素，每一个被安排在一3×3矩阵中(行向)。则AkTF+FTAk=0]]>k＝1，2，3，4这提供由m＞3视图线性求取F的静态最佳方法。
“表射偏振”点v′可由方程FTv′＝0或直接由张量如下这样求取三个矩阵E1、E2、E3中的每二个矩阵的相应列的交叉乘积提供一表射偏振线，即一与v′相一致的线。同样，如果以a1、a2、a3标示一个E矩阵的列和以b1、b2、b3标示另一E矩阵的列，则ai×bj-bj×aj即为一表射偏振线(i≠j)。综合张量即得到18个其交点定义表射偏振v′的表射偏程线。此交点可采用例如最小二乘方法来求得。
同样，此表射偏振点v′可如下述由A1…A4，四个最大基本元素链接矩阵来求取。假定a1、a2、a3为矩阵A1的列，b1、b2、b3为A2的列，则ai×bj-bj×aj为一表射偏振线，i≠j，而ai×bi提供另一组的三表射偏振线。由这样来选择四个矩阵A1…A4中的二个，就得到提供为解v′的冗余等式组的24个表射偏程线。
此3D景物的投影模型遵循F，而v′(例如3×4矩阵[(v′)F，v′])是由第二视图的3D座标到图象座标的摄象机变换。换句话说，对第一和第二视图各自的每一对匹配点p，p′均具有p′≅[v′]Fp+kv′,]]>这样，座标向量[x，y，l，k]即提供一景象的3D(投影)模型。
注1置换αi…也生成三个其他平面的三个单应性矩阵W1，W2，W3。单应性矩阵由最先的二第三视图得到，即对于来自与W1相关的平面共面的P的p，p″的W1p≈p″。从而与矩阵E1、E2、E3的应用相似，能获取与第一和第三视图之间的表射偏振几何条件相关的基本矩阵F″和表射偏振点v″。实现跨越3视图的图象的对应性。
由2D成象来处理3D景象中的一个重要步骤是就近取得跨越视图的可靠的(仅可能密集的)匹配点组。现有技术的当前状态是采用二视图的图象灰度间的相关技术连同几何图形信息，如表射偏振几何结构。
采用张量可以将一致性方法延伸来同时利用三个视图而无须假定摄象机定标。给定一组初始的跨越三个视图的对应点就得到此三线性张量。然后，此三线性张量提供用于采用相关方法来准确地定位匹配点的制约。这可按如下完成。
假定E1、E2、E3标明张量的设置a.j.和W1、W2、W3标明张量的设置ai..。首先将所有三个图象对准使得匹配点之间的位移仅仅只沿着表射偏振线。这按如下来完成。
假定E＝∑jpjEj，这里系数pj为对第一和第二图象中各自的所有初始匹配点pk，pk′的线性最小二乘方条件Epk≈pk′的解。同样，假定W＝∑iμiWi，这里系数μi为对第一和第三图象中各自的所有初始匹配点pk、pk″的线性最小二乘方条件Wpk≈pk″的解。以变换E-1“变形”图象二和以变换W-1“变形”图象三的处理产生二个标以2′和3′的新图象，从而使得图象1，2′间的匹配点沿着表射偏振线，图象1，3′间的匹配点沿着表射偏振线。
表射偏振线可由张量推导得，如这里所介绍的(优选实施例的详细说明中的节IIB)。假定W′和W″分别为1，2′和1，3′之间的表射偏振线的方向向量。因此，如果p＝(x，y)，p′＝(x′，y′)，p″＝(x″，y″)分别为图象1，2′，3′中的匹配点，则对于某些系数α，β，p′＝p+αW′及p″＝βW″。应指出，W′、W″、α、β为随位置(x，y)改变的量。
这里表明的三线性等式(3)-(6)提供对系数α、β的几何条件约束，亦即对α的每一任意值，β值即被唯一地确定，或反之。为了取得α，β两者，如下所述利用1，2′，3′的图象亮度。
系数α，β满足“恒定亮度方程”α(Ix,Iy)Tω′+It12′=0]]>β(Ix,Iy)Tω″+It13′=0]]>这里(Ix，Iy)是第一图象中位置p＝(x，y)处的梯度向量，It12′是图象1，2′间的时间导数，和It13′为图象1，3′的时间导数(所有导数均为离散近似)。由这里提出的三线性方程(3)-(6)可得到f(α，β)＝0，这里函数f( )的系数为张量的基元。由此而取得对α和β的最小二乘方解，满足几何条件和光度测定双方的约束。
此整个的过程可按拉普拉斯金字塔(Laplacian Pyramid)分层次地进行，如J.R.Bergen和R Hingorani的“基于运动的层次帧速率变换”(Technical report，David Sarnoff Research Center，1990)中所描述的。
现在参看图6，它是按照本发明一优选实施例构造和运行的3D景象重现设备的简化功能方框图，此设备的运行是由一3D景象的至少三个视图产生此3D景象的3D表示。
图6的设备包含为由至少3个相应的观察点提供一3D景象、即对象的至少3个相应的观察点提供一3D景物、即对象的至少3个数字图象的设备，例如一由3个不同观察点运行的CCD摄象机10，或者如一可以是机载的电影摄影机20结合一象Zeiss PS-1那样的用来作将由电影摄影机20所产生的图象数字化运行的扫描器30。
此至少3个数字视图被馈送到匹配点搜索器40，它由此3个数字视图识别至少7个而最好是多个匹配点的三组值。来自不同的视图的“匹配点”为对应于真实的3D世界中一单个位置的点(地点位置)。这些点例如说可以是作人工识别的。对于航空摄影可以采用市场供应的匹配点软件，例如执行相应功能的Inpho(Stuttgart)提供的Match-T程序包。
求取匹配点的普通方法在Wang，H.和Brady，J.M.的“棱角检测一些新的结果”(IEEE colloquium Digest of System Aspects of Machine Perception and vision”，Londone 1992，pp.1.1-1.4)中有介绍。
三线性张量计算单元50接收此匹配点三值组并计算表示三个视图之间的几何参数关系的三线性张量。
然后此三线性张量被用来产生3D景象的3D表示。按照本发明一实施例，此三线性张量由一表射偏振几何关系生成单元60利用来计算此三视图中二个的表射偏振几何表述。接着，3D表示产生单元70按单元60输出的表射偏振表述产生景象、即对象的3D表示，如下述文献中详细描述的a.Fangeras，O.D.“利用未作标定的立体设备能在三维空间看到什么？”，Proceedings of the EuropeanConfernce on computer vision，pp.563-578，Santa Margherita Ligure，Italy，June 1992。
b.Hartley，R.等“由未作标定的摄象机得的立体效果”，Proceedings IEEE Conf.on Computer Vision and Pattern Recognition，761-764，Champaign，IL，June 1992。
c.Shashua，A.“由未标定图象得的投影结构由运动和识别得的结构”，IEEE Transactions on PAMI，16(8)778-790，1994。
单元50和60的优选实现方案在附录A中以“C”计算机语言加以描述。单元50的优选实现方案由附录A的P1直至P4末尾中说明。单元60的优选实现方案则在附录A中由P4末尾到P8中部说明。有助于理解上述材料的子程序和统计过程列于附录A的P8以下。
包含有表示3D景象、即对象的至少一部分或一个方面的信息的3D表示可由一激光计算机打印机用来产生该对象(景象)的一新的视图。或者，通常的CAD(计算机辅助设计)软件与普通绘图机相结合也可被用来产生对象(景象)的新的视图。CAD软件还可在质量保障应用中用于比较2个CAD文件。
现在来参看图7说明一优选3D景象再现方法的简化流程图，按照本发明一优选实施例运行的这一方法有助于与图6的设备相结合。图7总体上是足以说明自身的。如果景象(对象)的3个视图是由一数字案卷取得或者是由例如一数字摄象机产生的，它们就可以是数字式的。如果它们不是数字的，就可对它们进行扫描或者进行数字化。
任何合适的一般的或其他格式均可被采用。例如Silicon Graphics′s Inventor格式一开始就可被用于3D表示。此Inventor格式可被变换成Postsript以便打印3D表示的一新视图。
如所示，景象的3D表示可用于进行广泛范围的活动中，例如景象、即对象的3D测量，如由打印输出来生成景象(对象)的新视图，和质量保证比较，后者是将所生成的对象(景象)的3D表示利用一普通方法与所希的对象(景象)或其所希望的3D表示进行比较。
现在参看按照本发明一优选实施例构成和运行的图象变换设备的简化功能方框图的图8，该设备作由一3D景象的至少二参考视图生成此3D景象的新视图的运行。图8的设备与图6的设备类似。但图8中一3D景象的新视图是直接仅由其至少二个参考视图产生的，最好无需中间产生一3D表示。
在图8中，有关此二参考视图和所希望的新视图的几何信息被用来产生表示三个视图之间的几何关系的三线性张量。一般，可例如人工地识别此新视图中至少7个三值组的位置，分别对应于至少二参考视图的每一个中的7个位置。通常至少某些有关此新视图的信息是可利用的。例如，如果希望根据相同地区的至少二新的基础视图更新一GIS(地理信息系统)的年前的视图的话。一般可能识别对应于该二个参考视图中的7个位置并被认为依然存在于此年前视图的即将被产生的当前型式中的此年前视图中的至少7个位置。
此新视图一般由根据所述第一和第二参考视图中的第一和第二相应位置和所述三线性张量分别计算多个各自对应于不同的所述第一和第二相应位置的新视图位置来产生。例如，匹配点搜索器140可由此二参考视图产生多对的匹配点，例如说1000对。对于最先的7对，用户可人工地指明此新型视图中的7个匹配点。
一旦由三线性张量计算单元150生成三线性张量后，就可将匹配点对8-1000的座标插入此三线性张量，如图9中所示，以便生成此新视图中的匹配点8-1000的座标。
这样生成的新视图可被例如与一年前的同一视图进行比较，来辨别这一年的过程内该景象中所发生的差异。
现在参看说明按照本发明一优选实施例运行的优选图象变换方法的简化流程的图9，此方法可用于与图8的设备相结合。图9通常是自说明性的。
在这里所介绍的任一实施例中，如果提供有多于3个的视图，可对每三个视图计算“中间”张量。然后可根据这些“中间”张量间的关系计算一具有代表性的张量。
图10为图2中也曾作图解的鞋的3D再现的图示。图10是由求得匹配点和再现它们的3D位置所生成的。接着将座标用CAD软件处理来生成图10中所示的表面。
图11为说明作工作质量保证的优选方法和设备的简化方框图。图11包含一3CCD摄象机210的阵列200，这些摄象机均对准一单个位置由此来生成该位置的三个透视图。这些摄象机被装附到一机器人臂212，因而能按照来自控制器224的适当指令作相对于沿着传送带220到达的工件214的运动。
当一工件进入CCD摄象机210的视野时，传送带通常会作暂停来使得基本上整个的工件表面区域均能由摄象机210成象。控制器224通过机器人臂212围绕对象移动摄象机阵列200，以使得一般几乎整个表面区域均被成象。例如可将摄象机阵列200围绕对象移动通过10个不同的位置，而在每一个位置3个CCD摄象机中的每一个均对工件摄象。所采用位置的数量取决于工件的复杂性。
这一过程由三个相应透视图生成多个三值组图象，每一三值组图象包含有该工件同一部分的三个数字图象。对各个图象三值组，可根据已知的机器人臂的位置和采用手-眼定标来计算阵列200的和每一摄象机210的对应位置。
各图象三值组由可能类似于图6的单元40、50、60、70和90的单元240、250、260、270和290分别进行处理。
由CAD S/W290按各个图象三值组(各取样位置)所产生的CAD模型信息被存贮在适当的存贮器300中。一计算单元310用来将多个对应于多个CCD摄象机阵列200的位置的取样地点汇集进一单一的座标系统。所需的座标变换由将确定CCD摄象机阵列的运动的反变换来进行计算。
计算单元320将单元310的输出与一参考CAD模型进行比较计算其间的差。这些差在计算单元330中与可接受的容差相比较。
图11的设备也适应于作反向工程或CAD(计算机辅助设计)中的3D数字化应用。
图12为说明为生成数字地形图的优选方法和设备的简化方框图，例如为了更新都市地图，为检查非法建筑，为提供汽车导航系统服务，或者甚至为制作诸如集成电路等的微观对象的图形。
机载CCD摄象机334飞行越过希望对其产生数字地形图的景象。此摄象机334由三个相应的透视图来生成此景象的3个数字图象。此3个图象由可能类似于图6的单元40、50、60和70的340、350、360和370进行处理。
由一表面插补单元380对3D显示生成单元370的输出进行一表面插补处置。一适宜的表面插补方法在Grimson，W.E.L.的“视觉表面插补的计算理论”(Proc of the Royal Soc.of London B.298，pp.395-427，1982)中有介绍。
最后由一正射照片生成单元390从表面插补单元380的输出去除透视失真，而得到一正射象片(平面地物图)。单元390的适当的运行方法在Slama，C.C.的“摄影测绘手则”(Amercan Society ofPhotogrammetry，1980)中有描述。
本发明还有助于作3D对象的显现，特别是对于娱乐、动画或教育上的应用。象图11的阵列200那样的摄象机陈列围绕一待加以显现的对象旋转而基本上将希望显示给用户的整个表面区域成象。例如摄象机可在围绕该对象的200个位置的每一个处对该对象成象。
于是可生成不是上述200个位置处的合成图象。用户可利用遥控操纵杆来指明所希望的位置。本发明的设备可被用来为该位置生成一合成图象。
普通的驾驶模拟游戏采用合成背景，但本发明可被用来提供带有真实背景的驾驶模拟游戏。为此，使至少3个最好5-10个摄象机的阵列在所希望的景象中运动以使得此摄象机阵列能由至少3个而最好更多些的不同透视图俘获基本上整个的景象。例如，此景象可由摄象机阵列的接近1000个位置中的每一个来俘获。然后按照本发明产生新的视图来适应用户例如以一控制杆指明新视图的需要。
现在根据一篇将在《IEEE，Transactions on PAMI》上公开的标题为“用于识别的代数函数”的文献说明有助于识别的代数函数。
用于识别的代数函数I.导言建立有关一3D景象的三个透视图之间的代数关系的通常结果和论证其通过对准对视觉识别的应用。一般表明，一景象的任何三个透视图能满足图象座标的一对三线性函数。在极限情况中，当所有三个视图均为正射的是，这些函数就成为线性的并减化到[38]所揭示的形式。利用此三线性结果就可处理一对象的视图(例如由二模型视图生成新视图)而无需恢复景象结构(量度的或非量度的)，摄象机变换，或者甚至表射偏振几何参数。而且，三线性函数可依靠线性方法7点的最小结构来求取。后者被证明是为对将一3D景象重新投影到对应于二参考视图之间的点的给定的任一新视图上的问题的通用线性解所需的最小结构上的新的下边界。前面的解依靠恢复表射偏振几何参数，这又需要对一线性解的8点的最小结构。
中心结果包含在Theorem(定理)1、2和3中。第一定理说明，由未标定的针孔摄象机得到的一固定的3D对象的各种视图ψ满足分类关系F(ψ，ψ1，ψ2)＝0，其中ψ1、ψ2为对象的二个任意视图，F具有一特殊的三线性形式。F的系数可线性地求得无须首先建立对象或摄象机运动的表射偏振几何关系的3D结构。为证明Theoreml所需的辅助定理可能自身是重要的，因为它们建立平面投影变换间的某些规则并导入新的视图不变量(Lemma4)。
Theorem2针对以最经济的方式获取三线性函数的系数的问题。已表明在三个视图之间的所有可能的三线性函数中，最多存在四个线性独立的这种函数。结果，这些函数的系数可由三个视图之间的7个对应点线性地得到。
Theorem3是Theorem1显然的推论但含有重大的实用特性。已表明如果由并行投影得到视图ψ1、ψ2，则F简化到一特殊的双线性形式，或者相当地，任何透视图ψ都可由二正射视图的有理线性函数得到。简化到双线性形式意味着如果被存贮在存贮器中的二参考视图(模型视图)为正射的就有可能采用较简单的识别方法。
这些结果可能有数个应用(在节VI中讨论)，但在整个这一说明中着重的一个是用于通过配准作3D对象的识别任务。作识别的配准方法([37，16]及其中的参考文献)是依据这样的结果，即经受3D刚性仿射或投影变换的对象的同样等级的视图(不计自身的闭塞)的俘获可借助存贮对象的3D模型，或者仅仅依靠存贮对象的至少二任意的“模型”视图一假定模型视图的对应问题可在一定程度上得到解决(参看[27，5，33])。识别期间，新输入视图与一特定候选对象的模型视图之间少量的对应点就足以将模型“再投影”到新观察位置上。如此被再投影的图象成功地与输入图象匹配即完成识别。这种利用有限数量的对应点由一组模型视图预测一新视图的问题被称之为再投影。
此再投影问题原则上可通过形状和摄象机运动的3D重现来处理。这包含为求取刚性摄象机运动参数和量度形状由运动方法得的经典结构[36，18，35，14，15]和为求取非量度结构的较近代的方法，即假定对象经受3D仿射或投影变换，或同等地，摄象机未经标定[17，25，39，10，13，30]。用于透视图的经典方法在图象测量中的误差情况，窄的视野，和内部摄象机标定方面不稳定是已经公知的[3，9，12]，因此，对于再投影的目的多半未得到实际上的应用。总体观念上的非量度方法还没有对真实图象作完全的测试，但到目前为止所建议的方法依赖于首先求取表射偏振几何参数-一当存在噪声时其不稳定也是公知的过程。
也已看到，表射偏振几何参数独自就足以依靠表射偏振线的相交[24，6，8，26，23，11]利用三视图之间的至少8个对应点来实现再投影。但这只有在三摄象机的中心是非共线的-除非这些中心远不是双线性否则会导致数量的不稳定性-才可能，而且三焦点平面上的任一对象点也不能被再投影。而且，在采用非量度再现方法时，获得表射偏振几何参数即使在利用数十个对应点而在现有技术方法的状态下充其量也不过是一灵敏的过程(见第V的详述和以模拟的与真实的图象进行的比较分析)。
因此，为了稳定目的，值得进行探索为实现再投影的更直接的手段。例如，替代形状和不变量的再现，可以建立单独以图象座标的函数表示的视图之间的直接关联，称之为“视图的代数函数”。[38]已得到正射情况下的这样的结果。它表明一对象的任意三个正射视图满足对应图象座标的线性函数-在此将仅表明更大组的代数函数的限制情况，它们一般均具有三线性形式。利用这些函数就可以处理对象的视图，例如生成新的视图，而无需作为中间步骤获取形状或摄象机几何参数，全部需要做的只是适当地组合二参考视图的图象座标。而且，采用这些函数，表射几何参数被缠绕在一起，不仅使得没有奇异性和点的最小组构上的下边界，而且在实验章节中还可看到图象测绘中存在差错时的更正确的性能。这一部分的成果(仅Theorem1)简要地在[31]作了介绍。
II.符号标记现考虑三维投影空间P3的对象空间和二维投影空间P2的图象空间。设P3为一组代表3D对象的点，和ψiP2表示以i作下标的的视图(随意的)。假定二摄象机其中心分别位于O，O′∈P3，表射偏振极即被定义为在线OO′与二图象平面的交点。因为图象平面是有限的，所以能不致丧失普遍性地给每一观测的图象点指派值1作为第三齐座标。就是说，如果(x，y)为某点的被观测的图象座标(相对于某一随意的原点-例如说图象的几何中心)，则P＝(x，y，1)就是指此图象平面的齐性座标。应指出，这一约定忽略掉其中内的点在这些视图中是无限的那些特殊的视图-这些奇异情况这里未加模拟。
由于这里将以最多每次三个视图讨论，所以相关的表射偏振极作如下指定假定v∈ψ1和v′∈ψ2为视图ψ1、ψ2之间的对应表射偏振极，假定v∈ψ1和v″∈ψ3为视图ψ1、ψ3之间的对应偏振极。同样，三个视图之间的对应图象点将被指定为P＝(x，y，1)，P′＝(x′，y′，1)和P″＝(x″，y″，1)。术语“图象座标”将是指P2的非齐性座标表示，例如用于三个对应点的(x，y)，(x′，y′)，(x″，y″)。
平面以πi表示并加以下标i，而如果仅讨论一个平面就是π。所有平面都被认为是随意的并且相互区分。符号≈是指明按比例相等，GLn代表n×n的矩阵组，和PGLn是按比例定义的组。
III.三线性公式这一说明书的中心结果在以下二定理中表明。本章其余部分用于证明这一结果及其蕴涵。
Theorem1(Trilinearity，三线性)假定ψ1、ψ2、ψ3为以一组3D中的点模拟的某一对象的三个随意的透视图。三个视图之间的三个对应点的图象座标(x，y)∈ψ1，(x′，y′)∈ψ2和(x″，y″)∈ψ3满足一对以下形式的三线性方程x″(α1x+α2y+α3)+x″x′(α4x+α5y+α6)+x′(c7x+α8y+α9)+α10x+α11y+α12＝0，和y″(β1x+β2y+β3)+y″x′(β4x+β5y+β6)+x′(β7x+β8y+β9)+β10x+β11y+β12＝0，式中系数αj，βj，j＝1，……12，对所有点均为固定的，被按整个比例唯一地确定，和αj＝βj，j＝1，……6。
下面的辅助定理被用作证明的一部分。
Lemmal(Auxiliary-Existence，辅助-存在)假定A∈PGL3为因某平面π产生的投影映射(单应性)ψ1→ψ2。假定A被定标得满足Po′≈APo+v′，其中Po∈ψ1和P′o∈ψ2为来自一随意的点Po
的对应点。这样，对于来自随意点P∈P3的对应点对P∈ψ1和P′∈ψ2，即具有p′≅Ap+kv′]]>系数K与ψ2无关，亦即对第二视图的选择是不变量。
此定理、其证明及理论和实际的涵义在[28，32]中有详细讨论。应指出，在单应性A是仿射的和表射偏振v′是在无穷远的线上时的特定情况，对应于由二正射视图来的仿射结构的构成[17]。概括地说，P3的表述Ro(座标四元组)可总被选择得使一随意平面π为无穷远处的平面。这样，一总体未定标摄象机运动产生能由仿射组的元件表明为与Ro相关的表述R。因此，标量K是一投射框架中的仿射不变量，并被叫做相对仿射不变量。二个这种不变量各自对应一不同的参考平面之比就是投影不变量[32]。在此无需讨论取得K的方法，所有需要作的只是应用已存在的与某参考平面相关连的相对仿射不变量K，其进而导致正射A。
Definition(定义)1，因同一平面π产生的由ψ1→ψ1单应性Ai∈PGL3被说成是标定兼容的，如果它们被定标得满足Lemma1的话，亦即，对于任何投影到P∈ψ1和P′∈ψi上的点P∈均存在一标量K，它对任何视图ψi满足p′≅Aip+kvi,]]>其中vi∈ψ1是带有ψ1的表射偏振极(随意标定的)。Lemma2(Auxiliary-Unqueness，辅助-单值性)假定A，A′∈PGL3分别为因平面π1、π2产生的二ψ1→ψ2的单应性，则存在一标量S，它满足等式A-sA′＝[αv′，βv′，γv′]，对于某一系数α，β，γ。
证明假定q∈ψ1为第一视图中的任一点。存在有一标量Sq，它满足v′≈Aq-SqAq。令H＝A-SqA′，即得到Hq≈v′。但如[29]中证明的，对于因任一平面引起的任何单应性ψ1→ψ2，Av≈v′。因此也存在Hv≈v′。仅在H为因经由平面(与投影中心O共面)产生的单应性时才会发生二分开的点q、v映射在同一点v″上，这样，对于所有p∈ψ1Hp≈v′，Sq为一固定的标量S，。后者也就意味着H为一其列为v′的倍数的矩阵。Lemma3(Auxiliary for Lemma4，对Lemma4的辅助)假定A，A′∈PGL3分别是因不同平面π1、π2引起的由ψ1→ψ2的单应性，和β，β′∈PGL3分别为因不同平面π1、π2引起的由ψ1→ψ3的单应性。则对某些T∈PGL3，A′＝AT，和B＝BCTC-1，其中Cv≈v。证明假定A＝A2-1A1，其中A1、A2分别为由ψ1ψ2到π1上的单应性。同样，B＝B2-1B1，其中B1、B2分别为由ψ1、ψ3到π1上的单应性。假定A1v＝(c1、c2、c3)T和假定C≈A1-1diag(c1、c2、c3)A1。则B1≈A1C-1，这样得到B≈B-12A1C-1。应指出，A1与B1之间的差是因表射偏振极v，v的不同位置引起的，它由C(Cv≈v)补偿。假定E1∈PGL3为由ψ1到π的单应性，和E2∈PGL3为由π2到π1的单应性。则利用正常的定标E1和E2，得到A′=A2-1E2E1=AA1-1E2E1=AT,]]>和用正常的定标C，得到B′=B2-1E2E1C-1=BCA1-1E2E1C-1=BCTC-1]]>Lemma4(Auxiliary-Uniqueness)、对定标兼容的单位性，Lemma2的标量s、α、β、γ为由ψ1、π1、π2指明的不变量。即就是，给定一随意的第三视图ψ3，假定B，B′分别为因π1、π2生成的由ψ1→ψ3的单应性。假定B为与A定标相兼容的，B′为与A′定标相兼容的。则，B-sB′＝[αv″，βv″，γv″]。证明首先证明s为一不变量，亦即，B-sB′为其列为v″的倍数的矩阵。按Lemma2和Lemma3存在有一其列为v′的倍数的矩阵H，一满足A′＝AT的矩阵T，和一使得I-sT＝A-1H的标量s。在两边乘以BC之后，和再予乘以C-1，就得到B-sBCTC-1＝BCA-1HC-1按Lemma3有B′＝BCTC-1。矩阵A-1H具有列为v的倍数的列(因为A-1v′≈v)，CA-1I是一其列为v的倍数的矩阵，和BCA-1H是其列为v″的倍数的矩阵。将BCA-1H乘以C-1并不改变其形式因为BCA-1HC-1的每一列仅只是BCA-1H的列的线性组合。结果，B-SB′就是一其列为v″的倍数的矩阵。
假定H＝A-sA′和H＝B-sB′。由于单应性是定标兼容的，按Lemma1，存在有与一随意的P∈ψ1相关的不变量K、K′，其中K是因π1产生的，K′是因π2产生的p′≈Ap+Kv′≈A′p+K′v′和p″≈Bp+Kv″≈B′p+k′v″。而按Lemma2具有Hp＝(sk′-k)v′和Hp＝(sk′-k)v″。由于p是随意的，这仅有在H中的v′的倍数的系数与H中的v″的倍数的系数相一致时才会发生。Proof of Theorem(定理的证明)Lemma1提供定理的存在部分，如下述。由于Lemma1适用于任何平面，选择平面π1并假定A、B分别为定标兼容的单应性→ψ2和ψ1→ψ3。则对每一点P∈ψ1在对应点为p′∈ψ2和p″∈ψ3时，存在有满足p′≈Ap+Kv′和p″≈Bp+kv″的标量k。可由二式提取k并得到
其中b1、b2、b3和a1、a2、a3为A和B的行向量，v′＝(v1′，v2″，v3′)，v″＝(v1″，v2″，v3″)。因为k的不变性，可使(1)的项与(2)的项相等以得到三视图之间图象座标的三线性函数。例如，使各等式中的起始二项相等就得到x″(v1′b3-v3″a1)Tp-x″x′(v3′b3-v3″a3)Tp+]]>x′(v3′b1-v1″a3)Tp-(v1′b1-v1″a1)Tp=0,----(3)]]>以同样方式，以使(1)的第一项与(2)的第二项相等之后得到y″(v1′b3-v3″a1)Tp-y″x′(v3′b3-v3″a3)Tp+]]>x′(v3′b3-v2″a3)Tp-(v1′b2-v2″a1)Tp=0.----(4)]]>此二式均为所希望的形式，其最初的六个系数在二等式中相同。
因为Lemma1适用于任何平面而出现唯一性问题。如果选择一带单应性A′、B′的不同平面，例如说π2，则必须证明新的单应性导致相同的系数(按整个比例)。(3)和(4)中括号内的项对某些i和j具有通用形式vi′bj-vj″ai。因而需要证明存在有一标量S能满足vj″(ai-sai′)=vi′(bj-sbj′)]]>但这是由Lemma2和4直接得出的。
此定理的直接蕴涵是仅由组合二模型视图(ψ1，ψ2)即可得到一新视图(ψ3)。组合系数αj和βj可一起作为一给定三视图之间的9个对应点的17个方程(24-6-1)的线性系统的解来求得(多于9个的点可被用于一最小二乘方解)。
在下一定理中得到为求解此二线性函数的系数所需点的数量的下边界。Theorem1证明的实体部分指明存在有该型式的9个三线性函数，其系数具有通用形式v1′bj-vj″ai。这样最多就有27个不同的系数(按统一的比例)，并因此如果此9个三线性函数中的多于二个为线性独立的，就可利用少于9个的点来求解该系数。下一定理表明最多4个三线性函数为线性独立的，并因而7个点即足以求解系数。
Theorem2存在有Theorem1中所述类型的9个不同的三线性形式，其中最多4个是线性独立的。此4个三线性形式的系数可利用三视图之间的7个对应点求得。
证明9个三线性形式的存在直接由(1)和(2)得到。假定αij＝vi′bj-vj″ai。9个形式如下(最先的2个为了方便是(3)和(4)的重复)x″α13Tp-x″x′α33Tp+x′α31Tp-α11Tp=0,]]>y″α13Tp-y″x′α33Tp+x′α32Tp-α12Tp=0,]]>x″α23Tp-x″y′α33Tp+y′α31Tp-α21Tp=0.----(5)]]>y″α23Tρ-y″y′α33Tp+y′α32Tp-α22Tp=0.----(6)]]>y″x′α31Tp-x″x′α32Tp+x″α12Tp-y″α11Tp=0.----(7)]]>y″y′α31Tp-x″y′α32Tp+x″α22Tp-y″α21Tp=0,----(8)]]>x″y′α13Tp-x″x′α23Tp+x′α21Tp-y′α11Tp=0,----(9)]]>y″y′α13Tp-y″x′α23Tp+x′α22Tp-y′α12Tp=0,----(10)]]>x″y′α12Tp-x″x′α22Tp+y″x′α21Tp-y″y′α11Tp=0,----(11)]]>对于一给定的三值组p，p′，p″，表上的最先4个函数产生一4×27矩阵。此矩阵的阶为4因为它含有4个正交列(与α11，α12，α21和α22相关连的列)，因此这些函数是线性独立的。由于具有27个系数，每一三值组p，p′，p″用于4个线性方程，这样三视图间的7个对应点即提供为线性求解系数的足够数量的方程(只要系统按一公共比例确定，7个点就产生能被用作一致性检查或用于得到最小二乘方解的二个附加的方程)。
其余的三线性形式是最先4个的线性扩展，如下
(7)＝y″(3)-x″(4)(8)＝y″(5)-x″(6)(9)＝y′(3)-x′(5)(10)＝y′(4)-x′(6)(11)＝y″y′(3)-x″y′(4)+x″x′(6)-y″x′(5)其中括号内的数表示各不同三线性函数的量程数。
共同地二定理为在给定二模型视图间的对应性p，p′求解一新视图中的位置x″，y″提供积极的手段。这一产生新视图的过程很容易完成而无需明确地了解结构、摄象机变换、或甚至表射偏振几何参数，而需要的对应点数量少于其他任一已知的方案。
对x″，y″的求解是唯一的，没有对所容许的摄象机变换的约束。不过有些摄象机组构需要与Theorem2的证明中所建议的不同的一组4个三线性函数。例如方程(5)、(6)、(9)和(10)组也是线性独立的。这样，如果v′1和v3′同时成为0，即v′≈(0，1，0)，则该组即可被代之以应用。同样，方程(3)，(4)，(9)和(10)为线性独立的，而在v≈(1，0，0)时应被加以应用。v″≈(1，0，0)和v″≈(0，1，0)时亦产生类似情况，可由上述讨论的6个函数中选择4个作为适当的依据来进行处理。应指出，这里没有涉及7点的奇异组构问题。例如很显然，如果此7点是共面的，则它们在三个视图间的对应性也许就不能为求取系数而产生唯一解。奇异表面的实质已就为寻求表射偏振几何参数所需的8点的情况进行了研究[19，14，22]。有关本文所提供的结果的同一实质是一未解决的问题。
消除对获取表射偏振几何参数的需要带来明显和巨大的优点。为更好理解这些优点，现在来粗略地考虑采用表射偏振几何参数作再投影的过程。表射偏振交叉方法可作如下的扼要说明(见[11])。假定F13和F23为满足p″F13p＝0和p″F23p′＝0的矩阵(当前词汇中“基本″矩阵[10])。则由p″与其表射偏振线的关联而得到p″≅F13p×F23p′.]]>因此，三个视图间的8个对应点就是以为二基本矩阵作线性求解，而后所有其他对象点均能再投影到第三视图上。等式(12)也为三线性形式，但不是Theorem1中引入的类型。不同点包含有(i)表射偏振交叉需要来自8点而不是7点的对应；(ii)p″的位置由一其在三摄象机中心共线的情况下是奇异的线交叉过程来求解；在三线性结果中p″的成份被分开求解而三个共线摄象机的情况是允许的；(iii)此表射偏振交叉过程是可分解的，即每次只用二个视图，而表射偏振几何参数在三线性结果中则是缠结在一起的不能独立获取。后者还意味着在存在噪声时的一种较好的数字运行方法，如下面将表明的，即使采用最小数量的所需的点，其性能也远远超过采用多得多的点的表射偏振交叉的性能。换句话说，由于避免了获取表射偏振几何参数的需要而得到巨大的实际利益，因为在协同透视图工作时表射偏振几何参数是最易于发生差错的组成部分。视图三线性函数的通用结果与对正射视图的“视图的线性组合”结果[38]之间的关系可由将p2中A和B设定为仿射的和v3′＝v3″＝0很容易看到。例如，(3)化简为v1′x″-v1″x′+(v1″a1-v1′b1)Tp=0,]]>这是一种α1x″+α2x′+α3x+α4y+α5＝0.的形式。如在透视情况中那样，每一点用于4个等式，但这里利用其中所有的4个来求取系数是不利的，因此可以仅利用4个等式中的二个和需要4个对应点来求取系数。这样，在所有三个视图均为正射的情况中，x″(y″)即被表明为一个二另外视图的图象座标的线性组合，如[38]中所揭示的。IV.双线性形式现在考虑正射地取一对象的二参考(模型)视图(采用一电视镜头将提供合理的近似性)，在识别期间对象的任一透视图均是容许时的情况。可容易地证明，三个视图将通过双线性函数(而不是三线性的)相关连Theorem3(Bilinearity，双线性)在Theorem1的条件中，如果由并行投影得到视图ψ1和ψ2，则一对Theorem1的三线性形式简化成如下的一对双线性方程x″(α1x+α2y+α3)+α4x″x′+α5x′+α6x+α7y+α8＝0.和 y″(β1x+β2y+β3)+β4y″x′+β5x′+β6x+β7y+β8＝0.其中αj＝βj；j＝1，…4。证明在这些条件下由Lemma1得到A为P2中仿射的和v3′＝0，因此(3)简化成x″(v1′b3-v3″a1)Tp+v3″x″x′-v1″x′+(v1″a1-v1′b1)Tp=0.]]>同样，(4)简化成y″(v1′b3-v3″a1)Tp+v3″y″x′-v2″x′+(v2″a1-v1′b2)Tp=0.]]>此二等式均是所希望的形式，其最先的4个系数在二等式间是相同的。
其余三线性形式均经过类似的简化，而Theorem2仍然成立，即仍然有4个线性独立的形式。结果得到21个按共同比例的系数(而不是27个)和每点的4个方程，从而5个对应点(而不是7个)就足够用于作线性求解。
三视图的双线性函数较通常的三线性函数具有二个优点。第一，如前面提到的，为求系数的解仅需要三视图间的5个对应点(不是7个)。第二，代数函数的方次越低，测量对应性中求解将越不易于出现差错。换句话说，很可能(虽然未必一定)较高阶项如等式3中的项x″x′x将对系统的整体误差敏感性的作用较大。
与认为所有视图均是正射的情况相比较，这一情况更接近真实。由于仅只取模型视图一次，所有要求采取特定方法、亦即采用电视镜头(假定处理对象识别而不是景象识别)就不是不合理的。如果满足这种要求，则就是普通的识别任务，因为此时的识别处理期间可以取任何透视图。V.试验数据进行本节所描述的试验是为了评估为再投影采用三线性结果与采用[38]的表射偏振交叉和线性组合结果相比较的实际作用(后者在这里仅表明了三线性结果的有限情况)。
节III中所述是由首先获取基本矩阵来实现的表射偏振交叉方法。虽然8个对应点足以作线性求解，实践中将采用多于8点来以线性的或非线性的二乘方方法求得基本矩阵。由于线性最小二乘方法仍然对图象噪声敏感，所以采用的是由T.Luong和L.Quan提供的[20]中描述的非线性实现方法(在[20]中建议有二种实现方法，是采用得到较佳结果的一种实现方式)。
第一个试验采用的模拟数据表明，即使正确地求取得表射偏振几何参数，采用免除线交叉处理的三线性结果仍然要好得多。第二个试验是对真实的图象组进行的，比较了各种方法的性能和实际达到合理的再投影效果所需的对应点的数量。A.计算机模拟利用了一被随机地置于Z轴座标在100单元至120单元和x，y轴座标随意地从-125至+125的范围内的46点的对象。焦点长为50单元和第一视图由fx/z，fy/z得到。第二视图(ψ2)是由围绕点(0，0，100)利用轴(0.14，0.7，0.7)和由0.3弧度的角度旋转所产生的。第三视图(ψ3)是由围绕轴(0，1，0)以同样的变换和角度旋转产生的。各种数量的随机噪声被作用到所有被再投影到第三视图上的点，但不作用到8或7个被用于求取参数的点(基本矩阵，或三线性系数)。噪声被随机地分开地加到各座标并带有由0.5至2.5象素误差的变化水平。进行了1000次如下的试验生成20个随机对象，对每一级误差各个对象进行10次模拟。采集了最大再投影误差(以象素为单位)和平均再投影误差(被作再投影的所有点的平均)。由对所有试样(它们中的200个)进行平均分别采集各个误差级别的这些数字和记录标准的偏离。由于没有误差被加到被用来确定表射偏振几何参数和三线性系数的8或7个点，所以仅进行了为得到基本矩阵或三线性系数所需的有关的线性方程系统的求解。
结果示于图1A-1B。右图表示测量最大投影误差得的各图象噪声水平的二种算法的性能。发现在所有的噪声水平下，三线性方法明显较佳而且还具有较小的标准偏离。对右边图中所示的平均再投影误差也同样。
性能中的这种差别是所预料的，因为三线性方法是一并考虑所有的三个视图，而不是每对分别考虑，并因而避免了线交叉。B.对真实图象的试验图2A-2C表示选择来作试验的对象的三个视图。此对象为带有为便于作对应性处理的附加结构的运动鞋。选取这一对象是因为其复杂性，即，它具有一自然对象的形状和不易于作参数描述(如平面或代数表面的集合)。应指出，这里所说明的情形是挑战性的因为被再投影的视图不是在二模型视图的中间，即应当预期到较之在中间情况对图象噪声更为敏感。人工地选择了帧面中的一个ψ1上的一组34点，它们的对应物自动地在此试验中所采用的所有其他帧面上得到。此对应过程是根据[7]中所描述的由粗到细的光流算法的执行。为达到分开的视图间的准确对应，采用了过渡的中间帧面和附加以连续帧面的位移。然后将全部位移场用来将第一帧面推(“翘曲”)向目的帧面这样来生成一合成图象。再次将光流作用到合成帧面与目的帧面之间并将最后得的位移加到早先得到的整个位移上。这一过程提供一密集的位移场，它随后被采样来获得最初在第一帧面中所选择的34点的对应物。这一过程的结果在图2A-2C中以显示中心围绕所计算得的对应点位置的方块来表示。可以看到，这种方式所取得的对应是合理的，和在大多数情况下达到子象素的准确度。依靠选择第一帧面中的点很容易进一步自动执行这一过程，为此空间导数Hessina矩阵可作良好地适应，类似于[4，7，34]的方案中所建议的置信度值，但这里并不企图这样来建议一完整的系统，仅只是测试三线性再投影法的性能和次其与表射偏振交叉和线性组合法的性能和比较。
三线性法需要三视图间至少7个对应点(这里需要26个方程，而7点提供28个方程)，而表射偏振交叉可以8点进行(原则上)。现在要解决的问题是为达到合理的性能实际需要多少点数(由于对应中的误差、透镜畸变和其他可能由针孔摄象机模型不适当地模拟的效果)。
三线性结果首先被以最小点数(7个)用于为系数求解，然后被应用于以8、9和10点的一线性最小二乘方解(应指出，一般利用SVD或Jacobi法替代线性最小二乘方可得到较佳的解，但这里并不想这样做)。结果如图3A-3B所示。7点提供的再投影具有最大误差为3.3象素和平均误差为0.98象素。利用10点的解取得最大误差为1.44和平均误差为0.44象素的改善。利用8和9点的性能合乎情理地在以上的性能之间。采用更多的点没有显著改善此结果，例如，当所有的34点均采用时最大误差降到1.14象素，平均误差为0.42象素。
接着应用表射偏振交叉法。为求取基本矩阵采取2种方法。一个方法是采用[20]的措施，另一个是利用来自一个平面(基础平面)的4个对应点。在前一情况中，为实现合理的结果需要比8点多得多的点。例如在利用全部34点时最大误差为43.4象素和平均误差为9.58象素。在后一情况中，首先求得因基础平面所引起的单应性B而后利用二附加点(暗盒上的点)得到表射偏振v″。这样就知道(见[28，21，32])F13＝[v″]B，其中[v″]是v″的反对称矩阵。求取F23采用了同样的过程。因此仅6点被用于再投影，但尽管如此结果还是稍有改善最大误差为25.7和平均误差为7.7象素。图4A-4B表明这些结果。
最后测试了利用线性组合方法再投影的性能。由于线性组合方法仅适用于正射视图，这里实际上是在透视情况下测试正射假设，或者换句话说，即三线性方程的高阶项(二线性和三线性)是否有意义。线性组合方法需要至少三视图间的4个对应点。应用此方法是利用4、10(为了与图3A-3B中所示的三线性情形比较)和全部34点(后两者利用线性最小二乘方)。结果显示在图5A-5C中。所有情况中的性能均大大差于采用三线性函数时，但比表射偏振交叉法为佳。VI.讨论已经看到，一般情况下一固定的3D对象的任一视图均能用二个参考视图表示为一三线性函数，或者在参考视图由平行投影所产生时表现为一双线性函数。这些函数提供为处理一景象的视图可选择的较之其他方法简单得多的措施。而且它们在理论上需要较少的对应点，而在实践中还少得多。试验结果表明三线性函数实践中还有利于产生较之表射偏振交叉或线性组合方法显著要好的性能。
一般二个视图使一“基本”矩阵(见[10].)能表述此二视图之间的表射偏振几何参数，且该矩阵的阵元受主体约束(矩阵的秩为2)。三线性结果(Theorem1，2)意味着三个视图可以有一带有27个不同的元素的张量。已看到在表射偏振约束不起作用，例如在三个摄象机沿着一直线(并非一不是普通的情况)时，此张量仍能起作用。本文中没有讨论7点的奇异组构情况(除7个共面点的明显的奇异组构外)。但再投影结果的健壮性可以指明任一这样的组构是非常少有或者不存在的。因而重要的是要研究这一课题，因为广泛相信表射偏振约束的数量上的不稳定是在于这种危形表面的存在。“三线性”张量的标记，其特性，对三视图的几何参数的关系，和对由多视图作3D重现的应用，构成了重要的未来的方向。
整个本说明中所强调的应用是通过对准作视觉识别。以新视图(ψ3)所需的最小点数(7)得到合理的性能，这可能是太多了，如果由试验点匹配的所有可能的组合来达到图象的模型匹配的话。在模型为正射的但新视图是透视的特殊情况下的双线性函数的存在从对点计数的观点上看较吸引人。这里已得到为识别透视图仅需要5个对应点(只要能满足模型是正射的要求)。没有进行以双线性函数的试验来检查在实践上需要多少点，但计划在未来进行。因为它们的简单性，人们可以推测这些代数函数将会在除视觉识别外的任务中获得应用，其中一些在下面讨论。
可能存在一些简单性成为主要重点的其他应用，而点数则较少考虑。考虑例如基于模型的压缩的应用。采用三线性函数需要17个参数来将一视图表示成完全对应中的二参考视图的函数(回想，为了将对应点数由9降到7采用了27个系数)。假定发送方和接收方都具有二参考视图并应用为获得二视图间的对应的相同算法。为发送第三视图(忽略可作分开处理的自闭问题)发送方可利用许多点对17个参数求解，但最后仅发送此17个参数。然后接收方仅只是按给于所接收参数的“三线性法”来组合此二参考视图。这显然属于点数不是要主要考虑的、而简单性和因计算的简化而得的健壮性(如上述)则是最重要的范畴。
关于图象编码，最近[1，2]推荐将图象分解为“层”的方法。此法中一系列视图被分成为区，每一区的运动被近似地以一2D仿射变换来描述。发送方发送第一个图象后仅跟随各后续帧的每一区的6个仿射参数。应用视图的代数函数有可能使这一方法更为有效，因为替代将景象分成为平面而可以试着将景象分成为对象，每一个带有描述其在随后的帧上的位移的17个参数。
应用的另外领域可以是计算机制图。再投影技术使得图象复制简化。给定某3D对象的二个经完成复制的视图，其他视图(仍然忽略自闭锁)就可仅以“组合”参考视图来加以复制。这里，对应点数仍然不作太多考虑。
参考文献[1]E.H.Adelson.Layerea representations for image coding.Technical Report 181，MediaLaboratory，Massachusetts Institute of Technology.1991.[2]E.H.Adelson and J.Y.A.Wang.Layered representation for motion analysis.In Proceed-ings IEEE Conf.on Computer Vision and Pattern Recognition，pages 361-366，NewYork.NY，June 1993.[3]G.Adiv.Innerent ambiguities in recovering 3D motion and structure from a noisyflow field.IEEE Transactions on Pattern Analysis and Machine Intelligence，PAMI-11(5)477-489，1989.[4]P.Anandan. A unified perspective on computational techniques for the measurementof visual motion.In Proceedings Image Understanding Workshop. pages 219-230，LosAngeles，CA，February 1987.Morgan Kaufmann，San Mateo，CA[5]I.A.Bachelder and S.Ullman.Contour matching using local affine transformations.InProceedings Image Understanding Workshop.Morgan Kaufmann，San Mateo，CA，1992.[6]E.B.Barrett，M.H.Brill，N.N.Haag，and P.M.Payton.Invariant linear methods inphotogrammetry and model-matching，In J.L.Mundy and A.Zisserman.editors.Ap-plications of invariances in computer vision. MIT Press，1992.[7]J.R.Bergen and R.Hingorani.Hierarchical motion-based frame rate conversion.Tech-nical report，David Sarnoff Research Center，1990.[8]S.Demey，A.Zisserman，and P.Beardsley.Affine and projective structure from motion.In Proceedings of the British Machine Vision Conference，October 1992.[9]R.Dutta and M.A.Synder.Robustness of correspondence based structure from motion.In Proceedings of the International Conference on Computer Vision，pages 106-110，Osaka，Japan，December 1990.[10]O.D.Faugeras.What can be seen in three dimensions with an uncalibrated stereo rig？In Proceedings of the European Conference on Computer Vision，pages 563-578. SantaMargherita Ligure.Italy.June 1992.[11]O.D.Faugeras and L.Robert.What can two images tell us about a third one？TechnicalReport INRIA，France.1993.[12]W.E.L.Grimson.Why stereo vision is not always about 3D reconstruction.A.I.MemoNo.1435，Artificial Intelligence Laboratory.Massachusetts Institute of Technology.July1993.[13]R.Hartley，R.Gupta，and T.Chang.Stereo from uncalibrated cameras.In ProceedingsIEEE Conf.on Computer Vision and Pattern Recognition，pages 761-764，Champaign，IL.，June 1992.[14]B.K.P.Horn.Relative orientation.International Journal of Computer Vision，459-78，1990.[15]B.K.P.Horn.Relative ocientation revisited.Journal of the Optical Society of America，
S1630-1638，1991.[16]D.P.Huttenlocher and S.Ullman.Recognizing solid objects by alignment with an image.
International Journal of Computer Vision，5(2)195-212，1990.[17]J.J.Koenderink and A.J.Van Doorn.Affine structure from motion.Journal of theOptical Society of America，8377-385，1991.[18]H.C.Longuet-Higgins.A computer algorithm for reconstructing a scene from two pro-jections.Nature，293133-135，1981.[19]H.C.Longuet-Higgins.the reconstruction of a scene from two projections-configura-tions that defeat the 8-point algorithm.In Proc.of the 1st Conf.AI applications，pages395-397，Denver，December 1984.[20]Q.T.Luong，R.Deriche，O.D.Faugeras，and T.Papadopoulo.On determining thefundamental matrixAnalysis of different methods and experimental results.TechnicalReport INRIA，France，1993.[21]Q.T.Luong and T.Vieville.Canonical representations for the geometries of multipleprojective views.Technical Report INRIA，France，fall 1993.[22]S.J.Maybank.The projective geometry of ambiguous surfaces.Proceedings of the RoyalSociety of London，3321-47，1990.[23]J.Mundy and A.Zisserman.Appendix-projective geometry for machine vision.InJ.Mundy and A.Zisserman，editors，Geometric invariances in computer vision. MITPress，Cambridge，1992.[24]J.L.Mundy，R.P.Welty，M.H.Brill，P.M.Payton，and E.B.Barrett.3D model align-ment without computing pose.In Proceedings Image Understanding Workshop.pages727-735.Morgan Kaufmann，San Mateo，CA，January 1992.[25]A.Shashua.Correspondence and affine shape from two orthographic viewsMotionand Recognition. A.I.Memo No.1327，Artificial Intelligence Laboratory，MassachusettsInstitute of Technology，December 1991.[26]A.Shashua.Geometry and Photometry in 3D visual recognition. PhD thesis.M.I.TArtificial Intelligence Laboratory，AI-TR-1401，November 1992.[27]A.Shashua.Illumination and view position in 3D visual recognition. In S.J.HansonJ.E.Moody and R.P.Lippmann，editors，Advances in Neural Information ProcessingSystems 4，pages 404-411.San Mateo，CAMorgan Kaufmann Publishers，1992.Pro-ceedings of the fourth annual conference NIPS，Dec.1991.Denver.CO[28]A.Shashua.On geometric and algebraic aspects of 3D affine and projective structuresfrom perspective 2D views. In The 2nd European Workshop on Invariants.Ponta Dela-gada. Azores，October 1993. Also in MIT AI memo No.1405，July 1993.[29]A.Shashua.Projective depthA geometric invariant for 3D reconstruction from twoperspective/orthographic views and for visual recognition. In Proceedings of the Inter-national Conference on Computer Vision，pages 583-590，Berlin，Germany，May 1993.[30]A.Shashua.Projective structure from uncalibrated imagesstructure from motion andrecognition.IEEE Transactions on Pattern Analysis and Machine Intelligence，1994.inpress.[31]A.Shashua.Trilinearity in visual recognition by alignment. In Proceedings of the Euro-pean Conference on Computer Vision，Stockholm，Sweden，May 1994.[32]A.Shashua and N.Navab.Relative affine structureTheory and application to 3dreconstruction from perspective views.In Proceedings IEEE Conf.on Computer Visionand Pattern Recognition，Seattle，Washington，1994.[33]A.Shashua and S.Toelg.The quadric reference surfaceApplications in registeringviews of complex 3d objects.In Proceedings of the European Conference on ComputerVision，Stockholm，Sweden，May 1994.[34]C.Tomasi and T.Kanade.Factoring image sequences into shape and motion.In IEEEWorkshop on Visval Motion，pages 21-29，Princeton，NJ，September 1991.[35]R.Y.Tsai and T.S.Huang.Uniqueness and estimation of three-dimensional motionparameters of rigid objects with curved surface.IEEE Transactions on Pattern Analysisand Machine Intelligence，PAMI-613-26.1984.[36]S.Ullman.The Interpretation of Visual Motion.MIT Press，Cambridge and London，1979.[37]S.Ullman.Aligning pictorial descriptionsan approach to object recognition.Cognition，32193-254，1989.Alsoin MIT AI Memo 931，Dec.1986.[38]S.Ullman and R.Basri.Recognition by linear combination of models.IEEE Transac-tions on Pattern Analysis and Machine Intelligence，PAMI-13992-1006.1991.Also inM.I.T AI Memo 1052，1989.[39]D.Weinshall.Model based invariants for 3-D vision.International Journal of ComputerVision，10(1)27-42，1993.
应当理解，这里所证明和描述的某些或所有的图象处理操作可以以软件来实现。换言之，这里所证明和描述的某些或所有的图象处理操作可以以ROM(只读存贮器)形式来实现。另外，此软件成分如果希望的话可以采用通常的技术以硬件来实现。
本发明具有非常广泛的应用，特别能应用于所有由2D到3D的技术是公认为有用的领域。本发明的应用至少包含下列这些包括由空中和卫星照相绘制地图和宇宙空间和船坞装配场地的座标测量的摄影测绘方面的应用，生产部件的座标测量，生产部件的自动光学检查，机器人单元对准，机器人轨迹识别，3D机器人反馈，景象的3D模拟，对象的3D模拟，反向工程，和3D数字化。
附录B为按照本发明一优选实施例构成和运行的3D再现设备的优选软件实现的列表。为根据附录B实现此3D再现设备，可将Math Soft提供的商业Maple2软件与IBM PCT或SUNI作站结合使用。可采用以下过程a.利用OPEN FILE命令将按附录B的列表生成的程序装载进Maple2；b.将光标移动到出现WITH(LINALG)命令的行以运行此子例程序；c.操作RETURN键直至光标到达文件的末尾；d.然后光标返回到文件的起始；e.按压RETURN直至到达随后的行来进行模拟，包括EVALN(PCOMP)然后系统在屏幕上打印出给定的3D世界与再现的3D世界间的比较。
应理解，附录B列表也有用于实现这里所证明和描述的图象变换设备和方法。此图象变换实施方案例如可基于与上述标题为“用于识别的代数函数”一节的等式3和4相结合使用的附录B的列表。
应理解附录中所描述的特定实施方案仅仅是企图用来对本发明作详尽的揭示而不是要对之作出限制。
应理解的是，为了清楚而在各个不同实施例的组成中说明的本发明的各种特点也可组合进一单个的实施例中来得到。相反，为了简要而在单个实施例的组成中说明的本发明的各种特点也可分开地或在任一适合的子组合中得到。
熟悉本技术领域的人们将会理解本发明并不限于以上已特别证明和描述的一切。而是，本发明的范围仅由后列的权利要求来确定。
BOOL ComputeAllTensors(IN int nTrials，TENSTensors[MAX_FRAMES][MAX_FRAMES][MAX_FRAMES]，OUT DMAT3 E_basis[MAX_FRAMES][MAX_FRAMES][4]，OUT int*size_E_basis){int rf，bf，i，j，k，l，m，n；int num_frames_for_i；int legalFrame[MAX_FRAMES]；double**rank4mat，*w，**v_tr，doublethresh4；BOOLfoundThresh；doublesum；int*dummy；/*for the call of findTensor*//*check the the number of frames is not too large*/if(NumberOfFrames＞MAX_FRAMES){return(FALSE)；}/*loop on the reference frame*/for(rf＝0；rf＜NumberOfFrames；rf++){/*consider only marked frame*/if(((specified_rf＜0)‖(rf＝specified_rf))&&Images_Flags[rf]){/*loop on the i′th frame*/for(i＝0；i＜NumberOfFrames；i++){/*consider only marked frame different of rf*/if(Images_Flags[i]&&(i！＝rf)){/*initialize num_frames_for_i*/(*size_E_basis)＝0；num_frames_for_i＝0；/*loop on the base frame*//**There is dependence between*Tensors[rf][i][bf]and Tensors[rf][bf][i]*so we avoid the computation of Tensors[rf][i][bf]*for bf＜i*/for(bf＝(i+1)；bf＜NumberOfFrames；bf++){/*consider only marked frame different of rf*/if(Images_Flags[bf]&&(bf！＝rf)){findTensor(Tensors[rf][i][bf]，rf，i，bf，0，dummy，0，nTrials)；legalFrame[num_frames_for_i]＝bf；num_frames for_i++；}}for(bf＝0；bf＜i；bf++){if(Images_Flags[bf]&&(bf！＝rf)){for(j＝0；j＜3；j++){for(k＝0；k＜3；k++){for(l＝0；l＜3；l++){Tensors[rf][i][bf][j][k][l]＝Tensors[rf][bf][i][k][j][l]；}}}legalFrame[num_frames_for_i]＝bf；num_frames_for_i++；}}if(num_frames_for_i＞1){/*there are num_frames_for_i Tensors for these*//*rf and i and when ordering the E0，E1，E2 to a*//*9-by-(3*num_frames_for_i)matrix its rank＝4*//*Note that E[j]＝Tensors[rf][i][bf][*][j][*]*/(*size_E_basis)＝4；/*if num_frames_for_i＞2 we consider the transposed*//*matrix for the svdcmp which requires rows＞＝cols*/if(num_frames_for_i＝2){rank4mat＝dmatrix(1，9，1，6)；w＝dvector(1，6)；v_tr＝dmatrix(1，6，1，6)；for(bf＝0；bf＜num_frames_for_i；bf++){for(j＝0；j＜3；j++){for(k＝0；k＜3；k++){for(l＝0；l＜3；l++){rank4mat[l+3*k+i][l+3*bf+j]＝Tensors[rf][i][legalFrame[bf]][k][j][l]；}}}}if(！svdcmp(rank4mat，9，6，w，v_tr)){free_dmatrix(rank4mar，1，9，1，6)；free_dvector(w，1，6)；free_dmatrix(v_tr，1，6，1，6)；return(FALSE)；}/*look for threshold for the largest 4 values*/foundThresh＝FALSE；thresh4＝0.02；while(！foundThresh){n＝0；for(k＝1；k＜＝6；k++){if(w[k]＞thresh4){n++；}}if(n＞4){thresh4*＝2；}else if (n＜4){thresh4/＝2；}else{foundThresh＝TRUE；}}/*the 4 columnsof rank4mat which correspond to*//*the 4 largest singular values form E_basis*//*Note that E_basis is[4][9]i.e.goes by rows*/n＝0；for(m＝1；m＝6；m++){if(w[m]＞thresh4){for(k＝0；k＜3；k++){for(l＝0；l＜3；l++){E_basis[rf][i][n][k][l]＝rank4mat[l+3*k+l][m]；}}n++；}}if(！ImproveE_basis(E_basis[rf][i])){free_dmatrix(rank4mat，1，9，1，6)；free_dvector(w，1，6)；free_dmatrix(v_tr，1，6，1，6)；return(FALSE)；}/*improve the tensors with the svd results*/for(bf＝0；bf＜num_frames_for_i；bf++){for(j＝0；j＜3；j++){for(k＝0；k＜3；k++){for(l＝0；l＜3；l++){sum＝0.0；n＝0；for(m＝1；m＜＝6；m++){if(w[m]＞thresh4){sum+＝E_basis[rf][i][n][k][l]*w[m]*v_tr[l+3*bf+j][m]；n++；}}Tensors[rf][i][legalFrame[bf]][k][j][l]＝sum；}}}}free_dmatrix(rank4mat，1，9，1，6)；free_dvector(w，1，6)；free_dmatrix(v_tr，1，6，1.6)；else{rank4mat＝dmatrix(1，3*num_frames_for_i，1，9)；w＝dvector(1，9)；v_tr＝dmatrix(1，9，1，9)；for(bf＝0；bf＜num_frames_for_i；bf++){for(j＝0；j＜3；j++){for(k＝0；k＜3；k++){for(l＝0；l＜3；l++){rank4mat[l+3*bf+j][l+3*k+l]＝Tensors[rf][i][legalFrame[bf]][k][j][l]；}}}}if(！svdcmp(rank4mat，3*num_frames_for_i，9，w，v_tr)){free_dmatrix(rank4mat，1，3*num_frames_for_i，1，9)；free_dvector(w，1，9)；free_dmatrix(v_tr，1，9，1，9)；return(FALSE)；}/*look for threshold for the largest 4 values*/foundThresh＝FALSE；thresh4＝0.02；while(！foundThresh){n＝0；for(k＝1；k＝9；k++){if(w[k]＞thresh4){n++；}}if(n＞4){thresh4*＝2；}else if(n＜4){thresh4/＝2；}else{foundThresh＝TRUE；}}/*the 4 rows of v which correspond to the*//*4 largest singular values form E_basis*/n＝0；for(m＝1；m＜＝9；m++){if(w[m]＞thresh4){for(k＝0；k＜3；k++){for(l＝0；l＜3；l++){E_basis[rf][i][n][k][l]＝v_tr[l+3*k+l][m]；}}n++；}}if(！ImproveE_basis(E_basis[rf][i])){free_dmatrix(rank4mat，1，3*num_frames_for_i，1，9)；free_dvector(w，1，9)；free_dmatrix(v_tr，1，9，1，9)；return(FALSE)；}/*improve the tensors with the svd results*/for(bf＝0；bf＜num_frames_for_l；bf++){for(j＝0；j＜3；j++){for(k＝0；k＜3；k++){for(l＝0；l＜3；l++){sum＝0.0；n＝0；for(m＝1；m＝9；m++){if(w[m]＞thresh4){sum+＝rank4mat[l+3*bf+j][m]*w[m]*E_basis[rf][i][n][k][l]；n++；}}Tensors[rf][i][legalFrame[bf]][k][j][l]＝sum；}}}}free_dmatrix(rank4mat，1，3*num_frames_for_i，1，9)；free_dvector(w，1，9)；free_dmatrix(v_tr，1，9，1，9)；}}else if(num_frames_for_i＝1){(*size_E_basis)＝3；/*keep the tensor as the first 3 of E_basis[rf][i]*/for(j＝0；j＜3；j++){for(k＝0；k＜3；k++){for(l＝0；l＜3；l++){E_basis[rf][i][j][k][l]＝Tensors[rf][i][legalFrame
][k][j][l]；}}}}}}}}return(TRUE)；}ComputeAllEpipoies(IN DMAT3E_basis[MAX_FRAMES][MAX_FRAMES][4]，IN int size_E_basis，IN BOOL removeOutliers，OUT DXYZ Epipoles[MAX_FRAMES][MAX_FRAMES]){int rf，i；doubleerror；int n_Eq；int row，i，j，k，m；double**Eqn；int nCols＝3，nSizes＝1，nGroups[1]，Sizes[1]；int minRows＝4，nTrials＝100，n_Eq_new；doublerow_sum；DMAT3 mat；DMAT3 mat3，Epimat；/*check the the number of frames is not too large*/if(NumberOfFrames＞MAX_FRAMES){return(FALSE)；}/*given a pair of images (rf，i)*//*there are 3*size_E_basis*(size_E_basis-1)+*//*6*size_E_basis*size_E_basis equations*/Eqn＝dmatrix(0，3*size_E_basis*(size_E_basis-1)+6*size_E_basis*size_E_basis-1，0，2)；/*loop on the reference frames*/for(rf＝0；rf＜NumberOfFrames；rf++){/*consider only marked frames*/if(((specified_rf＜0)‖(rf＝specified_rf))&& Images_Flags[rf]){/*loop on the i′th frame*/f or(i＝0；i＜NumberOfFrames；i++){/*consider only marked frames that differ from rf*/if(Images_Flags[i]&&(i！＝rf)){/*accumulate equations on Epipole*/n_Eq＝0；/*contribution of Tensors[rf][i][bf][*][j][*]*//*denoting E[j]＝Tensors[rf][i][bf][*][j][*]and*//*E[j][k]＝Tensors[rf][i][bf[*][j][k]then*/

/*the fundamental matrix.*//*In particular，transposed(E[j])*[v′]*E[i]*//*is a skew symmetric matrix and therefore*//*v′CROSS E[j][k]dot E[l][k]＝0 for j＜1*//*and also v′CROSS E[j][k]dot E[l][m]+*//*v′CROSS E[j][m] dot E[l][k]＝0 for j＜l，k＜m*//*Equivalently(being the same″volume″)*//*E[j][k]CROSS E[l][k]dot v′＝0 for j＜1 k＝0，1，2*//*and also E[j][k]CROSS E[l][m]dot v′+*//*E[j][m]CROSS E[l][k]dot v′＝0 for j＜1 k＜m*/for(row＝0；row＜(3*size_E_basis*(size_E_basis-1)+6*size_E_basis*size_E_basis)；row++){for(l＝0；l＜3；l++){Eqn[n_Eq+row][l]＝0.0；}}/*loop on E_basis[j]*/for(j＝0；j＜(size_E_basis-1)；j++){/*loop on E_basis[1]*/for(l＝(j+1)；l＜size_E_basis；l++){/*loop on the column*/for (k＝0；k＜3；k++){Eqn[n_Eq]
＝E_basis[rf][i][j][1][k]*E_basis[rf][i][l][2][k]-E_basis[rf][i][j][2][k]*E_basis[rf][i][l][1][k]；Eqn[n_Eq][1]＝E_basis[rf][i][l]
[k]*E_basis[rf][i][j][2][k]-E_basis[rf][i][l][2][k]*E_basis[rf][i][j]
[k]；Eqn[n_Eq][2]＝E_basis[rf][i][j]
[k]*E_basis[rf][i][l][1][k]-E_basis[rf][i][j][1][k]*E_basis[rf][i][l]
[k]；/*normalize the equation*/row_sum＝Eqn[n_Eq]
*Eqn[n_Eq]
+Eqn[n_Eq][1]*Eqn[n_Eq][1]+Eqn[n_Eq][2]*Eqn[n_Eq][2]；if(row_sum＞0.0){row_sum＝sqrt(row_sum)；Eqn[n_Eq]
/＝row_sum；Eqn[n_Eq][1]/＝row_sum；Eqn[n_Eq][2]/＝row_sum；n_Eq++；}}for(k＝0；k＜2；k++){for(m＝(k+1)；m＜3；m++){Eqn[n_Eq]
＝E_basis[rf][i][j][1][k]*E_basis[rf][i][l][2][m]-E_basis[rf][i][j][2][k]*E_basis[rf][i][l][1][m]+E_basis[rf][i][j][1][m]*E_basis[rf][i][l][2][k]-E_basis[rf][i][j][2][m]*E_basis[rf][i][l][1][k]；Eqn[n_Eq][1]＝E_basis[rf][i][l]
[m]*E_basis[rf[i][j][2][k]-E_basis[rf][i][l][2][m]*E_basis[rf][i][j]
[k]+E_basis[rf][i][l]
[k]*E_basis[rf][i][j][2][m]-E_basis[rf][i][l][2][k]*E_basis[rf][i][j]
[m]；Eqn[n_Eq][2]＝E_basis[rf][i][j]
[k]*E_basis[rf][i][l][1][m]-E_basis[rf][i][j][1][k]*E_basis[rf][i][l]
[m]+E_basis[rf][i][j]
[m]*E_basis[rf][i][l][1][k]-E_basis[rf][i][j][1][m]*E_basis[rf][i][l]
[k]；/*normalize the equation*/row_sum＝Eqn[n_Eq]
*Eqn[n_Eq]
+Eqn[n_Eq][1]*Eqn[n_Eq][1]+Eqn[n_Eq][2]*Eqn[n_Eq][2]；if(row_sum＞0.0){row_sum＝sqrt(row_sum)；Eqn[n_Eq]
/＝row_sum；Eqn[n_Eq][1]/＝row_sum；Eqn[n_Eq][2]/＝row_sum；n_Eq++；}}}}}/*additional size_E_basis*size_E_basis*6eqns*//*come from[v′]*E[rf][i][j]*E[i][rf][k]being*//*skew symmetric for any j，k.*/if(specified_rf＜0){/*loop on E_basis[rf][i][j]*/for(j＝0；j＜size_E_basis；j++){/*loop on E_basis[i][rf][l]*/for(l＝0；l＜size_E_basis；l++){matmulmat3(mat.E_basis[rf][i][j]，E_basis[i][rf][l])；/*ignore mat if it is a scalar matrix*/if((mat
[1]＞mm_epsilon)‖(mat
[1]＜-mm_epsilon)‖(mat
[2]＞mm_epsilon)‖(mat
[2]＜-mm_epsilon)‖(mat[1][2]＞mm_epsilon)‖(mat[1][2]＜-mm_epsilon)‖(mat[1]
＞mm_epsilon)‖(mat[1]
＜-mm_epsilon)‖(mat[2]
＞mm_epsilon)‖(mat[2]
＜-mm_epsilon)‖(mat[2][1]＞mm_epsilon)‖(mat[2][1]＜-mm_epsilon)){Eqn[n_Eq]
＝0.0；Eqn[n_Eq][1]＝mat[2]
；Eqn[n_Eq][2]＝-mat[1]
；/*normalize the equation*/row_sum＝Eqn[n_Eq]
*Eqn[n_Eq]
+Eqn[n_Eq][1]*Eqn[n_Eq][1]+Eqn[n_Eq][2]*Eqn[n_Eq][2]；if(row_sum＞0.0){row_sum＝sqrt(row_sum)；Eqn[n_Eq]
/＝row_sum；Eqn[n_Eq][1]/＝row_sum；Eqn[n_Eq][2]/＝row_sum；n_Eq++；}Eqn[n_Eq]
＝-mat[2][1]；Eqn[n_Eq][1]＝0.0；Eqn[n_Eq][2]＝mat
[1]；/*normalize the equation*/row_sum＝Eqn[n_Eq]
*Eqn[n_Eq]
+Eqn[n_Eq][1]*Eqn[n_Eq][1]+Eqn[n_Eq][2]*Eqn[n_Eq][2]；if(row_sum＞0.0){row_sum＝sqrt(row_sum)；Eqn[n_Eq]
/＝row_sum；Eqn[n_Eq][1]/＝Tow_sum；Eqn[n_Eq][2]/＝row_sum；n_Eq++；}Eqn[n_Eq]
＝mat[1][2]；Eqn[n_Eq][1]＝mat
[2]；Eqn[n_Eq][2]＝0.0；/*normalize the equation*/row_sum＝Eqn[n_Eq]
*Eqn[n_Eq]
+Eqn[n_Eq][1]*Eqn[n_Eq][1]+Eqn[n_Eq][2]*Eqn[n_Eq][2]；if(row_sum＞0.0){row_sum＝sqrt(row_sum)；Eqn[n_Eq]
/＝row_sum；Eqn[n_Eq][1]/＝row_sum；Eqn[n_Eq][2]/＝row_sum；n_Eq++；}Eqn[n_Eq]
＝-mat[2]
；Eqn[n_Eq][1]＝mat[2][1]；Eqn[n_Eq][2]＝mat
-mat[1][1]；/*normalize the equation*/row_sum＝Eqn[n_Eq]
*Eqn[n_Eq]
+Eqn[n_Eq][1]*Eqn[n_Eq][1]+Eqn[n_Eq][2]*Eqn[n_Eq][2]；if(row_sum＞0.0){row_sum＝sqrt(row_sum)；Eqn[n_Eq]
/＝row_sum；Eqn[n_Eq][1]/＝row_sum；Eqn[n_Eq][2]/＝row_sum；n_Eq++；}Eqn[n_Eq]
＝mat[1]
；Eqn[n_Eq][1]＝mat[2][2]-mat
；Eqn[n_Eq][2]＝-mat[1][2]；/*normaiize the equation*/row_sum＝Eqn[n_Eq]
*Eqn[n_Eq]
+Eqn[n_Eq][1]*Eqn[n_Eq][1]+Eqn[n_Eq][2]*Eqn[n_Eq][2]；if(row_sum＞0.0){row_sum＝sqrt[row_sum)；Eqn[n_Eq]
/＝row_sum；Eqn[n_Eq][1]/＝row_sum；Eqn[n_Eq][2]/＝row_sum；n_Eq++；}Eqn[n_Eq]
＝mat[1][1]-mat[2][2]；Eqn[n_Eq][1]＝-mat
[1]； Eqn[n_Eq][2]＝mat
[2]；/*normaiize the equation*/row_sum＝Eqn[n_Eq]
*Eqn[n_Eq]
+Eqn[n_Eq][1]*Eqn[n_Eq][1]+Eqn[n_Eq][2]*Eqn[n_Eq][2]；if(row_sum＞0.0){row_sum＝sqrt(row_sum)；Eqn[n_Eq]
/＝row_sum；Eqn[n_Eq][1]/＝row_sum；Eqn[n_Eq][2]/＝row_sum；n_Eq++；}}}}}if(removeOutliers){/*remove outliers*nGroups
＝n_Eq；Sizes
＝1；if(！removeEqOutliers(Eqn，nCols，nSizes，nGroups，Sizes，minRows，nTrials，SquaredThreshold，&n_Eq_new)){n_Eq＝n_Eq_new；}}if(n_Eq＞0){/*solve F by Jacobi*/find_sol_jacobi(Eqn.n_Eq，3，Epipoles[rf][i]，&error)；}}}}}free_dmatrix(Eqn，0.3*size_E_basis*(size_E_basis-1)+6*size_E_basis*size_E_basis-1，0，2)；return(TRUE)；}BOOL ImproveE_basis(INOUT DMAT3 E_basis[4]){inti，j，k，l，m；double**u，**v_tr，*w；double min_sv，sum；/*Note that each of the 3 4-by-3 matrices E_basis[*][*][j]is of*//*rank 2.After applying svd to each of them，ignore the smallest*//*singular value(which should be 0)in order to improve them*/u＝dmatrix(1，4，1，3)；v_tr＝dmatrix(1，3，1，3)，w＝dvector(1，3)；for(j＝0；j＜3；j++){for(k＝0；k＜4；k++){for(l＝0；l＜3；l++){u[1+k][1+l]＝E_basis[k][l][j]；}}if(！svdcmp(u，4，3，w，v_tr)){free_dmatrix(u，1，4，1，3)；free_dmatrix(v_tr，1，3，1，3)；free_dvector(w，1，3)；return(FALSE)；}/*find smallest singular value*min_sv＝w[1]；i＝1；for(l＝2；l＜＝3；l++){if(w[1]＜min_sv){min_sv＝w[1]；i＝1；}}/*replace the E_basis by the product ignoring w[i]*/for(k＝0；k＜4；k++){for(l＝0；l＜3；l++){sum＝0.0；for(m＝1；m＝3；m++){if(m！＝i){sum+＝u[1+k][m]*w[m]*v_tr[1+l][m]；}}E_basis[k][l][j]＝sum；}}}free_dmatrix(u，1，4，1，3)；free_dmatrix(v_tr，1，3，1，3)；free_dvector(w，1，3)；return(TRUE)；}FUNCTIONBOOLremoveEqOutliers /*remove outliers for homogeneous system*/(INOUT double* *Mat， /*the initial equations*/IN int nCols， /*there are nCols columns (unknowns)*/IN int nSizes， /**it is assumed that the equations are*organized in groups.These groups*have nSizes different sizes.*/IN int nGroups[]， /**there are nSizes nGroups which mark the*number of groups of each size*/IN int Sizes[]， /**there are nSizes Sizes which mark the*size of each group.*e.g.the first nGroups
are groups*each of size Sizes
，and so on*/IN int minRows， /**each iteration， try at least*minRows number of rows*/IN int nTrials， /*number of random trials to perform*/IN double SquaredThreshold，/*for boundary determination*/OUT int*outMatRows/**The number of rows in Mat to be*considered.These rows are the first*rows in the output Mat and do not*keep any grouping order*/)/DESCRIPTIONGiven the initial equations in Mat，get rid of theoutliers and keep the fitting equations as the leading rowsof the ontput Mat.Perform nTrials random trials，each time randomly choose groupsso that there are at least minRows equations.Solve this setof n_current_equations by Jacobi and check the solution againstALL equations.Keep track of the best solution(for which thenorm of the Mat*sol is minimal).Use the best solution to determine outliers(keeping theinitial partition of the equations into groups).Put the groups which fit the best solution as the leading*outMatRows rows of MatASSUMPTIONSThe equations are normalized so that there is no need toconsider the equation′s norm when performing the outliersremovalAUTHORTamir Shalom(April 10 1995)*/{int i，j，k，l；int iTrial；int iRows；int iGroup；int maxSize；int initialRows；/*how many initial rows*/int nTotalGroups； /*how many total groups*/int nAvailGroups； /*how many available groups*/double**subMat；int*from_row；/**starting row in Mat which*corresponds to Pth group*(among all nTotalGroups)*/int*to_row_plusl； /*ending row ...*/int nrandGroup；doubleerror，boundary；double*solution，*best_solution；doublesquared_norm，entry，min_error；/**Let maxSize be the maximal size of groups，thene*the maximal number of rows that will be used is*minRows+maxSize-1.*//*find maxSize，initialRows and nTotalGroups*/initialRows＝nGroups
*Sizes
；nTotalGroups＝nGroups
；maxSize＝Sizes
；for(i＝1；i＜nSizes；i++){initialRows+＝nGroups[i]*Sizes[i]；nTotalGroups+＝nGroups[i]；if(Sizes[i]＞maxSize){maxSize＝Sizes[i]；}}if(initialRows＜minRows){return(FALSE)；}/*allocate subMatrix*/subMat＝dmatrix(0，minRows+maxSize-1-1，0，nCols-1)；/*allocate solution and best_solution*/solution＝dvector(0，nCols-1)；best_solution＝dvector(0，nCols-1)；/*allocate from_row and to_row_plus1*/from_row＝ivector(0，nTotalGroups-1)；to_row_plus1＝ivector(0，nTotaiGroups-1)；/*perform the trials*/for(iTrial＝0；iTrial＜nTrials；iTrial++){/*initialize from_row and to_row_plus1*/iRows＝0；iGroup＝0；nAvailGroups＝nTotalGroups；for(i＝0；i＜nSizes；i++){for(j＝0；j＜nGroups[i]；j++){from_row[iGroup]＝iRows；iRows+＝Sizes[i]；to_row_plus1[iGroup]＝iRows；iGroup++；}}/*iRows denotes the number of already chosen equations*/iRows＝0；while(iPows＜minRows){nrandGroup＝randomn(nAvailGroups)；/*move the group to the subMat*/for(i＝from_row[nrandGroup]；i＜to_row_plus1[nrandGroup]；i++){for(j＝0；j＜nCols；j++){subMat[iRows][j]＝Mat[i][j]；}iRows++；}/*update nAvailGroups，from_row and to_row_plus1*/nAvailGroups--；for(i＝nrandGroup；i＜nAvailGroups；i++){from_row[i]＝from_row[i+1]；to_row_plus1[i]＝to_row_plus1[i+1]；}}/*solve the equations of subMat*/find_sol_jacobi(subMat.iRows，nCols.solution，&error)；/*estimate the solution against all equations*/squared_norm＝0.0；for(i＝0；i＜initialRows；i++){/**compute the Pth entry of the resuiting product of*Mat by solution*/entry＝0.0；for(j＝0；j＜nCols；j++){entry+＝Mat[i][j]*solution[j]；}squared_norm+＝entry*entry；}if((iTrial＝0)‖(squared_norm＜min_error)){min_error＝squared_norm；/*printf(″iTrial ％d error＝％g\n″，iTrial，squared_norm)；*//*keep the best solution so far*/for(j＝0；j＜nCols；j++){best_solution[j]＝solution[j]；}}}/*determine boundary for rejection of outliers*/boundary＝min_error*SquaredThreshold/(initialRows-1)；/*remove outlying groups of equations from Mat*/(*outMatRows)＝0；iRows＝0；iGroup＝0；for(i＝0；i＜nSizes；i++){for(j＝0；j＜nGroups[i]；j++){/*accumulate the norm of the whole group*/squared_norm＝0.0；for(k＝0；k＜Sizes[i]；k++){entry＝0.0；for(l＝0；l＜nCols；l++){entry+＝Mat[iRows+k][l]*best_solution[l]；}squared_norm+＝entry*entry；}WO 96/34365iGroup++；if(squared_norm＜Sizes[i]*boundary){/*add Sizes[i]equations*//*xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxprintf(″％d of kind ％d is taken\n″，j，i)；xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/if((*outMatRows)＝iRows){/*there is no need to copy the equations*//*yet(*outMatRows)should be updated*/(*outMatRows)+＝Sizes[i]；}else{for(k＝0；k＜Sizes[i]；k++){for(l＝0；l＜nCols；l++){Mat[(*outMatRows)][l]＝Mat[iRows+k][l]；}(*outMatRows)++；}}}iRows+＝Sizes[i]；}}free_dmatrix(subMat，0，minRows+maxSize-1-1，0，nCols-1)；free_dvector(solution，0，nCols-1)；free_dvector(best_solution，0，nCols-1)；free_ivector(from_row，0，nTotalGroups-1)；free_ivector(to_row_plus1，0，nTotalGroups-1)；printf(″after outliers removal only ％d out of ％d equations\n″，(*outMatRows)，initialRows)；return(TRUE)；}附件B


＞ TT[j+2，10]＝evalm(-p[1，i])＞ TT[j+2，11]＝evalm(-p(2，i])；＞ TT[j+2，12]＝-1＞ TT[j+2，16]＝evalm(ppp[1，i]*p[1，i])＞ TT[j+2，17]＝evalm[ppp[1，i]*p[2，i]) ＞ TT[j+2，18]＝evalm(ppp[1，i])＞ TT[j+2，19]＝evalm(pp[2，i]*p[1，i])＞ TT[j+2，20]＝evalm(pp[2，i]*p[2，i])＞ TT[j+2，21]＝evalm(pp[2，i])＞ TT[j+2，25]＝evalm(-ppp[1，i]*pp[2，i]*p[1，i])＞ TT[j+2，26]＝evalm(-ppp[1，i]*pp[2，i]*p[2，i])＞ TT[j+2，27]；＝evalm(-ppp[1，i]*pp[2，i])＞ TT[j+3，13]＝evalm(-p[1，i])＞ TT[j+3，14]＝evalm(-p[2，i])＞ TT[j+3，15]＝-1＞ TT[j+3，16]＝evalm(ppp[2，i]*p[1，i])＞ TT[j+3，17]＝evalm(ppp[2，i]*p[2，i])＞ TT[j+3，18]＝evalm(ppp[2，i])＞ TT[j+3，22]＝evalm(pp[2，i]*p[1，i])＞ TT[j+3，23]＝evalm(pp[2，i]*p[2，i])＞ TT[j+3，24]＝evalm(pp[2，i])＞ TT[j+3，25]＝evalm(-ppp[2，i]*pp[2，i]*p[1，i])＞ TT[j+3，26]＝evalm(-ppp[2，i]*pp[2，i]*p[2，i])＞ TT[j+3，27]＝evalm(-ppp[2，i]*pp[2，i])＞ j＝j+4＞ od＞ b＝-col(TT，27)＞ TT＝submatrix(TT.1..40，1..26)＞ sol＝leastsqrs(TT.b)＞ X＝vector(27，i-＞1)＞ for i from 1 to 26 do X[i]＝sol[i]；od＞ alpha＝array(1..3，1..3，1..3)＞ for i from 1 to 3 do＞ for j from 1 to 3 do＞ for k from 1 to 3 do＞ alpha[i，j，k]＝X[k+(i-1)*9+(j-1)*3]；＞ ododod＞ E1＝matrix(3，3)＞ E2＝matrix(3，3)＞ E3＝matrix(3，3)＞ for i from 1 to 3 do＞ for k from 1 to 3 do＞ E1[i，k]＝alpha[i，1，k]＞ E2[i，k]＝alpha[i，2，k]＞ E3[i，k]＝alpha[i，3，k]＞ odod＞＞ W1＝matrix(3，3)＞ W2＝matrix(3，3)＞ W3＝matrix(3，3)＞ for j from 1 to 3 do＞ for k from 1 to 3 do＞ W1[j，k]＝alpha[1，j，k]＞ W2[j，k]＝alpha[2，j，k]＞ W3[j，k]＝alpha[3，j，k]＞ odod＞＞ end

＞

＞ fnd_matrix＝proc(M1，M2，M3)＞ Test＝matrix(18，9，i-＞0)＞ for i from 1 to 3 do＞ Test[1，i]＝col(M1，1)[i]；Test[2，i+3]＝col(M1，2)[i]；Test[3，i+6]＝col(M1，3)[i]；＞ Test[4，i]＝col(M1，2)[i]；＞ Test[4，i+3]＝col(M1，1)[i]；＞ Test[5，i]＝col(M1，3)[i]；＞ Test[5，i+6]＝col(M1，1)[i]；＞ Test[6，i+3]＝col(M1，3)[i]；＞ Test[6，i+6]＝col(M1，2)[i]；＞ Test[7，i]＝col(M2，1)[i]；＞ Test[8，i+3]＝col(M2，2)[i]；＞ Test[9，i+6]＝col(M2，3)[i]；＞ Test[10，i]＝col(M2，2)[i]；＞ Test[10，i+3]＝col(M2，1)[i]；＞ Test[11，i]＝col(M2，3)[i]；＞ Test[11，i+6]＝col(M2，1)[i]；＞ Test[12，i+3]＝col(M2，3)[i]；＞ Test[12，i+6]＝col(M2，2)[i]；＞ Test[13，i]＝col(M3，1)[i]；＞ Test[14，i+3]＝col(M3，2)[i]；＞ Test[15，i+6]＝col(M3，3)[i]；＞ Test[16，i]＝col(M3，2)[i]；＞ Test[16，i+3]＝col(M3，1)[i]；＞ Test[17，i]＝col(M3，3)[i]；＞＞ Test[17，i+6]＝col(M3，1)[i]；＞ Test[18，i+3]＝col(M3，3)[i]；＞ Test[18，i+6]＝col(M3，2)[i]；＞ od＞ b＝-col(Test，9)＞ Test＝submatrix(Test.1..18，1..8)＞ sol＝lesstsqrs(Test，b)＞ X＝vector(9，i-＞1)＞ for i from 1 to 8 do X[i]＝sol[i]；od＞ F＝matrix(3，3)＞ for i from 1 to 3 do＞ for j from 1 to 3 do＞ F[i，j]＝X[i+(j-1)*3]；＞ odod＞ end

＞＞ skew_cross＝proc(w)＞ stack(evalm(crossprod([1，0，0]，w))；stack(evalm(crossprod(
，w))，evalm(cro＞ ssprod(
，w))))；＞ end＞

＞ get epipole from F^tv′＝0＞ vp_from_F＝proc(F)＞ G＝transpose(F)＞ b＝-col(G，3)＞ G＝submatrix(G，1..3.1..2)＞ sol＝leestsqrs(G，b)＞＞ vp_est＝vector(3，i-＞1)＞ vp_est[1]＝sol[1]vp_est[2]＝sol[2]＞ end


＞ get_Image_points＝proc()＞ randn＝rand(1..100)＞ R＝rot_mat(0.9，
)＞ Rp＝rot_mat(0.9，[2.3，0.7，1.1])＞ vpp＝vector(3，[1.1，3.3，1.0])＞ vp＝vector(3，
)＞ p＝matrix(4，10)＞ for i from 1 to 3 do＞ for j from 1 to 10 do＞ p[i，j]＝randn()/10；＞ odod＞ for j from 1 to 10 do＞ p[4，j]＝1；＞ od＞ A＝augment(R，vp)＞ B＝augment(Rp，vpp)＞ p＝evalm(submatrix(p，1..3.1..10))＞ pp＝evalm(A &*P)

＞ b＝-col(T，16)＞ T＝submatrix(T，1..15，1..15)＞ sol＝linsolve(T，b)＞ X＝vector(15，l-＞1)＞ for i from 1 to 15 do X[l]＝sol[i]；od＞ T＝matrix(4，4)＞ for i from 1 to 4 do＞ for j from 1 to 4 do＞ T[i，j]＝X[j+(i-1)*4]＞ odod＞ end＞

＞＞ getk＝proc(a，b，c)＞ evalm((transpose(crossprod(a，b))&*＞ crossprod(c，s))/(transpose(crossprod(c，a)) &*crossprod(c，a)))；＞ end＞

＞ rot_mat＝proc(angle，w)＞ w_unit＝evalm((1/norm(w，frobenius))*w)；＞ WW＝evalm(w_unit &*transpose(w_unit))＞ cos_mat＝evalm(cos(angle)*array(identtty，1..3，1.3))；＞ W_antisym＝stack(evalm(crossprod([1，0，0]，w_unlt))，stack(evalm(crossprod(
，w_unit))，evalm(crossprod(
，w_unit))))；＞ evalm(cos_mat+sin(angle)*W_antlsym+(1-cos(angle))*WW)；＞ end

＞ cross_mat＝proc(w)＞ evalm(stack(evalm(crossprod([1，0，0]，w_unit))，stack(evalm(crossprod(
，w＞ _unit))，evalm(crossprod(
，w_unit)))))；＞ end
权利要求
1.为由一3D对象的至少一2D投影生成关于此3D对象的信息的方法，此方法包括提供一3D对象的至少一2D投影；产生以αijk＝vi′bjk-vj″aik(i，j，k＝1，2，3)，描述的数的阵列，其中aij和bjk分别为矩阵A和B的阵元及vi′和vi″分别为向量v′和v″的元素，而所述矩阵和向量共同说明此3D对象的三个视图的摄象机参数；和利用此阵列生成关于此3D对象的信息。
2.为由一3D对象的具有n个对应点pi和pi′(i＝1……n)的第一和第二现有的视图产生此3D对象的一新视图的方法，此方法包括产生一张量αijk；及将点pi和pi′(i＝1……n)的值插入到三线性公式，由此来提取表示此新视图中一位置的x″，y″值；和根据插入步骤的结果生成新的视图。
3.为由一3D对象的至少一2D投影再现此3D对象的方法，此方法包括提供一3D对象的至少一2D投影；产生一以αijk＝ vi′bjk-vj″aik(i，j，k＝1，2，3)，描述的数的阵列，其中aij和bjk分别为矩阵A和B的阵元及vi′和vi″分别为向量v′和v″的元素，而所述矩阵和向量共同说明此3D对象的三个视图的摄象机参数；将此阵列置换成与3D空间中的三个对应平面相关的三个单应性矩阵；和利用此三个单应性矩阵再现此3D对象。
4.一光学识别方法，包括提供一3D对象的相互间存在一三线性关系的三个透视图；和利用视图之间的三线性关系以便由对准来进行光学识别。
5.按照权利要求4的方法，还包括再投影此3D对象。
6.按照权利要求1的方法，其中所述关于3D对象的信息包括3D对象的再现。
7.按照权利要求1的方法，其中所述关于3D对象的信息包括至少一个未作再现3D对象而产生的3D对象的新视图。
8.按照权利要求1的方法，其中所述至少一2D投影包括至少一航空照片。
9.按照权利要求1的方法，其中所述至少一2D投影包括至少一卫星照片。
10.按照权利要求1的方法，其中所述关于3D对象的信息包括至少一3D对象的座标。
11.按照权利要求1的方法，其中此3D对象包括宇宙空间对象。
12.按照权利要求1的方法，其中此3D对象包括很大的对象如一船舶。
13.按照权利要求1的方法，其中此3D对象包括不实在的对象。
14.为由一3D对象的至少一2D摄影产生关于此3D对象的信息的设备，该设备包括为提供一3D对象的至少一2D投影的设备；运行产生以αijk＝vi′bjk-vj″aik(i，j，k＝1，2，3)，描述的数的阵列的阵列发生器，其中aij和bjk分别为矩阵A和B的阵元及vi′和vi″分别为向量v′和v″的元素，而所述矩阵和向量共同说明此3D对象的三个视图的摄象机参数；和利用此阵列产生关于此3D对象的信息的阵列分析器。
15.为由一3D对象的具有n个对应点pi和pi′(i＝1……n)的第一和第二现有的视图产生此3D对象一新视图的设备，该设备包括为产生一张量αijk的设备；和为将点pi和pi′(i＝1……n)的值插入成为三线性形式的设备，由此提取表示新视图中一位置的x″，y″值；和为根据插入步骤的结果产生新视图的设备。
16.为由一3D对象的至少一2D投影再现此3D对象的设备，该设备包括为提供一3D对象的至少一2D投影的设备；运行来产生一以αijk＝vi′bjk-vj″aik(i，j，k＝1，2，3)，描述的数的阵列的阵列发生器，其中aij和bjk分别为矩阵A和B的阵元及vi′和vi″分别为向量v和v″的元素，而所述矩阵和向量共同说明3D对象的三个视图的摄象机参数；为将此阵列置换成为与3D空间中的三个对应平面相关的三个单应性矩阵运行的阵列置换器；和利用此三个单应性矩阵再现此3D对象运行的3D对象再现器。
17.视觉识别设备，包括为提供一3D对象的相互间存在有三线性关系的三个透视图的设备；和利用此视图间三线性关系以便由对准来进行视觉识别的设备。
18.按照权利要求1的方法，其中所述至少一2D投影包括三个2D投影。
19.按照权利要求14的设备，其中所述至少一2D投影包括三个2D投影。
20.按照权利要求1的方法，还包括利用所述方法的结果来进行进行生产部件的座标测量。
21.按照权利要求1的方法，还包括利用所述方法的结果来进行生产部件的自动光学检验。
22.按照权利要求1的方法，还包括利用所述方法的结果来进行机器人单元校准。
23.按照权利要求1的方法，其特征是还包括利用所述方法的结果来进行机器人轨迹识别。
24.按照权利要求1的方法，还包括利用所述方法的结果来进行3D机器人反馈。
25.按照权利要求1的方法，还包括利用所述方法的结果来进行景象的3D模拟。
26.按照权利要求1的方法，还包括利用所述方法的结果来进行对象的3D模拟。
27.按照权利要求1的方法，还包括利用所述方法的结果来进行反向工程。
28.按照权利要求1的方法，还包括利用所述方法的结果来进行3D数字化。
29.3D景象再现方法，用于由一3D景象的第一、第二和第三视图产生此3D景象的3D表示，此方法包括提供一3D景象的第一、第二和第三视图；利用关于所述第一、第二和第三视图的几何信息产生一表示第一、第二和第三视图之间的几何关系的三线性张量；和由所述三线性张量产生3D景象的3D表示。
30.按照权利要求29，其中所述产生3D表示的步骤包括按此三线线张量计算第一和第二视图的表射偏振几何表示；和由所述表射偏振几何表示生成所述3D表示。
31.图象变换设备，为由一3D景象的第一和第二参考视图产生此3D景象的新视图，此设备包括为提供一3D景象的第一和第二参考视图的设备；一三线性张量发生器，运行来分别利用关于所述第一参考视图、第二参考视图和新视图的几何信息产生表示第一、第二和新视图之间的几何关系的三线性张量；和一新视图发生器，运行来以根据所述第一和第二参考视图中的第一和第二对应位置和所述三线性张量计算多个各自分别对应于不同的所述第一和第二对应位置的新视图位置产生所述新视图。
32.3D再现设备，为由一3D景象的第一、第二和第三视图产生此3D景象的一3D表示，该设备包括为提供3D景象的第一、第二和第三视图的设备；一三线性张量发生器，运行来利用关于所述第一、第二和第三视图之间的几何关系来产生表示所述第一、第二和第三视图之间的几何关系的三线性张量；和3D景象表示发生器，运行来由所述三线性张量产生3D景象的3D表示。
33.按照权利要求29的方法，所述提供步骤包括提供3D景象的m＞3的视图，且所述利用步骤对m个视图中的多个3视图的子集的每一个进行，由此来产生多个中间张量，此方法还包括在所述产生步骤之前将此多个中间张量组合成最终的三线性张量，所述产生的步骤包括由所述最终的三线性张量生成3D景象的3D表示。
34.按照权利要求33的方法，其中所述利用步骤对第一和第二视图结合每一其余的视图进行，由此产生m-2个中间张量。
35.按照权利要求33的方法，其中所述组合步骤包括将每一中间张量排列在一对应的9×3(m-2))矩阵中；和将此最终的三线性张量定义作为所述矩阵的四个最大的基本成份。
36.匹配点求取设备，为寻求三个视图之间的匹配位置，此设备包括提供三视图之间的一初始匹配位置组；根据所述初始组产生表示所述三视图间的关系的三线性张量；和利用三线性张量产生一最终的匹配位置组。
37.按照权利要求36的设备，其中所述利用步骤包括对三视图的第一视图中多个位置的每一个由张量产生各自与其余二个视图相关的第一第二对应的表射偏振线；沿第一表射偏振线选择一第一候选匹配位置，并根据第一匹配位置计算沿第二表射偏振线的第二候选匹配位置；和重复选择步骤直至此二候选匹配位置与第一视图位置相符。
38.按照权利要求29的方法，还包括对3D表示进行表面插补，由此来产生作表面插补的输出；和由作表面插补的输出产生正射投影照片。
39.按照权利要求29的方法，其中所述提供的步骤包括顺序地在景象内的一序列位置上定位至少三个成象装置，在此序列内的每一位置利用每一个成象装置形成此景象的至少一部分的图象。
40.按照权利要求29的方法，还包括将3D景象的3D表示与所存贮的此3D景象的模型进行比较。
全文摘要
用于由一3D对象的至少一2D投影产生关于此对象的信息的方法。此方法包括提供一3D对象的至少一2D投影(40),产生一以aijk=vi′bik－vj″ajk(i,j,k=1,2,3)描述的数阵列(50、60),其中ajk和bjk分别为矩阵A和B的元而vi′和vj′分别为向量v′和v″的元素,而矩阵(50)和向量60一起表述此3D对象的三个视图(102)的摄象机参数,并利用此矩阵来产生有关此3D对象(70)的信息。
文档编号G01B11/00GK1198230SQ96194705
公开日1998年11月4日 申请日期1996年4月24日 优先权日1995年4月25日
发明者阿姆农·沙舒阿 申请人:考格内特恩斯有限公司