CN103136452A - Quick direct demodulation method applicable to spherical data - Google Patents

Abstract

The invention discloses a spherical direct demodulation method which includes the following steps: 1 determining size of a spherical window according to view size of a collimator, determining number of sampling points along any side in the spherical window to generate a discretization gridding of a reference window, utilizing the discretization gridding with identical size and shape with the reference window to cover the whole spherical face, utilizing the orthogonal pixel mode to conduct discretization on the spherical face and enabling the spherical face to be scattered into the gridding with equal-distance adjacent pixels, consistent pixel shape and pixel distribution orthogonality; 2 utilizing a clustering algorithm to group different modulation functions of azimuth angle alpha and dividing the modulation kernel function with the value of the azimuth angle located in the clustering range in one group; and 3 utilizing digital approximation based on clustering to achieve acceleration value calculation in the modulation algorithm.

Description

A kind of fast direct demodulation method that is applicable to sphere data

Technical field

The present invention relates to technical field of data processing, relate to particularly a kind of fast direct demodulation method that is applicable to sphere data.

Background technology

Optics collimator (hereinafter to be referred as collimating apparatus) is a kind of device that only allows incident angle to pass through with interior photon beam in certain limit.Therefore, when using collimating apparatus that light source is scanned, recording collimating apparatus, to pass through the photon line of collimating apparatus when difference is pointed to strong, can observation light source at the stream of collimating apparatus this locality by force in the distribution of diverse location, different directions.

The point spread function of collimating apparatus is to describe optical system to the function of the parsing of pointolite (hereinafter to be referred as point source).Desirable collimating apparatus only allows to pass through along the photon beam of its pointing direction, so in scanning process, the response to point source is still a point, and collimating apparatus that namely and if only if has photon beam to pass through during over against this point source, does not have photon beam to pass through when pointing to other directions.At this moment the point spread function of this collimating apparatus is the δ function of its sensing.

Actual collimating apparatus always allows the photon beam of incident angle in certain successive range to pass through, just generally maximum along the strong density of photon line of collimating apparatus pointing direction, the strong density of photon line of the incident angle scope edge that permission is passed through generally levels off to zero, so in scanning process, collimating apparatus is generally a therefrom hot spot of mind-set edge disperse to the response of point source.This is the origin of point spread function title just.

For example, Fig. 1 is the point spread function of a square pass collimating apparatus.X axle, y axle be square pass collimating apparatus respectively along square hole minor face and the inswept pointing direction of long side direction, unit is radian; Centre coordinate x=0, y=0, the expression collimating apparatus is pointed to the point source over against unit brightness, and this position is the reference position, and the coordinate of other positions i.e. the difference of the pointing direction of reference position therewith.Usually represent as shown in Figure 1 two-dimensional points spread function with P (x, y).

The point spread function of imaging system (as telescope etc.) is described with two-dimentional scalar-valued function traditionally.Yet, in fact the mode of using collimating apparatus scanning is to target (comprising source and on every side background) when observing, especially when the point spread function of collimating apparatus does not have rotational symmetry, observation data except with collimating apparatus point to relevant, also relevant with respect to the rotation (describing with the position angle) of self optical axis with collimating apparatus.

With Modulation equation, data that collimating apparatus scanning observation obtains and the relation between observed object used described.In more general signal system, Modulation equation be used for descriptive system input and output between relation.

Modulation equation comprises: input data, modulation kernel function, output data.In order to describe collimating apparatus scanning observation process, order input data are that the spatial brightness of observed object distributes, the output data are observation data, namely, given scanning position (collimating apparatus sensing, position angle) is strong (for example by the stream of collimating apparatus, can use the detector recording light quantum count such as photomultiplier, read simultaneously the time shutter to this position).As:

d(x，y，α)＝∫Ωp(x，y，α，x’，y’)f(x’，y’)dx’dy’+n(x，y，α)，(1)

Wherein x, y and x ', y ' are the rectangular coordinate of target on the observer plane, α is the position angle that relative its initial position of collimating apparatus rotates on optical axis, 0≤π＜α, f (x ', y ') be the Luminance Distribution of observed object on objective plane Ω, d (x, y, α) is the point on collimating apparatus sensing objective plane And the observation data the when position angle is α, that is strong by the stream of collimating apparatus, n (x, y, α) is noise item, p (x, y, α, x ', y ') is so-called modulation kernel function.The target brightness distribution function also is called image, is defined on objective plane, and coordinate is vector

The observation data function definition is at the collimating apparatus state space, and coordinate is vector

The modulation kernel function namely represents to be positioned at objective plane

Place's unit strength point source collimation device state space

The contribution of observation data.

