Muhieddine Khodari

Rachid Souktani

Olivier Le Coz

Dina Bedretdinova

Florence Figeac

Adrien Acquistapace

Pierre Francois Lesault

Julie Cognet

Anne Marie Rodriguez
email: anne-marie.rodriguez@inserm.fr.

René Yiou
email: rene.yiou@hmn.aphp.fr

Prof René Yiou

Monitoring of erectile and urethral sphincter dysfunctions in a rat model mimicking radical prostatectomy damage

Keywords: erectile dysfunction, urinary incontinence, radical prostatectomy, retrograde leak point pressure, arterial penile laser Doppler flow, durometry, apomorphine, post prostate cancer impairments erectile function RP: radical prostatectomy SUS: striated urethral sphincter USF: urethral sphincter function

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

Radical prostatectomy (RP), the reference treatment for organ-confined prostate cancer, continues to cause erectile dysfunction and urinary incontinence despite improvements in the surgical technique [START_REF] Glickman | Changes in continence and erectile function between 2 and 4 years after radical prostatectomy[END_REF][START_REF] Coelho | Retropubic, laparoscopic, and robot-assisted radical prostatectomy: a critical review of outcomes reported by high-volume centers[END_REF]. These impairments can severely impact quality of life [START_REF] Litwin | The UCLA Prostate Cancer Index: development, reliability, and validity of a health-related quality of life measure[END_REF][START_REF] Droupy | Quality of life after localized prostate cancer treatment[END_REF].

Post-RP impotence and urinary incontinence constitute a complex bifocal dysfunction resulting from a surgical injury. Post-RP erectile dysfunction results from injury to either the cavernous nerves [START_REF] Walsh | Impotence following radical prostatectomy: insight into etiology and prevention[END_REF] or the accessory pudendal arteries [START_REF] Benoit | Supra and infralevator neurovascular pathways to the penile corpora cavernosa[END_REF]. There is growing evidence that massive apoptosis of corporeal cells after RP plays a key role in the pathophysiology of post-RP erectile dysfunction [START_REF] User | Penile weight and cell subtype specific changes in a post-radical prostatectomy model of erectile dysfunction[END_REF][START_REF] Fall | Apoptosis and Effects of Intracavernous Bone Marrow Cell Injection in a Rat Model of Postprostatectomy Erectile Dysfunction[END_REF]. Urinary incontinence is chiefly related to direct damage to the striated urethral sphincter (SUS) [START_REF] Krane | Urinary incontinence after treatment for localized prostate cancer[END_REF] or its nerve supply [START_REF] Berry | Is there correlation of nerve-sparing status and return to baseline urinary function after robotassisted laparoscopic radical prostatectomy?[END_REF].

At present, there are no animal models replicating post-RP injuries and allowing monitoring of both urethral sphincter function (USF) and erectile function (EF) over time.

This absence of models impedes the development of novel therapeutic strategies.

Recently, cell-based treatment strategies were developed to repair post-RP damage to the penile or sphincter cells. For instance, adipose-derived stem cells were successfully used to treat animal models of either urinary incontinence [START_REF] Jack | Processed lipoaspirate cells for tissue engineering of the lower urinary tract: implications for the treatment of stress urinary incontinence and bladder reconstruction[END_REF] or erectile dysfunction [START_REF] Fandel | Recruitment of Intracavernously Injected Adipose-Derived Stem Cells to the Major Pelvic Ganglion Improves Erectile Function in a Rat Model of Cavernous Nerve Injury[END_REF][START_REF] Albersen | Multipotent stromal cell therapy for cavernous nerve injury-induced erectile dysfunction[END_REF].

However, these models only partially replicated the clinical situation, as most patients with post-RP urinary incontinence also have erectile dysfunction. Therefore, the development of new regenerative treatments for post-RP impairments would benefit from an animal model allowing assessments of both USF and EF.

Another limitation in earlier studies is the invasiveness of the functional methods used to evaluate USF or EF [START_REF] Fall | Apoptosis and Effects of Intracavernous Bone Marrow Cell Injection in a Rat Model of Postprostatectomy Erectile Dysfunction[END_REF][START_REF] Yiou | Restoration of functional motor units in a rat model of sphincter injury by muscle precursor cell autografts[END_REF][START_REF] Chung | Investigative models in erectile dysfunction: a state-of-the-art review of current animal models[END_REF][START_REF] Hr | A mouse model of cavernous nerve injuryinduced erectile dysfunction: functional and morphological characterization of the corpus cavernosum[END_REF]. The repeated use of these methods during longitudinal studies may induce additional damage and cause histopathological artifacts. Thus, the existing animal models are not suitable for longitudinal studies requiring multiple assessments in each animal.

Our objective here was to develop an improved animal model for investigating treatments for post-RP impairments such as cell-based strategies. We developed a rat model replicating the sphincter and erectile dysfunctions seen after RP and non-invasive methods allowing repeated measures of EF and USF in the same animal.

MATERIALS AND METHODS

Animals

All experiments were performed on 70 adult 6-week-old male Sprague Dawley rats according to institutional guidelines for animal care. The research protocol was approved by the local ethics committee (agreement number A-94028-245).

All surgeries were performed under general anesthesia with isoflurane gas (Baxter, France) inhalation (2.5% in 500 mL/min of air).

Experimental design

Preliminary experiments

In a first set of experiments (n=10), we sought to develop a non-invasive method for assessing USF. The reference method in rats is measurement of the antegrade leak point pressure (LPPa) [START_REF] Yiou | Restoration of functional motor units in a rat model of sphincter injury by muscle precursor cell autografts[END_REF][START_REF] Rodriguez | New objective measures to quantify stress urinary incontinence in a novel durable animal model of intrinsic sphincter deficiency[END_REF][START_REF] Hijaz | Animal models of female stress urinary incontinence[END_REF] defined as the lowest bladder pressure at which urine leakage occurs through the urethral meatus in the absence of detrusor contraction. LPPa corresponds to the maximal bladder capacity. In rats, LPPa measurement is an invasive procedure requiring surgical catheter insertion. Retrograde LLP (LPPr) has been measured in humans after RP [START_REF] Comiter | Retrograde leak point pressure for evaluating postradical prostatectomy incontinence[END_REF] and in female rats [START_REF] Rodriguez | New objective measures to quantify stress urinary incontinence in a novel durable animal model of intrinsic sphincter deficiency[END_REF] as the intraurethral pressure required to open the urethral sphincter during retrograde infusion of the distal urethra. This non-invasive method has not been investigated in male rats. In preliminary experiments, we measured LPPr in male rats as described below and we compared the LPPa and LPPr profiles. The aim was to determine whether LPPr could be used alone to study USF.

Main experiments

Next, we investigated EF and USF after an injury replicating the effects of RP. We randomly allocated 60 rats to either a sham procedure (n=30) or bilateral electrocautery of both sides of the SUS (n=30). USF and EF were evaluated at baseline and on postoperative days 7, 15, 30, 60, and 90 using LPPr measurement.

Rats in the sham-operated and injury group were sacrificed on either day 7 (n=10 including 5 sham operated and 5 injured rats) or day 30 (n=10) post-injury, after USF and EF assessments. The penises and urethras were harvested from each animal for histological analysis. Electrocautery-induced cell damage was evaluated by assessing apoptosis and fibrosis at both time points. All remaining rats were sacrificed on day 90 after the last functional evaluation.

Surgical procedures

A lower-abdominal midline incision was made and the SUS was approached under an operating microscope (Figure 1A). In the sham group, the incision was then closed. In the injury group, an electrocautery was applied using a surgical electrical blade connected to a electrosurgical unit (Erbotom T71, Erbe), for 1 second (50 W, 600 Hz, intensity 7) at each of four points along the SUS (Figure 1B). The cavernous nerves emerge from the lower part of the pelvic plexus and run in contact with the lateral side of the SUS, accompanied with penile blood vessels. Each electrocautery point was also applied on the penile neurovascular supply. Thus, the experimental injury caused bilateral damage to the SUS and penile neurovascular supply. The abdominal incision was then sutured.

Assessment of USF

In preliminary experiments, we investigated LPPa and LPPr in 10 male rats. The mean of four LPPa and LPPr measurements was computed for each animal as follows. A polyethylene catheter (PE20 Phymep, France) was inserted 0.5 cm into the distal urethra via the urethral meatus. The catheter was connected via a Statham pressure transducer amplifier [START_REF] Pogrund | Improved connectors for the Statham blood pressure transducer[END_REF] to a data acquisition board (Notochord Systems, France) and to a programmable infusion pump (Harvard Apparatus, Holliston, MA, USA). A 4-0 polyglactin suture was used to produce a watertight closure of the urethral meatus around the catheter. This catheter served to measure urethral pressure and retrograde bladder filling. The bladder was filled until it became palpable. A second supra-pubic polyethylene catheter (PE50) was then inserted transcutaneously into the bladder (Figure 1C) to measure bladder volume and pressure. The bladder was emptied through the supra-pubic catheter, and saline was infused at a constant rate. LPPa was defined as the bladder pressure when the first drop of urine appeared at the urethral meatus.

Next, we studied the pressure profile corresponding to urethral sphincter opening during retrograde infusion. The bladder was filled to half its maximum capacity using the bladder catheter. Laparotomy was performed to expose the bladder. Saline mixed with methylene blue was infused through the urethral catheter, and urethral pressure was recorded continuously. LPPr was defined as the urethral pressure when methylene blue appeared in the bladder (Figure 1D). After the first LPPr measurement, the bladder was emptied and three further measurements were performed. In each animal, we compared the profiles of four successive LPPr and LPPa values using a filling flow of 200 µL/min for each measurement.

Assessment of EF

To assess EF, apomorphine (Sigma Aldrich, France) was injected subcutaneously (0.05 mg/Kg). Apomorphine is a dopamine agonist, which can create a sexual impulse similar to psychological erection by activating the dopamine receptors in the central nervous system. Apomorphine-induced erection is primarily mediated via the sacral parasympathetic nerve system and the thoracolumbar sympathetic pathway [START_REF] Paick | The neural mechanism of apomorphine-induced erection: an experimental study by comparison with electrostimulation-induced erection in the rat model[END_REF].

Pharmacologically induced erections were monitored using two methods, laser Doppler and durometry, as follows.

Laser Doppler was performed to investigate erythrocyte movements due to arterial flow in the cavernous sinus (Figure 1 E), as previously described [START_REF] Mills | Androgenic maintenance of inflow and venoocclusion during erection in the rat[END_REF]. Briefly, the laser Doppler surface probe (Probe 411, Perimed, France) was maintained at the base of the penis after pulling back the foreskin using a micromanipulator. The probe is designed to assess erythrocyte movement at a depth of 1 mm. Concomitantly, a second laser Doppler probe (Probe 404, Perimed) was placed on the tail of the rat to monitor the effect of apomorphine injection on systemic arterial flow. Both probes were connected to the same transducer amplifier and data acquisition board (Periflux System 5000, Perimed). Results were reported as Perfusion Units.

Apomorphine-induced erections were also investigated using the Shore hardness measurement method consisting in measuring resistance to penetration of a tip. We used a durometer with a needle-shaped indenter (Digital Durometer Shore A, Mitutoyo America Corporation, Chicago, IL, USA) previously used for measuring the hardness of living tissues such as the skin [START_REF] Kissin | Durometry for the assessment of skin disease in systemic sclerosis[END_REF][START_REF] Falanga | Use of a durometer to assess skin hardness[END_REF][START_REF] Leblanc | Durometer measurements of skin induration in venous disease[END_REF] or penile prostheses [START_REF] Wang | Hardness evaluation of penile prostheses[END_REF] (Figure 1F). The force applied can range from 0 to 100 on the Shore scale indicated on the device, i.e., from 55 g to 855 g. When the indenter is pressed to the surface of the specimen, the balance between the resistance force of the specimen (reflecting hardness) and the force of a spring in the tester governs the depth of the indentation.

