DETECTION OF CANCER IN ANIMAL TISSUES: A WAVELET APPROACH

: Considering that the biospeckle laser is a dynamic interferometric phenomenon adopted as a tool to monitor changes in biological samples and that the temporal variation of speckle pattern depend on the activity level of the sample surface illuminated, this work proposes to analyse the time-varying scale-mixing matrix. Using two-dimensional scale-mixing wavelet transform several descriptive summaries varying on time are derived. These descriptors are signature of image regularity and fractality useful in tissue classiﬁcation. In this work we propose to verify the behavior of the energy-ﬂux between the scales, considering a set of 128 images to classifying cancer areas in images of an anaplastic mammary carcinoma in a female canine and in images of skin cancer in a cat obtained over time. The time-varying spectral slopes applied in the analysis of dissimilarities of tissues allowed to note that healthy area descriptors have lower values than cancer area descriptors, resulting in higher Hurst exponents. By using scaling properties of tissue images, we have captured information contained in the background tissue of images which is not utilized when only considering traditional morphological analysis.


Introduction
The research of cancer identification demands permanent efforts to scientists. The optical techniques are an actual alternative to achieve the diagnosis of tumors (KURACHI et al., 2008 ) which are also known as optical biopsy. The development of technologies associated with the dynamic laser speckle has offered alternatives to access the sensitive activities in animal and vegetable tissues (BRAGA et al., 2012).
The biospeckle laser is a dynamic interferometric phenomenon, which has been adopted as a sensitive tool to monitor changes in biological samples, and thus it has been applied in many areas, from medicine to agriculture, since the fact that it is a non-destructive technique is relevant in biological applications (BRAGA et al., 2008).
The multitude of applications is associated to a range of methods to illuminate, to assemble the images and to provide their analysis.
The static appearance of a speckle pattern is expressed by an image with clear and dark grains distributed over all the illuminated material. Since it is a dynamic phenomenon, the speckle pattern expresses the boiling effect with the grains changing their shape and level of lightness related to the level of movement of the scatterers of the coherent light. The scatterers in an inert material, for example, do not change, so much so, the boiling effect isn't observed and its intensity is linked to the level of activity. This measurement technique is identified as an alternative to that of a conventional mechanical stylus (DAINTY et al., 1975.) Several speckle pattern analysis methods have been investigated in literature to facilitate different application tasks including the speckle contrast (GEORGE, 1975), autocorrelation function (CHENG, 2004), etc. In addition to speckle pattern analysis, the temporal variations of speckle patterns are found to depend strongly on the activity level of the sample surface illuminated.
The irregular behavior of complex structures is difficult or impossible to quantify by standard modeling techniques; but when observations are inspected at different scales, there is in fact a regular relationship between the behavior among the scales. This phenomenon has been demonstrated in many medical images (JEON et al., 2014), leading to its diagnostic use as a tool capable of quantifying statistical similarity of data patterns at various scales (NICOLIS et al., 2011).
Wavelets are a mathematical tool for extracting information from different types of data. The Fourier analysis is suitable for stationary signals while the wavelet analysis is robust to non-stationary signals, because its transform has a local aspect. One key difference, however, is that the wavelet transform is localized in both frequency and time, while the standard Fourier transform is only localized in frequency. In other words, the Fourier transform tells us what frequencies are present in a signal, and the wavelet transform tells us what frequencies are present and at which locations. Wavelets are used for a variety of purposes, including the compression, denoising, and filtering of signals as well as measuring the degree of regularity of a signal. Often, discussions of wavelets focus on the study of one-dimensional signals, but wavelet techniques may also be applied to the study of multi-dimensional objects (NICOLIS et al., 2011;TONER, 2019). Statistical aspects of wavelets are discussed in Vidakovic (1999) and Morettin (2014).
Wavelet transforms lead to coefficients (numerical values) representing the nature of a given signal at different locations/resolutions. These coefficients may be used to form the wavelet-based spectra of the signal, showing the relationship between the resolution of the signal and the averaged magnitudes of the coefficients. By assessing the wavelet-based spectra, we may better understand the mathematical characteristics of the overall signal. If the energies (an engineering term for squared coefficients in the wavelet decomposition) decay regularly, this signifies scaling in the data, meaning all resolutions contribute to the overall observed phenomenon. In this case, a measure of regularity can be calculated as the rate of energy decay. More precisely, if the logarithms of average energies in different scales decay linearly with the scale index, then the slope of this decay is describing the regularity of the original signal/object. Thus the spectral slope of the wavelet-based spectra can precisely measure the degree of a signal's regularity. For details about wavelet-based spectra and their application to assessing regularity of signals/images we direct the reader to Jeon et al. (2014), Nicolis et al. (2011) andSáfadi et al. (2016).
Traditional 2-D wavelet transforms lead to three directional spectra: horizontal, vertical, and diagonal. These spectra are defined by directional hierarchies of detail coefficient. Marques and Sáfadi (2020) develop new applications for image analysis based on the non-decimated discrete wavelet transform. The behavior of (directional) Hurst exponents were verified, over time, using the nondecimated wavelet transform in 128 interferometric images of a canine anaplastic mammary carcinoma, obtained in regular time intervals by means of the dynamic biospeckle method. They concluded that the temporal analysis of the Hurst exponents is an alternative tool to distinguish between images of the healthy tissue and the carcinogenic one.
An ensemble of wavelet based-spectral tools for analysis of 2-D images was proposed on Ramirez-Cobo et al. (2011) where generalized the notion of 2-D spectra by generalizing the form of the 2-D wavelet transform and showed that the new spectra was capable of interfacing different scales when assessing the energy distribution of the image. Scale-mixing wavelet spectrum is applied to the analysis of time sequences of two-dimensional spatial rainfall radar images characterized by either convective or frontal systems. The scale-mixing (or covariance) matrix captures the "energy flux" between the scales. Jeon et al. (2014) generalized the scale-mixing wavelet spectra to the complex wavelet domain. In this domain, they estimated Hurst parameter and image phase and used them as discriminatory descriptors to classify mammographic images to benign and malignant. Namazi and Kiminezhadmalaie (2015) evaluated the Hurst exponent for 50 healthy and 50 lung cancer patients, the analysis found that DNA walks have smaller values of the Hurst exponent compared to normal DNA. They concluded that DNA damage is less predictable and more complex compared to normal DNA.
In this paper, we propose to find dissimilarities between cancer images and healthy images by analyzing the flow of energy between scales and over time. The methodology is based on a modification of a 2-D discrete scale-mixing wavelet transform proposed by Ramirez-Cobo et al. (2011) to classifying cancer areas in images of an anaplastic mammary carcinoma in a female canine and in images of skin cancer in a cat. We considered a set of 128 interferometric images of cancer tissues (breast and skin), obtained in regular time intervals in a rate of 0.08 s. Braga et al. (2012) showed that with these 128 images it is possible to separate different tissues in the same material through the frequency signature, and through the association of graphic and numerical results of the biospeckle laser images. The empirical spectral slopes densities and the Hurst exponent are used in the analysis of dissimilarities in tissues of distinct areas.
The rest of the article is organized as follows. In Section 2, concepts from 2-D scale-mixing non-decimated wavelet spectra are reviewed. In Section 3, the data and the methodology are described. In section 4 are presented the results. The conclusions are presented in Section 5.

