An Image Restoration Method under Complex Optical Imaging Conditions Based on Blind Restoration

Abstract

The present invention discloses an image restoration method under complex optical imaging conditions based on blind restoration, belonging to the field of image restoration. The present invention removes the noise brought by the imaging system based on the morphological filtering method, performs edge detection on the image spectrum of the blurred image under a point light source to remove the blur brought by the imaging system, determines the motion blur parameters based on the cepstrum method, uses the method of flow field data simulation to remove the pneumatic thermal radiation noise, and jointly uses the iterative blind convolution method and the image spectrum reconstruction method for simulation modeling to remove the image blur and image distortion caused by the aero-optical effect. The present invention can denoise, deblur, and correct the distortion of the degraded image under complex imaging conditions, including the preprocessing and image quality evaluation of the image, and combines two methods for flow field modeling. This method can solve the problem of drowning of useful information in the image under complex optical imaging conditions and facilitate the subsequent target detection and positioning of the image.

Description

Technical Field

[0001] The present invention belongs to the field of image restoration, and more specifically, relates to an image restoration method under complex optical imaging conditions based on blind restoration, which is used for image distortion removal, deblurring, and denoising under complex optical imaging conditions. Background Art

[0002] When an imaging detector or a target aircraft flies at high speed in the atmosphere, the flow field near the window becomes extremely complex. Due to its high Mach number and Reynolds number, a shock wave structure will be generated at the head of the aircraft, and phenomena such as a turbulent boundary layer, a shear layer, and a laminar flow field will also be generated near the window, causing thermal radiation interference and image transmission interference to the optical imaging detection system, and ultimately resulting in image offset, jitter, and a decrease in imaging intensity. This phenomenon is also called the aero-optical effect. In the engineering applications of weaponry and equipment in dense atmosphere, due to the existence of the aero-optical effect, in order to pursue imaging quality, the flight speed at the end of the aircraft is often severely restricted, resulting in a serious decline in penetration and strike effects. Therefore, the research on the aero-optical effect and image restoration under the influence of the aero-optical effect has great engineering significance.

[0003] The influence of the aero-optical effect on the image is manifested as vibration, blurring, phase distortion of the optical transmission path, aero-thermal radiation, and noise, etc. The phase distortion is caused by the optical path difference, that is, the plane wave is advanced or delayed by different numbers of wavelengths at different spatial positions, ultimately resulting in the offset of the line of sight or pointing. This part is mainly generated by the laminar part of the flow field. Vibration or jitter stems from the impact effect between the airflow and the aircraft platform, which can be represented by amplitude, frequency, and direction. It will cause the centroid position of the target to jump, resulting in inaccuracies in measurement or positioning. Blurring occurs when the light beam moves relative to the imaging focal plane during a relatively long exposure time interval, causing the recorded image to be blurred, ultimately resulting in attenuation of the detection object intensity, reduction of the detection distance, increase in measurement error, and even detection failure. According to the reasons for its degradation, the restoration of the aero-optical effect is divided into 5 parts, which are, in order, removing imaging system noise, removing imaging system blurring, removing motion blurring, removing aero-thermal radiation noise, and removing aero-optical effect distortion and blurring.

[0004] Currently, the domestic and foreign image restoration methods for degraded images under complex optical imaging conditions mainly use blind restoration methods and non-blind restoration methods, and a small part relies on deep learning and neural network technologies. At present, the main deficiencies of the image restoration method under the aero-optical effect include:

[0005] (1) Since the aero-optical effect will cause multiple types of image degradation, and noise, blurring, and distortion will interact with each other, a single algorithm cannot obtain a restored image with good effects;

[0006] (2) It is difficult to obtain actual flow field data as support and it is impossible to model complex optical imaging conditions;

[0007] (3) The restoration method based on deep learning requires a large amount of data. Since the cost of obtaining measured data is high and the data is not open, it is difficult to build a complete database, making it difficult to apply deep learning methods to restore degraded images under aero-optical effects.

[0008] (4) Existing restoration technologies are sensitive to noise, and the image degradation process will include detector noise and aerodynamic thermal radiation noise, which will reduce the accuracy of the restoration results and have limitations in practical applications. Summary of the invention

[0009] In view of the problems that the restoration of existing aero-optical effect degraded images is easily affected by noise, it is difficult to obtain accurate flow field measured data, and it is impossible to accurately judge the image pixel displacement offset, the present invention proposes an image restoration method under complex optical imaging conditions based on blind restoration to solve the problems of noise, blur and distortion of degraded images under aero-optical effect.

[0010] The basic principle of the image restoration method adopted in the present invention is: according to the mechanism of image degradation, the image degradation process is decomposed into five parts, namely, removing the noise caused by the imaging system, removing the blur caused by the imaging system, removing motion blur, removing aerodynamic thermal radiation noise, and removing aerodynamic optical effect distortion and blur. The degradation principle of each degraded part is studied, and partial algorithms are used in a targeted manner to reconstruct a clear image.

[0011] The method proposed in the present invention for the problem of degraded image restoration under complex optical imaging conditions can remove the noise, blur and distortion of the image under complex imaging conditions, and build a model based on the degradation mechanism to improve the signal-to-noise ratio and clarity of the image, while correcting the pixel displacement distortion in the image imaging process, and finally obtaining a clearer restored target image.

[0012] To achieve the above object, the technical solution adopted by the present invention is:

[0013] The present invention provides a method for restoring a degraded image under complex optical imaging conditions based on blind restoration, comprising the following steps:

[0014] Step 1: Calculate the peak signal-to-noise ratio (PSNR), average sharpness gradient, and coordinates of the maximum pixel value of the original degraded image g0(x,y);

[0015] Step 2: Perform morphological filtering on the degraded image to remove the noise introduced by the imaging system during the imaging process and the noise generated by part of the pneumatic thermal radiation effect, and generate a denoised image g1(x, y); calculate the peak signal-to-noise ratio (PSNR) of the image g1(x, y), and determine whether the increase in this value compared to the PSNR value obtained in Step 1 is greater than 5 dB. If so, proceed to the next step; otherwise, repeat Step 2.

[0016] Step 3: Simulate multiple point light source images passing through the imaging system, and average the multiple point light source images to reduce the influence of noise; perform edge detection on the averaged image after Fourier transform, calculate the estimated values of the blur radius r and standard deviation σ of the imaging system, and through Gaussian function modeling, obtain the imaging system blur point spread function H1(u, v):

