CN106949968B - 一种基于频谱能量的数字全息图谐波探测与消除方法 - Google Patents
一种基于频谱能量的数字全息图谐波探测与消除方法 Download PDFInfo
- Publication number
- CN106949968B CN106949968B CN201710179861.4A CN201710179861A CN106949968B CN 106949968 B CN106949968 B CN 106949968B CN 201710179861 A CN201710179861 A CN 201710179861A CN 106949968 B CN106949968 B CN 106949968B
- Authority
- CN
- China
- Prior art keywords
- harmonic
- spectrum
- frequency
- harmonic wave
- coordinate
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Fee Related
Links
- 238000001228 spectrum Methods 0.000 title claims abstract description 157
- 238000000034 method Methods 0.000 title claims abstract description 44
- 238000001514 detection method Methods 0.000 title claims abstract description 22
- 230000009466 transformation Effects 0.000 claims abstract description 38
- 238000012545 processing Methods 0.000 claims abstract description 8
- 238000004422 calculation algorithm Methods 0.000 claims description 20
- 230000000737 periodic effect Effects 0.000 claims description 9
- 230000003595 spectral effect Effects 0.000 claims description 8
- 101100173586 Schizosaccharomyces pombe (strain 972 / ATCC 24843) fft2 gene Proteins 0.000 claims description 3
- ODKSFYDXXFIFQN-UHFFFAOYSA-M argininate Chemical compound [O-]C(=O)C(N)CCCNC(N)=N ODKSFYDXXFIFQN-UHFFFAOYSA-M 0.000 claims description 3
- 239000004744 fabric Substances 0.000 claims description 2
- 210000000746 body region Anatomy 0.000 claims 1
- 238000003702 image correction Methods 0.000 claims 1
- 238000013461 design Methods 0.000 abstract description 7
- 238000001093 holography Methods 0.000 abstract description 4
- 238000012937 correction Methods 0.000 description 7
- 230000009467 reduction Effects 0.000 description 7
- 230000008030 elimination Effects 0.000 description 5
- 238000003379 elimination reaction Methods 0.000 description 5
- 238000002474 experimental method Methods 0.000 description 4
- 239000004615 ingredient Substances 0.000 description 4
- 238000005259 measurement Methods 0.000 description 4
- 238000004364 calculation method Methods 0.000 description 3
- 230000000694 effects Effects 0.000 description 3
- 241000208340 Araliaceae Species 0.000 description 2
- 235000005035 Panax pseudoginseng ssp. pseudoginseng Nutrition 0.000 description 2
- 235000003140 Panax quinquefolius Nutrition 0.000 description 2
- 238000005094 computer simulation Methods 0.000 description 2
- 238000001914 filtration Methods 0.000 description 2
- 235000008434 ginseng Nutrition 0.000 description 2
- 230000003287 optical effect Effects 0.000 description 2
- 238000010587 phase diagram Methods 0.000 description 2
- 230000008569 process Effects 0.000 description 2
- 238000003672 processing method Methods 0.000 description 2
- 230000007547 defect Effects 0.000 description 1
- 230000007812 deficiency Effects 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 238000003384 imaging method Methods 0.000 description 1
- 230000002452 interceptive effect Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000000149 penetrating effect Effects 0.000 description 1
- 230000010363 phase shift Effects 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01J—MEASUREMENT OF INTENSITY, VELOCITY, SPECTRAL CONTENT, POLARISATION, PHASE OR PULSE CHARACTERISTICS OF INFRARED, VISIBLE OR ULTRAVIOLET LIGHT; COLORIMETRY; RADIATION PYROMETRY
- G01J3/00—Spectrometry; Spectrophotometry; Monochromators; Measuring colours
- G01J3/28—Investigating the spectrum
- G01J3/45—Interferometric spectrometry
-
- G—PHYSICS
- G03—PHOTOGRAPHY; CINEMATOGRAPHY; ANALOGOUS TECHNIQUES USING WAVES OTHER THAN OPTICAL WAVES; ELECTROGRAPHY; HOLOGRAPHY
- G03H—HOLOGRAPHIC PROCESSES OR APPARATUS
- G03H1/00—Holographic processes or apparatus using light, infrared or ultraviolet waves for obtaining holograms or for obtaining an image from them; Details peculiar thereto
- G03H1/04—Processes or apparatus for producing holograms
- G03H1/16—Processes or apparatus for producing holograms using Fourier transform
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01J—MEASUREMENT OF INTENSITY, VELOCITY, SPECTRAL CONTENT, POLARISATION, PHASE OR PULSE CHARACTERISTICS OF INFRARED, VISIBLE OR ULTRAVIOLET LIGHT; COLORIMETRY; RADIATION PYROMETRY
- G01J3/00—Spectrometry; Spectrophotometry; Monochromators; Measuring colours
- G01J3/28—Investigating the spectrum
- G01J3/45—Interferometric spectrometry
- G01J2003/452—Interferometric spectrometry with recording of image of spectral transformation, e.g. hologram
Landscapes
- Physics & Mathematics (AREA)
- Spectroscopy & Molecular Physics (AREA)
- General Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- Holo Graphy (AREA)
Abstract
本发明属于数字全息技术领域，公开了一种基于频谱能量的数字全息图谐波探测与消除方法，包括图像的傅里叶变换、谐波频谱的寻找、读取谐波频谱的坐标和复数值、计算谐波在水平竖直两个方向上的频率和干扰强度系数以及初相位、利用参数构建谐波分布、图像校正六个步骤完成有谐波干扰的图像的处理。本发明不需要设计滤波器，只利用零级傅里叶变换谱和谐波谱的能量关系即可；无需获得谐波对应频谱的范围；不会对原物光的信息有任何影响；不需要对频谱图进行逆傅里叶变换；利用能量关系获得谐波参数后，构建谐波强度分布，直接相减消掉谐波，计算简单。
Description
技术领域
本发明属于数字全息技术领域，尤其涉及一种基于频谱能量的数字全息图谐波探测与消除方法。
背景技术
利用数字全息技术进行显微成像或精密测量时，数字全息图的质量对于还原像的质量或测量精度影响很大，激光散斑、参考光倾斜、谐波的存在会大幅降低成像质量或测量精度。数字全息需要物光与参考光干涉记录干涉图，相对于普通照相，记录光路中需要更多的光学元件，如分束器、合束器、空间滤波器等。光波经过这些器件时，其前后表面会发生反射、折射等现象。在棱镜或其他镜面之间多次反射的光波会发生干涉，给全息图附加周期性的正余弦分布的直条纹，被称为谐波。这种谐波是由于分束器两平面不严格平行，或某个或几个光学器件表面不平行造成的多次反射的具有一定夹角的平面波干涉引起的。这一干扰信息会严重干扰全息图的准确定位并造成全息再现时最终物光波的畸变或测量失败。要消除这些谐波的影响，需要对全息图进行预处理，消掉全息图中的这些谐波条纹。现有全息图预处理技术使用傅里叶变换在频谱面上进行滤波消掉谐波的频谱后再逆傅里叶变换获得无谐波的全息图。这一谐波消除方法需要设计相应的频谱滤波器进行滤波，去掉谐波频谱分量。这一滤波方法无法准确消除谐波成分，会剩余谐波分量，或者误伤有用的成分。而且，滤波器的滤波范围和滤波幅度也会影响滤波效果：滤波幅度较小时，不能完全除掉谐波信息，滤波幅度较大时，会增加错误信息。滤波范围小时，无法消除谐波影响；范围较大时，也会滤掉有用的物光频谱信息。滤波法的第二个技术缺陷是计算工作量大，滤波后还要进行一次逆傅里叶变换，获得没有谐波的全息图，两次傅里叶变换比一次傅里叶变换增大一倍的计算量。这也导致了滤波法的第三个缺陷：无法应用到连续的视频图像的谐波消除中。由于傅里叶变换计算量大，每幅图像都进行两次傅里叶变换需要大量的计算时间，即使是一次傅里叶变换，通常情况下也无法达到视频播放的速度。
综上所述，现有技术存在的问题是：
现有全息图预处理技术存在设计滤波器过程繁琐，成本高，寻找频谱困难；难以确定滤波的范围；难以确定滤波幅度；容易丢失有用的物光信息；需要对滤波后的频谱图进行逆傅里叶变换获得处理后的全息图，增加一倍的计算量；无法应用到在线连续图像处理中。
发明内容
针对现有技术存在的问题，本发明提供了一种基于频谱能量的数字全息图谐波探测与消除方法。
本发明是这样实现的，一种基于频谱能量的数字全息图谐波探测与消除方法，所述基于频谱能量的数字全息图谐波探测与消除方法包括图像的傅里叶变换、谐波频谱的寻找、读取谐波频谱的坐标和复数值、计算谐波在水平竖直两个方向上的频率和干扰强度系数及初相位、利用参数构建谐波分布、图像校正六个步骤完成有谐波干扰的图像的处理。
进一步，所述基于频谱能量的数字全息图谐波探测与消除方法具体包括以下步骤:
第一步，图像的傅里叶变换：使用计算机内的快速二维傅里叶变换算法完成对图像的傅里叶变换；谐波在原图像上是有一定偏置的正余弦周期性分布，其傅里叶频谱的模(强度)是三个对称分布的δ函数，其中一个在频谱面上的原点，另外两个关于原点对称分布；
第二步，谐波频谱的寻找：由于δ函数是无限小的点，在频谱面上的谐波的频谱是三个很亮的点，零频的点在原点，另两个点对称分布在原点的两侧；先分析出只包含该亮点的大体区域，再用依次求最大值的算法寻找；
第三步，读取谐波频谱的坐标和复数值：
用寻找最大值的方法得到谐波的非零频谱的水平和竖直方向的频谱坐标，分别设为u1和v1；再读出该像素频谱函数值的实部和虚部，分别设为Er和Ei，设其模为E1；另一非零频频谱坐标对称分布，频谱函数值相等，不用读取；同样，原点的频谱坐标为(0，0)，无需读取，但需要读取原点的频谱的函数值，这一函数值对应全息图的零级频谱，设为E0，其函数值为正实数；
第四步，计算谐波在水平竖直两个方向上的频率和干扰强度系数及初相位：
假设谐波频谱在频谱面上水平方向上的坐标为u1，在竖直方向上的频谱坐标为v1；求出水平方向的谐波频率和在竖直方向上的谐波频率fx和fy；
假设谐波和全息图零级项对应的频谱函数值的模分别为E1和E0，则由这两个参数求得谐波的干扰强度系数m；
假设谐波频谱对应函数值的实部和虚部分别设为Er和Ei；求出谐波对应的初相位
第五步，利用参数构建谐波分布：
利用得到的水平方向的谐波频率和在竖直方向上的谐波频率和干扰强度系数及初相位四个参数组建谐波的强度分布；
第六步，图像校正：
求出谐波的强度分布后，通过图像相减完全去掉强度图中的谐波成分，得到消除谐波影响后的干涉图。
进一步，所述图像的傅里叶变换具体包括：
设图像的强度分布为I(x,y)，使用计算机内的快速二维傅里叶变换算法完成对图像的傅里叶变换，即：
E(u,v)＝FFT2[I(x,y)] (1)
在这些频谱中，含有物光信息的频谱是连续的，且分布范围较大；谐波在图像上是有一定偏置的正余弦函数分布；谐波的傅里叶频谱的模是三个对称分布的δ函数，其中一个在频谱面上的原点，另外两个关于原点对称分布。
进一步，用依次求最大值的算法寻找频谱坐标的方法包括：
对全息图，零频分量的能量远大于高频分量；对数字全息图进行傅里叶变换后，零级分量最强，利用数值算法对整个频谱面复数组取模，再求最大值，只能得到零级频谱；
先将零级频谱赋值为0，再找最大值，找到两个谐波谱中的一个，记下其坐标和复数值；最后把零频的值再恢复到原值。
现有技术是通过人眼观察频谱位置，依据观察的位置设计滤波器，比较繁琐。
进一步，所述一定偏置来自于光的直透项或常数项。
进一步，所述计算谐波在水平竖直两个方向上的频率和干扰强度系数，以及初相位，具体包括：
假设谐波频谱在频谱面上水平方向上的坐标为u1，在竖直方向上的频谱坐标为v1，求出水平方向的谐波频率为：
fx＝u1 (2)
在竖直方向上的谐波频率为：
fy＝v1 (3)
假设谐波和全息图的零级项对应的频谱函数值的模分别为E1和E0，则由这两个参数求得谐波的干扰强度系数为：
m＝2E1/(E0-2E1) (4)。
假设谐波对应频谱的频谱函数值的实部与虚部分别为Er和Ei，则谐波的初相位是：
式中arg()表示对括号内的复数取幅角。
进一步，利用所述参数构建谐波分布，具体包括：
利用公式(2)、(3)、(4)、(5)得到的四个参数组建谐波的强度分布，方程如下：
式中<>符号表示对所有像素取平均值，坐标x和y以像素为坐标单位。
进一步，所述图像校正，具体包括：
由公式(6)求出谐波的强度分布后，通过图像相减完全去掉强度图I中的谐波成分：
I1就是消除谐波影响后的干涉图；
存在多个谐波或连续谱时，依次找到其频谱坐标或连续的频谱坐标范围，记录坐标和相应的频谱函数值，求出每一项谐波的干扰强度系数和初相位，构建每一项强度分布，最后依次从带有谐波的图像中减掉；
对于非谐波的周期性干扰，展开成周期性谐波叠加的形式进行求解消除非谐波的影响；对于非周期性非简谐波干扰，通过傅里叶级数展开成不同谐波的叠加来求解。
本发明的优点及积极效果为：
不需要设计滤波器，只利用零级傅里叶变换谱和谐波谱的能量关系和谐波谱的复数值即可。
本发明无需获得谐波对应频谱的范围。本发明能够精确消除谐波成分，不会对原物光的信息有任何影响。本发明不需要对滤波后的频谱图进行逆傅里叶变换，可以使计算速度提高近50％。本发明利用能量关系获得谐波参数后，构建谐波强度分布，直接相减消掉谐波，计算简单。如果是连续记录的图像，只对一幅图像进行一次傅里叶变换，求出谐波参数后，可以对多幅带有谐波的图像进行快速消除，需要的运算只是一次图像相减。
本发明计算出谐波参数后，反演出谐波强度分布，直接进行图像相减，完成谐波的消除。只需对一幅图像进行一次傅里叶变换即可，不需要对所有图像进行傅里叶变换和逆傅里叶变换，大大节省了计算时间。特别是在视频显示方面，滤波法需要对大量图片进行傅里叶变换和逆傅里叶变换，计算负担大，而本发明中提出的方法，只进行一次傅里叶变换后，直接进行图片相减即可，极大地减轻了视频处理的计算负担。本发明的方法不仅适用于全息图谐波的探测与消除，也适用于其它数字图像和视频的谐波探测与处理，处理方法相同。
附图说明
图1是本发明实施例提供的基于频谱能量的数字全息图谐波探测与消除方法流程图。
图2是本发明实施例提供的原始物光强度和相位分布图；
图中：(a)强度分布；(b)相位分布。
图3是本发明实施例提供的棱镜前后表面平行与不平行情况对比图；
图中：(a)平行时全息图；(b)不平行时全息图。
图4是本发明实施例提供的傅里叶频谱图。
图5是本发明实施例提供的校正后的全息图。
图6是本发明实施例提供的未校正还原像强度与相位分布图；
图中：(a)强度图；(b)相位图。
图7是本发明实施例提供的校正后还原像强度与相位分布图；
图中：(a)强度图；(b)相位图。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
本发明实施例提供的基于频谱能量的数字全息图谐波探测与消除方法包括图像的傅里叶变换、谐波频谱的寻找、读取谐波频谱的坐标和复数值、计算谐波在水平竖直两个方向上的频率和干扰强度系数及初相位参数、利用参数构建谐波分布、图像校正六个步骤完成有谐波干扰的图像的处理。
下面结合附图对本发明的应用原理作详细描述。
如图1所示，本发明实施例提供的基于频谱能量的数字全息图谐波探测与消除方法包括以下步骤：
S101:图像的傅里叶变换：使用计算机内的快速二维傅里叶变换算法完成对图像的傅里叶变换；谐波在图像上是有一定偏置(来自于光的直透项或常数项)的正余弦周期性分布，其傅里叶频谱的模是三个对称分布的δ函数，其中一个在频谱面上的原点，另外两个关于原点对称分布。
S102:谐波频谱的寻找：由于δ函数是无限小的点，在频谱面上的谐波的频谱的模是三个很亮的点，零频的点在原点，另两个点对称分布在原点的两侧；先观察出包含该亮点的大体区域，再用求最大值的算法寻找。
S103:读取谐波频谱的坐标和复数值：用依次寻找最大值的方法得到谐波的非零频谱的水平和竖直方向的频谱坐标，分别设为u1和v1；再读出该像素频谱函数值，其实部和虚部分别设为Er和Ei，则该频谱函数的模为复数值的模E1＝|Er+jEi|，也就是该频谱的能量或强度；另一非零频频谱坐标对称分布，频谱函数值相等，不用读取；同样，原点的频谱坐标为(0，0)，无需读取，但需要读取原点的频谱的函数值，这一函数值对应全息图的零级频谱，为实数，设为E0。
S104:计算谐波在水平竖直两个方向上的频率和干扰强度系数及初相位：假设谐波频谱在频谱面上水平方向上的坐标为u1，在竖直方向上的频谱坐标为v1；求出水平方向的谐波频率和在竖直方向上的谐波频率；假设谐波和全息图零级频谱对应的频谱函数值的模分别为E1和E0，则由这两个参数求得谐波的干扰强度系数m；假设谐波频谱对应函数值的实部和虚部分别设为Er和Ei；求出谐波对应的初相位
S105:利用参数构建谐波分布：利用得到的参数求出水平方向的谐波频率和在竖直方向上的谐波频率和干扰强度系数及初相位四个参数组建谐波的强度分布。
S106:图像校正：求出谐波的强度分布后，通过图像相减完全去掉强度图中的谐波成分，得到消除谐波影响后的干涉图。
用依次求最大值的算法寻找频谱坐标的方法包括：
对全息图，零频分量的能量远大于高频分量；对数字全息图进行傅里叶变换后，零级分量最强，利用数值算法对整个频谱面复数组取模再求最大值，只能得到零级频谱；
先将零级频谱赋值为0，再找最大值，找到两个谐波谱中的一个，记下其坐标和复数值；最后把零频的值再恢复到原值。
现有技术是通过人眼观察频谱位置，依据观察的位置设计滤波器，比较繁琐。
下面结合具体实施例对本发明的应用原理作进一步描述。
本发明实施例提供的基于频谱能量的数字全息图谐波探测与消除方法，包括：
第一步，图像的傅里叶变换：
设图像的强度分布为I(x,y)，使用计算机内的快速二维傅里叶变换算法完成对图像的傅里叶变换，即：
E(u,v)＝FFT2[I(x,y)] (1)
在这些频谱中，含有物光信息的频谱是连续的，且分布范围较大。谐波在图像上是有一定偏置量(来自于光的直透项或常数项)的正余弦周期性分布，其傅里叶频谱的模是三个对称分布的δ函数，其中一个在频谱面上的原点，另外两个关于原点对称分布。
第二步，谐波频谱的寻找。
由于δ函数是无限小的点，在频谱面上的谐波的频谱的模是三个很亮的点，零频的点在原点，另两个点对称分布在原点的两侧。先观察出只包含该亮点的大体区域，再用求最大值的算法寻找。
第三步，读取谐波频谱的坐标和复数值。
用寻找最大值的方法得到谐波的非零频谱的水平和竖直方向的频谱坐标，分别设为u1和v1。再读出该频谱函数值，其实部和虚部分别设为Er和Ei，则该频谱函数的模为复数值的模E1＝|Er+jEi|，也就是该频谱的能量或强度。另一非零频谱的坐标对称分布，频谱函数值相等，不用读取。同样，原点的频谱坐标为(0，0)，无需读取，但需要读取原点的频谱的函数值，这一函数值对应全息图的零级频谱，为实数，设为E0。
第四步，计算谐波在水平竖直两个方向上的频率和干扰强度系数及初相位参数。
假设谐波频谱在频谱面上水平方向上的坐标为u1，在竖直方向上的频谱坐标为v1，可求水平方向的谐波频率为：
fx＝u1 (2)
在竖直方向上的谐波频率：
fy＝v1 (3)
假设在谐波和全息图的零级项对应的频谱函数的模分别为E1和E0，则可以由这两个参数求得谐波的干扰强度系数为：
m＝2E1/(E0-2E1) (4)
假设谐波频谱复数值对应的实部和虚部分别为Er和Ei，则可以求出与谐波对应的初相位为：
式中“arg()”表示对复数取幅角。
第五步，利用参数构建谐波分布。
利用公式(2)、(3)、(4)、(5)得到的四个参数组建谐波的强度分布方程如下：
式中<>符号表示对所有像素取平均值，坐标x和y以像素为坐标单位。
第六步，图像校正。
由公式(6)求出谐波的强度分布后，就可以通过图像相减完全去掉强度图I中的谐波成分：
I1就是消除谐波影响后的干涉图。
下面结合模拟实验对本发明的应用效果作详细的描述。
在广义相移数字全息实验中，参考光一般是滤波准直后的平面光波，如果棱镜前后两表面不严格平行，透射过分束器棱镜时会形成等厚干涉斜条纹。随着参考光强度的增加，干涉条纹会越强，意味着在全息图的记录中干扰信息会更强，在全息图还原时干扰信息对还原结果的影响也会越大，尤其是单步还原和两步还原算法，会导致还原结果严重扭曲，甚至无法还原出结果。
本发明不需要更多的操作也不需要采集更多的数据，用两步还原算法所需的两幅全息图中的一幅就可以得到需要的信息，通过计算机运算就可以有效降低棱镜前后两表面不平行造成的误差。接下来的内容中，本发明将对这种新算法进行计算机模拟并比较有校正与无校正的两种情况下还原像的质量。
在计算机模拟实验中设置物平面上有一球面波。模拟的物平面和像平面的大小都是512×512像素，像素尺寸为15×15微米。如图2所示，(a)为物光强度图，振幅分布满足高斯分布，(b)为物光相位。原始物光经过菲涅尔衍射后与参考光在记录面上形成干涉全息图，见图3(a)。由于棱镜前后表面不平行，导致在全息图上形成与全息图无关的斜条纹，如图3(b)所示：
图3(a)中棱镜前后表面平行，可以看到在全息图上没有谐波。图3(b)中棱镜前后表面存在很小的倾斜角度，造成了全息图上干扰性的斜条纹。在模拟实验中，本发明设置的相移值为1.5rad、干扰强度系数m值设置为0.2、谐波初相位设置为3rad，对全息图图3(b)进行傅里叶变换，其频谱函数值的强度分布如图4所示。
图4中，频谱图上有两个极大值点。由于零级谱太强，在显示图4时，零级谱被我们去掉了。利用这两个点的坐标值和谐波谱的模及幅角便可以求出谐波的水平方向和竖直方向的空间频率，以及校正公式中需要的m值，初相位
m的计算值为0.2032，初相位的计算值为3.0018rad，这与本发明的设定值误差很小，将fx、fy、m、的值代入校正算法中，可以将有干扰信息的全息图中干扰信息有效的消除，图5给出了一幅校正后的全息图，完全消除了谐波的影响。下面本发明将对无校正时的还原图像(图6)与校正后的还原结果(图7)进行对比。
在图6与图7之间的对比中，可以看出未校正的还原物光波无论是强度分布还是相位分布都受到了很大影响，尤其是相位分布图发生了严重的扭曲，而在校正后的物光波还原图中，强度分布与相位分布都得到了很好的校正。
本发明不需要设计滤波器，只利用全息图零级傅里叶变换谱和谐波谱的能量关系即可。
本发明无需获得谐波对应频谱的范围。本发明不会对原物光的信息有任何影响。本发明不需要对滤波后的频谱图进行逆傅里叶变换。本发明利用能量关系获得谐波参数后，构建谐波强度分布，直接相减消掉谐波，计算简单。
本发明计算出谐波参数后，反演出谐波强度分布，直接进行图像相减，完成谐波的消除。只需对一幅图像进行一次傅里叶变换即可，不需要对所有图像进行傅里叶变换和逆傅里叶变换，大大节省了计算时间，至少节省50％以上的计算资源。特别是在视频显示方面，滤波法需要对大量图片进行傅里叶变换和逆傅里叶变换，计算负担大，而本发明中提出的方法，只进行一次傅里叶变换后，直接进行图片相减即可，极大地减轻了视频处理的计算负担。本发明的方法不仅适用于全息图像谐波的探测与消除，也适用于其它数字图像和视频谐波的探测与处理，处理方法相同。
存在多个谐波或连续谱时，依次找到其频谱坐标或连续的频谱坐标范围，记录坐标和相应的频谱函数值，求出每一项谐波的干扰强度系数和初相位，构建每一项强度分布，最后依次从带有谐波的图像中减掉即可。
对于非谐波的周期性干扰，可以展开成周期性谐波叠加的形式进行求解消除其影响。对于非周期性非简谐波干扰，可以通过傅里叶级数展开成不同谐波的叠加来求解。
以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。
Claims (6)
1.一种基于频谱能量的数字全息图谐波探测与消除方法，其特征在于，所述基于频谱能量的数字全息图谐波探测与消除方法包括图像的傅里叶变换、谐波频谱的寻找、读取谐波频谱的坐标和复数值、计算谐波频谱在水平竖直两个方向上的频率和在水平竖直两个方向上的干扰强度系数以及初相位、利用参数构建谐波分布、图像校正六个步骤完成有谐波干扰的图像的处理；
所述基于频谱能量的数字全息图谐波探测与消除方法具体包括以下步骤:
第一步，图像的傅里叶变换：使用计算机内的快速二维傅里叶变换算法完成对图像的傅里叶变换；谐波在图像上是有一定偏置的正余弦周期性分布的函数，其傅里叶频谱函数值的模是三个对称分布的δ函数，其中一个在频谱面上的原点，另外两个关于原点对称分布；
第二步，谐波频谱的寻找：由于δ函数是无限小的点，在频谱面上的谐波的频谱是三个很亮的点，零频的点在原点，另两个点对称分布在原点的两侧；先分析出只包含该亮点的大体区域，再用依次求最大值的算法寻找；
第三步，读取谐波频谱的坐标和复数值：
用寻找最大值的方法得到谐波的其中一个非零频谱的水平和竖直方向的频谱坐标，分别设为u1和v1；再读出对应该频谱坐标的频谱函数值，其实部和虚部分别设为Er和Ei，则该频谱的模为复数值的模E1＝|Er+jEi|，也就是该频谱的能量或强度；另一非零频谱坐标对称分布，频谱函数值相等，不用读取；同样，原点的频谱坐标为(0，0)，无需读取，但需要读取原点的频谱的函数值，这一函数值对应全息图的零级频谱，为正实数，设为E0；
第四步，计算谐波在水平竖直两个方向上的频率和干扰强度系数：
假设谐波频谱在频谱面上水平方向上的坐标为u1，在竖直方向上的频谱坐标为v1；求出水平方向的谐波频率和在竖直方向上的谐波频率fx和fy；
假设谐波对应的非零频谱和全息图的零级频谱的函数值的模分别为E1和E0，则由这两个参数求得谐波的干扰强度系数m；
假设谐波频谱对应函数值的实部和虚部分别设为Er和Ei；求出谐波对应的初相位
第五步，利用参数构建谐波分布：
利用得到的水平方向的谐波频率和在竖直方向上的谐波频率和干扰强度系数以及初相位四个参数组建谐波的强度分布；
第六步，图像校正：
求出谐波的强度分布后，通过图像相减完全去掉强度图中的谐波成分，得到消除谐波影响后的干涉图。
2.如权利要求1所述的基于频谱能量的数字全息图谐波探测与消除方法，其特征在于，所述图像的傅里叶变换具体包括：
设图像的强度分布为I(x,y)，使用计算机内的快速二维傅里叶变换算法完成对图像的傅里叶变换，即：
E(u,v)＝FFT2[I(x,y)] (1)
在这些频谱中，含有物光信息的频谱是连续的，且分布范围较大；谐波在原图像上是具有一定偏置的正余弦函数分布；谐波的傅里叶频谱的模是三个对称分布的δ函数，其中一个在频谱面上的原点，另外两个关于原点对称分布。
3.如权利要求1所述的基于频谱能量的数字全息图谐波探测与消除方法，其特征在于，用依次求最大值的算法寻找频谱坐标方法包括：
对全息图，零频分量的能量远大于高频分量；对数字全息图进行傅里叶变换后，零级分量最强，利用数值算法对整个频谱面的数组求最大值，只能得到零级频谱；
先将零频赋值为0，再找最大值，找到两个谐波谱中的一个，记下其坐标和复数值；最后把零频的值再恢复到原值。
4.如权利要求2所述的基于频谱能量的数字全息图谐波探测与消除方法，其特征在于，所述一定偏置来自于光的直透项或常数项。
5.如权利要求1所述的基于频谱能量的数字全息图谐波探测与消除方法，其特征在于，所述计算谐波频谱在水平竖直两个方向上的频率和在水平竖直两个方向上的干扰强度系数，以及初相位，具体包括：
假设谐波频谱在频谱面上水平方向上的坐标为u1，在竖直方向上的频谱坐标为v1，求出水平方向的谐波频率为：
fx＝u1 (2)
在竖直方向上的谐波频率为：
fy＝v1 (3)
假设谐波频谱的模为E1和全息图零级频谱的模为E0，，则由这两个参数求得谐波的干扰强度系数为：
m＝2E1/(E0-2E1) (4)
假设谐波频谱复数值对应的实部和虚部分别为Er和Ei，则可以求出与谐波对应的初相位为：
式中“arg()”表示对复数取幅角；
利用所述参数构建谐波分布，具体包括：
利用公式(2)、(3)、(4)、(5)得到的四个参数组建谐波的强度分布方程如下：
式中<>符号表示对所有像素取平均值，坐标x和y以像素为坐标单位。
6.如权利要求1所述的基于频谱能量的数字全息图谐波探测与消除方法，其特征在于，所述图像校正，具体包括：
由公式求出谐波的强度分布后，通过图像相减完全去掉强度图I中的谐波成分：
I1就是消除谐波影响后的干涉图；
存在多个谐波或连续谱时，依次找到其频谱坐标或连续的频谱坐标范围，记录坐标和相应的频谱函数值，求出每一项谐波的干扰强度系数和初相位，构建每一项强度分布，最后依次从带有谐波的图像中减掉；
对于非谐波的周期性干扰，展开成周期性谐波叠加的形式进行求解消除非谐波的影响；对于非周期性非简谐波干扰，通过傅里叶级数展开的方法展开成不同谐波的叠加来求解。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710179861.4A CN106949968B (zh) | 2017-03-23 | 2017-03-23 | 一种基于频谱能量的数字全息图谐波探测与消除方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710179861.4A CN106949968B (zh) | 2017-03-23 | 2017-03-23 | 一种基于频谱能量的数字全息图谐波探测与消除方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN106949968A CN106949968A (zh) | 2017-07-14 |
| CN106949968B true CN106949968B (zh) | 2019-01-22 |
Family
ID=59472861
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201710179861.4A Expired - Fee Related CN106949968B (zh) | 2017-03-23 | 2017-03-23 | 一种基于频谱能量的数字全息图谐波探测与消除方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN106949968B (zh) |
Families Citing this family (10)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108802115A (zh) * | 2018-06-12 | 2018-11-13 | 贵州省产品质量监督检验院 | 一种基于大数据的食品检测肉类含水量检测方法及系统 |
| CN109002382A (zh) * | 2018-07-13 | 2018-12-14 | 广东水利电力职业技术学院（广东省水利电力技工学校） | 一种服务器主板监测系统及监测方法、信息处理终端 |
| CN108983578A (zh) * | 2018-08-17 | 2018-12-11 | 重庆工业职业技术学院 | 一种基于多媒体的虚拟投影汽车维修展示方法及系统 |
| CN109223060A (zh) * | 2018-08-23 | 2019-01-18 | 荆门市第二人民医院 | 一种外科手术中辅助装置的控制系统及控制方法 |
| CN111025939A (zh) * | 2019-11-28 | 2020-04-17 | 徐州华邦专用汽车有限公司 | 一种用于汽车配件制造的切割系统 |
| CN111136494A (zh) * | 2019-12-24 | 2020-05-12 | 重庆兰羚天和科技股份有限公司 | 一种横切机板剪可转动刀头系统 |
| CN113218337B (zh) * | 2021-04-29 | 2022-08-19 | 苏州天准软件有限公司 | 条纹图像能量提取方法、存储介质和能量提取系统 |
| CN115395967B (zh) * | 2021-05-25 | 2024-08-09 | 瑞昱半导体股份有限公司 | 发射器及其校正方法 |
| CN114500311B (zh) * | 2022-02-21 | 2024-05-14 | 普源精电科技股份有限公司 | 分辨率带宽可调节的频谱分析方法、装置及计算机设备 |
| CN114998126B (zh) * | 2022-05-24 | 2024-05-24 | 中国人民解放军国防科技大学 | 一种低波段星载sar图像电离层幅度闪烁条纹去除方法 |
Family Cites Families (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102023311A (zh) * | 2010-08-10 | 2011-04-20 | 中国石油大学(华东) | 地层的品质因子谱及其求取方法 |
| EP3186616B1 (en) * | 2014-08-28 | 2019-12-18 | The Regents Of The University Of Colorado | Coherent diffractive imaging with arbitrary angle of incidence |
| CN106092158B (zh) * | 2016-08-19 | 2018-08-21 | 北京理工大学 | 物理参数估计方法、装置和电子设备 |
-
2017
- 2017-03-23 CN CN201710179861.4A patent/CN106949968B/zh not_active Expired - Fee Related
Also Published As
| Publication number | Publication date |
|---|---|
| CN106949968A (zh) | 2017-07-14 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN106949968B (zh) | 一种基于频谱能量的数字全息图谐波探测与消除方法 | |
| JP5130311B2 (ja) | 波面の位相情報を回復するためのシステムおよび方法 | |
| JP5467321B2 (ja) | ３次元形状計測方法および３次元形状計測装置 | |
| WO2010077675A2 (en) | Two grating lateral shearing wavefront sensor | |
| KR101899026B1 (ko) | 단일 생성 위상 천이 기법을 이용한 디지털 홀로그래픽 복원 장치 및 방법 | |
| Qiu et al. | Linear polarization demosaicking for monochrome and colour polarization focal plane arrays | |
| JP6628103B2 (ja) | デジタルホログラフィ記録装置、デジタルホログラフィ再生装置、デジタルホログラフィ記録方法、およびデジタルホログラフィ再生方法 | |
| KR20180036921A (ko) | 개선된 홀로그래픽 복원 장치 및 방법 | |
| US7924430B2 (en) | Optical heterodyne fourier transform interferometer | |
| Rastogi | Digital optical measurement techniques and applications | |
| CN109313009A (zh) | 用于确定输入射束簇的相位的方法 | |
| US7617060B1 (en) | Extracting higher order information from scene-based Shack-Hartmann wave-front sensing | |
| Fan et al. | Suppressing carrier removal error in the Fourier transform method for interferogram analysis | |
| JP4209023B2 (ja) | 画像計測システム及びその画像校正方法 | |
| Zhou et al. | Speckle-correlation-based non-line-of-sight imaging under white-light illumination | |
| Neuner et al. | Digital adaptive optical imaging for oceanic turbulence mitigation | |
| Xia et al. | Dual-wavelength high-speed digital holographic tomography system for asymmetric air-fluid three-dimensional visualization | |
| Bingleman et al. | DIC-based surface motion correction for ESPI measurements | |
| JPH05306916A (ja) | 縞位相分布解析方法および縞位相分布解析装置 | |
| JP6040469B2 (ja) | デジタルホログラフィ装置 | |
| JP2010002840A (ja) | ディジタルホログラフィ像再生方法及びプログラム | |
| Lekshmi et al. | Fried’s coherence length measurement of dynamic Kolmogorov type turbulence using the autocorrelation function | |
| JP4789799B2 (ja) | 干渉写真における搬送波の発生装置及び方法 | |
| EP4404004A1 (en) | Optical measurement system and optical measurement method | |
| JP3610244B2 (ja) | モアレ測定方法及びそれを用いたモアレ測定装置 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| CF01 | Termination of patent right due to non-payment of annual fee | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20190122 Termination date: 20200323 |