Real-time recognition and positioning method and system for breast duct inner wall lesion based on optical fiber imaging

Abstract

The invention provides a real-time breast duct inner wall lesion recognition and positioning method and system based on optical fiber imaging, and relates to the technical field of medical image processing, and the method comprises the steps: obtaining a breast duct inner wall image sequence, extracting a displacement vector field, decomposing the displacement vector field into a dominant frequency and residual components, analyzing the dominant frequency phase to obtain a tissue motion period, an abnormal displacement area is screened from the residual error to establish a tissue anomaly map, an area descriptor is constructed through local wavelet coefficient analysis, a lesion core area is determined through density clustering, and finally the range and the expansion direction of a lesion area are determined through a multi-scale radial basis function and isoline analysis. According to the invention, real-time accurate identification and positioning of the lesion of the inner wall of the breast duct can be realized.

Description

TECHNICAL FIELD

[0001] The present application relates to medical image processing technology, and in particular to a real-time identification and positioning method and system for breast duct wall lesions based on optical fiber imaging. BACKGROUND

[0002] Traditional breast examination methods include mammography, ultrasonic examination and magnetic resonance imaging, etc. Although these methods can detect larger lesions, the identification effect for early micro breast duct wall lesions is limited. With the development of medical technology, optical fiber imaging technology is introduced into the field of breast disease diagnosis. Through a minimally invasive method, a miniature optical fiber probe is introduced into the breast duct to obtain high-resolution images of the breast duct wall in real time. This technology can directly observe the morphology and lesion of the breast duct wall, providing a new means for the diagnosis of early breast disease. Currently, optical fiber ductoscopy has been preliminarily applied in clinical practice, and can obtain image sequences of the breast duct wall to assist doctors in diagnosis.

[0003] The existing ductoscopy technology still has many deficiencies in practical application. The image sequences obtained by the ductoscopy are large in quantity and complex. Relying entirely on the naked eye of the doctor to identify the lesion area not only has low efficiency, but also is easily affected by subjective factors, causing missed diagnosis or misdiagnosis. The breast duct wall tissue has natural movement and deformation caused by probe movement during the examination process, so that the position of the lesion area in the image sequence changes constantly, making it difficult to accurately track and locate. The morphology of the breast duct wall lesion is diverse and the boundary is blurred. In particular, the distinction between early micro lesions and normal tissues is relatively subtle. The traditional recognition method based on morphological features cannot effectively distinguish the lesion area, resulting in insufficient diagnostic accuracy.

[0004] Therefore, it is of important clinical value to develop a method that can automatically and accurately identify and locate the breast duct wall lesion, so as to improve the efficiency and accuracy of early diagnosis of breast disease. SUMMARY

[0005] The embodiment of the present application provides a real-time identification and positioning method and system for breast duct wall lesions based on optical fiber imaging, which can solve the problems in the prior art.

[0006] In a first aspect, the embodiment of the present application provides a real-time identification and positioning method for breast duct wall lesions based on optical fiber imaging, comprising: obtaining breast duct wall image sequences and dividing them into multiple subsequences according to a preset time interval; extracting the displacement vector field between adjacent images in each subsequence, decomposing the displacement vector field into a main frequency component and a residual component in the frequency domain, obtaining the tissue motion period through phase analysis of the main frequency component, screening abnormal displacement regions that do not conform to the tissue motion period from the residual component, merging the abnormal displacement regions in multiple subsequences, and establishing a tissue abnormality atlas; The tissue abnormality map is divided into blocks, and the local wavelet coefficients in each block area are calculated. The local wavelet coefficients are grouped and statistically analyzed to obtain internal structural features. The coefficient differences between adjacent blocks are calculated to obtain structural change features. The internal structural features and structural change features are combined to form a regional descriptor. The density clustering method is used to analyze the regional descriptor to determine the core area of ​​the lesion. A multi-scale radial basis function is constructed with the core area of ​​the lesion as the center. The response distribution of the regional descriptor under the multi-scale radial basis function is calculated, and the spatial gradient contour lines of the response distribution are extracted. The spatial range of the lesion area is determined based on the contour lines, and the evolution of the contour lines is tracked by optical flow field analysis to determine the direction of lesion expansion.

[0007] In an optional embodiment, Acquiring a sequence of images of the inner wall of the milk duct and dividing the sequence into multiple subsequences according to preset time intervals includes: Acquiring a continuous image sequence of the inner wall of the milk duct, performing time domain analysis on the continuous image sequence to extract a periodic signal of the image sequence, setting a time segmentation benchmark according to the periodic signal, and dividing the continuous image sequence into initial subsequences using an integer multiple of the time segmentation benchmark as a subsequence duration; An overlapping interval is set between adjacent sequences of the initial subsequence, and the duration of the overlapping interval is a preset proportion of the subsequence duration; a weight value is calculated based on the temporal distance between the image and the boundary of the overlapping interval, and the images in the overlapping interval are weighted to obtain a transition sequence in the overlapping interval. The transition sequence replaces the image sequence in the original overlapping interval to generate multiple temporally continuous image subsequences.

[0008] In an optional embodiment, Extract the displacement vector field between adjacent images in each subsequence and decompose the displacement vector field into the main frequency component and the residual component in the frequency domain, including: Obtaining the grayscale values ​​of adjacent image frames in a subsequence and constructing an image block centered on a pixel point; establishing grayscale constancy constraints and spatial consistency constraints within the image block to form a local grayscale change constraint matrix; substituting the local grayscale change constraint matrix into the optical flow equation to obtain a pixel-level displacement vector; and performing regional integration on the pixel-level displacement vector to obtain a displacement vector field; Decomposing the displacement vector field using wavelet basis functions to obtain displacement components of different scales; calculating the amplitude distribution of the displacement components of each scale, generating directional weight coefficients, fusing the displacement components according to the directional weight coefficients, and reconstructing the characteristics of the displacement vector field; performing two-dimensional Fourier transform on the displacement vector field features to obtain a frequency domain representation, counting energy distribution of each frequency component in the frequency domain representation, and calculating cumulative energy value of the frequency component; determining a frequency segmentation threshold value according to the cumulative energy value, dividing the frequency domain representation by the frequency segmentation threshold value, and reconstructing by inverse Fourier transform to obtain a main frequency component and a residual component.

[0009] In an alternative embodiment, obtaining a tissue motion cycle by phase analysis of the main frequency component, screening abnormal displacement regions from the residual component that do not conform to the tissue motion cycle, merging abnormal displacement regions in multiple sub-sequences, and establishing a tissue abnormality atlas including: performing Hilbert transform on the main frequency component to extract instantaneous phase, counting time sequence variation of the instantaneous phase, identifying extreme points in the phase sequence, calculating time interval sequence between adjacent extreme points, and determining the tissue motion cycle and its fluctuation range according to the time interval sequence; dividing grid points in the residual component, calculating corner point response values of the grid points, selecting grid points with corner point response values greater than a set response threshold value as feature points, calculating cross-correlation coefficients of the feature points between adjacent frames for feature point matching, and obtaining position displacement of the feature points based on the matching results; concatenating and filtering the position displacement along the time dimension to obtain a motion trajectory of the feature points, calculating tangent vectors and normal vectors of the motion trajectory, obtaining arc length variation and curvature variation of the trajectory, comparing the arc length variation and the curvature variation with the tissue motion cycle, marking trajectory points that exceed the cycle fluctuation range, and extracting a region containing the marked points as an abnormal displacement region; calculating motion direction, speed, and acceleration features of the abnormal displacement region, constructing motion features, calculating similarity between abnormal displacement regions according to the motion features, performing connected component analysis on abnormal displacement regions with similarity greater than a preset similarity threshold value and adjacent in space, and integrating connected components into a tissue abnormality atlas.

[0010] In an alternative embodiment, dividing the tissue abnormality atlas into blocks, calculating local wavelet coefficients in each block region, grouping and counting the local wavelet coefficients to obtain internal structure features, calculating coefficient difference degree between adjacent blocks to obtain structure change features, and combining the internal structure features and the structure change features to form a region descriptor including: performing recursive quadtree segmentation on the tissue abnormality atlas, calculating information entropy of pixel gray scale distribution in a block region, determining whether to continue to subdivide the current block according to a comparison result of the information entropy and a segmentation threshold value, until all blocks meet a stop condition, and obtaining an adaptive block atlas; Perform wavelet decomposition on each block area in the adaptive block atlas to obtain detail coefficient matrices in the horizontal, vertical and diagonal directions; calculate the grayscale gradient distribution of the block area, generate directional weight coefficients, calculate statistics on the detail coefficient matrix and weight it using the directional weight coefficients to construct the internal structure feature vector of the block; Determine block adjacency based on the adaptive block map, perform scale normalization on adjacent block pairs, calculate the difference of local wavelet coefficients, construct a structural change matrix based on the difference, calculate the propagation characteristics of the structural change matrix at the region boundary, and perform cumulative statistics on the propagation characteristics to generate structural change characteristics; The internal structure features and the structural change features are aligned by dimension, a covariance matrix is ​​calculated for the aligned features, feature orthogonal decomposition is performed based on the covariance matrix, feature components are sorted and screened according to the decomposition results, and the screened feature components are combined to construct a region descriptor.

[0011] In an optional embodiment, The density clustering method is used to analyze the regional descriptors to determine the core areas of the lesions, including: The region descriptor is normalized to a preset numerical range through maximum and minimum value mapping to obtain a normalized region descriptor. The search radius is determined according to the feature distribution of the normalized region descriptor. The neighbor set of each normalized region descriptor is obtained within the search radius, and the distance matrix is ​​calculated based on the neighbor set. Determining a bandwidth parameter of an adaptive kernel function according to the median of the distance matrix, calculating a local density value of a normalized region descriptor using the adaptive kernel function, determining a target region descriptor according to the local density value, calculating a minimum distance value from each normalized region descriptor to the target region descriptor as a relative distance, and constructing a density-distance distribution map according to the local density value and the relative distance; Extracting statistical features of the density-distance distribution map, setting an adaptive density threshold and a distance threshold based on the statistical features, and screening normalized region descriptors to obtain initial cluster centers; The density connectivity between the initial cluster centers is calculated, and the spatially adjacent initial cluster centers are merged based on the density connectivity. The cluster center with the maximum local density value and the strongest density connectivity after the merger is determined as the location of the lesion core area.

