Real-time recognition and positioning method and system for intraductal wall lesions of breast based on optical fiber imaging

Abstract

The application provides a breast duct inner wall lesion real-time identification and positioning method and system based on optical fiber imaging, relates to the technical field of medical image processing, and comprises the following steps: acquiring a breast duct inner wall image sequence, extracting a displacement vector field and decomposing the displacement vector field into a main frequency and a residual component, analyzing the main frequency phase to obtain a tissue motion cycle, screening an abnormal displacement area from the residual component to establish a tissue abnormality atlas, constructing a region descriptor through local wavelet coefficient analysis, determining a lesion core area by applying density clustering, and finally determining a lesion area range and an extension direction by using a multi-scale radial basis function and contour line analysis. The application can realize real-time accurate identification and positioning of breast duct inner wall lesions.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of medical image processing, 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-stage small 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. By minimally invasive way, 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 early-stage breast disease diagnosis. Currently, optical fiber ductoscopy has been preliminarily applied in clinical practice, which can obtain image sequences of the breast duct wall to assist doctors in diagnosis.

[0003] The existing ductoscopy technology still has many shortcomings 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 is low in 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-stage small lesions and normal tissues is relatively subtle. The traditional recognition method based on morphological features is difficult to 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 early-stage diagnosis efficiency and accuracy of breast diseases. 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:

[0007] obtaining breast duct wall image sequences and dividing them into multiple subsequences according to a preset time interval;

[0008] extracting a displacement vector field between adjacent images in each sub-sequence, decomposing the displacement vector field into a dominant frequency component and a residual component in a frequency domain, obtaining a tissue motion cycle through phase analysis of the dominant frequency component; screening abnormal displacement regions not conforming to the tissue motion cycle from the residual component, merging the abnormal displacement regions in multiple sub-sequences, and establishing a tissue abnormality atlas;

[0009] blocking the tissue abnormality atlas, calculating local wavelet coefficients in each block region, grouping and counting the local wavelet coefficients to obtain internal structure features, calculating coefficient differences between adjacent blocks to obtain structure change features, and combining the internal structure features and the structure change features to form a region descriptor;

[0010] analyzing the region descriptor using a density clustering method to determine a lesion core region, constructing a multi-scale radial basis function centered on the lesion core region, calculating a response distribution of the region descriptor under the multi-scale radial basis function, extracting a spatial gradient contour line of the response distribution, determining a spatial range of the lesion region based on the contour line, and determining a lesion extension direction through light flow field analysis and tracking of the contour line evolution.

[0011] In an optional embodiment,

[0012] obtaining a sequence of images of an inner wall of a breast duct and dividing the sequence into multiple sub-sequences according to a preset time interval includes:

[0013] obtaining a sequence of continuous images of an inner wall of a breast duct, performing time domain analysis on the sequence of continuous images, extracting a periodic signal of the sequence of images, setting a time segmentation reference according to the periodic signal, taking an integer multiple of the time segmentation reference as a sub-sequence length, and dividing the sequence of continuous images into initial sub-sequences.

[0014] setting an overlapping interval between adjacent sequences of the initial sub-sequences, the overlapping interval having a length of a preset proportion of the sub-sequence length; calculating a weight value according to a time distance between an image and a boundary of the overlapping interval, performing weighted processing on the image in the overlapping interval to obtain a transition sequence in the overlapping interval, replacing the image sequence in the original overlapping interval with the transition sequence, and generating multiple image sub-sequences that are continuous in time sequence.

[0015] In an optional embodiment,

[0016] extracting a displacement vector field between adjacent images in each sub-sequence, decomposing the displacement vector field into a dominant frequency component and a residual component in a frequency domain, obtaining a tissue motion cycle through phase analysis of the dominant frequency component; screening abnormal displacement regions not conforming to the tissue motion cycle from the residual component, merging the abnormal displacement regions in multiple sub-sequences, and establishing a tissue abnormality atlas;

[0017] Obtain the grayscale values ​​of adjacent image frames in the subsequence, and construct image blocks centered on pixels; establish grayscale constantness constraints and spatial consistency constraints within the image blocks to form a local grayscale change constraint matrix; substitute the local grayscale change constraint matrix into the optical flow equation to solve for pixel-level displacement vectors; and perform regional integration on the pixel-level displacement vectors to obtain a displacement vector field.

[0018] The displacement vector field is decomposed using wavelet basis functions to obtain displacement components at different scales; the amplitude distribution of each scale displacement component is calculated to generate directional weighting coefficients; the displacement components are fused based on the directional weighting coefficients to reconstruct the characteristics of the displacement vector field.

[0019] 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 components is calculated. A frequency segmentation threshold is determined based on the cumulative energy value. The frequency domain representation is divided by the frequency segmentation threshold and reconstructed by inverse Fourier transform to obtain the dominant frequency component and the residual component.

[0020] In one alternative embodiment,

[0021] The tissue motion cycle is obtained through phase analysis of the dominant frequency components; abnormal displacement regions that do not conform to the tissue motion cycle are screened from the residual components, and abnormal displacement regions from multiple subsequences are merged to establish a tissue anomaly map, including:

[0022] The Hilbert transform is applied to the main frequency component to extract the instantaneous phase. The temporal changes of the instantaneous phase are statistically analyzed to identify the extreme points in the phase sequence. The time interval sequence between adjacent extreme points is calculated, and the tissue movement cycle and its fluctuation range are determined based on the time interval sequence.

[0023] Divide the residual components into grid points, calculate the corner response values ​​of the grid points, select grid points whose corner response values ​​are greater than a set response threshold as feature points, calculate the cross-correlation coefficient of the feature points between adjacent frames to perform feature point matching, and obtain the position offset of the feature points based on the matching results.

[0024] The position offset is concatenated along the time dimension and filtered to obtain the motion trajectory of the feature point. The tangent vector and normal vector of the motion trajectory are calculated to obtain the arc length change and curvature change of the trajectory. The arc length change and curvature change are compared with the tissue motion cycle to mark the trajectory points that exceed the cycle fluctuation range. The area containing the marked points is extracted as the abnormal displacement area.

[0025] 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 domain analysis, and the connected domains are integrated into an organization abnormality map.

[0026] In an alternative embodiment,

[0027] The organization abnormality map is divided into blocks, the local wavelet coefficients in each block region are calculated, the local wavelet coefficients are grouped and counted to obtain internal structure characteristics, the coefficient difference between adjacent blocks is calculated to obtain structure change characteristics, and the internal structure characteristics and the structure change characteristics are combined to form a region descriptor including:

[0028] The organization abnormality map is subjected to recursive quadtree partitioning, the information entropy of the pixel gray scale distribution in the block region is calculated, and it is determined whether to continue to subdivide the current block according to the comparison result of the information entropy and the partition threshold, until all blocks meet the stop condition, to obtain an adaptive block map;

[0029] Each block region in the adaptive block map is subjected to wavelet decomposition to obtain detail coefficient matrices in horizontal, vertical and diagonal directions; the gray scale gradient distribution of the block region is calculated to generate direction weight coefficients, and the statistics of the detail coefficient matrices are calculated and weighted by the direction weight coefficients to construct an internal structure characteristic vector of the block;

[0030] Based on the adaptive block map, the block adjacency relationship is determined, the adjacent block pairs are subjected to scale normalization processing, the difference of the local wavelet coefficients is calculated, a structure change matrix is constructed according to the difference, the propagation characteristics of the structure change matrix at the region boundary are calculated, the propagation characteristics are accumulated and counted to generate a structure change characteristic;

[0031] The internal structure characteristics and the structure change characteristics are aligned in dimension, a covariance matrix is calculated for the aligned characteristics, a feature orthogonalization decomposition is performed based on the covariance matrix, the feature components are sorted and filtered according to the decomposition result, and the filtered feature components are combined to construct a region descriptor.

[0032] In an alternative embodiment,

[0033] The density clustering method is used to analyze the region descriptor to determine a lesion core region including:

[0034] The region descriptor is normalized to a preset numerical interval by maximum and 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 near neighbor set of each normalized region descriptor is obtained within the search radius, and a distance matrix is calculated based on the near neighbor set.

[0035] A bandwidth parameter of an adaptive kernel function is determined according to a median of the distance matrix, a 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, a 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;

[0036] Statistical features of the density-distance distribution map are extracted, an adaptive density threshold and a distance threshold are set based on the statistical features, and initial cluster centers are obtained by screening the normalized region descriptors;

[0037] Density connectivity between the initial cluster centers is calculated, spatially adjacent initial cluster centers are merged based on the density connectivity, and a cluster center with the largest local density value and the strongest density connectivity after merging is determined as the position of the lesion core region.

[0038] In an optional embodiment,

[0039] A multi-scale radial basis function is constructed with the lesion core region as the center, a response distribution of the region descriptor under the multi-scale radial basis function is calculated, a spatial gradient contour of the response distribution is extracted, a spatial range of the lesion region is determined based on the contour, and a lesion expansion direction is determined by tracking the evolution of the contour through an optical flow field analysis, including:

[0040] A reference scale value is calculated according to the spatial range of the lesion core region with the center position of the lesion core region as the reference point, an incremental scale sequence is constructed based on the reference scale value according to a preset ratio, and a multi-scale radial basis function is constructed using the incremental scale sequence;

[0041] A response distribution of the region descriptor under the multi-scale radial basis function is calculated, an anisotropic weight matrix is generated based on the gradient of the structural features of the region descriptor, and the anisotropic weight matrix is multiplied by the response distribution to obtain a modified response distribution;

[0042] A spatial gradient contour of the modified response distribution is extracted at each scale, a stability index is calculated by calculating the spatial overlap ratio of the contour at adjacent scales, a change index is extracted by extracting the response gradient of the contour neighborhood, and a consistency index is calculated by calculating the difference in morphological features between the contour and the lesion core region;

[0043] The spatial gradient contour is optimized based on the stability index, the change index, and the consistency index, and the largest closed contour after optimization is determined as the spatial range of the lesion region;

[0044] The optical flow field vector of the spatial range boundary of the lesion area is calculated, 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.

[0045] In a second aspect of the embodiment of the present application, a real-time identification and positioning system for lesions on the inner wall of a milk duct based on fiber imaging is provided, comprising:

[0046] A first unit is configured to 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;

[0047] A second unit is configured to extract a displacement vector field between adjacent images in each subsequence, decompose the displacement vector field into a principal frequency component and a residual component in a frequency domain, obtain a tissue motion cycle through phase analysis of the principal 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;

[0048] 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;

[0049] 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 the 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 the lesion area based on the contour line, and track evolution of the contour line by optical flow field analysis to determine an expansion direction of the lesion area.

[0050] In a third aspect of the embodiment of the present application, an electronic device is provided, comprising:

[0051] a processor;

[0052] a memory for storing processor-executable instructions;

[0053] The processor is configured to invoke the instructions stored in the memory to execute the method described above.

[0054] In a fourth aspect of the embodiment of the present application, 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.

[0055] In the embodiment, the breast duct wall lesion real-time identification and positioning method based on fiber imaging can dynamically analyze the breast duct wall image sequence, identify the abnormal motion area through frequency domain decomposition of the displacement vector field, effectively capture the micro change characteristics of the potential lesion site, and improve the early detection sensitivity of the breast duct lesion. Through local wavelet analysis and regional descriptor construction of the tissue abnormal atlas, accurate characterization of different types of breast duct lesions is realized, the density clustering technology is combined to automatically identify the lesion core area, the subjectivity of manual interpretation in the traditional method is avoided, and the accuracy and objectivity of lesion identification are significantly improved. The present application adopts multi-scale radial basis function and spatial gradient contour line technology, accurately depicts the spatial range and expansion direction of the lesion area, provides three-dimensional positioning information of the lesion for clinicians, helps to develop more accurate treatment plan, reduces unnecessary tissue damage, and has important clinical application value for early diagnosis and minimally invasive treatment of breast diseases. BRIEF DESCRIPTION OF DRAWINGS

[0056] Figure 1 A flowchart of the breast duct wall lesion real-time identification and positioning method based on fiber imaging of the embodiment of the present application is shown in

[0057] Figure 2 A feature point motion trajectory analysis and abnormality detection schematic diagram of the embodiment of the present application is shown in

[0058] Figure 3 A lesion analysis effect comparison diagram of the embodiment of the present application is shown in DETAILED DESCRIPTION

[0059] In order to make the purpose, technical scheme and advantages of the embodiment of the present application clearer, the technical scheme in the embodiment of the present application will be described clearly and completely in combination with the drawings in the embodiment of the present application. Obviously, the described embodiments are only part of the embodiments of the present application, not all the embodiments. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art without creative labor are within the scope of protection of the present application.

[0060] The technical scheme of the present application will be described in detail in the following specific embodiments. The following several specific embodiments can be combined with each other, and the same or similar concepts or processes can not be described in some embodiments.

[0061] Figure 1 A flowchart of the breast duct wall lesion real-time identification and positioning method based on fiber imaging of the embodiment of the present application is shown in Figure 1 As shown, the method comprises:

[0062] Obtaining a breast duct wall image sequence and dividing it into multiple subsequences according to a preset time interval;

[0063] extracting displacement vector fields between adjacent images in each sub-sequence, decomposing the displacement vector fields into principal frequency components and residual components in the frequency domain, obtaining a tissue motion cycle through phase analysis of the principal frequency components; screening abnormal displacement regions from the residual components that do not conform to the tissue motion cycle, merging abnormal displacement regions in multiple sub-sequences, and establishing a tissue abnormality atlas;

[0064] blocking the tissue abnormality atlas, calculating local wavelet coefficients in each block region, grouping and counting the local wavelet coefficients to obtain internal structure features, calculating coefficient differences between adjacent blocks to obtain structure change features, and combining the internal structure features and the structure change features to form a region descriptor;

[0065] analyzing the region descriptor using a density clustering method to determine a lesion core region, constructing a multi-scale radial basis function centered on the lesion core region, calculating a response distribution of the region descriptor under the multi-scale radial basis function, extracting spatial gradient contour lines of the response distribution, determining a spatial range of the lesion region based on the contour lines, and determining a lesion extension direction by tracking evolution of the contour lines through optical flow field analysis.

[0066] In an alternative embodiment, obtaining a sequence of images of the inner wall of the breast duct and dividing the sequence into multiple sub-sequences according to a preset time interval comprises:

[0067] obtaining a sequence of continuous images of the inner wall of the breast duct, performing time domain analysis on the sequence of continuous images, extracting a periodic signal of the sequence of images, setting a time segmentation reference according to the periodic signal, taking an integer multiple of the time segmentation reference as a sub-sequence length, and dividing the sequence of continuous images into initial sub-sequences.

[0068] setting an overlap interval between adjacent sequences of the initial sub-sequences, the overlap interval having a length of a preset proportion of the sub-sequence length, calculating a weight value according to a time distance between an image and a boundary of the overlap interval, performing weighted processing on images in the overlap interval to obtain a transition sequence in the overlap interval, replacing image sequences in the original overlap interval with the transition sequence, and generating multiple image sub-sequences that are continuous in time sequence.

[0069] In the present embodiment, the medical device first obtains a sequence of continuous images of the inner wall of the breast duct through a fiber optic endoscope. In a specific embodiment, the fiber optic endoscope has a diameter of 0.5 millimeters and can smoothly enter the inside of the breast duct for imaging. The endoscope collects images at a frequency of 30 frames per second, recording the state of the tissue of the inner wall of the breast duct. These collected continuous images constitute a complete image sequence, containing all information of the inner wall of the breast duct.

[0070] The acquired continuous image sequence is analyzed in time domain, and periodic signals in the image sequence are extracted by detecting features such as pixel brightness change and contour line change. In actual application, the imaging of the inner wall of the breast duct will be affected by physiological activities such as patient breathing and heartbeat, and periodic image changes will occur. For example, in one collection, the patient's breathing period was detected to be about 4 seconds, and the heartbeat period was about 0.8 seconds. These periodic changes will affect the stability of the image.

[0071] The time segmentation reference is set according to the extracted periodic signal. In actual operation, the breathing period can be selected as the main reference because breathing has a greater impact on tissue movement. Assuming that the detected breathing period is 4 seconds, 4 seconds can be set as the time segmentation reference. An integer multiple of the time segmentation reference is selected as the subsequence length, for example, 1 times the reference, that is, 4 seconds, is selected as the subsequence length. In this way, the continuous image sequence is segmented every 4 seconds to form multiple initial subsequences. For a collection frequency of 30 frames per second, each subsequence contains 120 images.

[0072] In order to ensure smooth transition between adjacent subsequences, an overlap interval is set between adjacent sequences of the initial subsequence. The length of the overlap interval is a preset proportion of the subsequence length, for example, 25% of the subsequence length, that is, 1 second of overlap time. In this way, the two adjacent subsequences have 1 second of overlap in time, corresponding to 30 overlapping images.

[0073] The weight value is calculated according to the time distance of the image and the boundary of the overlap interval. In the overlap interval, the image closer to the end of the previous subsequence has a higher weight value, and the image closer to the beginning of the next subsequence has a lower weight value, showing a gradual transition characteristic. Specifically, for the ith frame of image in the overlap interval, if the overlap interval has n frames of image, the weight value of the image can be calculated according to the distance relationship. For example, the weight of the ith frame of image in the overlap interval of the previous subsequence can be set as (n-i) / n, and the weight of the ith frame of image in the overlap interval of the next subsequence can be set as i / n. Taking 30 overlapping images as an example, the weight of the first frame of image in the previous sequence is 29 / 30, and the weight of the first frame of image in the next sequence is 1 / 30; the weight of the 15th frame of image in the previous sequence is 15 / 30, and the weight of the 15th frame of image in the next sequence is 15 / 30; the weight of the 30th frame of image in the previous sequence is 0 / 30, and the weight of the 30th frame of image in the next sequence is 30 / 30.