Laser Doppler and durometry were performed before and after apomorphine injection and the results were expressed as the variations in laser Doppler flow percentage and durometry, respectively.

Histological analysis

Penis and sphincter tissue samples were included in Tissue Tek (Thermo Scientific, France) and snap-frozen in isopentane precooled with liquid nitrogen.

To visualize the extent of fibrosis, serial cryo-sections (8 µm in thickness) of sphincters and penises were stained with Masson's Trichrome to identify collagen in blue, as previously reported [START_REF] Phull | Vulnerability of continence structures to injury by simulated childbirth[END_REF][START_REF] Valente | Larginine and phosphodiesterase (PDE) inhibitors counteract fibrosis in the Peyronie's fibrotic plaque and related fibroblast cultures[END_REF]. Apoptotic cells were detected in sphincters and penises using the TUNEL assay (Roche, France), which ensures enzymatic labeling of DNA strand breaks generated during apoptosis.

Statistical Analysis

GraphPad Prism (Version 5.00, USA) was used to perform the statistical analyses. The Kolmogorov-Smirnov test (with the Dallal-Wilkinson-Lilliefor P value) and Shapiro-Wilk normality test were performed to check the distribution of values. Values of LPPr, LPPa, and unit Shore A durometer hardness were normally distributed and were therefore described as mean±SD. Values of penile laser Doppler flow percentage variation were not normally distributed and were therefore described as the median with the 95% confidence interval (95% CI). P values less than 0.05 were considered statistically significant. Overall betweengroup comparisons were performed using analysis of variance (two-way ANOVA) with repeated measures by columns.

RESULTS

First set of experiments

Filling the bladder through the suprapubic catheter increased bladder pressure to a plateau corresponding to the maximal bladder capacity (Figure 2) and coinciding with urine leakage at the urethral meatus (LPPa). Retrograde bladder filling resulted in a rapid increase in the intraurethral pressure to a plateau phase that consistently coincided with the appearance of methylene blue in the bladder and was therefore taken as the LPPr.

Using a flow perfusion rate of 200 µL/min, we found similar values of LPPa (78±7.1 mmHg) and LPPr (76±8.9 mmHg) and similar profiles in all 10 rats (Figure 2). These findings validated the use of LPPr alone in the main experiments to assess the effect of electrocautery on USF.

Second set of experiments

Electrocautery caused alterations in both EF and USF that persisted throughout the 3 months of the study. 3 and4)

-Erectile function (Figures

In all the rats tested, apomorphine effects became detectable 3 minutes after the injection and lasted 30 minutes (Figure 3). Median variation in baseline arterial penile laser Doppler flow induced by apomorphine injection was +180% (95%CI, 156%-202%) overall, with no significant difference between sham rats (+174%; 95%CI, 149%-217%) and injured rats (+192%; 95%CI, 139%-209%; P>0.05). Apomorphine injection increased the mean baseline durometry value from 3.49±0.5 to 7.16±0.82 Shore A units overall, with no significant difference between sham rats (3.71±0.31) and injured rats (3.63± 0.27, P>0.05).

After electrocautery injury, significant declines occurred in the apomorphine-induced variations in arterial penile flow (p<0.05) and durometry values (p<0.01) and persisted throughout the 3-month study period (figure 4 and table 1). At the end of the study, no significant recovery of EF parameters was noted.

Peripheral arterial flow recorded at the surface of the tail (systemic arterial flow) was not significantly modified by apomorphine injection in any of the groups or at any of the time points (p>0.05).

-Urethral sphincter function

Mean baseline LPPr was 76.8±6.18 mmHg overall, with no significant difference between the sham group (76.87±6.89) and the injured group (76.80±5.6, P>0.05).

Electrocautery injury significantly decreased LPPr, by nearly 50% (Figure 5 and table 1). LPPr values remained unchanged in the sham group, with no significant differences across days 7, 15, 30, 60, and 90 (P>0.05).

-Histological analysis (Figure 6)

Electrocautery injury caused lasting loss of sphincter myofibers. Compared to sham animals, animals 30 days post-injury exhibited significant fibrosis replacing the sphincter myofibers, as assessed by Masson trichrome staining (Figures 6A and6B). TUNEL assays showed numerous apoptotic cells throughout the corpus cavernosum 7 days after electrocautery; no apoptotic cells were found in the sham animals (Figures 6C-6D). In addition, penile fibrosis was seen 30 days post-injury (Figures 6E-6F).

TUNEL assays were also performed on sphincter sections 7 and 30 days post-injury.

However, as nearly all the myofibers and other cells were destroyed immediately after electrocautery, very few apoptotic bodies were detectable 7 and 30 days later (data not shown).

DISCUSSION

We developed a male rat injury model replicating the two main functional impairments seen after RP, i.e., erectile dysfunction and urethral sphincter insufficiency. We also described non-invasive methods allowing repeated measurements of both EF and USF over time in the same animal.

Electrocautery of the lateral aspects of the SUS resulted in lasting USF impairment with LPPr decreases that persisted up to 3 months post-injury. Histological studies showed loss of sphincter myofibers, fibrosis, and no evidence of regeneration. We previously showed that SUS electrocautery in male rats also destroyed the sphincteric nerves [START_REF] Yiou | Restoration of functional motor units in a rat model of sphincter injury by muscle precursor cell autografts[END_REF][START_REF] Zini | The striated urethral sphincter of the pig shows morphological and functional characteristics essential for the evaluation of treatments for sphincter insufficiency[END_REF].

Several hypotheses have been put forward to explain post-RP stress urinary incontinence. including i) direct damage to the urinary sphincter [START_REF] Sasaki | Upfront transection and subsequent ligation of the dorsal vein complex during laparoscopic radical prostatectomy[END_REF][START_REF] Gerullis | Sphincter lesions after radical prostatectomy-evaluation and classification[END_REF]; ii) damage to its innervation [START_REF] Hansen | A neurophysiological study of patients undergoing radical prostatectomy[END_REF]; and iii) a decrease in functional urethral length [START_REF] Paparel | Recovery of urinary continence after radical prostatectomy: association with urethral length and urethral fibrosis measured by preoperative and postoperative endorectal magnetic resonance imaging[END_REF]. Direct damage to the sphincter has been hypothesized to result from accidental electrocoagulation during hemostasis of the dorsal vein complex [START_REF] Sasaki | Upfront transection and subsequent ligation of the dorsal vein complex during laparoscopic radical prostatectomy[END_REF] or sphincter transection [START_REF] Gerullis | Sphincter lesions after radical prostatectomy-evaluation and classification[END_REF]. To our knowledge, no prior studies have documented the exact nature of histopathological changes in the SUS after RP. MRI studies in patients with post-RP stress urinary incontinence showed evidence of periurethral fibrosis that probably involved the SUS [START_REF] Paparel | Recovery of urinary continence after radical prostatectomy: association with urethral length and urethral fibrosis measured by preoperative and postoperative endorectal magnetic resonance imaging[END_REF]. In females with intrinsic sphincter insufficiency, histopathological studies also revealed fibrosis associated with loss of sphincter myofibers [START_REF] Hale | Histologic analysis of needle biopsy of urethral sphincter from women with normal and stress incontinence with comparison of electromyographic findings[END_REF], as observed in our model.

Electrocautery also caused dramatic and persistent decreases in penile responses to apomorphine. Accidental electrocautery of the penile neurovascular bundles is a plausible cause of erectile dysfunction following RP. A previous study in dogs showed that the use of hemostatic energy sources near the prostate during neurovascular bundle dissection was associated with a significantly decreased erectile response to cavernous nerve stimulation [START_REF] Ong | Nerve sparing radical prostatectomy: effects of hemostatic energy sources on the recovery of cavernous nerve function in a canine model[END_REF]. Consequently, most urologists recommend avoiding electrocautery when a nervesparing RP is planned and dissect the penile neurovascular bundles using scissors and hemoclips [START_REF] Salomon | Radical prostatectomy by the retropubic, perineal and laparoscopic approach: 12 years of experience in one center[END_REF]. In keeping with other reports on cavernous-nerve-injury animal models [START_REF] User | Penile weight and cell subtype specific changes in a post-radical prostatectomy model of erectile dysfunction[END_REF][START_REF] Fall | Apoptosis and Effects of Intracavernous Bone Marrow Cell Injection in a Rat Model of Postprostatectomy Erectile Dysfunction[END_REF][START_REF] Hr | A mouse model of cavernous nerve injuryinduced erectile dysfunction: functional and morphological characterization of the corpus cavernosum[END_REF][START_REF] Klein | Apoptosis in the rat penis after penile denervation[END_REF], we found evidence of apoptosis and fibrosis in the cavernous tissue compared to sham animals. In a previous study [START_REF] Ong | Nerve sparing radical prostatectomy: effects of hemostatic energy sources on the recovery of cavernous nerve function in a canine model[END_REF], we showed that bilateral cavernous nerve ablation caused lasting EF impairment, although evidence of spontaneous nerve regeneration was noted, suggesting the presence of compensatory nerve pathways.

Our injury model is original because, in addition to direct damage to the sphincter, it causes indirect injury to the penis through destruction of its neurovascular supply traveling along the sides of the SUS. Electrocautery of the sphincter may not replicate the full spectrum of injuries that can occur during RP but produces some of the histopathological changes previously reported in the penis and sphincter after RP.

Several animal models have been used to replicate either urethral sphincter or erectile dysfunction. They include selective section [START_REF] User | Penile weight and cell subtype specific changes in a post-radical prostatectomy model of erectile dysfunction[END_REF], crushing [START_REF] Fandel | Recruitment of Intracavernously Injected Adipose-Derived Stem Cells to the Major Pelvic Ganglion Improves Erectile Function in a Rat Model of Cavernous Nerve Injury[END_REF][START_REF] Hr | A mouse model of cavernous nerve injuryinduced erectile dysfunction: functional and morphological characterization of the corpus cavernosum[END_REF][START_REF] Yamashita | Inhibition of interleukin-6 attenuates erectile dysfunction in a rat model of nerve-sparing radical prostatectomy[END_REF], or freezing [START_REF] El-Sakka | Effect of cavernous nerve freezing on protein and gene expression of nitric oxide synthase in the rat penis and pelvic ganglia[END_REF] of the cavernous nerves; and sphincter injury by pudendal nerve section [START_REF] Kwon | Neurologic recovery and improved detrusor contractility using muscle-derived cells in rat model of unilateral pelvic nerve transection[END_REF], sphincter laceration [START_REF] Praud | Improvement of urethral sphincter deficiency in female rats following autologous skeletal muscle myoblasts grafting[END_REF], electrocautery, [START_REF] Yiou | Restoration of functional motor units in a rat model of sphincter injury by muscle precursor cell autografts[END_REF] and vaginal distension [START_REF] Hijaz | Animal models of female stress urinary incontinence[END_REF][START_REF] Phull | Vulnerability of continence structures to injury by simulated childbirth[END_REF]. It can be assumed that some of these sphincter injuries -most of which were studied in females-would also induce erectile dysfunction if used in males. However, to our knowledge, combined evaluation of both EF and USF in the same animal has not been reported previously.

In addition to a new model of post-RP impairments, we report evaluation methods whose non-invasive nature is an advantage over previously used investigations. So far, the most widely used method for monitoring EF in rats has been measurement of intracavernosal pressure in response to cavernous nerve electrical stimulation [START_REF] Fall | Apoptosis and Effects of Intracavernous Bone Marrow Cell Injection in a Rat Model of Postprostatectomy Erectile Dysfunction[END_REF][START_REF] Fandel | Recruitment of Intracavernously Injected Adipose-Derived Stem Cells to the Major Pelvic Ganglion Improves Erectile Function in a Rat Model of Cavernous Nerve Injury[END_REF][START_REF] Bochinski | The effect of neural embryonic stem cell therapy in a rat model of cavernosal nerve injury[END_REF]. This is an invasive procedure that is not appropriate for longitudinal studies requiring repeated measurements over time in the same animals. Measurement of intracavernosal pressure requires a surgical approach to stimulate the cavernous nerves, as well as denudation of the penile skin to insert a pressure catheter. Repeated intracavernosal pressure measurements in the same animal are likely to damage the cavernous tissue and therefore to bias the results of subsequent histological or molecular assays. In addition, according to our experience, multiple laparotomies result in adhesions that complicate the approach to the cavernous nerves needed to perform electrostimulation.

