WO2021031711A1 - 一种识别位置谱的方法、装置以及计算机存储介质 - Google Patents

一种识别位置谱的方法、装置以及计算机存储介质 Download PDF

Info

Publication number
WO2021031711A1
WO2021031711A1 PCT/CN2020/099959 CN2020099959W WO2021031711A1 WO 2021031711 A1 WO2021031711 A1 WO 2021031711A1 CN 2020099959 W CN2020099959 W CN 2020099959W WO 2021031711 A1 WO2021031711 A1 WO 2021031711A1
Authority
WO
WIPO (PCT)
Prior art keywords
spectrum
position spectrum
extreme
extreme points
points
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.)
Ceased
Application number
PCT/CN2020/099959
Other languages
English (en)
French (fr)
Inventor
赵安江
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Raycan Technology Co Ltd
Original Assignee
Raycan Technology Co Ltd
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Raycan Technology Co Ltd filed Critical Raycan Technology Co Ltd
Priority to US17/597,642 priority Critical patent/US12164071B2/en
Priority to JP2021563135A priority patent/JP7198532B2/ja
Priority to FIEP20855620.9T priority patent/FI4020020T3/fi
Priority to EP20855620.9A priority patent/EP4020020B1/en
Publication of WO2021031711A1 publication Critical patent/WO2021031711A1/zh
Anticipated expiration legal-status Critical
Ceased legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01T—MEASUREMENT OF NUCLEAR OR X-RADIATION
    • G01T1/00—Measuring X-radiation, gamma radiation, corpuscular radiation, or cosmic radiation
    • G01T1/29—Measurement performed on radiation beams, e.g. position or section of the beam; Measurement of spatial distribution of radiation
    • G01T1/2914—Measurement of spatial distribution of radiation
    • G01T1/2985—In depth localisation, e.g. using positron emitters; Tomographic imaging (longitudinal and transverse section imaging; apparatus for radiation diagnosis sequentially in different planes, steroscopic radiation diagnosis)
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01T—MEASUREMENT OF NUCLEAR OR X-RADIATION
    • G01T1/00—Measuring X-radiation, gamma radiation, corpuscular radiation, or cosmic radiation
    • G01T1/16—Measuring radiation intensity
    • G01T1/20—Measuring radiation intensity with scintillation detectors
    • G01T1/202—Measuring radiation intensity with scintillation detectors the detector being a crystal
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01T—MEASUREMENT OF NUCLEAR OR X-RADIATION
    • G01T1/00—Measuring X-radiation, gamma radiation, corpuscular radiation, or cosmic radiation
    • G01T1/29—Measurement performed on radiation beams, e.g. position or section of the beam; Measurement of spatial distribution of radiation
    • G01T1/2914—Measurement of spatial distribution of radiation
    • G01T1/2992—Radioisotope data or image processing not related to a particular imaging system; Off-line processing of pictures, e.g. rescanners
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/23—Clustering techniques
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/70—Determining position or orientation of objects or cameras
    • G06T7/73—Determining position or orientation of objects or cameras using feature-based methods
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V10/00—Arrangements for image or video recognition or understanding
    • G06V10/20—Image preprocessing
    • G06V10/26—Segmentation of patterns in the image field; Cutting or merging of image elements to establish the pattern region, e.g. clustering-based techniques; Detection of occlusion
    • G06V10/267—Segmentation of patterns in the image field; Cutting or merging of image elements to establish the pattern region, e.g. clustering-based techniques; Detection of occlusion by performing operations on regions, e.g. growing, shrinking or watersheds
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V10/00—Arrangements for image or video recognition or understanding
    • G06V10/40—Extraction of image or video features
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
    • A61B6/02—Arrangements for diagnosis sequentially in different planes; Stereoscopic radiation diagnosis
    • A61B6/03—Computed tomography [CT]
    • A61B6/037—Emission tomography