[0074] The images in the overlapping interval are weighted and processed, and corresponding images in the front and rear sub-sequences are fused according to the weights to generate a transition sequence. During the weighting and processing, for each pixel point of an image, the pixel value of the corresponding image in the front sub-sequence and the pixel value of the corresponding image in the rear sub-sequence are calculated according to the respective weight values to obtain a weighted pixel value. For example, for the i-th image in the overlapping interval, if a pixel value of the corresponding image in the front sub-sequence is P1, a pixel value of the corresponding image in the rear sub-sequence at the same position is P2, the weight of the front sub-sequence is w1, and the weight of the rear sub-sequence is w2, then the weighted pixel value is P1 multiplied by w1 plus P2 multiplied by w2. Through such pixel-level weighted fusion, smooth transition of the images can be achieved.

[0075] The generated transition sequence is substituted for the image sequence in the original overlapping interval, thereby generating a plurality of time-sequentially continuous image sub-sequences. Finally, the entire continuous image sequence is divided into a plurality of sub-sequences with smooth transition characteristics, each sub-sequence has the same length, and smooth transition is achieved between adjacent sub-sequences through the weighting and processing, thereby avoiding abrupt phenomena that can be caused by simple splicing.

[0076] In actual applications, the system dynamically adjusts the parameters according to actual conditions. For example, for different patients, the respiratory cycle can be different, and the system adjusts the time segmentation reference according to the real-time detected respiratory cycle; for different examination sites, the proportion of the overlapping interval can also be adjusted as needed, and an appropriate value is usually selected between 15% and 35%.

[0077] In this embodiment, by acquiring a continuous image sequence and performing periodic signal analysis, this method can effectively eliminate the interference of physiological activities such as respiration and heartbeat on imaging, and improve the image stability. The adaptive adjustment of the time segmentation reference ensures that the sub-sequence division conforms to the physiological cycle characteristics, thereby ensuring the accuracy of image analysis. The overlapping interval setting and the weighting and processing technology achieve smooth transition between sub-sequences, avoid image abruptness caused by simple splicing, and enhance the sequence coherence. The pixel-level weighted fusion does not reduce the image quality in the transition area, and maintains the information integrity of the original image. This method can adapt to the physiological characteristic differences of different patients, and improves the adaptability and robustness through dynamic adjustment of parameters. The finally formed high-quality image sub-sequences lay a foundation for subsequent lesion recognition and positioning, significantly improve the detection rate and positioning accuracy of intraductal lesions, and enable doctors to make more accurate diagnosis and treatment planning, which has important value for improving the early diagnosis rate and treatment effect of breast diseases.

[0078] In an alternative embodiment, the displacement vector field between adjacent images in each sub-sequence is extracted, and the displacement vector field is decomposed into a main frequency component and a residual component in the frequency domain, which includes:

[0079] The gray scale values of adjacent image frames in the subsequence are acquired to construct an image block centered at a pixel point; a gray scale constancy constraint and a spatial consistency constraint are established within the image block to form a local gray scale variation constraint matrix, the local gray scale variation constraint matrix is substituted into an optical flow equation, and a pixel-level displacement vector is obtained by solving; the pixel-level displacement vector is regionally integrated to obtain a displacement vector field;

[0080] The displacement vector field is decomposed using a wavelet basis function to obtain displacement components of different scales; the amplitude distribution of each scale displacement component is calculated to generate a direction weight coefficient, and the displacement components are fused according to the direction weight coefficient to reconstruct the displacement vector field feature;

[0081] The displacement vector field feature is subjected to two-dimensional Fourier transform to obtain a frequency domain representation, the energy distribution of each frequency component in the frequency domain representation is counted, and the cumulative energy value of the frequency component is calculated; a frequency segmentation threshold is determined according to the cumulative energy value, the frequency domain representation is divided by the frequency segmentation threshold, and is reconstructed by inverse Fourier transform to obtain a main frequency component and a residual component.

[0082] In the specific embodiment, first, the gray scale values of adjacent image frames in the subsequence are acquired, and the color image is converted into a gray scale image to simplify the calculation. In actual application, the image of the inner wall of the breast duct obtained using an optical fiber endoscope is usually an RGB color image, which can be converted into a gray scale image according to the weight combination of R channel value multiplied by 0.299 plus G channel value multiplied by 0.587 plus B channel value multiplied by 0.114. For an original endoscope image with a resolution of 640×480, a gray scale image with the same resolution is obtained after conversion, and the gray scale value ranges from 0 to 255.

[0083] An image block is constructed centered at each pixel point, and the size of the image block is selected to be N×N pixels. In actual application, the value of N is usually set to an odd number between 5 and 15, for example, an image block of 11×11 can be selected. This size can contain enough local information without bringing too much computational burden. The system establishes a gray scale constancy constraint within each image block, which is based on the assumption that the gray scale value of the same tissue point in adjacent images should remain unchanged in a short time. For the inner wall tissue of the breast duct, this assumption is true in most cases, unless there is a significant change in illumination or tissue deformation.

[0084] Meanwhile, spatial consistency constraints are established within the image block, which are based on the assumption that the motion of neighboring pixels should be similar and not change drastically. This constraint is particularly suitable for breast duct wall tissue, as breast duct tissue usually moves as a whole, and the moving direction and amplitude of adjacent regions are often similar. Combining the gray constant constraint and the spatial consistency constraint forms a local gray change constraint matrix. For an 11x11 image block, this matrix describes the gray change relationship and spatial constraint relationship between 121 pixels in two adjacent frames.

[0085] The local gray 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 moving vector of the pixel can be obtained by solving the equation. In practical applications, an iterative solving method is used, usually setting the number of iterations to 3 to 5 times, and each iteration optimizes the results of the previous iteration 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.

[0086] 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, each grid contains 20x20 pixels. In each grid, the average value of all pixel displacement vectors 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.

[0087] 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 breast duct tissue.

[0088] The amplitude distribution of each scale displacement component is calculated to generate direction weight coefficients. Specifically, for each scale displacement component, the 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.

[0089] A two-dimensional Fourier transform is performed on the reconstructed displacement vector field to convert it from the spatial domain to the frequency domain representation. After the 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 ductal tissue, and the high frequency part corresponds to the local fine changes or noise.

[0090] The energy distribution of each frequency component in the statistical frequency domain representation is calculated, and the cumulative energy value of the frequency component is calculated. In practical applications, the proportion of the cumulative energy to the total energy is calculated in the order of frequency from low to high. For example, it can be found that the first 10% of the low frequency components have contained 85% of the total energy, indicating that the movement of the ductal tissue is mainly concentrated in the low frequency part, and the high frequency part has less energy.

[0091] The frequency segmentation threshold is determined according to the cumulative energy value. Generally, the frequency point at which the cumulative energy reaches 80% to 90% of the total energy can be taken as the segmentation threshold. In the above example, if the cumulative energy of the first 10% of the frequency components reaches 85% of the total energy, then 10% can be taken as the frequency segmentation threshold. With the frequency segmentation threshold as the boundary, the frequency domain representation is divided into a low frequency part and a high frequency part. The inverse Fourier transform is performed on the divided low frequency part and high frequency part respectively to convert it back to the spatial domain, and the main frequency component and the residual component of the displacement vector field are obtained. The main frequency component corresponds to the low frequency part and reflects the main movement mode of the ductal tissue, such as the overall movement caused by breathing; the residual component corresponds to the high frequency part and reflects the deformation of the local tissue or the characteristics of the potential lesion area.

[0092] In the identification of lesions in the inner wall of the breast duct, 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 the subsequent analysis focuses on the changes of the tissue itself rather than the artifacts caused by movement. For example, the system can use the displacement field calculated by the main frequency component to register adjacent images, eliminating the influence of overall movement caused by factors such as breathing, heartbeat, etc. The residual component contains the characteristic information of the potential lesion area, because the lesion tissue often shows different movement characteristics from the surrounding normal tissue. For example, the hardness of the lesion tissue such as intraductal papilloma is usually different from that of normal tissue, and it shows different deformation modes during the movement of the endoscope, and these differences will be reflected in the residual component.

[0093] In this embodiment, the micro-motion characteristics of the intraductal wall tissue are accurately captured by pixel-level displacement vector calculation and region integration, overcoming the problem that traditional methods are difficult to handle image changes caused by tissue motion. Wavelet basis function decomposition technology can effectively separate tissue motion of different scales, enabling the system to analyze both overall movement and local deformation characteristics. The introduction of directional weight coefficients enhances the expression of the main motion direction and suppresses the interference of the secondary direction, improving the accuracy of displacement field reconstruction. Frequency domain decomposition technology successfully separates tissue motion into main frequency components and residual components, enabling the system to distinguish between physiological motion and abnormal motion patterns related to lesions. This decomposition method significantly improves the sensitivity and specificity of lesion identification, reducing the false positive detection rate. At the same time, this technology can adapt to complex lighting conditions and tissue deformation within the duct, maintaining high robustness and providing reliable technical support for early and accurate identification and positioning of intraductal wall lesions.