Background
We describe briefly the non-decimated wavelet transform (NDWT), and provide motivation for the choice of this redundant transform over usual orthogonal discrete wavelet transform.
Any square-integrable L 2 (R) function f (x) can be represented in the wavelet domain as where c J0,k indicates coarse coefficients, d j,k detail coefficients, φ J0,k (x) scaling functions, and ψ jk (x) wavelet functions. We use different decomposing atom functions, as scaling and wavelet functions, depending on a version of wavelet transform. For standard discrete wavelet transform (DWT), the atoms are where x ∈ R, j is a resolution level, J 0 is the coarsest resolution level, and k is the location of an atom.
For NDWT, atoms are Notice that atoms in NDWT have a constant location shift k at all levels, which yields the maximal sampling rate at each level. Two types of coefficients, c J0,k and d j,k , capture coarse and detail fluctuations of an input signal, respectively. These are obtained as In a p-depth decomposition of an input signal of size m, a NDWT yields m × (p + 1) wavelet coefficients, while DWT yields m wavelet coefficients independent of p.
The redundant transform NDWT decreases the variance of the scaling estimators.
In wavelet analysis of signals/images the redundancy of non-decimated wavelet decompositions is often favored to minimality of standard orthogonal wavelet transforms because of their stationarity, smaller variance, and ease of use for signals of arbitrary size. Expanding on the 1-D definitions, we overview a 2-D NDWT of f(x,y), where (x, y) ∈ R 2 . Several versions of 2-D NDWT exist but we focus on the scale-mixing version. For the standard 2-D NDWT, the wavelet atoms are where J 0 is the coarsest decomposition level, and i ∈ {h, v, d} indicates the "orientation" of details coefficients as horizontal, vertical and diagonal (VIDAKOVIC, 1999) . As in the univariate case, any function f ∈ L 2 (R 2 ) can be represented as For the scale-mixing 2-D NDWT, the wavelet atoms are where J 01 and J 02 are coarsest levels, j 1 ≥ J 01 ; j 2 ≥ J 02 , and k = (k 1 , k 2 ). As a result, we obtain wavelet coefficients for f (x, y) from the scale-mixing NDWT as The coefficients (Equation (3)) will be denoted as c, h, v, and d-type coefficients. Notice that in the standard NDWT, we use common j to denote a scale, while in the scale-mixing NDWT, we use the pair (j 1 , j 2 ), which indicates that two scales are mixed. The d-coefficients correspond to decomposing atoms consisting of two wavelet functions, while the atoms of c, v or h-type coefficients contain at least one scaling function ( see Figure 1a). When an image possess a certain degree of smoothness, the coefficients corresponding to diagonal decomposition atoms (d-type coefficients) tend to be smaller in magnitude compared to the c, v or h-type coefficients. Figure 1a illustrates the tessellation of the scale-mixing 2-D NDWT (KANG, 2016) .
In practice, we operate with vectors and images as sampled versions of 1D and 2D functions. In this discrete case, the wavelet transforms (orthogonal and nondecimated) due to their linearity can be represented in form of matrix multiplication. Such wavelet-matrices can be constructed prior to implementation, which facilitates computational efficiency in repeated use.
A 2-D signal A of size [m x n] for p 1 and p 2 level decomposition along rows and columns, respectively, is obtained by NDWT matrix multiplication from the left and its transpose from the right. The transform results in a 2-D signal B of size (p 1 + 1)m × (p 2 + 1)n, Here p 1 , p 2 , m and n can take any integer value, and W  Matrix B will be called the scale-mixing (or covariance) wavelet transform of matrix A, and will be the basis for defining the scale-mixing spectra. The motivation for the name scale-mixing is that the atoms mix the scale indices thus capturing the "energy flux" between the scales. From the transformed images several spectral indices, describing innate regularity of image background, are derived. The descriptors involve spectral slopes which are directly connected with degree of image regularity.
The standard measure of regular scaling is the Hurst exponent. This measure can also be connected to the presence of long memory and fractality in signals and images and is viewed as an informative summary. Wavelet transforms are powerful tool in estimating the Hurst exponent and modeling statistical similarity at different scales.
For defining a wavelet spectrum, and subsequently for estimating Hurst exponent H, only detail wavelet coefficients are used. Empirically, we look at the levelwise average of squared detail coefficients (energies), where n j is the number of wavelet coefficients at level j. The relationship between average energy d 2 j and H is The spectral slope −(2H +1) is a function of Hurst exponent and is a signature of irregularity of a signal. For D-dimensional signals, the spectral slope becomes −(2H + D). Operationally, the spectral slopes are found by regressing log 2 d 2 j on the level j.
By inspecting the tessellation in Figure 1b, several details spaces can be identified. Analogous to the one-dimensional case, the scale-mixing spectra are defined in terms of the scale-mixing coefficients (Equation (3)) as where j ∈ Z and s ∈ Z is fixed. The empirical counterpart of Equation (5) iŝ for j ∈ Z.
As in Veitch and Abry (1999) and Nicolis et al. (2011), it is assumed here that the coefficients within and across the scales are uncorrelated. This assumption is reasonable; Flandrin(1992) showed that when the number of vanishing moments of the scaling function is large, the correlation between the coefficients within a scale decays exponentially fast, while the coefficients from different scales are almost uncorrelated  .
It can be seen that E(d 2 j,j+s,k ) = 2 −j(2H+2) V ψ,s (H); where V ψ,s (H) is a constant depending on ψ, H and s, but not on the scale j. By taking logarithms in Equation (7), for j ∈ Z, and thus the Hurst exponent can be estimated from the slope of the linear Equation (8). Finally, the empirical counterpart of (Equation (8)) is a regression defined on (j, log 2 (d 2 j,j+s,k )), j, s ∈ Z.
The slope is indeed the tangent of angle of a regression where log-average energies are regressed against the dyadic index of multiresolution subspaces.
Notice that the case s = 0 in Equation (6) corresponds to the traditional diagonal 2-D spectra (descriptors s 1 and s 2 ), while s = −1 corresponds to descriptor s 3 and s = 1 corresponds to s 4 . If the image were a monofractal (such as a fractional Brownian field) the two slopes (s 1 , s 2 ) would be identical and would lead to a global irregularity index of the image, such as Hurst exponent. The difference between the slopes will measure the deviation from the monofractality.
On the other hand, the two lateral slopes s 3 and s 4 measure directional violation of isotropy. For example, the slope s 3 is calculated from average energies in two multiresolution cells where horizontal and vertical resolutions differ. The average energies in corresponding cells measure the degree of horizontal or vertical anisotropy in the image.
If the image were fully isotropic, both s 3 and s 4 would be equal to s 1 and s 2 , that is, all four indices will be the same.
We draw the reader's attention to the fact that, in this work, the matrix B and the descriptors associated with it are indexed in time (index omitted throughout the text).