More general, Modulation equation (1) has been described input signal (occurring with forms such as target brightness distribution function, images) through system modulation functional core FUNCTION MODULATION, with the process of the formal output of observation data.

The demodulation problem, that is, the output data in given Modulation equation (1) and modulation kernel function are found the solution the problem of input signal.The integral equation problem of finding the solution shape such as equation (1) on mathematics is called Fredholm integral equation of the first kind.If the modulation kernel function of given input signal and system, the problem of finding the solution the output data is " direct problem ", and demodulation problem so is exactly a kind of indirect problem.Direct problem is corresponding with observation process, and indirect problem is processed with observation data, the target reduction process is corresponding.

Shown in Modulation equation (1), because the observation process of reality inevitably is subject to noise, indirect problem usually can not guarantee stability and the convergence of uniqueness of solution and numerical solution process.Therefore indirect problem is one of common type of numerical evaluation field ill-conditioning problem.

The direct demodulation method that the cautious a word used in place name of Chinese scholar Lee, Wu Mei propose provides framework flexibly for the numerical solution of demodulation problem.Directly demodulation method is by the iterative computation input signal, and the background constraint that the method is used is guaranteed the numerical procedure convergence with upper limit constraint; Per step iteration is avoided error with the iterative process accumulation and amplifies all from Modulation equation; The initial modulation equation is done integral transformation, improve the stability of numerical procedure.When using direct demodulation method, alternative manner commonly used is the RL iteration.

Although in general directly demodulation method is effectively, when the Spatial Dimension at the higher while objective function of observation data sampling rate place is higher, be again inapplicable under existing design conditions.This be due to, per step of RL iteration all must numerical evaluation Modulation equation (1) the right.The discrete form of Modulation equation (1) is:

$d_{i,j,k}={\mathrm{\Σ}}_{i^{\′}=1,j^{\′}=1}^{M,N}{p}_{i,j,k,i^{\′},j^{\′}}{f}_{i^{\′},j^{\′}}+n_{i,j,k}---\left(2\right)$

Wherein (Δ _x, Δ _yAnd Δ _αThe discrete sampling interval on respective dimensions).

Might as well suppose that at x, y and three dimension up-sampling rates of α be all N, at this moment the modulation meter calculator has the complexity of O (N5).Such complexity will be brought the time overhead that is difficult to bear in solving practical problems.

Yet, if point spread function has rotational symmetry, the value and the α that namely modulate kernel function are irrelevant, point spread function due to collimating apparatus has translation invariance again, modulation is calculated can be reduced to convolutional calculation, and then can accelerate with quick Fu's rank transformation, the reduced complexity of at this moment calculating is to O (N2log (N)).

In some occasions, in the whole day scope sky patrol that need to carry out as cosmology, astrophysics and uranology field, objective function is defined on sphere (as celestial sphere) coordinate system.When observation data has higher sampling rate, utilize existing design conditions just to be difficult to use direct demodulation method to process the data of modulation type recording geometry output.

The below explains respectively from two aspects.

1) sphere is different from the Topological property on plane.The properties that rectangular coordinate system on the plane possesses makes it become the coordinate system that numerical evaluation is carried out in natural being convenient to.For example: uniformly-spaced after discretize, all sampled point areas equate, shape is identical, can travel through all sampled points along orthogonal any both direction (being called row and row) respectively along the both direction of quadrature, etc.The Topological property of sphere makes can't construct the discretize grid with aforesaid properties on sphere.Therefore, the acceleration strategy that is applicable to the plane is not directly applied for sphere.

2) the modulation kernel function of recording geometry not necessarily has translation invariance and rotational symmetry.Above-mentioned problem can be evaded by global calculation being divided into a series of local calculation.Even if but on certain partial sphere, even this curved surface can be approximately the plane with any high precision, still need the modulation kernel function of recording geometry to possess some characteristics, can realize the acceleration strategy that chapters and sections are mentioned.That is, translation invariance and rotational symmetry.So-called translation invariance refers to:

p(x，y，α，x’，y’)＝p(x+r _x，y+r _y，α，x’+r _x，y’+r _y)

＝p(x-x’，y-y’，α)， (3)

So-called rotational symmetry refers to:

p(x，y，α，x’，y’)＝p(x，y，α’，x’，y’)

＝p(x，y，x’，y’)， (4)

Still the square pass collimating apparatus of introducing in the literary composition is as example, and this collimating apparatus does not just have rotational symmetry.Do not have equally honeycomb type collimating apparatus that also having of rotational symmetry use in real work etc.