[0094] In an alternative embodiment, the tissue motion cycle is obtained by phase analysis of the main frequency component; abnormal displacement regions that do not conform to the tissue motion cycle are selected from the residual component, and the abnormal displacement regions in multiple sub-sequences are merged to establish a tissue abnormality map, including:

[0095] Performing Hilbert transform on the main frequency component, extracting the instantaneous phase, and calculating the time interval sequence between adjacent extreme points in the phase sequence, the tissue motion cycle and its fluctuation range are determined according to the time interval sequence;

[0096] In the residual component, the grid points are divided, the corner point response values of the grid points are calculated, the grid points with corner point response values greater than a set response threshold are selected as feature points, the correlation coefficients between the feature points in adjacent frames are calculated for feature point matching, and the position offset of the feature points is obtained based on the matching results;

[0097] 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, and the arc length change and curvature change are compared with the tissue motion cycle to mark the trajectory points that exceed the cycle fluctuation range, and the region containing the marked points is extracted as an abnormal displacement region;

[0098] The motion direction, speed and acceleration characteristics of the abnormal displacement region are calculated to construct the motion characteristics, the similarity between the abnormal displacement regions is calculated according to the motion characteristics, the abnormal displacement regions with 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 abnormality map.

[0099] In implementation, the Hilbert transform is first 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 practical application, the Hilbert transform can be performed on the horizontal and vertical displacement components in the dominant frequency component respectively to obtain the instantaneous phase in two directions. For the breast duct wall tissue, the dominant frequency component usually presents periodic changes due to the influence of respiration and heartbeat, and these changes can be reflected by the trend of the instantaneous phase change.

[0100] The time sequence variation 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 period 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 periods of the patient.

[0101] The time interval sequence between adjacent extreme points is calculated, and the tissue motion period and its fluctuation range are determined according to the time interval sequence. In the breast duct wall imaging process, respiration and heartbeat are the two main periodic factors affecting tissue motion. By analyzing the distribution characteristics of the time interval sequence, these two periods can be identified. For example, for adult patients, the respiration period is usually in the range of 3 to 5 seconds, and the heartbeat period is usually in the range of 0.6 to 1.2 seconds. In statistical analysis, the histogram method or density estimation method can be used to find the peak value of the time interval distribution as the main period of tissue motion. At the same time, the standard deviation of each period component is calculated to determine the fluctuation range of the period. For example, the respiration period can be identified as 4.2 seconds with a fluctuation range of ±0.5 seconds, and the heartbeat period can be identified as 0.8 seconds with a fluctuation range of ±0.1 seconds.

[0102] Grid points are divided in the residual component, and the corner response value of the grid points is calculated. The residual component reflects the local irregular motion of the tissue, and may contain potential lesion information. In order to effectively analyze the residual component, the image can be uniformly divided into grids, for example, a grid point is set every 10 pixels. For each grid point, the gray level gradient covariance matrix of the surrounding area is calculated, and the corner response value is calculated based on the matrix. The corner response value reflects the richness of the texture in the local area, and the larger the value, the richer the feature information contained in the area.

[0103] The grid points with corner response values greater than a set response threshold are selected as feature points. The response threshold can be dynamically set according to the overall characteristics of the image, for example, the threshold can be set as the average value of the corner response values of all grid points plus 1.5 times the standard deviation. For the intraductal wall image, the feature points are usually distributed in the areas with rich tissue texture or obvious edges. In an intraductal wall image with a resolution of 640x480, 100 to 200 feature points can be selected for subsequent analysis.

[0104] The cross-correlation coefficients of the 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 result. In specific implementation, a fixed-size image block, for example, a 15x15-pixel region, can be extracted around the feature point, and the cross-correlation coefficient of the image block with the possible corresponding position in the next frame is calculated. The cross-correlation coefficient reflects the similarity of the two image blocks, and the value closer to 1 indicates a higher matching degree. The corresponding position of the feature point can be determined by searching for the position with the maximum cross-correlation coefficient in the next frame, and thus the position offset of the feature point is calculated. To improve the matching accuracy, a pyramid matching strategy can be used, that is, coarse matching is performed on low-resolution images first, and then accurate positioning is performed on high-resolution images.

[0105] The position offsets are concatenated along the time dimension and filtered to obtain the motion trajectory of the feature points. Since noise can be introduced in the image acquisition and feature point matching processes, the original offsets need to be filtered. Median filtering or Gaussian filtering methods can be used to eliminate abnormal jumps in the trajectory. The filtered trajectory reflects the true motion of the intraductal wall tissue and provides a reliable basis for subsequent analysis.

[0106] The tangent vector and normal vector of the motion trajectory are calculated to obtain the arc length variation and curvature variation of the trajectory. The tangent vector reflects the direction of the motion of the feature point, and the normal vector reflects the change in the direction of the motion. Based on the tangent vector, the arc length variation, that is, the distance moved by the feature point in unit time, can be calculated; based on the normal vector, the curvature variation, that is, the bending degree of the motion trajectory, can be calculated. Under the influence of breathing and heartbeat, the motion trajectory of normal breast duct tissue usually exhibits regular arc length and curvature variations.

[0107] The arc length variation and curvature variation are compared with the tissue motion period to mark out the trajectory points that exceed the period fluctuation range. The motion of normal tissue should be synchronized with the breathing and heartbeat period, and the arc length and curvature variations should be within the period fluctuation range. If the arc length variation or the curvature variation of the motion trajectory of a certain feature point at a certain time significantly exceeds the period fluctuation range, it may indicate that there is an abnormality in that region. For example, if the arc length variation period of a certain feature point is 4.3 seconds, which is consistent with the identified breathing period of 4.2 seconds ± 0.5 seconds, but suddenly becomes 6.2 seconds for a certain period of time, the trajectory of that period is marked as abnormal.

[0108] The region containing the marked points is extracted as an abnormal displacement region. A spatial clustering method can be used to gather feature points that are close in space and exhibit abnormalities at the same time, forming an abnormal displacement region. For example, a density clustering algorithm can be used, and a spatial distance threshold of 20 pixels is set. Abnormal 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 range, time span, and number of abnormal feature points contained is recorded.

[0109] The motion direction, speed, and acceleration characteristics of the abnormal displacement region are calculated 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 region; the speed can be obtained by calculating the displacement amount per unit time; and the acceleration can be obtained by calculating the rate of change of the speed. These characteristics together constitute the motion characteristic description of the abnormal displacement region. For example, a certain abnormal displacement region may exhibit motion in the direction of the nipple, with an average speed of 2.5 pixels / frame and an acceleration of 0.3 pixels / frame 2 .

[0110] The similarity between abnormal displacement regions is calculated based on the motion characteristics. The Euclidean distance or cosine similarity method can be used to calculate the distance between different abnormal displacement regions in the motion characteristic space. The smaller the distance, the higher the similarity. For example, the motion direction, speed, and acceleration can be normalized to construct a feature vector, and the Euclidean distance between the vectors can be calculated as a similarity index.

[0111] The abnormal displacement regions with a similarity greater than a preset similarity threshold and adjacent in space are subjected to connected component analysis. The preset similarity threshold can be set to 0.8 or higher to ensure that the merged regions have highly similar motion characteristics. Connected component analysis can use region growing or morphological processing methods to connect adjacent regions that meet the conditions to form larger connected regions.

[0112] The connected components are integrated into an organizational abnormality map. During the integration process, the spatial range, time span, and motion characteristic information of each connected component are retained to form a complete abnormality map. This map visually displays the regions in the breast duct wall tissue that may have lesions, providing an important diagnostic reference for doctors. For example, in the abnormality map, a certain part of the breast duct may exhibit a motion pattern different from that of normal tissue, which may be a manifestation of intraductal papilloma or intraductal carcinoma.

[0113] In the present embodiment, the instantaneous phase is extracted by Hilbert transform of the main frequency component, the scheme can accurately identify the periodic motion characteristics of the intraductal wall tissue, and overcome the defects that the prior art can not distinguish physiological motion from pathological motion only by image gray level change. The prior art usually uses fixed threshold or empirical parameters to judge tissue abnormalities, which is easily disturbed by physiological activities such as breathing and heartbeat, leading to misdiagnosis. The present scheme calculates the tissue motion cycle and fluctuation range, establishes a dynamic reference standard, and improves the accuracy of abnormal region detection. In the residual component, the feature points are selected by the corner response value and the cross-correlation coefficient is calculated for matching, which solves the problem of uneven illumination and large view angle change of the intraductal wall, and enhances the robustness of feature tracking. Comparing the changes of the arc length and curvature of the feature point trajectory with the tissue motion cycle can effectively identify abnormal regions that do not conform to the physiological cycle, and reduce the false positive detection rate. Through motion feature similarity calculation and connected component analysis, the abnormal displacement region is integrated to form a complete tissue abnormality atlas, which improves the detection rate and positioning accuracy of intraductal wall lesions, and provides reliable technical support for early clinical diagnosis.