Material and Methods
The biological materials chosen to be illuminated were neoplastic tissues from surgeries in a cat (skin) and in a dog (breast) and were analyzed before the histological exams. They were maintained cooled just after the surgeries and they were not fixed in phormaldehyde which is a routine step to implement histological exams. The samples were cut in pieces creating a flat surface with the neoplastic and the normal tissues naturally linked side by side. After the illumination, the neoplastic tissue was fixed in phormaldehyde 10%, and dehydrated through increasing concentrations of ehtyl alcohol, diaphanized in xylol and included in paraffin in order to be analyzed in histological exams. The cuts of 5 µm were pigmented using Hematoxylin-Eosin (HE) staining protocol, and the neoplastic diagnosis was conducted in accordance to histological and to pathological aspects.
Interferometric images of neoplastic tissues were obtained through an optical method using the interference pattern formed when a biological material is illuminated by a laser, called Biospeckle laser technique (BSL). A laser set of HeNe with 10mw, of 632nm, was opened and illuminated the sample directed by a mirror. The CCD camera and the computer were responsible for the assembling of the image with 486×640 pixels for the animal tissues summing 128 frames for neoplastic tissues in a rate of 0.08s.
In Figure 2 it is possible to observe the images from histological sections with the basophilic neoplastic tissues. The neoplastic tissues were histologically classified as basosquamous carcinoma in the feline's skin and anaplastic mammary carcinoma in a female canine. The anaplastic mammary carcinoma, besides its distinct cell origin, presented lower differentiation if compared to the basosquamous carcinoma. In turn, the anaplastic mammary carcinoma was plenty of fibrous stroma while it was rare in the basosquamous, and in addition the anaplastic mammary carcinoma had areas of inflammatory infiltrated dominated by neutrophils (characterized by whitish areas) which was either absent in the basosquamous carcinoma. The anaplastic mammary carcinoma is observed in the Figure 2a in dark situated bellow and the normal tissues in light on the top of the image. In Figure 2b, the basosquamous carcinoma can be seen in dark, and a normal tissue in the right bottom of the image with light (eosinophlic), in a format of a tail. On right are the illuminated images. To implement the 2D-non-decimated wavelet transform, we consider 4 subimages (64 x 64) in each of the 128 images, two from the cancer area and two from the healthy area. The location of healthy and cancerous tissues was based on research by Braga et al. (2012). The proposed methodology is evaluated through the empirical densities of spectral slopes and the Hurst exponent.
We use the Daubechies' db2 basis with 2 vanishing moments to perform the non-decimated wavelet decomposition with 3 levels of details. Selection of this wavelet was motivated by its good localization properties and the continuity of scaling and wavelet functions.