The methods used for determining erectile function in this study have not been validated with respect to conventional measurements of erectile responses (intracavernosal pressure). However, laser Doppler and durometry studies provided reproducible values that were dramatically altered after injury. In the clinical setting, Doppler sonography has provided invaluable insights into the pathophysiology of post-RP erectile dysfunction.

Mulhall et al. [START_REF] Mulhall | The hemodynamics of erectile dysfunction following nerve-sparing radical retropubic prostatectomy[END_REF][START_REF] Mulhall | Erectile dysfunction after radical prostatectomy: hemodynamic profiles and their correlation with the recovery of erectile function[END_REF] demonstrated a vascular mechanism involving arterial insufficiency or veno-occlusive dysfunction in the genesis of post-RP erectile dysfunction. Apoptosis of smooth muscle cells of the subalbugineal veins was also suggested in rat models of cavernous nerve injury [START_REF] User | Penile weight and cell subtype specific changes in a post-radical prostatectomy model of erectile dysfunction[END_REF][START_REF] Fall | Apoptosis and Effects of Intracavernous Bone Marrow Cell Injection in a Rat Model of Postprostatectomy Erectile Dysfunction[END_REF]. Arterial insufficiency can be caused by injury to the accessory pudendal arteries and to the vessels accompanying the cavernous nerves [START_REF] Droupy | Assessment of the functional role of accessory pudendal arteries in erection by transrectal color Doppler ultrasound[END_REF]. In our model we clearly identified and severed the vascular pedicles traveling with the cavernous nerves toward the corpus cavernous. Post-RP erectile dysfunction can be regarded as a complex vascular disease. Consequently, the development of methods for assessing penile vascularization in animal models of post-RP erectile dysfunction is essential to obtain further pathophysiological information and to evaluate novel therapeutic strategies. Our data do not allow us to determine whether the laser Doppler signal impairment in response to apomorphine preponderantly reflected arterial penile insufficiency, veno-occlusive dysfunction, or loss of neurological command. Further laser Doppler studies are ongoing in our laboratory to identify the exact mechanism involved.

Methods for investigating USF rely chiefly on LPPa measurement. A major challenge with these methods is standardization of the procedure used to increase the bladder pressure up to the point when urine leakage occurs. This pressure increase has been achieved using progressive manual pressure to the abdomen [START_REF] Brune | Comparison of alpha 1-adrenoceptor agonists in canine urethral pressure profilometry and abdominal leak point pressure models[END_REF], placement of the animal on a tilt table with a pressure clamp [START_REF] Chermansky | A model of intrinsic sphincteric deficiency in the rat: electrocauterization[END_REF], sneezing induced by pepper [START_REF] Lin | Effect of simulated birth trauma on the urinary continence mechanism in the rat[END_REF], and constant bladder infusion through a supra-pubic catheter [START_REF] Yiou | Restoration of functional motor units in a rat model of sphincter injury by muscle precursor cell autografts[END_REF][START_REF] Rodriguez | New objective measures to quantify stress urinary incontinence in a novel durable animal model of intrinsic sphincter deficiency[END_REF]. Such methods require evaluation of bladder filling before leak-point pressure measurement. In male rats, urethral catheterization up to the bladder is difficult to perform. Consequently, surgery is needed to insert the perfusion and pressure catheters into the bladder. Contrary to LPPa assessment, LPPr assessment does not require verification of bladder filling. In keeping with previous reports in female rats, LPPr has been found reproducible for assessing USF, with a strong correlation to LPPa [START_REF] Rodriguez | New objective measures to quantify stress urinary incontinence in a novel durable animal model of intrinsic sphincter deficiency[END_REF].

This animal model holds potential as a tool for developing new therapeutic strategies designed to correct post-RP impairments, including cell therapy strategies. We are currently testing various sources of progenitor cells to repair the damage observed in this model.

Figure 1 :Figure 2 :

 12 Figure 1: Surgical injury to the striated urethral sphincter and penile neurovascular

Figure 3 :

 3 Figure 3: Example of the laser Doppler response to apomorphine injection

Figure 4 :

 4 Figure 4: Results of erectile function assessments using laser Doppler and durometry

Figure 5 :

 5 Figure 5: Results of retrograde leak point pressure (LPPr) measurement to assess

Figure 6 :

 6 Figure 6: Histopathological changes in the penis and striated urethral sphincter after

Table legend Table 1 :

 legend1 Values of penile laser Doppler flow percentage variation, Shore A durometer study and LPPr values corresponding to figures 4 and 5.

Conclusion:

Electrocautery of the striated urethral sphincter caused severe and lasting impairment of EF and USF that could be monitored repeatedly using minimally invasive methods. This new animal model may hold potential for developing new treatments designed to correct post-RP impairments.

CONCLUSIONS

We developed a rat model that replicates the combined urological impairments seen in men after RP. The non-invasiveness of the evaluation methods allows repeated measurements in each animal. The absence of spontaneous regeneration after electrocautery injury is a key factor for long-term studies. Taken together, these characteristics should prove valuable for the development of future therapeutic strategies aimed at restoring both EF and USF after RP.


J.-F Bercher

Escort entropies and divergences and related canonical distribution

Keywords: Generalized Rényi and Tsallis entropies, q-gaussian distributions, Escort distributions PACS: 02.50.-r, 05.90.+m, 89.70.+c

We discuss two families of two-parameter entropies and divergences, derived from the standard Rényi and Tsallis entropies and divergences. These divergences and entropies are found as divergences or entropies of escort distributions. Exploiting the nonnegativity of the divergences, we derive the expression of the canonical distribution associated to the new entropies and a observable given as an escort-mean value. We show that this canonical distribution extends, and smoothly connects, the results obtained in nonextensive thermodynamics for the standard and generalized mean value constraints.

Introduction

Rényi and Tsallis entropies extend the standard Shannon-Boltzmann entropy, enabling to build generalized thermostatistics, that include the standard one as a special case. The thermodynamics derived from Tsallis entropy, the nonextensive thermodynamics, has received a high attention and there is a wide variety of applications where experiments, numerical results and analytical derivations fairly agree with the new formalisms [START_REF] Tsallis | Introduction to Nonextensive Statistical Mechanics[END_REF]. Some physical applications of the generalized entropies, -including statistics of cosmic rays, defect turbulence, optical lattices, systems with long-range interactions, superstatistical systems, etc., can be found in the recent review [START_REF] Beck | [END_REF] and references therein. In the extended thermodynamics, it has been found particularly useful to use generalized moments [3,4]. These moments are computed with respect to a deformed version of the density at hands, which is called its escort distribution [5,[START_REF] Beck | Thermodynamics of Chaotic Systems[END_REF]. Actually, several type of constraints [START_REF] Tsallis | [END_REF]8,9] have been used in order to derive the canonical distributions: a first type of constraints is expressed as standard linear mean values while a second type of constraints is given as generalized escort mean values. In the two cases, the related canonical distributions are expressed as a q-Gaussian distribution, but with two opposite exponents, and both solutions reduce to a standard Gaussian distribution in the q = 1 case. These q-Gaussian distributions can exhibit a power-law behavior, with a remarkable agreement with experimental data, see for instance [START_REF] Tsallis | Introduction to Nonextensive Statistical Mechanics[END_REF]10,11], and references therein. These distributions are also analytical solutions of actual physical problems, e.g. [12,13,14]. There has been numerous discussions regarding the choice of a 'correct' form of the constraints, either as a standard mean or an escort-average, and on the connections between the solutions and associated thermodynamics. In particular, dualities and equivalences between the two settings have been described. This is precisely the context of the present Letter, where we propose a simple connection between these different formulations and between the related canonical distributions. More specifically, we suggest a simple way to combine the originally distinct concepts of entropy and escort distributions into a single two-parameter (a, λ)entropy. Then, we propose to look at an associated extended maximum entropy problem. This approach includes the two aforementioned formulations as special cases. Exploiting the nonnegativity of the associated divergence, we derive the expression of the canonical distribution for an observable given as an escort-mean value. We show Email address: jf.bercher@esiee.fr (J.-F. Bercher)

Preprint submitted to Elsevier September 15, 2011 that this canonical distribution extends, and smoothly connects, the results obtained in nonextensive thermodynamics for the standard and generalized mean value constraints. We begin by recalling the context and definitions in section 2. Then, the combined escort divergences and entropies are introduced and discussed in section 3. The related maximum entropy problem and its solution are described in section 4. Finally, in section 5, we illustrate the results in the case of a two-level system.

Context and definitions

Main definitions

Let us recall that if f and g are two probability densities defined with respect to a common measure µ, then for a parameter q > 0, called the entropic index, the Tsallis information divergence is defined by

D (T ) q (f ||g) = 1 q -1 f (x) q g(x) 1-q dµ(x) -1 ,
and similarly, the Rényi information divergence is defined by

D (R) q (f ||g) = 1 q -1 log f (x) q g(x) 1-q dµ(x), (1)
provided, in both cases, that the integral is finite. By l'Hospital's rule, both Tsallis and Rényi information divergences reduce to the Kullback-Leibler information in the case q = 1. Associated to these divergences, we obtain the entropies when g(x) is uniform with respect to µ. For q = 1,

S q [f ] = 1 1 -q f (x) q dµ(x) -1 , (2)
is the Tsallis entropy, and

H q [f ] = 1 1 -q log f (x) q dµ(x), (3)
is the Rényi entropy. In both cases, we have

S 1 [f ] = H 1 [f ] = -f log f dµ(x)
, the Shannon entropy, for q = 1.

For the Lebesgue measure, the definitions reduces to the standard differential entropies. Finally, in the discrete case, the continuous sum is replaced by a discrete one which extends on a subset D of integers, and µ is a counting measure, usually taken as uniform. We will also denote

M q [f ] = f (x) q dµ(x) (4)
the integral of the density raised to the power q.

The escort distributions used in nonextensive statistics are defined as follows. If f (x) is an univariate probability density, then its escort f a (x) of order a, a ≥ 0, is given by

f a (x) = f (x) a f (x) a dµ(x) , (5)
provided that M a [f ] = f (x) a dµ(x) is finite. Given that f a (x) is the escort of f (x), we see that f (x) is itself the escort of order 1/a of f a (x). Accordingly, the (absolute) generalized a-moment of order p is defined by

m p,a [f ] = |x| p f a (x)dx = |x| p f (x) a dµ(x) f (x) a dµ(x) . (6)
Of course, standard moments are recovered in the case a = 1.

In the context of nonextensive statistical mechanics, the parameter a is taken equal to the parameter q of the q-entropy. We shall also point out that the 'deformed' information measure like the Rényi entropy (3) and the escort distribution (5) are originally two distinct concepts, as indicated here by the different notations a and q.

There is an interesting discussion on this point in [15]. It is a contribution of the present paper to indicate that it is still possible to adopt two different values for these parameters.

The maximum entropy problems

In nonextensive thermostatistics, the maximum q-entropy distributions, which are identified as q-gaussian distributions, can be obtained using a constraint expressed as a classical mean value; but it has also been found adequate to use constraints given as escort-mean values, with the same index q as in the entropy. This is discussed in [START_REF] Tsallis | [END_REF] or in [3,4]. In the first case, the maximum entropy problem stands as the variational problem

