Methods and systems for lesion characterisation

Abstract

Systems and methods for delineating tumour boundaries using a three-dimensional (3D) image stack, such as photoacoustic image slices, comprising a plurality of images. Multiple maximum intensity projections (MIPs) are derived from the stack, and level set segmentation is performed on those MIPs, to derive a boundary that bounds a region of interest (i.e., a tumour boundary). Structural and functional information is then derived for the region of interest, to analyse the tumour tissue.

Description

METHODS AND SYSTEMS FOR LESION CHARACTERISATIONTechnical Field

[0001] The present invention relates, in general terms, to the field of medical imaging, and more particularly, to methods and systems for precise lesion characterisation and three-dimensional (3D) visualisation via an automated photoacoustic imaging analysis technique.Background

[0002] Tumour boundary delineation is a critical aspect, especially in the field of oncological diagnostics. Within this domain, traditional methods, such as histopathology, have long been employed for tumour boundary assessment. However, there are several limitations associated with these traditional methods.

[0003] First, traditional methods tend to be reliant on invasive procedures. They may also lack real-time imaging capabilities, impeding the dynamic capture of changes in tumour boundaries crucial for adapting treatment strategies. Second, the two-dimensional nature of histopathology typically offers limited spatial representation, potentially leading to an incomplete understanding of three-dimensional tumour structures. The post-mortem or biopsy nature of specimens not only restricts real-time monitoring, but also renders the process time-consuming and labour-intensive. Subjectivity in histopathological image interpretation introduces inter-observer variability, posing challenges to standardisation and potentially resulting in inconsistent tumour boundary delineation among pathologists.

[0004] Existing medical imaging technologies play a crucial role in diagnosing and monitoring various diseases, particularly in the context of lesion detection. Despite the advancement of established imaging modalities, such as Optical Coherence Tomography (OCT) and Reflectance Confocal Microscopy (RCM),accurate lesion boundary delineation remains a persistent challenge. These imaging techniques provide inadequate resolution for fine structures, hindering precise lesion boundary definition and impacting diagnostic accuracy. Moreover, these modalities lack the incorporation of functional information, such as tissue chromophores or metabolic activity, thereby limiting their ability to comprehensively characterise lesions. The limited incorporation of functional features, coupled with sensitivity to noise and artefacts, presents challenges in achieving comprehensive lesion characterisation.

[0005] There is hence a need for an improved imaging technique, such as a photoacoustic technique coupled with a robust automated tumour boundary delineation algorithm, to overcome the above-mentioned limitations, and to enhance the precision and efficiency of lesion delineation for improved patient care and outcomes.Summary

[0006] Disclosed is a system for delineating tumour boundaries, comprising: a receiving module, for receiving an input of a three-dimensional (3D) image stack comprising a plurality of images; a preprocessing module, for generating a plurality of maximum intensity projections from the 3D image stack; a level set segmentation module, for initialising a set of boundary estimates for bounding a region of interest in each image, and iteratively performing level set segmentation to update the boundary estimates until convergence at a final boundary bounding the region of interest; and a post processing module, for extracting structural and functional information from the region of interest.

[0007] Also disclosed herein is a method for delineating tumour boundaries, comprising : receiving an input of a three-dimensional (3D) image stack comprising a plurality of images; generating a plurality of maximum intensity projections from the 3D image stack;initialising a set of boundary estimates for bounding a region of interest in each image, and iteratively performing level set segmentation to update the boundary estimates until convergence at a final boundary bounding the region of interest; and extracting structural and functional information from the region of interest.Brief description of the drawings

[0008] Embodiments of the present invention will be described with reference to the following figures, in which:

[0009] Figure 1 shows an image analysis workflow of the algorithm of the invention, from an initial 3D image stack obtained from Multispectral Optoacoustic Tomography (MSOT) to a final set of structural and functional information obtained after analysis.

[0010] Figure 2 shows a 3D reconstruction of segmented 2D basal cell carcinoma (BCC) image slices, to obtain a volume of the tumour, as compared to an original 3D render as collected by the MSOT.[Oil] Figure 3 shows a comparison of MSOT-obtained measurements with histology. Figure 3(a) illustrates the correlation of tumour width between MSOT and histology measurements. Figure 3(b) illustrates the correlation of tumour depth between MSOT and histology measurements.

