本发明涉及古生物化石技术领域，更为具体的是涉及一种古生物化石的三维重建方法。
背景技术：
化石是古生物学研究的基本材料，其主要信息都携带在化石的形状之中，化石标本是一类特殊的物体，是生物在死亡以后经过搬运或否，然后埋葬，在经石化过程而成。在埋藏和石化过程中，软组织腐烂降解，硬组织矿化成岩石，由于生物的结构很复杂，直到成为化石以后有的形态暴露在外面而一目了然，而有些结构形态被包藏在骨骼内部而无法直接看到。如果是数量较多的标本，可通过实体切片进行观察，但当标本数量很少，属于珍稀标本时，只能进行无创解剖，即通过透视方法获取化石标本内部的图像，并根据这些图像提供的三维信息来研究化石内部的构造形态。
目前，获取标本内部形态信息的途径主要通过CT扫描成像的方法，而CT扫描成像由一系列二维X光透视扫描图像组成，可以通过计算机程序进行三维图像的重建，在进行影像分割提取工作时，需要将每一张二维切片影像进行观察，鉴别不同的器官、组织或结构成分的界线。
但化石保存状况受诸多因素影响，由于化石标本在保存中经过了埋藏或者搬运，在软组织腐烂讲解以后，有些空穴会被围岩或矿物质填充，使得标本内部的密度分布发生与原结构不同的现象。有时化石标本内部或有变形，或因充填物的密度与骨骼接近而无法区分骨骼与充填物的界线，使得标本的影像色度分布复杂，难以使用现有的自动选择工具进行图形分割，只能使用手工分割，不仅使得提取的三维图形进度下降，并且分割的工作量很大。
技术实现要素：
为了解决上述问题，本发明提供了一种计算过程简单并能快速对化石内部复杂结构进行重建，且能准确区分骨骼和充填物界线的古生物化石的三维重建方法。
本发明采用的技术方案是：一种古生物化石的三维重建方法，包括以下步骤：
S1，向古生物化石投射线状激光进行CT扫描；
从至少两个角度连续采集射线状激光照射的古生物化石的图像信息，得二维CT扫描切片影像；
S2，将二维CT扫描切片影像转换为BMP格式；
S3，对BMP格式的二维CT扫描切片影像进行预处理；
所述预处理包括图像平滑、图像增强、图像插值、图像分割过程；
S4，对经预处理后的二维CT扫描切片影像进行压缩处理得压缩CT扫描切片数据；
S5，对压缩CT扫描切片数据进行图像处理；
S6，将经图像处理后的数据进行非线性优化与几何拟合得拟合数据；
S7，对拟合数据进行数据可视化计算，得3D可视化模型；
S8，对3D可视化模型进行动态模拟计算。
作为上述方案的进一步设置，所述图像平滑采用中值滤波法。
作为上述方案的进一步设置，所述图像增强采用拉普拉斯锐化法。
作为上述方案的进一步设置，所述图像插值采用线性加权平均法。
作为上述方案的进一步设置，所述图像分割包括图像二值化、边缘检测两个步骤。
作为上述方案的进一步设置，所述边缘检测采用Canny边缘算子。
作为上述方案的进一步设置，所述图像处理包括2D和3D图像处理。
作为上述方案的进一步设置，所述可视化模型进行动态模拟计算前先进行图像修复。
本发明计算过程简单、运算速度高，可对含有充填物的生物化石影像的边界进行识别与自动化分割。
附图说明
图1为本发明的结构示意图。
具体实施方式
下面结合附图及实施例对本发明做进一步描述。
如图1所示，一种古生物化石的三维重建方法，包括以下步骤：
S1，向古生物化石投射线状激光进行CT扫描；
从至少两个角度连续采集射线状激光照射的古生物化石的图像信息，得二维CT扫描切片影像；
如果被测物体较大，可采取分块多测站的形式对古生物化石进行扫描，然后进行量测或对应对点云进行拼接，为保证拼接的准确性，应使相邻测站扫描的影像具有一定的重叠度。
S2，将二维CT扫描切片影像转换为BMP格式；
S3，对BMP格式的二维CT扫描切片影像进行预处理；
所述预处理包括图像平滑、图像增强、图像插值、图像分割过程；
由于激光扫描所得到的数据是全离散的矢量距离点，也即是点云，在扫描的过程之中，扫描速度、设备精度、被测古生物化石的表面情况和操作者的熟练程度等都会对测量数据造成影响，使得得到的数据可能带有很多离散点和小振幅噪声，会影响重建后模型的质量。
因此，用计算机对图像做后处理，对获取的图像进行增强信噪比的工作，即滤除图形的噪声和干扰，突出感兴趣对象区域或边缘，便于边缘提取和图像分割。
1、图像平滑
对于噪声干扰的图像，可用图像平滑方法来滤除噪声。在频率域中，低频相应于空间的大物体，高频相应于细节和边缘，与噪声相混淆。一个较好的平滑方法或滤波方法应该是既能消掉这些噪声又不使图像的边缘轮廓和线条变模糊。
中值滤波在一定条件下可以克服线性滤波器如最小均方滤波、均值滤波等带来的图像细节模糊，而且对抑止图像中的脉冲干扰和椒盐噪声特别有效。因为这些噪声在图像中往往以孤立点的形式出现，与之对应的像素又比较少，所以采用中值滤波能有效地去除这些噪声，达到图像平滑的目的。
中值滤波法一般采用一个含有奇数个点的滑动窗口，对该滑动窗口内的诸像素灰度排列，用其中值代替窗口中心像素原来的灰度值，或者说用局部中值代替局部平均值。具体来说，就是假设有一个N维奇数的离散序列a1，a2，…，an，其中值为m，那么用它来代替窗口中心所对应像素的灰度。例如，取一维序列的元素为N＝5个，数值分别为80,90,200,110,1200计算的中值m＝110。将m＝110放在这5个数中间，对序列重新排序后得80,90,110,120,200。显然，如果200是一个噪声点，通过中值滤波后，图像质量就会有很大的改进。
一般情况下，序列为N的一维窗口的中值为：
median(x1，x2，...，xN)＝m (I)
如果N为偶数，则取两个中间数的平均值作为中值。
2、图像增强
图像增强的目的是为了改善图像的视觉效果，或者是为了更便于人或计算机的分析和处理。为了使一幅图像的边缘更加鲜明，有时也需要对图像进行尖锐化增强处理，因而图像尖锐化技术也常常用做图像增强的工具。在图像增强的过程中突出了一部分信息，同时可能会压制另一部分信息。
本发明的图像增强采用拉普拉斯锐化法，原理如下：
对于给定的二维离散图像∫(x，y)，可以计算出其一阶差分为：
二阶差分为：
从而根据拉普拉斯算子：
计算出：
对于因扩散引起的图像模糊，可以用式(3.12)进行锐化：
其中，k是与扩散有关的函数。如果令k＝1，则有：
g(x，y)＝5∫(x，y)-∫(x-1，y)-∫(x+1，y)-∫(x，y+1)-∫(x，y-1) (7)
3、图像插值
当断层图像间的距离比断层图像内像素间的距离大得多时，就需要用图像插值方法在原来的断层图像之间再插值生成一些中间断层图像。因为许多体数据三维重建和三维显示方法要求体数据是各向同性的。因而断层图像就要插值成各向同性，即经插值后的断层图像序列中断层间距等于断层图像内像素间距。然而，图像插值是一个具有很大任意性的问题。为了使图像插值成为一个确定的、可解的问题，通常引入下面三个约束条件：
(1)插值图像要与原始断层图像相似；
(2)插值图像与两个原始断层图像的相似度应该分别和它与这两个断层图像的距离成反比关系；
(3)插值图像序列应该呈现出从一幅原始断层图像到另一幅原始断层图像的渐变过程。
最简单的插值方法是对上下两个相邻的断层图像进行加权平均，产生一组插值图像。当断层间距与断层图像内像素间距相差不是很大时，采用这种插值方法是可行的。
但是，当断层间距与像素间距相差很大时，这种方式生成的插值图像模糊不清。其原因是，两个相邻断层图像中处于同一位置的像素不一定对应同一物质，他们的加权平均没有什么意义。对断层间距较大的插值，一种较好的方法是基于匹配的图像插值，这种方法能够更好的构造中间断层图像。关于断层图像的匹配是模式识别和计算机视觉领域所研究的主要问题之一。图像匹配的任务就是要在两幅图像之间建立对应关系，也就是寻找图像间的变换。传统的匹配方法一般包括特征选取、特征匹配和整幅图像之间匹配映射关系的确定三大步骤，匹配的可靠性和准确性很大程度取决于是否能够找到鲁棒的方法去完成特征提取和特征匹配。
本发明采用线性加权平均法来进行图像差值，线性加权平均的图像插值方法描述如下：
设Sk(i，j)、Sk+1(i，j)分别是第k层和第k+1层切片图像。他们之间的插值图像可表示为：
Sλ(i，j)＝(1-λ)·Sk(i，j)+λ·Sk+1(i，j) (8)
其中λ＝d1/(d1+d2)，d1、d2分别是插值图像到第k、k+1层图像的距离。显然，当λ＝0时，Sλ(i，j)＝Sk(i，j)，当λ＝1时，Sλ(i，j)＝Sk+1(i，j)。给出一组{λi|λi∈(0，1)，i＝1，2，...，n}，就相应的得n个插值图像。为得到等间距的插值图像，参数序列{λi|i＝1，2，...，n}应该取
4、图像分割
对CT图像进行基于边缘检测的图像分割处理，包括两个步骤：图像二值化和边缘检测。
(1)图像二值化
在进行图像处理时，由于普通计算机的内存和计算能力都是有限的，从而为了有效降低需要处理的数据量，同时清晰分辨图像边界轮廓，可采用图像二值化的方法对图像样本进行预处理。二值图像有以下几方面特点：
①计算二值图像特性的算法非常简单，容易理解和实现，并且计算速度很快。
②二值图像所需的内存少，对计算设备要求低。
③许多二值图像处理技术也可以用于灰度图像的处理。
二值化指将灰度图像转换为二值图像，二值图像指图像中的每个象素只取两个离散的值之一，即非此即他。用式(9)表示：
上式中，∫(x，y)表示一幅数字图像，x、y表示该图像中某象素的坐标值，而0和1表示该象素的象素取值。
由灰度级直方图(Grey Level Histogram)确定整体阈值
设规范化灰度值g范围为0≤g≤1，g＝0为最黑，g＝1为最白。M为灰度级数目，p(gk)为第k级灰度的概率。nk是在图像中出现的灰度级为k的次数，n为图像中像素的总数。则有：
p(gk)＝nk/n 0≤g≤1，k＝1，2，...，M (10)
通常称以p(gk)为纵坐标，以gk为横坐标的图形为灰度直方图。阈值应该取在两个峰值的波谷处，波谷越深陡二值化效果越好。
根据灰度直方图选取阈值T将图像分成C1和C2两个类，其中C1类像素灰度值比T值高，C2类像素灰度值比T值低，计算类间方差σ2B和类内方差σ2W：
类间方差：
类内方差：
其中，w1和w2分别是类C1和类C2的像素数，u1和u2分别是类C1和类C2中像素的灰度平均值，u是所有像素的灰度平均值，σ21和σ22分别是类C1和类C2的灰度方差。
使分离η(T)为最大值的阈值T即为最佳阈值，即：
(2)边缘提取
将CT图像二值化后，根据图像本身的情况可以对其进行边缘检测。其目的是把图像中人们感兴趣的部分分离出来，突出想要的目标，减少信息量。
图像的边缘是图像的最基本特征，所谓边缘是指其周围像素灰度有阶跃变化或屋顶变化的那些像素的集合，本发明采用Canny边缘算子来计算：
Canny检测主要检测阶跃性边缘，Canny算子的原理如下，
式(14)定义了G(X)：
分别对式(14)求x偏导，求y偏导得：
利用两个微分模板估计梯度矢量
其幅度边缘为强度梯度方向角为
对像素边缘强度和梯度方向角应用非最大值抑制技术，具体如下：
①计算梯度方向上与3×3边界相交的两个虚拟像素Q1和Q2。(i，j)是3×3模块的中心像素，Q1和Q2是(i，j)梯度方向上与3×3边界相交的虚拟像素。
②用所Q1在边框上的两个与Q1相邻像素的梯度强度线性内插Q1点的边缘强度用Q2所在边框上的两个与Q2相邻像素的梯度强度线性内插Q2点的边缘强度
③如果中心像素(i，j)的边缘强度同时大于和则保留此中心像素为候选边缘，否则极为非边缘像素；
应用Hysteris阈值判断是否是边缘像素，具体如下：
①对每个候选边缘，如果则标记为边缘像素；
②对剩下的每个候选边缘，如果在其3×3邻域，至少已有一个邻像素是边缘，那么标记为边缘；
③迭代进行第二步，知道没有新的边缘像素被标记出位置，则剩下的像素被标记为非边缘像素；经过以上步骤，图像就变成二值化的边缘图像。
S4，对经预处理后的二维CT扫描切片影像进行压缩处理得压缩CT扫描切片数据；
由于测量精度的提高，形成的数据可能会很大，而且其中会包括大量的冗余数据，必须对数据进行压缩处理，因此，采取先确定整个数据在整个三维空间的最外区域，然后，将整个区域细分，在每一个小区域内求小区域内的点云，然后根据点云曲率的特点压缩点云，主要步骤如下：
①确定原始点云数据最外区域；
②区域细化；
③细化区域内点云曲率的计算；
④计算小区域内每两个的点云曲率的差值；
⑤确定曲率差值的阈值，保留小阈值的点，删除大于阈值的亮点中曲率较小的点；
⑥压缩结果不满足要求时，循环步骤4和5，第二次压缩至数据压缩比在70％～80％。
S5，对压缩CT扫描切片数据进行图像处理；
由于投影数据形成的图像为2D图像，不能直观的显示古生物化石的效果，因此对滤波反投影数据进行2D和3D图像处理，使之可以形成3D图像
S6，将经图像处理后的数据进行非线性优化与几何拟合得拟合数据；
S7，对拟合数据进行数据可视化计算，得3D可视化模型；
由于发掘的古生物化石在形成过程中可能出现破碎、扭曲和不完整的情况，因此需要对可视化模型进行图像修复，图像修复采用Criminisi算法，其算法如下：
(1)读入可视化模型的图像与图像掩码；
(2)计算待修复区域的边界，初始化置信度和数据项；
(3)计算边界上像素点的置信度，数据项，以及邻域相关信息，根据优先权计算公式计算各点的优先权；
(4)以优先权最高的像素点为中心，形成初始待修复模板，再根据模板自适应大小方法，修改模板大小，若模板超过一定阈值大小，记录其边界；
(5)按照新的模板匹配方式，计算每个已知模板与待修复块的直方图相交距离和两者间的SSD值，以搜寻最佳匹配模板；
(6)将最佳匹配模板的信息拷贝到待修复模块中的未知部分，更新置信度和优先权等信息；
(7)若边界不为空，转到第三步继续执行；否则，转到下一步；
(8)根据第三步记录的较大模块的边界，对其采用Telea算法再进行一次快速修复，以消除边界效应，然后算法结束。
S8，对3D可视化模型进行动态模拟计算，使形成的古生物化石可视化模型可以以动态形式呈现。
上面结合附图对本发明的具体实施方式作了详细说明，但是本发明并不限于上述实施方式，在本领域普通技术人员所具备的知识范围内，还可以在不脱离本发明宗旨的前提下作出各种变化。