基于MRI的半月板三维重建系统的制作方法
本发明涉及医学影像技术领域，具体涉及一种基于mri的半月板三维重建系统。
背景技术：
半月板是膝关节保持稳定、进行正常功能的重要生理结构，半月板损伤是引起膝关节疼痛、功能障碍最常见的原因之一。核磁共振成像(mri)是半月板损伤最主要的无创性诊断方法，医生大多凭借经验由多幅二维图像去估计损伤的部位和形状。但二维断层图像表达的是某一个截面的解剖信息，而半月板是三维立体结构，损伤部位上有体部、边缘及前角、后角之分，损伤形态上有水平撕裂型、纵行撕裂型、放射状撕裂型、斜行撕裂型、复杂撕裂型之分。若单纯根据二维平面半月板内损伤的异常信号来分型，容易遗漏重要的临床信息。因此，如果能通过计算机三维重建技术显示半月板损伤的三维形态结构，并由计算机自动给出诊断，必然有利于降低人为因素导致的误诊率和漏诊率，确立更有效的治疗方案，提高医生的工作效率。
但半月板三维重建与损伤识别的前提，是将膝关节mri中的半月板从周围图像背景中提取出来，也即所谓的图像分割(segmentation)。如果没有经过半月板的图像分割，会导致得到的三维重建图像混合了膝关节韧带、软骨、肌肉等多种组织，严重影响计算机的识别与医生的观测。由于医学图像中的噪声和不均匀问题，膝关节半月板的识别特征不明显，与周围组织在空间中交错复杂，且图像灰度分布互相重叠，导致用经典的图像处理技术难以进行快速有效的半月板图像分割。
由于传统mri诊断半月板损伤的局限性，有学者尝试利用不同序列或成像技术加强其诊断能力。吴婧等通过对比mri质子压脂成像(fatsuppressed-protondensity-weightedimaging,fs-pdwi)和关节镜对膝关节半月板损伤的检查结果，认为fs-pdwi灵敏度和特异度均较高。刘新球等研究发现mri扫描参数的优化组合和选择正确的扫描方法、恰当的线圈及应用脂肪抑制技术更有助于显示正常膝关节解剖结构及病变的定位和诊断。nguyenjc等认为快速自旋回波(fse)3d序列能够提供高分辨率的半月板图像，降低临床的误诊率。赖清泉等则研究认为mri平扫与mri造影对半月板损伤的诊断准确性差异无统计学意义。choijy等通过对比研究发现，mri的utet2序列较普通t2序列可以更好地区分半月板的层状层结构，有利于发现早期的半月板损伤改变。calixtone等研究认为，高分辨率fset1ρ及t2序列能更有效地分辨膝关节骨性关节炎患者的早期半月板退变性损伤。周进等认为常规的mri序列首先由于其间隔扫描会导致病灶的漏诊，其次由于扫描角度的原因常造成膝关节各结构显示的不连续性，不利于诊断，第三由于常规序列扫描时间长，常导致急性外伤疼痛患者无法耐受而导致扫描失败，因此采用3d-cubepd序列进行扫描，但经比较发现与常规的2s-fsepdwi序列在显示半月板方面无统计学差异；而赵晓梅等研究也表明2d和3dmr序列在诊断内侧半月板和外侧半月板损伤的敏感度、特异度和准确度上并无显著差异。因此，尽管mri技术的不断提高在理论上有助于清晰准确的诊断半月板损伤，但单纯的二维图像信息始终无法直观地体现半月板损伤的部位及损伤形态。
目前较为常用的通用型影像学处理软件，如mimics等，存在一定的不足：第一，其处理过程有一定复杂性，必须由熟练掌握相关软件的人员进行操作；第二，由于其并非专用于膝关节mri三维重建，需要进行较多的参数调校，不利于临床的快速操作；第三，由于软件是由厂家开发完成，程序接口的开源性不强，不利于进一步的自主研究，如生物力学分析等。因此，有必要从临床实用性出发，为半月板损伤个性化定制一款三维诊断软件系统。
但半月板的图像分割是进行三维重建的基础，由于膝关节mri图像中各种解剖结构大小不一，在空间中有机交错，形态复杂，很难进行计算机自动分割。
人工智能(artificialintelligence)是研究计算机模拟、延伸和扩展人的智能的一门新的技术科学，近年来发展迅速。其中的模糊集合、蚁群算法等人工智能图像识别算法与人工神经网络识别技术，能较好地解决图像中的噪声和不均匀问题，快速准确地实现医学图像中感兴趣区域的图像分割。比较著名的应用有ibm的诊断系统，能快速有效地诊断肺癌、乳腺癌、结肠癌、直肠癌等肿瘤，并被浙江省中医院引入使用；阿里健康研发的医疗ai“doctoryou”，能在30分钟内阅完九千张ct影像，对肺部结节的诊断准确度达到90％以上。
随着近年来人工智能的发展，使计算机自动进行半月板的图像分割成为可能；而且经典图像处理技术与人工智能多种类型方法的结合与改进使用，相互弥补算法功能的缺陷，已是医学图像处理一个重要的发展方向。何璇研究并提出了一种基于多特征支持向量机和弹性自动区域生长的算法对膝关节mri图像中的软骨进行自动分割，并对分割出的软骨进行三维重建及定量分析，对基于人工智能技术的半月板分割算法具有一定的参考意义。zarandimh等利用2型模糊图像处理模型及神经感知网络技术，通过预处理、分割和分类的步骤来识别半月板损伤。因此，基于人工智能算法的半月板损伤三维诊断系统是可行的，而且也必将是未来用计算机自动诊断半月板损伤的研究趋势。
因此，需要提供一种基于mri的半月板三维重建系统，在人工智能算法的辅助下，快速准确地提取膝关节mri中的半月板数据，生成三维重建图像；并能自动识别半月板损伤的部位及形态，作出初步诊断，提高医生工作效率。
技术实现要素：
本发明目的是提供一种基于mri的半月板三维重建系统，其能够快速准确地提取膝关节mri中的半月板数据，生成交互式的三维重建图像；自动识别半月板损伤的部位及形态，作出初步诊断。
本发明通过以下技术方案来实现上述目的：
一种基于mri的半月板三维重建系统，其特征在于：该系统包括依次连接的原始图像获取模块、图像处理模块、3d图像输出模块、图像存储模块和损伤诊断模块。
进一步地，图像处理模块进一步包括区域提取组件、边缘提取组件、中心计算组件、阈值调整组件和卷积神经网络组件。
进一步地，原始图像获取模块在接收到序列切片图像后，将上述切片图像进行格式转换、灰度转换和去噪滤波。
进一步地，去噪滤波通过高斯滤波的方式进行。
进一步地，区域提取组件，用于在接收到序列切片图像后，对每一张切片图像进行图像序列灰度信息分析，并获得可能包含半月板目标的最小计算区域。
进一步地，区域提取组件用于计算得到所有切片图像的最小计算区域后，接收上述序列切片图像的最小计算区域的图像信息，并进一步确定半月板图像更准确的区域边缘信息。
进一步地，中心计算组件，用于接收边缘提取组件提供的图像信息和边缘图像函数信息，并根据上述信息计算每一张序列切片图像中的边缘函数图像围城的区域的中心oi(k,l)。
进一步地，阈值调整组件，用于接收中心计算组件传送的图像信息以及图像中心oi(k,l)的数据信息，并根据oi(k,l)数据信息制作切片时序-中心点x坐标和切片时序-中心点y坐标函数图；通过人工智能数据分析，筛选上述切片时序-中心点x坐标和切片时序-中心点y坐标函数图中的函数曲线突变点。
进一步地，如果不存在突变点，则不需要反馈控制边缘提取组件调整高低阈值；如果存在突变点，则需要反馈边缘提取组件需要调整相应图像的高低阈值。
一种所述的一种基于mri的半月板三维重建系统的一种基于mri的半月板三维重建方法，该方法包括：
s1.接收到上述序列切片图像后，将上述切片图像进行格式转换、灰度转换和去噪滤波；
s2.进行图像序列的图像处理：
s21.对每一张切片图像进行图像序列灰度信息分析，并获得可能包含半月板目标的最小计算区域；
s22.在区域提取组件计算得到所有切片图像的最小计算区域后，接收上述序列切片图像的最小计算区域的图像信息，并进一步确定半月板图像更准确的区域边缘信息；
s23.接收边缘提取组件提供的图像信息和边缘图像函数信息，并根据上述信息计算每一张序列切片图像中的边缘函数图像围城的区域的中心oi(k,l)；
s24.接受序列扫描图像的边界以及修正边界，并通过卷积神经网络进行三维图像的重建。
本发明的有益效果在于：快速准确地提取膝关节mri中的半月板数据，生成三维重建图像；并能自动识别半月板损伤的部位及形态，作出初步诊断，提高医生工作效率，降低误诊率和漏诊率。
附图说明
图1是本发明提供的基于mri的半月板三维重建系统的结构示意图；
具体实施方式
以下将对本发明的基于mri的半月板三维重建系统作进一步的详细描述。
下面将参照附图对本发明进行更详细的描述，其中表示了本发明的优选实施例，应该理解本领域技术人员可以修改在此描述的本发明而仍然实现本发明的有益效果。因此，下列描述应当被理解为对于本领域技术人员的广泛知道，而并不作为对本发明的限制。
为了表述清楚，不描述实际实施例的全部特征。在下列描述中，不详细描述公知的功能和结构，因为它们会使本发明由于不必要的细节而混乱。应当认为在任何实际实施例的开发中，必须做出大量实施细节以实现开发者的特定目标。
为使本发明的目的、特征更明显易懂，下面结合附图对本发明的具体实施方式作进一步的说明。需要说明的是，附图均采用非常简化的形式且均使用非精准的比例，仅用于方便、清晰地辅助说明本发明实施例的目的。
实施例1基于mri的半月板三维重建系统
本发明提供的基于mri的半月板三维重建系统包括：依次连接的原始图像获取模块、图像处理模块、3d图像输出模块、图像存储模块和损伤诊断模块。
其中，原始图像获取模块的输入端连接到核磁共振(mri)设备，并接受核磁共振设备发送的半月板影像。上述影像为通过核磁共振扫描设备获得的序列切片图像。
原始图像获取模块在接收到上述序列切片图像img1，img2，……，imgi后，将上述切片图像进行格式转换、灰度转换和去噪滤波。优选地，去噪滤波通过高斯滤波的方式进行。
原始图像获取模块的输出端连接到图像处理模块，上述模块用于对上述的序列切片图像img1，img2，……，imgi进行图像处理以形成半月板的3d图像。
图像处理模块进一步包括区域提取组件、边缘提取组件、中心计算组件、阈值调整组件和卷积神经网络组件。
其中，区域提取组件，用于在接收到序列切片图像img1，img2，……，imgi后，对每一张切片图像进行图像序列灰度信息分析，并获得可能包含半月板目标的最小计算区域，减小目标搜索范围，同时降低误提取概率。
具体地，上述灰度信息分析可以通过：1)标引切片图像中的灰度等级；2)将相邻的灰度等级相同的像素区域连通为大片区域；3)将切片图像中最外围的灰度等级最高(即亮度最低)的部分排除，剩余的部分可以保留为最小计算区域。可替换地，步骤3)还可以为将切片图像中最外围的灰度等级最高(即亮度最低)的部分排除后，在剩余区域中的矩阵像素中的每个像素的灰度和边缘特征进行聚类特征矩阵，并计算显著性因子，并根据显著性因子对像素聚类，像素聚类后的外部边缘则形成最小计算区域的外围轮廓。
边缘提取组件，用于在区域提取组件计算得到所有切片图像的最小计算区域后，接收上述序列切片图像img1，img2，……，imgi的最小计算区域的图像信息，并进一步确定半月板图像更准确的区域边缘信息。
具体地，每一个切片图像img1，img2，……，imgi都由x*y个像素点阵列构成，每一个像素点的坐标为(x，y)。图像函数f(x,y)在x，y轴两个方向的偏导数分别为：
hx(x,y)＝(f(x,y+1)-f(x,y)+f(x+1,y+1)-f(x+1,y))/2
hy(x,y)＝(f(x,y)-f(x+1,j)+f(x,y+1)-f(x+1,y+1))/2
图像灰度的阶度幅值为d(x,y)＝(hx2(x,y)+hy2(x,y))1/2；
阶度方向a(x,y)＝arctg(hx(x,y)/hy(x,y))
利用梯度方向a(x,y)将梯度角离散为圆周的四个扇区，分别编号为0到3，在梯度幅值图像中以点(x,y)为中心n×n(n＝3-6)的邻域内沿梯度方向与领域两个像素点相比较，如果(x,y)处幅度值小则标记为非边缘点，即得到二级边缘图像n(x,y)。
随后设定高低阈值，具体地，每一个切片图像img1，img2，……，imgi的上述高低阈值的计算方式为：
1)对于序列图片中的每一张，选择整个最小计算区域的像素中的最高灰度和最低灰度的平均值；
2)计算灰度高于上述平均值的所有像素的平均灰度g1以及低于上述平均值的所有像素的平均灰度g2；
3)获得中位阈值t＝(g1+g2)/2；
4)高阈值为高于中位阈值t的所有像素的平均灰度，低阈值为低于中位阈值t的所有像素的平均灰度。
当像素点灰度大于高阈值时，则该点为边缘点；当像素点小于低阈值时，该点一定不是边缘点；对于灰度在高低阈值之间的点，当其周围所有的像素点内存在大于高阈值的像素，则该点为边缘点，进而获得边缘图像函数e(x,y)。每一张序列切片图像img1，img2，……，imgi均具有边缘图像函数ei(x,y)。
中心计算组件，用于接收边缘提取组件提供的图像信息和边缘图像函数信息，并根据上述信息计算每一张序列切片图像img1，img2，……，imgi中的边缘函数图像围城的区域的灰度中心oi(k,l)。
具体地，该中心的计算方式为：
其中，当(x,y)位于目标区域i(x,y)＝1；当(x,y)位于背景区域。
阈值调整组件，用于接收中心计算组件传送的图像信息以及图像中心oi(k,l)的数据信息，并根据oi(k,l)数据信息制作切片时序-中心点x坐标和切片时序-中心点y坐标函数图。即函数图以序列图片的编号为横轴，以图片边缘区域内的中心点的x或y坐标为纵轴。通过人工智能数据分析，筛选上述切片时序-中心点x坐标和切片时序-中心点y坐标函数图中的函数曲线突变点。具体的数据分析方法可以为ks检验、小波分析法、t检验或者mann-kendal检验法等本领域公知的检测方法，此处不再赘述。
如果不存在突变点，则不需要反馈控制边缘提取组件调整高低阈值；如果存在突变点，则需要反馈边缘提取组件需要调整相应图像的高低阈值。
这是由于对于半月板图像的时序扫描以及人体组织的生长都是连续的，因此对于临近的扫面图像在灰度阈值选定合理的情况下，连续图像之间的灰度中心点应该也是连续变化，而非突变的。突变的形成主要是由于单张扫描图像的边界提取时高低阈值的选择不恰当，从而造成了边界偏移，从而形成了灰度中心点的突变。因此，在出现连续序列图片之间的灰度中心点的突变时，需要及时调整突变点对应的图片的高低阈值以及边界。
具体的阈值调整方式为：
1)通过人工智能数据分析，计算突变点偏离的方向α(+x,+y,-x,-y)和步长δ；具体的数据分析方法可以为ks检验、小波分析法、t检验或者mann-kendal检验法等。变点偏离的方向α(+x,+y,-x,-y)和步长δ具体可以为现有偏离数据点与数据分析中拟合的数据线上相应序列号的数据点之间的差距；
2)调取突变点对应的切片图像的最小计算面积图像以及现有的灰度中心点信息，将灰度中心点沿突变偏离方向的相反方向移动步长δ，获得修正灰度中心点o’(x,y)；
3)计算灰度偏差值δg为沿偏离方向α(+x,+y,-x,-y)相反的方向上的距离步长δ的所有像素的灰度差的加权平均值。修正高低阈值分别为原高低阈值加灰度偏差值δg；
4)边缘提取组件根据重新调整的高低阈值，进一步确定修正边界。
卷积神经网络组件，用于接受序列扫描图像的边界以及修正边界，并通过卷积神经网络进行三维图像的重建。
3d图像输出模块，用于将图像处理模块处理的图像进行三维显示，具体地，使用vtk中的gpu体绘制方法快速生成交互式的三维重建图像，可以进行旋转、缩放等操作，并用方向控件指示三维空间上的方向，以便观察。
图像存储模块，用于存储图像处理模块传送的处理图像。图像存储模块还可以预存半月板损伤图像信息，以供自动诊断中的人工智能学习。
损伤诊断模块，用于基于人工智能算法的半月板损伤诊断。
利用人工智能中神经网络的图像识别技术，识别半月板三维数据中的损伤部位与形态，并给出诊断结论。
实施例2基于mri的半月板三维重建方法
本发明提供的基于mri的半月板三维重建方法包括：
s1.接收到上述序列切片图像img1，img2，……，imgi后，将上述切片图像进行格式转换、灰度转换和去噪滤波。优选地，去噪滤波通过高斯滤波的方式进行。
s2.进行图像序列的图像处理：
s21.对每一张切片图像进行图像序列灰度信息分析，并获得可能包含半月板目标的最小计算区域，减小目标搜索范围，同时降低误提取概率。
具体地，上述灰度信息分析可以通过：1)标引切片图像中的灰度等级；2)将相邻的灰度等级相同的像素区域连通为大片区域；3)将切片图像中最外围的灰度等级最高(即亮度最低)的部分排除，剩余的部分可以保留为最小计算区域。可替换地，步骤3)还可以为将切片图像中最外围的灰度等级最高(即亮度最低)的部分排除后，在剩余区域中的矩阵像素中的每个像素的灰度和边缘特征进行聚类特征矩阵，并计算显著性因子，并根据显著性因子对像素聚类，像素聚类后的外部边缘则形成最小计算区域的外围轮廓。
s22.在区域提取组件计算得到所有切片图像的最小计算区域后，接收上述序列切片图像img1，img2，……，imgi的最小计算区域的图像信息，并进一步确定半月板图像更准确的区域边缘信息。
具体地，每一个切片图像img1，img2，……，imgi都由x*y个像素点阵列构成，每一个像素点的坐标为(x，y)。图像函数f(x,y)在x，y轴两个方向的偏导数分别为：
hx(x,y)＝(f(x,y+1)-f(x,y)+f(x+1,y+1)-f(x+1,y))/2
hy(x,y)＝(f(x,y)-f(x+1,j)+f(x,y+1)-f(x+1,y+1))/2
图像灰度的阶度幅值为d(x,y)＝(hx2(x,y)+hy2(x,y))1/2；
阶度方向a(x,y)＝arctg(hx(x,y)/hy(x,y))
利用梯度方向a(x,y)将梯度角离散为圆周的四个扇区，分别编号为0到3，在梯度幅值图像中以点(x,y)为中心n×n(n＝3-6)的邻域内沿梯度方向与领域两个像素点相比较，如果(x,y)处幅度值小则标记为非边缘点，即得到二级边缘图像n(x,y)。
随后设定高低阈值，具体地，每一个切片图像img1，img2，……，imgi的上述高低阈值的计算方式为：
1)对于序列图片中的每一张，选择整个最小计算区域的像素中的最高灰度和最低灰度的平均值；
2)计算灰度高于上述平均值的所有像素的平均灰度g1以及低于上述平均值的所有像素的平均灰度g2；
3)获得中位阈值t＝(g1+g2)/2；
4)高阈值为高于中位阈值t的所有像素的平均灰度，低阈值为低于中位阈值t的所有像素的平均灰度。
当像素点灰度大于高阈值时，则该点为边缘点；当像素点小于低阈值时，该点一定不是边缘点；对于灰度在高低阈值之间的点，当其周围所有的像素点内存在大于高阈值的像素，则该点为边缘点，进而获得边缘图像函数e(x,y)。每一张序列切片图像img1，img2，……，imgi均具有边缘图像函数ei(x,y)。
s23.接收边缘提取组件提供的图像信息和边缘图像函数信息，并根据上述信息计算每一张序列切片图像img1，img2，……，imgi中的边缘函数图像围城的区域的中心oi(k,l)。
具体地，该中心的计算方式为：
其中，当(x,y)位于目标区域i(x,y)＝1；当(x,y)位于背景区域。
接收中心计算组件传送的图像信息以及图像中心oi(k,l)的数据信息，并根据oi(k,l)数据信息制作切片时序-中心点x坐标和切片时序-中心点y坐标函数图。即函数图以序列图片的编号为横轴，以图片边缘区域内的中心点的x或y坐标为纵轴。通过人工智能数据分析，筛选上述切片时序-中心点x坐标和切片时序-中心点y坐标函数图中的函数曲线突变点。具体的数据分析方法可以为ks检验、小波分析法、t检验或者mann-kendal检验法等本领域公知的检测方法，此处不再赘述。
如果不存在突变点，则不需要反馈控制边缘提取组件调整高低阈值；如果存在突变点，则需要反馈边缘提取组件需要调整相应图像的高低阈值。
这是由于对于半月板图像的时序扫描以及人体组织的生长都是连续的，因此对于临近的扫面图像在灰度阈值选定合理的情况下，连续图像之间的灰度中心点应该也是连续变化，而非突变的。突变的形成主要是由于单张扫描图像的边界提取时高低阈值的选择不恰当，从而造成了边界偏移，从而形成了灰度中心点的突变。因此，在出现连续序列图片之间的灰度中心点的突变时，需要及时调整突变点对应的图片的高低阈值以及边界。
通过序列图像之间的灰度中点位置的比较，进一步保证了序列扫描图像中的边缘确定不会犹豫阈值选择不恰当出现突然的偏移，从而造成3d图像的错误。
具体的阈值调整方式为：
1)通过人工智能数据分析，计算突变点偏离的方向α(+x,+y,-x,-y)和步长δ；具体的数据分析方法可以为ks检验、小波分析法、t检验或者mann-kendal检验法等。变点偏离的方向α(+x,+y,-x,-y)和步长δ具体可以为现有偏离数据点与数据分析中拟合的数据线上相应序列号的数据点之间的差距；
2)调取突变点对应的切片图像的最小计算面积图像以及现有的灰度中心点信息，将灰度中心点沿突变偏离方向的相反方向移动步长δ，获得修正灰度中心点o’(x,y)；
3)计算灰度偏差值δg为沿偏离方向α(+x,+y,-x,-y)相反的方向上的距离步长δ的所有像素的灰度差的加权平均值。修正高低阈值分别为原高低阈值加灰度偏差值δg；
4)边缘提取组件根据重新调整的高低阈值，进一步确定修正边界。
s24.接受序列扫描图像的边界以及修正边界，并通过卷积神经网络进行三维图像的重建。
实施例3准确性验证
膝关节半月板损伤80例。
纳入标准：a.近1月内有膝关节疼痛，膝关节负重时扭转或者下蹲后疼痛加重；b.内或外侧关节间隙压痛；c.mcmurray试验及apley试验阳性；d.常规mri检查见半月板内与关节相通的高信号征象；e.关节镜手术中见半月板损伤。
排除标准：a.合并类风湿性关节炎；b.合并强直性脊柱炎；c.合并痛风性关节炎；d.合并化脓性关节炎；e.未行膝关节mri检查者；f.无关节镜手术指征者。
半月板损伤分类：按损伤部位分类：a.前角，b.体部，c.后角，d.多发性。按损伤形态分类：a.水平撕裂型(平行于胫骨平台，将半月板分为上下两部分)，b.纵形撕裂型(沿着半月板前后轴方向垂直于胫骨平台的撕裂)，c.放射状撕裂型(垂直于半月板前后轴方向垂直于胫骨平台的撕裂)，d.斜形撕裂型(斜于半月板前后轴方向垂直于胫骨平台的撕裂)，e.复杂撕裂型(包含2种以上上述撕裂)。
本项目软件诊断结果：利用本项目软件读取80例患膝mri的dicom数据，作出半月板损伤部位与形态的诊断。
关节镜诊断结果：对80例患膝行关节镜手术，记录关节镜下的损伤部位与形态诊断。
统计分析：使用spss17.0软件，kappa值为一致性指标，以关节镜诊断为标准，验证软件诊断的准确性。经验证，利用本发明提供的基于mri的半月板三维重建系统进行的诊断的准确性在85.9％。
本发明所提供的基于mri的半月板三维重建系统能够快速准确地提取膝关节mri中的半月板数据，生成三维重建图像；并能自动识别半月板损伤的部位及形态，作出初步诊断，提高医生工作效率。
以上显示和描述了本发明的基本原理、主要特征和优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其效物界定。
- 还没有人留言评论。精彩留言会获得点赞！