[0114] Figure 2 For the feature point motion trajectory analysis and abnormality detection schematic diagram of the embodiment of the present application, as shown in Figure 2 , the figure directly shows the process of feature point motion trajectory analysis and abnormality detection. Four feature points are clearly marked in the figure, which are located at coordinates (7.5, 11.0), (15.0, 20.0), (25.0, 15.0) and (31.0, 25.0), and 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 region. By calculating the tangent vector and normal vector of the trajectory, analyzing the changes of arc length and curvature, the present technical scheme successfully identifies multiple abnormal displacement points, such as the feature point (7.5, 11.0) which has an abnormality of 3.0 deviation in its third trajectory, and the feature point (15.0, 20.0) which has an abnormality of 5.0 deviation in its third trajectory. These abnormal points are all beyond the fluctuation range of the normal motion cycle of the tissue, and by marking and connecting these abnormal points, a clear abnormal displacement region is formed. The present technical scheme has high precision and high sensitivity in trajectory analysis, and can accurately capture the tiny abnormal motion pattern.

[0115] In an alternative embodiment, the tissue abnormality atlas is divided into blocks, the local wavelet coefficients in each block region are calculated, the local wavelet coefficients are grouped and counted to obtain internal structure features, the coefficient difference degree between adjacent blocks is calculated to obtain structure change features, and the internal structure features and structure change features are combined to form a region descriptor, including:

[0116] The tissue abnormality map is subjected to recursive quadtree partitioning, information entropy of pixel gray scale distribution in the partitioned region is calculated, and whether to continue to subdivide the current block is determined according to a comparison result of the information entropy and a partition threshold, until all blocks meet a stop condition, and an adaptive partition map is obtained;

[0117] Each partitioned region in the adaptive partition map is subjected to wavelet decomposition, and detail coefficient matrices in horizontal, vertical and diagonal directions are obtained; a gray scale gradient distribution of the partitioned region is calculated, a direction weight coefficient is generated, a statistic of the detail coefficient matrix is calculated and weighted by the direction weight coefficient, and an internal structure feature vector of the partition is constructed;

[0118] Based on the adaptive partition map, a partition adjacency relationship is determined, adjacent partition pairs are subjected to scale normalization processing, a difference degree of local wavelet coefficients is calculated, a structure change matrix is constructed according to the difference degree, a propagation feature of the structure change matrix at a region boundary is calculated, the propagation feature is accumulated and counted, and a structure change feature is generated;

[0119] The internal structure feature and the structure change feature are aligned in dimension, a covariance matrix of the aligned features is calculated, feature orthogonalization decomposition is performed based on the covariance matrix, feature components are sorted and screened according to a decomposition result, and a region descriptor is constructed by combining the screened feature components.

[0120] For example, first, the tissue abnormality map is subjected to recursive quadtree partitioning, and the structural differences of the intraductal wall tissue are fully reflected by the adaptive partitioning. 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 partition. In the recursive partitioning process, each partition is equally divided into four sub-blocks, forming a quadtree structure. For each partitioned 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 partition 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.

[0121] The adaptive segmentation process is continued until all the blocks meet the stopping condition, and the final adaptive segmentation map is obtained. In practical applications, the normal tissue regions of the breast duct wall usually form larger blocks, while the lesion regions are often segmented into multiple small blocks due to high texture complexity.

[0122] Each block region in the adaptive segmentation map is subjected to wavelet decomposition to obtain the 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 practical applications, multi-level decomposition can be performed, usually 2 to 3 levels, to capture texture features at 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 breast duct wall tissue in different directions, which are of great significance for identifying the lesion tissue of the breast duct wall.

[0123] The gray level gradient distribution of the block region is calculated to generate the direction weight coefficient. The gray level gradient reflects the direction and intensity of pixel value change, which can be obtained by calculating the difference between adjacent pixels. For each block region, the gradients in horizontal and vertical directions are calculated, and then the gradient amplitude and direction are obtained. The distribution of gradient direction is calculated, 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.

[0124] The statistical quantities are calculated for the detail coefficient matrix and weighted by the directional weight coefficients to construct the internal structure feature vector of the block. The statistical quantities can include mean, variance, skewness, kurtosis, etc. These statistical quantities 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 quantities of the detail coefficients in each direction are calculated, and then combined by weighting according to the directional weight coefficients 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 area is usually larger than that of the normal tissue, reflecting the irregularity of the texture in the lesion area.

[0125] The block adjacency relationship is determined based on the adaptive block atlas, and the scale normalization processing is performed on the adjacent block pairs. The adjacency relationship can be determined by analyzing the spatial position of the block. If two blocks share a boundary, they are considered to be adjacent. Since the adaptive segmentation results in different sizes of adjacent 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 blocks of different sizes to a uniform size, such as 16x16 pixels.

[0126] The difference degree of the local wavelet coefficients is calculated, and the structure change matrix is constructed according to the difference degree. For each pair of adjacent blocks, the difference in the wavelet coefficients 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 blocks in the horizontal, vertical, and diagonal directions can be calculated to form a difference vector. The difference vectors of all adjacent 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 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.

[0127] The propagation feature of the structure change matrix at the region boundary is calculated, and the propagation feature is accumulated to generate the structure change feature. The propagation feature reflects how the structural difference spreads from one region to the surrounding region. By treating the blocks as nodes and the adjacency relationship as edges, the propagation of the difference degree from one block to other blocks can be calculated by graph theory methods. 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. This 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 intraductal wall lesions, this feature helps to distinguish the lesion boundary and the internal region, and improves the positioning accuracy.

[0128] The internal structure features and the structure change features are aligned in dimension, and a covariance matrix is calculated for the aligned features. Since the internal structure features and the structure change features can be different in dimension, a dimension alignment process is needed. The two types of features 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 a covariance matrix of the vector is calculated. The covariance matrix reflects the correlation between the feature components, which helps to eliminate redundant information.

[0129] Based on the covariance matrix, feature orthogonalization decomposition is performed, and the feature components are sorted and selected according to the decomposition results. The 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 size of the feature value reflects the amount of information contained in the principal component. The principal components are sorted according to the size of the feature value, and the principal components with larger contribution rates are retained, while 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.

[0130] The selected feature components are combined to construct a region descriptor. The region descriptor is a compact representation of the abnormal region features 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 intraductal wall lesions can be achieved.

[0131] In this embodiment, the recursive quadtree partitioning dynamically determines the block size according to the information entropy, so that the complex regions can obtain more detailed blocks, and the simple regions can maintain the integrity. The wavelet decomposition captures the multi-scale texture features, and the weighted processing of the directional weight coefficients enhances the sensitivity to the structure changes in different directions of the intraductal wall. The calculation of the local wavelet coefficient difference between blocks realizes the accurate expression of the boundary features of the tissue, so that the boundary features of the lesions and normal tissues can be effectively extracted. The propagation feature analysis of the structure change matrix can accurately describe the transition characteristics between the internal structure of the lesion region and the surrounding normal tissue, and improve the recognition ability of the fuzzy region of the lesion boundary. The feature orthogonalization decomposition eliminates the redundant information and reduces the calculation complexity, and the constructed region descriptor has high discriminative ability, so that the system can realize accurate identification and positioning of different types of intraductal wall lesions, and provide reliable technical support for early clinical diagnosis.

[0132] In an alternative embodiment, a density clustering method is used to analyze the region descriptor to determine the lesion core region, including:

[0133] The region descriptor is normalized to a preset numerical interval by maximum minimum value mapping to obtain a normalized region descriptor, a search radius is determined according to a 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;

[0134] A bandwidth parameter of an adaptive kernel function is determined according to a median of the distance matrix, a local density value of the normalized region descriptor is calculated by using the adaptive kernel function, a target region descriptor is determined according to the local density value, a 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;

[0135] Statistical features of the density-distance distribution map are extracted, an adaptive density threshold and a distance threshold are set based on the statistical features, and an initial clustering center is obtained by screening the normalized region descriptor;

[0136] Density connectivity between the initial clustering centers is calculated, spatially adjacent initial clustering centers are merged based on the density connectivity, and a clustering center with the maximum local density value and the strongest density connectivity after merging is determined as the position of the lesion core region.

[0137] Exemplarily, first, the region descriptor is 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 may have large differences, 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. During normalization, 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 in the maximum minimum value range. For example, if the minimum value of a certain dimension feature is 10, the maximum value is 50, 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, and the normalized region descriptor in a unified range can be obtained, which is convenient for subsequent analysis.

[0138] The search radius is determined according to the feature distribution of the normalized region descriptors. The search radius is a key parameter of density clustering, and affects the calculation accuracy of local density. An excessively large search radius will cause over-smoothing and loss of local features, while an excessively small search radius will cause unstable density estimation. In the identification of intraductal wall lesions, the search radius can be dynamically determined according to the distribution characteristics of the normalized region descriptors. 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 can be selected as the search radius. In practical applications, for a data set containing 100 normalized region descriptors, the calculated search radius may be 0.15, indicating that about 15% of the sample pairs in the feature space have a distance less than this value.