[0012] In an optional embodiment, A multi-scale radial basis function is constructed with the lesion core area as the center. The response distribution of the regional descriptor under the multi-scale radial basis function is calculated. The spatial gradient contour of the response distribution is extracted. The spatial range of the lesion area is determined based on the contour. The evolution of the contour is tracked by optical flow field analysis to determine the direction of lesion expansion. Taking a center position of the lesion core region as a reference point, a reference scale value is calculated according to a spatial range of the lesion core region, an incremental scale sequence is constructed according to a preset ratio based on the reference scale value, and a multi-scale radial basis function is constructed by using the incremental scale sequence; A response distribution of the region descriptor under the multi-scale radial basis function is calculated, an anisotropic weight matrix is generated based on a structural feature gradient of the region descriptor, the anisotropic weight matrix is multiplied by the response distribution to obtain a modified response distribution; A spatial gradient contour line of the modified response distribution is extracted at each scale, a stability index is calculated according to a spatial overlap ratio of the contour lines at adjacent scales, a change index is extracted according to a response gradient of a neighborhood of the contour line, and a consistency index is calculated according to a morphological feature difference between the contour line and the lesion core region; The spatial gradient contour line is optimized based on the stability index, the change index and the consistency index, and a maximum closed contour line after optimization is determined as a spatial range of the lesion region; An optical flow field vector is calculated for a boundary of the spatial range of the lesion region, principal component analysis is performed on the optical flow field vector, and a feature vector direction corresponding to a maximum eigenvalue obtained by the principal component analysis is determined as an extension direction of the lesion region.

[0013] In a second aspect of the embodiment, a real-time identification and positioning system for intraductal wall lesions based on optical fiber imaging is provided, and the system comprises: A first unit is configured to acquire an intraductal wall image sequence and divide the sequence into a plurality of subsequences according to a preset time interval; A second unit is configured to extract a displacement vector field between adjacent images in each subsequence, decompose the displacement vector field into a main frequency component and a residual component in a frequency domain, obtain a tissue motion cycle through phase analysis of the main frequency component, screen abnormal displacement regions that do not conform to the tissue motion cycle from the residual component, merge the abnormal displacement regions in the plurality of subsequences, and establish a tissue abnormality atlas; A third unit is configured to block the tissue abnormality atlas, calculate local wavelet coefficients in each block region, group and statistically analyze the local wavelet coefficients to obtain internal structure features, calculate a coefficient difference degree between adjacent blocks to obtain structure change features, and combine the internal structure features and the structure change features to form a region descriptor; A fourth unit is configured to analyze the region descriptor by using a density clustering method to determine a lesion core region, construct a multi-scale radial basis function with the lesion core region as a center, calculate a response distribution of the region descriptor under the multi-scale radial basis function, extract a spatial gradient contour line of the response distribution, determine a spatial range of a lesion region based on the contour line, and determine a lesion extension direction by tracking evolution of the contour line through optical flow field analysis.

[0014] According to a third aspect of an embodiment of the present invention, an electronic device is provided, including: processor; a memory for storing processor-executable instructions; The processor is configured to call the instructions stored in the memory to execute the aforementioned method.

[0015] According to a fourth aspect of an embodiment of the present invention, a computer-readable storage medium is provided, on which computer program instructions are stored. When the computer program instructions are executed by a processor, the method described above is implemented.

[0016] In this embodiment, a real-time identification and location method for ductal lesions based on fiber optic imaging can dynamically analyze ductal lesion image sequences, identify abnormal motion regions through frequency domain decomposition of the displacement vector field, effectively capture subtle changes in the characteristics of potential lesions, and improve the sensitivity of early detection of ductal lesions. By performing local wavelet analysis and regional descriptor construction on tissue abnormality maps, accurate characterization of different types of ductal lesions is achieved. Combined with density clustering technology, the core area of ​​the lesion is automatically identified, avoiding the subjectivity of manual interpretation in traditional methods and significantly improving the accuracy and objectivity of lesion identification. The present invention uses multi-scale radial basis functions and spatial gradient contouring technology to accurately depict the spatial extent and expansion direction of the lesion area, providing clinicians with three-dimensional location information of the lesion, helping to formulate more precise treatment plans and reduce unnecessary tissue damage. It has important clinical application value for the early diagnosis and minimally invasive treatment of breast diseases. BRIEF DESCRIPTION OF THE DRAWINGS

[0017] Figure 1 Schematic diagram of the process of a method for real-time identification and positioning of milk duct inner wall lesions based on optical fiber imaging according to an embodiment of the present invention; Figure 2 Schematic diagram of feature point motion trajectory analysis and anomaly detection according to an embodiment of the present invention; Figure 3 This is a comparison chart of the lesion analysis effects of the embodiments of the present invention. DETAILED DESCRIPTION

[0018] To make the objectives, technical solutions, and advantages of the embodiments of the present invention more clear, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts shall fall within the scope of protection of the present invention.

[0019] The following specific embodiments are used to describe the technical solution of the present invention in detail. The following specific embodiments can be combined with each other, and the same or similar concepts or processes may not be described in detail in some embodiments.

[0020] Figure 1 FIG. 1 is a flow chart of a method for real-time identification and positioning of milk duct inner wall lesions based on optical fiber imaging according to an embodiment of the present invention. Figure 1 As shown, the method includes: Acquire a sequence of images of the inner wall of the milk duct and divide the sequence into multiple subsequences according to preset time intervals; The displacement vector field between adjacent images in each subsequence is extracted and decomposed into a primary frequency component and a residual component in the frequency domain. The tissue motion cycle is obtained by phase analysis of the primary frequency component. Abnormal displacement regions that do not conform to the tissue motion cycle are screened out from the residual component, and the abnormal displacement regions in multiple subsequences are merged to establish a tissue abnormality map. The tissue abnormality map is divided into blocks, and the local wavelet coefficients in each block area are calculated. The local wavelet coefficients are grouped and statistically analyzed to obtain internal structural features. The coefficient differences between adjacent blocks are calculated to obtain structural change features. The internal structural features and structural change features are combined to form a regional descriptor. The density clustering method is used to analyze the regional descriptor to determine the core area of ​​the lesion. A multi-scale radial basis function is constructed with the core area of ​​the lesion as the center. The response distribution of the regional descriptor under the multi-scale radial basis function is calculated, and the spatial gradient contour lines of the response distribution are extracted. The spatial range of the lesion area is determined based on the contour lines, and the evolution of the contour lines is tracked by optical flow field analysis to determine the direction of lesion expansion.

[0021] In an optional embodiment, acquiring a milk duct inner wall image sequence and dividing it into a plurality of subsequences according to preset time intervals includes: Acquiring a continuous image sequence of the inner wall of the milk duct, performing time domain analysis on the continuous image sequence to extract a periodic signal of the image sequence, setting a time segmentation benchmark according to the periodic signal, and dividing the continuous image sequence into initial subsequences using an integer multiple of the time segmentation benchmark as a subsequence duration; An overlapping interval is set between adjacent sequences of the initial subsequence, and the duration of the overlapping interval is a preset proportion of the subsequence duration; a weight value is calculated based on the temporal distance between the image and the boundary of the overlapping interval, and the images in the overlapping interval are weighted to obtain a transition sequence in the overlapping interval. The transition sequence replaces the image sequence in the original overlapping interval to generate multiple temporally continuous image subsequences.

[0022] In this embodiment, the medical device first uses a fiber optic endoscope to acquire a continuous sequence of images of the inner wall of the milk duct. In one specific embodiment, the fiber optic endoscope has a diameter of 0.5 mm, allowing for smooth entry into the milk duct for imaging. The endoscope captures images at a rate of 30 frames per second, recording the tissue state of the milk duct inner wall. These captured continuous images constitute a complete image sequence, encompassing comprehensive information about the inner wall of the milk duct.

[0023] Time-domain analysis is performed on the acquired continuous image sequence, detecting features such as pixel brightness changes and contour line changes to extract periodic signals from the image sequence. In practice, imaging of the inner wall of the milk duct is affected by physiological activities such as the patient's breathing and heartbeat, resulting in periodic image changes. For example, in one acquisition, the patient's breathing cycle was detected to be approximately 4 seconds, and the heart cycle was approximately 0.8 seconds. These periodic changes can affect image stability.

[0024] The time segmentation benchmark is set based on the extracted periodic signal. In practice, the respiratory cycle is the primary reference, as respiration significantly influences tissue movement. Assuming a detected respiratory cycle of 4 seconds, 4 seconds can be used as the time segmentation benchmark. The subsequence duration is determined by using integer multiples of this time segmentation benchmark. For example, a 1x benchmark, or 4 seconds, is used as the subsequence duration. This allows the continuous image sequence to be segmented every 4 seconds, forming multiple initial subsequences. For an acquisition frequency of 30 frames per second, each subsequence contains 120 frames.

[0025] To ensure smooth transitions between adjacent subsequences, an overlap interval is set between adjacent sequences within the initial subsequence. The duration of the overlap interval is a preset proportion of the subsequence duration, for example, 25% of the subsequence duration, or 1 second. This way, two adjacent subsequences overlap for 1 second, corresponding to 30 overlapping frames.

[0026] Weights are calculated based on the temporal distance between the image and the boundary of the overlapping interval. Within the overlapping interval, images near the end of the previous subsequence have higher weights, while images near the beginning of the next subsequence have lower weights, resulting in a gradual transition. Specifically, for the i-th frame within the overlapping interval, if there are n frames in the overlapping interval, the weight of that image can be calculated based on the distance relationship. For example, the weight of the i-th frame of the previous subsequence within the overlapping interval can be set to (ni) / n, and the weight of the i-th frame of the next subsequence within the overlapping interval can be set to i / n. Taking 30 overlapping images as an example, the weight of the first frame is 29 / 30 for the previous sequence and 1 / 30 for the next sequence; the weight of the 15th frame is 15 / 30 for the previous sequence and 15 / 30 for the next sequence; and the weight of the 30th frame is 0 / 30 for the previous sequence and 30 / 30 for the next sequence.

[0027] The images within the overlapping interval are weighted, and the corresponding images in the preceding and following subsequences are fused according to the weights to generate a transition sequence. During weighted processing, each pixel in the image is calculated based on the pixel values ​​of the corresponding images in the preceding and following subsequences and their respective weights to obtain a weighted pixel value. For example, for the i-th frame of the image within the overlapping interval, if the pixel value of the corresponding image in the preceding subsequence is P1, and the pixel value at the same position in the corresponding image in the following subsequence is P2, and the weight of the preceding subsequence is w1 and the weight of the following subsequence is w2, then the weighted pixel value is P1 multiplied by w1 plus P2 multiplied by w2. This pixel-level weighted fusion achieves a smooth image transition.

[0028] The resulting transition sequence replaces the original image sequence within the overlapping interval, generating multiple temporally continuous image subsequences. Ultimately, the entire continuous image sequence is divided into multiple subsequences with smooth transition characteristics. Each subsequence has the same length, and adjacent subsequences are weighted to achieve smooth transitions, avoiding the sudden changes that may be caused by simple splicing.

[0029] In practice, the system dynamically adjusts parameters based on actual conditions. For example, the respiratory cycle may vary from patient to patient, so the system adjusts the time segmentation benchmark based on the real-time detected respiratory cycle. The overlap ratio can also be adjusted based on the specific examination site, typically selecting an appropriate value between 15% and 35%.

