Detection and identification method for visible moving foreign matters and bubbles in liquid medicine
Technical Field
The invention belongs to the technical field of image processing and automation, and relates to a method for detecting and identifying moving visible foreign matters and bubbles in liquid medicines.
Background
At present, liquid medicines widely used in our country include injection, infusion solution, oral liquid, eye drops and the like, and are usually filled in a glass bottle or plastic bottle mode, and due to the reasons of poor filtration, poor container cleaning, collision during packaging and the like in the production process of the medicines, visible foreign matters such as glass chips, suspended matters, hairs and the like exist in the liquid medicine, and in order to detect whether the visible foreign matters are mixed in the filled liquid medicine, the visible foreign matter detection needs to be carried out on each bottle of the medicines. The Chinese pharmacopoeia carries out clarity detection and insoluble particle detection on the detection of the medicines, namely the regulations from 1985, and strictly controls foreign matters and particles in the medicines. 24.8.1999, the national drug administration formally issued a notice on the implementation of < good manufacturing practice > regulations, clearly requiring that drug manufacturing enterprises must pass GMP (good) certification by the end of 2000. The "chinese pharmacopoeia" 2005 edition was executed from 7/1 in 2005, and the original "clarity inspection rules and judgment standards" were revised to "inspection method for visible foreign matter", which further makes clear that drugs such as injections and eye drops must be produced under the conditions conforming to the Good Manufacturing Practice (GMP), and that products should be inspected one by an appropriate method before shipment and rejected at the same time, and the visible foreign matter refers to an insoluble substance existing in the injections and eye drops and visually observed under the specified conditions, and the particle size or length of the insoluble substance is usually greater than 50 um. The visual foreign matter inspection method includes a light inspection method and a light scattering method, and the light inspection method is generally used. At present, the detection of visible foreign matters in domestic medicines mostly adopts a manual lamp detection method, the medicine bottle container is manually and lightly rotated or turned, the visible foreign matters possibly existing in the medicine liquid move (the medicine liquid does not generate bubbles) and are judged by visual inspection through human eyes, and if the visible foreign matters exist, the visible foreign matters are manually removed. The problems of low efficiency, high omission factor, low precision and the like which are inevitable during manual detection, poor repeatability and consistency of detection results, and harm to the health of lamp detection personnel caused by long-time detection work. Along with the popularization of GMP authentication, the automation and intelligence level of Chinese medicine production equipment and detection equipment are greatly improved, a semi-automatic lamp inspection machine and a full-automatic lamp inspection machine for detecting visible foreign matters in medicines are provided, wherein the semi-automatic lamp inspection machine adopts simple mechanical transmission and an optical system to realize the semi-automation of manual lamp inspection, the detection device can reduce manual labor to a certain extent, and the detection precision and the accuracy are not improved through lamp inspection industrial analysis and judgment. The full-automatic lamp inspection machine is only capable of being produced by a few companies in a few countries such as Germany, Japan, Italy and the like, but the application of the full-automatic lamp inspection machine in foreign countries is still few in China, and the following problems exist due to the full-automatic lamp inspection machine: the filtration and packaging materials in the domestic pharmaceutical process are very different from those in foreign countries, so that the detection quality is low. However, with the continuous improvement of the domestic pharmaceutical environment, the advanced medicine detection equipment in foreign countries will tend to change the detection status of the domestic pharmaceutical market. In the face of the impact of the foreign advanced medicine equipment sets, the research and development of the medicine detection equipment with independent intellectual property rights and the core detection and identification method have very important values.
In the process of visual detection of foreign matters in medical liquid, a specially-made mechanical device is needed to drive the foreign matters precipitated at the bottom of the liquid medicine to the upper part of the liquid so as to facilitate imaging of a camera, but because the movement of the mechanical device often causes bubbles and noise to exist in an imaged image, the technical problems in identification of foreign matters targets and bubbles in the image are mainly as follows:
1. foreign matters are different in size, the detection resolution ratio is up to 50um, the foreign matter target does not have information such as size, shape and texture, and the traditional image processing method cannot be applied;
2. due to the fact that scale embossments and spots exist on the surface of the bottle, the background of an imaged image is complex, and noise is various;
3. the signal-to-noise ratio in the image is low, and the foreign object target is often submerged in the background noise;
4. the detection speed of the production line is high in requirement, and the real-time performance of the identification algorithm is greatly required.
Disclosure of Invention
The invention aims to solve the technical problem of providing a detection and identification method for detecting moving visible foreign matters and bubbles in liquid medicines, which is used for rapidly detecting unqualified products with foreign matters in an automatic liquid medicine detection line, is matched with an imaging hardware system to carry out signal analysis and processing, and provides a reliable detection method for medicine detection equipment, thereby improving the detection precision and repeatability of the visible foreign matters in the medicines, thoroughly solving the problem of high false detection rate in the detection and meeting the performance requirements of the domestic existing lamp detection system.
The technical scheme adopted by the invention for solving the technical problems is as follows:
a method for detecting and identifying visible foreign bodies and bubbles moving in liquid medicines is characterized by comprising the following steps:
1) acquiring a plurality of frames of continuous images of a current liquid medicine to be detected;
2) image preprocessing: performing top-hat (top-hat) morphological filtering on the acquired image by using a 7 x 7 circular template as a structural element to remove various kinds of bright noise and dark noise smaller than the structural element to obtain a de-noised image;
3) extracting a moving target: according to the motion characteristic that the foreign object is represented as interframe on the time domain, a novel three-frame image difference algorithm is utilized to extract the motion object in the image, and the operation process is as follows: performing differential processing on the first frame and the second frame, performing differential processing on the second frame and the third frame, and then enhancing a target image through image accumulation to obtain an enhanced image;
4) and (3) moving object segmentation: applying an improved two-dimensional maximum entropy threshold segmentation algorithm to the enhanced image, and selecting a threshold by using a two-dimensional histogram consisting of gray levels and neighborhood average gray levels to obtain a segmented image;
5) tracking a moving target: when 3 continuous frames of possible targets in the segmented image are detected to appear in a certain neighborhood, the motion track is determined to be started; updating the target state by adopting Kalman filtering based on a uniform linear motion model according to the regularity of the motion of a foreign object target and a bubble target in the image, and generating an elliptic tracking gate; when only one measured value falls into the elliptic tracking gate, directly updating the target track; when more than one measured value falls into the elliptic tracking gate, the target track is updated by the measured value closest to the further predicted value; when the target is lost, extrapolating the possible positions of the targets of several frames after prediction by adopting the filtering value of the previous moment, and keeping stable memory tracking capability;
6) image identification and judgment: and obtaining a target motion track through a tracking result, and distinguishing a foreign matter target and a bubble target by utilizing the continuity and the directionality of the target track, thereby finally judging whether the current liquid medicine is qualified.
The step 2) is as follows:
assuming that the gray function of the image is f (x, y), the structural elements are B (x, y), and (x, y) are pixel coordinates, 4 basic operations are defined as follows:
and (3) corrosion: (f Θ B) (x, y) ═ min { f (x + i, y + j) -B (i, j) | (x + i, y + j) ∈ Df;(i,j)∈DB};
Expansion: <math><mrow> <mrow> <mo>(</mo> <mi>f</mi> <mo>⊕</mo> <mi>B</mi> <mo>)</mo> </mrow> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>=</mo> <mi>max</mi> <mrow> <mo>{</mo> <mi>f</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>-</mo> <mi>i</mi> <mo>,</mo> <mi>y</mi> <mo>-</mo> <mi>j</mi> <mo>)</mo> </mrow> <mo>+</mo> <mi>B</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>)</mo> </mrow> <mo>|</mo> <mrow> <mo>(</mo> <mi>x</mi> <mo>-</mo> <mi>i</mi> <mo>,</mo> <mi>y</mi> <mo>-</mo> <mi>j</mi> <mo>)</mo> </mrow> <mo>∈</mo> <msub> <mi>D</mi> <mi>f</mi> </msub> <mo>;</mo> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>)</mo> </mrow> <mo>∈</mo> <msub> <mi>D</mi> <mi>B</mi> </msub> <mo>}</mo> </mrow> <mo>;</mo> </mrow></math>
and (3) closed operation: <math><mrow> <mi>f</mi> <mo>·</mo> <mi>B</mi> <mo>=</mo> <mrow> <mo>(</mo> <mi>f</mi> <mo>⊕</mo> <mi>B</mi> <mo>)</mo> </mrow> <mi>ΘB</mi> <mo>;</mo> </mrow></math>
in the formula, DfAnd DBAre the domain definitions of the functions f and B,
introducing a top-hat transform filter operator: hat (f) ═ f- (foB);
in the formula (foB), the structural element B performs a gray-scale opening operation on the image f, and a 7 × 7 circular template is selected as the structural element to perform top-hat morphological filtering on the original image, so as to obtain an image containing foreign objects and bubbles, i.e., a denoised image.
The step 3) is as follows:
aiming at the denoised image, carrying out differential processing on a first frame and a second frame in each continuous three-frame image according to an image sequence, carrying out differential processing on the second frame and a third frame, and then enhancing a target image according to the following formula to obtain an enhanced image:
d(x,y)=|f1(x,y)-f2(x,y)|×|f2(x,y)-f3(x,y)|;
wherein f is1(x,y),f2(x,y),f3(x, y) are respectively the first frame, the second frame and the third frame of image after the high-cap filtering, | f1(x,y)-f2(x, y) | is a differential image of the first frame and the second frame, | f2(x,y)-f3And (x, y) | is a differential image of the second frame and the third frame. .
The step 4) comprises the following steps:
(a) segmenting the image into a series of sub-images;
(b) calculating the threshold value of each sub-image by a two-dimensional maximum entropy algorithm;
(c) segmenting the image according to the threshold value of each sub-image;
the method for solving the threshold value comprises the following steps: setting the size of an image as M multiplied by N, wherein the image has L gray levels, the neighborhood gray level average value of each pixel point is also L gray levels, and the total gray level is L multiplied by L; making a two-dimensional gray level histogram on a two-dimensional plane, forming a gray level binary value by the gray level value and the neighborhood gray level average value, recording the gray level binary value as the coordinate value of each pixel point as (i, j), wherein i and j respectively represent the gray level of the pixel point and the neighborhood gray level average value, and the frequency of the occurrence of the gray level combination is fijThe pixel point function value is a joint probability density function pij:
<math><mrow><msub><mi>p</mi><mi>ij</mi></msub><mo>=</mo><mfrac><msub><mi>f</mi><mi>ij</mi></msub><mrow><mi>M</mi><mo>×</mo><mi>N</mi></mrow></mfrac><mo>,</mo><mi>i</mi><mo>,</mo><mi>j</mi><mo>=</mo><mn>0,1</mn><mo>,</mo><mo>.</mo><mo>.</mo><mo>.</mo><mi>L</mi><mo>-</mo><mn>1</mn><mo>;</mo></mrow></math>
Probability distribution p for object O with (s, t) as segmentation thresholdij( i 1, 2,. s; j 1, 2,. t) performing a normalization operation; the entropy associated with the probability distribution of the object O after the normalization operation is defined as:
<math><mrow><mi>H</mi><mrow><mo>(</mo><mi>O</mi><mo>)</mo></mrow><mo>=</mo><mo>-</mo><munderover><mi>Σ</mi><mrow><mi>i</mi><mo>=</mo><mn>0</mn></mrow><mi>s</mi></munderover><munderover><mi>Σ</mi><mrow><mi>j</mi><mo>=</mo><mn>0</mn></mrow><mi>t</mi></munderover><mfrac><msub><mi>p</mi><mi>ij</mi></msub><msub><mi>P</mi><mi>st</mi></msub></mfrac><mi>ln</mi><mfrac><msub><mi>p</mi><mi>ij</mi></msub><msub><mi>P</mi><mi>st</mi></msub></mfrac><mo>=</mo><mo>-</mo><mfrac><mn>1</mn><msub><mi>P</mi><mi>st</mi></msub></mfrac><munderover><mi>Σ</mi><mrow><mi>i</mi><mo>=</mo><mn>0</mn></mrow><mi>s</mi></munderover><munderover><mi>Σ</mi><mrow><mi>j</mi><mo>=</mo><mn>0</mn></mrow><mi>t</mi></munderover><mrow><mo>(</mo><msub><mi>p</mi><mi>ij</mi></msub><mi>ln</mi><msub><mi>p</mi><mi>ij</mi></msub><mo>-</mo><msub><mi>p</mi><mi>ij</mi></msub><mi>ln</mi><msub><mi>P</mi><mi>st</mi></msub><mo>)</mo></mrow><mo>=</mo><mi>ln</mi><msub><mi>P</mi><mi>st</mi></msub><mo>+</mo><mfrac><msub><mi>H</mi><mi>st</mi></msub><msub><mi>P</mi><mi>st</mi></msub></mfrac><mo>;</mo></mrow></math>
wherein, <math><mrow><msub><mi>P</mi><mi>st</mi></msub><mo>=</mo><munderover><mi>Σ</mi><mrow><mi>i</mi><mo>=</mo><mn>0</mn></mrow><mi>s</mi></munderover><munderover><mi>Σ</mi><mrow><mi>j</mi><mo>=</mo><mn>0</mn></mrow><mi>t</mi></munderover><msub><mi>p</mi><mi>ij</mi></msub><mo>,</mo></mrow></math> <math><mrow><msub><mi>H</mi><mi>st</mi></msub><mo>=</mo><mo>-</mo><munderover><mi>Σ</mi><mrow><mi>i</mi><mo>=</mo><mn>0</mn></mrow><mi>s</mi></munderover><munderover><mi>Σ</mi><mrow><mi>j</mi><mo>=</mo><mn>0</mn></mrow><mi>t</mi></munderover><msub><mi>p</mi><mi>ij</mi></msub><mi>ln</mi><msub><mi>p</mi><mi>ij</mi></msub><mo>;</mo></mrow></math>
entropy associated with the probability distribution of background B is
<math><mrow><mi>H</mi><mrow><mo>(</mo><mi>B</mi><mo>)</mo></mrow><mo>=</mo><mo>-</mo><munderover><mi>Σ</mi><mrow><mi>i</mi><mo>=</mo><mi>s</mi><mo>+</mo><mn>1</mn></mrow><mrow><mi>L</mi><mo>-</mo><mn>1</mn></mrow></munderover><munderover><mi>Σ</mi><mrow><mi>j</mi><mo>=</mo><mi>t</mi><mo>+</mo><mn>1</mn></mrow><mrow><mi>L</mi><mo>-</mo><mn>1</mn></mrow></munderover><mfrac><msub><mi>p</mi><mi>ij</mi></msub><mrow><mn>1</mn><mo>-</mo><msub><mi>P</mi><mi>st</mi></msub></mrow></mfrac><mi>ln</mi><mfrac><msub><mi>p</mi><mi>ij</mi></msub><mrow><mn>1</mn><mo>-</mo><msub><mi>P</mi><mi>st</mi></msub></mrow></mfrac><mo>=</mo><mo>-</mo><mfrac><mn>1</mn><mrow><mn>1</mn><mo>-</mo><msub><mi>P</mi><mi>st</mi></msub></mrow></mfrac><munderover><mi>Σ</mi><mrow><mi>i</mi><mo>=</mo><mi>s</mi><mo>+</mo><mn>1</mn></mrow><mrow><mi>L</mi><mo>-</mo><mn>1</mn></mrow></munderover><munderover><mi>Σ</mi><mrow><mi>j</mi><mo>=</mo><mi>t</mi><mo>+</mo><mn>1</mn></mrow><mrow><mi>L</mi><mo>-</mo><mn>1</mn></mrow></munderover><mo>[</mo><msub><mi>p</mi><mi>ij</mi></msub><mi>ln</mi><msub><mi>p</mi><mi>ij</mi></msub><mo>-</mo><msub><mi>p</mi><mi>ij</mi></msub><mi>ln</mi><mrow><mo>(</mo><mn>1</mn><mo>-</mo><msub><mi>P</mi><mi>st</mi></msub><mo>)</mo></mrow><mo>]</mo><mo>;</mo></mrow></math>
Wherein, <math><mrow><msub><mi>H</mi><mrow><mi>L</mi><mo>-</mo><mn>1</mn><mi>L</mi><mo>-</mo><mn>1</mn></mrow></msub><mo>=</mo><mo>-</mo><munderover><mi>Σ</mi><mrow><mi>i</mi><mo>=</mo><mn>0</mn></mrow><mrow><mi>L</mi><mo>-</mo><mn>1</mn></mrow></munderover><munderover><mi>Σ</mi><mrow><mi>j</mi><mo>=</mo><mn>0</mn></mrow><mrow><mi>L</mi><mo>-</mo><mn>1</mn></mrow></munderover><msub><mi>p</mi><mi>ij</mi></msub><mi>ln</mi><msub><mi>p</mi><mi>ij</mi></msub><mo>;</mo></mrow></math>
the criterion function Ψ (s, t) is taken as the sum of H (O), H (B), i.e.
<math><mrow><mi>Ψ</mi><mrow><mo>(</mo><mi>s</mi><mo>,</mo><mi>t</mi><mo>)</mo></mrow><mo>=</mo><mi>H</mi><mrow><mo>(</mo><mi>O</mi><mo>)</mo></mrow><mo>+</mo><mi>H</mi><mrow><mo>(</mo><mi>B</mi><mo>)</mo></mrow><mo>=</mo><mi>ln</mi><msub><mi>P</mi><mi>st</mi></msub><mo>+</mo><mfrac><msub><mi>H</mi><mi>st</mi></msub><msub><mi>P</mi><mi>st</mi></msub></mfrac><mo>+</mo><mi>ln</mi><mrow><mo>(</mo><mn>1</mn><mo>-</mo><msub><mi>P</mi><mi>st</mi></msub><mo>+</mo><mfrac><mrow><msub><mi>H</mi><mrow><mi>L</mi><mo>-</mo><mn>1</mn><mi>L</mi><mo>-</mo><mn>1</mn></mrow></msub><mo>-</mo><msub><mi>H</mi><mi>st</mi></msub></mrow><mrow><mn>1</mn><mo>-</mo><msub><mi>P</mi><mi>st</mi></msub></mrow></mfrac><mo>)</mo></mrow><mo>;</mo></mrow></math>
The gray-scale value (s, t) that maximizes Ψ (s, t) is the optimum threshold value t that is obtained*I.e. by
<math><mrow> <msup> <mi>t</mi> <mo>*</mo> </msup> <mo>=</mo> <mi>Arg</mi> <mrow> <mo>{</mo> <munder> <mi>Max</mi> <mrow> <mi>s</mi> <mo>,</mo> <mi>t</mi> <mo>∈</mo> <mi>G</mi> </mrow> </munder> <mi>Ψ</mi> <mrow> <mo>(</mo> <mi>s</mi> <mo>,</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>}</mo> </mrow> <mo>;</mo> </mrow></math>
Wherein Arg { } denotes an inverse function.
As a refinement, the step 4) further comprises: a genetic algorithm is introduced to optimize a segmentation threshold value of a two-dimensional maximum entropy method so as to shorten the image processing time; the method comprises the following specific steps:
initializing a population: randomly generating and selecting uniform distribution;
population scale: through analysis, the population scale is popsize 20, and the maximum evolution generation number is gmax=100;
Fitness function: directly using the objective function formula Ψ (s, t) as a fitness function;
and (3) encoding: since the gray value of the image is between 0 and 255, and an 8-bit gray image is adopted, the threshold parameter satisfies s which is more than or equal to 0 and t which is less than or equal to 255, the individual is coded into a 16-bit binary code, the high 8 bits represent s, and the low 8 bits represent t;
and (3) decoding: decoding the upper 8 bits of the 16-bit binary code into an s value and the lower 8 bits into a t value;
selecting: adopting a proportion selection operator;
and (3) a crossover operator: adopting double-point crossing, wherein two randomly generated crossing points are respectively positioned at the front 8 bits and the back 8 bits of the coding string; performing cross operation according to the cross probability Pc; selecting Pc as 0.8 in the early stage of searching, and selecting Pc as 0.6 in the later stage of searching;
mutation operator: the mutation operation is bit-wise negation; defining the mutation probability Pm as a parabola-like mutation operator:
<math><mrow> <msub> <mi>P</mi> <mi>m</mi> </msub> <mo>=</mo> <msub> <mi>P</mi> <mi>m</mi> </msub> <mrow> <mo>(</mo> <mi>g</mi> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mi>P</mi> <mrow> <mi>m</mi> <mo>_</mo> <mi>max</mi> </mrow> </msub> <mo>-</mo> <mfrac> <mrow> <mn>4</mn> <mo>×</mo> <mrow> <mo>(</mo> <msub> <mi>P</mi> <mrow> <mi>m</mi> <mo>_</mo> <mi>max</mi> </mrow> </msub> <mo>-</mo> <msub> <mi>P</mi> <mrow> <mi>m</mi> <mo>_</mo> <mi>min</mi> </mrow> </msub> <mo>)</mo> </mrow> </mrow> <msubsup> <mi>g</mi> <mi>max</mi> <mn>2</mn> </msubsup> </mfrac> <msup> <mrow> <mo>(</mo> <mi>g</mi> <mo>-</mo> <mfrac> <msub> <mi>g</mi> <mi>max</mi> </msub> <mn>2</mn> </mfrac> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mrow></math>
wherein g is evolution algebra, and g belongs to [1, g ∈max],Pm_min=Pb,Pm_max=10Pm_min
PbFor the basic variation probability, the estimation is performed as follows:
<math><mrow> <msub> <mi>P</mi> <mi>b</mi> </msub> <mo>≈</mo> <mfrac> <mn>175</mn> <mrow> <mi>popsize</mi> <mo>×</mo> <msqrt> <mi>bits</mi> </msqrt> </mrow> </mfrac> <mo>;</mo> </mrow></math>
wherein bits represents the number of bits of the individual; when popsize is 20 and bits is 16, Pb≈0.02;
When P is presentb≈0.02,gmaxWhen 100, we can get:
Pm=Pm(g)=0.2-0.00008×(g-50)2;
termination criteria: when the algorithm is executed to the maximum algebra or the highest fitness value in the group is stable, the algorithm stops running, and the individual with the highest fitness value is the solution; at the termination of the algorithm, the individual with the highest fitness value is taken as the best threshold.
The step 5) comprises the following steps:
a single-frame detection result can be determined to be obtained through an improved two-dimensional maximum entropy threshold segmentation algorithm, and a target is approximated by a uniform linear motion model; the tracking system model is as follows: let the state vector be <math><mrow> <mi>X</mi> <mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo> </mrow> <mo>=</mo> <msup> <mrow> <mo>(</mo> <mi>x</mi> <mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo> </mrow> <mo>,</mo> <mover> <mi>x</mi> <mo>·</mo> </mover> <mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo> </mrow> <mo>,</mo> <mi>y</mi> <mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo> </mrow> <mo>,</mo> <mover> <mi>y</mi> <mo>·</mo> </mover> <mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo> </mrow> <mo>)</mo> </mrow> <mi>T</mi> </msup> <mo>,</mo> </mrow></math> The measurement vector is Z (k) ═ x (k), y (k)]TThe system equation can be written as
X(k+1)=ΦX(k)+Gw(k);
Z(k+1)=HX(k)+v(k);
Wherein <math><mrow> <mi>Φ</mi> <mo>=</mo> <mrow> <mfenced open='[' close=']'> <mtable> <mtr> <mtd> <mn>1</mn> </mtd> <mtd> <mi>T</mi> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mn>1</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mn>1</mn> </mtd> <mtd> <mi>T</mi> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mn>1</mn> </mtd> </mtr> </mtable> </mfenced> <mo>,</mo> </mrow> </mrow></math>
w (k), v (k) is zero-mean white Gaussian noise with variance Q (k) and R (k), respectively; wherein T is a sampling period; track start: if 3 consecutive frames of the detected possible target appear in the 3 × 3 neighborhood, the track starts;
tracking gates and data association: assume that the state vector at time k-1 is estimated as
The one-step prediction of the state vector is:
<math><mrow> <mover> <mi>X</mi> <mo>^</mo> </mover> <mrow> <mo>(</mo> <mi>k</mi> <mo>|</mo> <mi>k</mi> <mo>-</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>=</mo> <mi>Φ</mi> <mover> <mi>X</mi> <mo>^</mo> </mover> <mrow> <mo>(</mo> <mi>k</mi> <mo>-</mo> <mn>1</mn> <mo>|</mo> <mi>k</mi> <mo>-</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>;</mo> </mrow></math>
the innovation of the system is:
the innovation variance matrix is:
S(k)=HP(k|k-1)HT+R(k);
wherein P (k | k-1) is the one-step prediction variance;
let the norm of the innovation vector d (k) be:
g(k)=dT(k)S-1(k)d(k);
wherein g (k) is subject to
Distribution, M1 is the measurement dimension;
defining a tracking gate in the measurement space, so that the measurements are distributed within the tracking gate with a certain probability;
<math><mrow> <mover> <mi>V</mi> <mo>~</mo> </mover> <mrow> <mo>(</mo> <mi>γ</mi> <mo>)</mo> </mrow> <mo>=</mo> <mrow> <mo>[</mo> <mi>Z</mi> <mo>:</mo> <mi>g</mi> <mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo> </mrow> <mo>≤</mo> <mi>γ</mi> <mo>]</mo> </mrow> <mo>,</mo> </mrow></math> wherein γ can pass through x2Checking the difference of the distribution table;
if only one echo exists in the tracking gate after a certain frame is processed, the target track is directly updated; if more than one echo exists in the tracking gate, the target track is updated through a measured value which is closest to the one-step predicted value;
the track of the target is updated by a standard Kalman filtering algorithm
X(k|k-1)=ΦX(k-1|k-1)
P(k|k-1)=ΦP(k-1|k-1)ΦT+GQ(k-1)GT
K(k)=P(k|k-1)HT[HP(k|k-1)HT+R]-1;
X(k|k)=X(k|k-1)+K(k)[Z(k)-H(k)X(k|k-1)]
P(k|k)=[I-K(k)H]P(k|k-1)
After the target is detected and tracked, if the target disappears momentarily, the next possible position of the target can be predicted according to the previous position information and motion state of the target, and when the target appears again, the target can still be stably tracked without losing the target, and the specific process is as follows:
assume that the state vector at time k is estimated asThen the predicted value of the target at k + n frame is:
<math><mrow> <mover> <mi>X</mi> <mo>^</mo> </mover> <mrow> <mo>(</mo> <mi>k</mi> <mo>+</mo> <mi>n</mi> <mo>|</mo> <mi>k</mi> <mo>)</mo> </mrow> <mo>=</mo> <msup> <mi>Φ</mi> <mi>n</mi> </msup> <mover> <mi>X</mi> <mo>^</mo> </mover> <mrow> <mo>(</mo> <mi>k</mi> <mo>|</mo> <mi>k</mi> <mo>)</mo> </mrow> <mo>;</mo> </mrow></math> selecting n<6。
The step 6) is as follows: for each assumed track obtained by target tracking, whether the track is formed by foreign matters or noise or bubbles is distinguished according to the following 3 principles, wherein the judgment criteria are as follows:
(a) because the image is shot when the medicine bottle is still after being turned over by 180 degrees, the visible foreign matters move downwards, the centroids and the ordinate of the image should be sequentially enlarged, and the upper left corner is taken as the origin of coordinates;
(b) the bubbles move upwards in the opposite direction to the foreign bodies;
(c) the track curve formed by the foreign matters is smooth, and the track formed by the noise is disordered;
therefore, if the centroids of the detected motion tracks become smaller in sequence, the motion tracks are indicated to be generated by bubbles; and if the track with smooth motion track and sequentially enlarged centroid longitudinal coordinates is detected, indicating that foreign matters exist, and judging that the current liquid medicine is unqualified.
The invention has the following beneficial effects:
the invention adopts mathematical form filtering to obtain the preprocessed image, and the algorithm can be realized in parallel by hardware, thereby greatly improving the processing speed. In the extraction of the moving target, the defects that the bottle side wall possibly has difference and the detection effect of simple sequence image difference on the small foreign object target is poor are well overcome by using a novel three-frame image difference algorithm, the signal-to-noise ratio of an output image is obviously improved, and the difficulty of a subsequent segmentation detection algorithm is greatly simplified. The improved two-dimensional maximum entropy threshold motion target segmentation algorithm fully utilizes the gray distribution information of target pixels in the image and the related information among the pixels, improves the anti-noise capability of threshold segmentation, optimizes the parameters of the two-dimensional maximum entropy method through a genetic algorithm, reduces the image processing time, and greatly improves the real-time property and the robustness of the algorithm. According to the motion characteristics of the moving target in the image, the target state is updated by adopting Kalman filtering based on a uniform linear motion model, and the calculated amount of filtering is simplified under the condition of ensuring the tracking precision. The invention fully simplifies the calculation under the condition of ensuring the accuracy of the algorithm, improves the real-time performance of the algorithm, greatly improves the accuracy of detecting and identifying foreign matters and bubbles in the liquid medicine, reduces the missing rate and the false detection rate of the medicine, can be applied to various detection systems of medicine, beverage, wine and the like, and has wide market prospect and practical application value.
Drawings
FIG. 1 is a general flow chart of a method in accordance with the present invention;
FIG. 2 is an original image of a liquid drug within a vial taken in accordance with the present invention;
FIG. 3 is an enhanced image obtained after target extraction in the present invention;
FIG. 4 is an image obtained after segmentation of a moving object according to the present invention;
fig. 5 is a motion trail diagram obtained by tracking a moving object in the invention.
Description of reference numerals: 1-bubble trajectory, 2-foreign body trajectory.
Detailed Description
The following will be described in further detail with reference to fig. 1 to 5 and specific examples.
Example 1:
as shown in fig. 1, the present invention provides a method for identifying motion-visible foreign objects and bubbles for liquid medicine detection. The specific implementation details of each part are as follows:
1. top-hat morphological filtering process
And filtering the image noise of the acquired image by adopting gray morphological filtering. The most basic operations of gray scale morphology are erosion, dilation, open and close operations. If the gray scale function of the image is f (x, y) and the structural element is B (x, y), 4 basic operations are defined as follows; (ii) a Where (x, y) is the pixel coordinate;
and (3) corrosion: (f Θ B) (x, y) ═ min { f (x + i, y + j) -B (i, j) | (x + i, y + j) ∈ Df;(i,j)∈DB};
Expansion: <math><mrow> <mrow> <mo>(</mo> <mi>f</mi> <mo>⊕</mo> <mi>B</mi> <mo>)</mo> </mrow> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>=</mo> <mi>max</mi> <mrow> <mo>{</mo> <mi>f</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>-</mo> <mi>i</mi> <mo>,</mo> <mi>y</mi> <mo>-</mo> <mi>j</mi> <mo>)</mo> </mrow> <mo>+</mo> <mi>B</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>)</mo> </mrow> <mo>|</mo> <mrow> <mo>(</mo> <mi>x</mi> <mo>-</mo> <mi>i</mi> <mo>,</mo> <mi>y</mi> <mo>-</mo> <mi>j</mi> <mo>)</mo> </mrow> <mo>∈</mo> <msub> <mi>D</mi> <mi>f</mi> </msub> <mo>;</mo> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>)</mo> </mrow> <mo>∈</mo> <msub> <mi>D</mi> <mi>B</mi> </msub> <mo>}</mo> </mrow> <mo>;</mo> </mrow></math>
and (3) closed operation: <math><mrow> <mi>f</mi> <mo>·</mo> <mi>B</mi> <mo>=</mo> <mrow> <mo>(</mo> <mi>f</mi> <mo>⊕</mo> <mi>B</mi> <mo>)</mo> </mrow> <mi>ΘB</mi> <mo>;</mo> </mrow></math>
in the formula, DfAnd DBThe domains of the functions f and B, respectively. As a result of the erosion of the gray image, a portion darker than the background is expanded, and a portion brighter than the background is contracted; as a result of expanding the gray image, the portion brighter than the background is expanded, and the portion darker than the background is contracted; open operation graphThe image contour becomes smooth, breaking up narrow discontinuities and eliminating thin protrusions; the closed operation also smoothes the contour of the image, but in contrast to the open operation, it eliminates narrow discontinuities and long narrow gaps, eliminates small holes, and fills up cracks in the contour lines.
The Top-hat transform filter operator based on gray level morphology is defined as:
[0118]Hat(f)=f-(fοB); (1)
[0119]wherein B is a structural element, (f ° B) is the structural element B, and the gray level on operation is performed on the image f, and the structure thereof is the estimated background, and hat (f) is obtained after subtraction, that is, the target image from which noise is filtered, and of course, the filtered image also contains some single-point impulse noise, which can be filtered in the following processing. Because the noise in the acquired image is a small bright spot, a 7 × 7 circular template is selected as a structural element to perform Top-hat morphological filtering on the original image, so as to obtain a denoised image containing a foreign object and bubbles, as shown in fig. 3.
2. Moving object extraction based on three-frame image difference algorithm
And aiming at the images subjected to Top-hat filtering denoising, performing differential processing on the first frame and the second frame according to the image sequence, performing differential processing on the second frame and the third frame, and then enhancing the target image according to the following formula:
d(x,y)=|f1(x,y)-f2(x,y)|×|f2(x,y)-f3(x,y)| (2)
wherein f is1(x,y),f2(x,y),f3(x, y) are respectively the images of the first frame, the second frame and the third frame after Top-hat filtering, | f1(x,y)-f2(x, y) | is a differential image of the first frame and the second frame, | f2(x,y)-f3(x, y) | is a differential image of the second frame and the third frame;
after the three-frame image difference algorithm processing, only the pixel point position corresponding to the moving target in d (x, y) is not zero. However, in practice, due to the existence of various noises, the positions of many pixel points outside the moving object in the difference image are not zero, which requires the subsequent processing to effectively identify the foreign object information.
3. Moving object segmentation
The size of the image is M multiplied by N, the image has L gray levels, the average value of the neighborhood gray level of each pixel point is also the L gray levels, and the total gray level is L multiplied by L. Making a two-dimensional gray histogram on a plane, forming a gray binary value by a gray value and a neighborhood gray average value, and marking as (i, j), wherein i and j respectively represent the gray value of a pixel point and the neighborhood gray average value, and the frequency of the occurrence of the gray combination is fijThe pixel point function value is a joint probability density function pij。
<math><mrow><msub><mi>p</mi><mi>ij</mi></msub><mo>=</mo><mfrac><msub><mi>f</mi><mi>ij</mi></msub><mrow><mi>M</mi><mo>×</mo><mi>N</mi></mrow></mfrac><mo>,</mo><mi>i</mi><mo>,</mo><mi>j</mi><mo>=</mo><mn>0,1</mn><mo>,</mo><mo>.</mo><mo>.</mo><mo>.</mo><mi>L</mi><mo>-</mo><mn>1</mn><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>3</mn><mo>)</mo></mrow></mrow></math>
If (s, t) is used as the segmentation threshold, the probability distribution p corresponding to the object OijThe sum of ( i 1, 2.. s; j 1, 2.. t) should be 1, and thus a normalization operation is required.
The entropy associated with the probability distribution of the object O after the normalization operation is defined as:
<math><mrow><mi>H</mi><mrow><mo>(</mo><mi>O</mi><mo>)</mo></mrow><mo>=</mo><mo>-</mo><munderover><mi>Σ</mi><mrow><mi>i</mi><mo>=</mo><mn>0</mn></mrow><mi>s</mi></munderover><munderover><mi>Σ</mi><mrow><mi>j</mi><mo>=</mo><mn>0</mn></mrow><mi>t</mi></munderover><mfrac><msub><mi>p</mi><mi>ij</mi></msub><msub><mi>P</mi><mi>st</mi></msub></mfrac><mi>ln</mi><mfrac><msub><mi>p</mi><mi>ij</mi></msub><msub><mi>P</mi><mi>st</mi></msub></mfrac><mo>=</mo><mo>-</mo><mfrac><mn>1</mn><msub><mi>P</mi><mi>st</mi></msub></mfrac><munderover><mi>Σ</mi><mrow><mi>i</mi><mo>=</mo><mn>0</mn></mrow><mi>s</mi></munderover><munderover><mi>Σ</mi><mrow><mi>j</mi><mo>=</mo><mn>0</mn></mrow><mi>t</mi></munderover><mrow><mo>(</mo><msub><mi>p</mi><mi>ij</mi></msub><mi>ln</mi><msub><mi>p</mi><mi>ij</mi></msub><mo>-</mo><msub><mi>p</mi><mi>ij</mi></msub><mi>ln</mi><msub><mi>P</mi><mi>st</mi></msub><mo>)</mo></mrow><mo>=</mo><mi>ln</mi><msub><mi>P</mi><mi>st</mi></msub><mo>+</mo><mfrac><msub><mi>H</mi><mi>st</mi></msub><msub><mi>P</mi><mi>st</mi></msub></mfrac><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>4</mn><mo>)</mo></mrow></mrow></math>
wherein, <math><mrow><msub><mi>P</mi><mi>st</mi></msub><mo>=</mo><munderover><mi>Σ</mi><mrow><mi>i</mi><mo>=</mo><mn>0</mn></mrow><mi>s</mi></munderover><munderover><mi>Σ</mi><mrow><mi>j</mi><mo>=</mo><mn>0</mn></mrow><mi>t</mi></munderover><msub><mi>p</mi><mi>ij</mi></msub><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>5</mn><mo>)</mo></mrow></mrow></math>
<math><mrow><msub><mi>H</mi><mi>st</mi></msub><mo>=</mo><mo>-</mo><munderover><mi>Σ</mi><mrow><mi>i</mi><mo>=</mo><mn>0</mn></mrow><mi>s</mi></munderover><munderover><mi>Σ</mi><mrow><mi>j</mi><mo>=</mo><mn>0</mn></mrow><mi>t</mi></munderover><msub><mi>p</mi><mi>ij</mi></msub><mi>ln</mi><msub><mi>p</mi><mi>ij</mi></msub><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>6</mn><mo>)</mo></mrow></mrow></math>
similarly, the entropy associated with the probability distribution of background B is
<math><mrow><mi>H</mi><mrow><mo>(</mo><mi>B</mi><mo>)</mo></mrow><mo>=</mo><mo>-</mo><munderover><mi>Σ</mi><mrow><mi>i</mi><mo>=</mo><mi>s</mi><mo>+</mo><mn>1</mn></mrow><mrow><mi>L</mi><mo>-</mo><mn>1</mn></mrow></munderover><munderover><mi>Σ</mi><mrow><mi>j</mi><mo>=</mo><mi>t</mi><mo>+</mo><mn>1</mn></mrow><mrow><mi>L</mi><mo>-</mo><mn>1</mn></mrow></munderover><mfrac><msub><mi>p</mi><mi>ij</mi></msub><mrow><mn>1</mn><mo>-</mo><msub><mi>P</mi><mi>st</mi></msub></mrow></mfrac><mi>ln</mi><mfrac><msub><mi>p</mi><mi>ij</mi></msub><mrow><mn>1</mn><mo>-</mo><msub><mi>P</mi><mi>st</mi></msub></mrow></mfrac><mo>=</mo><mo>-</mo><mfrac><mn>1</mn><mrow><mn>1</mn><mo>-</mo><msub><mi>P</mi><mi>st</mi></msub></mrow></mfrac><munderover><mi>Σ</mi><mrow><mi>i</mi><mo>=</mo><mi>s</mi><mo>+</mo><mn>1</mn></mrow><mrow><mi>L</mi><mo>-</mo><mn>1</mn></mrow></munderover><munderover><mi>Σ</mi><mrow><mi>j</mi><mo>=</mo><mi>t</mi><mo>+</mo><mn>1</mn></mrow><mrow><mi>L</mi><mo>-</mo><mn>1</mn></mrow></munderover><mo>[</mo><msub><mi>p</mi><mi>ij</mi></msub><mi>ln</mi><msub><mi>p</mi><mi>ij</mi></msub><mo>-</mo><msub><mi>p</mi><mi>ij</mi></msub><mi>ln</mi><mrow><mo>(</mo><mn>1</mn><mo>-</mo><msub><mi>P</mi><mi>st</mi></msub><mo>)</mo></mrow><mo>]</mo><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>7</mn><mo>)</mo></mrow><mo>;</mo></mrow></math>
Wherein, <math><mrow><msub><mi>H</mi><mrow><mi>L</mi><mo>-</mo><mn>1</mn><mi>L</mi><mo>-</mo><mn>1</mn></mrow></msub><mo>=</mo><mo>-</mo><munderover><mi>Σ</mi><mrow><mi>i</mi><mo>=</mo><mn>0</mn></mrow><mrow><mi>L</mi><mo>-</mo><mn>1</mn></mrow></munderover><munderover><mi>Σ</mi><mrow><mi>j</mi><mo>=</mo><mn>0</mn></mrow><mrow><mi>L</mi><mo>-</mo><mn>1</mn></mrow></munderover><msub><mi>p</mi><mi>ij</mi></msub><mi>ln</mi><msub><mi>p</mi><mi>ij</mi></msub><mo>;</mo></mrow></math>
provided that p is present in the region (i ═ s + 1.. gtl-1 and j ═ 1, 2.. t) and in the region (i ═ 1, 2.. s and j ═ t + 1.., L-1)ijH (b) defined above can be established with 0, which is assumed to contribute to saving computation time. Since in many cases the probability distribution of the off-diagonal is negligible, it is reasonable to make this assumption.
The criterion function Ψ (s, t) is taken as the sum of H (O), H (B), i.e.
<math><mrow><mi>Ψ</mi><mrow><mo>(</mo><mi>s</mi><mo>,</mo><mi>t</mi><mo>)</mo></mrow><mo>=</mo><mi>H</mi><mrow><mo>(</mo><mi>O</mi><mo>)</mo></mrow><mo>+</mo><mi>H</mi><mrow><mo>(</mo><mi>B</mi><mo>)</mo></mrow><mo>=</mo><mi>ln</mi><msub><mi>P</mi><mi>st</mi></msub><mo>+</mo><mfrac><msub><mi>H</mi><mi>st</mi></msub><msub><mi>P</mi><mi>st</mi></msub></mfrac><mo>+</mo><mi>ln</mi><mrow><mo>(</mo><mn>1</mn><mo>-</mo><msub><mi>P</mi><mi>st</mi></msub><mo>+</mo><mfrac><mrow><msub><mi>H</mi><mrow><mi>L</mi><mo>-</mo><mn>1</mn><mi>L</mi><mo>-</mo><mn>1</mn></mrow></msub><mo>-</mo><msub><mi>H</mi><mi>st</mi></msub></mrow><mrow><mn>1</mn><mo>-</mo><msub><mi>P</mi><mi>st</mi></msub></mrow></mfrac><mo>)</mo></mrow><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>8</mn><mo>)</mo></mrow><mo>;</mo></mrow></math>
The gray-scale value (s, t) that maximizes Ψ (s, t) is the optimum threshold value t that is obtained*I.e. by
<math><mrow><msup><mi>t</mi><mo>*</mo></msup><mo>=</mo><mi>Arg</mi><mo>{</mo><munder><mi>Max</mi><mrow><mi>s</mi><mo>,</mo><mi>t</mi><mo>∈</mo><mi>G</mi></mrow></munder><mi>Ψ</mi><mrow><mo>(</mo><mi>s</mi><mo>,</mo><mi>t</mi><mo>)</mo></mrow><mo>}</mo><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>9</mn><mo>)</mo></mrow></mrow></math>
Wherein Arg { } denotes an inverse function.
Here, the optimal solution for the two parameter problem (s, t) is optimized using a genetic algorithm. The genetic algorithm is used as an efficient parallel global search method for solving problems, and is mainly characterized in that information exchange between a group search strategy and individuals in a group can automatically acquire and accumulate knowledge about a search space in the search process, and the search process is controlled in a self-adaptive mode to obtain an optimal solution or an approximately optimal solution. Genetic algorithms exploit simple coding techniques and propagation mechanisms to represent complex and nonlinear problems that are difficult to solve with traditional search methods.
Initializing a population: randomly generated and uniformly distributed.
Population scale: through analysis, the population scale is popsize 20, and the maximum evolution generation number is gmax=100。
Fitness function: the objective function (6) is used directly as the fitness function.
And (3) encoding: as the gray value of the image is between 0 and 255, and an 8-bit gray image is adopted for carrying out a simulation experiment, the threshold parameter satisfies s which is more than or equal to 0 and t which is less than or equal to 255, the individual code is a 16-bit binary code, s is represented by high 8 bits, and t is represented by low 8 bits.
And (3) decoding: the upper 8 bits of the 16-bit binary code are decoded to the s value and the lower 8 bits are decoded to the t value.
Selecting: adopting a proportion selection operator;
and (3) a crossover operator: two-point crossing is adopted, and two randomly generated crossing points are respectively positioned at the front 8 bits and the rear 8 bits of the coding string. And performing the crossing operation according to the preset crossing probability Pc. Too high a cross probability means that the individual updates faster, a larger solution space can be achieved, and the probability of obtaining a non-optimal solution can be reduced, but the destructiveness to the existing better model also increases. The crossover probability is too low and the search may become sluggish due to the reduction in search range. In the early stage of the search, Pc is 0.8, and in the later stage of the search, Pc is 0.6.
Mutation operator: because a binary coding mode is adopted, the mutation is bit-wise negation. Mutation controls the proportion of new genes introduced into the population. If the mutation probability is too low, some useful genes cannot be introduced, and if the mutation probability is too high, i.e., if the random variation is too much, the offspring may lose the good characteristics of parents. In this embodiment, the mutation probability Pm is defined as a parabola-like mutation operator:
<math><mrow> <msub> <mi>P</mi> <mi>m</mi> </msub> <mo>=</mo> <msub> <mi>P</mi> <mi>m</mi> </msub> <mrow> <mo>(</mo> <mi>g</mi> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mi>P</mi> <mrow> <mi>m</mi> <mo>_</mo> <mi>max</mi> </mrow> </msub> <mo>-</mo> <mfrac> <mrow> <mn>4</mn> <mo>×</mo> <mrow> <mo>(</mo> <msub> <mi>P</mi> <mrow> <mi>m</mi> <mo>_</mo> <mi>max</mi> </mrow> </msub> <mo>-</mo> <msub> <mi>P</mi> <mrow> <mi>m</mi> <mo>_</mo> <mi>min</mi> </mrow> </msub> <mo>)</mo> </mrow> </mrow> <msubsup> <mi>g</mi> <mi>max</mi> <mn>2</mn> </msubsup> </mfrac> <msup> <mrow> <mo>(</mo> <mi>g</mi> <mo>-</mo> <mfrac> <msub> <mi>g</mi> <mi>max</mi> </msub> <mn>2</mn> </mfrac> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>10</mn> <mo>)</mo> </mrow> </mrow></math>
wherein g is evolution algebra, and g belongs to [1, g ∈max],Pm_min=Pb,Pm_max=10Pm_min;
PbFor the basic mutation probability, it can be estimated as follows:
<math><mrow> <msub> <mi>P</mi> <mi>b</mi> </msub> <mo>≈</mo> <mfrac> <mn>175</mn> <mrow> <mi>popsize</mi> <mo>×</mo> <msqrt> <mi>bits</mi> </msqrt> </mrow> </mfrac> <mo>;</mo> </mrow></math>
where bits represents the number of bits of an individual. When popsize is 20 and bits is 16, Pb≈0.02;
When P is presentb≈0.02,gmaxWhen 100, we can get:
Pm=Pm(g)=0.2-0.00008×(g-50)2;
termination criteria: when the algorithm is executed to the maximum algebra or the highest fitness value in the group is stable, the algorithm stops running, and the individual with the highest fitness value is the solution. At the termination of the algorithm, the individual with the highest fitness value is taken as the best threshold.
The basic steps of the partitioning two-dimensional maximum entropy threshold segmentation algorithm are as follows:
(1) segmenting the image into a series of sub-images;
(2) calculating a threshold value for each sub-image (each threshold value is calculated by a two-dimensional maximum entropy algorithm);
(3) and segmenting the image according to the threshold value of each sub-image.
The segmented image is shown in fig. 4.
4. Moving object tracking
A single-frame detection result can be determined to be obtained through an improved two-dimensional maximum entropy threshold segmentation algorithm, and a target can be approximated by a uniform linear motion model because foreign matters and bubbles move at a low speed in water.
Tracking system model: let the state vector be <math><mrow> <mi>X</mi> <mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo> </mrow> <mo>=</mo> <msup> <mrow> <mo>(</mo> <mi>x</mi> <mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo> </mrow> <mo>,</mo> <mover> <mi>x</mi> <mo>·</mo> </mover> <mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo> </mrow> <mo>,</mo> <mi>y</mi> <mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo> </mrow> <mo>,</mo> <mover> <mi>y</mi> <mo>·</mo> </mover> <mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo> </mrow> <mo>)</mo> </mrow> <mi>T</mi> </msup> <mo>,</mo> </mrow></math> The measurement vector is Z (k) ═ x (k), y (k)]TThe system equation can be written as:
X(k+1)=ΦX(k)+Gw(k) (11)
Z(k+1)=HX(k)+v(k) (12)
wherein <math><mrow> <mi>Φ</mi> <mo>=</mo> <mrow> <mfenced open='[' close=']'> <mtable> <mtr> <mtd> <mn>1</mn> </mtd> <mtd> <mi>T</mi> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mn>1</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mn>1</mn> </mtd> <mtd> <mi>T</mi> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mn>1</mn> </mtd> </mtr> </mtable> </mfenced> <mo>,</mo> </mrow> </mrow></math>
w (k), v (k) are zero-mean white Gaussian noise with variance Q (k) and R (k), respectively. Where T is the sampling period.
Track start: if a detected possible target appears in the 3 x 3 neighborhood for 3 consecutive frames, the track starts.
Tracking gates and data association: assume that the state vector at time k-1 is estimated as
The one-step prediction of the state vector is:
<math><mrow> <mover> <mi>X</mi> <mo>^</mo> </mover> <mrow> <mo>(</mo> <mi>k</mi> <mo>|</mo> <mi>k</mi> <mo>-</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>=</mo> <mi>Φ</mi> <mover> <mi>X</mi> <mo>^</mo> </mover> <mrow> <mo>(</mo> <mi>k</mi> <mo>-</mo> <mn>1</mn> <mo>|</mo> <mi>k</mi> <mo>-</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>13</mn> <mo>)</mo> </mrow> </mrow></math>
the innovation of the system is:
the innovation variance matrix is:
S(k)=HP(k|k-1)HT+R(k) (15)
where P (k | k-1) is the one-step prediction variance.
Let the norm of the innovation vector d (k) be:
g(k)=dT(k)S-1(k)d(k); (16)
wherein g (k) is subject to
Distribution, M1 is the measurement dimension.
An ellipsoid (also called tracking gate) is defined in the measurement space such that the measurements are distributed with a certain probability within the tracking gate
<math><mrow> <mover> <mi>V</mi> <mo>~</mo> </mover> <mrow> <mo>(</mo> <mi>γ</mi> <mo>)</mo> </mrow> <mo>=</mo> <mrow> <mo>[</mo> <mi>Z</mi> <mo>:</mo> <mi>g</mi> <mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo> </mrow> <mo>≤</mo> <mi>γ</mi> <mo>]</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>17</mn> <mo>)</mo> </mrow> </mrow></math>
Wherein γ can pass through x2And (6) checking a distribution table.
If only one echo exists in the tracking gate after a certain frame is processed, the target track is directly updated; if there is more than one echo in the tracking gate, the target trajectory is updated by the measurement closest to the one-step predicted value.
The track of the target is updated by a standard Kalman filtering algorithm
X(k|k-1)=ΦX(k-1|k-1)
P(k|k-1)=ΦP(k-1|k-1)ΦT+GQ(k-1)GT
K(k)=P(k|k-1)HT[HP(k|k-1)HT+R]-1 (18)
X(k|k)=X(k|k-1)+K(k)[Z(k)-H(k)X(k|k-1)]
P(k|k)=[I-K(k)H]P(k|k-1)
Fig. 5 shows a motion trail diagram obtained by target tracking.
After the target is detected and tracked, if the target disappears momentarily due to some optical factors, the next possible position of the target can be predicted according to the previous position information and motion state of the target, and when the target appears again, the target can still be stably tracked without losing the target. This is done to greatly reduce the possibility of misinterpreting noise as a foreign object. Because the movement speed of the foreign matters in the liquid is slow, the filtering value of the state vector of the target can be obtained through the target tracking based on the uniform linear motion model, and the possible position of the target can be obtained through extrapolation prediction.
Assume that the state vector at time k is estimated as
Then the predicted value of the target at k + n frame is:
<math><mrow> <mover> <mi>X</mi> <mo>^</mo> </mover> <mrow> <mo>(</mo> <mi>k</mi> <mo>+</mo> <mi>n</mi> <mo>|</mo> <mi>k</mi> <mo>)</mo> </mrow> <mo>=</mo> <msup> <mi>Φ</mi> <mi>n</mi> </msup> <mover> <mi>X</mi> <mo>^</mo> </mover> <mrow> <mo>(</mo> <mi>k</mi> <mo>|</mo> <mi>k</mi> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>19</mn> <mo>)</mo> </mrow> </mrow></math>
when applied to actual detection, the past more total tracking data is less and less relevant to the future situation along with the increase of time and the change of possible object motion, and the prediction precision is reduced along with the increase of n. Here n <6 is generally chosen.
5. Object recognition and determination
In this step, some image recognition techniques are mainly used, and for each assumed track obtained by target tracking, whether the track is formed by foreign matters or noise or bubbles is distinguished according to the following 3 principles, which are as follows:
(1) since the image is taken when the vial is stationary after being turned 180 degrees, and the visible foreign object is moving downward, the centroid ordinate of the image should be sequentially enlarged (with the upper left corner as the origin of coordinates);
(2) the bubbles move upwards in the opposite direction to the foreign bodies;
(3) the trace curve formed by the foreign matter is smooth, and the trace formed by the noise is disordered.
Therefore, if the centroids of the detected motion tracks become smaller in sequence, the motion tracks are indicated to be generated by bubbles; if the track that the motion track is smooth and the centroid vertical coordinate is sequentially enlarged is detected, the existence of foreign matters is indicated, the judgment result is recorded, and a rejection signal is sent to the system. The trajectories of foreign substances and air bubbles and the recognition results are shown in fig. 5.