Results
In this section, we will consider time-varying subfigures of healthy and cancerous areas. For each subfigure, scale-mixing matrix wavelet and several spectral slopes are derived. The dissimilarity analysis is performed through the empirical density of the spectral slopes and the Hurst exponent.

Feline tissues
Initially, we considered two sub-images from the feline image (Figure 2b), the first in the healthy area and the second in the cancer area. The empirical densities for each spectral slope are shown in Figure 3. Note that the distributions for each descriptor are similarly in the two areas. In addition, s 2 and s 4 are bimodal with greater variability in healthy areas. The mean and the confidence interval for each descriptors are shown on Table 1, note that the values for s 1 , s 2 , s 3 and s 4 are distinct indicating anisotropy in the image. For the slopes s 1 (traditional) and s 3 (mixture of scales) whose distributions are unimodal, the confidence intervals indicate different descriptors when comparing images of cancer and healthy tissues. The Hurst exponent was obtained using the mode instead of the mean on Equation (8), because we have asymmetric and bimodal distributions. Note that healthy tissues have greater Hurst exponents than cancer tissue, indicating that healthy areas are more irregular than cancer areas. Considering, as proposed, 2 sub-images from the cancer area and 2 from the healthy area of the feline skin image, it is observed (Figure 4) that the empirical densities for the descriptors of the cancer areas are closer than the densities of healthy areas. Indicating that healthy areas are more irregular than cancer areas.