[0012] Figure 4 is a schematic of a computer system for implementing the workflow of Figure 1.Detailed Description

[0013] The following detailed description is merely exemplary in nature and is not intended to limit the invention or the application and uses of the invention. Furthermore, there is no intention to be bound by any theory presented in the preceding background of the invention or the following detailed description.

[0014] Disclosed herein are systems and methods for delineating tumour boundaries. Advantageously, the claimed invention provides a multi-step segmentation process for multi-modal photoacoustic images, an in-depth structural feature extraction, capturing intricate details within lesions, functional information extraction (for e.g., variations in tissue vascularity and metabolic activity), a three-dimensional (3D) volumetric reconstruction of organs and tumours, and an automation of algorithmic processes for the application of the method in photoacoustic imaging.

[0015] Figure 1 shows an image analysis workflow 100 of the algorithm of the invention, for delineating tumour boundaries. The workflow 100 progresses from an initial 3D image stack to a final set of structural and functional information obtained after analysis. In step 110, a three-dimensional (3D) image stack is collected. The image stack may be obtained from Multispectral Optoacoustic Tomography (MSOT) techniques, 3D optoacoustic tomography or another 3D imaging modality. The 3D image stack may be 100 x 100 x 100 pixels, or of other dimensions. It may contain a single feature or a plurality of features.

[0016] In step 120, image pre-processing is carried out. A series of maximum intensity projections (MIPs) is generated for each anatomical plane (step 122), or a single projection may be generated per plane. In some embodiments, three anatomical planes are used, namely the coronal plane, also known as the frontal or Y-X plane, the transverse plane, also known as the axial or X-Z plane, and the sagittal plane, also known as the lateral or Y-Z plane.

[0017] As presently used, MIP involves parallel ray tracing using rays normal to the imaging plane. Each ray passes through one or more images in the 3D stack and the highest intensity encountered by that ray will appear in the projection. The visualisation of relevant features is enhanced during the generation of maximum intensity projections.

[0018] Embodiments of the methods described herein make use of radiopacity changes, or emission / absorption spectra changes, resulting from chromophores - particularly, endogenous chromophores. Chromophoresimprove the visibility of the features sought to be identified - e.g., tumour tissue. Moreover, since different tissue types and structures produce different endogenous chromophores, detecting particular emission spectra or radiopacity characteristics can facilitate tissue differentiation. When detecting, for example, signals modified by chromophores from tumour tissue, the emission spectra of tumour cells will modify the reflected optoacoustic signal, resulting in greater signal intensity for tumour cells when the optoacoustic signal is tuned to the relevant emission spectra.

[0019] Since the intention of maximum intensity projection is to emphasise features relevant for tumour boundary delineation, selection of contrast medium may be based on the particular anatomical features that medium preferentially enters and the imaging modality.

[0020] Preprocessing the 3D image stack may also involve noise reduction. Noise reduction may be performed on the maximum intensity projections (step 124). This may be done by applying a median filter, a non-linear digital filtering technique to remove noise. Noise reduction may also be performed by averaging or Gaussian filters. Noise reduction serves to refine the image and improve the outcomes of subsequent steps.

[0021] Thresholding (step 126) is also performed during image pre-processing. Thresholding aims to eliminate portions of the image which fall below a predetermined percentage threshold relative to the maximum intensity, whether from 0 intensity to the maximum intensity, or from the minimum intensity and the maximum intensity, of pixels in the projections generated at step 122. This may also eliminate irrelevant information such as stray rays or pixel-to-pixel variation. Thresholding may result in the image being partitioned into a foreground and a background. Particularly where contrast- enhanced features are being identified, relevant image features should be significantly more prominent (i.e., correspond to higher intensity pixels or voxels) than otherwise confounding features. Therefore, for all voxels shown in the projections, there may be multiple peaks. For example, in a graph showing a number of voxels of each intensity from lowest intensity to highest intensity, there may be a peak in the voxels that have not absorbed contrast medium and another peak in the voxels that have absorbed contrast medium.The threshold may therefore be set between the two peaks, to remove a large number of voxels that are not of interest.