[0030] In this embodiment, by acquiring a continuous image sequence and performing periodic signal analysis, this method can effectively eliminate interference from physiological activities such as breathing and heartbeat on imaging, thereby improving image stability. Adaptive adjustment of the time segmentation benchmark ensures that the subsequence division conforms to the characteristics of the physiological cycle, ensuring the accuracy of image analysis. Overlapping interval setting and weighted processing technology achieve smooth transitions between subsequences, avoiding image mutations caused by simple splicing and enhancing sequence coherence. Pixel-level weighted fusion ensures that the image quality in the transition region does not degrade, maintaining the information integrity of the original image. This method can adapt to the differences in physiological characteristics of different patients and improves adaptability and robustness through dynamic parameter adjustment. The resulting high-quality image subsequence lays the foundation for subsequent lesion identification and localization, significantly improving the detection rate and localization accuracy of intraductal lesions, enabling doctors to make more accurate diagnoses and treatment plans, and is of great value in improving the early diagnosis rate and treatment effectiveness of breast diseases.

[0031] In an optional embodiment, extracting the displacement vector field between adjacent images in each subsequence and decomposing the displacement vector field into a primary frequency component and a residual component in the frequency domain includes: Obtaining the grayscale values ​​of adjacent image frames in a subsequence and constructing an image block centered on a pixel point; establishing grayscale constancy constraints and spatial consistency constraints within the image block to form a local grayscale change constraint matrix; substituting the local grayscale change constraint matrix into the optical flow equation to obtain a pixel-level displacement vector; and performing regional integration on the pixel-level displacement vector to obtain a displacement vector field; Decomposing the displacement vector field using wavelet basis functions to obtain displacement components of different scales; calculating the amplitude distribution of the displacement components of each scale, generating directional weight coefficients, fusing the displacement components according to the directional weight coefficients, and reconstructing the characteristics of the displacement vector field; A two-dimensional Fourier transform is performed on the displacement vector field characteristics to obtain a frequency domain representation, the energy distribution of each frequency component in the frequency domain representation is statistically analyzed, and the cumulative energy value of the frequency component is calculated; a frequency segmentation threshold is determined based on the cumulative energy value, the frequency domain representation is divided with the frequency segmentation threshold as the boundary and reconstructed through an inverse Fourier transform to obtain the main frequency component and the residual component.

[0032] In a specific embodiment, the grayscale values ​​of adjacent image frames in a subsequence are first obtained, and the color image is converted to a grayscale image to simplify calculations. In practical applications, images of the inner wall of milk ducts acquired using a fiber optic endoscope are typically RGB color images. These images can be converted to grayscale images by weighting the R channel value by 0.299, the G channel value by 0.587, and the B channel value by 0.114. For an original endoscopic image with a resolution of 640×480, the conversion yields a grayscale image of the same resolution, with a grayscale value range of 0 to 255.

[0033] An image block is constructed with each pixel as the center, and the block size is selected to be N×N pixels. In practical applications, the value of N is typically set to an odd number between 5 and 15, for example, an 11×11 image block can be selected. This size allows for sufficient local information without incurring excessive computational burden. The system establishes a grayscale constancy constraint within each image block. This constraint is based on the assumption that the grayscale value of the same tissue point in adjacent images should remain unchanged over a short period of time. For milk duct lining tissue, this assumption holds true in most cases, unless there is significant illumination change or tissue deformation.

[0034] At the same time, a spatial consistency constraint is established within the image block. This constraint is based on the assumption that the motion of adjacent pixels should be similar and not drastically change. This constraint is particularly applicable to milk duct lining tissue, as milk duct tissue typically moves as a whole, with adjacent regions often moving in similar directions and magnitudes. The grayscale constrains and spatial consistency constraints are combined to form a local grayscale variation constraint matrix. For an 11×11 image block, this matrix describes the grayscale variation relationships and spatial constraints of 121 pixels between two adjacent frames.

[0035] The local gray level change constraint matrix is substituted into the optical flow equation for solving. The optical flow equation expresses the relationship between the change of image gray value over time and spatial gradient, and the movement vector of the pixel can be obtained by solving the equation. In practical applications, an iterative solving method is used, and the number of iterations is usually set to 3 to 5 times. Each iteration uses the result of the last iteration for optimization to obtain a more accurate solution. In this way, the horizontal and vertical displacements of each pixel in the image are calculated to form a pixel-level displacement vector.

[0036] In order to reduce the influence of noise and improve stability, the pixel-level displacement vector is integrated by region. The specific implementation is to divide the image into a grid, for example, a 640x480 image is divided into a 32x24 grid, and each grid contains 20x20 pixels. In each grid, the average value of the displacement vectors of all pixel points is calculated as the overall displacement vector of the grid. Through this regional integration, a more stable displacement vector field is formed for subsequent analysis.

[0037] The displacement vector field is decomposed using wavelet basis functions to obtain displacement components of different scales. In practical applications, Haar wavelet, Daubechies wavelet or Symlet wavelet can be selected as the basis function. For example, using a 4th order Daubechies wavelet for 3 layer decomposition, the displacement vector field can be decomposed into a low frequency approximation component and 3 high frequency detail components of different scales. These components reflect the characteristics of large-scale overall movement and small-scale local deformation of the breast duct tissue.

[0038] The amplitude distribution of each scale displacement component is calculated to generate direction weight coefficients. Specifically, for each scale displacement component, its energy distribution in the horizontal, vertical and diagonal directions is calculated. For example, if the horizontal direction displacement energy of a certain scale accounts for 60%, the vertical direction accounts for 30%, and the diagonal direction accounts for 10%, the corresponding direction weight coefficients are set to 0.6, 0.3 and 0.1 respectively. The system fuses the displacement components according to these direction weight coefficients to reconstruct the displacement vector field features. In the fusion process, displacement components of different directions are combined according to the weight, highlighting the contribution of the main motion direction and suppressing the interference of the secondary direction.

[0039] The reconstructed displacement vector field features are subjected to two-dimensional Fourier transform to convert them from spatial domain to frequency domain representation. After Fourier transform, the component distribution of the displacement vector field at different frequencies is obtained. In the frequency domain representation, the low frequency part corresponds to the overall slow movement of the breast duct tissue, and the high frequency part corresponds to the local fine changes or noise.

[0040] The energy distribution of each frequency component in the frequency domain representation is statistically analyzed, and the cumulative energy of the frequency components is calculated. In practice, the proportion of the cumulative energy to the total energy is calculated in ascending order of frequency. For example, the first 10% of low-frequency components may contain 85% of the total energy, indicating that the movement of the milk duct tissue is mainly concentrated in the low-frequency region, while the high-frequency region has less energy.

[0041] The frequency segmentation threshold is determined based on the cumulative energy value. Typically, the frequency point where the cumulative energy reaches 80% to 90% of the total energy is used as the segmentation threshold. In the above example, if the cumulative energy of the top 10% of frequency components reaches 85% of the total energy, 10% can be set as the frequency segmentation threshold. Using this frequency segmentation threshold as the boundary, the frequency domain representation is divided into low-frequency and high-frequency components. An inverse Fourier transform is performed on the divided low-frequency and high-frequency components, and then converted back to the spatial domain to obtain the primary frequency component and residual component of the displacement vector field. The primary frequency component corresponds to the low-frequency component and reflects the primary movement pattern of the milk duct tissue, such as the overall movement caused by breathing. The residual component corresponds to the high-frequency component and reflects the deformation of the local tissue or the characteristics of potential lesions.

[0042] In the identification of lesions on the inner wall of milk ducts, the main frequency component and the residual component have different application values. The main frequency component can be used to compensate for the overall tissue movement, ensuring that subsequent analysis focuses on the changes in the tissue itself rather than the artifacts caused by the movement. For example, the system can use the displacement field calculated by the main frequency component to align adjacent images and eliminate the influence of overall movement caused by factors such as breathing and heartbeat. The residual component contains the characteristic information of the potential lesion area, because the lesion tissue often exhibits different motion characteristics from the surrounding normal tissue. For example, the hardness of lesion tissue such as intraductal papilloma is usually different from that of normal tissue, and it exhibits different deformation patterns during the movement of the endoscope. These differences will be reflected in the residual component.

[0043] In this embodiment, pixel-level displacement vector calculation and regional integration are used to accurately capture the subtle motion characteristics of the ductal lining tissue, overcoming the difficulty traditional methods have in processing image changes caused by tissue motion. Wavelet basis function decomposition technology effectively separates tissue motion at different scales, enabling the system to simultaneously analyze global movement and local deformation characteristics. The introduction of directional weight coefficients enhances the characteristic representation of the primary motion direction, suppresses interference from secondary directions, and improves the accuracy of displacement field reconstruction. Frequency domain decomposition technology successfully separates tissue motion into primary and residual components, enabling the system to distinguish physiological motion from abnormal motion patterns associated with lesions. This decomposition method significantly improves the sensitivity and specificity of lesion identification and reduces the false positive detection rate. Furthermore, this technology can adapt to complex lighting conditions and tissue deformation within the duct, maintaining high robustness, providing reliable technical support for the early and accurate identification and localization of ductal lining lesions.

[0044] In an optional embodiment, the tissue motion cycle is obtained by phase analysis of the dominant frequency component; the abnormal displacement regions not conforming to the tissue motion cycle are screened out from the residual component, the abnormal displacement regions in multiple sub-sequences are merged, and the tissue abnormal atlas is established, comprising: The Hilbert transform is performed on the dominant frequency component to extract the instantaneous phase, the time sequence change of the instantaneous phase is counted, the extreme points in the phase sequence are identified, the time interval sequence between adjacent extreme points is calculated, and the tissue motion cycle and its fluctuation range are determined according to the time interval sequence; The grid points are divided in the residual component, the corner point response value of the grid points is calculated, the grid points with a corner point response value greater than a set response threshold are selected as feature points, the cross-correlation coefficient between adjacent frames of the feature points is calculated for feature point matching, and the position offset of the feature points is obtained based on the matching result; The position offset is concatenated along the time dimension and filtered to obtain the motion trajectory of the feature points, the tangent vector and normal vector of the motion trajectory are calculated, the arc length change and curvature change of the trajectory are obtained, the arc length change and curvature change are compared with the tissue motion cycle, and the trajectory points exceeding the cycle fluctuation range are marked, and the region containing the marked points is extracted as an abnormal displacement region; The motion direction, speed and acceleration characteristics of the abnormal displacement region are calculated, the motion characteristics are constructed, the similarity between the abnormal displacement regions is calculated according to the motion characteristics, the abnormal displacement regions with a similarity greater than a preset similarity threshold and adjacent in space are subjected to connected component analysis, and the connected components are integrated into a tissue abnormal atlas.

[0045] In implementation, first, the Hilbert transform is performed on the dominant frequency component of the displacement vector field to extract its instantaneous phase information. The Hilbert transform can convert a real signal into an analytic signal, and the instantaneous phase can be obtained through the phase angle of the analytic signal. In actual application, the Hilbert transform can be performed on the horizontal and vertical displacement components in the dominant frequency component to obtain the instantaneous phase in two directions. For the breast duct wall tissue, due to the influence of respiration and heartbeat, the dominant frequency component usually presents periodic changes, which can be reflected by the change trend of the instantaneous phase.