[0017]

[0018] Use the imaging system blur point spread function to perform least squares filtering on the image g1(x, y) obtained in Step 2 to obtain an image G2(u, v) that removes the blur of the imaging system

[0019]

[0020] where G1(u, v) is the Fourier spectrum of the image g1(x, y), H1 * (u, v) represents the complex conjugate of H1(u, v), η1 represents the noise level of the image. Perform inverse Fourier transform on G2(u, v), that is, obtain the restored image g2(x, y); calculate the radiation intensity of the pneumatic thermal radiation effect according to the flow field data near the detection window, and subtract this matrix from the image g2(x, y) to obtain a restored image g3(x, y) that removes the pneumatic thermal radiation;

[0021] Step 4: Generate the cepstrum of the restored image g3(x, y) obtained in Step 3, and determine the image motion blur angle and image motion blur scale; synthesize the image motion blur angle and image motion blur scale to obtain the motion blur point spread function H2(u, v); use the motion blur point spread function H2(u, v) to perform least squares filtering on the image g3(x, y) that removes the blur of the imaging system obtained in Step 3 to generate an image G4(u, v) that removes the motion blur:

[0022]

[0023] where \(G_3(u, v)\) is the Fourier transform of the image \(g_3(x, y)\), and \(\eta_2\) represents the noise level of the image; perform the inverse Fourier transform on \(G_4(u, v)\) to obtain the restored image \(g_4(x, y)\); calculate the average gradient of sharpness of the image \(g_4(x, y)\) and check if there is an improvement of more than 5% compared to the value in step 1. If so, proceed to the next step; otherwise, return to step 3 and appropriately increase the values of the noise levels \(\eta_1\) and \(\eta_2\) of the image.

[0024] Step 5: According to the exposure time of the imaging system, divide the degraded images under the complex imaging system into long-exposure degraded images and short-exposure degraded images. The exposure time of the imaging system for long-exposure degraded images is at the second level, and the exposure time of the imaging system for short-exposure degraded images is at the 0.1-second level. Classify the degraded images accordingly. Short-exposure degraded images execute step 6, and long-exposure degraded images execute step 7.

[0025] Step 6: Obtain the spectrum of the image \(g_4(x, y)\) obtained in step 4, perform iterative blind deconvolution on the initial spectrum to estimate the point spread function \(H_3(u, v)\), and at the same time perform spatial domain constraint, support domain restriction, energy constraint, and energy redistribution. After a certain number of iterations, obtain the estimated point spread function \(H_3(u, v)\); use \(H_3(u, v)\) to perform least squares filtering on the image \(g_4(x, y)\) obtained in step 4 to obtain the final restored image \(g_5(x, y)\) under the complex optical imaging conditions; calculate the average gradient of sharpness of the image \(g_5(x, y)\) and the coordinates of the maximum point of the image pixel values. Check if the average gradient of sharpness has an improvement of more than 5% compared to the value in step 4 and if the coordinates of the maximum point of the image pixel values are different from the values in step 1. If all are satisfied, proceed to the next step; otherwise, re-execute step 6 and appropriately increase the value of the noise level \(\eta_3\) of the image.

[0026] Step 7: The mathematical form of the point spread function in the complex imaging environment under long-exposure conditions is where \(\alpha\) and \(\beta\) are parameters to be determined, and \(u\) and \(v\) are frequency coordinates; represent the image degradation process as

[0027] \(O(u, v)=R(u, v)S(u, v)+N(u, v)\)