[0139] The neighbor set of each normalized region descriptor is obtained within the search radius range, 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 metric methods such as Euclidean distance, Manhattan distance, or Mahalanobis distance. In the identification of intraductal wall lesions, Euclidean distance is commonly used, 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 sufficiently large value, indicating that there is no direct connection between them.

[0140] The bandwidth parameter of the adaptive kernel function is determined according to 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. Commonly used kernel functions include Gaussian kernel, Epanechnikov kernel, etc. The bandwidth parameter controls the shape of the kernel function, directly affecting the smoothing degree of density estimation. In the identification of intraductal wall lesions, an adaptive bandwidth strategy can be used to determine the bandwidth parameter according to the statistical characteristics of the distance matrix. Specifically, the median of the effective distance values (distance less than the search radius) in the distance matrix can be calculated as the bandwidth parameter. This method can automatically adjust the bandwidth according to 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.

[0141] The local density value of the normalized region descriptor is calculated using the adaptive kernel function. The local density value reflects the degree of aggregation of sample points in the feature space, and is the basis of density clustering. For each normalized region descriptor, the contribution of its neighbor points is calculated, and all contributions are accumulated to obtain the local density value. Specifically, the Gaussian kernel function can be used, and for two sample points with a distance of d, the contribution value is exp(-d 2 / h2 ), where h is a bandwidth parameter. In the identification of intraductal lesions, the descriptors of normal tissue regions usually have a scattered distribution and a low local density value; while the descriptors of lesion regions often form a dense distribution and have a high local density value due to the high similarity of features.

[0142] Descriptors with high local density values are selected as target region descriptors. Target region descriptors are the descriptors with high local density values, which may correspond to the feature centers of lesion regions. According to the distribution characteristics of local density values, descriptors with density values greater than a certain threshold can be selected as target region descriptors. For example, the top 20% of descriptors in terms of local density values, or the descriptors with density values greater than the average density plus one standard deviation, can be selected. In the identification of intraductal lesions, target region descriptors usually correspond to typical features of lesion tissue and are important for the determination of lesion core regions.

[0143] The minimum distance value of each normalized region descriptor to a target region descriptor is calculated as the relative distance. The relative distance reflects the proximity of the sample point to the high-density region and is an important basis for judging whether the sample point is a cluster center. For each normalized region descriptor, the distance to all target region descriptors is calculated, and the minimum value is taken as the relative distance. In the identification of intraductal lesions, the descriptors of the lesion edge region usually have moderate local density values and small relative distances, while the descriptors of the lesion core region often have high local density values and large relative distances.

[0144] A density-distance distribution map is constructed according to the local density values and relative distances. The density-distance distribution map is a two-dimensional scatter plot, with the horizontal axis representing the local density value and the vertical axis representing the relative distance. Each normalized region descriptor corresponds to a point in the map. This visualization method can intuitively show the clustering characteristics of sample points. In an ideal case, cluster centers appear as points with high density and high distance, which are clearly distinguished from other sample points. In the identification of intraductal lesions, the density-distance distribution map can help identify the feature centers of the lesion core region.

[0145] Statistical features of the density-distance distribution map are extracted, and adaptive density and distance thresholds are set based on the statistical features. Statistical features can include the mean, variance, skewness of the density distribution, the mean, variance, skewness of the distance distribution, and the correlation coefficient of density and distance, etc. These features reflect the overall distribution characteristics of the data set and are helpful in determining appropriate thresholds. In practical applications, the density threshold can be set to the mean density plus the density standard deviation multiplied by a coefficient, and the distance threshold can be set to the mean distance plus the distance standard deviation multiplied by a coefficient. The value of the coefficient can be adjusted according to actual needs, usually between 0.5 and 2. For example, if the mean density is 0.3, the standard deviation is 0.1, and the coefficient is 1.5, then the density threshold is 0.3 + 0.1 x 1.5 = 0.45.

[0146] The normalized region descriptors are screened to obtain initial cluster centers. The screening condition is that the descriptors with local density value greater than the density threshold and relative distance greater than the distance threshold are selected as initial cluster centers. These descriptors have both high local density, indicating that the region samples where they are located are clustered, and large relative distance, indicating that they are significantly distinguished from other high-density regions. In the identification of intraductal wall lesions, the initial cluster centers usually correspond to different types or different parts of lesion features. For example, for intraductal wall images containing intraductal papilloma, multiple initial cluster centers may be identified, corresponding to different parts or different growth stages of the lesion.

[0147] 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 density connectivity refers to the density continuity between two cluster centers, reflecting whether they belong to the same cluster. For two initial cluster centers, the minimum density value of all points on the path between them is calculated as the density connectivity. If the density connectivity is greater than a set threshold, it indicates that there is a density continuous path between the two centers, which may belong to the same lesion region and should be merged. The density connectivity threshold can be set to a certain proportion of the mean density of all samples, for example, 80% of the mean density. Merging spatially adjacent initial cluster centers can obtain a more compact cluster structure. In the identification of intraductal wall lesions, the merging process helps to eliminate the situation that the same lesion region is identified multiple times, and improves the accuracy of lesion positioning.

[0148] The cluster center with the largest local density value and the strongest density connectivity after merging is determined as the lesion core region position. The density connectivity refers to the degree of density connection of the cluster center to other sample points, which can be measured by calculating the number of density connected paths from the cluster center to other sample points. The stronger the density connectivity, the more concentrated the feature distribution of the region where the cluster center is located, and the more likely it is the core region of the lesion. In practical applications, the local density value and the density connectivity can be considered comprehensively, for example, the weighted sum of the two is calculated, and the weight can be adjusted according to actual needs. The cluster center with the largest weighted sum is selected as the lesion core region position. In the identification of intraductal wall lesions, the lesion core region usually corresponds to the most typical and concentrated part of the lesion tissue, which is the key to accurate positioning and diagnosis.

[0149] In this embodiment, the maximum-minimum value mapping normalization technique eliminates the dimensional differences between the dimensions of the regional descriptor, improving the accuracy of distance calculation in the feature space. The adaptive search radius determination strategy dynamically adjusts the parameters according to the feature distribution, enhancing the algorithm's adaptability to lesions of different density distributions. The adaptive kernel function bandwidth parameter based on the median of the distance matrix makes the density estimation more accurate, effectively dealing with complex situations such as uneven lighting and tissue deformation in the intraductal breast wall image. The joint analysis method of local density value and relative distance successfully distinguishes the core area and the edge area of the lesion, reducing the false positive rate. The statistical feature analysis of the density-distance distribution graph realizes the adaptive adjustment of the threshold parameter, making the algorithm maintain stable performance among different patients and different lesion types. The density connectivity calculation and spatial proximity cluster center merging mechanism effectively solve the problem of multiple identifications of the same lesion, improving the positioning accuracy. This method has excellent recognition ability for irregular-shaped and various-sized intraductal breast wall lesions, providing reliable technical support for accurate positioning of early micro-lesions, and significantly improving the sensitivity and specificity of intraductal breast wall lesion diagnosis.

[0150] In an alternative embodiment, a multi-scale radial basis function is constructed centered on the core area of the lesion, 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 lesion extension direction is determined by tracking the evolution of the contour through optical flow field analysis, including:

[0151] Taking the center position of the core area of the lesion as the reference point, calculating the reference scale value based on the spatial range of the core area of the lesion, constructing an increasing scale sequence based on the reference scale value according to a preset ratio, and constructing a multi-scale radial basis function using the increasing scale sequence;

[0152] Calculating the response distribution of the regional descriptor under the multi-scale radial basis function, generating an anisotropic weight matrix based on the gradient of the structural features of the regional descriptor, and multiplying the anisotropic weight matrix and the response distribution to obtain a modified response distribution;

[0153] Extracting the spatial gradient contour of the modified response distribution at each scale, calculating the spatial overlap ratio of the contour at adjacent scales to obtain a stability index, extracting the response gradient of the contour neighborhood to obtain a change index, and calculating the morphological feature difference between the contour and the core area of the lesion to obtain a consistency index;

[0154] Optimizing the spatial gradient contour based on the stability index, the change index, and the consistency index, and determining the largest closed contour after optimization as the spatial range of the lesion area;

[0155] The optical flow field vector of the spatial range boundary of the lesion area is calculated, 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 extension direction of the lesion area.

[0156] Exemplarily, a reference scale value is calculated according to the spatial range of the lesion core area with the center position of the lesion core area as a reference point. In actual application, the lesion core area is a feature aggregation area identified from the region descriptor by a density clustering method. The center position of the lesion core area can be obtained by calculating the average coordinates of all points in the area, and is taken as the reference point of the radial basis function. The spatial range of the lesion core area can be determined by calculating the extension distance of the area 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 area can be taken as the reference scale value. For example, if the lesion core area presents an approximate circular distribution in space with an average radius of 5 pixels, the reference scale value can be set to 5.