[0046] The time sequence change of the instantaneous phase is counted, and the extreme points in the phase sequence are identified. The extreme points in the phase sequence correspond to the turning points of the tissue motion direction, and are an important basis for cycle judgment. In the breast duct wall image sequence obtained by the optical fiber endoscope, the extreme points can be identified by differentiating the phase curve and finding the zero points. In actual operation, the three-point central difference method can be used to calculate the derivative of the phase curve, and then the positions of the sign change of the derivative are detected and marked as extreme points. For example, in a 30-second breast duct wall image sequence, about 75 phase extreme points can be detected, corresponding to the respiration and heartbeat cycles of the patient.

[0047] The time interval sequence between adjacent extreme points is calculated, and the tissue motion cycle and its fluctuation range are determined based on this time interval sequence. During imaging of the milk duct wall, respiration and heartbeat are the two primary periodic factors influencing tissue motion. These two periods can be identified by analyzing the distribution characteristics of the time interval sequence. For example, in adult patients, the respiratory cycle typically ranges from 3 to 5 seconds, and the cardiac cycle typically ranges from 0.6 to 1.2 seconds. During statistical analysis, a histogram or density estimation method can be used to identify the peak of the time interval distribution as the primary period of tissue motion. Simultaneously, the standard deviation of each periodic component is calculated to determine the period's fluctuation range. For example, the respiratory cycle might be identified as 4.2 seconds with a fluctuation range of ±0.5 seconds, and the cardiac cycle might be identified as 0.8 seconds with a fluctuation range of ±0.1 seconds.

[0048] Grid points are assigned to the residual component, and corner response values ​​are calculated for each grid point. The residual component reflects irregular local tissue motion and may contain potential pathological information. To effectively analyze the residual component, the image can be evenly divided into a grid, for example, with a grid point every 10 pixels. For each grid point, the grayscale gradient covariance matrix of the surrounding area is calculated, and the corner response value is calculated based on this matrix. The corner response value reflects the texture richness of the local area, with larger values ​​indicating richer feature information in the area.

[0049] Grid points whose corner response values ​​exceed a set response threshold are selected as feature points. The response threshold can be set dynamically based on the overall image characteristics. For example, the threshold can be set to the mean of all grid point corner response values ​​plus 1.5 times the standard deviation. For milk duct wall images, feature points are typically distributed in areas with rich tissue texture or distinct edges. In a 640×480 resolution milk duct wall image, 100 to 200 feature points may be selected for subsequent analysis.

[0050] The cross-correlation coefficients of feature points between adjacent frames are calculated for feature point matching, and the position offset of the feature points is obtained based on the matching results. In specific implementation, a fixed-size image block, such as a 15×15 pixel area, can be extracted around the feature point, and the cross-correlation coefficient between this image block and the possible corresponding position in the next frame can be calculated. The cross-correlation coefficient reflects the degree of similarity between two image blocks, with values ​​closer to 1 indicating a higher degree of match. The corresponding position of the feature point can be determined by searching for the position with the largest cross-correlation coefficient in the next frame, thereby calculating the position offset of the feature point. To improve matching accuracy, a pyramid matching strategy can be adopted, first performing coarse matching on a low-resolution image, followed by precise positioning on a high-resolution image.

[0051] The position offsets are concatenated along the time dimension and filtered to obtain the motion trajectory of the feature points. Because noise may be introduced during image acquisition and feature point matching, the raw offsets need to be filtered. Median filtering or Gaussian filtering can be used to eliminate abnormal jumps in the trajectory. The filtered trajectory reflects the actual movement of the milk duct lining tissue, providing a reliable basis for subsequent analysis.

[0052] The tangent vector and normal vector of the motion trajectory are calculated to determine the change in arc length and curvature of the trajectory. The tangent vector reflects the direction of motion of the feature point, while the normal vector reflects the change in direction of motion. Based on the tangent vector, the change in arc length (the distance the feature point moves per unit time) can be calculated; based on the normal vector, the change in curvature (the degree of curvature of the motion trajectory) can be calculated. Under the influence of respiration and heartbeat, the motion trajectory of normal milk duct tissue typically exhibits regular changes in arc length and curvature.

[0053] Compare changes in arc length and curvature with the tissue motion cycle, and mark trajectory points that exceed the cyclic fluctuation range. Normal tissue movement should be synchronized with the respiratory and cardiac cycles, and its arc length and curvature changes should be within the cyclic fluctuation range. If the arc length change or curvature change of a feature point's motion trajectory at a specific moment significantly exceeds the cyclic fluctuation range, it may indicate an abnormality in that area. For example, if the arc length change cycle of a feature point is 4.3 seconds, which is consistent with the identified respiratory cycle of 4.2 seconds ± 0.5 seconds, but suddenly changes to 6.2 seconds within a certain period of time, this segment of the trajectory is marked as abnormal.

[0054] Regions containing marked points are extracted as abnormal displacement regions. Spatial clustering methods can be used to cluster similarly located and abnormal feature points to form abnormal displacement regions. For example, a density clustering algorithm can be used, with a spatial distance threshold of 20 pixels. Anomalous feature points with a distance less than the threshold are grouped into the same region. For each abnormal displacement region, information such as its spatial extent, time span, and the number of abnormal feature points contained is recorded.

[0055] Calculate the motion direction, velocity, and acceleration characteristics of the abnormal displacement area to construct the motion characteristics. The motion direction can be obtained by calculating the average direction of the displacement vector of the feature points in the area; the velocity can be obtained by calculating the displacement per unit time; and the acceleration can be obtained by calculating the rate of change of the velocity. These features together constitute the motion feature description of the abnormal displacement area. For example, an abnormal displacement area may appear to be moving toward the nipple, with an average velocity of 2.5 pixels / frame and an acceleration of 0.3 pixels / frame. 2 .

[0056] Calculate the similarity between abnormal displacement regions based on motion characteristics. Methods such as Euclidean distance or cosine similarity can be used to calculate the distance between different abnormal displacement regions in the motion feature space. Smaller distances indicate higher similarity. For example, the motion direction, velocity, and acceleration can be normalized to construct feature vectors, and the Euclidean distance between these vectors can be calculated as a similarity metric.

[0057] Connected domain analysis is performed on spatially adjacent abnormal displacement regions with a similarity greater than a preset similarity threshold. The preset similarity threshold can be set to 0.8 or higher to ensure that the merged regions have highly similar motion characteristics. Connected domain analysis can use region growing or morphological processing methods to connect adjacent regions that meet the conditions to form larger connected regions.

[0058] Connected domains are integrated into a tissue abnormality atlas, preserving the spatial extent, temporal span, and motion characteristics of each connected domain during the integration process to form a complete abnormality atlas. This atlas visually displays areas of potential pathology within the milk duct lining, providing doctors with important diagnostic insights. For example, the abnormality atlas may reveal a location within a milk duct that consistently exhibits a motion pattern that differs from normal tissue, potentially indicating a pathology such as intraductal papilloma or intraductal carcinoma.

[0059] In this embodiment, by performing a Hilbert transform on the dominant frequency component to extract the instantaneous phase, this approach accurately identifies the periodic motion characteristics of the milk duct inner wall tissue, overcoming the drawback of existing techniques that rely solely on image grayscale changes to distinguish physiological from pathological motion. Existing techniques typically use fixed thresholds or empirical parameters to determine tissue abnormalities, which can be easily affected by physiological activities such as breathing and heartbeat, leading to misdiagnosis. This approach establishes a dynamic reference standard by calculating the tissue motion period and fluctuation range, improving the accuracy of abnormal area detection. A method of selecting feature points from the residual components based on corner response values ​​and calculating cross-correlation coefficients for matching addresses the issues of uneven illumination and large viewing angle variations within the milk duct inner wall, enhancing the robustness of feature tracking. Comparing the arc length and curvature changes of the feature point trajectories with the tissue motion period effectively identifies abnormal areas that do not conform to the physiological cycle, reducing the false positive detection rate. By calculating the similarity of motion features and analyzing connected domains, abnormal displacement areas are integrated to form a complete tissue abnormality map, which enables the identification of lesion areas to be upgraded from single-point judgment to overall regional analysis. This significantly improves the detection rate and positioning accuracy of lesions on the inner wall of the milk duct, providing reliable technical support for early clinical diagnosis.

[0060] Figure 2 Schematic diagram of feature point motion trajectory analysis and anomaly detection according to an embodiment of the present invention. Figure 2As shown in the figure, the process of feature point motion trajectory analysis and anomaly detection is intuitively demonstrated. The figure clearly marks four feature points, located at coordinates (7.5, 11.0), (15.0, 20.0), (25.0, 15.0) and (31.0, 25.0), as well as their motion trajectories. The blue line represents the normal trajectory, the red line represents the abnormal trajectory, and the red circle marks the abnormal displacement area. By calculating the tangent vector and normal vector of the trajectory and analyzing the arc length change and curvature change, this technical solution successfully identified multiple abnormal displacement points, such as the feature point (7.5, 11.0) with an abnormal offset of 3.0 in its third segment of the trajectory, and the feature point (15.0, 20.0) with an abnormal offset of 5.0 in its third segment of the trajectory. These abnormal points all exceed the fluctuation range of the normal motion cycle of the tissue. By marking and connecting these abnormal points, a clear abnormal displacement area is formed. This figure demonstrates the high precision and high sensitivity of this technical solution in trajectory analysis, which can accurately capture tiny abnormal motion patterns.

[0061] In an optional embodiment, the tissue abnormality map is divided into blocks, the local wavelet coefficients in each block area are calculated, the local wavelet coefficients are grouped and statistically analyzed to obtain internal structural features, the coefficient differences between adjacent blocks are calculated to obtain structural change features, and the internal structural features and structural change features are combined to form a regional descriptor, including: Perform recursive quadtree segmentation on the tissue abnormality map, calculate the information entropy of the pixel grayscale distribution in the block area, and determine whether to continue subdividing the current block based on the comparison result of the information entropy and the segmentation threshold until all blocks meet the stopping condition, thereby obtaining an adaptive block map; Perform wavelet decomposition on each block area in the adaptive block atlas to obtain detail coefficient matrices in the horizontal, vertical and diagonal directions; calculate the grayscale gradient distribution of the block area, generate directional weight coefficients, calculate statistics on the detail coefficient matrix and weight it using the directional weight coefficients to construct the internal structure feature vector of the block; Determine block adjacency based on the adaptive block map, perform scale normalization on adjacent block pairs, calculate the difference of local wavelet coefficients, construct a structural change matrix based on the difference, calculate the propagation characteristics of the structural change matrix at the region boundary, and perform cumulative statistics on the propagation characteristics to generate structural change characteristics; The internal structure features and the structural change features are aligned by dimension, a covariance matrix is ​​calculated for the aligned features, feature orthogonal decomposition is performed based on the covariance matrix, feature components are sorted and screened according to the decomposition results, and the screened feature components are combined to construct a region descriptor.