Canine tissues
The analysis for canine tissues was divided into two parts, in the first we tried to choose regions of the healthy area as homogeneous as possible, Figure 5, and in the second we considered sub-images of whitish regions associated with infections, Figure 6.
When considering healthy areas, but with more evident infection (whitish areas), Figure 6, we observed that the empirical densities of the descriptors of these regions are closer, indicating greater regularity. Statistics of the descriptors for the region with cancer and healthy are shown in Table 2. Note that, as in Table 1   healthy area descriptors have lower values than cancer area descriptors, resulting in higher Hurst exponents. It is interesting to note that despite being considered a healthy area because it does not contain cancer, it is an area of great infection and even so the methodology used was efficient.
In Figures 4 and 6, where we have a basosquamous carcinome and an anaplastic carcinome tissues, respectively, we see that the densities for healthy areas had lower values than the densities for cancer areas, mainly observed in s 1 (diagonal) and s 3 (scale-mixing) descriptors.
Note that healthy area descriptors have lower values than cancer area descriptors, resulting in higher Hurst exponents. A similar result was obtained by Marquez and Sáfadi (2020) using directional Hurst exponents.
Our results are in agreement with those found by Namazi and Kiminezhadmalaie (2015), whose results for damage DNA (lung cancer) in humans showed higher values of the Hurst exponent when compared with normal DNA.

Conclusions
The methodology makes no a priori assumptions about the morphology of a cancer, but rather detects it by departure from normal background. By using scalemixing spectra properties of tissue images, we have captured information contained in the background tissue of images which is not utilized when only considering traditional morphological analysis.
The time-varying spectral slopes applied in the analysis of dissimilarities of tissues allowed to note that healthy area descriptors have lower values than cancer area descriptors, resulting in higher Hurst exponents.
The method used in this research can be applied for analysis of other types of cancer. The ability to apply the methodology in the animal tissue can be considered a potential usage in this research field.
We suggest that this type of classifier may be used in conjunction with other methodologies in order to improve the detection. If the tool proves robust on further investigation, it may be useful in outlining cases that require heightened scrutiny or even addition of supplemental screening modalities.
Sensitivity and specificity analysis will consider in future work.