[0157] An incremental scale sequence is constructed according to a preset ratio based on the reference scale value, and the multi-scale radial basis function is constructed 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 a 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 area. 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 common forms include Gaussian function, multi-quadratic function, etc. In the identification of intraductal wall lesions, a Gaussian radial basis function can be used, and the function value decreases with the increase of the distance between the point and the reference point. Specifically, for a distance d and a scale s, the function value can be represented as exp(-d 2 / s 2 ). The set of radial basis functions constructed by different scales forms a multi-scale radial basis function, which can capture the feature distribution of the lesion area from different spatial scales.

[0158] The region descriptor describes the response distribution of the region under the multi-scale radial basis function. The region descriptor contains the structural features of the abnormal region of the inner wall of the breast duct. When combined with the multi-scale radial basis function, the distribution of the features at different spatial scales can be obtained. In specific implementation, for the region descriptor of each spatial position, the product of the region descriptor and the radial basis function is calculated to obtain the response value. The response values of all positions are organized into a response distribution map, which reflects the distribution intensity of the features in space. In the identification of lesions in the inner wall of the breast duct, the response value of the lesion region is usually higher than that of the surrounding normal tissue, and the highlight region is shown in the response distribution map.

[0159] An anisotropic weight matrix is generated based on the gradient of the structural features of the region descriptor, and the modified response distribution is obtained by multiplying the anisotropic weight matrix and the response distribution. The anisotropic weight matrix is used to adjust the weight of the response distribution in different directions, so that the modified response distribution better reflects the directional features of the tissue structure. The gradient of the structural features reflects the trend of the features in space, which can be obtained by calculating the partial derivative of the region 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 decomposed to obtain the principal direction and the degree of anisotropy, and the anisotropic weight matrix is generated accordingly. In the identification of lesions in the inner wall of the breast duct, the lesion tissue often shows obvious structural features in certain directions. Through the adjustment of the anisotropic weight matrix, the response in these directions can be enhanced, and the identification accuracy of the lesion boundary can be improved.

[0160] The spatial gradient contour of the modified response distribution is extracted at each scale. The spatial gradient contour refers to the curve with equal gradient amplitude of the response distribution, which usually corresponds to the boundary region of the response distribution. The edge detection and contour extraction method can be used to extract the spatial gradient contour. In the identification of lesions in the inner wall of the breast duct, the gradient amplitude map of the modified response distribution can be calculated first, and then threshold segmentation is applied on the gradient amplitude map to extract the point set with gradient amplitude equal to a certain value, and these points are connected into a closed curve, which is the spatial gradient contour. Usually, multiple contours corresponding to different gradient values can be extracted to form a contour set. For example, the 25%, 50% and 75% quantile points of the gradient amplitude can be selected as the threshold to extract three groups of contours.

[0161] 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.

[0162] 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.

[0163] 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.

[0164] The stability index, the change index, and the consistency index are used to optimize the spatial gradient contour lines. The maximum closed contour line after optimization is determined as the spatial range of the lesion region. The optimization process can use a weighted scoring method to consider the values of the three indexes comprehensively. For example, the three indexes can be normalized to the interval [0, 1], and the weights w1, w2, and w3 are set. 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 region. In the identification of intraductal wall lesions, the spatial range of the lesion region usually appears as a closed boundary, which contains the entire region of the lesion tissue.

[0165] The spatial range boundary of the lesion region is calculated to obtain the optical flow field vector. 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 region. In a continuous sequence of intraductal wall images, the optical flow field vector is extracted for the determined boundary of the lesion region. The Lucas-Kanade method or the Horn-Schunck method can be used to calculate the motion vector based on the image gray gradient and the 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 the boundary optical flow field. For example, in a sequence of 100 images, the optical flow field can be calculated every 5 frames to obtain the motion trend of the lesion boundary at different times.

[0166] The principal component analysis is performed on the optical flow field vector, and the direction of the eigenvector corresponding to the largest eigenvalue obtained by the principal component analysis is determined as the expansion direction of the lesion region. The principal component analysis is a dimension reduction method that can be used to extract the main variation direction of the data. The principal component analysis is performed on the data matrix composed of the boundary optical flow field vector to calculate the eigenvalues and eigenvectors of the covariance matrix. The eigenvector represents the main direction of data variation, and the eigenvalue represents the variation degree in that direction. The eigenvector corresponding to the largest eigenvalue is selected, and its direction is the main expansion direction of the lesion region. In the identification of intraductal wall lesions, the identification of the lesion expansion direction is of great significance for evaluating the development trend of the lesion and formulating a treatment plan. For example, if the largest eigenvector obtained by the principal component analysis is at 45 degrees, it indicates that the lesion region mainly expands to the right and upward direction.

[0167] In this embodiment, an adaptive scale sequence is constructed using the core region of the lesion as a reference point, overcoming the limitation of traditional fixed-scale methods in adapting to lesions of different sizes. The introduction of anisotropic weighting matrices effectively enhances the structural features of lesion tissue in specific directions, improving the accuracy of lesion boundary identification. Multi-scale spatial gradient contour extraction technology can capture lesion boundary information from different spatial scales, while the comprehensive evaluation mechanism of stability, change, and consistency indices effectively selects the optimal contour lines, 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 region, providing clinicians with a basis for predicting lesion development.

[0168] Figure 3 This is a comparison chart of the lesion analysis effects in an embodiment of the present invention, such as... Figure 3 As shown in the figure, this comparison illustrates the performance evaluation results of three different lesion analysis methods. The multi-scale radial basis function method performed excellently across all five evaluation indicators, particularly achieving high scores of 92.3% and 91.2% in boundary accuracy and clinical applicability, respectively, significantly outperforming traditional methods. The traditional gradient method performed moderately, with all indicators evenly distributed between 72% and 85%. While the region growing method showed good computational efficiency (90.6%), it was relatively weaker in other indicators. This comparative analysis confirms that the lesion analysis method based on multi-scale radial basis functions has significant advantages in accurately capturing lesion boundaries and adapting to the clinical environment. In particular, through the comprehensive evaluation of stability, change, and consistency indicators, it can more accurately determine the spatial extent and expansion direction of the lesion region.

[0169] A second aspect of the present invention provides a real-time identification and localization system for lesions on the inner wall of the mammary duct based on fiber optic imaging, the system comprising:

[0170] The first unit is used to acquire the image sequence of the inner wall of the mammary duct and divide it into multiple sub-sequences according to a preset time interval;

[0171] The second unit is used to extract the displacement vector field between adjacent images in each subsequence, decompose the displacement vector field into a main frequency component and a residual component in the frequency domain, obtain the tissue motion cycle through phase analysis of the main frequency component, screen out abnormal displacement regions that do not conform to the tissue motion cycle from the residual component, merge the abnormal displacement regions in multiple subsequences, and establish a tissue abnormality map.

[0172] The third unit is used to divide the tissue abnormality map into blocks, calculate the local wavelet coefficients in each block, group and statistically analyze the local wavelet coefficients to obtain internal structural features, calculate the coefficient difference between adjacent blocks to obtain structural change features, and combine the internal structural features and structural change features to form a region descriptor.

[0173] A fourth unit is configured to analyze the region 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 region 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 determine a lesion expansion direction by tracking evolution of the contour line through an optical flow field analysis.

[0174] In a third aspect, the present application provides an electronic device, comprising:

[0175] a processor;

[0176] a memory for storing processor-executable instructions;

[0177] The processor is configured to invoke the instructions stored in the memory to execute the method described above.

[0178] In a fourth aspect, the present application provides a computer-readable storage medium having stored thereon computer program instructions, which, when executed by a processor, implement the method described above.

[0179] 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.

[0180] 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 for 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 localization of lesions on the inner wall of mammary ducts based on fiber optic imaging, characterized in that, include: Acquire a sequence of images of the inner wall of the mammary ducts and divide it into multiple sub-sequences according to a preset time interval; The displacement vector field between adjacent images in each subsequence is extracted, and the displacement vector field is decomposed into a main frequency component and a residual component in the frequency domain. The tissue motion cycle is obtained by phase analysis of the main frequency component. Abnormal displacement regions that do not conform to the tissue motion cycle are screened out from the residual component. Abnormal displacement regions in multiple subsequences are merged to establish a tissue abnormality map. The tissue anomaly map is divided into blocks, the local wavelet coefficients in each block are calculated, the local wavelet coefficients are grouped and statistically analyzed to obtain internal structural features, the coefficient difference between adjacent blocks is calculated to obtain structural change features, and the internal structural features and structural change features are combined to form a region descriptor. Density clustering was used to analyze the region descriptors to determine the core region of the lesion. A multi-scale radial basis function was constructed with the core region of the lesion as the center. The response distribution of the region descriptors under the multi-scale radial basis function was calculated. Spatial gradient contour lines of the response distribution were extracted. The spatial range of the lesion region was determined based on the contour lines. The direction of lesion expansion was determined by tracking the evolution of the contour lines through optical flow field analysis.