H (1) q (m) = max f {H q [f ] : m p,1 [f ] = m and m 0,1 [f ] = 1} (7)
which consists in finding a distribution with maximum Rényi entropy on the set of all probability distributions with a fixed moment of order p: m p,1 [f ] = m. The value of the maximum entropy obtained for a given moment is denoted H q (m) -the use of the square brackets and parenthesis distinguishes between the functions of the state and the functions of the observable. Similarly, for a constraint given as an escort-mean value, we have

H (q) q (m) = max f {H q [f ] : m p,q [f ] = m and m 0,1 [f ] = 1} . ( 8
)
Note that although the problems above have been written for the Rényi entropy, the very same problems can be written for Tsallis entropy; actually the optimum distributions are identical and the optimum values of both maximum entropy problems are easily related.

In the case of the standard mean constraint, the density that achieves the maximum entropy is given by the generalized q-Gaussian f

(x) = 1 Z(β) (1 -(q -1) β|x| p ) 1 q-1 + (9)
while in the case of the escort mean constraint, the solution has an exponent with opposite sign

f (x) = 1 Z(β) (1 -(1 -q) β|x| p ) 1 1-q + . ( 10
)
In both cases, we used the notation (x) + = max {x, 0} , and β is a positive parameter chosen so as to satisfy the constraint. The limit case q → 1 gives f (x) ∝ exp (-β|x| p ) . The partition functions Z(β) used to normalize f (x) as a density will be given below.

Dualities

There has been a great deal of discussions concerning the relationships between the two problems ( 7) and ( 8) with standard and escort mean value constraints, and between their solutions (9) and (10). Dualities between the solutions have been introduced: the q → 2 -q transformation [16,17] enables to obtain the density (10) from the maximization of the entropy S 2-q [f ] subject to standard constraints rather than the maximization of S q [f ] subject to escort generalized constraints. The q → 1/q duality rests on the duality between the original density and its escort: as already noted, if f q is the escort of order q of f , then f is itself the escort of order 1/q of f q . This has been mentioned in [START_REF] Tsallis | [END_REF], the related duality discussed in [18] and in [19]. It is very easy to obtain the q → 1/q duality as the consequence of the symmetry above for escort distributions: it suffices to note that for the Rényi entropy, we always have the equality

H 1 q [f q ] = H q [f ].
Therefore, noting that we also have m p,q [f ] = m p,1 [f q ], we see that the problem

H (q) q (m) = max f {H q [f ] : m p,q [f ] = m and m 0,1 [f ] = 1} , (11)
is exactly equivalent to

H (q) q (m) = max f H 1 q [f q ] : m p,1 [f q ] = m and m 0,1 [f ] = 1 = H (1) 1 q (m) (12)
where the maximum entropy is obtained for the density (10). Therefore, it is always possible to swap the standard and generalized constraints, provided that entropic index q is changed into 1/q. We shall also mention that these dualities have been linked to combinatorial considerations in [20].

The (a, λ) divergences and entropies

In complement to the dualities mentioned above, we will show that there is a continuum of q-Gaussians, solutions of an extended maximum entropy problem, that smoothly connects the problems [START_REF] Tsallis | [END_REF] and (8) and their solutions (9) and (10). The basic idea is to mix the concepts of q-entropies and escort distributions into a single quantity. This leads us to a simple extension of the Rényi (Tsallis) information divergence and entropy to a two parameters case. Interestingly, the generalized (a, λ)-Rényi information divergence that emerges in (14) has been mentioned [21, eq. ( 49)] as a possible extension of a new class of (α, β)-divergences, generalizing a family of Gamma-divergences recently introduced in [22]. Definition 1. Given two parameters a and λ > 0 and two densities f and g, we call (a, λ)-Rényi information divergence the standard Rényi information divergence with index q = a/λ between the escort distributions like (5), of order λ, f λ and g λ , associated to f and g :

D (R) a,λ (f ||g) = D (R) a λ (f λ ||g λ ). (13)
In the developed form, this gives

D (R) a,λ (f ||g) = 1 a -λ log f (x) a g(x) λ-a dµ(x) λ f (x) λ dµ(x) a g(x) λ dµ(x) λ-a , (14)
provided that the different integrals exist. Similarly, we call (a, λ)-Tsallis information divergence the quantity

D (T ) a,λ (f ||g) = 1 a -λ f (x) a g(x) λ-a dµ(x) λ f (x) λ dµ(x) a g(x) λ dµ(x) λ-a -1 (15)
provided that the different integrals exist.

These definitions include the Rényi (Tsallis) information divergence D (.)

a (f ||g) in the case λ = 1. In the case λ = a, again by l'Hospital's rule, (14) becomes D (.) λ,λ (f ||g) = D(f λ ||g λ ), the Kullback-Leibler divergence between the escort distributions of order λ, f λ and g λ , associated to f and g. Finally, in the case a = λ = 1, the generalized Rényi (Tsallis) information divergence reduces to the standard Kullback-Leibler divergence.

Associated to this (a, λ)-Rényi information divergences, we can also define a two parameters Rényi (Tsallis) entropy.

Definition 2. Given a and λ > 0 and a density f , we call (a, λ)-Rényi entropy the entropy of order q = a/λ of the escort f λ :

H a,λ [f ] = H a λ [f λ ] = 1 λ -a log f (x) a dµ(x) λ f (x) λ dµ(x) -a , (16)
provided that the different integrals exist. Similarly, the (a, λ)-Tsallis entropy is given by

S a,λ [f ] = S a λ [f λ ] = 1 λ -a f (x) a dµ(x) λ f (x) λ dµ(x) -a -1 , (17)
Again, for λ = 1 we get the standard Rényi (Tsallis) entropy of order a, for a = 1 we get the standard Rényi (Tsallis) entropy of order λ, for a = λ it reduces to the Shannon entropy of the escort distribution, and to the standard Shannon entropy in the case a = λ = 1.

In the following, we will need the fact that the divergence D (.) a,λ (f ||g) is always nonnegative: Proposition 3. if f and g are two densities such that the involved integrals are finite, with a and λ > 0, then D (.) a,λ (f ||g) ≥ 0, with equality if and only if f = g.

Proof.

Obviously, once we have realized that D (.) a,λ (f ||g) in ( 14) is the Rényi (Tsallis) divergence between escort distributions, this is a direct consequence of the nonnegativity of the Rényi (Tsallis) divergence, which follows from Jensen inequality. Alternatively, this could also be derived from the Hölder inequality.

The (a, λ) maximum entropy problem

With the previous definitions at hands, we can now consider the following extended maximum entropy problem:

H a,λ (m) = max f {H a,λ [f ] : m p,a [f ] = m and m 0,1 [f ] = 1} (18)
which includes the previous problems as particular cases. The usual procedure for handling such variational problem is the technique of Lagrangian multipliers, e.g. [START_REF] Tsallis | [END_REF]8]. However, even though the objective functional is strictly concave (here for a ≤ λ), the constraints set is not convex and the uniqueness and nature of the maximum can not be guaranteed. The situation is still more involved in the continuous setting considered here which requires the results from calculus of variations [START_REF] Gelfand | Calculus of Variations[END_REF][START_REF] Giaquinta | Calculus of Variations: The Lagrangian formalism[END_REF]. Therefore, we propose to derive the solution as consequence of the nonnegativity of D a,λ (f ||g), without recourse to the technique of Lagrange multipliers. We obtain that the density f (x) that achieves the maximum in the right hand side of ( 18) is a generalized Gaussian, as given by the following proposition.

Proposition 4. The density that achieves the maximum of the (a, λ)-Rényi or Tsallis entropy, subject to the escort mean constraint m p,a [f ] = m, is the generalized Gaussian

G β (x) = 1 Z(β) (1 -(λ -a) β|x| p ) 1 λ-a + for λ = a 1 Z(β) exp (-β|x| p ) if λ = a (19)
where β is a positive parameter. The partition function Z(β) is given by

Z(β) = 2 p (β) -1 p ×          (a -λ) -1 p B 1 p , -1 (λ-a) -1 p for a > λ > a -p (λ -a) -1 p B 1 p , 1 (λ-a) + 1 for λ > a Γ 1 p and if λ = a (20)
where B(x, y) is the Beta function. Finally, the a-moment of order p is given by

m p,a [G β ] = 1 βp(λ-a)( a λ-a + 1 p +1)
for λ -a = 0

1 βpa for λ = a. (21)
For λ > a, the density has a compact support, while for λ ≤ a it is defined on the whole real axis and behaves as a power distribution for |x| → ∞.

Proof. The expressions of the partition function and of the a-moment are obtained by direct calculation. For the main result, let us consider the generalized Gaussian (19), and denote A(β) = 1/Z(β) the inverse of the partition function. We directly have

f a G λ-a β dµ(x) = A(β) λ-a M a [f ] × (1 -(λ -a)β|x| p ) + f a M a [f ] dµ(x) (22)
≤ A(β) λ-a (1 -(λ -a)β m p,a [f ]) M a [f ], (23)
where we have exhibited the escort distribution f a , where m p,a [f ] denotes the generalized a-moment, and where the inequality in [START_REF] Gelfand | Calculus of Variations[END_REF] results from the fact that for λ > a, the support of (1 -(λ -a)β|x| p ) + can be smaller than the support of f. From ( 23) we also immediately get with

f = G β that M λ [G β ] = A(β) λ-a (1 -(λ -a)β m a,p [G β ]) M a [G β ]. (24)
Therefore, for all distributions f with a given a-moment m p,a [f ] = m and for β such that the generalized Gaussian has the same moment m p,a [G β ] = m, then the combination of ( 23) and [START_REF] Giaquinta | Calculus of Variations: The Lagrangian formalism[END_REF] gives

f a G λ-a β dµ ≤ M λ [G β ] M a [G β ] M a [f ].
Finally, the (a, λ)-Rényi divergence can be expressed as

D (R) a,λ (f ||G β ) = log    f a G λ-a β dµ(x) λ M λ [f ] a M λ [G β ] λ-a    1 a-λ ≤ log M a [G β ] -λ M λ [G β ] a M a [f ] -λ M λ [f ] a 1 a-λ = H a,λ (G β ) -H a,λ (f ),
by the definitions ( 14) and ( 16) of the divergence and entropy. By the nonnegativity of the divergence, we obtain that

H a,λ (G β ) ≥ H a,λ (f )
for all distributions f with a given a-moment m p,a [f ] = m p,a [G β ] = m, and with equality iff f = G β . Using the (a, λ)-Tsallis divergence rather than Rényi's, the same result also directly follows for the (a, λ)-Tsallis entropy.

From this result, we see that the problem (18) and its solution indeed interpolates between the maximum q-entropy problem with standard constraint [START_REF] Tsallis | [END_REF] and with generalized q-moment constraint (8).

• For a = q and λ = 1, the (q, 1)-Rényi entropy is the standard Rényi entropy with index q, whose maximum subject to a q-moment constraint is attained for the generalized Gaussian (10) with exponent 1/(1 -q).

• For a = 1 and λ = q, the (1, q)-Rényi entropy is the Rényi entropy with index q, whose maximum subject to a standard moment constraint is the generalized Gaussian (9) with exponent 1/(q -1).

It shall be also mentioned that we still get a q → 1/q duality result in this setting. From the definition (16), we always have the identity

H a,λ [f ] = H a λ [f λ ] = H λ a [f a ].
In the other hand, the generalized a-moment m p,a [f ] can be considered as the standard moment of density f a , m p,a [f ] = m p,1 [f a ] but also as the escort moment of order a/λ of the density f λ , i.e. m p,a [f ] = m p, a λ [f λ ]. Therefore, the problem (18) can be recasted in the two equivalent forms

H a,λ (m) = max f H λ a [f a ] : m p,1 [f a ] = m and m 0,1 [f ] = 1 (25)
which is a problem with respect to the density f a under a standard mean constraint, or