[0022] In step 130, level set segmentation is performed. A number of initialisations (step 132) of the cell boundary estimate are first performed, based on the maximum intensity of the image. In some embodiments, a series of five initialisations are performed. In other embodiments, two, three, four, or six or more initialisations may be performed, and the number of initialisations can depend on the particular application or desired accuracy of the level set segmentation process. An initialisation step serves to guide the workflow toward potential boundary locations within the 3D image stack. This provides a series of diverse starting points for adaptability to different or varying lesion structures.

[0023] In step 134, level set segmentation is performed by applying iterations of the level set function. In effect, a curve is adjusted in the image until it meets some optimal condition corresponding to when the curve coincides with the boundary of the region of interest - e.g., a tumour. The curve encloses the region of interest, isolating or distinguishing it from the background. The boundary of the region of interest corresponds to the zero-level set, which is where the signed distance function of points on the evolved boundary (i.e., the boundary after a series of iterations of the level set function) is zero. The curve may be defined by one or more criteria, including that the curve should be smooth, evolves according to a local mean curvature (e.g., using mean curvature flow) and its progress is damped by intensities of pixels in the region of interest. Evolution of the boundary towards the optimal boundary (i.e., when the boundary coincides with that of the tumour or region of interest) may be guided by an energy minimisation process and an implicit representation of the evolving boundary, to precisely capture the true boundaries of the segmented region. Level set segmentation may comprise a progressive evaluation of the differences among neighbouring pixels to find object boundaries. This may involve considering the image gradients, edges, or region intensities, to determine the object boundaries. The cell boundary estimate is also updated.

[0024] In some embodiments, when performing step 130, the method may involve calculating the intensity gradient across the MIP (or across each slice if 3D reconstruction is intended). The boundaries are thus initialized using the intensity gradient from the contrast-enhanced images. Specifically, the segmentation process begins by dividing the image into a plurality of subdomains corresponding to 0%, 20%, 40%, 60%, and 80% of the maximum intensity within the images. Each value corresponds to an initialisation of the level set function, as discussed below. The segmentation process is based on dividing the domain into N = 5 disjoint subdomains (or another number of subdomains, as appropriate). Each subdomain is assigned an initial constant corresponding to specific intensity thresholds at 0%, 20%, 40%, 60%, and 80% of the maximum intensity within the images. These constants serve as initial values to guide the evolution process. This initialization ensures that the level set model is robust to variability in signal intensity and enables it to accurately delineate regions of interest.

[0025] Per step 134, the level set function iteratively evolves to minimize the energy function and accurately converge on the true boundary of the lesion (which may be referred to as the "cell boundary" which, in this case, is the boundary between the tumour tissue and healthy or surrounding tissue rather than the boundary of an anatomical 'cell' per se). This dynamic approach adjusts for the spatial variations in intensity and morphology without relying on fixed bounding boxes.

[0026] To accurately define the lesion in 3D, segmentation is performed directly on the MIPs - e.g. the MIP for each anatomical plane. As discussed above, the MIP highlights regions of interest across the entire image stack. This approach significantly streamlines the process by providing a comprehensive view of the lesion without segmenting each slice individually. For analyses requiring detailed 3D reconstructions, such as volume estimation, segmentation is performed on a slice-by-slice basis. This ensures accurate capture of lesion morphology for volumetric measurements while maintaining overall consistency across slices.

[0027] The dimensions of the estimated boundaries are then dynamically determined during the evolution of the level set function. The evolving surfaceiteratively adapts to the morphological features of the lesion, ensuring an accurate fit to the true boundary. The energy function drives this evolution, with convergence indicating that the estimated boundary aligns closely with the actual lesion boundary.

[0028] Segmentation performed directly on the MIPs (the three, orthogonal 2D projections) delineates the boundaries of lesions - e.g., pigmented BCC lesions. In some embodiments, the level set segmentation algorithm undergoes multi-partitioning and multi-labelling, based on a convex relaxation of the continuous Potts model, iteratively using the curve of the level set function over time to compute the optimal partitions and accurately adapt to the object boundaries. "Multi-labelling" refers to the process of assigning distinct labels to different partitions / subdomains in the domain of interest during the segmentation process. In the context of the level set algorithm, this is achieved by evolving the level set function iteratively over time. The level set function is used to represent multiple partitions implicitly, where each label corresponds to a specific partition.fi = \I - Ci\r(3)