[0062] For example, the tissue abnormality map is first subjected to recursive quadtree segmentation, and the structural differences of the intraductal wall tissue are fully reflected by adaptive block division. In actual application, the initial tissue abnormality map can be set as a root node, and the entire image region can be used as an initial block. In the recursive segmentation process, each block is equally divided into four sub-blocks to form a quadtree structure. For each block region, the information entropy of the internal pixel gray scale distribution is calculated, and the information entropy reflects the complexity of the pixel value distribution in the region. When calculating the information entropy, the frequency of each gray scale value in the region is first counted, and then the entropy value is calculated according to the frequency. For example, for a block with a size of 64x64 pixels, if the gray scale distribution is relatively uniform, the calculated information entropy can be close to the maximum value 8 (for an 8-bit gray scale image); if the gray scale changes in the region are small, the information entropy can be only 2 to 3.

[0063] Whether to continue to subdivide the current block is determined according to the comparison result of the information entropy and the preset segmentation threshold. If the information entropy is greater than the segmentation threshold, it indicates that the pixel distribution in the region is complex and needs to be further subdivided; if the information entropy is less than or equal to the segmentation threshold, it indicates that the pixel distribution in the region is relatively uniform and the segmentation can be stopped. In the processing of the intraductal wall image, the segmentation threshold can be set to 4.5, which is verified by multiple experiments and can capture sufficient detail changes while maintaining the integrity of the region. In addition to the information entropy, a minimum block size limit can also be set, for example, 8x8 pixels, to prevent over-segmentation. The segmentation process continues until all blocks meet the stopping condition, and finally an adaptive block map is obtained. In actual application, the normal tissue region of the intraductal wall usually forms a larger block, while the lesion region is often divided into multiple small blocks due to high texture complexity.

[0064] Each block region in the adaptive block map is subjected to wavelet decomposition to obtain detail coefficient matrices in horizontal, vertical and diagonal directions. Wavelet decomposition can select Haar wavelet or Daubechies wavelet as the basis function to perform two-dimensional discrete wavelet transform on each block. In actual application, multi-level decomposition can be performed, and usually 2 to 3 levels are selected to capture texture features of different scales. After wavelet decomposition, each block region will generate a set of detail coefficient matrices, including horizontal, vertical and diagonal directions. These coefficient matrices reflect the texture characteristics of the intraductal wall tissue in different directions, and are of great significance for identifying the lesion tissue of the intraductal wall.

[0065] The gray scale gradient distribution of the sub-block region is calculated to generate the direction weight coefficient. The gray scale gradient reflects the direction and intensity of the pixel value change, which can be obtained by calculating the difference between adjacent pixels. For each sub-block region, the gradients in the horizontal and vertical directions are calculated, and then the gradient amplitude and direction are obtained. The distribution of the gradient direction is counted, and the histogram method can be used to divide the gradient direction into several intervals, for example, 360 degrees are divided into 8 intervals, and the cumulative amplitude of the gradient in each interval is calculated. According to the proportion of the gradient amplitude in each direction to the total amplitude, the direction weight coefficient is generated. For example, if the gradient amplitude in the horizontal direction accounts for 40%, the vertical direction accounts for 35%, and the diagonal direction accounts for 25%, the corresponding direction weight coefficients can be set to 0.4, 0.35, and 0.25.

[0066] The statistical quantity of the detail coefficient matrix is calculated and weighted by the direction weight coefficient to construct the internal structure feature vector of the sub-block. The statistical quantity can include mean, variance, skewness, kurtosis, etc., which describe the characteristics of the coefficient distribution from different angles. For example, the variance reflects the dispersion degree of the coefficient distribution, the skewness reflects the asymmetry of the distribution, and the kurtosis reflects the sharpness of the distribution. The statistical quantity of the detail coefficient in each direction is calculated, and then combined by weighting according to the direction weight coefficient to form the internal structure feature vector. In the analysis of the intraductal wall image, the normal tissue and the lesion tissue usually show obvious differences in these statistical quantities. For example, the wavelet coefficient variance of the intraductal papilloma region is usually larger than that of the normal tissue, reflecting the irregularity of the texture in the lesion region.

[0067] The adjacency relationship of the sub-block is determined based on the adaptive sub-block atlas, and the scale normalization processing is performed on the adjacent sub-block pair. The adjacency relationship can be determined by analyzing the spatial position of the sub-block, and if two sub-blocks share a boundary, they are considered to be adjacent. Since the adaptive segmentation results in different sizes of adjacent sub-blocks, scale normalization processing is needed to make the statistical features comparable. The normalization method can use bilinear interpolation or nearest neighbor interpolation to adjust the different size sub-blocks to a uniform size, for example, to 16x16 pixels.

[0068] The difference degree of the local wavelet coefficient is calculated, and the structure change matrix is constructed according to the difference degree. For each pair of adjacent sub-blocks, the difference in the wavelet coefficient in each direction is calculated. The difference degree can be measured by Euclidean distance, Manhattan distance, or correlation coefficient, etc. For example, the difference between the mean values of the wavelet coefficients of the two sub-blocks in the horizontal, vertical, and diagonal directions can be calculated to form a difference vector. The difference vectors of all adjacent sub-block pairs are organized into a structure change matrix, which reflects the spatial change trend of the intraductal wall tissue structure. In the intraductal wall image, the difference degree between the adjacent sub-blocks in the normal tissue region is usually small, while the difference degree at the boundary between the lesion region and the surrounding normal tissue is usually large.

[0069] The propagation feature of the structure change matrix at the region boundary is accumulated to generate a structure change feature. The propagation feature reflects how the structure difference spreads from one region to the surrounding regions. By regarding the blocks as nodes and the adjacency relationship as edges, the spread of the difference degree from one block to other blocks can be calculated by a graph theory method. For example, the block pair with the largest difference degree can be selected as the starting point to analyze how the difference degree changes with the increase of the distance between the blocks. Such analysis helps to identify the boundary and internal structure characteristics of the lesion region. The cumulative statistics of the propagation feature can obtain a feature vector reflecting the overall structure change. In the identification of the lesion region of the intraductal wall, this feature helps to distinguish the lesion boundary and the internal region, and improves the positioning accuracy.

[0070] The internal structure feature and the structure change feature are aligned in dimension, and a covariance matrix is calculated for the aligned features. Since the internal structure feature and the structure change feature can have different dimensions, dimension alignment processing is needed. The features of the two types can be made compatible in dimension by feature dimension reduction or padding, etc. After alignment, the two types of features are combined into a unified feature vector, and the covariance matrix of the vector is calculated. The covariance matrix reflects the correlation between the feature components, which helps to eliminate redundant information.

[0071] Based on the covariance matrix, feature orthogonalization decomposition is performed, and the feature components are sorted and selected according to the decomposition result. Orthogonalization decomposition can use principal component analysis method to convert the original features into a set of orthogonal principal components. Each principal component corresponds to a feature value, and the feature value size reflects the amount of information contained in the principal component. The principal components are sorted according to the feature value size, and the principal components with larger contribution rates are retained and the principal components with smaller contribution rates are discarded. In practical applications, a set of principal components with a cumulative contribution rate of 95% can be selected. This selection method not only retains the main information, but also reduces the feature dimension and improves the calculation efficiency.

[0072] The selected feature components are combined to construct a region descriptor. The region descriptor is a compact representation of the features of the abnormal region of the intraductal wall, which contains the key information of the internal structure and the structure change. In practical applications, the region descriptor can be used for lesion type identification and severity assessment. For example, different types of lesions such as intraductal papilloma and intraductal carcinoma usually exhibit different feature patterns in the region descriptor. By comparing with the pre-established lesion feature library, automatic identification and classification of the intraductal wall lesion can be realized.

[0073] In the embodiment, the recursive quadtree partitioning dynamically determines the block size according to the information entropy, so that the complex region is finely partitioned, the simple region remains integrity, the wavelet decomposition captures the multi-scale texture features, and the weighted processing of the directional weight coefficient enhances the sensitivity to the changes of the tissue structure in different directions of the inner wall of the breast duct. The calculation of the local wavelet coefficient difference degree between the blocks realizes the accurate expression of the boundary features of the tissue, so that the boundary features of the lesions and normal tissues are effectively extracted. The propagation feature analysis of the structure change matrix can accurately describe the transition characteristics of the internal structure of the lesion region and the surrounding normal tissue, and improves the recognition ability of the fuzzy region of the lesion boundary. The feature orthogonalization decomposition eliminates the redundant information and reduces the calculation complexity, the constructed region descriptor has high discrimination ability, so that the system can realize the accurate recognition and positioning of different types of lesions of the inner wall of the breast duct, and provides reliable technical support for the early clinical diagnosis.

[0074] In an optional implementation, the density clustering method is used to analyze the region descriptor to determine the lesion core region, including: The region descriptor is normalized to a preset numerical interval through the maximum-minimum value mapping to obtain a normalized region descriptor, a search radius is determined according to the feature distribution of the normalized region descriptor, a neighbor set of each normalized region descriptor is obtained within the search radius, and a distance matrix is calculated based on the neighbor set; A bandwidth parameter of an adaptive kernel function is determined according to the median of the distance matrix, the local density value of the normalized region descriptor is calculated using the adaptive kernel function, a target region descriptor is determined according to the local density value, the minimum distance value of each normalized region descriptor to the target region descriptor is calculated as a relative distance, and a density-distance distribution map is constructed according to the local density value and the relative distance; Statistical features of the density-distance distribution map are extracted, adaptive density and distance thresholds are set based on the statistical features, the normalized region descriptor is screened to obtain an initial clustering center; The density connectivity between the initial clustering centers is calculated, the initial clustering centers that are spatially adjacent are merged based on the density connectivity, and the clustering center with the largest local density value and the strongest density connectivity after the merging is determined as the position of the lesion core region.

[0075] Exemplarily, the region descriptor is first normalized to a preset numerical interval by maximum-minimum value mapping to obtain a normalized region descriptor. The region descriptor contains the structural features of the abnormal region of the intraductal wall tissue, and the feature value range of different dimensions can have a large difference, which needs to be normalized to eliminate the dimension influence. The maximum-minimum value mapping is a commonly used normalization method, which maps the original feature value to a preset numerical interval. In actual application, [0, 1] or [-1, 1] can be selected as the preset numerical interval. In the normalization process, the maximum value and the minimum value of each dimension feature of the region descriptor are calculated, and then the current feature value is mapped to the preset interval according to the relative position of the current feature value in the maximum-minimum value range. For example, if the minimum value of a certain dimension feature is 10, the maximum value is 50, and the current value is 30, and the preset interval is [0, 1], then the normalized value is (30-10) / (50-10)=0.5. The region descriptor extracted from the intraductal wall image is normalized to obtain a normalized region descriptor in a unified range, which is convenient for subsequent analysis.

[0076] The search radius is determined according to the feature distribution of the normalized region descriptor. The search radius is a key parameter of density clustering, which affects the calculation accuracy of local density. An excessively large search radius will cause excessive smoothing and loss of local features, and an excessively small search radius will cause unstable density estimation. In the recognition of intraductal wall lesions, the search radius can be dynamically determined according to the distribution characteristics of the normalized region descriptor. An effective method is to calculate the statistical distribution of the distance between samples in the feature space, and select a specific percentile of the distance distribution as the search radius. For example, the Euclidean distance between all sample pairs can be calculated, and the 15th percentile of the distance distribution is selected as the search radius. In actual application, for a data set containing 100 normalized region descriptors, the calculated search radius can be 0.15, indicating that about 15% of the sample pairs in the feature space have a distance less than this value.