H a,λ (m) = max f H a λ [f λ ] : m p, a λ [f λ ] = m and m 0,1 [f ] = 1 (26)
which is a problem with respect to f λ with a generalized moment constraint.

The case of a two-level system

Finally, we close this paper with the example of a two-level system, with eigenenergies 0 and 1. Although the general (a, λ) maximum entropy problem (18) usually do not lead to closed-form expressions for the entropies H a,λ (m), see e.g. [START_REF] Bercher | [END_REF], we still can get an explicit solution for the simple two-level system. For this system, the measure at hands charges the two levels, with µ = δ 0 + δ 1 where δ x denotes the Dirac mass at x. Then, the maximum entropy distribution is the discrete distribution with a probability p for the excited state and (1 -p) for the ground state, with

p = (1 + β) 1 λ-a 1 + (1 + β) 1 λ-a .
In turn, the (a, λ)-entropy is then given by

H a,λ [G β ] = 1 λ -a log [p a + (1 -p) a ] λ p λ + (1 -p) λ -a ,
and the internal energy, m, which is computed with respect to the escort distribution of order a is simply

m = p a p a + (1 -p) a .
With a little algebra, it is possible to express p as a function of m:

p = m 1 a m 1 a + (1 -m) 1 a
, and then we can get the expression of the (a, λ)-entropy:

H a,λ (m) = a a -λ log m λ a + (1 -m) λ a
.

From this expression, it is easy to derive an expression of the inverse temperature as the derivative with respect to m, or the expression of the related heat capacity. We also see that the subsequent entropy only depends on the ratio λ/a. Finally, for λ = a, l'Hospital's rule gives

H a,a (m) = -m log [m] -(1 -m) log (1 -m) ,
which is the Fermi-Dirac entropy.

Conclusions

In this Letter, we have suggested a possible extension of the standard q-divergences and entropies by considering a combination of the concepts of q-divergences and of escort distributions. This leads to two-parameter information measures that recover the classical ones as special cases, and coincide with some recently proposed measures [21]. Further work should examine the general properties of these information measures. We have introduced a general maximum entropy problem and derived the expression of the canonical distribution maximizing the (a, λ)-entropy under a generalized a-moment constraint. This was obtained as a consequence of the nonnegativity of the related information divergence, without recourse to the technique of Lagrange multipliers. This canonical distribution is a two-parameter version of the generalized q-Gaussian distribution, extending the versatility of this distribution. It includes the standard solutions of the nonextensive thermostatistics as particular cases. Actually, the proposed generalized maximum entropy problem includes the standard approaches of nonextensive thermostatistics and provides a continuum of problems and solutions between them. In addition, this approach suggests that it is possible to adopt different indexes for the q-entropy and for the escort-mean constraint.


Abdourrahmane M Atto
email: abdourrahmane.atto@univ-savoie.fr

Yannick Berthoumieu
email: yannick.berthoumieu@ims-bordeaux.fr

Philippe Bolon
email: philippe.bolon@univ-savoie.fr

2-Dimensional Wavelet Packet Spectrum for Texture Analysis

Keywords: 2-D Wavelet packet transforms, Random fields, Spectral analysis, Spectrum estimation, Similarity measurements. I

The paper derives a 2-Dimensional spectrum estimator from some recent results on the statistical properties of wavelet packet coefficients of random processes. It provides an analysis of the bias of this estimator with respect to the wavelet order. The paper also discusses the performance of this wavelet based estimator, in comparison with the conventional 2-D Fourier-based spectrum estimator on texture analysis and content based image retrieval. It highlights the effectiveness of the wavelet based spectrum estimation.

that a large amount of frequencies contributes in the corresponding wavelet variances.

The present paper investigates the capability of 2-D wavelet packets to yield an accurate Power Spectral Density (PSD) estimator for characterizing second order statistical properties of random fields. In this paper, we consider the standard Fourier spectrum as benchmark and seek for estimating this spectrum by using suitable wavelets. This positions the Shannon wavelet at the focus of the paper: as the decomposition level tends to infinity, the bias of the Shannon wavelet spectrum estimator tends to 0. The bias of an arbitrary wavelet spectrum estimator then relates to the closeness of the wavelet under consideration to the Shannon wavelet. This closeness is measured through a parameter called wavelet order.

The results derived in the paper are in continuation of [START_REF] Atto | Wavelet packets of fractional brownian motion: Asymptotic analysis and spectrum estimation[END_REF] which introduces PSD estimation by using the statistical properties of the wavelet packet coefficients of random processes. The paper focuses on presenting the specificities that follow by dimensionality increasing from [START_REF] Atto | Wavelet packets of fractional brownian motion: Asymptotic analysis and spectrum estimation[END_REF]:

• the analytic form of the wavelet packet based PSD,

• the singular paths of fractional Brownian fields (Section III), • the 2-D wavelet packet PSD estimator (Section IV). Furthermore, for textured image analysis, the paper provides experimental results for evaluating the relevance of spectrum estimation by 1) comparing PSD estimated from Fourier and wavelet packet methods and 2) performing content based image retrieval associated with spectral similarity measurements in the Fourier and wavelet packet domains (Section V).

II. 2-D WAVELET PACKETS

This section provides a brief introduction to 2-D wavelet packet transforms. Advanced concepts and algorithms concerning 1D and 2-D wavelet packet analysis can be found in [START_REF] Wickerhauser | Adapted Wavelet Analysis from Theory to Software[END_REF], [START_REF] Mallat | A wavelet tour of signal processing[END_REF], [START_REF] Daubechies | Ten lectures on wavelets[END_REF]. The reader is also invited to refer to [START_REF] Craigmile | Asymptotic decorrelation of between-scale wavelet coefficients[END_REF], [START_REF] Kato | On the spectral density of the wavelet transform of fractional brownian motion[END_REF], [START_REF] Averkamp | Some distributional properties of the continuous wavelet transform of random processes[END_REF], [START_REF] Cambanis | On the continuous wavelet transform of second-order random processes[END_REF], [START_REF] Flandrin | Wavelet analysis and synthesis of fractional brownian motion[END_REF] (wavelets) and [START_REF] Atto | Wavelet packets of fractional brownian motion: Asymptotic analysis and spectrum estimation[END_REF], [START_REF] Atto | Wavelet packets of nonstationary random processes: Contributing factors for stationarity and decorrelation[END_REF] (wavelet packets) for more details on the statistical properties of wavelet transforms, when the decomposition relates to a random process.

A. Introduction on 2-D wavelet packets

We consider the 2-D separable wavelet packet decomposition in a continuous time signal setting for presenting theoretical results (see [START_REF] Mallat | A wavelet tour of signal processing[END_REF] for the connection between the continuous and the discrete wavelet transforms). In this decomposition, the wavelet paraunitary filters H 0 (low-pass, scaling filter) and H 1 (high-pass, wavelet filter) are used to split the input functional space U = W 0,0 ⊂ L 2 (R 2 ) into orthogonal subspaces (subbands):

W 0,0 = W 1,[0,0] ⊕ W 1,[1,0] ⊕ W 1,[0,1] ⊕ W 1,[1,1] , W j ,[n 1 ,n 2 ] = W j +1,[2n 1 ,2n 2 ] ⊕ W j +1,[2n 1 +1,2n 2 ] ⊕ . . . ⊕W j +1,[2n 1 ,2n 2 +1] ⊕ W j +1,[2n 1 +1,2n 2 +1]
for j 1. The subbands involved in this quad-tree splitting trick follow from a separable product W j ,[n 1 ,n 2 ] = W j ,n 1 ⊗ W j ,n 2 of 1D wavelet packet subspaces W j ,n 1 and W j ,n 2 . We assume that U is the subspace spanned by an orthonormal system

τ [k 1 ,k 2 ] Φ : (k 1 , k 2 ) ∈ Z 2 with Φ ∈ L 2 (R 2
) and where τ is the shift operator:

τ [k 1 ,k 2 ] f (t 1 , t 2 ) = f (t 1 -k 1 , t 2 -k 2 ). Denote F f , the Fourier transform of f ∈ L 1 (R 2 ) ∪ L 2 (R 2 ), with F f (ω 1 , ω 2 ) = R 2 f (t 1 , t 2 )e -i (ω 1 t 1 +ω 2 t 2 ) dt 1 dt 2 .
Since U does not necessarily depend on the wavelet filters, we assume in the following that U is the space of functions having their Fourier transforms with supports in [-π, π] × [-π, π]. In this respect, Φ is the sinc function:

Φ(t 1 , t 2 ) = sinc(t 1 , t 2 ) = sin(πt 1 ) sin(πt 2 ) π 2 t 1 t 2 .
The wavelet packet subband W j ,[n 1 ,n 2 ] is, by construction, the closure of the L 2 (R 2 ) subspace spanned by the orthornormal system of wavelet packet functions:

τ [2 j k 1 ,2 j k 2 ] W j ,[n 1 ,n 2 ] : (k 1 , k 2 ) ∈ Z 2 .
Due to separability, the wavelet packet function W j ,[n 1 ,n 2 ] satisfies:

F W j ,[n 1 ,n 2 ] (ω 1 , ω 2 ) = F W j ,n 1 (ω 1 )F W j ,n 2 (ω 2 ). ( 1
)
By decomposing the frequency indices n 1 , n 2 in the form

n 1 = j =1 1 2 j -, n 2 = j =1 2 2 j -, (2)
with ( 1 , 2 ) ∈ {0, 1} 2 for every ∈ {1, 2, . . . , j }, we have

F W j ,[n 1 ,n 2 ] = H j ,[n 1 ,n 2 ] F Φ, (3)
where the multiscale wavelet packet filter H j ,[n 1 ,n 2 ] applied to obtain the subband ( j , n) wavelet packet coefficients is

H j ,[n 1 ,n 2 ] (ω 1 , ω 2 ) = 2 i =1 H j ,n i (ω i ) with H j ,n i (ω) = 2 j /2 j =1 H i (2 -1 ω) (4) for i ∈ {0, 1}, i = 1, 2
, where H 0 and H 1 are the scaling and the wavelet filter introduced above. Assume that the scaling filter is with order r : H 0 ≡ H r 0 , where r is the largest non-negative integer such that [START_REF] Burrus | Introduction to Wavelets and Wavelet Transforms: A Primer[END_REF]

H r 0 (ω) = 1 + e -i ω 2 r Q(e i ω ) ( 5
)
and Q is a polynomial, with the following equality holding almost everywhere (a.e.):

lim r →∞ H r 0 (ω) a.e. = H S 0 (ω), (6)
filter H S 0 denoting the scaling filter associated with the Shannon wavelet. Then the 1D multiscale filters H r j ,n i i =1,2 have very tight supports when r is large (see [14, Figure 1] for illustration). These multiscale filters tend to the Shannon ideal filters H S j ,[n 1 ,n 2 ] when r increases (see Eq. ( 4) and Eq. ( 6)). The consequence is that the spectrum of the input random field is analyzed within a frequency window with small size. Indeed, the restriction of the support of

H S j ,[n 1 ,n 2 ]
to [-π, π]×[-π, π] is with size 2 -j ×2 -j . This makes possible a local analysis of the spectrum, as performed in Section IV below.

B. Wavelet packet paths

This section presents a specific wavelet packet path description derived from the binary sequence approach of [START_REF] Meyer | Wavelets, Algorithms and Applications[END_REF] for representing nested wavelet packet subspaces. This description is suitable for establishing asymptotic properties of 2-D wavelet packets with respect to the increase of the decomposition level.

For the sake of simplifying notation, 2-D wavelet packets are commonly presented through a quad-tree structure where a pair [n 1 , n 2 ] ∈ 0, 1, . . . , 2 j -1 2 of frequency indices is set in a one to one correspondence with an index n ∈ 0, 1, . . . , 4 j -1 . This correspondence leads to a single frequency indexing of the wavelet packet subbands. More precisely, this correspondence follows by associating a quaternary index µ ∈ {0, 1, 2, 3}