To sum up, the prior art modulation that also do not have at present a kind of method can accelerate on sphere is calculated.Existing unique similar techniques only can be accelerated the direct demodulation of the translation invariant Modulation equation on the plane and calculate.This acceleration is to replace translation invariant modulation (in fact the definition of convolution is exactly translation invariant modulation) with convolution, replaces general modulation matrix with the Toeplitz matrix, thereby utilizes fast fourier transform to accelerate.The shortcoming of this technology is obvious, i.e. narrow application range.

Summary of the invention

In order to overcome defects of the prior art, the present invention proposes the direct demodulation method of a kind of sphere.

The direct demodulation method of sphere of the present invention comprises step: step 1, according to collimating apparatus visual field size, determine the sphere window size and determine the interior quantity along any limit sampled point of sphere window, the discretize grid of generating reference window, use the discretize grid identical with the reference windows size, that shape is identical to cover whole sphere, use accurate quadrature pixelation mode to carry out discretize to sphere, with sphere discretely turn to that neighbor is equidistant, primitive shape is consistent, the grid of pixel distribution quadrature; Step 2 uses clustering algorithm that parameter alpha different modulating function in position angle is divided into groups, and the value of position angle parameter alpha is divided into one group in the cluster scope with interior modulation kernel function; Step 3 is used the acceleration numerical evaluation that realizes modulation operation based on the numerical value approximation computation of cluster.

Preferably, step 2 further comprises: step 210 is decomposed into the modulation kernel function form that is rotated conversion by the kernel function that does not contain azimuthal variations, thereby the modulation operation of azimuthal variation is decomposed into modulation operation and the twiddle operation at constant bearing angle; Step 220 according to the requirement to demodulation computing degree of accuracy, is determined the radius of above-mentioned cluster; Step 230 according to given cluster radius, is determined each cluster centre and cluster member, makes the position angle of each member in given cluster and the position angle of cluster centre differ in cluster radius; Step 240 is carried out numerical value with step 210 to the rotation of each kernel function and is approached, and namely is approximately the rotation to cluster centre.

Preferably, step 210 further comprises: step 211, and with former modulation kernel function

D (x, y, α)=∫ Ω p (x, y, α, x ', y ') f (x ', y ') dx ' dy '+n (x, y, α) is decomposed into by the position angle:

$p(x,y,\mathrm{\α},x^{\′},y^{\′})={\∫}_0^{2\mathrm{\π}}\mathrm{\Ω}\·p(x,y,\mathrm{\α},x^{\′},y^{\′})f(x^{\′},y^{\′}){\mathrm{dx}}^{\′}{\mathrm{dy}}^{\′}+n(x,y,\mathrm{\α}),,$

Wherein, p is the modulation kernel function, and x, y, x ' and y ' are respectively the azimuthal coordinates variablees, and α and α ' are azimuthal variations, and δ is the finite point impulse function, and Ω refers to azimuthal span; Step 212 according to the decomposition result of step 211, with former modulation Kernel Function Transformation is:

$d(x,y,\mathrm{\α})={\∫}_{\mathrm{\Ω}}{\∫}_0^{2\mathrm{\π}}p(x,y,{\mathrm{\α}}^{\′},x^{\′},y^{\′})f(x^{\′},y^{\′})\mathrm{\δ}(\mathrm{\α}-{\mathrm{\α}}^{\′}){\mathrm{dx}}^{\′}{\mathrm{dy}}^{\′}{\mathrm{d\α}}^{\′}+n(x,y,\mathrm{\α})$

$={\∫}_0^{2\mathrm{\π}}(p\left({\mathrm{\α}}^{\′}\right)*f)(x,y)\mathrm{\δ}(\mathrm{\α}-{\mathrm{\α}}^{\′}){\mathrm{d\α}}^{\′}+n(x,y,\mathrm{\α});$ Wherein, d is the output data, and f is the input data, and n is noise item.

Preferably, step 220 further comprises: described cluster is spherical cluster, and wherein cluster radius c passes through

Determine, wherein ε is the given positional precision ε of user, and the unit of ε is radian, and D is farthest two dot spacings in the point spread function field of definition, and unit is also radian.