[0077] The neighbor set of each normalized region descriptor is obtained within the search radius, and the distance matrix is calculated based on the neighbor set. For each normalized region descriptor, the distance between it and other descriptors is calculated, and if the distance is less than the search radius, the corresponding descriptor is included in the neighbor set. The distance calculation can use Euclidean distance, Manhattan distance or Mahalanobis distance and other measurement methods. In the recognition of intraductal wall lesions, Euclidean distance is a commonly used choice, which reflects the straight-line distance between two points in the feature space. Based on the neighbor set, a distance matrix is constructed to record the distance relationship between all normalized region descriptors. For sample pairs outside the neighbor set, the distance can be set to infinity or a large enough value, indicating that there is no direct connection between them.

[0078] The bandwidth parameter of the adaptive kernel function is determined based on the median of the distance matrix. The kernel function is an important tool for density estimation, used to calculate the local density of sample points. Common kernel functions include the Gaussian kernel and the Epanechnikov kernel. The bandwidth parameter controls the shape of the kernel function and directly affects the smoothness of the density estimate. In the identification of lesions on the inner wall of milk ducts, an adaptive bandwidth strategy can be adopted to determine the bandwidth parameter based on the statistical properties of the distance matrix. Specifically, the median of the valid distance values ​​(distances less than the search radius) in the distance matrix can be calculated and used as the bandwidth parameter. This method can automatically adjust the bandwidth based on the sparsity of the data distribution, improving the accuracy of density estimation. In practical applications, if the median of the distance matrix is ​​0.08, the bandwidth parameter can be set to 0.08.

[0079] Adaptive kernel function is used to calculate the local density value of normalized region descriptor. The local density value reflects the degree of clustering of sample points in the feature space and is the basis of density clustering. For each normalized region descriptor, the contribution of its neighboring points to it is calculated, and all contributions are accumulated to obtain the local density value. Specifically, the Gaussian kernel function can be used. For two sample points with a distance d, the contribution value is exp(-d 2 / h 2 ), where h is the bandwidth parameter. In the identification of milk duct lesions, descriptors in normal tissue regions are typically dispersed and have low local density values. However, descriptors in lesion regions, due to their high feature similarity, are often densely distributed and have high local density values.

[0080] Determine the target region descriptor based on the local density value. Target region descriptors are those with high local density values, which may correspond to the characteristic center of the lesion area. Based on the distribution characteristics of the local density values, descriptors with density values ​​greater than a certain threshold can be selected as target region descriptors. For example, descriptors with local density values ​​ranking in the top 20% or descriptors with density values ​​greater than the mean density plus one standard deviation can be selected. In the identification of milk duct lining lesions, target region descriptors often correspond to typical features of the lesion tissue and are important for determining the core area of ​​the lesion.

[0081] The minimum distance between each normalized region descriptor and the target region descriptor is calculated as the relative distance. Relative distance reflects the proximity of a sample point to a high-density region and is an important indicator for determining whether a sample point is a cluster center. For each normalized region descriptor, its distance to all target region descriptors is calculated, and the minimum value is taken as the relative distance. In the identification of milk duct lesions, descriptors for the lesion edge typically have medium local density values ​​and small relative distances, while descriptors for the lesion core often have high local density values ​​and large relative distances.

[0082] A density-distance distribution map is constructed based on local density values ​​and relative distances. A density-distance distribution map is a two-dimensional scatter plot with local density values ​​on the horizontal axis and relative distances on the vertical axis. Each normalized region descriptor corresponds to a point in the map. This visualization method can intuitively demonstrate the clustering characteristics of sample points. Ideally, cluster centers appear as points with high density and high distance in the map, clearly distinguishing them from other sample points. In the identification of milk duct lesions, density-distance distribution maps can help identify the characteristic centers of the lesion core.

[0083] Extract statistical features from the density-distance distribution map and use them to set adaptive density and distance thresholds. Statistical features can include the mean, variance, and skewness of the density distribution; the mean, variance, and skewness of the distance distribution; and the correlation coefficient between density and distance. These features reflect the overall distribution characteristics of the dataset and help determine appropriate thresholds. In practice, the density threshold can be set as the density mean plus the density standard deviation multiplied by a coefficient, and the distance threshold as the distance mean plus the distance standard deviation multiplied by a coefficient. The coefficient value can be adjusted based on actual needs and is typically between 0.5 and 2. For example, if the density mean is 0.3, the standard deviation is 0.1, and the coefficient is 1.5, then the density threshold is 0.3 + 0.1 × 1.5 = 0.45.

[0084] Normalized region descriptors are screened to obtain initial cluster centers. The screening criteria are: descriptors with local density values ​​greater than a density threshold and relative distances greater than a distance threshold are selected as initial cluster centers. These descriptors have both high local density, indicating that samples in their region are clustered, and large relative distances, indicating that they are clearly distinguishable from other high-density regions. In the identification of lesions in the inner wall of milk ducts, initial cluster centers typically correspond to lesion features of different types or locations. For example, for an image of the inner wall of a milk duct containing an intraductal papilloma, multiple initial cluster centers may be identified, corresponding to different locations or growth stages of the lesion.

[0085] Calculate the density connectivity between the initial cluster centers, and merge the spatially adjacent initial cluster centers based on the density connectivity. Density connectivity refers to the density continuity between two cluster centers, reflecting whether they belong to the same cluster. For two initial cluster centers, calculate the minimum density value of all points on the path between them and use it as the density connectivity. If the density connectivity is greater than the set threshold, it indicates that there is a density-continuous path between the two centers, which may belong to the same lesion area and should be merged. The density connectivity threshold can be set to a certain proportion of the mean density of all samples, such as 80% of the density mean. Merging spatially adjacent initial cluster centers can obtain a more compact cluster structure. In the identification of lesions on the inner wall of the milk duct, the merging process helps to eliminate the situation where the same lesion area is identified multiple times, thereby improving the accuracy of lesion localization.

[0086] The cluster center with the largest local density value and the strongest density connectivity after merging is determined as the location of the core area of ​​the lesion. Density connectivity refers to the degree of density connection between the cluster center and other sample points, which can be measured by calculating the number of density connectivity paths from the cluster center to other sample points. The stronger the density connectivity, the more concentrated the characteristic distribution of the area where the cluster center is located, and the more likely it is to be the core area of ​​the lesion. In practical applications, the local density value and density connectivity can be considered comprehensively, for example, by calculating the weighted sum of the two, and the weights can be adjusted according to actual needs. The cluster center with the largest weighted sum is selected as the location of the core area of ​​the lesion. In the identification of lesions on the inner wall of the milk duct, the core area of ​​the lesion usually corresponds to the most typical and concentrated part of the lesion tissue, and is the key to accurate positioning and diagnosis.

[0087] In this embodiment, the maximum-minimum mapping normalization technique eliminates dimensional differences between dimensions of the region descriptor, improving the accuracy of distance calculation in the feature space. An adaptive search radius determination strategy dynamically adjusts parameters based on the feature distribution, enhancing the algorithm's adaptability to lesions with varying density distributions. The adaptive kernel function bandwidth parameter, determined based on the median of the distance matrix, enables more accurate density estimation, effectively addressing complex conditions such as uneven illumination and tissue deformation in ductal wall images. A joint analysis of local density values ​​and relative distances successfully distinguishes lesion core and edge regions, reducing the false positive rate. Statistical analysis of the density-distance distribution map enables adaptive adjustment of threshold parameters, ensuring algorithm stability across different patients and lesion types. Density connectivity calculation and the merging of spatially adjacent cluster centers effectively address the problem of multiple identifications of the same lesion, improving localization accuracy. This method demonstrates excellent recognition of irregularly shaped and variably sized ductal wall lesions, providing reliable technical support for the precise localization of early, subtle lesions and significantly improving the sensitivity and specificity of ductal wall lesion diagnosis.

[0088] In an optional embodiment, a multi-scale radial basis function is constructed with the lesion core area as the center, the response distribution of the regional descriptor under the multi-scale radial basis function is calculated, the spatial gradient contour of the response distribution is extracted, the spatial range of the lesion area is determined based on the contour, and the evolution of the contour is tracked by optical flow field analysis to determine the direction of lesion expansion, including: Taking the center position of the lesion core area as the reference point, calculating the reference scale value according to the spatial range of the lesion core area, constructing an increasing scale sequence according to a preset ratio based on the reference scale value, and constructing a multi-scale radial basis function using the increasing scale sequence; Calculating the response distribution of the region descriptor under the multi-scale radial basis function, generating an anisotropic weight matrix based on the structural feature gradient of the region descriptor, and multiplying the anisotropic weight matrix by the response distribution to obtain a modified response distribution; extracting spatial gradient isograms of the corrected response distribution at each scale, calculating a stability index by the proportion of spatial overlap of the isograms at adjacent scales, extracting a change index by the response gradient of the neighborhood of the isograms, calculating a consistency index by the difference of morphological features of the isograms and the lesion core region; selecting the spatial gradient isograms based on the stability index, the change index and the consistency index, and determining the maximum closed isograms after the selection as the spatial range of the lesion region; calculating the optical flow field vector on the boundary of the spatial range of the lesion region, and determining the expansion direction of the lesion region by the eigenvector direction corresponding to the maximum eigenvalue obtained by the principal component analysis of the optical flow field vector.

[0089] Exemplarily, the center position of the lesion core region is taken as the reference point, and the reference scale value is calculated according to the spatial range of the lesion core region. In actual application, the lesion core region is a feature aggregation region identified from the region descriptors by the density clustering method. The center position of the lesion core region can be obtained by calculating the average coordinates of all points in the region, and taken as the reference point of the radial basis function. The spatial range of the lesion core region can be determined by calculating the expansion distance of the region in each direction. The reference scale value is an important parameter for constructing the multi-scale radial basis function, and directly affects the shape and coverage of the function. In the identification of intraductal wall lesions, the average radius of the lesion core region can be taken as the reference scale value. For example, if the lesion core region presents an approximate circular distribution in space, and the average radius is 5 pixels, the reference scale value can be set to 5.