µ = 2 1 + 2 =        0 if ( 1 , 2 ) = (0, 0) 1 if ( 1 , 2 ) = (0, 1) 2 if ( 1 , 2 ) = (1, 0) 3 if ( 1 , 2 ) = (1, 1) (7)
to the binary indices ( 1 , 2 ) ∈ {0, 1} 2 and by defining:

n = n( j ) = j =1 µ 4 j -. ( 8
)
The above change of variable is a univoque correspondence between (n 1 , n 2 ) ∈ 0, 1, . . . , 2 j -1 2 defined by Eq. ( 2)

and n ∈ 0, 1, . . . , 4 j -1 defined by Eq. ( 8). This correspondence makes it possible to characterize the 4-band tree structure of the 2-D wavelet packet decomposition with a single shift parameter n, instead of the more intricate tree description with double indexed frequency indices. This characterization also makes it possible to use a straightforward path based description of the wavelet packet tree by associating an infinite quaternary sequence µ = µ 1 , µ 2 , µ 3 , . . . to a wavelet packet path P ≡ P µ and vice versa or, equivalently, an infinite sequence of frequency indices n = (n(1), n [START_REF] Nicolis | 2d wavelet-based spectra with applications[END_REF], n(3), . . .), those associated with µ to P ≡ P n from Eq. [START_REF] Daubechies | Ten lectures on wavelets[END_REF].

Let ( 1 k ) k∈N and ( 2 k ) k∈N be the binary sequences associated with µ from the correspondence given by Eq. [START_REF] Mallat | A wavelet tour of signal processing[END_REF]. Then, by taking into account Eq. ( 2), frequency index n P ( j ) is also associated with a pair of frequency indices (n 1 [P ]( j ), n 2 [P ]( j )).

In the rest of the paper, we will use notation n, n P ( j ) in general whereas notation n 1 , n 2 , n 1 [P ]( j ), n 2 [P ]( j ) will be used to specify the spatial location of the frequency support of the wavelet packet subband under consideration.

To end this section, it is worth mentioning that some specific paths will present singular behavior, depending on the input random field: the wavelet coefficients of certain non-stationary random fields on the subbands associated with these singular paths will remain non-stationary. As a matter of example,

• The separable fractional Brownian field analyzed in Section III-B1 admits frequency indices n( j ) such that n 1 ( j ) = 0 (resp. n 2 ( j ) = 0) for every j as singular frequency indices. The set of (singular) paths associated with these frequency indices will be denoted by

P n/n 1 =0 (resp. P n/n 2 =0
).

• The isotropic fractional Brownian field analyzed in Section III-B2 admits a unique singular path: the approximation path denoted by P 0 and associated with frequency indices n P 0 ( j ) = 0 for every j .

III. 2-D WAVELET PACKETS OF SECOND ORDER RANDOM FIELDS

A. The case of stationary random fields

Let X = X (t 1 , t 2 ) be a second-order zero-mean real valued WSS random field. We assume that X is continuous in quadratic mean.

The autocorrelation function of

X is R(t 1 , t 2 , s 1 , s 2 ) = E [X (t 1 , t 2 )X (s 1 , s 2 )]. We assume that R satisfies R 4 R(t 1 , t 2 , s 1 , s 2 )τ [2 j k 1 ,2 j k 2 ] W j ,n (t 1 , t 2 ) τ [2 j k 1 ,2 j k 2 ] W j ,n (s 1 , s 2 ) dt 1 dt 2 ds 1 ds 2 < ∞.
Under WSS assumption, this autocorrelation function reduces to a two variable function by using the convention:

R(t 1 , t 2 , s 1 , s 2 ) ≡ R(t 1 -s 1 , t 2 -s 2 ) = R(ρ 1 , ρ 2 )
where

ρ i = t i -s i for i = 1, 2.
We assume that γ = F R exists,

γ(ω 1 , ω 2 ) = R 2 R(ρ 1 , ρ 2 )e -i (ω 1 ρ 1 +ω 2 ρ 2 ) dρ 1 dρ 2 . ( 9
)
Function γ is the PSD of the WSS random field X . Under assumptions mentioned above, the coefficients of X on subband W j ,n define a discrete second order random field

c j ,n [k 1 , k 2 ]= R 2 X (t 1 , t 2 )τ 2 j [k 1 ,k 2 ] W j ,n (t 1 , t 2 ) dt 1 dt 2 . ( 10
)
Random field c j ,n inherits the WSS property of X : the autocorrelation function R j ,n of c j ,n ,

R j ,n [k 1 , k 2 , 1 , 2 ] = R 4 R(t 1 , t 2 , s 1 , s 2 )τ 2 j [k 1 ,k 2 ] W j ,n (t 1 , t 2 ) τ 2 j [ 1 , 2 ] W j ,n (s 1 , s 2 ) dt 1 dt 2 ds 1 ds 2 (11) reduces to R j ,n [k 1 , k 2 , 1 , 2 ] = R j ,n [k 1 -1 , k 2 -2 ] = R j ,n [m 1 , m 2 ]. Moreover, if we assume that PSD γ of X is bounded on the support of F W j ,n , then function R j ,n [m 1 , m 2 ] has the form R j ,n [m 1 , m 2 ] = 1 4π 2 R 2 γ(ω 1 , ω 2 )|F W j ,n (ω 1 , ω 2 )| 2 e i 2 j (m 1 ω 1 +m 2 ω 2 ) dω 1 dω 2 .
(12)

B. The case of non-stationary random fields: application to separable and isotropic fractional Brownian fields

Associating a wavelet packet based PSD to a non-WSS random field involves finding the WSS subbands (if any) associated with the decomposition of this field. Specifically, we consider in this section, the wavelet packet decompositions of some separable and isotropic constructions of random fields from the 1D fractional Brownian motion. Such a random field is called a fractional Brownian field (fBf ).

1) Separable fractional Brownian field:

This fBf, hereafter denoted by X s α = X s α (t 1 , t 2 )
, where α is the Hurst parameter, follows from a separable 2-D extension (cartesian product) of a 1D fractional Brownian motion. Its autocorrelation function is

R s (t 1 , t 2 , s 1 , s 2 ) = σ 2 4 |t 1 | 2α + |s 1 | 2α -|t 1 -s 1 | 2α × |t 2 | 2α + |s 2 | 2α -|t 2 -s 2 | 2α
where σ is a constant that can be tuned so as to control the variance of the fBf at a specific spatial coordinates. Assume that the wavelet function W j ,n = W j ,n (t 1 , t 2 ) is with compact support and has (at least) one vanishingmoment in each variable t 1 , t 2 . This assumption is satisfied for the frequency indices n such that n 1 = 0 and n 2 = 0. The corresponding subbands (4 j -2 j +1 + 1 subbands at decomposition level j ) are hereafter called detail wavelet packet subbands since their generating functions have vanishing moments in each direction.

The following Proposition 1 thus concerns the 4 j -2 j +1 +1 detail wavelet packet subbands. In these subbands, the wavelet packet moments vanish and the autocorrelation function R s j ,n of c s j ,n reduces to

R s j ,n [k 1 , k 2 , 1 , 2 ] = σ 2 4 R 4 |t 1 -s 1 | 2α |t 2 -s 2 | 2α τ [2 j k 1 ,2 j k 2 ] W j ,n (t 1 , t 2 ) τ [2 j 1 ,2 j 2 ] W j ,n (s 1 , s 2 ) dt 1 dt 2 ds 1 ds 2 . ( 13
)
Assume furthermore that function W j ,n is with compact support or has sufficiently fast decay, and that the Fourier transform of W j ,n has a bounded derivative in the neighborhood of the origin, then we can state the following Proposition 1.

Proposition 1: The discrete random field c s j ,n representing the coefficients of X s α on W j ,n is WSS when the frequency index n is such that n 1 = 0 and n 2 = 0. Its

autocorrelation function R s j ,n [k 1 , k 2 , 1 , 2 ] = R s j ,n [k 1 -1 , k 2 - 2 ] = R s j ,n [m 1 , m 2 ], is: R s j ,n [m 1 , m 2 ] = 1 4π 2 R 2 γ s (ω 1 , ω 2 )|F W j ,n (ω 1 , ω 2 )| 2 e i 2 j (m 1 ω 1 +m 2 ω 2 ) dω 1 dω 2 , (14)
where

γ s (ω 1 , ω 2 ) = σ 2 Γ 2 (2α + 1) sin 2 (πα) |ω 1 | 2α+1 × |ω 2 | 2α+1 (15)
and the Gamma function is defined by

Γ(u) = R + x u-1 e -x dx.
From Proposition 1, we derive that stationarity holds true for the 4 j -2 j +1 + 1 detail wavelet subbands. Let us focus on the frequency axes, that is when n 1 = 0 and n 2 = 0, or vice versa. For such an index n, we have

R s j ,n [k 1 , k 2 , 1 , 2 ] = Λ s j ,n [k 1 , k 2 , 1 , 2 ] + 1 4π 2 R 2 γ s (ω 1 , ω 2 )|F W j ,n (ω 1 , ω 2 )| 2 e i 2 j (m 1 ω 1 +m 2 ω 2 ) dω 1 dω 2 (16)
where, as above,

m i = k i -i for i = 1, 2 and Λ s j ,n [k 1 , k 2 , 1 , 2 ] = M 2α j ,0,k 2 × 2 R s j ,n 2 [m 2 ] if n 1 = 0, n 2 = 0 M 2α j ,0, 1 × 1 R s j ,n 1 [m 1 ] if n 2 = 0, n 1 = 0 and M 2α j ,n i ,k i is the 2α fractional moment of the 1D wavelet function τ 2 j k i W j ,n i : M 2α j ,0,k i = R |t | 2α τ 2 j k i W j ,0 (t ) dt , ( 17
)
with 1 R s j ,n 1 (resp. 2 R s j ,n 2
) denoting the autocorrelation function of a 1D section of an fBf (thus a fractional Brownian motion),

1 R s j ,n 1 [m 1 ] = 1 2π R γ * (ω 1 )|F W j ,n 1 (ω 1 )| 2 e i 2 j m 1 ω 1 dω 1 when n 1 = 0 and 2 R s j ,n 2 [m 2 ] = 1 2π R γ * (ω 2 )|F W j ,n 2 (ω 2 )| 2 e i 2 j m 2 ω 2 dω 2 ,
when n 2 = 0. In these equations, function γ * is given by

γ * (ω) = σ 2 Γ(2α + 1) sin(πα) |ω| 2α+1 . ( 18
)
One can notice the presence of an additive nonstationary term Λ s j ,n [k 1 , k 2 , 1 , 2 ] in Eq. ( 16). This term is a 2α centered fractional moment of the scaling function W j ,0 . This fractional moment is centered at -2 j k 2 when n 1 = 0 (resp. at -2 j 1 when n 2 = 0). The height of this fractional moment depends on the scaling function associated with the wavelet packet decomposition.

Note also that when both n 1 = n 2 = 0, the nonstationarity in wavelet coefficients is more intricate, mainly because the analyzing function has no vanishing moments in neither of the variables t 1 , t 2 .

2) Isotropic fractional Brownian field: This random field, denoted by

X i α = X i α (t 1 , t 2 )
is such that any of its 1D section by a line segment issued from the origin is a fractional Brownian motion. Its autocorrelation function is defined by

R i (t 1 , t 2 , s 1 , s 2 ) = σ 2 2 (t 2 1 + t 2 2 ) α + (s 2 1 + s 2 2 ) α - σ 2 2 (t 1 -s 1 ) 2 + (t 2 -s 2 ) 2 α .
By following a reasoning analog to that of Section III-B1 above, we have ) α W j ,n (t 1 , t 2 ) (compact support or sufficiently fast decay for W j ,n ) and provided that n = 0 (the wavelet function W j ,n then have at least one vanishing moment in some direction).