Preferably, step 230 further comprises: step 231, all azimuthal variations are carried out quicksort, and it is arranged as: { α according to the numerical values recited ascending order ₀, α ₁, α ₂α _K-1, wherein K is azimuthal quantity; Step 232 makes position angle index k=0, i.e. α _k=α ₀, cluster index l=0, cluster centre Cluster member's number K ₀=1, cluster ordinal number g ₀If=0:a) k＜K carry out next step, otherwise finish to return; B) order Carry out next step; C) if α _k+1-α _k＜c makes g _k+1=l, K _l=K _l+ 1; Otherwise, order

L=l+1,

Make k=k+1, and return to (a), by single ergodic, obtain and azimuthal variations packet marking array { the cluster centre array of gk} and each grouping one to one

This result that obtains is used for numerical value and approaches modulation calculating, through the modulation meter that approaches at last according to K position angle rotation modulation kernel function one by one, after obtaining the above results, the modulation kernel function that belongs to a cluster together the angle according to its cluster centre is rotated.

Preferably, step 240 further comprises: step 241, and in this step, the is olation of discrete form is:

$d_{i,j,k}=\underset{i^{\′}=1,j^{\′}=1}{\overset{M,N}{\mathrm{\Σ}}}{k}^{\′}{P}_{i,j,k^{\′},i^{\′},j^{\′}}{\mathrm{\δ}}_{k-k^{\′}}+n_{i,j,k}$

$={\mathrm{\Σ}}_{k^{\′}=1}{(p_{k^{\′}}*f)}_{i,j}{\mathrm{\δ}}_{k-k^{\′}}+n_{i,j,k},k^{\′}=\mathrm{0,1,2},...,K,---\left(7\right)$

Wherein k is the position angle index, and K is the position angle number;

Step 242, order

$d_{i,j,k}\≈\underset{l}{\mathrm{\Σ}}{(p_l*f)}_{i,j}{\mathrm{\δ}}_{|{\mathrm{\α}}_k-{\stackrel{\·}{\mathrm{\α}}}_l|>c}+n_{i,j,k},l=\mathrm{0,1,2},...,L,---\left(8\right)$

Wherein

For the Kronecker symbol of promoting, work as subscript

Be true time as the logical expression value, the value of this symbol is 1, otherwise the value of this symbol is 0, α _kBe azimuthal k the sampled point of discretize,

α _kThe cluster centre value of place cluster, l is the cluster index, and L is total cluster number, and c is cluster radius.

Preferably, step 3 further comprises: the modulation kernel function position angle of same cluster is consistent in given accuracy, according to the cyclic convolution theorem, realizes with fast fourier transform.

Utilize method of the present invention, the modulation that can accelerate on sphere is calculated, and can realize that reconstruction obtains the higher-dimension image to the low-dimensional observation data.

Description of drawings

Fig. 1 is the two-dimensional points spread function of optics collimator in prior art;

Fig. 2 is the fast direct demodulation method process flow diagram that the present invention is applicable to sphere data;

Fig. 3 is that the method according to this invention is carried out the part first meridian window schematic diagram that the sphere discretize obtains;

Fig. 4 is that the method according to this invention is carried out the two poles of the earth that the sphere discretize obtains and the window schematic diagram of equatorial zone;

Fig. 5 is that the method according to this invention is carried out the Partial Window schematic diagram along given latitude that the sphere discretize obtains;

Fig. 6 is that the method according to this invention is carried out the hemispherical discretize grid covering schematic diagram that the sphere discretize obtains;

Fig. 7 is that the discretize grid that the method according to this invention is carried out the sphere that the sphere discretize obtains covers schematic diagram;

Fig. 8 is that the method according to this invention is based on the discretization scheme schematic diagram of dimetric projection.

Embodiment

For making the purpose, technical solutions and advantages of the present invention clearer, below in conjunction with specific embodiment, and with reference to accompanying drawing, the present invention is described in more detail.

The present invention proposes the direct demodulation method of a kind of sphere, can be used for the reconstruction of low-dimensional observation data and obtains the higher-dimension image, investigates over the ground etc. such as CT tomoscan or remote sensing satellite.

Relate to two special difficult points for the computing on sphere, the one, can not be with the pixel of identical shaped, formed objects with any sampling interval discretize; The 2nd, to not having the function of rotational symmetry, can't realize the translation to the everywhere continuous of this function on sphere, that is, before and after mobile, the function value is not only relevant with mobile front and back Anchor displacement, also the orientation angles of function is relevant, and can't keep orientation angles constant and make this function cover sphere.

Singularity based on the sphere computing, the ultimate principle of the direct demodulation method of sphere that the present invention proposes mainly comprises following three steps: step 1, use accurate quadrature pixelation method to carry out discretize to sphere, with sphere discretely turn to that neighbor is approximate equidistantly, the grid of approximate consistent, the pixel distribution nearly orthogonal of primitive shape; Step 2 uses clustering algorithm that the position angle different modulating function of parameter a is divided into groups, and the modulation kernel function that the value of position angle parameter a is approaching is divided into one group; Step 3 is used the acceleration numerical evaluation that realizes modulation operation based on the approximate treatment of cluster.

Wherein the target of step 1 is to set up nearly orthogonal, uniform pixel grid, realizes the discretize of sphere.The target of step 2 is that orientation angles parameter a is separated from the modulation kernel function, and the movement of modulation kernel function on sphere is decomposed into continuous local translation and rotation.Translation can be converted into convolution, and rotary manipulation can be independent.After step 3 occurs in step 1,2, can calculate a small amount of rotation on the basis of step 1 and 2, then the consistent orientation angles a of pairing approximation takes to be similar to.

Fig. 2 is the process flow diagram of the direct demodulation method of sphere of the present invention, and with reference to Fig. 2, the method comprising the steps of:

Step 1, the sphere discretize.

In this step, use square partial sphere projecting method to carry out the discretize of sphere, according to collimating apparatus visual field size, determine the sphere window size and determine the interior quantity along any limit sampled point of sphere window, the discretize grid of generating reference window uses the discretize grid identical with the reference windows size, that shape is identical to cover whole sphere; Use accurate quadrature pixelation mode to carry out discretize to sphere, with sphere discretely turn to that neighbor is equidistant, primitive shape is consistent, the grid of pixel distribution quadrature, this step further comprises:

Step 101 is selected sphere window size and sampling rate.

In this step, according to collimating apparatus visual field size, select the sphere window size, as, σ * σ, wherein σ is the arc length on unit sphere.Limit the even number that σ only gets π/one, for example: Deng.Sampling rate is the interior quantity along any limit sampled point of sphere window, as N.

Step 102, the discretize grid of generating reference window.

A) intersection point that is centered close to unit sphere and x axle of definition reference windows, that is, (1,0,0), the local ideal plane rectangular coordinate system of window x axle is parallel to global coordinate system xOy plane, and the y axle is parallel to global coordinate system z axle.

B) get and be parallel to yOz plane, the plane x=1 tangent with unit sphere, choose the length of side be on this plane

Square S, its center also is positioned at (1,0,0).

C) with the uniform discrete N that turns to of the area in S ²The little square that individual area equates, shape is identical, each little foursquare center is the sampled point at this place, claims that these little squares are pixel.The position of pixel namely is positioned at the wherein position of the sampled point of the heart.

D) with the pixel on the plane along unit sphere radially projecting to sphere, namely obtain the pixel grid on sphere, as shown in Figure 8.At this moment the elemental area approximately equal on sphere, shape approximation are the sphere square.So-called sphere square, namely the arc between adjacent 2 is orthodrome, and arc length equates, and the arc that intersects is orthogonal at intersection point place tangent line.

Step 103 uses the discretize grid identical with the reference windows size, that shape is identical to cover whole sphere.

A) hypothesis M is any positive integer, at first along the semi arch of x＞0 one side spheres and zOx Plane intersects, take σ as the interval, respectively to two ends rotary reference window, obtain 2M+1 and reference windows shape with big or small all unanimously, the different window (comprising reference windows itself) in position only, be called first meridian window (if specialize abstract unit ball discussed herein as an example of the earth example, get again (0,0,1) be arctic point, (1,0,0) locating longitude is 0, and the window that herein produces is all along the first meridian).As shown in Figure 3.The position of regulation window is with the positional representation at its center.

B) along " equator ", take σ as the interval, rotation obtains 4M equator window, as shown in Figure 4.

C) under the line and between two-stage, also have 2M-2 first meridian window.Respectively take this 2M-2 window as reference, along the parallel at place separately with

Be interval (θ _wLatitude for window :)

M=1 ..., M-1 generates

Individual window.As shown in Figure 5.

Finally obtain

Individual window.As Fig. 6, shown in Figure 7.

According to top step, sphere can be divided on the one hand a series of windows of complete covering sphere, can turn to suitable grid with above-mentioned window is discrete on the other hand.

For the observation data collection { d (x, y, α) } of determining, if given its Theory of Spherical Frame Vectors There is unique azimuthal variations α in Partial Variable, claims such observation data to be distributed on curved surface, that is, and and imagination Be Theory of Spherical Frame Vectors, and α is " radial height ", to be distributed in thickness be on 0 closure, continuous curve surface to observation data.

Like this, the position angle parameter that can also defined function α (x, y) represents observation data on given spherical surface position.

In conjunction with top discretization scheme, arbitrarily local observation data and position angle parameter thereof the form with the plane epigraph can be showed.

Step 2 uses clustering algorithm that parameter alpha different modulating function in position angle is divided into groups, and the value of position angle parameter alpha is divided into one group in the cluster scope with interior modulation kernel function.This step further comprises:

Step 201, the decomposition of modulation kernel function and modulation operation.

The variable α value of modulation kernel function p (x, y, α, x ', y ') in formula (1) is different, means due to rotation, and the position angle of modulating system is different.Former modulation kernel function is decomposed by the position angle:

$p(x,y,\mathrm{\α},x^{\′}{,y}^{\′})={\∫}_0^{2\mathrm{\π}}p(x,y,{\mathrm{\α}}^{\′},x^{\′},y^{\′})\mathrm{\δ}(\mathrm{\α}-{\mathrm{\α}}^{\′}){\mathrm{d\α}}^{\′}.---\left(5\right)$

Therefore modulation methods formula (1) is transformed to:

$d(x,y,\mathrm{\α})-{\∫}_{\mathrm{\Ω}}{\∫}_0^{2\mathrm{\π}}p(x,y,{\mathrm{\α}}^{\′},x^{\′},y^{\′})\∫(x^{\′},y^{\′})\mathrm{\δ}(\mathrm{\α}-{\mathrm{\α}}^{\′}){\mathrm{dx}}^{\′}{\mathrm{dy}}^{\′}{\mathrm{d\α}}^{\′}+n(x,y,\mathrm{\α})$

$={\∫}_0^{2\mathrm{\π}}(p\left({\mathrm{\α}}^{\′}\right)*f)(x,y)\mathrm{\δ}(\mathrm{\α}-{\mathrm{\α}}^{\′}){\mathrm{d\α}}^{\′}+n(x,y,\mathrm{\α}),---\left(6\right)$

(p (α ') * f wherein) (x, y) is p (α ') and the convolution of f,, Modulation equation is decomposed into a series of convolution equations that is.

Step 202, the cluster of discrete form and approximate mode.

In this step, contrast (6), the is olation of discrete form is:

$d_{i,j,k}=\underset{i^{\′}=1,j^{\′}=1}{\overset{M,N}{\mathrm{\Σ}}}{k}^{\′}{P}_{i,j,k^{\′},i^{\′},j^{\′}}{\mathrm{\δ}}_{k-k^{\′}}+n_{i,j,k}$

$={\mathrm{\Σ}}_{k^{\′}=1}{(p_{k^{\′}}*f)}_{i,j}{\mathrm{\δ}}_{k-k^{\′}}+n_{i,j,k},k^{\′}=\mathrm{0,1,2},...,K,---\left(7\right)$

Be approximately:

$d_{i,j,k}\≈\underset{l}{\mathrm{\Σ}}{(p_l*f)}_{i,j}{\mathrm{\δ}}_{|{\mathrm{\α}}_k-{\stackrel{\·}{\mathrm{\α}}}_l|>c}+n_{i,j,k},l=\mathrm{0,1,2},...,L,---\left(8\right)$

Wherein For the Kronecker symbol of promoting, work as subscript

Be true time as the logical expression value, the value of this symbol is 1, otherwise the value of this symbol is 0, α _kBe azimuthal k the sampled point of discretize,

α _kThe cluster centre value of place cluster.

The c that deserves to be called in the face formula is cluster radius, visible approximate process that is the process of k ' value degeneracy, and c is larger, and the degeneracy degree is higher, and k ' is different, and value is fewer, and is approximate more coarse, and calculated amount is also just less.

So-called cluster does not namely wait value but approaching azimuthal variations is approximately certain approximate value.Adopt the approximate cluster member of center of gravity (that is mean value of cluster member) of cluster.

Step 203, clustering parameter is determined.

In this step, the clustering method of employing only has a parameter c to be determined.

Suppose the positional precision that requires be ε (unit: radian), within the point spread function range of definition distance distance between two points farthest be D (unit: radian), the cluster radius c of use should be:

$c=\frac{\mathrm{\ϵ}}{2D},---\left(9\right)$

Like this, the interior site error of same cluster set is no more than

The purpose of hard clustering radius c is that accuracy requirement is still satisfied in the calculating after guaranteeing to be similar to.

Step 204, the design of cluster approximate data.

The target of clustering algorithm is: 1) azimuthal variations of given observation data is divided grouping, make every prescription parallactic angle variable maximum value and minimal value differ and be no more than cluster size c; 2) the azimuthal center of gravity of calculating group is as cluster member's approximate value.

The cluster approximate data realizes in accordance with the following steps:

1) all azimuthal variations are carried out quicksort, it is arranged according to the numerical values recited ascending order, as

{α ₀，α ₁，α ₂，...，α _K-1}；

2) make k=0, l=0, a ^- ₀=0, K ₀=1, g ₀=0:

(a) if k＜K carry out next step; Otherwise finish to return;

(b) order

Carry out next step;

(c) if α _k+1-α _k＜c makes g _k+1=l, K _l=K _l+ 1; Otherwise, order L=l+1,

Make k=k+1, and return to (a).

Like this, by above-mentioned single ergodic, just obtain and azimuthal variations packet marking array { the cluster centre array of gk} and each grouping one to one

