CN104820980B - Adaptive high-precision MTF measurement methods - Google Patents
Adaptive high-precision MTF measurement methods Download PDFInfo
- Publication number
- CN104820980B CN104820980B CN201510178462.7A CN201510178462A CN104820980B CN 104820980 B CN104820980 B CN 104820980B CN 201510178462 A CN201510178462 A CN 201510178462A CN 104820980 B CN104820980 B CN 104820980B
- Authority
- CN
- China
- Prior art keywords
- line
- mtf
- precision
- area
- image
- 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.)
- Active
Links
Landscapes
- Image Processing (AREA)
- Image Analysis (AREA)
Abstract
The present invention provides a kind of adaptive high-precision MTF measurement methods, including steps are as follows:Step 1:Detection image edge pixel point grey scale change;Step 2:The line for building image edge contour supports area;Step 3:Calculate the property parameters that the line based on resize-window supports area;Step 4:Calculate the high-precision MTF based on Kalman filtering.The present invention proposes a kind of adaptive sword side extracting method being based on " gray feature+horizontal line direction ", realize the improvement to recognition status, realize the extracted in self-adaptive of sword side target, the present invention proposes a kind of method of the line spread function optimal value estimation based on Kalman filtering simultaneously, realizes the calculating of high-precision MTF.
Description
Technical field
The invention belongs to aerospace optical remote sensing technical fields, more particularly to a kind of adaptive high-precision MTF measurement methods.
Background technology
Modulation transfer function (Modulation Transfer Function, be abbreviated as MTF) is the one of remote sensor performance
A comprehensive evaluation index, value directly reflect image different frequency attenuation.It can based on quality of image degeneration theory
Know, if the MTF of imaging system can accurately measure to obtain, real image can be by high quality resume.In view of shadow
The MTF of picture can be by atmospheric environment when satellite platform unstability, satellite load aging, image signal-to-noise ratio and satellite imagery etc.
The influence of factor, therefore the mtf value of satellite imaging equipment is one with the changed function of shooting time and camera site,
In order to obtain the accurate MTF of image, most simple most intuitive method is to calculate mtf value one by one to different images.At this stage
MTF extracting methods rely primarily on the ground reference artificially chosen, and the data volume drastically expanded makes the place of this choosing method
Reason efficiency faces the challenge, while artificially choosing ground reference can not make the high-precision of remote sensing image restore realization completely certainly
Dynamicization processing.Therefore it needs to consider the method that adaptive ground reference is chosen, improves the efficiency of processing.On the other hand, due to
The distribution of MTF is extremely sensitive to noise, how in the presence of noise, the optimal estimation of image MTF is obtained, to be promoted
The recovery accuracy of image also determines the key factor of the quality of image after restoring.Therefore, adaptive ground reference is studied
Object extracts and high-precision MTF measurement methods, for improving the precision of image restoration, reducing labor intensity, has important
Theoretical and practical significance.
MTF measurement methods based on remote sensing images, common method is recognition status and impulse method both at home and abroad at present, due to sword
Side target is easily chosen, energy foot, therefore uses in actual production more frequent.However as previously mentioned, current sword
Side method is larger to the dependence for artificially choosing target, cannot be satisfied the Time Inconsistency and space inconsistency of MTF estimations, together
When due to current recognition status it is insufficient to the rejection ability of noise, once to the estimation of MTF, there may be deviations, are based on this, need
Processing is optimized to current recognition status.
The present invention is namely based on the optimization of recognition status, wherein optimization includes mainly two aspects:First is sword side target
It is automatic to choose;Second is the optimal estimation of line spread function.
Invention content
The technical problem to be solved in the present invention is:A kind of measurement sides adaptive high-precision MTF for image procossing are provided
Method optimizes existing recognition status, improves treatment effeciency and graphical quality to image.
The technical scheme is that:
A kind of adaptive high-precision MTF measurement methods, including steps are as follows:Step 1:Detection image edge pixel point gray scale
Variation;Step 2:The line for building image edge contour supports area;Step 3:Calculate the attribute that the line based on resize-window supports area
Parameter;Step 4:Calculate the high-precision MTF based on Kalman filtering.
Further:Detection image edge pixel point grey scale change in step 1, including:Step 11:Calculate each pixel
Gradient magnitude and horizontal line direction;Step 12:Each pixel of image is encoded according to coding rule.
Further, the line that image edge contour is built in step 2 supports area, including:Step 21:It is neighbouring to geometric position
And divide same line of the encoded radio variation less than 1 supports area;Step 22:Area is supported to optimize place on the line that upper step divides
Reason, excluding gross error point.
Further, the property parameters that the line based on resize-window supports area are calculated in step 3, including, step 31:It calculates
Line supports the inclined direction of fitting a straight line in area;Step 32:Image in rectangle frame is calculated using an expansible rectangle frame
Peak factor value KU and crooked angle value SK;Step 33:Reject the line support area for being unsatisfactory for peak value degree and skewness requirement.
Further, the high-precision MTF based on Kalman filtering is calculated in step 4, including:Step 41:Line spread function
Extraction;Step 42:Line spread function Gauss curve fitting;Step 43:High-precision wire spread function based on Kalman filtering calculates;Step
Rapid 43:MTF is calculated and optimization.
The advantages of the present invention over the prior art are that:
The present invention proposes a kind of adaptive sword side extracting method being based on " gray feature+horizontal line direction ", Yi Zhongji
In the method that the best sword side target of resize-window is chosen automatically, a kind of optimal line spread function estimation based on Kalman filtering
Method realizes the improvement to recognition status, realizes the automatic optimal estimation chosen with line spread function of sword side target.
Description of the drawings
Fig. 1 is the flow diagram of the adaptive high-precision MTF measurement methods of the present invention;
Fig. 2 is the MTF curve calculation process schematic diagram of the present invention;
Fig. 3 is the high-precision MTF calculation process schematic diagrames based on Kalman filtering of the present invention.
Specific implementation mode
Since satellite image type of ground objects is various, complex distribution, it is larger that land-based target target chooses difficulty automatically, therefore this hair
The needs of the characteristics of bright combination MTF is calculated and ground automation production, it is little to construct a kind of calculation amount, can be adaptively
Detect continuous marginal point, and edge positioning can reach the Linear feature extraction algorithm of sub-pixel.In order to meet MTF
The intensity profile of the linear edge both sides of the ground reference of calculating should be more uniform feature, introduce KURTOSIS values (peaks
Value coefficient value) determine ground reference candidate region, then according to the intensity profile of the atural object of candidate region linear edge both sides
Situation finally determines whether to calculate desired ground reference to meet MTF.In terms of line spread function optimization, the present invention uses
Gaussian function fitting line spread function, while using the optimal estimation of Kalman filtering estimation line spread function, it is final to calculate life
At MTF curve.
As shown in Figure 1, the realization process of the present invention is:
Step 1:Image border wire-frame image vegetarian refreshments grey scale change detection:
Edge contour is the set that all gray levels change region the most violent in image, therefore the grey scale change of image
Information, i.e. Grad are the key factors of edge extracting, and the effect of pixel gray level variation detection is exactly to calculate in image one by one
The gradient magnitude of pixel, by setting threshold value appropriate, the smaller value of rejecting grey scale change retains the larger value of grey scale change,
Be typically chosen radiation quantizating index 30% to 50% is used as Grads threshold.Pixel gray level variation detection is as follows:
1, gradient magnitude and the horizontal line direction of each pixel are calculated.
Each pixel gradient magnitude and horizontal line direction are calculated according to formula (1) and formula (2):
Wherein,
I (x, y)-be coordinate be (x, y) at gray value, x=1,2,3 ..., N, y=1,2,3 ..., M, N, M are respectively
The length and width of image;
α-coordinate is the horizontal line direction of the pixel i (x, y) at (x, y);
G (x, y)-coordinate is the gradient magnitude of the pixel i (x, y) at (x, y);
2, each pixel of image is encoded according to coding rule,
The Grad calculated according to upper step and horizontal line direction, encode each pixel value, coding rule is as follows:
1) 30% point for being less than quantizating index for gradient magnitude is assigned a value of 0, otherwise is assigned a value of 1.
2) plane space is divided into 8 parts, is -180 ° to -135 ° respectively, -135 ° to -90 °, -90 ° to -45 °, -
45 ° to 0 °, 0 ° to 45 °, 45 ° to 90 °, 90 ° to 135 °, 135 ° to 180 °, each direction scope is divided into 1,2,3,4,5,
6,7,8 coding indicates that the horizontal line of pixel drops into corresponding direction scope just assigns its corresponding volume of range
Code, and the comprehensive coding result of acquisition that the encoded radio is multiplied with the coding result in step (1).
3) calculation code more than 1 point and it is 2 points adjacent between the angle that constitutes, the point for angle equal to 180 ° increases
0.5 encoded radio, with its prominent importance in similar point.
Step 2:Image border contour line supports area's structure:
Line supports that area refers to that geometric position is neighbouring, and the variation of pixel encoded radio is less than 1, and gradient magnitude is more than quantization resolution
The set of the pixel of maximum value 30%.Since atural object edge has continuity, it can be initially believed that these points are to belong to
The point at same edge.Line supports the specific construction method in area as follows:
1, divide same line of neighbouring to the geometric position and encoded radio variation less than 1 supports area;
2, area is supported to optimize processing, excluding gross error point on the line that upper step divides.
When due to coding, the encoded radio for the point that pixel angle is 180 ° is increased 0.5, is encoded to then just will appear
5.5 and 6.5 point is divided into the same line and supports in area, causes subsequent edges positioning to generate error, it is therefore desirable to line branch
It holds area and optimizes processing, weed out rough error point, specific processing method is as follows:
First, linear fit is carried out to all pixels point by the way of least square, fitting formula is as follows:
Yi=Xi·a+b (3)
The form that above-mentioned formula is write as to matrix is as follows:
I.e.:
Y=XA (5)
So non trivial solution is:
A=(XTPX)-1XTPY (6)
Wherein, YiFor the coordinate of pixel in the Y direction, XiIt is pixel in the coordinate of X-direction, i=1,2,3 ..., K, K is
Line supports the number of pixel in area, and a, b respectively represent the slope and intercept of fitting a straight line.
Secondly, it calculates line and supports all the points and fitting a straight line in areaBetween orthogonal distance d, and will be away from
It is rejected as rough error point from the point more than threshold value.
X herein andIndicate the coordinate value of corresponding point in fitting a straight line.
3, it supports area to carry out linear fit on the line of excluding gross error point, is tentatively set to the position where edge line, wherein
Fitting formula is shown in formula (3).
Step 3:Image border line supports area's property parameters to calculate:
After to supporting area to carry out edge line positioning, needs to calculate and support some attribute characteristics in area to determine the branch
It holds in area with the presence or absence of there is the sword edge image for meeting MTF calculating, specifically calculates and judge as follows:
1, the inclined direction that line supports fitting a straight line in area is calculated, retains to tilt and supports area to be used for down for 8 ° to 20 ° of line
The line of one step supports area's attribute to calculate;
2, using an expansible rectangle frame calculate in rectangle frame peak factor value (KURTOSIS) KU of image and
Crooked angle value SK.
The experiment proves that the size of sword edge image for MTF calculating in the presence of influence, the experimental results showed that when scratch figure
Height be less than 15 pixel sizes when, the stability of mtf value is poor, and when scratch figure height be more than 17 pixel sizes when,
MTF calculated values are more stable, and variation is little.Therefore the minimum 17*17 of size for choosing expansible rectangle frame is up to line support
The size in area.
After selected rectangle frame, peak factor value KU and crooked angle value SK are calculated according to the following formula to image in rectangle frame.
Wherein, DiFor the gray value of ith pixel in rectangle frame,For the mean value of image in rectangle frame, n is picture in rectangle
The number of vegetarian refreshments, sd are the standard deviation of all pixels point gray value in rectangle frame.
3, the line support area for being unsatisfactory for peak value degree and skewness requirement is rejected.
Peak value degree is the statistic for describing all value distributional patterns in totality and delaying degree suddenly, KU>3 are known as with transition
Kurtosis.It should be the kurtosis for having transition, i.e. KU that MTF calculating, which requires sword edge image,>3;The degree of bias is similar with kurtosis, it is also description
The statistic of data distribution form describes the symmetry of certain overall value distribution.The degree of bias is close to 0 it is believed that distribution is
Symmetrically.For sword edge image, its distribution should be close to 0.So rejecting KU<3 and the degree of bias>1 line is supported
Area, while the sword edge image met the requirements is obtained for calculating mtf value.
Step 4:High-precision MTF based on Kalman filtering is calculated
Mainly there are four steps for high-precision MTF calculating based on Kalman filtering:Line spread function extracts;Line spread function
Gauss curve fitting;High-precision wire spread function based on Kalman filtering calculates;MTF is calculated and optimization, as shown in Fig. 3.
1, line spread function extracts:
The step of line spread function extracts is as follows:
1) the optimal sword edge image extracted to step 3 is handled, and the position of marginal point is recalculated using Fermi's equation,
And the marginal point extracted is handled using least square, obtain the marginal position of sub-pixel;
2) using marginal position as foundation, resampling is carried out to the row/column data of image, and obtained result is merged, from
And obtain edge spread function;
3) line spread function is obtained to edge spread function derivation.
2, line spread function Gauss curve fitting:
The expression formula of Gaussian curve is as follows:
Wherein, μ is mean value, and σ is standard deviation, and e indicates natural logrithm.
Gaussian curve approximation key is to determine the dispersion of the center and Gaussian Profile of Gaussian Profile.To Gaussian function
Number carries out linearization process, and processing formula is shown below:
Enabling ln (y)=U, x=V, then above formula can be simplified as:
U=a0V2+a1V+a2 (12)
A0, a1, a2 indicate fitting coefficient to be asked respectively herein.At this point,
In formula, ui=ln (yi), yiFor the gray value corresponding to line expansion curve;vi=xi, xiIt is right for line expansion curve
Answer position;N ' is the sampled point number of line expansion curve.
So, the dispersion of the center and Gaussian Profile of Gaussian Profile is calculated using following formula:
Finally, the result of calculating is substituted into formula 10, the line expansion curve y after being fittedi;
3, the high-precision wire spread function based on Kalman filtering generates
Kalman filtering is substantially an optimization autoregression data processing algorithm, and basic principle is according to certain moment
Obtained system mode value and actual measured value is come the true estimation of etching system when fitting this, in of the invention, precise Gaussian
Fitting result is state component, and the line spread function surveyed is actual measured value.
The core of Kalman filtering is five fundamental formulars, it is possible, firstly, to which the process model by system predicts next shape
The system of state, it is k to enable current system mode, then:
X (k | k-1)=AX (k-1 | k-1)+BU (k) (15)
In above formula, A, B indicate coefficient to be asked, X (k | k-1) is using laststate prediction as a result, X (k-1 | k-1)
It is the optimal result of laststate, U (k) is the controlled quentity controlled variable of present status, if not provided, it can be assumed that U (k) is 0.
When updating system model, corresponding system noise variance is also required to do synchronized update, and more new formula is:
P (k | k-1)=AP (k-1 | k-1) A '+Q (16)
In above formula, P (k | k-1) is X (k | k-1) corresponding noise variance of state, and P (k-1 | k-1) is X (k-1 | k-1) shape
The corresponding noise variance of state.A ' is the transposed matrix of A, and Q is the noise variance for including in systematic procedure, due to the band in system
Noise has been removed before recovery, thus herein it is contemplated that noise be white Gaussian noise, wrap in systematic procedure
The noise variance contained may be considered 0.
After the prediction result for obtaining standing state, current state can be obtained in conjunction with the measured value under present status
Optimum estimation value X (k | k):
X (k | k)=X (k | k-1)+Gain (k) (Z (k)-HX (k | k-1)) (17)
Wherein H is state equation, and in conjunction with the computation model of MTF, H herein indicates as follows:
What U was indicated is the secondary reduced equation of Gauss curve fitting, the value for being different sample point of the different subscripts expression of V.
Gain (k) is kalman gain
Wherein, linear equation H ' be H inverse, equation of higher order H ' be H transposed matrix in mathematical computations, in order to obtain
Optimal estimation of the data in function convergence needs to enable Kalman filtering constantly to run and goes down until system convergence, and constantly
The more noise variance of new state k so that noise variance tends to be minimum, and variance more new formula is as follows:
P (k | k)=(I-Gain (k) H) P (k | k-1) (20)
Wherein, I is unit matrix.
4, MTF curve calculates and optimizes
The high-precision wire expansion curve obtained to upper step carries out discrete Fourier transform, takes the mould of each component after transformation to be
The mtf value of each frequency, and on the basis of first mtf value, final MTF curve is obtained as normalized.
Finally, by Frequency point to make normalized on the basis of frequency, then cutoff frequency is 1, is taken at frequency 0~1
Mtf value constitute optical system MTF curve.
The content category ability that description in the present invention is not described in detail is in the known technology of technical staff.
Claims (4)
1. a kind of adaptive high-precision MTF measurement methods, which is characterized in that including steps are as follows:
Step 1:The grey scale change of detection image edge pixel point;
Step 2:The line for building image edge contour supports area;
Step 3:Calculate the property parameters that the line based on resize-window supports area;Including:
Step 31:Calculate the inclined direction that line supports fitting a straight line in area;
Step 32:The peak factor value KU of image and crooked angle value in rectangle frame are calculated using an expansible rectangle frame
SK;
Step 33:Reject the line support area for being unsatisfactory for peak value degree and skewness requirement;
Step 4:Calculate the high-precision MTF based on Kalman filtering.
2. adaptive high-precision MTF measurement methods as described in claim 1, it is characterised in that:Detection figure described in step 1
As edge pixel point grey scale change, including:
Step 11:Calculate gradient magnitude and the horizontal line direction of each pixel;
Step 12:Each pixel of image is encoded according to coding rule.
3. adaptive high-precision MTF measurement methods as described in claim 1, it is characterised in that:Structure figures described in step 2
As the line support area of edge contour, including:
Step 21:Divide same line of and encoded radio neighbouring to the geometric position variation less than 1 supports area;
Step 22:Area is supported to optimize processing, excluding gross error point on the line that step 21 divides.
4. adaptive high-precision MTF measurement methods as described in claim 1, it is characterised in that:Calculating base described in step 4
In the high-precision MTF of Kalman filtering, including:
Step 41:Line spread function extracts;
Step 42:Line spread function Gauss curve fitting;
Step 43:High-precision wire spread function based on Kalman filtering calculates;
Step 43:MTF is calculated and optimization.
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201510178462.7A CN104820980B (en) | 2015-04-15 | 2015-04-15 | Adaptive high-precision MTF measurement methods |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201510178462.7A CN104820980B (en) | 2015-04-15 | 2015-04-15 | Adaptive high-precision MTF measurement methods |
Publications (2)
Publication Number | Publication Date |
---|---|
CN104820980A CN104820980A (en) | 2015-08-05 |
CN104820980B true CN104820980B (en) | 2018-07-24 |
Family
ID=53731265
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201510178462.7A Active CN104820980B (en) | 2015-04-15 | 2015-04-15 | Adaptive high-precision MTF measurement methods |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN104820980B (en) |
Families Citing this family (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN110460765B (en) * | 2018-05-02 | 2020-12-18 | 展讯通信(上海)有限公司 | Method and device for acquiring spatial frequency response curve and electronic equipment |
CN110807768A (en) * | 2019-10-29 | 2020-02-18 | 核工业北京地质研究院 | Remote sensing image quality evaluation method based on MTF |
CN110954303B (en) * | 2019-11-20 | 2021-05-18 | 东南大学 | MTF automatic measurement and calculation method based on high-resolution remote sensing image reference |
CN115544473B (en) * | 2022-09-09 | 2023-11-21 | 苏州吉弘能源科技有限公司 | Photovoltaic power station operation and maintenance terminal login control system |
Citations (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN101793599A (en) * | 2010-03-29 | 2010-08-04 | 中国科学院对地观测与数字地球科学中心 | MTF (Modulation Transfer Function) parameter testing method under condition of nonideal target |
CN101980293A (en) * | 2010-09-02 | 2011-02-23 | 北京航空航天大学 | Method for detecting MTF of hyperspectral remote sensing system based on edge image |
CN102779333A (en) * | 2012-07-10 | 2012-11-14 | 武汉大学 | Optical image restoration method based on Kalman filter |
JP2014203162A (en) * | 2013-04-02 | 2014-10-27 | 日本放送協会 | Inclination angle estimation device, mtf measuring apparatus, inclination angle estimation program and mtf measurement program |
-
2015
- 2015-04-15 CN CN201510178462.7A patent/CN104820980B/en active Active
Patent Citations (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN101793599A (en) * | 2010-03-29 | 2010-08-04 | 中国科学院对地观测与数字地球科学中心 | MTF (Modulation Transfer Function) parameter testing method under condition of nonideal target |
CN101980293A (en) * | 2010-09-02 | 2011-02-23 | 北京航空航天大学 | Method for detecting MTF of hyperspectral remote sensing system based on edge image |
CN102779333A (en) * | 2012-07-10 | 2012-11-14 | 武汉大学 | Optical image restoration method based on Kalman filter |
JP2014203162A (en) * | 2013-04-02 | 2014-10-27 | 日本放送協会 | Inclination angle estimation device, mtf measuring apparatus, inclination angle estimation program and mtf measurement program |
Also Published As
Publication number | Publication date |
---|---|
CN104820980A (en) | 2015-08-05 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN107038717B (en) | A method of 3D point cloud registration error is automatically analyzed based on three-dimensional grid | |
CN104820980B (en) | Adaptive high-precision MTF measurement methods | |
CN106023298B (en) | Point cloud Rigid Registration method based on local Poisson curve reestablishing | |
KR101475382B1 (en) | Method for extracting self adaptive window fourie phase of optical three dimensionl measurement | |
CN110443836A (en) | A kind of point cloud data autoegistration method and device based on plane characteristic | |
CN110084840A (en) | Point cloud registration method, device, server and computer-readable medium | |
CN107481274B (en) | Robust reconstruction method of three-dimensional crop point cloud | |
CN105631939B (en) | A kind of three-dimensional point cloud distortion correction method and its system based on curvature filtering | |
CN103646395B (en) | A kind of High-precision image method for registering based on grid method | |
CN110060342B (en) | Three-dimensional curved surface fitting method | |
CN111524103A (en) | Circular tunnel central axis extraction method based on three-dimensional laser point cloud | |
CN111260668A (en) | Power line extraction method, system and terminal | |
CN109100719A (en) | Combine plotting method with the topographic map of optical image based on satellite-borne SAR image | |
CN115564926A (en) | Three-dimensional patch model construction method based on image building structure learning | |
CN112381862A (en) | Full-automatic registration method and device for CAD (computer-aided design) model and triangular mesh | |
CN113096181A (en) | Method and device for determining pose of equipment, storage medium and electronic device | |
CN114299242A (en) | Method, device and equipment for processing images in high-precision map and storage medium | |
CN108205645A (en) | A kind of reference map quality evaluating method of heterologous image matching system | |
CN113409332A (en) | Building plane segmentation method based on three-dimensional point cloud | |
Zhu et al. | Triangulation of well-defined points as a constraint for reliable image matching | |
CN112288784A (en) | Descriptor neighborhood self-adaptive weak texture remote sensing image registration method | |
CN116756836A (en) | Tunnel super-undermining volume calculation method, electronic equipment and storage medium | |
CN107945218B (en) | Edge large distortion image matching method | |
CN113963114A (en) | Discrete boundary point tracking method based on polygon growth | |
Kang et al. | Point cloud smooth sampling and surface reconstruction based on moving least squares |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
EXSB | Decision made by sipo to initiate substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |