近视度数与眼轴长之间的关系(数值模拟版)

其实这又是一篇optiland的教程。
更新：第一版的顶点距离设定错了，设定成了1e-4，因为是从我其他程序里copy出来的，那个程序里在后续位置设定了框架镜。这一版修订成12mm，所以每毫米眼轴增长对应的近视度数结果略有差异。

一个常见的临床问题：患者眼轴增长了1mm，近视度数会增长多少？我们可以实际算一算。注意，我们这里使用的仍然是Gullstand模型眼，其实它并不适用于近视，所以只是一个半定量的计算，主要目的还是为了熟悉optiland的使用。

建立模型眼

当然我们先要导入必要的库

try:
    from optiland.optic import Optic
except:
    ! pip install optiland[torch]
    from optiland.optic import Optic
import numpy as np
import matplotlib.pyplot as plt
from optiland.materials import AbbeMaterial
from optiland import analysis, optimization, solves
from optiland.mtf import FFTMTF, SampledMTF
from optiland.solves.quick_focus import QuickFocusSolve

我们按照之前的教程：用optiland建立一个Gullstrand模型眼建立一个Gullstrand模拟眼。为了后续方便，这次使用python的“类”，也就是把上次一条一条写的模拟眼设定，给放在一起，以后方便调用。一个类里面，首先要有一个初始化的函数__init__

class Model_Eye(Optic):
    def __init__(self,
                 power=0.0,
                 pupil_size=4.0
                 ):
        super().__init__()

        # wavelength
        self.wavelength_F= 0.48613
        self.wavelength_d= 0.58756
        self.wavelength_C= 0.65627
        self.add_wavelength(value = self.wavelength_F)
        self.add_wavelength(value = self.wavelength_d)
        self.add_wavelength(value = self.wavelength_C)
        # material
        material_aqueous = AbbeMaterial(n=1.336981, abbe=52.658991)
        material_cornea = AbbeMaterial(n=1.376981, abbe=56.279936)
        material_lens = AbbeMaterial(n=1.419976, abbe= 51.226142)
        material_vitreous = AbbeMaterial(n=1.335982, abbe=53.342173)

        # aperture
        self.set_aperture(aperture_type="float_by_stop_size", value = pupil_size )

        # FOV
        self.set_field_type(field_type="angle")
        self.add_field(y=0)

        f = 1000/power if power != 0 else np.inf

        surface_index = 0
        self.add_surface(index = surface_index , radius = np.inf, thickness =np.inf, comment = 'Object')
        surface_index += 1
        self.add_surface(index = surface_index, comment = 'Phoropter',
                        surface_type = 'paraxial',
                        thickness = 12, f = f
                        )

        surface_index += 1
        self.add_surface(index = surface_index , comment = 'Cornea',
                        surface_type='standard',
                        radius = 7.70, thickness = 0.55,
                        material=material_cornea,
                        aperture=8)
        surface_index += 1
        self.add_surface(index = surface_index , comment = 'Aqueous',
                        surface_type='standard',
                        radius = 6.80, thickness = 3.6,
                        material=material_aqueous,
                        aperture=8)
        surface_index += 1
        self.add_surface(index = surface_index, comment = 'Pupil',
                        surface_type='standard',
                        radius = np.inf, thickness = 0,
                        material=material_aqueous,
                        aperture = pupil_size,
                        is_stop=True)
        surface_index += 1
        self.add_surface(index = surface_index , comment = 'Lens front',
                        surface_type='standard',
                        radius = 10.0, thickness = 0.546,
                        material=AbbeMaterial(n=1.386, abbe=50.23),
                        aperture=8)
        surface_index += 1
        self.add_surface(index = surface_index , comment = 'Lens front core',
                        surface_type='standard',
                        radius = 7.91, thickness = 2.419,
                        material=AbbeMaterial(n=1.406, abbe=50.23),
                        aperture=8)
        surface_index += 1
        self.add_surface(index = surface_index , comment = 'Lens back core',
                        surface_type='standard',
                        radius = -5.76, thickness = 0.635,
                        material=AbbeMaterial(n=1.386, abbe=50.23),
                        aperture=8)
        surface_index += 1
        self.add_surface(index = surface_index , comment = 'Vitreous',
                        surface_type='standard',
                        radius = -6.0,
                        thickness = 17.185,
                        material=material_vitreous,
                        aperture=8)
        surface_index += 1
        self.add_surface(index = surface_index , comment = 'Retina',
                        surface_type='standard',
                        radius = -12.0, thickness =0,
                        aperture=8)