2. The method according to claim 1, characterized in that, The image sequence of the inner wall of the mammary duct is obtained and divided into multiple subsequences according to a preset time interval, including: A continuous image sequence of the inner wall of the mammary duct is obtained, and a time-domain analysis is performed on the continuous image sequence to extract the periodic signal of the image sequence. A time segmentation benchmark is set according to the periodic signal, and the continuous image sequence is divided into initial subsequences with an integer multiple of the time segmentation benchmark as the 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 duration of the subsequence. A weight value is calculated based on the time 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 in the frequency domain into a main frequency component and a residual component, including: Obtain the grayscale values ​​of adjacent image frames in the subsequence, and construct image blocks centered on pixels; establish grayscale constantness constraints and spatial consistency constraints within the image blocks to form a local grayscale change constraint matrix; substitute the local grayscale change constraint matrix into the optical flow equation to solve for pixel-level displacement vectors; and perform regional integration on the pixel-level displacement vectors to obtain a displacement vector field. The displacement vector field is decomposed using wavelet basis functions to obtain displacement components at different scales; the amplitude distribution of each scale displacement component is calculated to generate directional weighting coefficients; the displacement components are fused based on the directional weighting coefficients to reconstruct 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 components is calculated. A frequency segmentation threshold is determined based on the cumulative energy value. The frequency domain representation is divided by the frequency segmentation threshold and reconstructed by inverse Fourier transform to obtain the dominant frequency component and the residual component.

4. The method according to claim 1, characterized in that, The tissue motion cycle is obtained through phase analysis of the dominant frequency components; abnormal displacement regions that do not conform to the tissue motion cycle are screened from the residual components, and abnormal displacement regions from multiple subsequences are merged to establish a tissue anomaly map, including: The Hilbert transform is applied to the main frequency component to extract the instantaneous phase. The temporal changes of the instantaneous phase are statistically analyzed to identify the extreme points in the phase sequence. The time interval sequence between adjacent extreme points is calculated, and the tissue movement cycle and its fluctuation range are determined based on the time interval sequence. Divide the residual components into grid points, calculate the corner response values ​​of the grid points, select grid points whose corner response values ​​are greater than a set response threshold as feature points, calculate the cross-correlation coefficient of the feature points between adjacent frames to perform feature point matching, and obtain the position offset of the feature points based on the matching results. The position offset is concatenated along the time dimension and filtered to obtain the motion trajectory of the feature point. The tangent vector and normal vector of the motion trajectory are calculated to obtain the arc length change and curvature change of the trajectory. The arc length change and curvature change are compared with the tissue motion cycle to mark the trajectory points that exceed the cycle fluctuation range. The area containing the marked points is extracted as the abnormal displacement area. The motion direction, velocity, and acceleration characteristics of the abnormal displacement region are calculated to construct motion features. The similarity between abnormal displacement regions is calculated based on the motion features. Connectivity analysis is performed on abnormal displacement regions with similarity greater than a preset similarity threshold and spatially adjacent, and the connected regions are integrated into an organization anomaly map.

5. The method according to claim 1, characterized in that, The tissue anomaly map is divided into blocks, and the local wavelet coefficients within each block are calculated. These local wavelet coefficients are then grouped and statistically analyzed to obtain internal structural features. The coefficient differences between adjacent blocks are calculated to obtain structural change features. Finally, the internal structural features and structural change features are combined to form a region descriptor, including: The abnormal tissue map is recursively segmented into a quadtree, and the information entropy of the pixel grayscale distribution within the segmented region is calculated. Based on the comparison between the information entropy and the segmentation threshold, it is determined whether to continue subdividing the current block until all blocks meet the stopping condition, thus obtaining an adaptive segmented map. Wavelet decomposition is performed on each block region in the adaptive block map to obtain the detail coefficient matrix in the horizontal, vertical and diagonal directions; the gray-level gradient distribution of the block region is calculated to generate directional weight coefficients; statistics are calculated on the detail coefficient matrix and weighted by the directional weight coefficients to construct the internal structural feature vector of the block. Based on the adaptive block map, the block adjacency relationship is determined, the scale normalization process is performed on the adjacent block pairs, the difference degree of local wavelet coefficients is calculated, the structure change matrix is ​​constructed according to the difference degree, the propagation characteristics of the structure change matrix at the region boundary are calculated, the propagation characteristics are accumulated and statistically analyzed to generate the structure change characteristics. The internal structural features and structural change features are aligned by dimension. The covariance matrix is ​​calculated for the aligned features. Based on the covariance matrix, the features are orthogonally decomposed. The feature components are sorted and filtered according to the decomposition results. The filtered feature components are combined to construct a region descriptor.

6. The method according to claim 1, characterized in that, Density clustering was used to analyze the region descriptors to determine the core lesion regions, including: The region descriptor is normalized to a preset numerical range by mapping the maximum and minimum values ​​to obtain the normalized region descriptor. The search radius is determined according to the feature distribution of the normalized region descriptor. The nearest neighbor set of each normalized region descriptor is obtained within the search radius. The distance matrix is ​​calculated based on the nearest neighbor set. The bandwidth parameter of the adaptive kernel function is determined based on the median of the distance matrix. The local density value of the normalized region descriptor is calculated using the adaptive kernel function. The target region descriptor is determined based on the local density value. The minimum distance from each normalized region descriptor to the target region descriptor is calculated as the relative distance. A density-distance distribution map is constructed based on the local density value and the relative distance. Statistical features of the density-distance distribution map are extracted, and adaptive density and distance thresholds are set based on the statistical features. Normalized region descriptors are then filtered to obtain initial cluster centers. Calculate the density connectivity between the initial cluster centers, merge spatially adjacent initial cluster centers based on the density connectivity, and determine the location of the core lesion region as the merged cluster center with the largest local density value and the strongest density connectivity.

7. The method according to claim 1, characterized in that, A multi-scale radial basis function is constructed centered on the core region of the lesion. The response distribution of the region descriptor under the multi-scale radial basis function is calculated, and spatial gradient contour lines of the response distribution are extracted. Based on the contour lines, the spatial extent of the lesion region is determined, and the direction of lesion expansion is determined by tracking the evolution of the contour lines through optical flow field analysis. Using the center of the core lesion region as a reference point, a reference scale value is calculated based on the spatial range of the core lesion region. An incremental scale sequence is constructed based on the reference scale value according to a preset ratio, and a multi-scale radial basis function is constructed using the incremental scale sequence. The response distribution of the region descriptor under multi-scale radial basis functions is calculated, and an anisotropic weight matrix is ​​generated based on the structural feature gradient of the region descriptor. The anisotropic weight matrix is ​​multiplied by the response distribution to obtain the modified response distribution. Spatial gradient contour lines of the corrected response distribution are extracted at each scale. The spatial overlap ratio of contour lines at adjacent scales is calculated to obtain a stability index. The response gradient of the neighborhood of the contour lines is extracted to obtain a change index. The morphological feature difference between the contour lines and the core area of ​​the lesion is calculated to obtain a consistency index. Based on the stability index, change index and consistency index, the spatial gradient contour lines are optimized, and the optimized maximum closed contour line is determined as the spatial range of the lesion area. The optical flow field vector is calculated for the spatial boundary of the lesion area. Principal component analysis is performed on the optical flow field vector, and the direction of the eigenvector corresponding to the largest eigenvalue obtained from the principal component analysis is determined as the expansion direction of the lesion area.

8. A real-time identification and localization system for ductal wall lesions based on fiber optic imaging, used to implement the method of any one of claims 1-7, characterized in that, include: The first unit is used to acquire the image sequence of the inner wall of the mammary duct and divide it into multiple sub-sequences according to a preset time interval; The second unit is used to extract the displacement vector field between adjacent images in each subsequence, decompose the displacement vector field into a main frequency component and a residual component in the frequency domain, obtain the tissue motion cycle through phase analysis of the main frequency component, screen out abnormal displacement regions that do not conform to the tissue motion cycle from the residual component, merge the abnormal displacement regions in multiple subsequences, and 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, group and statistically analyze the local wavelet coefficients to obtain internal structural features, calculate the coefficient difference between adjacent blocks to obtain structural change features, and combine the internal structural features and structural change features to form a region descriptor. The fourth unit is used to analyze the region descriptors using density clustering to determine the core region of the lesion. A multi-scale radial basis function is constructed with the core region of the lesion as the center. The response distribution of the region descriptors under the multi-scale radial basis function is calculated. The spatial gradient contour lines of the response distribution are extracted. The spatial range of the lesion region is determined based on the contour lines. The direction of lesion expansion is determined by tracking the evolution of the contour lines through optical flow field analysis.

9. An electronic device, characterized in that, include: processor; Memory used to store processor-executable instructions; The processor is configured to invoke 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 the processor, they implement the method described in any one of claims 1 to 7.

Images