R i j ,n [k 1 , k 2 , 1 , 2 ] = - σ 2 2 R 4 (t 1 -s 1 ) 2 + (t 2 -s 2 ) 2 α τ [2 j k 1 ,2 j k 2 ] W j ,n (t 1 , t 2 ) τ [2 j 1 ,2 j 2 ] W j ,n (
Furthermore, if F W j ,n (ω 1 , ω 2 )/(ω 2 1 +ω 2 2 ) 3/4 is bounded in the neighborhood of the origin, then we have the following result:

Proposition 2: The random field c i j ,n = (c i j ,n [k 1 , k 2 ]) [k 1 ,k 2 ]∈Z×Z representing subband W j ,n coeffi- cients of X i α is WSS for n = 0, with autocorrelation function R i j ,n [k 1 , k 2 , 1 , 2 ] ≡ R i j ,n [k 1 -1 , k 2 -2 ] = R i j ,n [m 1 , m 2 ]
given by:

R i j ,n [m 1 , m 2 ] = 1 4π 2 R 2 γ i (ω 1 , ω 2 )|F W j ,n (ω 1 , ω 2 )| 2 e i 2 j (m 1 ω 1 +m 2 ω 2 ) dω 1 dω 2 , (20)
where

γ i (ω 1 , ω 2 ) = 2 -(2α+1) π 2 σ 2 sin(πα)Γ 2 (1 + α) 1 ω 2 1 + ω 2 2 α+1 . (21)
From Propositions 1 and 2, it follows that:

Remark 1: Except for a few frequency indices given in Propositions 1 and 2, the autocorrelation function of the wavelet packet coefficients of separable and isotropic fBfs can be written in the integral form involved in Eq. [START_REF] Cambanis | On the continuous wavelet transform of second-order random processes[END_REF].

From now on, when no confusion is possible, the autocorrelation close form given by Eq. ( 12) will also be used to denote the autocorrelation function of separable fBf (see Eq. ( 14)) or isotropic fBf (see Eq. ( 20)).

More precisely, when the autocorrelation function of the wavelet packet coefficients of an arbitrary, non-stationary random field can be written in the form given by Eq. ( 12) for a specific function γ, this function will be called the spectrum of the random field under consideration.

The following section addresses the estimation of spectrum γ. The method proposed relies on the frequency partitioning induced by the Shannon wavelet packets.

IV. 2-D WAVELET PACKET BASED SPECTRUM ESTIMATION

The spectrum estimation method presented in this section follows from the asymptotic analysis of the autocorrelation functions of the 2-D wavelet packet coefficients. This asymptotic analysis is performed with respect to the wavelet order r and the wavelet decomposition level j . When r increases, the asymptotic behavior of the sequence of wavelet functions is driven by the Shannon wavelet functions (see Eq. ( 6)). In this respect, we consider the Shannon wavelets in Section IV-A below and derive asymptotic results with respect to the wavelet decomposition level. The case of a wavelet with order r will be considered as an approximation of the Shannon limiting behavior when r is large enough.

A. Autocorrelation functions of the Shannon wavelet packets and asymptotic analysis

When the Shannon ideal paraunitary filters are used, the Fourier transform of the Shannon wavelet packet function W S j ,n satisfies:

F W S j ,n (ω 1 , ω 2 ) = 2 j 1l ∆ j ,G(n 1 ) ×∆ j ,G(n 2 ) (ω 1 , ω 2 ), (22)
where 1l K denotes the indicator function of the set K and

∆ j ,G(n) = ∆ - j ,G(n) ∪ ∆ + j ,G(n) , with ∆ - j ,G(n) = -∆ + j ,G(n) , ∆ + j ,G(n) = G(n)π 2 j , (G(n) + 1)π 2 j . ( 23
)
The function G in Eq. ( 23) is a permutation (see [START_REF] Wickerhauser | Adapted Wavelet Analysis from Theory to Software[END_REF], [START_REF] Mallat | A wavelet tour of signal processing[END_REF] for details) defined by G(0) = 0 and

G(2 + ) = 3G( ) + -2 G( ) + 2 , ( 24
)
where z is the largest integer less than or equal to z and ∈ {0, 1}. From Eqs. ( 12), ( 15), ( 21), [START_REF] Kwitt | Lightweight detection of additive watermarking in the DWT-domain[END_REF], the autocorrelation functions of the Shannon wavelet packet coefficients of the WSS random field described above, the separable fBf (with n 1 = 0 and n 2 = 0 when separable fBf is considered) and the isotropic fBf (assumes (n 1 , n 2 ) = 0 for the isotropic fBf ) can be written in the form

R j ,n [m 1 , m 2 ] = 2 2 j π 2 ∆ + j ,G(n 1 ) ×∆ + j ,G(n 2 ) cos (2 j m 1 ω 1 ) cos (2 j m 2 ω 2 )
γ(ω 1 , ω 2 ) dω 1 dω 2 (25) where γ is the spectrum of the input random field: γ is given by Eq. ( 9) for a WSS random field and γ is given by Eq. (15) (resp. Eq. ( 21)) for the separable (resp. isotropic) fBf.

The asymptotic analysis of these autocorrelations when j tends to infinity is performed by considering a path P of the 2-D wavelet packet decomposition tree. We assume that

• P = P 0 in the case of the isotropic fBf and

• P ∉ P n/n 1 =0 ∪ P n/n 2 =0 in the case of the separable fBf.
Thus, the sequence of autocorrelation functions R S j ,n P ( j ) j 1 of Shannon wavelet packets in path P is such that every R S j ,n P ( j ) , for j 1, has the form given by Eq. (25).

Let

ω 1 [P ] = lim j →+∞ G(n 1 [P ]( j ))π 2 j , ( 26
) ω 2 [P ] = lim j →+∞ G(n 2 [P ]( j ))π 2 j . ( 27
)
Assume that (ω 1 [P ], ω 2 [P ]) is a continuity point of the spectrum γ. Then, the sequence of autocorrelation functions R S j ,n P ( j ) satisfies:

lim j →+∞ R S j ,n [m 1 , m 2 ] = γ(ω 1 [P ], ω 2 [P ])δ[m 1 , m 2 ] (28) uniformly in (m 1 , m 2 ) ∈ Z 2 , where δ[•, •] is defined by δ[k 1 , k 2 ] = 1 if k 1 = k 2 = 0, 0 otherwise.

B. 2-D Wavelet packet based spectrum estimation

The following provides a non-parametric method for estimating spectrum γ of 2-D random fields on the basis of the convergence criteria given by Eqs. ( 6) and (28).

From Eq. ( 28), it follows that γ(ω 1 [P ], ω 2 [P ]) = lim j →+∞ R S j ,n [0, 0] so that the continuity points of spectrum γ can be estimated by subband variances (values R S j ,n [0, 0] n ), provided that the Shannon wavelet is used and j is large enough. Furthermore, we can derive from the convergence criteria given by Eq. ( 6), several spectrum estimators by considering wavelets with finite orders r (Shannon wavelet corresponds to r = +∞ ≡ S), the accuracy of the spectrum estimation being dependent on the wavelet order as shown in Proposition 3 below.

Assuming a uniform sampling (regularly spaced frequency plane tiling), the method applies upon the following steps.

2-D Wavelet Packet Spectrum

1) Define a frequency grid composed with frequency points (

p 1 π 2 j , p 2 π 2 j ) for p 1 , p 2 ∈ 0, 1, . . . , 2 j -1 (natural ordering).
2) Compute, from Eqs. ( 2), ( 7) and ( 8) the index n ∈ 0, 1, . . . , 4 j -1 (corresponding to the wavelet packet ordering) associated with (p 1 , p 2 ). 3) Set, for any pair (p 1 , p 2 ) given in step 1) and the corresponding n obtained from step 2),

γ p 1 π 2 j , p 2 π 2 j = var[c r j ,n ]. ( 29
)
From a practical point of view, the estimator of γ

p 1 π 2 j , p 2 π 2 j
given by Eq. ( 29) assumes i ) approximating the Shannon wavelet by a suitable wavelet and i i ) estimating the variance of c j ,n .

Approximating the Shannon wavelet can be performed by using Daubechies wavelets, depending to the Daubechies filter order r (see Eq. ( 5)), the larger the order, the more relevant the approximation. In practice, Daubechies filters with r = 7, 8 yield relevant approximations of the Shannon filters.

Estimating the variance Var[c j ,n ] of c j ,n can be performed by using the standard unbiased variance estimator (also called unbiased sample variance). In the Section dedicated to the experimental results, var[c j ,n ] will denote this unbiased estimator when formula of Eq. ( 29) is under consideration.

Assuming that an unbiased variance estimator is used, the bias of the above sample spectrum estimator depends on j , r and the spectrum γ shape. More precisely, we have: Proposition 3: Let X be a zero-mean random field having Shannon wavelet packet autocorrelation function of the form given by Eq. ( 25). Let p 1 , p 2 ∈ 0, 1, . . . , 2 j -1 .

Let n be the frequency index associated with p 1 , p 2 in the sense of the correspondence ( j , n)

≡ ( j , [n 1 , n 2 ]) = j , [G -1 (p 1 ),G -1 (p 2 )] given in Section II. Assume that γ is continuous on ∆ + j ,p 1 × ∆ + j ,p 2
and differentiable on the interior of ∆ + j ,p 1

× ∆ + j ,p 2
. Assume that the standard unbiased estimator of variance is used in Eq. ( 29). We have:

Bias γ p 1 π 2 j , p 2 π 2 j πC ∇γ e 2 j + C γ 4π 2 ε r j ,n , (30)
where constants e =

2 + ln(1 + 2) /3 < 1, C ∇γ = sup ∆ + j ,p 1 ×∆ + j ,p 2 ∇γ , C γ = sup ∆ + j ,p 1 ×∆ + j ,p 2 γ and ε r j ,n R 2 |F W r j ,n (ω 1 , ω 2 )| 2 -|F W S j ,n (ω 1 , ω 2 )| 2 dω 1 dω 2 .
The proof sketch is the following.

Proof: The bias of the estimator given by Eq. ( 29) is:

Bias γ p 1 π 2 j , p 2 π 2 j = E var[c r j ,n ] -γ p 1 π 2 j , p 2 π 2 j = Bias var[c r j ,n ] + Var[c r j ,n ] -γ p 1 π 2 j , p 2 π 2 j . (31)
Thus, if we consider the standard unbiased estimator of variance (Bias var[c r j ,n ] = 0), we obtain:

Bias γ p 1 π 2 j , p 2 π 2 j Var c r j ,n -Var c S j ,n + Var[c S j ,n ] -γ p 1 π 2 j , p 2 π 2 j . ( 32
)
On one hand, we have

Var[c r j ,n ] -Var[c S j ,n ] C γ 4π 2 × R 2 |F W r j ,n (ω 1 , ω 2 )| 2 -|F W S j ,n (ω 1 , ω 2 )| 2 dω 1 dω 2 . ( 33
)
On the other hand, we have

Var[c S j ,n ] -γ p 1 π 2 j , p 2 π 2 j = 2 2 j π 2 ∆ + j ,p 1 ×∆ + j ,p 2 γ(ω 1 , ω 2 ) -γ p 1 π 2 j , p 2 π 2 j dω 1 dω 2 2 2 j π 2 ∆ + j ,p 1 ×∆ + j ,p 2 γ(ω 1 , ω 2 ) -γ p 1 π 2 j , p 2 π 2 j dω 1 dω 2 . ( 34
)
Since γ is continuous on

∆ + j ,p 1 × ∆ + j ,p 2
and differentiable on the interior of ∆ + j ,p 1 × ∆ + j ,p 2 , we have, from the mean value theorem:

γ(ω 1 , ω 2 ) -γ p 1 π 2 j , p 2 π 2 j sup (u 1 ,u 2 )∈∆ + j ,p 1 ×∆ + j ,p 2 ∇γ(u 1 , u 2 ) × ω 1 - p 1 π 2 j 2 + ω 2 - p 2 π 2 j 2 (35) for any (ω 1 , ω 2 ) ∈ ∆ + j ,p 1 × ∆ + j ,p 2
, with

∆ + j ,p 1 ×∆ + j ,p 2 ω 1 - p 1 π 2 j 2 + ω 2 - p 2 π 2 j 2 dω 1 dω 2 = e π 2 j 3 .
(36) Proposition 3 then follows from Eqs. (32), ( 33), (34), ( 35) and (36).

From Proposition 3, we derive that the bias of the estimator given by Eq. ( 29) depends on the decomposition level and wavelet order used. This bias tends to 0 when both j and r tends to infinity. Notice that ε r j ,n tends to 0 as r tends to infinity, ∀ j , n, due to the convergence of the filter H r j ,n to the Shannon filter H S j ,n as r tends to infinity.

V. EXPERIMENTAL RESULTS

This section provides experimental results on spectral analysis of textures. Section V-A addresses the quality of the wavelet packet spectrum estimator to capture textural information. Section V-B provides quantitative evaluation of the wavelet packet spectrum estimator in the framework of content based image retrieval.

A. Spectral analysis and spectral texture contents

Wavelet packet spectra 1 of some texture images are provided in Figures 1 and2. The wavelet packet spectra have been computed from the method given in Section IV-B, where the decomposition level is 6 and the Daubechies wavelet with order r = 7 is used. Spectra computed from the Fourier transform are also given in this figure, for comparison purpose.

From a visual analysis of images given in Figures 1 and2 (by focusing on image interpretation without watching spectra given in the same figures), one can remark that most of these textures exhibit non-overlapping textons replicating repeatedly: thus, coarsely, we can distinguish several frequencies having significant variance contributions (from a theoretical consideration), when the texture does not reduce to the replications of a single texton.

In addition, when these textons occupy approximately the same spatial area (see for instance "Fabric" textures in Figure 2), the frequencies with high variance contributions (peak in the spectrum) are close in terms of their spatial location (from a theoretical consideration). 1 Colors represented in Figures 1 and2 are simulated from a light source in order to ease 3-D visualization: red color [value 1] corresponds to fully illuminated shapes whereas blue color [value 0] is associated to shaded areas, green color corresponds to value 0.5. The above heuristics, issued from visual image analysis, are confirmed by considering the wavelet packet spectra (see for instance spectra of "Fabric" textures in Figure 2), whereas, in most cases, the two dimensional discrete Fourier transform exhibits only one peak.

One can highlight that the poorness of the Fourier spectra is not due to a lack of resolution in the sampling step of the Fourier transform. This poorness can be explained by noting that Fourier transform is sensitive to global spatial regularity. In contrast wavelet packets can capture local spatial regularity and lead to a more informative spectrum estimator when several frequencies contribute in texture variance distribution. Note that when we have some a priori on the spectrum shape, we can apply a non-uniform sampling scheme by simply focusing on the particular tree structure associated with the subbands-of-interest. These subbands are those associated with wavelet packet functions having tight Fourier transform support across the frequencies where the spectrum exhibits sharp components.

B. Content based image retrieval by using spectral texture features

Section V-A have shown qualitatively that wavelet packets are relevant in analyzing textures with rich spectral content. In this section, we evaluate the contribution of wavelet packets in Content Based Image Retrieval (CBIR) from a database composed with such textures.

The database under consideration is composed with M = 94 images associated with the 6 texture classes given in Figures 1 and2. These M images are obtained from the splitting of every 'large' texture image into 16 nonoverlapping images. An image from the database thus has 160 × 160 pixels when it pertains to one of the classes "D3", "D10", "D17" or "D87" and 128×128 pixels for classes "Fabric.09" and "Fabric.11" (the larger images have 640×640 pixels in the Brodatz album and 512 × 512 pixels in the VisTeX album).

In order to compare a query image I q with an arbitrary image I , 1 q, M , we use the Log-Spectral distance and the symmetric version of the Kullback-Leibler divergence. For images I 1 , I 2 with sample spectra γ 1 , γ 2 , their Log-Spectral distance is given by

LS(γ 1 , γ 2 ) ∝ q, log γ 2 1 (ω q , ω ) γ 2 2 (ω q , ω ) 1/2 (37)
and their symmetric Kullback-Leibler divergence is defined by

K L(γ 1 , γ 2 ) = K L(γ 1 ||γ 2 ) + K L(γ 2 ||γ 1 ) (38)
where

K L(γ 1 ||γ 2 ) ∝ q, γ 1 (ω q , ω ) log γ 1 (ω q , ω ) γ 2 (ω q , ω )
Table I provides the CBIR retrieval rates per texture and the average CBIR retrieval rate over the 6 texture classes. In this table, the above spectral similarity measurements are computed upon discrete Fourier, wavelet and wavelet packet spectra. The Daubechies wavelet of order r = 7 and j = 4 have been used to compute wavelet and wavelet packet spectra from Eq. ( 29).

This table emphasizes

• the suitability of wavelet based approaches, in comparison with the Fourier approach, • the sensitivity of the wavelet transform to the similarity measures considered, which may follow from the coarse-irregular spectrum sampling induced by using wavelet method. • the good performance of the wavelet packet spectral analysis for the Log-Spectral and Kullback-Leibler similarity measurements. One can note that the wavelet transform is a particular case of the wavelet packet transform so that the wavelet spectrum can be seen through the wavelet packet spectrum by zooming on the neighborhood of the zero frequency. This wavelet spectrum thus follows from a non-uniform sub-sampling of the wavelet packet spectrum (non-regular frequency grid). From similarity evaluations, this sampling scheme can be interpreted as a binary sparse weighting, where the non-zero values are concentrated in low frequencies. This sampling scheme consequently penalizes highly medium and high frequencies. This strategy is suitable for textures with spectral supports localized in the neighborhood of zero (up to 90 % of textures from natural images). However, this strategy can justify the contrasted results obtained by using the wavelet transform for the tests performed upon the database of Table I: the spectral content of these textures is not restricted to low frequencies.

In addition to the above binary weighting regarding spectra sample points, the log function involved in Log-Spectral and Kullback-Leibler measures impacts differently the sample spectra γ 1 (ω q , ω ), γ 2 (ω q , ω ) under comparison. Indeed, for large values of γ 1 (ω q , ω )/γ 2 (ω q , ω ), the log function imposes a higher penalty in Kullback-Leibler measurements than in Log-Spectral measurements.

VI. CONCLUSION

Two issues are addresses in the paper: (i ) associating a PSD to the wavelet packet coefficients of random fields and (i i ) estimating this PSD from the statistical properties of the wavelet packet coefficients.

Issue (i ) has been addressed for stationary random fields and some non-stationary fractional Brownian fields. Issue (i i ) has been tackled from asymptotic properties of the Shannon wavelet packets and the spectrum estimation method proposed is more effective for wavelet filters with large order.

The PSD estimation method derived in the paper has shown relevancy in estimating the PSD of texture images. This method outperforms Fourier based spectral estimation when the textures under consideration have substantial spectral content (when spectrum do not reduces to a peak at single frequency).

Prospects regarding this work may concern:

(p 1 ) The extension of the method to adaptive sampling of the PSD: a best basis can be derived in order to achieve an efficient strategy. Indeed, the experimental results have shown that for many textures, a large frequency domain is unoccupied whereas estimating several spectrum sample points in this domain could increase estimation errors due to the decimation step involved by the wavelet packet decomposition. (p 2 ) The extension of the method with respect to nonorthogonal wavelet packet transforms: the redundancy of these transforms will probably benefit to the spectrum estimation. Note that the spectrum needs to be defined with respect to the transform and its estimation cannot reduce to variances per subbands unless any given subband has been processed accordingly.

Some applications of the method are:

(a 1 ) Hurst parameter estimation for self-similar medical images, see for instance [START_REF] Nicolis | 2d wavelet-based spectra with applications[END_REF].

(a 2 ) texture modeling by using Wold decompositions [START_REF] Liu | Periodicity, directionality, and randomness: Wold features for image modeling and retrieval[END_REF], [START_REF] Stitou | Three-dimensional textured image blocks model based on wold decomposition[END_REF]: estimation the poles of the spectrum is necessary to determine the spectral singularities involved in the deterministic texture contribution. These poles are associated with peaks of the spectrum and their number, as well as their location determines the accuracy of the modeling. (a 3 ) Spread-spectrum image watermarking [START_REF] Bors | Image watermarking using block site selection and DCT domain constraints[END_REF], [START_REF] Wang | Wavelet tree quantization for copyright protection watermarking[END_REF], [START_REF] Gkizeli | Optimal signature design for spread-spectrum steganography[END_REF], [START_REF] Kwitt | Lightweight detection of additive watermarking in the DWT-domain[END_REF], [START_REF] Cui | Adaptive multiwavelet-based watermarking through jpw masking[END_REF], textural change detection by using spectral distances [START_REF] Vu | S3: a spectral and spatial measure of local perceived sharpness in natural images[END_REF], among other references. We conclude by noting that the wavelet packet decomposition tree can be chosen adaptively with respect to the spectrum shape in order to makes possible non-uniform sampling schemes: the analysis performed in the paper suggests estimating the spectrum from 1) a full wavelet packet decomposition especially for applications involving images with rich spectral contents (for instance when the input image has many spectral singularity points, as in the case of the multifractional brownian fields) and 2) a sub-tree of the full wavelet packet tree, this sub-tree being deleafed in frequency domains associated with less spectral information (in particular, a wavelet spectrum follows from a non-uniform sampling grid associated with a tree that expands mainly in low frequencies).

FTFig. 1 .

 1 Fig. 1. Textures images and their spectra γ computed by using discrete Fourier and wavelet packet transforms. Abscissa of the spectra images consist of a regular grid over [0, π/2] × [0, π/2].

Remark 2 :

 2 In the above experimental results, we have considered a regular frequency grid over [0, π] × [0, π]. Consequently, we have performed a uniform sampling of the spectrum of the input texture.

FTFig. 2 .

 2 Fig. 2. Textures images and their spectra γ computed by using discrete Fourier and wavelet packet transforms. Abscissa of the spectra images consist of a regular grid over [0, π/2] × [0, π/2].

  s 1 , s 2 ) dt 1 dt 2 ds 1 ds 2

						(19)
	under	integrability	condition	for	the	term
	(t 2 1 + t 2 2

TABLE I RETRIEVAL

 I RESULTS PER TEXTURE AND AVERAGE RETRIEVAL RESULTS FOR THE 6 TEXTURE CLASSES GIVEN IN FIGURES 1 AND 2. EXPERIMENTAL RESULTS ARE PERFORMED WITH THE LOG-SPECTRAL SPECTRAL DISTANCE AND THE KULLBACK-LEIBLER DIVERGENCE. SPECIFICALLY, FT, WT AND WPT DESIGNATE THE DISCRETE FOURIER, WAVELET AND WAVELET PACKET TRANSFORMS, RESPECTIVELY.

	Measure:	Log-Spectral	Kullback-Leibler
	Texture	FT	WT	WPT	FT	WT	WPT
	D3	44.5	88.7	99.2	39.5	58.2	91.0
	D10	71.5	82.8	90.2	90.6	58.2	85.6
	D17	93.0	100	99.6	100	99.6	100
	D87	77.3	100	96.1	87.9	87.5	99.2
	Fabr.09	84.4	100	83.6	82.8	80.5	98.1
	Fabr.11	69.5	95.7	77.0	74.6	73.4	97.3

ACKNOWLEDGEMENT

The authors address special thanks to the anonymous reviewers and the associate editor for their constructive comments.