[0029] The energy function in Equation 1 is composed of the sum of multipartition and the boundary lengths of these partitions - boundary lengths are calculated by counting boundary pixels. The sum of multi-partition is the summation of energy contributions from each of the N disjoint partitionsin the segmentation process. The boundary lengths d0zare calculated, or approximated, by counting the pixels along the edges separating adjacent partitions. This discrete approximation serves as a regularization term in the energy function, encouraging the formation of smooth and well-defined boundaries between the partitions. In the formulae above, I represent the input image, which is divided into N disjoint partitions <pt, <ptdenote the characteristic functions of the disjoint subdomains with constrains listed inEquation 2, c = {c1;c2, ..cN} are the constants for the initialization. The level set function ft = | / - ct\rcan have coefficients y = 1 or 2 and 2 can be empirically determined or manually adjusted to find the optimal partitions. The boundary lengths of the subdomains are given by d0(. By minimizing the energy function, the level set function f (Equation 3) evolves to fit the optimal object boundary. The energy function is considered 'minimized' if the change of the output of the function between successive iterations is below a predetermined threshold.

[0030] The segmentation process begins with initializing an estimate of the cell boundary into N = 5 disjoint subdomains, with initial constants c, incorporating intensity thresholds corresponding to 0%, 20%, 40%, 60%, and 80% of the maximum intensity within the images. The level set function ftis used to iterate the evolving surface of the tumour to separate it from the background signals. Subsequently, it is iteratively updated to minimize the energy function that quantifies the difference between the evolving boundary and the actual boundary of the lesion. Thus, evolution involves region-based partitioning, which constrains the intensity differences within each partition (region) and boundary regularization in which the boundary lengths of the partitions are minimized, to encourage smooth and well-defined boundaries.

[0031] In step 136, the system is configured to check for boundary convergence. The convergence criteria ensures that the algorithm stops evolving when an optimal boundary has been attained. If boundary convergence has not yet been attained, the level set function is performed, and the cell boundary estimate is updated. Steps 132, 134 and 136 may be repeated until boundary convergence is attained.

[0032] More specifically, convergence of the boundary estimate is checked at each iteration against the convergence of the energy function. Once convergence is achieved, the final cell boundary estimate is obtained, providing an accurate representation of the BCC lesion morphology within the MSOT images. In this instance, convergence is the point at which the energy function stabilizes - i.e., there is minimal change between iterations where the change in the energy function. This stabilization may be found where the change between successive iterations drops below a predetermined threshold.That threshold can be manually defined to be a level of change below which further iteration will yield insignificant benefit.

[0033] Change below the predetermined threshold indicated that the level set function has evolved to its optimal configuration, effectively delineating the boundary of the region of interest. In this context, "optimal" may not be perfectly optimal, but is sufficiently optimal that further evolution of the boundary will yield insignificant benefit in delineating the tumour boundary.

[0034] As mentioned above, boundary convergence occurs at the zero level set of the energy function, where the signed distance function equals zero. When equal to zero, the signed distance function shows that movement of the boundary in any direction will result in a worse approximation of the boundary of the region of interest. This ensures that the boundary separates the region of interest from the background, with no significant deviation along the boundary points.

[0035] Since multiple initializations are used, and in the circumstance that the slices collectively contain a single continuous tumour mass, convergence carries the further property that all initial estimates have evolved to collectively define a consistent boundary that encloses the same region. The energy minimization process ensures stability, preventing the presence of multiple distinct "optimal" boundaries (i.e., instability). This stability is inherent to the iterative refinement, as the energy function inherently balances regional homogeneity and boundary smoothness, driving all partitions toward a unified solution.

[0036] The level set segmentation is designed to delineate tumour boundaries rather than individual cell boundaries. While the initialization process begins with a coarse estimate of the lesion boundary, the iterative refinement ensures that the final segmentation represents the tumour morphology accurately. The coarse estimate can essentially be arbitrary since convergence of the boundary on the line of minimum energy will occur from any starting point. This focus on tumour boundaries enables subsequent analyses, such as measuring tumour size, volume, and texture features.

[0037] The level set segmentation method integrates both region partitioning and constrained boundary information, providing robust and precise boundary detection. This approach is particularly valuable for accurately quantifying tumour characteristics within MSOT images. The iterative refinement process ensures reproducibility and precision in segmentation.

[0038] At the end of the level set segmentation, a final cell boundary estimate is produced (step 138). A final cell boundary is produced for each anatomical plane and comprises one intensity plot for the respective anatomical plane. The "intensity plot" is the set of pixel intensities of pixels within the MIP, that lie inside the final boundary for each anatomical plane. This plot represents the spatial distribution of the intensity values, derived from the MIP of the respective anatomical plane, and provides information about the signal characteristics within the delineated region of interest - i.e., within the boundary of the respective MIP.

[0039] In step 140, data post-processing is carried out. This involves extracting both structural and functional information from the segmented regions. Other information, such as texture features and tumour volume may also be extracted during data post-processing.

[0040] With regard to tumour volume, an estimate may be made from segmentation performed on MIPs as set out above. This can be sufficient for some analyses, and some tumour types. However, when segmentation is performed on MIP images, the resulting boundary reflects the projection of the highest intensity signals through the 3D image stack for each anatomical plane. This approach is efficient for estimating lesion boundaries but has inherent limitations, as regions within the defined boundary may not fully correspond to tumour tissue.

[0041] For detailed 3D reconstructions, segmentation is conducted on individual slices, generating one boundary per slice. This method ensures accurate lesion morphology in all three anatomical planes. To handle planes perpendicular to the slices, boundaries are interpolated between adjacent slices, maintaining continuity and ensuring the 3D boundary reflects the true lesion morphology while minimizing the inclusion of non-tumour regions.

[0042] MIP-based segmentation provides a rapid and practical solution for initial boundary estimation, while slice-by-slice segmentation offers higher precision for applications such as volumetric analysis. Combining these methods balances efficiency and accuracy, allowing segmentation to be tailored to specific analytical needs.

[0043] Structural information extraction may comprise calculating the width and depth of the feature of interest (step 142a). A first MIP can be used to estimate the tumour width along two axes and the depth will be estimated from the other two MIPs, and the greatest of the two estimates of the depth is taken to avoid underestimation of tumour depth. Knowledge of the width and depth of the segmented region may allow for determination of a maximum feret diameter, or it may allow for determination of the length of the major and minor axes.

[0044] It may also comprise extracting the texture of the feature of interest (step 142b). This may be done using a Gray Level Co-occurrence Matrix (GLCM), which quantifies the spatial relationships between pixel intensities in an image. The GLCM method quantifies textural patterns present in the image, with all background pixels (i.e., pixels outside the boundary determined at step 130) being set to zero and GLCM being applied only to non-zero pixels, and may reveal information about the homogeneity, contrast and randomness of the underlying tissue structures. By considering the spatial relationships between pixels, the GLCM enhances the algorithm's ability to discriminate between different tissue textures within the segmented regions, offering a more comprehensive understanding of the internal composition of lesions. In this context, "textures" refer to the spatial relationships between pixel intensities and their surroundings. By analyzing features such as homogeneity, contrast, and entropy, this approach has the potential to differentiate between various subtypes of tumour tissue (e.g., basal cell carcinoma) based on their distinct textural patterns. This may also result in more accurate diagnosis and treatment planning, in medical imaging applications.

[0045] Each segmented 2D layer may be reconstructed into a 3D volumetric representation of the lesion, for visualisation by a user. This facilitatespreoperative mapping and enhances the precision of surgical planning. Figure 2 shows a 3D reconstruction of segmented 2D basal cell carcinoma (BCC) image slices, to obtain a volume of the tumour, as compared to an original 3D render as collected by the MSOT.

[0046] Functional information may also be extracted. This may include extraction of tissue chromophore information (step 144). Tissue chromophore information is obtained using spectral unmixing, involving separating and quantifying tissue chromophores (e.g., haemoglobin, melanin) shown in photoacoustic images by analyzing their distinct absorption spectra across multiple wavelengths. Identifying each absorption spectrum enables the visualization of chromophore distributions within the boundary identified at step 138, providing insights into tissue composition. In view of the present teachings, it will be evident to the skilled person how to perform spectral unmixing from images obtained from a photoacoustic system. Thus, the method 100 can involve extracting functional information by spectral unmixing of absorption spectra in the photoacoustic images, isolating one or more tissue chromophore channels or spectra, and calculating the average intensity within the segmentation mask or boundary. This ensures an incorporation of functional features, such as variations in tissue vascularity or metabolic activity, enhancing the algorithm's capacity to comprehensively characterise lesions.

[0047] The workflow set out above is embodied by the system 400 in Figure 4. An imager 402 uses one or more imaging modalities to capture the 3D image stack. In some embodiments, the 3D image stack is a series of 2D images, or "slices", that are stacked atop one another to form the stack and the depth of the stack (i.e., passing through each image, perpendicular to a plane of the image) provides the third dimension. When performing operations on the stack, spacing between successive slices may be based on the position, through the patient, at which the image was taken.

[0048] The imager 402 passes the captured 3D stack to the system 400 for processing. The preprocessor 404 receives the 3D image stack and applies MIP module 406 to perform MIP on the 3D image stack. The preprocessor 404 may also use one or both of Noise reducer 408 for denoising the images inthe 3D image stack and / or denoising the MIPs, and the thresholder 410 for thresholding voxels to distinguish between voxels of interest (or "foreground") and voxel that are not of interest (or "background").

[0049] After preprocessing, the preprocessed 3D image stack is sent to the level set segmenter 412 for boundary determination. The level set segmenter 412 uses the cell boundary initializer 414 to generate initial estimates of the boundary of the region of interest. The boundary updater 416 applies a level set function iteration to the initial estimates in the first instance, and updated estimates for each iteration, to determine an updated estimate for each initial estimate of the boundary. To check whether the boundary has been optimized, a convergence checker 418 checks whether the updated boundaries resulting from all estimates have converged to a common boundary.

[0050] After defining the boundaries, a post-processor 420 processes the bounded images - i.e., the region of each image, within the respective boundary. The post- processor 420 extracts structural and textural features. It uses a width and depth calculator 422 for large structural elucidation, texture extractor 424 for applying a texture algorithm (e.g, Gray Level Co-occurrence Matrix) to extract textural features and a chromophore extractor 426 that extracts chromophore features useful in understanding the function of tissues in the region or volume of interest.

[0051] All of the above operations are governed by a program code 432, stored in memory 430 and executed by one or more processors 436.

[0052] While the system 400 has been shown with a specific structure, some modules referenced in Figure 4 may be implemented in hardware, software, firmware or a combination thereof, and may be in a single server or distributed computing environment. All such hardware variations are intended to fall within the scope of the present teachings.Example 1

[0053] To validate the algorithm, a cohort of 31 clinical subjects diagnosed with pigmented basal cell carcinoma (BCC) were studied. Figure 3 shows acomparison of MSOT-obtained measurements with histology. Figure 3(a) illustrates the correlation of tumour width between MSOT and histology measurements. Figure 3(b) illustrates the correlation of tumour depth between MSOT and histology measurements. A robust correlation of 0.84 (100% correlation having a value of 1) for the tumour width and a correlation of 0.81 for the tumour depth were obtained, for studies carried out using the systems and methods of the present invention. This signifies a strong alignment between the algorithm's prediction and the clinically observed dimensions of pigmented BCC lesions, where a correlation value of 1 would indicate a 100% alignment between the algorithm's prediction and the clinically observed dimensions.

[0054] Traditional methods for measuring tumour width and depth may be performed by image processing programs, such as ImageJ, which involves a user manually drawing a line across an observed tumour on the image, for determining the greatest tumour width. This may be subject to inter-observer variability.

[0055] It will be appreciated that may further modifications and permutations of various aspects of the described embodiments are possible. Accordingly, the described aspects are intended to embrace all such alterations, modifications and variations that fall within the spirit and the scope of the appended claims.

[0056] Throughout this disclosure, unless the context requires otherwise, the word "comprise" and variations such as "comprises" and "comprising" will be understood to imply the inclusion of a stated integer or step or group of integers or steps, but not the exclusion of any other integer or step or group of integers or steps.

[0057] The reference to any prior art in this disclosure is not, and should not be taken as, an acknowledgement or any form of suggestion that the prior art forms part of the common general knowledge.

Claims

CLAIMS1. A system for delineating tumour boundaries, comprising: a receiving module, for receiving an input of a three-dimensional (3D) image stack comprising a plurality of images; a preprocessing module, for generating a plurality of maximum intensity projections from the 3D image stack; a level set segmentation module, for initialising a set of boundary estimates for bounding a region of interest in each image, and iteratively performing level set segmentation to update the boundary estimates until convergence at a final boundary bounding the region of interest; and a post processing module, for extracting structural and functional information from the region of interest.

2. The system of claim 1, wherein each of the plurality of maximum intensity projections is generated for a respective anatomical plane.

3. The system of claim 1 or claim 2, wherein the preprocessing module further comprises a median filter applied to the 3D image stack, for noise reduction.

4. The system of any one of claim 1 to claim 3, wherein the level set segmentation module initializes the boundary estimates using a maximum intensity pixel, with each boundary estimate differing from all other boundary estimates.

5. The system of any one of claim 1 to claim 4, wherein level set segmentation is guided by an energy minimisation process.

6. The system of any one of claim 1 to claim 5, wherein level set segmentation stops when the updated boundary estimates converge on a common boundary.

7. The system of any one of claim 1 to claim 6, wherein extracting structural information comprises one or more than one of: extracting width and depth measurements, determining a maximum feret diameter, and extracting a length of major and minor axes of the region of interest.

8. The system of any one of claim 1 to claim 7, wherein extracting structural information comprises analysing texture features, using Gray Level Cooccurrence Matrix (GLCM) which quantifies the spatial relationships between pixel intensities in a range.

9. The system of any one of claim 1 to claim 8, wherein functional information comprises information relating to tissue chromophores, which are isolated by spectral unmixing of the input.

10. The system of any one of claim 1 to claim 9, wherein functional information comprises variations in tissue vascularity and / or metabolic activity.

11. A method for delineating tumour boundaries, comprising: receiving an input of a three-dimensional (3D) image stack comprising a plurality of images; generating a plurality of maximum intensity projections from the 3D image stack; initialising a set of boundary estimates for bounding a region of interest in each image, and iteratively performing level set segmentation to update the boundary estimates until convergence at a final boundary bounding the region of interest; and extracting structural and functional information from the region of interest.

12. The method of claim 11, wherein each of the plurality of maximum intensity projections is generated for a respective anatomical plane.

13. The method of claim 11 or claim 12, comprising preprocessing the 3D image stack, before generating the plurality of maximum intensity projections, by applying a median filter to reduce noise in the 3D image stack.

14. The method of any one of claim 11 to claim 13, wherein the boundary estimates are initialized using a maximum intensity pixel, with each boundary estimate differing from all other boundary estimates.

15. The method of any one of claim 11 to claim 14, wherein level set segmentation is guided by an energy minimisation process.

16. The method of any one of claim 11 to claim 15, wherein level set segmentation stops when the updated boundary estimates converge on a common boundary.

17. The method of any one of claim 11 to claim 16, wherein extracting structural information comprises one or more than one of: extracting width and depth measurements, determining a maximum feret diameter, and extracting a length of major and minor axes of the region of interest.

18. The method of any one of claim 11 to claim 17, wherein extracting structural information comprises analysing texture features, using Gray Level Cooccurrence Matrix (GLCM) which quantifies the spatial relationships between pixel intensities in a range.

19. The method of any one of claim 11 to claim 18, wherein functional information comprises information relating to tissue chromophores, which are isolated by spectral unmixing of the input.

20. The method of any one of claim 11 to claim 19, wherein functional information comprises variations in tissue vascularity and / or metabolic activity.

Images