通过对焦计算玻璃体腔厚度

这就是上一次建立模拟眼的过程。
然后，我们加入一个对焦的函数，让它来修改玻璃体厚度，这也在上一讲中讲过。这个函数也是放在类的定义中的。

class Model_Eye(Optic):
    def __init__(self,
                 power=0.0,
                 pupil_size=4.0
                 ):
        ......
    def def quick_focus(self):
        ......

具体内容也和上次差不多，只是定位surface_number的时候，用了整个总面数来推导，在optiland中是从0开始计数的，所以物体面是第0面，比如一共10个面，那么视网膜面就是第9号面，而玻璃体面就是再之前一个面，是总面数-2。

    def quick_focus(self):
        problem = optimization.OptimizationProblem()
        problem.add_variable(self, "thickness", surface_number = self.surface_group.num_surfaces-2 , min_val=10, max_val=30)

        input_data = {
            "optic": self,
            "surface_number": -1,
            "Hx": 0,
            "Hy": 0,
            "num_rays":16,
            "wavelength": self.wavelength_d, 
            "distribution": "uniform",
        }
        problem.add_operand(
            operand_type = "rms_spot_size",
            target=0.0,
            weight =10,
            input_data=input_data,
        )

        optimizer = optimization.OptimizerGeneric(problem)
        optimizer.optimize()

        # Update the vitreous humor thickness with the optimized value
        optimized_thickness = problem.variables[0].value
        self.surface_group.surfaces[8].thickness = optimized_thickness

计算眼轴长

在计算眼轴长的时候，可以把视网膜面的位置（position）减去角膜面的位置。当然也可以直接把各个面的厚度加起来，但没有这样来得优雅。@property是一个“装饰器”，这样读取眼轴长的时候，就可以写成Model_Eye.axial_length


    @property
    def axial_length(self):
        AL=self.surface_group.positions[-1]-self.surface_group.positions[-1-7]
        return AL

整合在一起

把模拟眼这个“类”整合在一起的结果就是这样的：