Because general observation process position angle is continually varying, proper if cluster size c chooses, array

Scale will be much smaller than array { α _k.

In superincumbent step, K is azimuthal sum, for numbering one by one at K position angle, K azimuthal subscript (numbering) be 0,1,2 ..., K-1, this K numbering is exactly their index, k represents with variable.With lower target K, as K ₀, be azimuthal quantity that different clusters comprises.Subscript represents different clusters.Cluster ordinal number g _kRefer to azimuth angle alpha _kBelong to which cluster.Due to cluster also from 0 open numbering, so the value of cluster ordinal number is also the number since 0.

Step 3, the direct demodulation method approximate based on cluster accelerates.

The front has been introduced method how to utilize cluster to be similar to original discrete form Modulation equation (2) has been approximately formula (8).The item p that comprises in formula (8) _lBe the position angle parameter is existed

Being similar to of modulation kernel function in scope.If observation data is distributed on curved surface, namely there are unique observation data and position angle parameter thereof on sphere everywhere, aforementioned position angle packet marking { g _kAlso can write { g _{I, j}Form.The brief note Kronecker symbol

Be δ _{L, i, j}, that is, if under be designated as i, the observation data position angle of j is in grouping l, this symbol value is 1, otherwise value is 0. like this, formula (8) also can further be reduced to:

d _i，j≈∑δ _l，i，j(p _l*f) _i，j+n _i，j，l＝0，1，2，...，L (10)

Because pl is exactly in fact discrete form

So { pl} is exactly in fact that system point spread function p (x-x ', y-y ', α=0) is around central rotation

Conversion afterwards.For the modulation kernel function of discrete form, because it exists with the form of matrix, can calculate postrotational new matrix by the common method of interpolation for two-dimensional matrix.Convolutional calculation in formula uses fast fourier transform to realize.So far, possessed the full terms that carries out direct demodulation method acceleration realization on sphere.

Iterative after acceleration is:

P ' wherein _lP _lTransposition,

The demodulation result of t step iteration,

That is the iteration initial value is made as observation data, With Respectively Fourier transform and inverse transformation.

About other aspects of direct demodulation method, as consistent with original direct demodulation method in the applying method of background constraint, the applying method of upper limit constraint, and computational complexity is negligible a small amount of than the iterative part that provides above.

The present invention is mainly applied widely with respect to the advantage of prior art, can be applicable to the demodulation that do not have Rotational Symmetry and comprise the modulated process of rotary manipulation, and can according to degree of accuracy require flexible accelerate and precision between equilibrium point.

Above-described specific embodiment; purpose of the present invention, technical scheme and beneficial effect are further described; institute is understood that; the above is only specific embodiments of the invention; be not limited to the present invention; within the spirit and principles in the present invention all, any modification of making, be equal to replacement, improvement etc., within all should being included in protection scope of the present invention.

Claims (7)

1. direct demodulation method of sphere, the method comprising the steps of:

Step 1, according to collimating apparatus visual field size, determine the sphere window size and determine the interior quantity along any limit sampled point of sphere window, the discretize grid of generating reference window, use the discretize grid identical with the reference windows size, that shape is identical to cover whole sphere, use accurate quadrature pixelation mode to carry out discretize to sphere, with sphere discretely turn to that neighbor is equidistant, primitive shape is consistent, the grid of pixel distribution quadrature;

Step 2 uses clustering algorithm that parameter alpha different modulating function in position angle is divided into groups, and the value of position angle parameter alpha is divided into one group in the cluster scope with interior modulation kernel function;

Step 3 is used the acceleration numerical evaluation that realizes modulation operation based on the numerical value approximation computation of cluster.

2. method according to claim 1, is characterized in that, step 2 further comprises:

Step 210 is decomposed into the modulation kernel function form that is rotated conversion by the kernel function that does not contain azimuthal variations, thereby the modulation operation of azimuthal variation is decomposed into modulation operation and the twiddle operation at constant bearing angle;

Step 220 according to the requirement to demodulation computing degree of accuracy, is determined the radius of above-mentioned cluster;

Step 230 according to given cluster radius, is determined each cluster centre and cluster member, makes the position angle of each member in given cluster and the position angle of cluster centre differ in cluster radius;

Step 240 is carried out numerical value with step 210 to the rotation of each kernel function and is approached, and namely is approximately the rotation to cluster centre.

3. method according to claim 2, is characterized in that, step 210 further comprises:

Step 211 is with former modulation kernel function

D (x, y, α)=∫ Ω p (x, y, α, x ', y ') f (x ', y ') dx ' dy '+n (x, y, α) is decomposed into by the position angle:

$p(x,y,\mathrm{\α},x^{\′},y^{\′})={\∫}_0^{2\mathrm{\π}}\mathrm{\Ω}\·p(x,y,\mathrm{\α},x^{\′},y^{\′})f(x^{\′},y^{\′}){\mathrm{dx}}^{\′}{\mathrm{dy}}^{\′}+n(x,y,\mathrm{\α})$

Wherein, p is the modulation kernel function, and x, y, x ' and y ' are respectively the azimuthal coordinates variablees, and α and α ' are azimuthal variations, and δ is the finite point impulse function, and Ω refers to azimuthal span;

Step 212 according to the decomposition result of step 211, with former modulation Kernel Function Transformation is:

$d(x,y,\mathrm{\α})={\∫}_{\mathrm{\Ω}}{\∫}_0^{2\mathrm{\π}}p(x,y,{\mathrm{\α}}^{\′},x^{\′}{,y}^{\′})f(x^{\′},y^{\′})\mathrm{\δ}(\mathrm{\α}-{\mathrm{\α}}^{\′}){\mathrm{dx}}^{\′}{\mathrm{dy}}^{\′}{\mathrm{d\α}}^{\′}+n(x,y,\mathrm{\α})$

$={\∫}_0^{2\mathrm{\π}}(p\left({\mathrm{\α}}^{\′}\right)*f)(x,y)\mathrm{\δ}(\mathrm{\α}-{\mathrm{\α}}^{\′}){\mathrm{d\α}}^{\′}+n(x,y,\mathrm{\α});$

Wherein, d is the output data, and f is the input data, and n is noise item.

4. method according to claim 3, is characterized in that, step 220 further comprises:

Described cluster is spherical cluster, and wherein cluster radius c passes through

Determine, wherein ε is the given positional precision ε of user, and the unit of ε is radian, and D is farthest two dot spacings in the point spread function field of definition, and unit is also radian.

5. method according to claim 4, is characterized in that, step 230 further comprises:

Step 231 is carried out quicksort to all azimuthal variations, it is arranged as: { α according to the numerical values recited ascending order ₀, α ₁, α ₂α _K-1, wherein K is azimuthal quantity;

Step 232 makes position angle index k=0, i.e. α _k=α ₀, cluster index l=0, cluster centre

Cluster member's number K ₀=1, cluster ordinal number g ₀=0:

If a) k＜K carry out next step, otherwise finish to return;

B) order Carry out next step;

C) if α _k+1-α _k＜c makes g _k+1=l, K _l=K _l+ 1; Otherwise, order

L=l+1,

Make k=k+1, and return to (a),

By single ergodic, obtain and azimuthal variations packet marking array { g one to one _kAnd the cluster centre array of each grouping

This result that obtains is used for numerical value and approaches modulation calculating, through the modulation meter that approaches at last according to K position angle rotation modulation kernel function one by one, after obtaining the above results, the modulation kernel function that belongs to a cluster together the angle according to its cluster centre is rotated.

6. method according to claim 5, is characterized in that, step 240 further comprises:

Step 241, in this step, the is olation of discrete form is:

$d_{i,j,k}=\underset{i^{\′}=1,j^{\′}=1}{\overset{M,N}{\mathrm{\Σ}}}{k}^{\′}{P}_{i,j,k^{\′},i^{\′},j^{\′}}{\mathrm{\δ}}_{k-k^{\′}}+n_{i,j,k}$

$={\mathrm{\Σ}}_{k^{\′}=1}{(p_{k^{\′}}*f)}_{i,j}{\mathrm{\δ}}_{k-k^{\′}}+n_{i,j,k},k^{\′}=\mathrm{0,1,2},...,K,---\left(7\right)$

Wherein k is the position angle index, and K is the position angle number;

Step 242, order

$d_{i,j,k}\≈\underset{l}{\mathrm{\Σ}}{(p_l*f)}_{i,j}{\mathrm{\δ}}_{|{\mathrm{\α}}_k-{\stackrel{\·}{\mathrm{\α}}}_l|>c}+n_{i,j,k},l=\mathrm{0,1,2},...,L,---\left(8\right)$

Wherein For the Kronecker symbol of promoting, work as subscript

Be true time as the logical expression value, the value of this symbol is 1, otherwise the value of this symbol is 0, α _kBe azimuthal k the sampled point of discretize,

α _kThe cluster centre value of place cluster, l is the cluster index, and L is total cluster number, and c is cluster radius.

7. method according to claim 6, is characterized in that, step 3 further comprises: the modulation kernel function position angle of same cluster is consistent in given accuracy, according to the cyclic convolution theorem, realizes with fast fourier transform.

Images