[0090] Based on the reference scale value, an incremental scale sequence is constructed according to a preset ratio, and the multi-scale radial basis function is constructed by using the incremental scale sequence. The incremental scale sequence refers to a series of scale values from small to large, which is used to construct radial basis functions of different sizes. The incremental ratio can be set in the form of geometric progression, for example, increasing by 50% each time, that is, the scale sequence is the reference scale value multiplied by 1, 1.5, 2.25, 3.375, etc. In actual application, 5 to 8 incremental scales can be selected to cover from the lesion core to the possible lesion edge region. For each scale value, a corresponding radial basis function is constructed. The radial basis function is a function with distance as the independent variable, and the commonly used forms include Gaussian function, multi-quadratic function, etc. In the identification of intraductal wall lesions, the Gaussian type radial basis function can be used, and the function value decreases with the increase of the distance from the point to the reference point. Specifically, for the distance d and the scale s, the function value can be expressed as exp(-d 2 / s 2 The set of radial basis functions constructed by different scales forms the multi-scale radial basis function, which can capture the feature distribution of the lesion region from different spatial scales.

[0091] The response distribution of the region descriptor under a multi-scale radial basis function is calculated. The region descriptor incorporates the structural features of abnormal regions of the milk duct lining tissue. Combining this with the multi-scale radial basis function reveals the distribution of these features at different spatial scales. Specifically, for each spatial location, the region descriptor is multiplied by the radial basis function to obtain a response value. The response values ​​at all locations are organized into a response distribution map, which reflects the spatial distribution intensity of the features. In milk duct lining lesion identification, the response value of the lesion region is typically higher than that of the surrounding normal tissue, appearing as a highlighted area in the response distribution map.

[0092] An anisotropic weight matrix is ​​generated based on the structural feature gradient of the regional descriptor, and the anisotropic weight matrix is ​​multiplied by the response distribution to obtain a modified response distribution. The anisotropic weight matrix is ​​used to adjust the weights of the response distribution in different directions so that the modified response distribution better reflects the directional characteristics of the tissue structure. The structural feature gradient reflects the changing trend of the feature in space and can be obtained by calculating the partial derivatives of the regional descriptor in the horizontal and vertical directions. Based on the gradient information, a second-order tensor is constructed, which contains the components of the gradient in each direction. The tensor is subjected to eigendecomposition to obtain the main direction and degree of anisotropy, and the anisotropic weight matrix is ​​generated based on this. In the identification of lesions on the inner wall of the milk duct, the diseased tissue often exhibits obvious structural features in certain directions. By adjusting the anisotropic weight matrix, the response in these directions can be enhanced, improving the recognition accuracy of the lesion boundary.

[0093] The spatial gradient contour lines of the modified response distribution are extracted at each scale. Spatial gradient contour lines refer to curves with equal response distribution gradient amplitudes, which usually correspond to the boundary areas of the response distribution. Edge detection and contour extraction methods can be used to extract spatial gradient contour lines. In the identification of lesions on the inner wall of the milk duct, the gradient amplitude map of the modified response distribution can be calculated first, and then threshold segmentation is applied to the gradient amplitude map to extract the set of points with gradient amplitudes equal to specific values. These points are connected into closed curves, which are the spatial gradient contour lines. Generally, contour lines corresponding to multiple different gradient values ​​can be extracted to form a set of contour lines. For example, the 25%, 50%, and 75% quantiles of the gradient amplitude can be selected as thresholds to extract three groups of contour lines.

[0094] The stability index is obtained by calculating the spatial overlap ratio of the contours at adjacent scales. The stability index reflects the stability degree of the contours at scale change, and is an important basis for evaluating the quality of the contours. For the contours at two adjacent scales, the spatial overlap area ratio of the contours is calculated. The larger the ratio is, the more stable the contours are. In the identification of the lesions in the inner wall of the breast duct, the stable contours usually correspond to the real tissue boundary, while the unstable contours may be caused by noise or artifacts. For example, if the overlap area of the contours at two adjacent scales is 80 pixels, and the union area is 100 pixels, the stability index is 0.8, indicating that the contours are relatively stable at scale change.

[0095] The change index is obtained by extracting the response gradient of the neighborhood of the contours. The change index reflects the change degree of the response distribution near the contours, and is an important basis for judging whether the contours are located in the boundary region. For each contour, the response gradient in the neighborhood of the contour is extracted, and the average or median of the gradient amplitude is calculated as the change index. In the identification of the lesions in the inner wall of the breast duct, the contours located at the boundary of the lesions usually have a larger change index, because the response value in the boundary region changes significantly; while the contours located in the uniform region have a relatively small change index. For example, if the average amplitude of the response gradient in the neighborhood of a contour is 0.4, which is higher than the overall average value of 0.2, it indicates that the contour may be located in the boundary region.

[0096] The consistency index is obtained by calculating the difference of the morphological features between the contours and the core region of the lesions. The consistency index reflects whether the shape of the contours is consistent with the morphological features of the core region of the lesions, and is an important basis for screening suitable contours. The morphological features can include circularity, long axis to short axis ratio, direction consistency, etc. For each contour, the difference between the morphological features of the contour and the morphological features of the core region of the lesions is calculated. The smaller the difference is, the higher the consistency index is. In the identification of the lesions in the inner wall of the breast duct, the real lesion boundary usually has similar morphological features with the core region, such as both showing circular or elliptical shape. For example, if the circularity of the core region of the lesions is 0.85, and the circularity of a contour is 0.82, the morphological difference between them is small, and the consistency index is high.

[0097] Based on the stability index, change index and consistency index, the spatial gradient contour line is optimized, and the largest closed contour line after optimization is determined as the spatial range of the lesion area. The optimization process can adopt a weighted scoring method to comprehensively consider the values ​​of the three indicators. For example, the three indicators can be normalized to the interval [0, 1], and the weights are set as w1, w2 and w3, and the weighted sum is calculated as the comprehensive score of the contour line. The contour line with the highest score is selected as the optimal contour line. If there are multiple contour lines with similar scores, the largest closed contour line can be selected as the spatial range of the lesion area. In the identification of lesions on the inner wall of the milk duct, the spatial range of the lesion area is usually manifested as a closed boundary that includes the entire area of ​​the lesion tissue.

[0098] The optical flow field vector is calculated for the spatial range boundary of the lesion area. The optical flow field vector reflects the motion trend of the pixel points in the image sequence and can be used to analyze the expansion direction of the lesion area. In the continuous image sequence of the inner wall of the milk duct, the optical flow field vector is extracted for the determined lesion area boundary. The optical flow field calculation can use the Lucas-Kanade method or the Horn-Schunck method to solve the motion vector based on the image grayscale gradient and time derivative. In practical applications, multiple points can be uniformly sampled on the lesion boundary, and the optical flow vectors of these points can be calculated to form a boundary optical flow field. For example, in a sequence containing 100 frames of images, the optical flow field can be calculated every 5 frames to obtain the motion trend of the lesion boundary at different times.

[0099] Principal component analysis is performed on the optical flow field vectors, and the direction of the eigenvector corresponding to the maximum eigenvalue obtained from the principal component analysis is determined as the expansion direction of the lesion area. Principal component analysis is a dimensionality reduction method that can be used to extract the main direction of change in data. Principal component analysis is performed on the data matrix composed of boundary optical flow field vectors, and the eigenvalues ​​and eigenvectors of its covariance matrix are calculated. The eigenvector represents the main direction of data change, and the eigenvalue represents the degree of change in that direction. The eigenvector corresponding to the maximum eigenvalue is selected, and its direction is the main expansion direction of the lesion area. In the identification of lesions on the inner wall of the milk duct, the identification of the lesion expansion direction is important for assessing the development trend of the lesion and formulating treatment plans. For example, if the direction of the maximum eigenvector obtained by principal component analysis is 45 degrees, it indicates that the lesion area is mainly expanding in the upper right direction.

[0100] In the embodiment, the lesion core region is used as a reference point to construct an adaptive scale sequence, which overcomes the limitations of the traditional fixed scale method in adapting to different sizes of lesions. The introduction of the anisotropic weight matrix effectively enhances the structural characteristics of the lesion tissue in a specific direction, improving the accuracy of lesion boundary recognition. The multi-scale spatial gradient contour extraction technology can capture lesion boundary information from different spatial scales, and the comprehensive evaluation mechanism of stability index, change index and consistency index effectively selects the optimal contour, reducing the interference of noise and artifacts. The extended direction determination method combining optical flow field vector analysis and principal component analysis can accurately capture the dynamic evolution trend of the lesion area, providing a basis for predicting the development of the lesion for clinicians.

[0101] Figure 3 For the lesion analysis effect comparison chart of the embodiment of the application, as shown in Figure 3 , the chart compares and displays the performance evaluation results of three different lesion analysis methods. The multi-scale radial basis function method generally performs well in the five evaluation indexes, especially in the boundary accuracy and clinical applicability, reaching high scores of 92.3% and 91.2%, respectively, which is significantly better than the traditional method. The traditional gradient method performs moderately, and each index is uniformly distributed between 72% and 85%. Although the region growing method performs well in computational efficiency (90.6%), it is relatively weak in other indexes. The comparative analysis confirms that the lesion analysis method based on multi-scale radial basis function has significant advantages in accurately capturing the lesion boundary and adapting to the clinical environment, especially through the comprehensive evaluation of stability index, change index and consistency index, which can more accurately determine the spatial range and expansion direction of the lesion area In a second aspect of the embodiment of the application, a real-time identification and positioning system for intraductal wall lesions based on optical fiber imaging is provided, which comprises: A first unit for acquiring an intraductal wall image sequence and dividing it into multiple subsequences according to a preset time interval; A second unit for extracting the displacement vector field between adjacent images in each subsequence, decomposing the displacement vector field into a main frequency component and a residual component in the frequency domain, obtaining the tissue motion period through phase analysis of the main frequency component, screening abnormal displacement regions that do not conform to the tissue motion period from the residual component, merging the abnormal displacement regions in multiple subsequences, and establishing a tissue abnormality atlas; A third unit for blocking the tissue abnormality atlas, calculating the local wavelet coefficients in each block region, grouping and statistically analyzing the local wavelet coefficients to obtain internal structure features, calculating the coefficient difference between adjacent blocks to obtain structure change features, and combining the internal structure features and structure change features to form a region descriptor; The fourth unit is configured to analyze the regional descriptor by using a density clustering method to determine a lesion core region, to construct a multi-scale radial basis function centered on the lesion core region, to calculate a response distribution of the regional descriptor under the multi-scale radial basis function, to extract a spatial gradient contour line of the response distribution, to determine a spatial range of the lesion region based on the contour line, and to track evolution of the contour line by using an optical flow field analysis to determine a lesion expansion direction.

[0102] In a third aspect, an electronic device is provided, comprising: a processor; a memory for storing processor-executable instructions; wherein the processor is configured to invoke the instructions stored in the memory to perform the method described above.

[0103] In a fourth aspect, a computer-readable storage medium is provided, which stores computer program instructions, and the computer program instructions are executed by a processor to implement the method described above.

[0104] The present application can be a method, apparatus, system and / or computer program product. The computer program product can include a computer-readable storage medium having stored thereon computer-readable program instructions that, when executed by a computer, cause the computer to carry out various aspects of the present application.

[0105] Finally, it should be noted that: the above embodiments are only used to illustrate the technical solutions of the present application, and not to limit them; although the present application has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand: it can still modify the technical solutions recorded in the foregoing embodiments, or make equivalent replacement to part or all of the technical features; and these modifications or replacements do not make the essence of the corresponding technical solutions deviate from the scope of the technical solutions of the embodiments of the present application.

Claims

1. A method for real-time identification and positioning of milk duct inner wall lesions based on optical fiber imaging, characterized in that: include: Acquire a sequence of images of the inner wall of the milk duct and divide the sequence into a plurality of subsequences according to a preset time interval; The displacement vector field between adjacent images in each subsequence is extracted and decomposed into a primary frequency component and a residual component in the frequency domain. The tissue motion cycle is obtained by phase analysis of the primary frequency component. Abnormal displacement regions that do not conform to the tissue motion cycle are screened out from the residual component, and the abnormal displacement regions in multiple subsequences are merged to establish a tissue abnormality map. The tissue abnormality map is divided into blocks, and the local wavelet coefficients in each block area are calculated. The local wavelet coefficients are grouped and statistically analyzed to obtain internal structural features. The coefficient differences between adjacent blocks are calculated to obtain structural change features. The internal structural features and structural change features are combined to form a regional descriptor. The density clustering method is used to analyze the regional descriptor to determine the core area of ​​the lesion. A multi-scale radial basis function is constructed with the core area of ​​the lesion as the center. The response distribution of the regional descriptor under the multi-scale radial basis function is calculated, and the spatial gradient contour lines of the response distribution are extracted. The spatial range of the lesion area is determined based on the contour lines, and the evolution of the contour lines is tracked by optical flow field analysis to determine the direction of lesion expansion.

2. The method according to claim 1, characterized in that Acquiring a sequence of images of the inner wall of the milk duct and dividing the sequence into multiple subsequences according to preset time intervals includes: Acquiring a continuous image sequence of the inner wall of the milk duct, performing time domain analysis on the continuous image sequence to extract a periodic signal of the image sequence, setting a time segmentation benchmark according to the periodic signal, and dividing the continuous image sequence into initial subsequences using an integer multiple of the time segmentation benchmark as a subsequence duration; An overlapping interval is set between adjacent sequences of the initial subsequence, and the duration of the overlapping interval is a preset proportion of the subsequence duration; a weight value is calculated based on the temporal distance between the image and the boundary of the overlapping interval, and the images in the overlapping interval are weighted to obtain a transition sequence in the overlapping interval. The transition sequence replaces the image sequence in the original overlapping interval to generate multiple temporally continuous image subsequences.

3. The method according to claim 1, characterized in that Extract the displacement vector field between adjacent images in each subsequence and decompose the displacement vector field into the main frequency component and the residual component in the frequency domain, including: Obtaining the grayscale values ​​of adjacent image frames in a subsequence and constructing an image block centered on a pixel point; establishing grayscale constancy constraints and spatial consistency constraints within the image block to form a local grayscale change constraint matrix; substituting the local grayscale change constraint matrix into the optical flow equation to obtain a pixel-level displacement vector; and performing regional integration on the pixel-level displacement vector to obtain a displacement vector field; Decomposing the displacement vector field using wavelet basis functions to obtain displacement components of different scales; calculating the amplitude distribution of the displacement components of each scale, generating directional weight coefficients, fusing the displacement components according to the directional weight coefficients, and reconstructing the characteristics of the displacement vector field; A two-dimensional Fourier transform is performed on the displacement vector field characteristics to obtain a frequency domain representation, the energy distribution of each frequency component in the frequency domain representation is statistically analyzed, and the cumulative energy value of the frequency component is calculated; a frequency segmentation threshold is determined based on the cumulative energy value, the frequency domain representation is divided with the frequency segmentation threshold as the boundary and reconstructed through an inverse Fourier transform to obtain the main frequency component and the residual component.

4. The method according to claim 1, wherein The tissue motion cycle is obtained by phase analysis of the main frequency component; abnormal displacement areas that do not conform to the tissue motion cycle are screened out from the residual components, and abnormal displacement areas in multiple subsequences are merged to establish a tissue abnormality map, including: Performing a Hilbert transform on the main frequency component to extract the instantaneous phase, counting the temporal changes of the instantaneous phase, identifying extreme points in the phase sequence, calculating the time interval sequence between adjacent extreme points, and determining the tissue motion cycle and its fluctuation range based on the time interval sequence; Dividing the residual components into grid points, calculating the corner point response values ​​of the grid points, selecting the grid points whose corner point response values ​​are greater than a set response threshold as feature points, calculating the mutual correlation coefficients of the feature points between adjacent frames to perform feature point matching, and obtaining the position offset of the feature points based on the matching results; The position offsets are concatenated and filtered along the time dimension to obtain the motion trajectory of the feature points, the tangent vector and normal vector of the motion trajectory are calculated, and the arc length change and curvature change of the trajectory are obtained. The arc length change and curvature change are compared with the tissue motion cycle, and the trajectory points that exceed the cycle fluctuation range are marked. The area containing the marked points is extracted as the abnormal displacement area; The motion direction, velocity, and acceleration characteristics of the abnormal displacement area are calculated to construct motion features. The similarity between the abnormal displacement areas is calculated based on the motion features. The abnormal displacement areas whose similarity is greater than a preset similarity threshold and whose spatial positions are adjacent are analyzed for connected domains, and the connected domains are integrated into a tissue abnormality map.

5. The method according to claim 1, wherein The tissue abnormality map is divided into blocks, the local wavelet coefficients in each block area are calculated, the local wavelet coefficients are grouped and statistically analyzed to obtain internal structural features, the coefficient differences between adjacent blocks are calculated to obtain structural change features, and the internal structural features and structural change features are combined to form a regional descriptor including: Perform recursive quadtree segmentation on the tissue abnormality map, calculate the information entropy of the pixel grayscale distribution in the block area, and determine whether to continue subdividing the current block based on the comparison result of the information entropy and the segmentation threshold until all blocks meet the stopping condition, thereby obtaining an adaptive block map; Perform wavelet decomposition on each block area in the adaptive block atlas to obtain detail coefficient matrices in the horizontal, vertical and diagonal directions; calculate the grayscale gradient distribution of the block area, generate directional weight coefficients, calculate statistics on the detail coefficient matrix and weight it using the directional weight coefficients to construct the internal structure feature vector of the block; Determine block adjacency based on the adaptive block map, perform scale normalization on adjacent block pairs, calculate the difference of local wavelet coefficients, construct a structural change matrix based on the difference, calculate the propagation characteristics of the structural change matrix at the region boundary, and perform cumulative statistics on the propagation characteristics to generate structural change characteristics; The internal structure features and the structural change features are aligned by dimension, a covariance matrix is ​​calculated for the aligned features, feature orthogonal decomposition is performed based on the covariance matrix, feature components are sorted and screened according to the decomposition results, and the screened feature components are combined to construct a region descriptor.

6. The method according to claim 1, characterized in that The density clustering method is used to analyze the regional descriptors to determine the core areas of the lesions, including: The region descriptor is normalized to a preset numerical range through maximum and minimum value mapping to obtain a normalized region descriptor. The search radius is determined according to the feature distribution of the normalized region descriptor. The neighbor set of each normalized region descriptor is obtained within the search radius, and the distance matrix is ​​calculated based on the neighbor set. Determining a bandwidth parameter of an adaptive kernel function according to the median of the distance matrix, calculating a local density value of a normalized region descriptor using the adaptive kernel function, determining a target region descriptor according to the local density value, calculating a minimum distance value from each normalized region descriptor to the target region descriptor as a relative distance, and constructing a density-distance distribution map according to the local density value and the relative distance; Extracting statistical features of the density-distance distribution map, setting an adaptive density threshold and a distance threshold based on the statistical features, and screening normalized region descriptors to obtain initial cluster centers; The density connectivity between the initial cluster centers is calculated, and the spatially adjacent initial cluster centers are merged based on the density connectivity. The cluster center with the maximum local density value and the strongest density connectivity after the merger is determined as the location of the lesion core area.

7. The method according to claim 1, characterized in that A multi-scale radial basis function is constructed with the lesion core area as the center. The response distribution of the regional descriptor under the multi-scale radial basis function is calculated. The spatial gradient contour of the response distribution is extracted. The spatial range of the lesion area is determined based on the contour. The evolution of the contour is tracked by optical flow field analysis to determine the direction of lesion expansion. Taking the center position of the lesion core area as the reference point, calculating the reference scale value according to the spatial range of the lesion core area, constructing an increasing scale sequence according to a preset ratio based on the reference scale value, and constructing a multi-scale radial basis function using the increasing scale sequence; Calculating the response distribution of the region descriptor under the multi-scale radial basis function, generating an anisotropic weight matrix based on the structural feature gradient of the region descriptor, and multiplying the anisotropic weight matrix by the response distribution to obtain a modified response distribution; Extract the spatial gradient contour lines of the modified response distribution at each scale, calculate the spatial overlap ratio of the contour lines at adjacent scales to obtain a stability index, extract the response gradient of the contour line neighborhood to obtain a change index, and calculate the difference in morphological characteristics between the contour line and the lesion core area to obtain a consistency index; Optimizing the spatial gradient contour line based on the stability index, the change index, and the consistency index, and determining the optimal maximum closed contour line as the spatial range of the lesion area; An optical flow field vector is calculated for the spatial range boundary of the lesion area, a principal component analysis is performed on the optical flow field vector, and the direction of the eigenvector corresponding to the maximum eigenvalue obtained by the principal component analysis is determined as the expansion direction of the lesion area.

8. A real-time identification and positioning system for milk duct inner wall lesions based on optical fiber imaging, used to implement the method according to any one of claims 1 to 7, characterized in that: include: The first unit is used to obtain a sequence of images of the inner wall of the milk duct and divide the sequence into multiple subsequences according to preset time intervals; The second unit is used to extract the displacement vector field between adjacent images in each subsequence, decompose the displacement vector field into a primary frequency component and a residual component in the frequency domain, and obtain the tissue motion cycle through phase analysis of the primary frequency component. Abnormal displacement areas that do not conform to the tissue motion cycle are screened out from the residual component, and abnormal displacement areas in multiple subsequences are merged to establish a tissue abnormality map. The third unit is used to divide the tissue abnormality map into blocks, calculate the local wavelet coefficients in each block area, perform group statistics on the local wavelet coefficients to obtain internal structural features, calculate the coefficient differences between adjacent blocks to obtain structural change features, and combine the internal structural features and structural change features to form a regional descriptor; The fourth unit is used to analyze the regional descriptor using the density clustering method to determine the core area of ​​the lesion, construct a multi-scale radial basis function with the core area of ​​the lesion as the center, calculate the response distribution of the regional descriptor under the multi-scale radial basis function, extract the spatial gradient contour lines of the response distribution, determine the spatial range of the lesion area based on the contour lines, and track the evolution of the contour lines through optical flow field analysis to determine the direction of lesion expansion.

9. An electronic device, characterized in that: include: processor; a memory for storing processor-executable instructions; The processor is configured to call the instructions stored in the memory to execute the method according to any one of claims 1 to 7.

10. A computer-readable storage medium having computer program instructions stored thereon, characterized in that: When the computer program instructions are executed by a processor, the method according to any one of claims 1 to 7 is implemented.

Images