class Model_Eye(Optic):
    def __init__(self,
                 power=0.0,
                 pupil_size=4.0
                 ):
        super().__init__()

        # wavelength
        self.wavelength_F= 0.48613
        self.wavelength_d= 0.58756
        self.wavelength_C= 0.65627
        self.add_wavelength(value = self.wavelength_F)
        self.add_wavelength(value = self.wavelength_d)
        self.add_wavelength(value = self.wavelength_C)
        # material
        material_aqueous = AbbeMaterial(n=1.336981, abbe=52.658991)
        material_cornea = AbbeMaterial(n=1.376981, abbe=56.279936)
        material_lens = AbbeMaterial(n=1.419976, abbe= 51.226142)
        material_vitreous = AbbeMaterial(n=1.335982, abbe=53.342173)

        # aperture
        self.set_aperture(aperture_type="float_by_stop_size", value = pupil_size )

        # FOV
        self.set_field_type(field_type="angle")
        self.add_field(y=0)

        f = 1000/power if power != 0 else np.inf
        surface_index = 0
        self.add_surface(index = surface_index , radius = np.inf, thickness =np.inf, comment = 'Object')
        surface_index += 1
        self.add_surface(index = surface_index, comment = 'Phoropter',
                        surface_type = 'paraxial',
                        thickness = 1e-4, f = f
                        )

        surface_index += 1
        self.add_surface(index = surface_index , comment = 'Cornea',
                        surface_type='standard',
                        radius = 7.70, thickness = 0.55,
                        material=material_cornea,
                        aperture=8)
        surface_index += 1
        self.add_surface(index = surface_index , comment = 'Aqueous',
                        surface_type='standard',
                        radius = 6.80, thickness = 3.6,
                        material=material_aqueous,
                        aperture=8)
        surface_index += 1
        self.add_surface(index = surface_index, comment = 'Pupil',
                        surface_type='standard',
                        radius = np.inf, thickness = 0,
                        material=material_aqueous,
                        aperture = pupil_size,
                        is_stop=True)
        surface_index += 1
        self.add_surface(index = surface_index , comment = 'Lens front',
                        surface_type='standard',
                        radius = 10.0, thickness = 0.546,
                        material=AbbeMaterial(n=1.386, abbe=50.23),
                        aperture=8)
        surface_index += 1
        self.add_surface(index = surface_index , comment = 'Lens front core',
                        surface_type='standard',
                        radius = 7.91, thickness = 2.419,
                        material=AbbeMaterial(n=1.406, abbe=50.23),
                        aperture=8)
        surface_index += 1
        self.add_surface(index = surface_index , comment = 'Lens back core',
                        surface_type='standard',
                        radius = -5.76, thickness = 0.635,
                        material=AbbeMaterial(n=1.386, abbe=50.23),
                        aperture=8)
        surface_index += 1
        self.add_surface(index = surface_index , comment = 'Vitreous',
                        surface_type='standard',
                        radius = -6.0,
                        thickness = 17.185,
                        material=material_vitreous,
                        aperture=8)
        surface_index += 1
        self.add_surface(index = surface_index , comment = 'Retina',
                        surface_type='standard',
                        radius = -12.0, thickness =0,
                        aperture=8)

    def quick_focus(self):
        problem = optimization.OptimizationProblem()
        problem.add_variable(self, "thickness", surface_number = self.surface_group.num_surfaces-2 , min_val=10, max_val=30)

        input_data = {
            "optic": self,
            "surface_number": -1,
            "Hx": 0,
            "Hy": 0,
            "num_rays":16,
            "wavelength": self.wavelength_d, 
            "distribution": "uniform",
        }
        problem.add_operand(
            operand_type = "rms_spot_size",
            target=0.0,
            weight =10,
            input_data=input_data,
        )

        optimizer = optimization.OptimizerGeneric(problem)
        optimizer.optimize()

        # Update the vitreous humor thickness with the optimized value
        optimized_thickness = problem.variables[0].value
        self.surface_group.surfaces[8].thickness = optimized_thickness
    def reset_phoropter(self,power=0.0):
        self.surface_group.surfaces[1].f = 1000/self.power if power != 0 else np.inf

    @property
    def axial_length(self):
        AL=self.surface_group.positions[-1]-self.surface_group.positions[-1-7]
        return AL

计算近视度数与眼轴长之间的关系

当这个模拟眼设计好了以后，我们就可以代入不同的近视度数，计算其眼轴长，得到结果了。注意这并非一个严谨的测试，因为Gullstrand模型眼并非适用于近视。我们只是假定了除了玻璃体厚度，眼球的任何结构与近视发展都无关，这个假设其实并不成立。更合适的方法是使用近视专用的模型眼，比如Atchison模型眼，但它需要使用梯度折射率来描述晶状体，当前的optiland还不支持。

下面这一段代码，一口气画出了从-10D到+5D的度数与眼轴的关系，

power_list = np.linspace(-10,5,16)
AL_list = np.zeros(len(power_list))
for power in power_list:
    my_eye = Model_Eye(power = power, pupil_size= 4)
    my_eye.quick_focus()    
    AL_list[power_list==power]=my_eye.axial_length
plt.plot(AL_list, power_list,)
plt.scatter(AL_list, power_list)
plt.ylabel('Power (dB)')
plt.xlabel('Axial Length (mm)')
plt.show()

An image to describe post

可以看出几乎是线性的。

我们做一个简单的线性拟合：

# Fit a linear equation to the data
coefficients = np.polyfit(AL_list,power_list, 1)
# Print the coefficients (slope and intercept)
print(f"Slope: = {coefficients[0]:.2f}D")
x = np.linspace(21,29,100)
y = coefficients[0]*x+coefficients[1]
plt.plot(x,y, "--")
plt.scatter(AL_list, power_list)
plt.ylabel('Power (dB)')
plt.xlabel('Axial Length (mm)')
plt.show()

An image to describe post

计算得到斜率 = -2.67，也就是1mm眼轴增长，大致对应的是-2.67D近视度数。