[0028] where \(O(u, v)\) is the spectrum of the degraded image, \(R(u, v)\) is the point spread function, \(S(u, v)\) is the spectrum of the clear image, and \(N(u, v)\) is the noise; substitute \(G_4(u, v)\) into \(O(u, v)\), take the modulus normalization of each component of the above formula, divide both sides by \(\max(|G_4(u, v)|)\), approximate \(|N'(u, v)| \ll 1\) and take \(u = 0\), and calculate to obtain \(-\alpha|v|\) 2β \(=\ln|G_4'(0, v)|-\ln|S'(0, v)|\); use an isosceles triangle to reconstruct the spectrum of the clear image, and thus calculate \(-\alpha|v|\) 2βThe image is fitted to obtain parameters α and β and a point spread function H4(u, v); the point spread function H4(u, v) is used to perform least squares filtering on the image g4(x, y) to obtain the final clearly restored image g6(x, y) under complex optical imaging conditions; calculate the average gradient of sharpness of the image g6(x, y) and the coordinates of the point with the maximum pixel value of the image. Whether the average gradient of sharpness is increased by more than 5% compared to the value in step 4 and whether the coordinates of the point with the maximum pixel value of the image are different from the value in step 1. If all are met, proceed to the next step; otherwise, re-execute step 7 and at the same time appropriately increase the value of the noise level η4 of the image.

[0029] Step 8: Calculate the comprehensive image quality evaluation index. When the comprehensive image quality parameter is greater than a certain threshold, it is considered that the image restoration effect is good and the image restoration is completed; otherwise, return to step 2, appropriately increase the noise level of the image in the formula, and recalculate.

[0030] Compared with the prior art, the advantages of the present invention are as follows:

[0031] (1) The present invention adopts a method of collecting degraded images under a point light source to calculate the blur parameters of the imaging system, greatly simplifies the process of modeling the blur point spread function, and improves the accuracy of the parameters at the same time;

[0032] (2) The present invention divides the degraded image into long-exposure degraded images and short-exposure degraded images, and proposes different restoration methods for different image types, improving the accuracy of restoration;

[0033] (3) The present invention guides the selection of parameters in the image restoration process through negative feedback of the image quality evaluation index, enabling the image processing and image quality detection to form a closed loop, and the restoration result can be directly applied to subsequent target detection, greatly improving the practicality of the technology;

[0034] (4) The present invention adopts a method of fusing iterative blind deconvolution and image filtering, reduces the influence on the estimated value of the clear image due to the modeling deviation of the point spread function, and improves the accuracy of image restoration;

[0035] (5) The present invention adopts a variety of image quality evaluation indexes to evaluate the blur degree, noise degree and distortion degree of the image, overcoming the defect of only evaluating the image quality in a single dimension in the prior art;

[0036] (6) The present invention adopts a method of decomposing the image degradation process, decomposing the degradation process into noise, imaging system blur, motion blur, blur and distortion based on the aero-optical effect, overcoming the problem of restoring the image with a single model in the prior art. Description of the Drawings

[0037] Figure 1It is the flowchart of the image restoration method under complex optical imaging conditions provided by the embodiments of the present invention.

[0038] Figure 2(a) is the denoising result diagram of the image degraded by uniform pneumatic effect provided by the embodiments of the present invention.

[0039] Figure 2(b) is the denoising result diagram of the image degraded by random pneumatic effect provided by the embodiments of the present invention.

[0040] Figure 3(a) is the result diagram of removing the blur of the imaging system for the image degraded by uniform pneumatic effect provided by the embodiments of the present invention.

[0041] Figure 3(b) is the result diagram of removing the blur of the imaging system for the image degraded by random pneumatic effect provided by the embodiments of the present invention.

[0042] Figure 4(a) is the result diagram of removing motion blur for the image degraded by uniform pneumatic effect provided by the embodiments of the present invention.

[0043] Figure 4(b) is the result diagram of removing motion blur for the image degraded by random pneumatic effect provided by the embodiments of the present invention.

[0044] Figure 5(a) is the result diagram of removing the blur and distortion of pneumatic optics for the short-exposure image degraded by uniform pneumatic effect provided by the embodiments of the present invention.

[0045] Figure 5(b) is the result diagram of removing the blur and distortion of pneumatic optics for the short-exposure image degraded by random pneumatic effect provided by the embodiments of the present invention.

[0046] Figure 6(a) is the result diagram of removing the blur and distortion of pneumatic optics for the long-exposure image degraded by uniform pneumatic effect provided by the embodiments of the present invention.

[0047] Figure 6(b) is the result diagram of removing the blur and distortion of pneumatic optics for the long-exposure image degraded by random pneumatic effect provided by the embodiments of the present invention. Detailed implementation manners

[0048] In order to make the objectives, technical solutions and advantages of the present invention clearer, the present invention will be described in detail below with reference to specific embodiments. Specific embodiments are described below to simplify the present invention. However, it should be recognized that the present invention is not limited to the described embodiments, and various modifications of the present invention are possible without departing from the basic principles, and these equivalent forms also fall within the scope defined by the appended claims of this application.

[0049] The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings in the embodiments of the present invention. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without creative efforts shall fall within the protection scope of the present invention.

[0050] AsFigure 1 As shown in the figure, an image restoration method under complex optical imaging conditions based on blind restoration provided by the present invention mainly includes the following steps:

[0051] Step 1: Calculate the peak signal-to-noise ratio (PSNR), the average gradient of sharpness, and the coordinates of the maximum point of the image pixel values of the original degraded image g0(x, y);

[0052] Step 2: Perform morphological filtering on the degraded image to remove the noise introduced by the imaging system during the imaging process and the noise generated by part of the pneumatic thermal radiation effect, and generate a denoised image g1(x, y); calculate the peak signal-to-noise ratio (PSNR) of the image g1(x, y), and determine whether the improvement amount of this value compared with the PSNR value obtained in Step 1 is greater than 5 dB. If so, proceed to the next step; otherwise, repeat Step 2;

[0053] Step 3: Simulate multiple point light source images passing through the imaging system, and average the multiple point light source images to reduce the influence of noise; perform edge detection after Fourier transform on the averaged image, calculate the estimated values of the blur radius r and the standard deviation σ of the imaging system, and obtain the blur point spread function H1(u, v) of the imaging system through Gaussian function modeling:

[0054]

[0055] Use the blur point spread function of the imaging system to perform least-squares filtering on the image g1(x, y) obtained in Step 2 to obtain an image G2(u, v) that removes the blur of the imaging system

[0056]

[0057] where G1(u, v) is the Fourier spectrum of the image g1(x, y), H1 * (u, v) represents the complex conjugate of H1(u, v), η1 represents the noise level of the image. Perform inverse Fourier transform on G2(u, v) to obtain the restored image g2(x, y); calculate the radiation intensity of the pneumatic thermal radiation effect according to the flow field data near the detection window, and subtract this matrix from the image g2(x, y) to obtain a restored image g3(x, y) that removes the pneumatic thermal radiation;

[0058] Step 4: Generate the cepstrum of the restored image g3(x, y) obtained in Step 3, and determine the image motion blur angle and the image motion blur scale; synthesize the image motion blur angle and the image motion blur scale to obtain a motion blur point spread function H2(u, v); use the motion blur point spread function H2(u, v) to perform least-squares filtering on the image g3(x, y) that removes the blur of the imaging system obtained in Step 3 to generate an image G4(u, v) that removes the motion blur:

[0059]

[0060] where G3(u, v) is the Fourier transform of the image g3(x, y), and η2 represents the noise level of the image; perform the inverse Fourier transform on G4(u, v) to obtain the restored image g4(x, y); calculate the average gradient of sharpness of the image g4(x, y), and check if there is an improvement of more than 5% compared to the value in step 1. If so, proceed to the next step; otherwise, return to step 3 and appropriately increase the values of the noise levels η1 and η2 of the image.

[0061] Step 5: According to the exposure time of the imaging system, divide the degraded images under complex imaging systems into long-exposure degraded images and short-exposure degraded images. The exposure time of the imaging system for long-exposure degraded images is at the second level, and the exposure time of the imaging system for short-exposure degraded images is at the 0.1-second level. Based on this, classify the degraded images. Short-exposure degraded images execute step 6, and long-exposure degraded images execute step 7;

[0062] Step 6: Obtain the spectrum of the image g4(x, y) obtained in step 4, perform iterative blind deconvolution on the initial spectrum to estimate the point spread function H3(u, v), and at the same time perform spatial domain constraint, support domain restriction, energy constraint, and energy redistribution. After a certain number of iterations, obtain the estimated point spread function H3(u, v); use H3(u, v) to perform least squares filtering on the image g4(x, y) obtained in step 4 to obtain the final restored image g5(x, y) under complex optical imaging conditions; calculate the average gradient of sharpness of the image g5(x, y) and the coordinates of the maximum point of the image pixel values. Check if the average gradient of sharpness has an improvement of more than 5% compared to the value in step 4 and if the coordinates of the maximum point of the image pixel values are different from the values in step 1. If all are satisfied, proceed to the next step; otherwise, re-execute step 6 and appropriately increase the value of the noise level η3 of the image.

[0063] Step 7: The mathematical form of the point spread function in a complex imaging environment under long-exposure conditions is where α and β are parameters to be determined, and u and v are frequency coordinates; represent the image degradation process as

[0064] O(u, v) = R(u, v)S(u, v) + N(u, v)

[0065] where O(u, v) is the spectrum of the degraded image, R(u, v) is the point spread function, S(u, v) is the spectrum of the clear image, and N(u, v) is the noise; substitute G4(u, v) into O(u, v), take the modulus normalization of each component of the above formula, divide both sides by max(|G4(u, v)|), approximate |N′(u, v)| << 1 and take u = 0, and calculate to obtain -α|v| 2β= ln|G4′(0, v)| - ln|S′(0, v)|; Reconstruct the clear image spectrum using an isosceles triangle, and thus calculate -α|v| 2β of the image, and obtain parameters α and β and the - point spread function H4(u, v) after fitting; Use the point spread function H4(u, v) to perform least - squares filtering on the image g4(x, y) to obtain the final clear restored image g6(x, y) under complex optical imaging conditions; Calculate the average gradient of sharpness of the image g6(x, y) and the coordinates of the point with the maximum pixel value of the image. Whether the average gradient of sharpness has an improvement of more than 5% compared to the value in step 4 and whether the coordinates of the point with the maximum pixel value are different from the value in step 1. If all are met, proceed to the next step; Otherwise, re - execute step 7, and at the same time appropriately increase the value of the noise level η4 of the image;

[0066] Step 8: Calculate the comprehensive image quality evaluation index. When the comprehensive image quality parameter is greater than a certain threshold, it is considered that the image restoration effect is good, and the image restoration is completed; Otherwise, return to step 2, appropriately increase the noise level of the image in the formula, and recalculate.

[0067] The following introduces the specific implementation methods.

[0068] In a specific embodiment of the present invention, step 1 is: Calculate the image quality evaluation parameters of the original degraded image g0(x, y), and the main steps are:

[0069] 1 - 1: Calculate the peak signal - to - noise ratio PSNR of the original degraded image g0(x, y)

[0070]

[0071]

[0072] where M and N represent the length and width of the image, g(x, y) represents the pixel value of the point (x, y) in the distorted image, g′(x, y) represents the pixel value of the reference image, and L is the maximum value of pixel gray level, generally 255 is used.

[0073] 1 - 2: Calculate the average gradient of sharpness of the original degraded image g0(x, y)

[0074]

[0075] where k x (x, y) and k y (x, y) represent the image gradients calculated by the Sobel operator in the x and y directions at the pixel point (x, y), and l(x, y) represents the image gradient calculated by the Laplacian operator at the pixel point (x, y).

[0076] 1 - 3: Calculate the coordinates (x max1 , y max1 ) of the point with the maximum pixel value in the original degraded image.

[0077] In this embodiment, step 2 is: denoise the original degraded image g0(x, y), and the main steps are as follows:

[0078] 2 - 1: Perform an erosion operation on the original degraded image g0(x, y) to smooth the noise points and make the image noise into blocks of the same level.

[0079] 2 - 2: Perform morphological reconstruction on the eroded image to extract the connected regions with the same features in the original image. Use the eroded image and the original image as the marker and mask respectively to constrain the transformation process and implement the opening operation of the image to obtain the denoised image g1(x, y). Process the uniformly aerodynamically degraded image and the randomly aerodynamically degraded image respectively, and the results are shown in Figures 2(a) and 2(b) respectively. This attached figure is the restored image after denoising the degraded image.

[0080] 2 - 3: Calculate the peak signal - to - noise ratio PSNR of the image g1(x, y), and determine whether the improvement amount of this value compared to the PSNR value obtained in step 1 - 1 is greater than 5 dB. If so, proceed to the next step; otherwise, repeat step 2.

[0081] In this embodiment, step 3 is: regard the point - spread function generated by the imaging system as a Gaussian function, defined as

[0082]

[0083] where σ is the standard deviation of the normal distribution, and (u, v) are the frequency coordinates of the image. The standard deviation σ determines the smoothness of the point - spread function. The larger σ is, the greater the influence of distant pixels on the central pixel, and the smoother the image. Since the point - spread function H1(u, v) of the detector imaging system does not vary according to the imaging scenario, the parameters of the point - spread function H1(u, v) can be determined through experiments. Step 3 specifically includes the following steps:

[0084] 3 - 1: Obtain multiple degraded images after passing through the imaging system under the same point light source. These images contain imaging system noise and blur. First, perform weighted averaging on multiple images to suppress noise, and then perform Fourier transform on the noise - suppressed image to obtain the spectrum of the degraded image.

[0085] 3 - 2: Perform edge detection on the spectrum of the degraded image. Apply the Sobel operator approximation to find the point with the maximum image gradient to locate the edge. Set the edge points to 1 and the non - edge points to 0. After determining the edge points, calculate the radius size of the non - zero region to obtain the parameter r.

[0086] 3-3: Generally, the standard deviation σ ranges from 0 to 8. Use the parameters r and σ to perform Wiener filtering on the image g1(x, y) obtained in step 2-2, calculate the PSNR of the filtered image, and take the standard deviation σ corresponding to the optimal PSNR value as the estimated value of the Gaussian kernel standard deviation σ to obtain the point spread function H1(u, v) of the imaging system.

[0087] 3-4: Perform constrained least squares filtering on the image g1(x, y) obtained in step 2-2 to obtain the image G2(u, v) that removes the blur of the imaging system.

[0088]

[0089] where G1(u, v) is the Fourier spectrum of the image g1(x, y), H1 * (u, v) represents the complex conjugate of H1(u, v), and η1 represents the noise level of the image, which is selected according to empirical values. Perform inverse Fourier transform on G2(u, v) to obtain the restored image g2(x, y).

[0090] 3-5: According to the radiation data near the detection window, calculate the radiation intensity of the aerodynamic heat radiation effect, calculate the radiation intensities of CO2 at 2.7 μm and 4.3 μm and H2O at 1.9 μm, 2.7 μm, and 6.3 μm respectively, and obtain a radiation matrix of the same size as the degraded image after interpolation. Subtract this matrix from the image g2(x, y) obtained in step 3-4 to obtain the restored image g3(x, y) that removes the aerodynamic heat radiation. For the uniformly aerodynamically degraded image, a restored image that completely filters out the aerodynamic heat radiation can be obtained; for the randomly aerodynamically degraded image, since the image has been denoised in step 2, this processing also has a certain filtering effect on the aerodynamic random radiation, so no other processing is performed, and the restored images are shown in Figures 3(a) and 3(b). This attached figure is the restored image after removing the detector blur and the aerodynamic heat radiation from the degraded image.

[0091] In this embodiment, step 4 is: perform motion blur removal processing on the image, and the main steps are as follows:

[0092] 4-1: Since the exposure time of the CCD on the imaging surface or the near-infrared imaging surface is extremely short, it can be approximately considered that there is relative uniform motion between the aircraft and the imaging detector during this exposure period. Therefore, the point spread function H2(u, v) of the motion blur is a straight line.

[0093] 4-2: The cepstrum of the image a(x, y) is defined as A(p, q) = F -1{log|A(u,v)|}, where 1(A,v) is the Fourier transform of a(x,y), (x,y) are the image coordinate points, (u,v) are the image spectrum coordinate points, and (p,q) are the cepstrum coordinate points of the image. Calculate the cepstrum of the image g3(x,y) obtained in step 3-5. After centering the low-frequency components in the cepstrum, perform edge detection. The edge points are assigned a value of 1, and the other points are assigned a value of 0.

[0094] 4-3: Perform a Radon transform on the edge detection result obtained in step 4-2 within the range of 0-180°. The angle corresponding to the maximum value of the Radon curve is the image motion blur angle.

[0095] 4-4: Rotate the image cepstrum obtained in step 4-2 clockwise by an angle equal to the image motion blur angle. Calculate the average value of the pixel values in the column direction of the rotated image. The column number corresponding to the first negative column average is the motion blur scale.

[0096] 4-5: Draw the motion blur function H2(u,v) based on the blur scale and blur angle obtained in steps 4-3 and 4-4; perform constrained least squares filtering on the image g3(x,y) obtained in step 3-5 to obtain the image G4(u,v)

[0097]

[0098] where G3(u,v) is the Fourier transform of the image g3(x,y), and η2 represents the noise level of the image, which is selected based on empirical values. Perform an inverse Fourier transform on G4(u,v) to obtain the restored image g4(x,y). The filtering results are shown in Figures 4(a) and 4(b). This figure is the restored image after removing motion blur.

[0099] 4-6: Calculate the average gradient of the sharpness of the image g4(x,y) obtained in step 4-5. Check if there is an improvement of more than 5% compared to the value in step 1-2. If so, proceed to the next step; otherwise, return to step 3-4, appropriately increase the values of the image noise levels η1 and η2, and recalculate.

[0100] In this embodiment, step 5 is: According to the exposure time of the imaging system, divide the degraded image under the complex imaging system into long-exposure images and short-exposure images. The exposure time of the long-exposure degraded image imaging system is at the second level, and the exposure time of the short-exposure degraded image imaging system is at the 0.01-0.1 second level. Classify the degraded images accordingly. The short-exposure degraded images execute step 6, and the long-exposure degraded images execute step 7.

[0101] In this embodiment, step 6 is as follows: Under short exposure conditions, it is impossible to model the point spread function in a complex imaging environment. Therefore, an iterative blind deconvolution method is used to determine the point spread function H3(u, v).

[0102] 6-1: Obtain the spectrum C(u, v) of the motion-blur-removed image g4(x, y) obtained in steps 4-5. Let Fw(u, v) be the spectrum of the clear image, with the initial value equal to C(u, v), and the initial value of the point spread function be a matrix Hw(u, v) with the same size as the image spectrum, and the initial value in the matrix is 0.

[0103] 6-2: Iteratively calculate the spectrum H(u, v) of the point spread function

[0104]

[0105] where θ1 is an iterative parameter, which is selected to be a smaller value at a higher signal-to-noise ratio for better recovery of image details, and a larger value at a lower signal-to-noise ratio; noiselevel is determined using prior knowledge of the noise level. Since the image has been denoised in step 2, it is generally taken as 0.001 - 0.01.

[0106] 6-3: Perform an inverse Fourier transform on H(u, v) obtained above to get h(x, y), and perform spatial domain constraints on it, that is

[0107] hw(x, y) = h(x, y)

[0108] hw(x, y) = 0 when hw(x, y) < 0

[0109] Perform support region limitation and energy constraint on hw(x, y) obtained above

[0110] hw(x, y) = hw(x, y) * d(x, y)

[0111] hw(x, y) = hw(x, y) / ∑hw(x, y)

[0112] where the size of the support region d(x, y) is the same as the size of the point spread function, and the value is fixed at 100. Perform a Fourier transform on the point spread function hw(x, y) obtained above to get a new estimated value Hw(u, v) of the point spread function after one iteration.

[0113] 6-4: Iteratively calculate the spectrum F(u, v) of the clear image

[0114]

[0115] Where θ2 is an iterative parameter, which is chosen to be a smaller value at higher signal-to-noise ratios for better recovery of image details and a larger value at lower signal-to-noise ratios. The obtained F(u, v) is subjected to inverse Fourier transform to obtain f(x, y), and spatial domain constraints are imposed on it.

[0116] fw(x, y) = f(x, y)

[0117] fw(x, y) = 0 when fw(x, y) < 0

[0118] Energy redistribution is performed on the obtained fw(x, y).

[0119] E = ∑|fw(x, y) - f(x, y)|

[0120]

[0121] Where E represents the energy contained in the image, and m and n are the dimensions of the image g4(x, y) in step 6-1. Based on the above fw(i, j), its Fourier transform can be performed to obtain the estimated value Fw(u, v) of the clear image spectrum after one iteration.

[0122] 6-5: Repeat steps 6-2 to 6-4 for iteration. The number of iterations is chosen to be 300 times to obtain the estimated value H3(u, v) of the final point spread function. Constrained least squares filtering is performed on the image g4(x, y) in step 6-1.

[0123]

[0124] Where G4(u, v) is the Fourier spectrum of g4(x, y), and η3 represents the noise level of the image, which is chosen according to empirical values. Performing inverse Fourier transform on G5(u, v) gives the restored image g5(x, y), as shown in Figures 5(a) and 5(b). This attached figure is the restored image after removing the blur and pixel distortion caused by the aero-optical effect from the degraded image.

[0125] 6-6: Calculate the average gradient of sharpness and the coordinates of the maximum point of the image pixel values of the image g5(x, y) obtained in step 6-5. Whether the average gradient of sharpness has an increase of more than 5% compared to the value in step 4-6 and whether the coordinates of the maximum point of the image pixel values are different from the values in step 1-3. If all are satisfied, proceed to the next step; otherwise, return to step 6-5, appropriately increase the value of the image noise level η3, and recalculate.

[0126] Step 7: Under long exposure conditions, the point spread function in a complex imaging environment has the characteristic of Gaussian distribution, and its mathematical form is

[0127]

[0128] Among them, α and β are parameters to be determined, and u and v are frequency coordinates.

[0129] 7-1: For the image degradation caused by complex imaging conditions such as aero-optics, the degradation process is expressed as the following model:

[0130] O(u, v) = R(u, v)S(u, v) + N(u, v)

[0131] Among them, O(u, v) is the spectrum of the degraded image, R(u, v) is the point spread function, S(u, v) is the spectrum of the clear image, and N(u, v) is the noise. Substitute the image spectrum G4(u, v) obtained in step 4-5 into O(u, v), take the modulus of each component of the above formula for normalization, and divide both sides by max(|G4(u, v)|) to obtain

[0132] |G4′(u, v)| = |R′(u, v)||S′(u, v)| + |N′(u, v)|

[0133] Since the image has been denoised in step 2, it can be approximately considered that |N′(u, v)| << 1, and the above formula can be approximated as |G4′(u, v)| ≈ |R′(u, v)||S′(u, v)|. Substitute And take the logarithm of both ends:

[0134] ln|G4′(u, v)| = -α(u 2 + v 2 ) β + ln|S′(u, v)|

[0135] Since the Gaussian distribution has rotational symmetry, only considering any single line passing through the origin of the Fourier plane can reflect the frequency domain characteristics of the entire image. Select u = 0:

[0136] ln|G4′(0, v)| = -α(v 2 ) β + ln|S′(0, v)|

[0137] -α|v| 2β = ln|G4′(0, v)| - ln|S′(0, v)|

[0138] According to the degraded image spectrum information G4′(0, v) and the clear image spectrum information S′(0, v), the parameters α and β can be calculated.

[0139] 7-2: According to the description in step 7-1, the degraded image spectrum G4′(0, v) is equivalent to the product of the clear image spectrum S′(0, v) and the Gaussian distribution. Since the Gaussian distribution is in the form of high in the middle and low on both sides, and its spectrum is close to 1 at the middle, the spectrum of the central part of the degraded image changes little before and after the multiplication, while the part far from the center will decrease rapidly, and the overall spectrum moves down. At the same time, according to the statistics of the spectral information of a large number of clear images, the best fitting curve of the natural image spectrum has a similar slope in the logarithmic coordinate system, and 80% of it is distributed between -0.8 and -1.5. Combining the above two assumptions, an isosceles triangle is used to reconstruct the clear image spectrum, that is

[0140]

[0141] where the size of the image spectrum is 2M*2N, and a and b represent the slope and intercept of one side of the isosceles triangle, which are parameters to be determined

[0142] 7-3: In the blurred image spectrum G4′(0, v), determine two points A(x1, y1) and B(x2, y2), where

[0143] x1 = n1

[0144] x2 = N - n2

[0145]

[0146]

[0147] where 0 < n1, n2 < 10, and the parameters -1 < ε1 < 0 < ε2 < 1 are used for compensation. Combining the coordinates of two points on the straight line, according to the straight line equation y = ax + b, calculate the parameters a and b, so as to obtain the clear image spectrum ln|S′(0, v)| in step 7-2, and then calculate -α|v| according to the formula in step 7-1 2β curve

[0148] 7-4: The error of the calculated -α|v| 2β will increase with the increase of |v|, and the image of -α|v| 2β will show a decreasing-increasing-decreasing-increasing pattern in the range of -N < v < N; according to the expression form of -α|v| 2β it is concluded that the parameters α and β in the increasing-decreasing range within a certain threshold ε at the origin are correct. Therefore, determine the frequency threshold ε, and perform a second-order fitting on the image and the expression -α|v| 2β in the range of -ε < v < ε to obtain the parameters α and β, and finally calculate the point spread function

[0149] 7-5: Perform constrained least squares filtering on the motion-blur-removed image g4(x, y) obtained in step 4-5

[0150]

[0151] where η4 represents the noise level of the image, which is selected according to empirical values; G4(u, v) is the Fourier transform of the image g4(x, y). Perform inverse Fourier transform on G6(u, v) to obtain the restored image g6(x, y). The results are shown in Figures 6(a) and 6(b). This attached figure is the restored clear image after removing the blur and pixel distortion caused by the aero-optical effect from the degraded image

[0152] 7-6: Calculate the average gradient of sharpness and the coordinates of the maximum point of the image pixel values of the image g6(x, y) obtained in step 7-5. Check if the average gradient of sharpness has an increase of more than 5% compared to the value in step 4-6 and if the coordinates of the maximum point of the image pixel values are different from those in step 1-3. If all conditions are met, proceed to the next step; otherwise, return to step 7-5, appropriately increase the value of the image noise level η4, and recalculate

[0153] In this embodiment, step 8 is: Calculate the comprehensive image quality evaluation parameter. The main steps are as follows

[0154] 8-1: Calculate the SSIM index of the final restored clear images g5(x, y) and g6(x, y)

[0155] SSIM = [L((g(x, y), g′(x, y))] α ·[C((g(x, y), g′(x, y))] β ·[S((g(x, y), g′(x, y))] γ

[0156] where L((g(x, y), g′(x, y)), C((g(x, y), g′(x, y)), and S((g(x, y), g′(x, y)) are the gray-scale similarity, contrast similarity, and structure similarity between the distorted image and the reference image respectively, and α, β, γ are weights, usually taking α = β = γ = 1

[0157] 8-2: Calculate the MSSIM index of the final restored clear images g5(x, y) and g6(x, y)

[0158]

[0159] where M represents the number of times of low-pass filtering and 1 / 2 downsampling of the reference image and the distorted image, usually taking α M = β j = γ j and

[0160] 8-3: Calculate the GSSIM index of the finally restored clear images g5(x, y) and g6(x, y)

[0161] GSSIM = [L((g(x, y), g′(x, y))] α ·[C((g(x, y), g′(x, y))] β ·[G((g(x, y), g′(x, y))] γ where G((g(x, y), g′(x, y)) is the structural similarity based on the gradient component.

[0162] 8-4: Calculate the comprehensive image quality parameter

[0163]

[0164] where τ1, τ2, and τ3 are weights, and generally τ1 = τ2 = τ3 = 1. When the comprehensive image quality parameter is greater than a certain threshold, it is considered that the image restoration effect is good and the image restoration is completed. This threshold generally takes a value greater than 0.95; otherwise, return to step 2, appropriately increase the noise level of the image in the formula, and recalculate.

[0165] The method of the embodiment of the present invention, compared with the prior art, fully considers the mechanism of image degradation, applies multiple algorithms for different degradation methods, and can obtain a better restoration effect than the image restoration using a single algorithm. At the same time, restoration methods are provided for both short-exposure images and long-exposure images, increasing the practicability of the method. The whole process from image preprocessing to final image quality evaluation can be realized, and the accuracy and quality of the restoration meet the requirements.

[0166] The attached drawings description shown in the embodiment of the present invention can make the purpose, technical solution and advantages of the present invention introduced more clearly. It should be noted that the specific embodiments described here are only used to explain the present invention and are not used to limit the present invention. Any equivalent replacement, improvement, etc. made within the method ideas and principles provided by the present invention shall be included within the protection scope of the present invention.

Claims

1. An image restoration method under complex optical imaging conditions based on blind restoration, characterized in that, It includes the following steps: Step 1: Calculate the peak signal-to-noise ratio (PSNR), the average gradient of sharpness, and the coordinates of the maximum point of the image pixel values of the original degraded image \(g_0(x,y)\); Step 2: Perform morphological filtering on the degraded image to remove the noise introduced by the imaging system during the imaging process and the noise generated by part of the aerodynamic heat radiation effect, and generate a denoised image \(g_1(x,y)\); Calculate the PSNR of the image \(g_1(x,y)\), and determine whether the improvement in this value compared to the PSNR value obtained in Step 1 is greater than 5 dB. If so, proceed to the next step; otherwise, repeat Step 2; Step 3: Simulate multiple point light source images passing through the imaging system, and average the multiple point light source images to reduce the influence of noise; Perform edge detection after Fourier transform on the averaged image, calculate the estimated values of the blur radius \(r\) and the standard deviation \(\sigma\) of the imaging system, and obtain the blur point spread function \(H_1(u,v)\) of the imaging system through Gaussian function modeling: Use the blur point spread function of the imaging system to perform least-squares filtering on the image \(g_1(x,y)\) obtained in Step 2 to obtain an image \(G_2(u,v)\) that removes the blur of the imaging system where G1(u, v) is the Fourier spectrum of the image g1(x, y), H1 * (u, v) represents the complex conjugate of H1(u, v), η1 represents the noise level of the image. Performing the inverse Fourier transform on G2(u, v), the restored image g2(x, y) can be obtained. According to the flow field data near the detection window, calculate the radiation intensity of the aerodynamic heat radiation effect. Subtract this matrix from the image g2(x, y) to obtain the restored image g3(x, y) without aerodynamic heat radiation; Step 4: Generate the cepstrum of the restored image \(g_3(x,y)\) obtained in Step 3, and determine the image motion blur angle and the image motion blur scale; Synthesize the image motion blur angle and the image motion blur scale to obtain a motion blur point spread function \(H_2(u,v)\); Use the motion blur point spread function \(H_2(u,v)\) to perform least-squares filtering on the image \(g_3(x,y)\) that removes the blur of the imaging system obtained in Step 3 to generate an image \(G_4(u,v)\) that removes motion blur: where \(G_3(u,v)\) is the Fourier transform of the image \(g_3(x,y)\), and \(\eta_2\) represents the noise level of the image; Perform inverse Fourier transform on \(G_4(u,v)\) to obtain the restored image \(g_4(x,y)\); Calculate the average gradient of sharpness of the image \(g_4(x,y)\), and determine whether there is an improvement of more than 5% compared to the value in Step 1. If so, proceed to the next step; otherwise, return to Step 3, and at the same time appropriately increase the values of the noise levels \(\eta_1\) and \(\eta_2\) of the image; Step 5: According to the exposure time of the imaging system, divide the degraded images under the complex imaging system into long-exposure degraded images and short-exposure degraded images. The exposure time of the imaging system for long-exposure degraded images is at the second level, and the exposure time of the imaging system for short-exposure degraded images is at the 0.1-second level. Based on this, classify the degraded images. The short-exposure degraded images execute Step 6, and the long-exposure degraded images execute Step 7; Step 6: Obtain the spectrum of the image \(g4(x,y)\) obtained in Step 4, perform iterative blind deconvolution on the initial spectrum to estimate the point spread function \(H3(u,v)\), and at the same time perform spatial domain constraint, support domain limitation, energy constraint, and energy redistribution. After a certain number of iterations, obtain the estimated point spread function \(H3(u,v)\); use \(H3(u,v)\) to perform least squares filtering on the image \(g4(x,y)\) obtained in Step 4 to obtain the final restored image \(g5(x,y)\) under complex optical imaging conditions; calculate the average gradient of sharpness of the image \(g5(x,y)\) and the coordinates of the point with the maximum pixel value of the image. Whether the average gradient of sharpness has an increase of more than 5% compared to the value in Step 4 and whether the coordinates of the point with the maximum pixel value are different from the value in Step 1. If all are satisfied, proceed to the next step; otherwise, re - execute Step 6, and at the same time appropriately increase the value of the noise level \(\eta3\) of the image; Step 7: The mathematical form of the point spread function in a complex imaging environment under long exposure conditions is where α and β are parameters to be determined, and u and v are frequency coordinates; the image degradation process is expressed as O(u, v) = R(u, v)S(u, v) + N(u, v) where \(O(u, v)\) is the spectrum of the degraded image, \(R(u, v)\) is the point spread function, \(S(u, v)\) is the spectrum of the clear image, and \(N(u, v)\) is the noise; substituting \(G4(u, v)\) into \(O(u, v)\), taking the modulus normalization of each component of the above formula, dividing both sides by \(\max(|G4(u, v)|)\), approximating \(|N′(u, v)| \ll 1\) and taking \(u = 0\), after calculation, we get \(-α|v|\) 2β = \(\ln|G4′(0, v)|-\ln|S′(0, v)|\); using an isosceles triangle to reconstruct the spectrum of the clear image, thus calculating the image of \(-α|v|\) 2β The image is fitted to obtain the parameters \(α\) and \(β\) and the point spread function \(H4(u, v)\); using the point spread function \(H4(u, v)\) to perform least squares filtering on the image \(g4(x, y)\) to obtain the final clear restored image \(g6(x, y)\) under complex optical imaging conditions; calculating the average gradient of the sharpness of the image \(g6(x, y)\) and the coordinates of the maximum point of the image pixel values, whether the average gradient of the sharpness is increased by more than 5% compared to the value in step 4 and whether the coordinates of the maximum point of the image pixel values are different from the value in step 1, if all are satisfied, then proceed to the next step; otherwise, re - execute step 7, and at the same time appropriately increase the value of the noise level \(η4\) of the image. Step 8: Calculate the comprehensive image quality evaluation index. When the comprehensive image quality parameter is greater than a certain threshold, it is considered that the image restoration effect is good and the image restoration is completed; otherwise, return to Step 2, appropriately increase the noise level of the image in the formula, and recalculate.

2. The image restoration method under complex optical imaging conditions based on blind restoration according to claim 1, wherein In Step 3, the point - source image passing through the imaging system is a blurred image under the point - source, which contains a certain amount of imaging system noise. Therefore, multiple images contaminated by noise and blur under the same point - source are acquired, and the multiple images are first weighted and averaged to suppress noise.

3. The image restoration method under complex optical imaging conditions based on blind restoration according to claim 1, wherein, The specific steps for determining the image motion blur angle in Step 4 are as follows: a) Calculate the cepstrum matrix; b) Center the low - frequency components; c) Perform edge detection, with edge points taking the value of 1 and the remaining points taking the value of 0; d) Perform Radon transform in the range of 0 - 180°; e) Detect the maximum value of the Radon curve, and the angle corresponding to the maximum value is the image motion blur angle.

4. The image restoration method under complex optical imaging conditions based on blind restoration according to claim 1, characterized in that, The specific steps for determining the image motion blur scale in Step 4 are as follows: a) Rotate the image cepstrum clockwise by the image motion blur angle; b) Calculate the sum of pixel values in the column direction of the image; c) Calculate the average value of pixel values in each column; d) With the right - hand direction as the positive x - axis, the column number where the first column average value is negative is the motion blur scale.

5. The image restoration method under complex optical imaging conditions based on blind restoration according to claim 1, wherein, The calculation process of estimating the point spread function \(H3(u,v)\) by the iterative blind deconvolution method in Step 6 is as follows: 6.1 Obtain the spectrum \(C(u,v)\) of \(g4(x,y)\). The spectrum of the to - be - restored clear image is \(Fw(u,v)\), and its initial value is \(C(u,v)\); the to - be - determined point spread function is \(Hw(u,v)\), and the initial value in the matrix is 0; 6.2 Iteratively calculate the spectrum \(H(u,v)\) of the point spread function where \(\theta\) is the iteration parameter and noiselevel is the noise level; 6.3 Perform inverse Fourier transform on the point spread function \(H(u,v)\) to obtain \(h(x,y)\), perform spatial domain constraint, support domain limitation, and energy constraint on it, and after Fourier transform, obtain the estimated value \(Hw(u,v)\) of the point spread function after one complete iteration; 6.4 Iteratively calculate the spectrum \(F(u,v)\) of the clear image where \(\theta2\) is the iteration parameter and noiselevel is the noise level; 6.5 Perform the inverse Fourier transform on the obtained F(u, v) to get f(x, y), perform spatial domain constraint and energy redistribution on it, and after Fourier transform, obtain the estimated value Fw(u, v) of the clear image spectrum after one complete iteration; 6.6 Repeat steps 6.2 to 6.5 for iteration until the set number of iterations is reached to obtain the estimated value H3(u, v) of the final point spread function.

6. The image restoration method under complex optical imaging conditions based on blind restoration according to claim 5, wherein, The steps for performing spatial domain constraint on h(x, y) in step 6 are: hw(x, y) = h(x, y) hw(x, y) = 0 when hw(x, y) < 0 The steps for support region restriction are hw(x, y) = hw(x, y) * d(x, y) where the size of the support region d(i, j) is the same as the size of the PSF; the steps for energy constraint are hw(x, y) = hw(x, y) / ∑hw(x, y).

7. The image restoration method under complex optical imaging conditions based on blind restoration according to claim 5, characterized in that, The steps for performing spatial domain constraint on f(x, y) in step 6 are: fw(x, y) = f(x, y) fw(x, y) = 0 when fw(x, y) < 0 The steps for energy redistribution are: where E represents the energy contained in the image, and m and n are the sizes of the image g4(x, y).

8. The image restoration method under complex optical imaging conditions based on blind restoration according to claim 1, characterized in that The process of reconstructing the clear image spectrum using an isosceles triangle in step 7 is: a) Let where the size of the image spectrum is 2M * 2N, and a and b represent the slope and intercept of one side of the isosceles triangle, which are parameters to be determined; b) In the blurred image spectrum G4′(0, v), determine two points A(x1, y1) and B(x2, y2), where: x1 = n1 x2 = N - n2 where 0 < n1, n2 < 10, -1 < ε1 < 0 < ε2 < 1; combining the coordinates of two points on the straight line, according to the straight line equation y = ax + b, calculate the parameters a and b, so as to obtain the clear image spectrum ln|S′(0, v)|, and calculate the curve of -α|v| 2β curve; c) Quadratically fit the image and the expression -α|v| in the range -ε < v < ε 2β to obtain the parameters α and β, and finally calculate the point spread function 9. The image restoration method under complex optical imaging conditions based on blind restoration according to claim 1, wherein, The definition of the comprehensive image quality parameter in step 8 is where τ1, τ2, and τ3 are weights, and take τ1 = τ2 = τ3 = 1; SSIM is defined as: SSIM = [L((g(x,y), g′(x,y))] α ·[C((g(x,y), g′(x,y))] β ·[S((g(x,y), g′(x,y))] γ where L((g(x,y), g′(x,y)), C((g(x,y), g′(x,y)), and S((g(x,y), g′(x,y)) are the luminance similarity, contrast similarity, and structure similarity between the reference image and the distorted image, respectively, and α, β, γ are weights with α = β = γ = 1; MSSIM is defined as: where M represents the number of times of low-pass filtering and 1 / 2 downsampling on the reference image and the distorted image, and taking α M = β j = γ j and GSSIM is defined as: GSSIM = [L((g(x,y), g′(x,y))] α ·[C((g(x,y), g′(x,y))] β ·[G((g(x,y), g′(x,y))] γ where G((g(x, y), g′(x, y)) is the structural similarity of the reference image and the distorted image based on the gradient component.

Images