Definitions

  • the present invention relates to the field of image processing, and more specifically to a method, device and computer storage medium for identifying a position spectrum.
  • PET Positron Emission Tomography
  • the basic principle of PET is to label a molecular probe with a radionuclide that can generate positrons and inject it into the organism.
  • the positrons generated by the decay of the radionuclide collide with the negative electrons in the organism and annihilate, and at the same time release a pair of energy
  • All 511keV gamma photons with the opposite direction of motion are converted into electrical signals by using a position-sensitive radiation detector arranged around the biological body to obtain the energy, location and time information of the annihilation event, and further pass the annihilation
  • the position of the response line of the annihilation event is obtained, and the distribution of radionuclides in the organism is obtained through the two-dimensional or three-dimensional tomographic reconstruction algorithm, so as to observe the physiological and biochemical changes in the organism in vitro.
  • PET imaging equipment can obtain the location information of the annihilation event through the radiation detector at the front end.
  • the radiation detector usually adopts a structure in which a crystal array and a photoelectric conversion device (such as a photomultiplier tube or a silicon photomultiplier tube) are coupled, that is, the crystal ( Also known as scintillation crystal) is cut into crystal strips of a certain specification, and then several crystal strips are arranged according to a certain rule to form a crystal array (crystal array can also be called a module), and then the crystal array is coupled with the photoelectric conversion device to form a radiation detector .
  • a large number of gamma photons fly at different angles.
  • the crystal strips When gamma photons are incident on a crystal strip in the crystal array, the crystal strips interact with the gamma photons to produce fluorescence, which is then transferred to the optoelectronics
  • the multiplier device is converted into a corresponding electrical signal output by the photomultiplier device.
  • the coordinates (x, y) of the incident position of the gamma photon can be calculated, and the coordinate system XOY can be established according to the cross section of the crystal array. Finally, according to the position coordinate information, it is judged from which crystal bar the ⁇ photon is incident.
  • the corresponding relationship between the position coordinates and the crystal bar is a linear relationship, that is, when the dynamic range of the position coordinates (x, y) is divided in proportion to the size of the crystal bar, each range is The coordinate area covered by the cross section of a crystal bar corresponds linearly.
  • the linear correspondence relationship always exhibits nonlinear characteristics. .
  • a two-dimensional position scatter plot or a two-dimensional position histogram, which records the number of incident gamma photons measured at each coordinate position.
  • the existing methods for establishing position spectra are all based on experimentally measured scatter plots. Their general idea is to first determine the boundary of each crystal strip, and then use the boundary information to establish the position spectrum. Rogers et al. used the centroid of each small area as the peak point of the area, and then found the local minimum points around it, and used these local minimum points as the boundaries of the corresponding crystal strips.
  • the advantage of this method is that the crystal There will be no detection dead zone on the bar, there is a crystal bar corresponding to any position coordinate, there is no flicker case abandoned, and the sensitivity of the detection system will not be reduced.
  • This method is intuitive and easy to understand, but lacks a rigorous theoretical basis. Stronger and Johnson et al. gave a method to establish a position spectrum for a crystal array based on the Gaussian mixture distribution model (GMM). This method improves the accuracy of the established position spectrum, but it still exists in specific implementation. Where further improvement is needed, for example, multiple use of the mean value deletion method when searching for local peak points is likely to cause accidental deletion.
  • GMM Gaussian mixture distribution model
  • the number of crystal bars is used as the loop termination condition, and the local peak points found may be There are false local peak points, so the true local peak points are missed.
  • semi-automatic or automated location spectrum segmentation methods can also be used, such as using watershed algorithm to obtain extreme points of crystal locations, and unrecognized location points can be found by straight line fitting, and the location spectrum can be segmented.
  • this method also has the problems that the edge of the position spectrum is blurred and the deformation is serious and difficult to accurately identify.
  • the method for realizing position spectrum segmentation in the prior art has poor applicability, and has poor recognition effect on the position spectrum with blurred edges and severe deformation, and the accuracy is low; at the same time, the prior art can only deal with certain characteristics. Effective recognition of the position spectrum cannot solve the problem of different characteristics and large differences in the position spectrum generated by different radiation detectors.
  • the purpose of the present invention is to provide a method, a device and a computer storage medium for identifying a position spectrum, so as to solve at least one of the above-mentioned problems in the prior art.
  • the technical solution of the present invention is to provide a method for identifying the position spectrum.
  • the position spectrum is a two-dimensional position distribution map of photons.
  • the position spectrum contains the position information of all photons extracted from the single event output by the detector. The method includes the following steps:
  • Step S1 preprocessing the initial position spectrum to obtain the first position spectrum
  • Step S2 Perform feature extraction on the first position spectrum to obtain multiple second position spectra
  • Step S3 Perform an extreme value search on the second position spectrum to obtain N extreme value points, the extreme value points form a third position spectrum, where N is the number of scintillation crystals in the detector;
  • Step S4 Globally number the extreme points in the third position spectrum to form a sixth position spectrum
  • Step S5 cluster the single events according to the extreme points in the sixth position spectrum to form a final position spectrum, and complete the segmentation of the position spectrum.
  • the specific steps of the preprocessing are: selecting the size and the first variance of the first Gaussian template, and comparing the initial position spectrum with the first Gaussian template.
  • the template undergoes convolution processing to obtain the first position spectrum.
  • the specific steps of feature extraction include:
  • Step S21 Determine a set of second variances ⁇ 1 , ⁇ 2 ,..., ⁇ n , where n is a natural number, and generate n second Gaussian templates according to each of the second variances ⁇ n and the same size, denoted as g ;
  • Step S22 Find the Hessian matrix after convolution of each of the second Gaussian templates and the first position spectrum
  • Step S23 Calculate the determinants of n Hessian matrices: Obtain n second position spectra.
  • the specific steps of obtaining the Hessian matrix include:
  • the second-order partial derivative of the second Gaussian template with respect to x is calculated according to the different second variances:
  • the second step is to calculate the second-order partial derivative of the second Gaussian template with respect to y according to the different second variances:
  • the third step is to calculate the partial derivative of the second Gaussian template with respect to x according to the different second variances, and then calculate the partial derivative with respect to y:
  • the above three partial derivatives are normalized, that is, multiplied by ⁇ 2 , and then the three partial derivatives are convolved with the first position spectrum respectively to obtain n L xx , L xy , L yy , thus obtaining n Hessian matrices, where g is a second Gaussian template, and x and y are orthogonal coordinate axes corresponding to the cross-section of the scintillation crystal in the detector.
  • step S3 it is necessary to perform morphological expansion processing on a plurality of the second position spectra to obtain the expansion position spectra before performing the extreme value search.
  • the specific step of the extreme value search is to search for the largest N values among the plurality of expansion position spectra to obtain the third position spectra.
  • the specific steps of global numbering include:
  • Step S41 Obtain the position coordinates and position response radius of all the extreme points
  • Step S42 Extract the area where the position coordinates of the extreme point are within the position response radius
  • Step S43 Calculate the minimum value of the corresponding dimensional coordinates in the row/column direction in each of the regions;
  • Step S44 Sort the N minimum values, the extreme points corresponding to the m ⁇ (n-1)+1 to m ⁇ n minimum values belong to the nth row/column, and m is the flicker on the row/column
  • n is a natural number
  • Step S45 Numbering the rows/columns allocated by the extreme points according to the actual geometric order of the scintillation crystal array, and displaying the numbers around the extreme points to form the sixth position spectrum.
  • the event clustering means that after obtaining the position coordinates of the extreme points in the position spectrum, clustering all single events according to the position information to generate a crystal comparison table and Obtain the segmentation result of the position spectrum.
  • the specific steps of event clustering include:
  • Step S51 randomly select one of the position coordinates of the single event, and calculate the distance between the position coordinates and the position coordinates of the N identified extreme points;
  • Step S52 Select the extreme point where the shortest distance is located among the N distances in the step S51 as the crystal cluster of the position coordinates in the step S51;
  • Step S53 Traverse all the position coordinates of the single event, obtain all the crystal clusters corresponding to the position coordinates, and distinguish different categories by a broken line to form the final position spectrum.
  • the method further includes the following steps:
  • Step S31 Eliminate invalid extreme points in the third position spectrum to form a fourth position spectrum
  • Step S32 Distribute rows and columns of all extremum points after the elimination to form a fifth position spectrum
  • Step S33 Predict the missing extreme points in the fifth position spectrum.
  • the specific steps of removing invalid extreme points in the third position spectrum include:
  • Step S311 Extract the location distribution area images of all extreme points according to the position coordinates of the extreme points and the second variance
  • Step S312 Determine an initial segmentation threshold T, mark the ratio of the number of pixels smaller than the initial segmentation threshold T in the position distribution area image to the number of all pixels as w 0 , and mark the average gray as u 0 ; record the position distribution area image
  • the ratio of the number of pixels greater than the initial segmentation threshold T to the number of all pixels is denoted as w 1
  • the average gray level is denoted as u 1 ;
  • Step S314 Traverse the value of the initial segmentation threshold T so that the inter-class variance k is the largest, thereby obtaining a binarized image of the location distribution area image of the scintillation crystal;
  • Step S315 Obtain a neighborhood in the center of the binarized image and calculate the pixel sum. If the value of the pixel sum in the neighborhood is 0, the extreme point is determined as an invalid extreme point.
  • the specific steps of rank allocation may include:
  • Step S321 sort the position coordinates of the extreme points identified in step S31 according to the size of the corresponding dimensional coordinates;
  • Step S322 Extract the m ⁇ (n-1)+1 to m ⁇ nth extreme points as the fuzzy row set of the nth row/column, where m is the number of crystals on each row, and n is a natural number;
  • Step S323 randomly select p extreme points from the fuzzy row set of the nth row/column as the point set S participating in the fitting, and the set of remaining extreme points is recorded as the complement R, and p is the sample required for fitting the curve number;
  • Step S324 Perform quadratic curve fitting on the extreme points in the point set S;
  • Step S325 Obtain and fit the error sum E of the quadratic curve for all extreme points in the complementary set R;
  • Step S326 Repeat the process of step S323-step S325 until the minimum value of the error and E obtained within the number of iterations or the error and E is less than the set minimum tolerance threshold, and the distribution curve of the nth row/column extreme point is obtained;
  • Step S327 Assign the identified extreme points to the distribution curve of each row/column according to the principle of row-column allocation.
  • the rank allocation principle includes the nearest allocation principle and the conflict allocation principle.
  • the principle of nearby allocation includes: the closest allocation according to the distance between the extreme point and the intersection of the distribution curve; the intersection of a distribution curve occupies at most one extreme point, and one extreme point can only be allocated Give a point of intersection.
  • the conflict allocation principle includes: finding the corresponding multiple extreme points through the intersection of the distribution curve; finding the two closest intersection points for each corresponding extreme point; judging each extreme value The size of the sum of the distances between a point and the two closest intersections. The extreme point with the largest sum of distances is assigned to the intersection of the conflicting distribution curves.
  • the specific step of predicting the missing extreme point is: find the coordinate of the intersection of the unassigned extreme point, and assign the position coordinates of the intersection to the current row and The missing extreme points on the current column.
  • the present invention also provides a device for identifying a position spectrum.
  • the device includes: a preprocessing module, a feature extraction module, an extremum search module, a global numbering module, and an event clustering module.
  • the preprocessing module receives the output from the detector An initial position spectrum and preprocessing the initial position spectrum to form a first position spectrum;
  • the feature extraction module receives the first position spectrum and performs feature extraction on the first position spectrum to obtain a second position spectrum;
  • the extreme value search module receives the second position spectrum and performs an extreme value search on the second position spectrum to obtain a third position spectrum;
  • the global numbering module performs a global number on the third position spectrum to form a second position spectrum Six position spectrum;
  • the event clustering module clusters single events in the sixth position spectrum to form a final position spectrum, and completes the segmentation of the position spectrum.
  • the device further includes: an invalid value elimination module, a rank assignment module, and a missing point prediction module.
  • the invalid value elimination module eliminates invalid extreme points in the third position spectrum and forms a first Four-position spectrum; the row and column allocation module allocates extreme points in the fourth position spectrum to form a fifth position spectrum; the missing point prediction module performs operations on the extreme points missing in the fifth position spectrum prediction.
  • the present invention also provides a computer storage medium with a computer program stored on the computer storage medium, and when the computer program is executed, the method as recorded in any of the above embodiments is implemented.
  • any position spectrum can be processed automatically according to needs, and the position spectrum can also be divided and post-processed according to needs, which can be highly efficient and accurate It realizes the recognition of the position spectrum with fuzzy edges and severe deformation, and supports the effective recognition of the position spectrum generated by detectors of different structures.
  • Fig. 1 is a schematic diagram of the steps of a method for identifying a position spectrum according to an embodiment of the present invention
  • FIG. 2 is a schematic diagram of the effect of the method for identifying position spectra according to an embodiment of the present invention, where the initial position spectrum (left) is obtained by a SiPM (silicon photomultiplier) detector, and the first position spectrum can be obtained after filtering.
  • FIG. 3 is a schematic diagram of a second position spectrum according to the method for identifying a position spectrum in FIG. 2, wherein one of the second position spectra is obtained after feature extraction on the first position spectrum;
  • FIG. 4 is a schematic diagram of a position spectrum after morphological expansion of the second position spectrum shown in FIG. 2 according to the method for identifying a position spectrum according to an embodiment of the present invention
  • FIG. 5 is a schematic diagram of a fourth position spectrum according to the method for identifying a position spectrum of FIG. 4, wherein the third position spectrum is obtained after an extreme value search is performed on the second position spectrum;
  • Fig. 6 is a schematic diagram of a final position spectrum after segmentation according to an embodiment of the present invention, wherein the initial position spectrum is obtained by a SiPM detector;
  • Fig. 7 is a schematic diagram of the effect of a method for identifying a position spectrum according to another embodiment of the present invention, in which the initial position spectrum (left) is acquired by a PSPMT (position sensitive photomultiplier tube) detector, and the initial position spectrum can be filtered after filtering. Get the first position spectrum (right);
  • PSPMT position sensitive photomultiplier tube
  • FIG. 8 is a schematic diagram of a second position spectrum according to the method for identifying a position spectrum in FIG. 6, wherein one of the second position spectrums is obtained after feature extraction on the first position spectrum;
  • FIG. 9 is a schematic diagram of a position spectrum after morphological expansion of the second position spectrum shown in FIG. 7 according to the method for identifying a position spectrum according to an embodiment of the present invention.
  • FIG. 10 is a schematic diagram of a third position spectrum obtained after an extreme value search is performed on the position spectrum shown in FIG. 8 in a method for identifying a position spectrum according to an embodiment of the present invention
  • FIG. 11 is a schematic diagram of the fourth position spectrum of the method for identifying position spectrum according to an embodiment of the present invention, wherein the fourth position spectrum is obtained after invalid value elimination of the third position spectrum, and the white circle marks the recognition error that needs to be eliminated Extreme point
  • FIG. 12 is a schematic diagram of an allocation conflict in a method for identifying a position spectrum according to an embodiment of the present invention.
  • FIG. 13 is a schematic diagram of a result of an allocation conflict according to the method for identifying a position spectrum of FIG. 11;
  • FIG. 14 is a schematic diagram of a fifth position spectrum after ranks are allocated according to an embodiment of the present invention, wherein the initial position spectrum is obtained by a SiPM detector;
  • 15 is a schematic diagram of a fifth position spectrum after ranks are allocated according to another embodiment of the present invention, wherein the initial position spectrum is acquired by a PSPMT detector;
  • 16 is a schematic diagram of a sixth position spectrum with a global number in the method for identifying a position spectrum according to an embodiment of the present invention
  • FIG. 17 is a schematic diagram of a final position spectrum after segmentation according to another embodiment of the present invention, wherein the initial position spectrum is acquired by a PSPMT detector;
  • Fig. 18 is a schematic structural diagram of an apparatus for identifying a position spectrum according to an embodiment of the present invention.
  • connection/connection may include electrical and/or mechanical physical connection/connection.
  • including/comprising refers to the existence or addition of features, steps or components/parts, but does not exclude the existence or addition of one or more other features, steps or components/parts.
  • the term “and/or” as used herein includes any and all combinations of one or more of the associated listed items.
  • the method for identifying position spectrum may include the following steps:
  • Step S1 preprocessing the initial position spectrum to obtain the first position spectrum
  • Step S2 Perform feature extraction on the first position spectrum to obtain multiple second position spectra
  • Step S3 Perform an extreme value search on the second position spectrum to obtain N extreme value points, which form a third position spectrum, where N is the number of scintillation crystals in the detector;
  • Step S4 Globally number extreme points in the third position spectrum to form a sixth position spectrum
  • Step S5 cluster the single events according to the extreme points in the sixth position spectrum to form a final position spectrum, and complete the segmentation of the position spectrum.
  • the initial position spectrum can be obtained by equipment or methods commonly used in the art, such as obtaining the initial position spectrum by SiPM detector and multi-threshold sampling method, or by PSPMT detector and multi-threshold sampling method.
  • the threshold sampling method obtains the initial position spectrum, that is, the initial position spectrum in this application may be the position spectrum obtained by any method and device in the prior art.
  • the position spectrum is a two-dimensional position distribution map of all photons formed by the position information extracted from the single event data output by the PET detector, and each pixel value in the image represents the two-dimensional coordinate
  • the light spots in the image indicate the position distribution of photons generated by different crystal strips. The larger the pixel value, the more photons that fall into the position.
  • the initial position spectrum in this application can be a position spectrum with sharp edges and small deformation, or a position spectrum with fuzzy edges and severe deformation. In the art, the sharpness of the edge of the position spectrum and the magnitude of the deformation can be judged by the experience of those skilled in the art, which are not difficult to determine in the art, and will not be repeated here.
  • the specific method for preprocessing the initial position spectrum can adopt two-dimensional Gaussian filtering. Specifically, those skilled in the art can select the size and first variance of the first Gaussian template according to actual experience or randomly. , And then convolve the initial position spectrum with the first Gaussian template to obtain the first position spectrum.
  • the selection of the first Gaussian template, the size and the first variance of the first Gaussian template, and the convolution operation are all conventional technical choices of those skilled in the art, which are not difficult to determine in the art, and will not be repeated here.
  • each spot in the initial position spectrum can be made more uniform and obvious, that is, the shape and distribution of the spots in the initial position spectrum are closer to a two-dimensional Gaussian function, as shown in Figure 2.
  • step S2 since the main feature of the first position spectrum lies in the local light spot, in order to be able to correctly identify the light spot in a wider range, it is necessary to perform feature extraction on the light spot in the first position spectrum. Therefore, based on The derivative method is used to judge the similarity between the light spot and the differential operator according to the light spot distribution characteristics.
  • the distribution feature of the light spot refers to the distribution of the light spot on the two-dimensional coordinate, which can be the arrangement shape of the light spot, etc., such as high middle count and low surrounding count, that is, the middle of the light spot is clear and the surrounding is weak.
  • the specific steps of feature extraction may include:
  • Step S21 Determine a set of second variances ⁇ 1 , ⁇ 2 ,..., ⁇ n (where n is a natural number), and generate a two-dimensional Gaussian template according to each second variance ⁇ n and the same size, that is, generate n
  • the second Gaussian template is denoted as g, and an orthogonal coordinate system XOY corresponding to the image of the position spectrum is established according to the cross section of the scintillation crystal in the detector;
  • Step S22 Find the solution after convolution of each second Gaussian template and the first position spectrum, that is, find the Hessian matrix:
  • Step S23 Calculate the determinant of n Hessian matrices, namely Obtain n second position spectra, as shown in Figure 3.
  • step S21 those skilled in the art can select the specific value of the second variance based on actual experience or randomly, which is not difficult to determine in the art, and will not be repeated here.
  • the n second position spectra can be morphologically expanded to obtain n dilated position spectra, as shown in Figure 4; the extremum search is to search for the largest value among the n dilated position spectra.
  • the position coordinates and the second variance value corresponding to the N values are obtained at the same time to form the third position spectrum.
  • N is the number of corresponding scintillation crystals in the PET detector, and the position corresponding to the N maximum value is Identified N extreme points, the expansion position spectrum where each extreme point is located corresponds to the size of the second Gaussian template variance; usually 3 times the second Gaussian template variance is taken as the third position spectrum of each crystal Radius, as shown in Figure 5.
  • step S4 can be performed.
  • the global numbering can be performed by sorting the neighborhood of extreme points, and the specific steps can include:
  • Step S41 Obtain the position coordinates and position response radius of all extreme points, where the position coordinates of the extreme points are consistent with the originally established XOY coordinate system, and the position response radius can be selected according to the actual cross-sectional size of the crystal bar, such as usual Take 3 times the second Gaussian template variance as the position response radius of each extreme point;
  • Step S42 Extract the area where the position coordinates of the extreme points are within the radius of the position response. If the number of extreme points is N, then the number of extracted areas is also N;
  • Step S43 Calculate the minimum value of the corresponding coordinates in each area in the row/column direction; specifically, when judging rows, it is necessary to judge the smallest y value among all points in the entire area; when judging columns, it needs to judge The smallest x value among all points in the entire area;
  • Step S44 Sort the N minimum values, then the extreme points corresponding to the m ⁇ (n-1)+1 to m ⁇ n minimum values belong to the nth row/column, and m is the row/column The number of crystals, n is a natural number;
  • Step S45 Number the rows/columns assigned by the extreme points in the actual geometric order of the crystal array, and display the numbers around the extreme points to form a sixth position spectrum, as shown in Figure 6, where the numbers refer to the global Numbered result.
  • event clustering means that after obtaining the position coordinates of the extreme points in the position spectrum, all single events need to be clustered according to the position information to generate a crystal lookup table. Since the annihilation of the positron and the electron produces a pair of photons with the same energy and opposite directions, each detected photon event contains the annihilation position, time and energy information. In order to accurately determine the position information, it is necessary to analyze all the single photons. The events are clustered to determine which crystal bar each photon event originates from, and each photon event corresponds to the spot of the position spectrum. Specifically, event clustering may include the following steps:
  • Step S51 arbitrarily select the position coordinates of a single event, and calculate the distance between the position coordinates and the position coordinates of the N recognized extreme points, where N is the number of crystals;
  • Step S52 selecting the crystal with the shortest distance among the N distances in step S51 as the crystal cluster of the position coordinates in step S51;
  • Step S53 Traverse the position coordinates of all single events to obtain crystal clusters of all position coordinates, thereby forming a final position spectrum, as shown in FIG. 6.
  • the position of the extreme point of the identified crystal can also be used to initialize the Gaussian function mean in GMM (Gaussian Mixture Model), initialize the initial neuron weights in SOM (Self-Organizing Mapping Algorithm), and initialize The center of the gathering point in the KNN (K-nearest neighbor algorithm) algorithm can efficiently obtain single-event crystal clusters, and the specific steps will not be repeated.
  • GMM Global System Mixture Model
  • SOM Self-Organizing Mapping Algorithm
  • KNN Know-nearest neighbor algorithm
  • the position coordinate of the single event is essentially the pixel number of the position spectrum.
  • the light spot represents the crystal.
  • the total pixel value is 256 ⁇ 256, indicating that the photon may hit 256 ⁇ 256 positions. Therefore, for a position spectrum of 256 ⁇ 256 pixels,
  • step S1-step S5 can process and identify the position spectrum with higher accuracy, faster speed, higher efficiency, and is especially suitable for the position spectrum with neatly arranged and clear distribution positions, as shown in Figure 2.
  • Shown is the initial position spectrum acquired by SiPM (silicon photomultiplier tube) detector.
  • SiPM silicon photomultiplier tube
  • PSPMT Part Sensitive Photomultiplier Tube
  • the steps of the method for identifying a position spectrum provided by the present invention may further include the following steps:
  • Step S1 preprocessing the initial position spectrum to obtain the first position spectrum
  • Step S2 Perform feature extraction on the first position spectrum to obtain multiple second position spectra
  • Step S3 Perform an extremum search on the second position spectrum to obtain N extremum points, where N is the number of scintillation crystals in the detector; wherein, after the extremum search, steps S31 to S33 are required;
  • Step S31 Eliminate invalid extreme points in the third position spectrum to obtain a fourth position spectrum
  • Step S32 Distribute ranks and columns of all extremum points after being eliminated to form a fifth position spectrum
  • Step S33 Predict the missing extreme points to form a third position spectrum
  • Step S4 Globally number extreme points in the third position spectrum to form a sixth position spectrum
  • Step S5 cluster single events according to the extreme points in the sixth position spectrum to form a final position spectrum, and complete the segmentation of the position spectrum.
  • step S1 and step S2 are similar or the same as in the previous embodiment, and will not be repeated here.
  • the arrangement is disorderly ,
  • the distribution is fuzzy, the effect after preprocessing is shown in Figure 7; the position spectrum after feature extraction is shown in Figure 8, and the effects of morphological expansion and extreme value search are shown in Figures 9 and 10, respectively.
  • the specific steps of removing invalid extreme points may include:
  • Step S311 Extract the location distribution area images of all extreme points according to the position coordinates of the extreme points identified by the extreme search and the second variance, as shown in FIG. 10;
  • Step S312 Determine an initial segmentation threshold T, mark the ratio of the number of pixels smaller than the initial segmentation threshold T in the position distribution area image to the number of all pixels as w 0 , and mark the average gray as u 0 ; record the position distribution area image
  • the ratio of the number of pixels greater than the initial segmentation threshold T to the number of all pixels is denoted as w 1
  • the average gray level is denoted as u 1
  • the average gray level of the overall position distribution area image is denoted as u;
  • Step S314 Traverse the value of the initial segmentation threshold T to maximize the inter-class variance k, thereby obtaining a binarized image of the location distribution area image of the crystal;
  • Step S315 Obtain a neighborhood in the center of the binarized image and calculate the pixel sum. If the value of the pixel sum in the neighborhood is 0, the extreme point is determined to be an incorrectly identified extreme point, that is, an invalid extreme point. As shown in Figure 11, the white circle marks the extreme point of the recognition error.
  • step S32 the specific steps of row and column allocation can be performed according to the rows or columns of the crystal arrangement.
  • row allocation is taken as an example.
  • the specific steps of row and column allocation may include:
  • Step S321 sort the position coordinates of the extreme points identified in step S31 according to the size of the corresponding dimensional coordinates;
  • Step S322 If sorting according to the x coordinate, extract the m ⁇ (n-1)+1 to m ⁇ nth extreme points as the fuzzy row set of the nth row, where m is the number of crystals on each row, and n is Natural number;
  • Step S323 randomly select p extreme points from the fuzzy line set of the nth row as the point set S participating in the fitting, and the set of remaining extreme points is recorded as the remainder R, and p is the number of samples required to fit the curve;
  • Step S324 Perform quadratic curve fitting based on the least square method on the extreme points in the point set S;
  • Step S325 Obtain and fit the error sum E of the quadratic curve for all extreme points in the complementary set R;
  • Step S326 Repeat the process from step S323 to step S325 until the minimum value of the error and E obtained in the number of iterations or the error and E is less than the set minimum tolerance threshold.
  • the fitted quadratic curve is the nth row Distribution curve of extreme points;
  • Step S327 Assign the identified extreme points to the distribution curve of each row/column according to the principle of row-column allocation.
  • the minimum tolerance threshold of error and E can be determined through multiple experiments, and those skilled in the art can easily determine the range of the minimum tolerance threshold through multiple test results, which will not be repeated here.
  • the essence of the rank allocation principle is to allocate extreme points on the position spectrum to the intersection of the row and column distribution curves.
  • the rank allocation principle includes the nearest allocation principle and the conflict allocation principle.
  • the principle of nearest distribution includes: first, the closest distribution is based on the distance between the extreme point and the intersection of the distribution curve; second, the intersection of a distribution curve occupies at most one extreme point, and one extreme point can only be assigned to one Intersection.
  • the two nearest intersections of extreme point A are a and c
  • the two nearest intersections of extreme point B are c and d. ;
  • the extreme point with the largest sum of distances is assigned to the intersection of the conflicting distribution curve, for example, the extreme point A and The distances between the intersection points c and a are l A1 and l A2 respectively, the distances between the extreme point B and the intersection points c and d are l B1 and l B2 respectively , and the extreme point A and the two closest intersection points c and a
  • the sum of the distances between l A1 + l A2 is less than the sum of the distances between extreme point B and the two closest intersection points c and d l B1 + l B2 , so extreme point B is assigned to intersection c, extreme point A is allocated to another intersection a, which is closest to it, and the allocation result is shown in Figure 13.
  • step S33 the prediction of the missing extreme points is performed based on the results of the distribution curve intersection row and column allocation. After the row and column allocation, all the identified extreme points have and only one intersection of the distribution curve corresponds to it. Missing point prediction is to find the intersection coordinates of the unassigned extreme points, and assign the intersection coordinates to the missing extreme points on the current row and the current column.
  • step S4 can be carried out directly according to the result of rank assignment, here No longer.
  • the present invention provides method operation steps as described in the above-mentioned embodiments or flowcharts, more or fewer operation steps may be included in the method based on conventional or no creative labor. In steps where there is no necessary causality logically, the execution order of these steps is not limited to the execution order provided in the embodiment of the present invention.
  • any position spectrum can be processed according to needs, and the position spectrum can be divided and post-processed according to needs.
  • the position spectrum For example, for the position spectrum with blurred edges and severe sexuality, you can directly follow the steps S1-step S2-step S3-step S31, step S32, step S33-step S4-step S5 for segmentation, for neatly arranged and clearly distributed position spectrum, you can follow steps S1-step S2-step S3-step S4-step S5 performs division.
  • the method provided by the present invention can realize the identification of the position spectrum with fuzzy edges and severe deformation with high efficiency and high accuracy, and support the effective identification of the position spectrum generated by detectors of different structures.
  • the present invention also provides a device for identifying position spectrum based on the above method.
  • the device includes a preprocessing module 10, a feature extraction module 20, an extremum search module 30, a global numbering module 40 and The event clustering module 50, wherein the preprocessing module 10 receives the initial position spectrum output by the detector and preprocesses the initial position spectrum to form a first position spectrum; the feature extraction module 20 receives the first position spectrum and analyzes the first position The spectrum performs feature extraction to obtain the second position spectrum; the extreme value search module 30 receives the second position spectrum and performs extreme value search on the second position spectrum to obtain the third position spectrum; the global numbering module 40 performs global numbering on the third position spectrum To form a sixth position spectrum; the event clustering module 50 clusters single events in the sixth position spectrum to form a final position spectrum, and completes the segmentation of the position spectrum.
  • the initial position spectrum can be obtained by equipment or methods commonly used in the art, such as obtaining the initial position spectrum by a SiPM detector and a multi-threshold sampling method, or obtaining the initial position spectrum by a PSPMT detector and a traditional time interval sampling method.
  • the preprocessing module 10 may preprocess the initial position spectrum by using two-dimensional Gaussian filtering.
  • the preprocessing module 10 may set the size of the first Gaussian template and the first party based on actual experience or randomly set by those skilled in the art. Then, the initial position spectrum is convolved with the first Gaussian template to obtain the first position spectrum.
  • each spot in the initial position spectrum can be made more uniform and obvious, that is, the shape and distribution of the spots in the initial position spectrum are closer to a two-dimensional Gaussian function, as shown in Figure 2.
  • the feature extraction module 20 may adopt a derivative-based differentiation method to determine the light spot and the light spot distribution characteristics.
  • the distribution feature of the light spot refers to the distribution of the light spot on the two-dimensional coordinate, which can be the arrangement shape of the light spot, such as high middle count and low surrounding count, that is, the middle light spot is clear and the surrounding light spot is weak.
  • the feature extraction module 20 can perform feature extraction according to the following specific steps:
  • Step S21 preset a set of second variances ⁇ 1 , ⁇ 2 ,..., ⁇ n (where n is a natural number), and generate a two-dimensional Gaussian template according to each second variance ⁇ n and the same size, that is, generate n
  • a second Gaussian template denoted as g;
  • Step S22 Find the solution after convolution of each second Gaussian template and the first position spectrum, that is, find the Hessian matrix:
  • Step S23 Calculate the determinant of n Hessian matrices, namely Obtain n second position spectra, as shown in Figure 3.
  • the extremum search module 30 can perform morphological expansion operations on the n second position spectra to obtain n expansion position spectra, and then the extremum search module 30 searches for the largest N values in the n expansion position spectra as extreme points to obtain
  • the third position spectrum, N is the number of corresponding scintillation crystals in the PET detector.
  • the global numbering module 40 can be used for global numbering.
  • the global numbering can be passed through the extreme point adjacent
  • the specific steps can include:
  • Step S41 Obtain the position coordinates and position response radius of all extreme points, where the position coordinates of the extreme points are consistent with the originally established XOY coordinate system, and the position response radius can be selected according to the actual cross-sectional size of the crystal bar, such as usual Take 3 times the second Gaussian template variance as the position response radius of each extreme point;
  • Step S42 Extract the area where the position coordinates of the extreme points are within the radius of the position response. If the number of extreme points is N, then the number of extracted areas is also N;
  • Step S43 Calculate the minimum value of the corresponding coordinates in each area in the row/column direction; specifically, when judging rows, it is necessary to judge the smallest y value among all points in the entire area; when judging columns, it needs to judge The smallest x value among all points in the entire area;
  • Step S44 Sort the N minimum values, then the extreme points corresponding to the m ⁇ (n-1)+1 to m ⁇ n minimum values belong to the nth row/column, and m is the row/column The number of crystals, n is a natural number;
  • Step S45 Number the rows/columns assigned by the extreme points in the actual geometric order of the crystal array, and display the numbers around the extreme points to form a sixth position spectrum, as shown in Figure 6, where the numbers refer to the global Numbered result.
  • the event clustering module 50 clusters all single events according to the location information, and generates a crystal lookup table.
  • the event clustering module 50 can process according to the following steps:
  • Step S51 arbitrarily select the position coordinates of a single event, and calculate the distance between the position coordinates and the position coordinates of the N recognized extreme points, where N is the number of crystals;
  • Step S52 selecting the crystal with the shortest distance among the N distances in step S51 as the crystal cluster of the position coordinates in step S51;
  • Step S53 Traverse the position coordinates of all single events to obtain crystal clusters of all position coordinates, thereby forming a final position spectrum, as shown in FIG. 6.
  • the position spectrum with disorderly arrangement and fuzzy distribution such as the initial position spectrum obtained by PSPMT (Position Sensitive Photomultiplier Tube) detector shown in Fig. 7, it passes through the preprocessing module 10, the feature extraction module 20, and the extreme value search module After the processing of 30, many wrong extreme points may be obtained. Therefore, it is necessary to use the invalid value elimination module 60, the rank allocation module 70 and the missing point prediction module 80 after the extreme value search module to improve the fuzzy distribution of crystals. The effect of identification and prediction.
  • the processing steps of the invalid value elimination module 60, the rank assignment module 70, and the missing point prediction module 80 are the same as the invalid value elimination step, rank assignment step, and missing point prediction step in the second embodiment described above, and will not be repeated here.
  • the invalid value elimination module 60 performs adaptive error filtering based on the features of the wrong recognition area, which can solve the problem that the position spectrum with fuzzy distribution will get wrong extreme points after feature extraction, and can effectively eliminate most of the wrong extreme points. Value point location.
  • the rank assignment module 70 can make all the recognized extreme points have and only one intersection point of the fitted curve corresponds to it.
  • the missing point prediction module 80 can obtain the coordinates of the intersection of the unassigned extreme point, and assign the coordinates of the intersection to the missing extreme points on the current row and the current column.
  • any position spectrum can be processed according to needs, and the position spectrum can be divided and post-processed according to needs. For example, for a neatly arranged and clearly distributed position spectrum, it can be directly pre-processed. Processing module, feature extraction module, extreme value search module, global numbering module and event clustering module for segmentation; for location spectrum with fuzzy edges and severe deformation, invalid value removal module, rank assignment module and missing point prediction module can be additionally used segmentation.
  • the device provided by the present invention can realize the identification of the position spectrum with fuzzy edges and severe deformation with high efficiency and high accuracy, and support the effective identification of the position spectrum generated by detectors of different structures.
  • the present invention also provides a computer storage medium that stores a computer program.
  • the steps in the above method embodiments can be implemented, such as: preprocessing the initial position spectrum to obtain the first A position spectrum; feature extraction is performed on the first position spectrum to obtain multiple second position spectra; the extreme value search is performed on the second position spectrum to obtain N extreme value points, and the third position spectrum is formed by these extreme value points, where , N is the number of scintillation crystals in the detector; global number the extreme points in the third position spectrum to form the sixth position spectrum; cluster single events according to the extreme points in the sixth position spectrum to form the final Position spectrum, complete the segmentation of position spectrum.
  • invalid extreme points in the third position spectrum can be eliminated to obtain the fourth position spectrum, and all extreme points after the elimination can be allocated to form the fifth position spectrum.
  • the computer storage medium of the present invention may include any entity or device capable of carrying computer program code, such as ROM/RAM, magnetic disk, optical disk, flash memory and other memories, which will not be listed here.
  • the devices, modules, etc. described in the foregoing embodiments may be specifically implemented by computer chips and/or entities, or implemented by products with certain functions.
  • the functions are divided into various modules and described separately.
  • the functions of each module can be integrated into the same or multiple computer chips.

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Theoretical Computer Science (AREA)
  • Health & Medical Sciences (AREA)
  • High Energy & Nuclear Physics (AREA)
  • Molecular Biology (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Data Mining & Analysis (AREA)
  • Multimedia (AREA)
  • Crystallography & Structural Chemistry (AREA)
  • Chemical & Material Sciences (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Evolutionary Biology (AREA)
  • Evolutionary Computation (AREA)
  • General Engineering & Computer Science (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Artificial Intelligence (AREA)
  • Measurement Of Radiation (AREA)
  • Image Processing (AREA)
  • Nuclear Medicine (AREA)
  • Image Analysis (AREA)

Abstract

一种识别位置谱的方法、装置以及计算机存储介质，方法包括：对初始位置谱进行预处理以得到第一位置谱(S1)；对第一位置谱进行特征提取，得到多个第二位置谱(S2)；对第二位置谱进行极值搜索，得到N个极值点并形成第三位置谱(S3)；对第三位置谱中的极点值进行全局编号，形成第六位置谱(S4)；根据第六位置谱中的极值点对单事件进行聚类以获得最终位置谱(S6)。装置包括预处理模块(10)、特征提取模块(20)、极值搜索模块(30)、全局编号模块(40)和事件聚类模块(50)。计算机存储介质的程序被执行时实现识别位置谱的方法。识别位置谱的方法可根据需要全自动化的对任意的位置谱进行处理，可高效率、高准确度地实现对边缘模糊、形变严重的位置谱的识别，支持对不同结构探测器生成的位置谱进行有效识别。

Description

一种识别位置谱的方法、装置以及计算机存储介质 技术领域
本发明涉及图像处理领域，更具体地涉及一种识别位置谱的方法、装置以及计算机存储介质。
背景技术
正电子发射断层成像(Positron Emission Tomography，以下简称PET)是一种非侵入式的造影方法。PET的基本原理是将能够产生正电子的放射性核素标记于分子探针上并注入生物体内部，放射性核素衰变产生的正电子与生物体内的负电子碰撞发生湮灭，同时释放出一对能量均为511keV且运动方向相反的γ光子，通过采用环绕生物体布置的位置灵敏型辐射探测器将接收的γ光子转换为电信号，从而获得湮灭事件的能量、位置和时间信息，进一步地通过湮灭符合技术，得到湮灭事件所在响应线的位置，并通过二维或三维断层重建算法获得放射性核素在生物体内部的分布，从而实现在体外观测生物体内的生理、生化变化过程。
PET成像设备通过位于前端的辐射探测器可以获取湮灭事件的位置信息，辐射探测器通常采用晶体阵列与光电转换器件(比如光电倍增管或者硅光电倍增管)耦合的结构形式，即先将晶体(也称为闪烁晶体)切割成一定规格的晶体条，再由若干个晶体条按一定规则排列组成晶体阵列(晶体阵列也可称为模块)，之后将晶体阵列与光电转换器件耦合形成辐射探测器。当进行辐射探测时，数量众多的γ光子以不同角度飞行，当γ光子入射到晶体阵列中的某一根晶体条上时，晶体条与γ光子相互作用产生荧光，这些荧光随之传递至光电倍增器件并通过光电倍增器件转换为相应的电信号输出，利用该电信号，可以计算出γ光子入射位置的坐标(x，y)，坐标系XOY可以根据晶 体阵列的横截面建立。最后根据位置坐标信息判断γ光子是从哪个晶体条入射的。
理想情况下，位置坐标和晶体条之间的对应关系为线性关系，也就是说，当按照晶体条的尺寸大小成比例的划分位置坐标(x，y)的动态变化范围时，每个范围与一个晶体条的截面所覆盖的坐标区域线性对应。然而在实际应用中，由于受到辐射探测器空间响应非线性、晶体条制作规格不完全一致、物理特性存在一定的差异以及康普顿散射等因素的影响，使得线性对应关系总是呈现非线性特征。为确保对应关系的准确性，就必须建立位置谱，在位置谱中为晶体条和位置坐标(x，y)定义准确的对应关系。为了满足PET成像设备的高位置分辨率的要求，这一对应关系就必须非常准确，因此，怎样建立合理正确的位置谱就成为本领域关注的焦点。
建立位置谱时，通常需要用到二维位置散点图，或者称为二维位置柱状图，该图记录了每个坐标位置上测量到的入射γ光子的事例数。现有的几种建立位置谱的方法都是基于实验测量的散点图，它们的总体思路都是先确定出每个晶体条的边界，然后利用边界信息建立位置谱。Rogers等人用每个小区域的质心点作为该区域的峰值点，然后找出其周围的局部极小值点，以这些局部极小值点作为对应晶体条的边界，该方法的优点在于晶体条上不会出现探测死区，任意位置坐标都有晶体条与之对应，没有被舍弃的闪烁事例，不会降低探测系统的灵敏度。这种方法比较直观易于理解，但缺乏严谨的理论基础。Stronger和Johnson等人以高斯混合分布模型(Gaussian mixture model，简称GMM)为理论基础给出了为晶体阵列建立位置谱的方法，该方法提高了所建位置谱的精度，然而具体实现时依然存在需要进一步完善的地方，比如，在寻找局部峰值点时多次使用平均值删除法容易造成误删除，同时，在寻找局部峰值点时采用晶体条数目为循环中止条件，找到的局部峰值点中可能存在伪局部峰值点，从而漏选真正的局部峰值点。现有技术中还可以通过半自动或者自动化的位置谱分割的方法，比如使用分水岭算法来获得晶体位置的极值点，通过直线拟合的方式来查找未识别的位置点，并分割出位置谱。但是该方 法也存在位置谱边缘模糊、形变严重难以准确识别的问题。
总之，现有技术中实现位置谱分割的方法适用性较差，对于边缘模糊、形变严重的位置谱的识别效果不佳，准确度较低；同时，现有技术只能针对具有某种特征的位置谱进行有效识别，无法解决不同的辐射探测器生成的位置谱的特征不一、差异较大的问题。
发明内容
本发明的目的是提供一种识别位置谱的方法、装置以及计算机存储介质，从而解决上述现有技术中的至少一种问题。
为了解决上述技术问题，本发明的技术方案是提供一种识别位置谱的方法，位置谱为光子的二维位置分布图，位置谱包含从探测器输出的单事件中提取所有光子的位置信息，该方法包括以下步骤：
步骤S1：对初始位置谱进行预处理以得到第一位置谱；
步骤S2：对所述第一位置谱进行特征提取，得到多个第二位置谱；
步骤S3：对所述第二位置谱进行极值搜索，得到N个极值点，所述极值点形成第三位置谱，其中，N为所述探测器中闪烁晶体的数量；
步骤S4：对所述第三位置谱中的所述极值点进行全局编号，形成第六位置谱；
步骤S5：根据所述第六位置谱中的所述极值点对所述单事件进行聚类，形成最终位置谱，完成位置谱的分割。
根据本发明的一个实施例，在所述步骤S1中，所述预处理的具体步骤为：选取第一高斯模板的尺寸大小和第一方差，将所述初始位置谱与所述第一高斯模板进行卷积处理以得到第一位置谱。
根据本发明的一个实施例，在所述步骤S2中，特征提取的具体步骤包括：
步骤S21：确定一组第二方差σ 1,σ 2,…,σ n，n为自然数，根据每一个所述第二方差σ n和相同的尺寸大小生成n个第二高斯模板，记为 g；
步骤S22：求每一个所述第二高斯模板与第一位置谱卷积后的Hessian矩阵，
Figure PCTCN2020099959-appb-000001
，其中，x、y为所述位置谱的图像中对应的正交坐标轴；
步骤S23：计算n个Hessian矩阵的行列式：
Figure PCTCN2020099959-appb-000002
得到n个第二位置谱。
根据本发明的一个实施例，在所述步骤S22中，求Hessian矩阵的具体步骤包括：
第一步，根据不同的所述第二方差计算所述第二高斯模板关于x的二阶偏导：
Figure PCTCN2020099959-appb-000003
第二步，根据不同的所述第二方差计算所述第二高斯模板关于y的二阶偏导：
Figure PCTCN2020099959-appb-000004
第三步，根据不同的所述第二方差计算所述第二高斯模板关于x的偏导后再计算关于y的偏导：
Figure PCTCN2020099959-appb-000005
第四步，对上述三种偏导进行规范化，即乘上σ 2，之后将三种偏导分别与第一位置谱进行卷积操作，得到n个L xx，L xy，L yy，从而获得n个Hessian矩阵，其中g为第二高斯模板，x、y为所述探测器中所述闪烁晶体的横截面对应的正交坐标轴。
根据本发明的一个实施例，在上述步骤S3中，进行极值搜索之 前需要对多个所述第二位置谱进行形态学膨胀处理以得到膨胀位置谱。
根据本发明的一个实施例，在上述步骤S3中，极值搜索的具体步骤是在多个所述膨胀位置谱中搜索最大的N个值以得到所述第三位置谱。
根据本发明的一个实施例，在所述步骤S4中，全局编号的具体步骤包括：
步骤S41：获取所有所述极值点的位置坐标和位置响应半径；
步骤S42：提取所述极值点的所述位置坐标在所述位置响应半径范围内的区域；
步骤S43：计算每个所述区域内在行/列方向上的相应维度坐标的最小值；
步骤S44：对N个最小值进行排序，第m×(n-1)+1到第m×n个最小值对应的极值点隶属于第n行/列，m为行/列上的闪烁晶体数量，n为自然数；
步骤S45：根据所述极值点所分配的行/列按闪烁晶体阵列实际的几何顺序进行编号，将编号显示在极值点周围以形成所述第六位置谱。
根据本发明的一个实施例，在所述步骤S5中，事件聚类是指获得位置谱中极值点的所述位置坐标后，对所有单事件按照位置信息进行聚类，生成晶体对照表并获得位置谱的分割结果。
根据本发明的一个实施例，事件聚类的具体步骤包括：
步骤S51：任意选取一个所述单事件的位置坐标，计算所述位置坐标与N个所识别出的所述极值点的所述位置坐标之间的距离；
步骤S52：在所述步骤S51的N个距离中选取最短的距离所在的极值点作为所述步骤S51中所述位置坐标的晶体聚类；
步骤S53：遍历所有的所述单事件的位置坐标，获得所有所述位置坐标对应的晶体聚类，将不同类别用折线区分开来以形成所述最终位置谱。
根据本发明的一个实施例，该方法还包括以下步骤：
步骤S31：对所述第三位置谱中无效的极值点进行剔除以形成第四位置谱；
步骤S32：对剔除后的所有极值点进行行列分配以形成第五位置谱；
步骤S33：对所述第五位置谱中缺失的极值点进行预测。
根据本发明的一个实施例，在所述步骤S31中，对所述第三位置谱中无效的极值点进行剔除的具体步骤包括：
步骤S311：根据所述极值点的位置坐标和第二方差大小提取出所有极值点的位置分布区域图像；
步骤S312：确定一个初始分割阈值T，将位置分布区域图像中小于初始分割阈值T的像素点的数目占所有像素点数目的比例记为w 0，平均灰度记为u 0；将位置分布区域图像中大于初始分割阈值T的像素点的数目占所有像素点数目的比例记为w 1，平均灰度记为u 1；
步骤S313：按照公式k＝w 0×w 1×(u 0-u 1) 2计算类间方差k；
步骤S314：遍历初始分割阈值T的值，使得类间方差k最大，从而获得闪烁晶体的位置分布区域图像的二值化图像；
步骤S315：取得二值化图像中心一个邻域并计算像素和，若该邻域内的像素和的值为0，则认定该极值点为无效的极值点。
根据本发明的一个实施例，在所述步骤S32中，行列分配的具体步骤可以包括：
步骤S321：将步骤S31中识别的极值点的位置坐标按照对应的维度坐标的大小分别进行排序；
步骤S322：抽取第m×(n-1)+1到第m×n个极值点作为第n行/列的模糊行集，m为每一行上的晶体数量，n为自然数；
步骤S323：从第n行/列的模糊行集中随机抽取p个极值点作为参与拟合的点集S，剩余极值点的集合记为余集R，p为拟合曲线所需的样本数；
步骤S324：对点集S中的极值点进行二次曲线拟合；
步骤S325：对余集R中的所有极值点求取与拟合的二次曲线的误差和E；
步骤S326：重复步骤S323-步骤S325的过程，直至在迭代次数内获得的误差和E的最小值或者误差和E小于设定的最小容忍阈值，得到第n行/列极值点的分布曲线；
步骤S327：将已经识别的极值点按照行列分配原则分配至每一行/列的分布曲线上。
根据本发明的一个实施例，所述行列分配原则包括就近分配原则和冲突分配原则。
根据本发明的一个实施例，所述就近分配原则包括：按照极值点和分布曲线交点之间的距离最近分配；一个分布曲线的交点处最多占据一个极值点，一个极值点仅能分配给一个交点。
根据本发明的一个实施例，所述冲突分配原则包括：通过分布曲线的交点找到对应的多个极值点；对每个对应的极值点分别寻找最近的两个交点；判断每个极值点距离两个最近的交点的距离之和的大小，距离之和最大的极值点分配至有冲突的分布曲线的交点。
根据本发明的一个实施例，在所述步骤S33中，对缺失的极值点进行预测的具体步骤是：找到未分配极值点的交点坐标，并且将该交点的位置坐标赋予给当前行和当前列上缺失的极值点。
本发明还提供一种识别位置谱的装置，所述装置包括：预处理模块、特征提取模块、极值搜索模块、全局编号模块以及事件聚类模块，所述预处理模块接收通过探测器输出的初始位置谱并对所述初始位置谱进行预处理以形成第一位置谱；所述特征提取模块接收所述第一位置谱并对所述第一位置谱进行特征提取以得到第二位置谱；所述极值搜索模块接收所述第二位置谱并对所述第二位置谱进行极值搜索以得到第三位置谱；所述全局编号模块对所述第三位置谱进行全局编号以形成第六位置谱；所述事件聚类模块对所述第六位置谱中的单事件进行聚类以形成最终位置谱，完成位置谱的分割。
根据本发明的一个实施例，所述装置进一步包括：无效值剔除模 块、行列分配模块以及缺失点预测模块，所述无效值剔除模块剔除所述第三位置谱中无效的极值点并形成第四位置谱；所述行列分配模块对所述第四位置谱中的极值点进行分配以形成第五位置谱；所述缺失点预测模块对所述第五位置谱中缺失的极值点进行预测。
本发明还提供了一种计算机存储介质，该计算机存储介质上存储有计算机程序，该计算机程序被执行时实现如上述任一实施例中记载的方法。
通过本发明提供的位置谱分割的方法、装置以及计算机存储介质，可以根据需要全自动化的对任意的位置谱进行处理，也可以根据需要对位置谱进行划分后处理，可以高效率、高准确度地实现对边缘模糊、形变严重的位置谱的识别，同时支持对不同结构探测器生成的位置谱进行有效识别。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明中记载的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
图1是根据本发明的一个实施例的识别位置谱的方法的步骤示意图；
图2是根据本发明的一个实施例的识别位置谱的方法的效果示意图，其中，初始位置谱(左)通过SiPM(硅光电倍增管)探测器获取，初始位置谱经过滤波处理后可以得到第一位置谱(右)；
图3是根据图2的识别位置谱的方法的第二位置谱的示意图，其中通过对第一位置谱进行特征提取后得到的其中一个第二位置谱；
图4是根据本发明一个实施例的识别位置谱的方法对图2所示的第二位置谱进行形态学膨胀后的位置谱示意图；
图5是根据图4的识别位置谱的方法的第四位置谱的示意图，其中通过对第二位置谱进行极值搜索后得到第三位置谱；
图6是根据本发明一个实施例的经过分割后的最终位置谱的示意图，其中初始位置谱通过SiPM探测器获取；
图7是根据本发明另一个实施例的识别位置谱的方法的效果示意图，其中，初始位置谱(左)通过PSPMT(位置灵敏型光电倍增管)探测器获取，初始位置谱经过滤波处理后可以得到第一位置谱(右)；
图8是根据图6的识别位置谱的方法的第二位置谱的示意图，其中通过对第一位置谱进行特征提取后得到其中一个第二位置谱；
图9是根据本发明一个实施例的识别位置谱的方法对图7所示的第二位置谱进行形态学膨胀后的位置谱示意图；
图10是根据本发明一个实施例的识别位置谱的方法对图8所示的位置谱进行极值搜索后得到的第三位置谱的示意图；
图11是根据本发明一个实施例的识别位置谱的方法的第四位置谱的示意图，其中通过对第三位置谱进行无效值剔除后得到第四位置谱，白色圆圈标记的是识别错误需要剔除的极值点；
图12是根据本发明一个实施例的识别位置谱的方法的分配冲突的示意图；
图13是根据图11的识别位置谱的方法的分配冲突后的结果示意图；
图14是根据本发明一个实施例的经过行列分配后的第五位置谱的示意图，其中初始位置谱通过SiPM探测器获取；
图15是根据本发明另一个实施例的经过行列分配后的第五位置谱示意图，其中初始位置谱通过PSPMT探测器获取；
图16是根据本发明一个实施例的识别位置谱的方法经过全局编号的第六位置谱的示意图；
图17是根据本发明另一个实施例的的经过分割后的最终位置谱的示意图，其中初始位置谱通过PSPMT探测器获取；
图18是根据本发明一个实施例的识别位置谱的装置的结构示意图。
具体实施方式
以下结合具体实施例，对本发明做进一步说明。应理解，以下实施例仅用于说明本发明而非用于限制本发明的范围。
需要说明的是，当部件/零件被称为“设置在”另一个部件/零件上，它可以直接设置在另一个部件/零件上或者也可以存在居中的部件/零件。当部件/零件被称为“连接/联接”至另一个部件/零件，它可以是直接连接/联接至另一个部件/零件或者可能同时存在居中部件/零件。本文所使用的术语“连接/联接”可以包括电气和/或机械物理连接/联接。本文所使用的术语“包括/包含”指特征、步骤或部件/零件的存在，但并不排除一个或更多个其它特征、步骤或部件/零件的存在或添加。本文所使用的术语“和/或”包括一个或多个相关所列项目的任意的和所有的组合。
除非另有定义，本文所使用的所有的技术和科学术语与属于本发明的技术领域的技术人员通常理解的含义相同。本文中所使用的术语只是为了描述具体实施例的目的，而并不是旨在限制本发明。
另外，在本发明的描述中，术语“第一”、“第二”等仅用于描述目的和区别类似的对象，两者之间并不存在先后顺序，也不能理解为指示或暗示相对重要性。此外，在本发明的描述中，除非另有说明，“多个”的含义是两个或两个以上。
如图1所示，本发明提供的识别位置谱的方法可以包括以下步骤：
步骤S1：对初始位置谱进行预处理以得到第一位置谱；
步骤S2：对第一位置谱进行特征提取，得到多个第二位置谱；
步骤S3：对第二位置谱进行极值搜索，得到N个极值点，这些极值点形成第三位置谱，其中，N为探测器中闪烁晶体的数量；
步骤S4：对第三位置谱中的极值点进行全局编号，形成第六位置谱；
步骤S5：根据第六位置谱中的极值点对单事件进行聚类，形成最终位置谱，完成位置谱的分割。
更具体地，在上述步骤S1之前，本领域技术人员可以通过本领域常用的设备或者方法获取初始位置谱，比如通过SiPM探测器和多 阈值采样方法获取初始位置谱，或者通过PSPMT探测器和多阈值采样方法方法获取初始位置谱，即，本申请中的初始位置谱可以是现有技术中通过任何方法和设备所获取的位置谱。本领域技术人员公知的是，位置谱是通过从PET探测器输出的单事件数据中提取的位置信息所形成的所有光子的二维位置分布图，图像中每一个像素值表示该二维坐标上落入光子的数目，图像中的光斑表示了不同晶体条所产生的光子的位置分布，像素值越大表示落入该位置的光子数量越多。本申请中的初始位置谱可以是边缘清晰、形变较小的位置谱，也可以是边缘模糊、形变严重的位置谱。在本领域中，位置谱的边缘的清晰度和形变的大小可以通过本领域技术人员的经验进行判断，其并不是本领域中难以确定的，在此不再赘述。
在上述步骤S1中，对初始位置谱进行预处理的具体方法可以采用二维高斯滤波，具体地，本领域技术人员可以根据实际经验或者随机选定第一高斯模板的尺寸大小和第一方差，然后将初始位置谱与第一高斯模板进行卷积处理，得到第一位置谱。第一高斯模板、第一高斯模板的尺寸大小和第一方差的选择以及卷积操作均属于本领域技术人员的常规技术选择，其并不是本领域中难以确定的，在此不再赘述。通过预处理可以使得初始位置谱中的每一个光斑变得更加均匀和明显，即初始位置谱中的光斑的形状和分布更加接近二维高斯函数，如图2所示。
在上述步骤S2中，由于第一位置谱的主要特征在于局部的光斑，为了能够在更广的范围内正确的识别光斑，需要对第一位置谱中的光斑进行特征提取，因此，可以采用基于求导的微分方法，针对光斑分布特征以判断光斑与微分算子的相似性。光斑的分布特征指光斑在二维坐标上的分布，可以是光斑的排布形状等，比如中间计数高、四周计数低，即光斑中间清晰、四周较弱等特征。特征提取的具体步骤可以包括：
步骤S21：确定一组第二方差σ 1,σ 2,…,σ n(其中n为自然数)，根据每一个第二方差σ n和相同的尺寸大小生成一个二维高斯模板，即生成n个第二高斯模板，记为g，同时根据探测器中闪烁晶体的横截 面建立与位置谱的图像对应的正交坐标系XOY；
步骤S22：求每一个第二高斯模板和第一位置谱卷积后的解，即求Hessian矩阵：
首先，根据不同的第二方差计算第二高斯模板关于x的二阶偏导：
Figure PCTCN2020099959-appb-000006
其次，根据不同的第二方差计算第二高斯模板关于y的二阶偏导：
Figure PCTCN2020099959-appb-000007
再次，根据不同的第二方差计算第二高斯模板关于x的偏导后再计算关于y的偏导：
Figure PCTCN2020099959-appb-000008
最后，对上述三种偏导进行规范化，即乘上σ 2，之后将三种偏导分别与第一位置谱进行卷积操作，得到n个L xx，L xy，L yy，从而可以获得n个Hessian矩阵，即
Figure PCTCN2020099959-appb-000009
步骤S23：计算n个Hessian矩阵的行列式，即
Figure PCTCN2020099959-appb-000010
得到n个第二位置谱，如图3所示。
在上述步骤S21中，本领域技术人员可以根据实际经验或者随机选定第二方差的具体值，这并不是本领域中难以确定的，在此不再赘述。
在上述步骤S3中，进行极值搜索之前可以对n个第二位置谱进行形态学膨胀操作得到n个膨胀位置谱，如图4所示；极值搜索即在n个膨胀位置谱中搜索最大的N个值，同时获得N个值对应的位置坐标和第二方差的值，从而以形成第三位置谱，N为PET探测器中对应的闪烁晶体的数目，N个最大值对应的位置即识别出的N个极 值点，每个极值点所在的膨胀位置谱对应着第二高斯模板方差的大小；通常取3倍的第二高斯模板方差作为每个晶体在第三位置谱上的半径，如图5所示。
对于排列整齐、分布清晰的位置谱，当经过步骤S1-S3的处理后，即可进行步骤S4。在步骤S4中，全局编号可以通过极值点邻域排序的方式进行，具体步骤可以包括：
步骤S41：获取所有极值点的位置坐标和位置响应半径，其中，极值点的位置坐标与原始建立的XOY坐标系保持一致，位置响应半径可以根据实际晶体条的横截面尺寸选取，比如通常取3倍的第二高斯模板方差作为每个极值点的位置响应半径；
步骤S42：提取极值点的位置坐标在位置响应半径范围内的区域，若极值点个数为N，则提取的区域个数也为N；
步骤S43：计算每个区域内在行/列方向上对应坐标的最小值；具体地，比如对行进行判断时，需要判断整个区域内所有点中最小的y值；对列进行判断时，需要判断整个区域内所有点中最小的x值；
步骤S44：对N个最小值进行排序，则第m×(n-1)+1到第m×n个最小值对应的极值点隶属于第n行/列，m为行/列上的晶体数量，n为自然数；
步骤S45：根据极值点所分配的行/列按晶体阵列实际的几何顺序进行编号，将编号显示在极值点周围以形成第六位置谱，如图6所示，其中的数字为经过全局编号的结果。
在上述步骤S5中，事件聚类是指在获得位置谱中极值点的位置坐标后，需要对所有单事件按照位置信息进行聚类，生成晶体对照表(crystal lookup table)。由于正电子与电子湮灭后产生一对能量相同、方向相反的光子，每一个被探测到的光子事件都包含了湮灭的位置、时间和能量信息，为了准确的判断位置信息，就需要对所有单事件进行聚类，判断每一个光子事件来源于哪一根晶体条，将每一个光子事件与位置谱的光斑对应起来。具体地，事件聚类可以包括以下步骤：
步骤S51：任意选取一个单事件的位置坐标，计算该位置坐标与N个所识别出的极值点的位置坐标之间的距离，N为晶体数量；
步骤S52：在步骤S51的N个距离中选取最短的距离所在的晶体作为步骤S51中位置坐标的晶体聚类；
步骤S53：遍历所有单事件的位置坐标，获得所有位置坐标的晶体聚类，从而形成最终位置谱，如图6所示。
此外，在上述步骤S5中，还可以采用已识别晶体的极值点位置来初始化GMM(高斯混合模型)中的高斯函数均值、初始化SOM(自组织映射算法)中的初始神经元权值、初始化KNN(K-近邻算法)算法中的聚集点中心以高效率地获得单事件的晶体聚类，具体步骤不再赘述。
在上述步骤S51中，单事件的位置坐标本质上是位置谱的像素编号。光斑代表着晶体，比如对于13×13排列的晶体条，其总共的像素值为256×256，表示光子可能会打中的256×256个位置，因此，对于256×256像素大小的位置谱，对于第一行第一列的单事件位置，其位置坐标为(1，1)，对于第i行第j列的单事件位置，其位置坐标为(i，j)。
在上述实施例中，步骤S1-步骤S5能够以更高的准确度、更快的速度、更高的效率处理、识别位置谱，尤其适用于排列整齐、分布位置清晰的位置谱，比如图2所示的通过SiPM(硅光电倍增管)探测器获取的初始位置谱。但是，对于对于排列杂乱、分布模糊的位置谱，比如图7所示的通过PSPMT(位置灵敏型光电倍增管)探测器获取的初始位置谱，按照上述步骤S1-S3进行处理之后，可能会得到较多错误的极值点，因此，需要在极值搜索之后进行无效值剔除、行列分配和缺失点预测以提高对模糊分布的晶体进行识别和预测的效果。
根据本发明的另一个实施例，本发明提供的识别位置谱的方法的步骤还可以包括以下步骤：
步骤S1：对初始位置谱进行预处理以得到第一位置谱；
步骤S2：对第一位置谱进行特征提取，得到多个第二位置谱；
步骤S3：对第二位置谱进行极值搜索，得到N个极值点，N为探测器中闪烁晶体的数量；其中，在极值搜索后需要进行步骤S31- 步骤S33；
步骤S31：对第三位置谱中无效的极值点进行剔除，得到第四位置谱；
步骤S32：对剔除后的所有极值点进行行列分配，形成第五位置谱；
步骤S33：对缺失的极值点进行预测以形成第三位置谱；
步骤S4：对第三位置谱中的极值点进行全局编号，形成第六位置谱；
步骤S5：根据对第六位置谱中的极值点对单事件进行聚类，形成最终位置谱，完成位置谱的分割。
上述步骤S1和步骤S2中的具体方法与上一实施例中相似或者相同，在此不再赘述，比如，对于通过PSPMT(位置灵敏型光电倍增管)探测器获取的初始位置谱，其排列杂乱、分布模糊，经过预处理后的效果如图7所示；进行特征提取后的位置谱如图8所示，进行形态学膨胀和极值搜索后的效果分别如图9、图10所示。
在上述步骤S31中，对无效的极值点进行剔除的具体步骤可以包括：
步骤S311：根据极值搜索所识别出的极值点的位置坐标和第二方差大小提取出所有极值点的位置分布区域图像，如图10所示；
步骤S312：确定一个初始分割阈值T，将位置分布区域图像中小于初始分割阈值T的像素点的数目占所有像素点数目的比例记为w 0，平均灰度记为u 0；将位置分布区域图像中大于初始分割阈值T的像素点的数目占所有像素点数目的比例记为w 1，平均灰度记为u 1；将整体位置分布区域图像的平均灰度记为u；
步骤S313：按照公式k＝w 0×w 1×(u 0-u 1) 2计算类间方差k；
步骤S314：遍历初始分割阈值T的值，使得类间方差k最大，从而获得晶体的位置分布区域图像的二值化图像；
步骤S315：取得二值化图像中心一个邻域并计算像素和，若该邻域内的像素和的值为0，则认定该极值点为错误识别的极值点，即无效的极值点，如图11所示，白色圆圈标记的是识别错误的极值点。
在上述步骤S32中，行列分配的具体步骤可以按照晶体排列的行或者列进行，在此以行分配举例，行列分配的具体步骤可以包括：
步骤S321：将步骤S31中识别的极值点的位置坐标按照对应的维度坐标的大小分别进行排序；
步骤S322：若按照x坐标排序，则抽取第m×(n-1)+1到第m×n个极值点作为第n行的模糊行集，m为每一行上的晶体数量，n为自然数；
步骤S323：从第n行的模糊行集中随机抽取p个极值点作为参与拟合的点集S，剩余极值点的集合记为余集R，p为拟合曲线所需的样本数；
步骤S324：对点集S中的极值点进行基于最小二乘法的二次曲线拟合；
步骤S325：对余集R中的所有极值点求取与拟合的二次曲线的误差和E；
步骤S326：重复步骤S323-步骤S325的过程，直至在迭代次数内获得的误差和E的最小值或者误差和E小于设定的最小容忍阈值，此时拟合的二次曲线即为第n行极值点的分布曲线；
步骤S327：将已经识别的极值点按照行列分配原则分配至每一行/列的分布曲线上。
在上述步骤S323-步骤S325中，基于最小二乘法进行曲线拟合是本领域技术人员容易实现的，在此不再赘述。
在上述步骤S326中，误差和E的最小容忍阈值可以通过多次试验确定，本领域技术人员通过多次的试验结果，可以很容易的确定最小容忍阈值的范围，在此不再赘述。
在上述步骤S327中，行列分配原则的本质是将位置谱上的极值点分配给行、列分布曲线的交点，行列分配原则包括就近分配原则和冲突分配原则。
对于就近分配原则，包括：第一、按照极值点和分布曲线交点之间的距离最近分配；第二、一个分布曲线的交点处最多占据一个极值点，一个极值点仅能分配给一个交点。
对于冲突分配原则，具体地，对于一个交点对应多个极值点的情况，如图12所示，空心圆表示极值点，实心圆表示分布曲线的交点，当多个极值点与某个分布曲线的交点的距离相同时，该交点便会对应多个极值点。如图12中l A1＝l B1，则表示分布曲线的交点c同时对应极值点A和B，l为极值点与分布曲线交点之间的距离。此时，冲突分配原则按照以下步骤进行：
第一、通过分布曲线的交点找到对应的多个极值点，比如，交点c对应的两个极值点A和B；
第二、对每个对应的极值点分别寻找最近的两个交点，比如，极值点A最近的两个交点分别为a和c，极值点B最近的两个交点分别为c和d；
第三、判断每个极值点距离两个最近的交点的距离之和的大小，距离之和最大的极值点分配至有冲突的分布曲线的交点，比如，图12中极值点A与交点c、a之间的距离分别为l A1、l A2，极值点B与交点c、d之间的距离分别为l B1、l B2，极值点A与两个最近的交点c、a之间的距离之和l A1+l A2小于极值点B与两个最近的交点c、d之间的距离之和l B1+l B2，因此极值点B分配至交点c，极值点A分配至与其距离最近的另外一个交点a，分配结果如图13所示。
对于通过PSPMT(位置灵敏型光电倍增管)探测器获取的初始位置谱，其经过极值搜索和无效值剔除后的位置谱如图11所示，经过行列分配之前的位置谱如图14所示，可以看到其中一些极值点共享了相同的分布曲线交点；经过行列分配之后的位置谱如图15所示，可以看到经过行列分配之后所有的极值点和分布曲线的交点一一对应。
进一步地，由于经过错误识别和剔除以后，最终识别出的晶体的数目或者极值点的数量要少于目标数目，因此需要进一步对缺失晶体的位置进行预测，即进行步骤S33。在上述步骤S33中，对缺失的极值点进行预测是根据分布曲线交点行列分配的结果进行的，经过行列分配之后，所有识别的极值点都有且只有一个分布曲线的交点与之对应，缺失点预测即找到未分配极值点的交点坐标，并且将该交点坐标 赋予给当前行和当前列上缺失的极值点。
本领域技术人员需要注意的是，对于排列杂乱、分布模糊的位置谱，由于其需要进行无效值剔除、行列分配以及缺失点预测的步骤，步骤S4可以根据行列分配后的结果直接进行，在此不再赘述。
虽然本发明提供了如上述实施例或流程图所述的方法操作步骤，但基于常规或者无需创造性的劳动在所述方法中可以包括更多或者更少的操作步骤。在逻辑性上不存在必要因果关系的步骤中，这些步骤的执行顺序不限于本发明实施例提供的执行顺序。
通过本发明提供的位置谱分割的方法，可以根据需要对任意的位置谱进行处理，也可以根据需要对位置谱进行划分后处理，比如对于边缘模糊、性变严重的位置谱，可以直接按照步骤S1-步骤S2-步骤S3-步骤S31、步骤S32、步骤S33-步骤S4-步骤S5进行分割，对于排列整齐、分布清晰的位置谱，可以按照步骤S1-步骤S2-步骤S3-步骤S4-步骤S5进行分割。总之，通过本发明提供的方法，可以高效率、高准确度地实现对边缘模糊、形变严重的位置谱的识别，同时支持对不同结构探测器生成的位置谱进行有效识别。
进一步地，如图18所示，本发明还提供了一种基于上述方法的识别位置谱的装置，该装置包括预处理模块10、特征提取模块20、极值搜索模块30、全局编号模块40和事件聚类模块50，其中，预处理模块10接收通过探测器输出的初始位置谱并对初始位置谱进行预处理以形成第一位置谱；特征提取模块20接收第一位置谱并对第一位置谱进行特征提取以得到第二位置谱；极值搜索模块30接收第二位置谱并对第二位置谱进行极值搜索以得到第三位置谱；全局编号模块40对第三位置谱进行全局编号以形成第六位置谱；事件聚类模块50对第六位置谱中的单事件进行聚类以形成最终位置谱，完成位置谱的分割。
初始位置谱可以通过本领域常用的设备或者方法获取，比如通过SiPM探测器和多阈值采样方法获取初始位置谱，或者通过PSPMT探测器和传统的时间间隔采样方法获取初始位置谱。
更具体地，预处理模块10可以采用二维高斯滤波对初始位置谱 进行预处理，预处理模块10根据本领域技术人员通过实际经验或者随机设定的第一高斯模板的尺寸大小和第一方差，然后将初始位置谱与第一高斯模板进行卷积处理，得到第一位置谱。通过预处理可以使得初始位置谱中的每一个光斑变得更加均匀和明显，即初始位置谱中的光斑的形状和分布更加接近二维高斯函数，如图2所示。
进一步地，为了能够在更广的范围内正确的识别光斑，需要对第一位置谱中的光斑进行特征提取，特征提取模块20可以采用基于求导的微分方法，针对光斑分布特征以判断光斑与微分算子的相似性。光斑的分布特征指光斑在二维坐标上的分布，可以是光斑的排布形状等，比如中间计数高、四周计数低，即中间光斑清晰、四周光斑较弱等特征。特征提取模块20可以按照下述具体步骤进行特征提取：
步骤S21：预设一组第二方差σ 1,σ 2,…,σ n(其中n为自然数)，根据每一个第二方差σ n和相同的尺寸大小生成一个二维高斯模板，即生成n个第二高斯模板，记为g；
步骤S22：求每一个第二高斯模板和第一位置谱卷积后的解，即求Hessian矩阵：
首先，根据不同的第二方差计算第二高斯模板关于x的二阶偏导：
Figure PCTCN2020099959-appb-000011
其次，根据不同的第二方差计算第二高斯模板关于y的二阶偏导：
Figure PCTCN2020099959-appb-000012
再次，根据不同的第二方差计算第二高斯模板关于x的偏导后再计算关于y的偏导：
Figure PCTCN2020099959-appb-000013
最后，对上述三种偏导进行规范化，即乘上σ 2，之后将三种偏导分别与第一位置谱进行卷积操作，得到n个L xx，L xy，L yy，从而可以 获得n个Hessian矩阵，即
Figure PCTCN2020099959-appb-000014
步骤S23：计算n个Hessian矩阵的行列式，即
Figure PCTCN2020099959-appb-000015
得到n个第二位置谱，如图3所示。
极值搜索模块30可以对n个第二位置谱进行形态学膨胀操作得到n个膨胀位置谱，然后极值搜索模块30在n个膨胀位置谱中搜索最大的N个值作为极值点以得到第三位置谱，N为PET探测器中对应的闪烁晶体的数目。
对于排列整齐、分布清晰的位置谱，当经过预处理模块10、特征提取模块20、极值搜索模块30的处理后，即可通过全局编号模块40进行全局编号，全局编号可以通过极值点邻域排序的方式进行，具体步骤可以包括：
步骤S41：获取所有极值点的位置坐标和位置响应半径，其中，极值点的位置坐标与原始建立的XOY坐标系保持一致，位置响应半径可以根据实际晶体条的横截面尺寸选取，比如通常取3倍的第二高斯模板方差作为每个极值点的位置响应半径；
步骤S42：提取极值点的位置坐标在位置响应半径范围内的区域，若极值点个数为N，则提取的区域个数也为N；
步骤S43：计算每个区域内在行/列方向上对应坐标的最小值；具体地，比如对行进行判断时，需要判断整个区域内所有点中最小的y值；对列进行判断时，需要判断整个区域内所有点中最小的x值；
步骤S44：对N个最小值进行排序，则第m×(n-1)+1到第m×n个最小值对应的极值点隶属于第n行/列，m为行/列上的晶体数量，n为自然数；
步骤S45：根据极值点所分配的行/列按晶体阵列实际的几何顺序进行编号，将编号显示在极值点周围以形成第六位置谱，如图6所示，其中的数字为经过全局编号的结果。
最后，事件聚类模块50对所有单事件按照位置信息进行聚类， 生成晶体对照表(crystal lookup table)。事件聚类模块50可以按照以下步骤处理：
步骤S51：任意选取一个单事件的位置坐标，计算该位置坐标与N个所识别出的极值点的位置坐标之间的距离，N为晶体数量；
步骤S52：在步骤S51的N个距离中选取最短的距离所在的晶体作为步骤S51中位置坐标的晶体聚类；
步骤S53：遍历所有单事件的位置坐标，获得所有位置坐标的晶体聚类，从而形成最终位置谱，如图6所示。
对于对于排列杂乱、分布模糊的位置谱，比如图7所示的通过PSPMT(位置灵敏型光电倍增管)探测器获取的初始位置谱，经过预处理模块10、特征提取模块20、极值搜索模块30的处理后，可能会得到较多错误的极值点，因此，需要在极值搜索模块之后采用无效值剔除模块60、行列分配模块70和缺失点预测模块80处理以提高对模糊分布的晶体进行识别和预测的效果。无效值剔除模块60、行列分配模块70和缺失点预测模块80的处理步骤与上述第二实施例中的无效值剔除步骤、行列分配步骤和缺失点预测步骤相同，在此不再赘述。无效值剔除模块60针对错误识别区域的特征来进行适应性的错误筛除工作，可以解决分布模糊的位置谱经过特征提取会得到错误的极值点的问题，能够有效剔除大部分识别错误的极值点位置。行列分配模块70可以使所有识别的极值点都有且只有一个拟合曲线的交点与之对应。缺失点预测模块80可以获得未分配极值点的交点坐标，同时将该交点坐标赋予给当前行和当前列上缺失的极值点。
通过本发明提供的位置谱分割的装置，可以根据需要对任意的位置谱进行处理，也可以根据需要对位置谱进行划分后处理，比如，对于排列整齐、分布清晰的位置谱，可以直接通过预处理模块、特征提取模块、极值搜索模块、全局编号模块和事件聚类模块进行分割；对于边缘模糊、形变严重的位置谱，可以额外采用无效值剔除模块、行列分配模块和缺失点预测模块进行分割。总之，通过本发明提供的装置，可以高效率、高准确度地实现对边缘模糊、形变严重的位置谱的识别，同时支持对不同结构探测器生成的位置谱进行有效识别。
本发明还提供了一种计算机存储介质，该计算机存储介质存储有计算机程序，该计算机程序被处理器执行时可以实现上述方法实施例中的步骤，比如：对初始位置谱进行预处理以得到第一位置谱；对第一位置谱进行特征提取，得到多个第二位置谱；对第二位置谱进行极值搜索，得到N个极值点，通过这些极值点形成第三位置谱，其中，N为探测器中闪烁晶体的数量；对第三位置谱中的极值点进行全局编号，形成第六位置谱；根据第六位置谱中的极值点对单事件进行聚类，形成最终位置谱，完成位置谱的分割。另外，该计算机程序被执行时也可以实现对第三位置谱中无效的极值点进行剔除以得到第四位置谱以及对剔除后的所有极值点进行行列分配以形成第五位置谱。
关于该实施例的详细描述，可以参照上述实施例中对位置谱的分割方法的详细描述，在此不再赘叙。
本发明的计算机存储介质可以包括能够携带计算机程序代码的任何实体或装置，比如ROM/RAM、磁盘、光盘、闪存等存储器，在此不再一一列举。
上述实施例阐明的装置、模块等，具体可以由计算机芯片和/或实体实现，或者由具有某种功能的产品来实现。为了描述的方便，描述以上装置时以功能分为各种模块分别描述。当然，在实施本申发明的实施例时可以把各模块的功能集成在同一个或多个计算机芯片中实现。
以上所述的，仅为本发明的较佳实施例，并非用以限定本发明的范围，本发明的上述实施例还可以做出各种变化。即凡是依据本发明申请的权利要求书及说明书内容所作的简单、等效变化与修饰，皆落入本发明专利的权利要求保护范围。本发明未详尽描述的均为常规技术内容。

Claims (19)

  1. 一种识别位置谱的方法，所述位置谱为光子的二维位置分布图，所述位置谱包含从探测器输出的单事件中提取的所有光子的位置信息，其特征在于，所述方法包括以下步骤：
    步骤S1：对初始位置谱进行预处理以得到第一位置谱；
    步骤S2：对所述第一位置谱进行特征提取，得到多个第二位置谱；
    步骤S3：对所述第二位置谱进行极值搜索，得到N个极值点，所述极值点形成第三位置谱，其中，N为所述探测器中闪烁晶体的数量；
    步骤S4：对所述第三位置谱中的所述极值点进行全局编号，形成第六位置谱；
    步骤S5：根据所述第六位置谱中的所述极值点对所述单事件进行聚类，形成最终位置谱，完成位置谱的分割。
  2. 根据权利要求1所述的识别位置谱的方法，其特征在于，在所述步骤S1中，所述预处理的具体步骤为：选取第一高斯模板的尺寸大小和第一方差，将所述初始位置谱与所述第一高斯模板进行卷积处理以得到第一位置谱。
  3. 根据权利要求1所述的识别位置谱的方法，其特征在于，在所述步骤S2中，特征提取的具体步骤包括：
    步骤S21：确定一组第二方差σ 1,σ 2,…,σ n，n为自然数，根据每一个所述第二方差σ n和相同的尺寸大小生成n个第二高斯模板，记为g；
    步骤S22：求每一个所述第二高斯模板与所述第一位置谱卷积后的Hessian矩阵，
    Figure PCTCN2020099959-appb-100001
    其中，x、y为所述位置谱的图像中对应的正交坐标轴；
    步骤S23：计算n个Hessian矩阵的行列式：
    Figure PCTCN2020099959-appb-100002
    得到n个第二位置谱。
  4. 根据权利要求3所述的识别位置谱的方法，其特征在于，在所述步骤S22中，求Hessian矩阵的具体步骤包括：
    第一步，根据不同的所述第二方差计算所述第二高斯模板关于x的二阶偏导：
    Figure PCTCN2020099959-appb-100003
    第二步，根据不同的所述第二方差计算所述第二高斯模板关于y的二阶偏导：
    Figure PCTCN2020099959-appb-100004
    第三步，根据不同的所述第二方差计算所述第二高斯模板关于x的偏导后再计算关于y的偏导：
    Figure PCTCN2020099959-appb-100005
    第四步，对上述三种偏导进行规范化，即乘上σ 2，之后将三种偏导分别与所述第一位置谱进行卷积操作，得到n个L xx，L xy，L yy，从而获得n个Hessian矩阵。
  5. 根据权利要求1所述的识别位置谱的方法，其特征在于，在上述步骤S3中，进行极值搜索之前需要对多个所述第二位置谱进行形态学膨胀处理以得到膨胀位置谱。
  6. 根据权利要求5所述的识别位置谱的方法，其特征在于，在上述步骤S3中，极值搜索的具体步骤是在多个所述膨胀位置谱中搜索最大的N个值以得到所述第三位置谱。
  7. 根据权利要求1所述的识别位置谱的方法，其特征在于，在所述步骤S4中，全局编号的具体步骤包括：
    步骤S41：获取所有所述极值点的位置坐标和位置响应半径；
    步骤S42：提取所述极值点的所述位置坐标在所述位置响应半径范围内的区域；
    步骤S43：计算每个所述区域内在行/列方向上的相应维度坐标的最小值；
    步骤S44：对N个最小值进行排序，第m×(n-1)+1到第m×n个最小值 对应的极值点隶属于第n行/列，m为某一行/列上的闪烁晶体数量，n为自然数；
    步骤S45：根据所述极值点所分配的行/列按闪烁晶体阵列实际的几何顺序进行编号，将编号显示在极值点周围以形成所述第六位置谱。
  8. 根据权利要求7所述的识别位置谱的方法，其特征在于，在所述步骤S5中，事件聚类是指获得位置谱中极值点的所述位置坐标后，对所有单事件按照位置信息进行聚类，生成晶体对照表并获得位置谱的分割结果。
  9. 根据权利要求8所述的识别位置谱的方法，其特征在于，事件聚类的具体步骤包括：
    步骤S51：任意选取一个所述单事件的位置坐标，计算所述位置坐标与N个所识别出的所述极值点的所述位置坐标之间的距离；
    步骤S52：在所述步骤S51的N个距离中选取最短的距离所在的极值点作为所述步骤S51中所述位置坐标的晶体聚类；
    步骤S53：遍历所有的所述单事件的位置坐标，获得所有所述位置坐标对应的晶体聚类，将不同类别用折线区分开来以形成所述最终位置谱。
  10. 根据权利要求3所述的识别位置谱的方法，其特征在于，所述步骤S3中还包括以下步骤：
    步骤S31：对所述第三位置谱中无效的极值点进行剔除以形成第四位置谱；
    步骤S32：对剔除后的所有极值点进行行列分配以形成第五位置谱；
    步骤S33：对所述第五位置谱中缺失的极值点进行预测以形成所述第三位置谱。
  11. 根据权利要求10所述的识别位置谱的方法，其特征在于，在所述步骤S31中，对所述第三位置谱中无效的极值点进行剔除的具体步骤包括：
    步骤S311：根据所述极值点的位置坐标和第二方差大小提取出所有极值点的位置分布区域图像；
    步骤S312：确定一个初始分割阈值T，将位置分布区域图像中小于初始分割阈值T的像素点的数目占所有像素点数目的比例记为w 0，平均灰度记为u 0；将位置分布区域图像中大于初始分割阈值T的像素点的数目占所有像 素点数目的比例记为w 1，平均灰度记为u 1；
    步骤S313：按照公式k＝w 0×w 1×(u 0-u 1) 2计算类间方差k；
    步骤S314：遍历初始分割阈值T的值，使得类间方差k最大，从而获得闪烁晶体所产生单事件的位置分布区域图像的二值化图像；
    步骤S315：取得二值化图像中心一个邻域并计算像素和，若该邻域内的像素和的值为0，则认定该极值点为无效的极值点。
  12. 根据权利要求10所述的识别位置谱的方法，其特征在于，在所述步骤S32中，行列分配的具体步骤可以包括：
    步骤S321：将步骤S31中识别的极值点的位置坐标按照对应的维度坐标的大小分别进行排序；
    步骤S322：抽取第m×(n-1)+1到第m×n个极值点作为第n行/列的模糊行集，m为每一行上的晶体数量，n为自然数；
    步骤S323：从第n行/列的模糊行集中随机抽取p个极值点作为参与拟合的点集S，剩余极值点的集合记为余集R，p为拟合曲线所需的样本数；
    步骤S324：对点集S中的极值点进行二次曲线拟合；
    步骤S325：对余集R中的所有极值点求取与拟合的二次曲线的误差和E；
    步骤S326：重复步骤S323-步骤S325的过程，直至在迭代次数内获得的误差和E的最小值或者误差和E小于设定的最小容忍阈值，得到第n行/列极值点的分布曲线；
    步骤S327：将已经识别的极值点按照行列分配原则分配至每一行/列的分布曲线上。
  13. 根据权利要求10所述的识别位置谱的方法，其特征在于，所述行列分配原则包括就近分配原则和冲突分配原则。
  14. 根据权利要求13所述的识别位置谱的方法，其特征在于，所述就近分配原则包括：按照极值点和分布曲线交点之间的距离最近分配；一个分布曲线的交点处最多占据一个极值点，一个极值点仅能分配给一个交点。
  15. 根据权利要求13所述的识别位置谱的方法，其特征在于，所述冲突分配原则包括：通过分布曲线的交点找到对应的多个极值点；对每个对应的极值点分别寻找最近的两个交点；判断每个极值点距离两个最近的交点的距 离之和的大小，距离之和最大的极值点分配至有冲突的分布曲线的交点。
  16. 根据权利要求15所述的识别位置谱的方法，其特征在于，在所述步骤S33中，对缺失的极值点进行预测的具体步骤是：找到未分配极值点的交点坐标，并且将该交点的位置坐标赋予给当前行和当前列上缺失的极值点。
  17. 一种识别位置谱的装置，其特征在于，所述装置包括：
    预处理模块，所述预处理模块接收通过探测器输出的初始位置谱并对所述初始位置谱进行预处理以形成第一位置谱；
    特征提取模块，所述特征提取模块接收所述第一位置谱并对所述第一位置谱进行特征提取以得到第二位置谱；
    极值搜索模块，所述极值搜索模块接收所述第二位置谱并对所述第二位置谱进行极值搜索以得到第三位置谱；
    全局编号模块，所述全局编号模块对所述第三位置谱进行全局编号以形成第六位置谱；以及
    事件聚类模块，所述事件聚类模块对所述第六位置谱中的单事件进行聚类以形成最终位置谱，完成位置谱的分割。
  18. 根据权利要求17所述的识别位置谱的装置，其特征在于，所述装置进一步包括：
    无效值剔除模块，所述无效值剔除模块剔除所述第三位置谱中无效的极值点并形成第四位置谱；
    行列分配模块，所述行列分配模块对所述第四位置谱中的极值点进行分配以形成第五位置谱；以及
    缺失点预测模块，所述缺失点预测模块对所述第五位置谱中缺失的极值点进行预测。
  19. 一种计算机存储介质，其特征在于，所述计算机存储介质上存储有计算机程序，所述计算机程序被执行时实现如权利要求1-16任一项所述方法的步骤。
PCT/CN2020/099959 2019-08-19 2020-07-02 一种识别位置谱的方法、装置以及计算机存储介质 Ceased WO2021031711A1 (zh)

Priority Applications (4)

Application Number Priority Date Filing Date Title
US17/597,642 US12164071B2 (en) 2019-08-19 2020-07-02 Method and apparatus for identifying location spectrum, and computer storage medium
JP2021563135A JP7198532B2 (ja) 2019-08-19 2020-07-02 位置スペクトルを識別するための方法、デバイス、及びコンピュータ記憶媒体
FIEP20855620.9T FI4020020T3 (fi) 2019-08-19 2020-07-02 Menetelmä ja laite paikannusspektrin tunnistamiseksi ja tietokoneen tallennusväline
EP20855620.9A EP4020020B1 (en) 2019-08-19 2020-07-02 Method and apparatus for identifying location spectrum, and computer storage medium

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
CN201910764960.8 2019-08-19
CN201910764960.8A CN110471102B (zh) 2019-08-19 2019-08-19 一种识别位置谱的方法、装置以及计算机存储介质

Publications (1)

Publication Number Publication Date
WO2021031711A1 true WO2021031711A1 (zh) 2021-02-25

Family

ID=68511128

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/CN2020/099959 Ceased WO2021031711A1 (zh) 2019-08-19 2020-07-02 一种识别位置谱的方法、装置以及计算机存储介质

Country Status (6)

Country Link
US (1) US12164071B2 (zh)
EP (1) EP4020020B1 (zh)
JP (1) JP7198532B2 (zh)
CN (1) CN110471102B (zh)
FI (1) FI4020020T3 (zh)
WO (1) WO2021031711A1 (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116989664A (zh) * 2023-09-27 2023-11-03 板石智能科技(深圳)有限公司 一种基于光谱共焦位移传感器的谱峰峰值计算方法及系统
CN118799591A (zh) * 2024-09-10 2024-10-18 山东科技大学 基于点云数据的露天采场坡面顶底线自动提取方法

Families Citing this family (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110471102B (zh) 2019-08-19 2021-06-01 苏州瑞派宁科技有限公司 一种识别位置谱的方法、装置以及计算机存储介质
CN113554673B (zh) * 2020-04-26 2024-08-20 中国石油化工股份有限公司 一种基于随钻电成像图像的自动识别裂缝的方法及系统
CN111938684B (zh) * 2020-07-08 2024-01-02 南昌大学 高速低时空复杂度的pet晶体位置谱数据分割方法
JP7509015B2 (ja) * 2020-12-04 2024-07-02 株式会社島津製作所 表面分析装置
CN112690818B (zh) * 2020-12-29 2022-06-17 赛诺联合医疗科技（北京）有限公司 Pet探测器的晶体位置查找表的校正方法及pet系统
JP7570374B2 (ja) * 2022-06-14 2024-10-21 アンリツ株式会社 イベント検出装置及びイベント検出方法
AU2024259579A1 (en) 2023-04-21 2025-11-06 Chevron U.S.A. Inc. Membrane preconcentration of carbon dioxide from exhaust gas sources
CN118761940B (zh) * 2024-09-04 2024-11-29 南京大学 一种卫星下传数据的快速解析方法

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101996325A (zh) * 2010-09-08 2011-03-30 北京航空航天大学 一种改进的提取图像中特征点的方法
CN103914860A (zh) * 2013-01-05 2014-07-09 苏州瑞派宁科技有限公司 晶体条位置查找表生成方法及装置
WO2015010393A1 (zh) * 2013-07-25 2015-01-29 苏州瑞派宁科技有限公司 针对全数字pet系统的在线能量符合方法及系统
CN104732541A (zh) * 2015-03-27 2015-06-24 北京永新医疗设备有限公司 Pet晶体位置查找表的获取方法和装置
CN106526651A (zh) * 2015-09-14 2017-03-22 上海联影医疗科技有限公司 一种探测器晶体位置表的建立方法及系统
CN110471102A (zh) * 2019-08-19 2019-11-19 苏州瑞派宁科技有限公司 一种识别位置谱的方法、装置及计算机存储介质

Family Cites Families (15)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6169285B1 (en) * 1998-10-23 2001-01-02 Adac Laboratories Radiation-based imaging system employing virtual light-responsive elements
US7026623B2 (en) * 2004-01-07 2006-04-11 Jacob Oaknin Efficient single photon emission imaging
US8586932B2 (en) * 2004-11-09 2013-11-19 Spectrum Dynamics Llc System and method for radioactive emission measurement
US9316743B2 (en) * 2004-11-09 2016-04-19 Biosensors International Group, Ltd. System and method for radioactive emission measurement
US8606349B2 (en) * 2004-11-09 2013-12-10 Biosensors International Group, Ltd. Radioimaging using low dose isotope
US8452064B2 (en) * 2009-12-07 2013-05-28 General Electric Company Apparatus and methods for geometric calibration of positron emission tomography systems
CN103177252B (zh) * 2013-03-04 2017-04-12 苏州瑞派宁科技有限公司 一种自动识别并分割位置谱的方法及装置
US9364192B2 (en) * 2013-03-15 2016-06-14 Siemens Medical Solutions Usa, Inc. Error estimates in quantitative functional imaging
US9928437B2 (en) * 2015-04-29 2018-03-27 Shanghai United Imaging Healthcare Co., Ltd. Method and system for crystal identification
CN106725560B (zh) * 2015-11-19 2021-01-12 上海联影医疗科技股份有限公司 光传感器的性能检测方法和医学成像设备
CN106295639A (zh) * 2016-08-01 2017-01-04 乐视控股（北京）有限公司 一种虚拟现实终端以及目标图像的提取方法和装置
CN107080551B (zh) * 2017-05-25 2023-08-22 苏州瑞派宁科技有限公司 一种三维异质pet系统
WO2018218441A1 (zh) * 2017-05-27 2018-12-06 北京悦琦创通科技有限公司 谱图分析方法、装置和设备及计算机可读存储介质
CN107238854A (zh) * 2017-07-25 2017-10-10 苏州瑞派宁科技有限公司 一种数字pet探测器的增益校正装置
WO2019084778A1 (en) * 2017-10-31 2019-05-09 Shenzhen United Imaging Healthcare Co., Ltd. System and method for pet correction

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101996325A (zh) * 2010-09-08 2011-03-30 北京航空航天大学 一种改进的提取图像中特征点的方法
CN103914860A (zh) * 2013-01-05 2014-07-09 苏州瑞派宁科技有限公司 晶体条位置查找表生成方法及装置
WO2015010393A1 (zh) * 2013-07-25 2015-01-29 苏州瑞派宁科技有限公司 针对全数字pet系统的在线能量符合方法及系统
CN104732541A (zh) * 2015-03-27 2015-06-24 北京永新医疗设备有限公司 Pet晶体位置查找表的获取方法和装置
CN106526651A (zh) * 2015-09-14 2017-03-22 上海联影医疗科技有限公司 一种探测器晶体位置表的建立方法及系统
CN110471102A (zh) * 2019-08-19 2019-11-19 苏州瑞派宁科技有限公司 一种识别位置谱的方法、装置及计算机存储介质

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
CHEN XIN , CHEN YUANBAO , NIU MING , WU ZHONGYI , WANG YUMBO , XIE QINGGUO: "Fast online crystal identification algorithm in positron emission tomography", CHINA SCIENCEPAPER, vol. 8, no. 4, 15 April 2013 (2013-04-15), pages 282 - 286, XP055782082, ISSN: 2095-2783 *
See also references of EP4020020A4 *

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116989664A (zh) * 2023-09-27 2023-11-03 板石智能科技(深圳)有限公司 一种基于光谱共焦位移传感器的谱峰峰值计算方法及系统
CN116989664B (zh) * 2023-09-27 2023-12-19 板石智能科技(深圳)有限公司 一种基于光谱共焦位移传感器的谱峰峰值计算方法及系统
CN118799591A (zh) * 2024-09-10 2024-10-18 山东科技大学 基于点云数据的露天采场坡面顶底线自动提取方法

Also Published As

Publication number Publication date
JP7198532B2 (ja) 2023-01-04
US12164071B2 (en) 2024-12-10
EP4020020A4 (en) 2023-10-11
EP4020020B1 (en) 2025-11-26
JP2022534357A (ja) 2022-07-29
FI4020020T3 (fi) 2026-02-23
EP4020020A1 (en) 2022-06-29
US20220163686A1 (en) 2022-05-26
CN110471102B (zh) 2021-06-01
CN110471102A (zh) 2019-11-19

Similar Documents

Publication Publication Date Title
WO2021031711A1 (zh) 一种识别位置谱的方法、装置以及计算机存储介质
Cao et al. A novel attention-guided convolutional network for the detection of abnormal cervical cells in cervical cancer screening
CN114930353A (zh) 受引导的缺陷发现的特性化系统及方法
US8787620B2 (en) Automated crystal identification achieved via watershed segmentation
CN110770752A (zh) 多尺度特征融合网络结合定位模型的害虫自动计数方法
Fang et al. An automatic method for counting wheat tiller number in the field with terrestrial LiDAR
WO2016090148A1 (en) Analysis and screening of cell secretion profiles
CN102981179B (zh) 用于闪烁探测器的位置表生成方法
CN109741325B (zh) 一种布线垂直度智能检测方法
Nandy et al. Automatic segmentation and supervised learning‐based selection of nuclei in cancer tissue images
Chen et al. An intelligent vision recognition method based on deep learning for pointer meters
CN109492706A (zh) 一种基于循环神经网络的染色体分类预测装置
US11804029B2 (en) Hierarchical constraint (HC)-based method and system for classifying fine-grained graptolite images
CN113537026A (zh) 建筑平面图中的图元检测方法、装置、设备及介质
Chen et al. Locating crop plant centers from UAV-based RGB imagery
CN113096080A (zh) 图像分析方法及系统
Di Cataldo et al. ANAlyte: A modular image analysis tool for ANA testing with indirect immunofluorescence
CN110309858A (zh) 基于判别学习的细粒度图像分类算法
Hu et al. Automatic detection of tuberculosis bacilli in sputum smear scans based on subgraph classification
Wang et al. GrainNet: efficient detection and counting of wheat grains based on an improved YOLOv7 modeling
Faezi et al. Multi-Spectral Source-Segmentation Using Semantically-Informed Max-Trees
US11067705B1 (en) CNN-based abnormality score map for SPECT gamma camera flood analysis
CN111723737B (zh) 一种基于多尺度匹配策略深度特征学习的目标检测方法
Gao et al. Spatio-temporal processing for automatic vehicle detection in wide-area aerial video
CN113780168B (zh) 一种高光谱遥感影像端元束自动提取方法

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 20855620

Country of ref document: EP

Kind code of ref document: A1

ENP Entry into the national phase

Ref document number: 2021563135

Country of ref document: JP

Kind code of ref document: A

NENP Non-entry into the national phase

Ref country code: DE

ENP Entry into the national phase

Ref document number: 2020855620

Country of ref document: EP

Effective date: 20220321

WWG Wipo information: grant in national office

Ref document number: 2020855620

Country of ref document: EP