\section{Introduction}
Recent practices of medicine strongly rely on the visualization of inner organs
in order to study the relationship of anatomic structure and physiological process
of living beings.
In this context, vibro-acoustography (VA)~\cite{fatemi98} has become a medical imaging tool for investigating
calcifications in human arteries~\cite{greenleaf98,fatemi00,fatemi03},
microcalcifications in the breast~\cite{fatemi01,wold04}, liver tumors~\cite{alizad04},
bone structure~\cite{calle01}, among others applications~\cite{urban:2011}.

The VA technique employs two co-focused ultrasound beams in megahertz range with slightly different
frequencies, typically in the kilohertz range, to scan biological tissues
(see Fig.~\ref{fig:figure1} for the system description).
The VA image formation process is related to the nonlinear interaction of the primary ultrasound waves (incident and scattered),
which gives rise to a secondary acoustic wave at the difference-frequency~\cite{silva:1326,silva:2011}.
The generated difference-frequency wave is detected by a hydrophone (typical bandwidth $1$--$\unit[70]{kHz}$).
The acquired signal is processed and used to form VA images.

The interaction of the primary incident ultrasound waves gives rise to a highly collimated difference-frequency beam,
known as the parametric array effect~\cite{westervelt:535}.
It has been recently shown that the difference-frequency signal
detected in the forward scattering direction is mostly due to the nonlinear interaction of the primary scattered
waves~\cite{silva:2011}.
Hence, the parametric array may behave as a background
signal present in the difference-frequency wave detected to form VA images.

The VA method can be extended to generate several images at once using a
multifrequency arrangement also known as multifrequency VA~\cite{urban06}.
In this configuration, $N$ ultrasound beams (more than two) are produced and
they overlap in the system focal region.
Due to the nonlinear interaction of the ultrasound waves,
$N(N-1)/2$ difference-frequency components arise in the medium and they are used
to form multiple images with just one scanning process.
\begin{figure}[htb!]
\centering
 \includegraphics[width=\linewidth]{figure1.pdf}
 \caption{Description of a typical VA imaging system.
 The central- and the difference-frequency are $\omega_0$ and $\Delta \omega$, respectively.}
 \label{fig:figure1}
\end{figure}

The spatial resolution of VA described in terms of the system point spread function (PSF)
is determined by the transducer geometrical properties and the ultrasound driving frequencies.
Beamforming has a great impact on the system resolution.
Some schemes to generate two or more ultrasound beams for VA  have been proposed based on a confocal transducer~\cite{fatemi99},
$x$-focal (crossing beams)~\cite{chen04},
sector array~\cite{glauber04}, linear array~\cite{glauber04,urban11}, and reconfigurable array~\cite{kamimura:}.
The VA resolution does not to correspond to the difference-frequency beam waist
due to the parametric array effect.
Indeed, the system PSF has a narrower lateral width than that of a parametric beam produced by the VA confocal transducer.
To verify this, we have computed the parametric beam with the theory presented in Ref.~\cite{ding:2759}
and found that this beam does not fit the PSF waist shown in Ref.~\cite{fatemi99}.
However, for the sake of brevity, these results will not be included in this article.

The study presented here is limited to VA systems based on the confocal arrangement.
Specifically, a confocal VA system composed by a two-element confocal transducer
with $\unit[44]{mm}$-diameter, focused at $\unit[70]{mm}$, and
operating with frequencies around $\unit[3.1]{MHz}$, achieves a spatial resolution
($\unit[-6]{dB}$ beamwidth)
of approximately \unit[$0.8\times 0.8$]{mm} in the focal plane and \unit[$16$]{mm} axially.
With this resolution only images formed at the focal plane are feasible.
Furthermore, compared to conventional pulse-echo ultrasound images, VA images generated by the confocal transducer do
not exhibit speckle noise.

Clinical applications of VA will demand systems based on linear array transducers because the ultrasound beams
can be electronically steered over the region to be imaged.
It has been shown through computational simulations that a typical linear array
with $\unit[32]{mm}$-aperture operating at $\unit[3]{MHz}$ can produce
VA images with $\unit[-6]{dB}$ of lateral resolution of $\unit[1.9]{mm}$ when
focused at $\unit[5]{cm}$ in depth~\cite{glauber04}.
However, the system PSF is not symmetric in the imaging plane as in other VA beamforming schemes.
The resolution is poorer in the elevation when compared to the lateral direction (roughly threefold larger).
Hence, linear array systems (1.25 and 1.5D) may produce distorted VA images.
Recently, a VA system based on the General Electric Vivid 7 ultrasound scanner with the GE 7L and 10L
linear array transducers was designed and assessed~\cite{urban11}. This system was set to operate at
nearly $\unit[5]{MHz}$, with the difference-frequency being $\unit[51.54]{kHz}$.
When the incident beams were focused at $\unit[25]{mm}$ in depth, the spatial resolution
achieved by the system at $\unit[-6]{dB}$ was
$\unit[0.573]{mm}$ laterally, $\unit[1.44]{mm}$ in elevation,  and $\unit[2.20]{mm}$ axially.

Limitations in the spatial resolution of VA systems degrade the yielded images, reducing
their utility to medical imaging applications.
The need for restoration of VA images in either elevation (with linear arrays) or axial directions
has prompted us to assess the efficiency of some restoration algorithms applied to these images.
In doing so, three deconvolution frequency domain filters, namely constrained least squares (CLS),
geometric mean, and Wiener filters~\cite{gonzalezcap5,hunt73} based on the theoretical PSF of the system are
employed to restore simulated VA images~\cite{perciano08}.
In addition, a phantom composed by three stretched $\unit[0.5]{mm}$ wires was imaged by a confocal VA system.
The VA line spread function (LSF) was estimated from the image and used to restore the image through the CLS filter.
Restoration with the theoretical LSF yielded results with no visual information of the wires.
The main reason for this problem lies on the lack of the phase information in the theoretical LSF.
This is illustrated by restoring the actual VA image with a combined LSF, in which the magnitude and the phase come from
the theoretical and the experimental LSF, respectively.
In this case, the CLS filter restoration is capable to restore most of the information such as position, shape and contrast
of the wires in the image.

This article is organized as follows.
A model for  vibro-acoustic image formation based on a three-dimensional PSF composed by the product
of the pressure amplitudes of the incident beams is presented in Sec.~\ref{sec:va}.
The general aspects of the restoration filters (Wiener, CLS, and geometric mean) are briefly described in
Sec.~\ref{sec:restoration}.
The computational and experimental results are explained in Sec.~\ref{sec:results}.
In Sec.~\ref{sec:conclusions}, the main conclusions of this study are outlined.

\section{Image formation}
\label{sec:va}
The VA technique employs two ultrasound beams at frequencies $\omega_1$ and $\omega_2$ to form
an image of the region-of-interest.
Hence the primary incident pressure upon this region is given by
$$
p(\mathbf{r},t) = p_1(\mathbf{r}) e^{-i\omega_1 t} + p_2(\mathbf{r})
e^{-i\omega_2 t},
$$
where $\mathbf{r}$ is the position vector indicating the observation point, $i$ is the imaginary unit,
$t$ is the time,
and $p_1$ and $p_2$ are the complex pressure amplitudes.

Based on the radiation force concept, a 2D image formation model was proposed in which
the system PSF is the product of the incident pressure fields in the transducer
focal plane~\cite{fatemi99}.
An extension of this model to account for depth-of-field effects in the VA image formation was presented in
Ref.~\cite{glauber06}.
According to these image formation models, the three-dimensional PSF of a VA system is given by
\begin{equation}
\label{psf1}
 h(\mathbf{r})  = p_1^*(\mathbf{r})p_2(\mathbf{r}),
\end{equation}
where the symbol * means complex conjugation.
It is worthy noting that the PSF of a VA system as given in Eq.~(\ref{psf1}) is a complex function.
Despite early success on explaining VA as a radiation force imaging method,
a recent study~\cite{silva:2011} has shown that the difference-frequency signal used to form VA images
depends on the nonlinear interaction of primary incident and scattered waves at the frequencies $\omega_1$ and $\omega_2$.
Both radiation force phenomenon and the nonlinear interaction of the waves scale quadratically with the primary
pressure amplitudes.
Hence, the description of the difference-frequency pressure magnitude based upon the radiation force and the nonlinear
interaction models are alike, but the phase information due to each phenomenon is different~\cite{silva:2011}.
Therefore, as the deconvolution of VA images needs to be done with the PSF (magnitude and phase).
Thus, we can expect that the deconvolution using the PSF given in Eq.~(\ref{psf1}) may not be able
to fully recover the visual information in VA images.

To compute the PSF given in Eq.~(\ref{psf1}) one has
to calculate the incident pressure field amplitudes $p_1$ and $p_2$ in the linear approximation for the VA transducer.
This has been done using the Field II simulation program~\cite{jensen:262}.
Even though the fields from Field II are linear the combination through equation 2 makes the problem inherently nonlinear.
The nonlinearity of the VA beamforming arises because of the interaction of the two beams,
not the nonlinear propagation of the pressure from the transducer.
The simulations were performed considering
the transducer and operation frequencies used in the experiments to be described later.
Hence, a two-element spherical confocal transducer with a
$\unit[44]{mm}$ diameter and focus
at $\unit[70]{mm}$ was considered.
The transducer inner radius is $\unit[14.8]{mm}$, while the outer ring ranges
from $15.2$ to $\unit[22]{mm}$.
The transducer elements were driven independently by sinusoidal signals
at $\unit[3.075]{MHz}$ and $\unit[3.125]{MHz}$.
The magnitudes of the simulated PSF laterally and in depth are presented in Fig.~\ref{fig:psf}.
\begin{figure}[hbt!]
 \centering
 \mbox{\subfigure{\label{fig:figure2}\includegraphics[width = 0.7\linewidth,angle=-90]{figure2.pdf}}}
 \mbox{\subfigure{\label{fig:figure3}\includegraphics[width = 0.44\linewidth,angle=-90]{figure3.pdf}}}
\caption{Magnitude of the PSF in the focal plane and in depth for a two-element confocal transducer
with inner and outer radii being respectively $14.8$ and $\unit[22]{mm}$.
The transducer focal distance is $\unit[70]{mm}$, which corresponds to zero in the depth scale.}
\label{fig:psf}
\end{figure}

The VA image of an object can be obtained by convolving the function which represents the object with the
system PSF.
Accordingly, the image function is given by~\cite{glauber06}
\begin{equation}
g(\mathbf{r})=f(\mathbf{r})\star h(\mathbf{r}) +
n(\mathbf{r}),
\label{eq1}
\end{equation}
where $f$ is the function that represents the object, $n$ represents the additive Gaussian noise, and
the symbol $\star$ denotes the spatial convolution.
The additive Gaussian noise is the only source of noise present in VA images
due to thermal fluctuations present in the electronic devices of the VA system.
The noise level of VA images has been assessed.
Thus, a white noise with
$\unit[20]{dB}$ signal-to-noise ratio (SNR) was chosen to be added to all simulated images here.
This level was considered by visual inspection of the VA simulated and actual images.

In the current arrangement of VA, real-time images cannot be done because the system requires
four or five difference-frequency cycles to form a pixel.
At $\unit[50]{kHz}$, this corresponds to $\unit[80]{\mu s}$ per pixel.
To form images with $100 \times 100$ pixels and show them in $24$ frames per
second, one needs $\unit[4.1]{\mu s}$ for each pixel.

\section{Restoration algorithms}
\label{sec:restoration}

Among different filtering schemes to perform image deconvolution, we choose the Wiener, the  constrained least-squares (CLSs), and the geometric mean (GM) filters~\cite{gonzalezcap5}.
This choice is based on the following digital image processing aspects.
Firstly, the chosen filters operate in the frequency domain, which is
computationally less expensive than to perform the filtering in the spatial domain.
Secondly, inasmuch as the VA image formation is based on a spatial
invariant model with additive noise as given in Eq.~(\ref{eq1}),
one can use inverse filter methods for deconvolution.
In this scheme, the Wiener filter removes the image degradation without division by small values, which can lead to numerical overflow.
However, this filter reduces the noise level by smoothing the image without sharpening the edges.
On the other hand, the CLS filter reduces the image degradation
and it allows the control of sharpening the edges through a regularization parameter.
The CLS and the Wiener filters were used for the deconvolution of
ultrasound images~\cite{yeoh2006}.
The GM filter combines the features of the Wiener and the CLS filters depending on the adjustment of two parameters.
The criteria used here to set the filters parameters will be discussed later.

To  describe the Wiener, the CLS, and the GM filters,
the image formation model in Eq.~(\ref{eq1}) is written in
the frequency domain as
\begin{equation}
 G({\bf k}) = F({\bf k})H({\bf k}) + N({\bf k}),
\label{eq3}
\end{equation}
where  ${\bf k}$ is the reciprocal of ${\bf r}$, and $F$, $G$, $H$, and $N$ are the Fourier transforms of the functions $f$,
$g$, $h$, and $n$, respectively.

\subsection{Wiener filter}

The Wiener filter is designed to minimize the mean squared
error between the filtered and the ideal images.
In the frequency domain, this filter is given by~\cite{gonzalezcap5}
\begin{equation}
F = \left[\frac{H^*}{|H|^2 + S_n/S_f}\right]G,
\end{equation}
where $S_n$ and $S_f$ are the power spectrum density of the noise and of the ideal image, respectively.

\subsection{Constrained least-squares filter}
In the frequency domain, the constrained least-squares
(CLSs) filter is expressed as~\cite{gonzalezcap5}
\begin{eqnarray}
F = \left[\frac{H^*}{|H|^2+\gamma|P|^2}\right]G,
\end{eqnarray}
where $\gamma$ is the regularization parameter and $P$ is the Fourier Transform of the discrete Laplacian operator in three-dimensions.

\subsection{Geometric mean filter}
The geometric mean (GM) filter is given in the frequency
domain by~\cite{gonzalezcap5}
\begin{eqnarray}
F =
\left\{\left[\frac{H^*}{|H|^2}\right]^\alpha\left[\frac{H^*}{|H|^2+\gamma
S_n/S_f}\right]^{1-\alpha}\right\}G,
\end{eqnarray}
where $\alpha \in [0, 1]$ and $\gamma>0$.

\subsection{Filter parameters estimation}
\label{subsec:parameters}
In performing the VA image deconvolutions, a trial and error
approach is adopted.
In doing so, the best visual and quantitative results
for the restored images are sought.
The quantitative analysis is performed using three indices of digital image processing,
namely~\cite{wang02,WangBovikLu:ICASSP:02} improvement in the signal-to-noise ratio (ISNR),
mean squared error (MSE), and universal image quality index (UIQI).
For the last two indices, the imaginary part of the restored image is ignored,
because  it approaches to zero after the filtering process.

\section{Results and discussion}
\label{sec:results}

\subsection{Computational simulations}
\label{subsec:simulations}

To perform the computational simulation of VA images, the PSF is represented by a
complex matrix whose dimensions are $256\times 256\times 512$.
This corresponds to a region of \unit[$32$]{mm} laterally and \unit[$64$]{mm} in depth.
Two computational phantoms are devised with three and five spherical inclusions
as shown in Fig.~\ref{fig:figure4}.
The phantoms have dynamic range of $\unit[48]{dB}$, i.e., a pixel is
represented with $\unit[8]{bits}$ depth.
Later, we present two computer phantoms with inclusions having up to four contrast levels.
The purpose of using such phantoms is to demonstrate that the deconvolution of VA images can recover most information regarding to
the inclusions morphology.
Thus, the adopted dynamic range, which represents inclusions with $256:1$ contrast ratio, is sufficient for this proposal.
In turn, the inclusions can be regarded as soft tissue with different mechanical properties compared to the background.

\begin{figure}[hbt!]
 \centering
 \includegraphics[width=\linewidth]{figure4.pdf}
 \caption{Phantoms used in the computational restoration.
(a) Phantom 1 seen in the \unit[6]{mm}-transverse plane.
(b) Phantom 1 seen in-depth with three inclusions with \unit[$4$]{mm}
 radius each and  $\unit[36]{dB}$ (left),
$\unit[0]{dB}$ (center), and $\unit[42]{dB}$ (right) contrasts.
(c) Phantom 2 seen in the \unit[8]{mm}-transverse plane.
(d) Phantom 1 seen in-depth with two additional inclusions
of \unit[$2$]{mm} radius and $\unit[46]{dB}$ contrast.
The dotted lines in (b) and (d) indicates the position of the images (a) and (b).
}
 \label{fig:figure4}
\end{figure}
Phantom 1 is composed by three inclusions having same radius of  \unit[$4$]{mm}.
The inclusions have contrasts of $\unit[36]{dB}$ (left),
$\unit[0]{dB}$ (center),  and $\unit[42]{dB}$ (right).
Phantom 2 is similar to phantom 1 with two additional spheres are included with \unit[$2$]{mm} radius each and $\unit[46]{dB}$
contrast.
The dimensions of the inclusions are typically found in VA images of mass legions and calcifications.

The VA computational images are generated by convolving the system PSF with the function
which represents the phantom, as given in Eq.~\eqref{eq1}.
A Gaussian white noise of $\unit[20]{dB}$ signal-to-noise ratio (SNR) was added
to the generated images.
The simulated images are presented in Fig.~\ref{fig:figure5}.
It is not possible to recognize any inclusion from the axial image information.
\begin{figure}[hbt!]
 \centering
 \includegraphics[width=\linewidth]{figure5.pdf}
 \caption{VA images obtained by computational simulations with white noise of
 $\unit[20]{dB}$ SNR. (a) Image
 of phantom 1 seen in the \unit[6]{mm}-transverse plane.
 (b) Phantom 1 seen along depth. (c) Image of phantom 2
 at the \unit[8]{mm}-transverse plane.
 (d) Phantom 2 shown along depth.
The dotted lines in (b) and (d) indicates the position of the images (a) and (b).}
 \label{fig:figure5}
\end{figure}

The simulated images are restored with the CLS, geometric mean and Wiener filters.
The obtained results for both phantoms are shown in Figs.~\ref{fig:figure6} and~\ref{fig:figure7}.
It is clear that all details in the axial images are recovered by the filters analyzed here.
\begin{figure*}[htb!]
\centering
 \includegraphics[width=0.56\linewidth,angle=-90]{figure6.pdf}
 \caption{Image restoration results for phantom 2 on both
 the $\unit[6]{mm}$-transverse plane and in-depth.
The filter parameters are $\gamma= 0.2$ (CLS), and $\alpha=0.5$ and $\gamma=1$ (geometric mean).
 The CLS filter yielded the best boarder preservation,
 while the Wiener filter yielded the most blurred image.
The dotted line indicates where the respective transverse image is placed in depth.}
 \label{fig:figure6}
\end{figure*}
\begin{figure*}[htb!]
\centering
 \includegraphics[width=0.59\linewidth,angle=-90]{figure7.pdf}
 \caption{Image restoration results for phantom 2 on both
 the $\unit[8]{mm}$-transverse plane and in-depth.
The filter parameters are $\gamma= 0.2$ (CLS), and $\alpha=0.5$ and $\gamma=1$ (geometric mean).
 The image restored with the Wiener filter is more blurred than the other results.
The CLS filter yielded the best edge preservation among the other filters. }
 \label{fig:figure7}
\end{figure*}

In Fig.~\ref{fig:figure8}, the profile of phantom 2  and the corresponding restored images at $\unit[$8$]{mm}$ in depth are displayed in normalized scale.
The VA image profile in Fig. 7(b) shows a noisy mixture
of two gray levels present in phantom 2.
After the restoration process, the gray levels resemble
that of phantom 2.
\begin{figure*}[hbt!]
 \centering
 \subfigure[Phantom 2]{\includegraphics[width=0.28\linewidth,angle=-90]{figure8a.pdf}\label{fig:figure8a}}
 \subfigure[VA image]{\includegraphics[width=0.28\linewidth,angle=-90]{figure8b.pdf}\label{fig:figure8b}}\\
\subfigure[Wiener filter]{\includegraphics[width=.28\linewidth,angle=-90]{figure8c.pdf}\label{fig:figure8c}}
\subfigure[CLS filter]{\includegraphics[width=.28\linewidth,angle=-90]{figure8d.pdf}\label{fig:figure8d}}
\subfigure[GM filter]{\includegraphics[width=.28\linewidth,angle=-90]{figure8e.pdf}\label{fig:figure8e}}
 \caption{Gray levels obtained for a transverse line at $\unit[8]{mm}$ in depth of phantom 2.
 \subref{fig:figure8a} Transverse line for phantom 2.
 \subref{fig:figure8b} Transverse line for the simulated VA image of phantom 2.
 Transverse lines for the restored images using~\subref{fig:figure8c} Wiener filter,
 \subref{fig:figure8d} CLS filter and~\subref{fig:figure8e} GM filter.
}
 \label{fig:figure8}
\end{figure*}

The quantitative analysis is performed using the measures defined in Sec.~\ref{subsec:parameters}.
Deconvolution restoration results for phantoms 1 and 2 are presented in Table~\ref{tab:table1}.
The values indicate good results for the three restoration algorithms
because the ISNR gives \unit[60]{dB} improved with
respect to the raw image.
Moreover, the MSE is lower than $0.02$ and the UIQI is higher than $0.77$ (when this index values one, it indicates
``perfect'' restoration).
The best restoration result was achieved by the Wiener filter for both phantoms.
We observe that the resulting images are smoother than those obtained
with the two other filters.
In other words, the restored Wiener image is more blurred and
the object borders are mixed with the background.
Furthermore,  the constrained least-squares filtered images display the best border preservation because
it is based on the Laplacian operator, which enhances the object borders.
\begin{table*}[ht!]
\caption{
Quantitative results based on the ISNR (improvement in signal-to-noise ratio),
MSE (mean squared error), and UIQI (universal image quality index) for
phantoms 1 and 2.
}
\label{tab:table1}
\begin{center}
\begin{tabular}{l|c|c|c||l|c|c}
\toprule
& \multicolumn{3}{c||}{\small{Phantom 1}} & \multicolumn{3}{c}{\small{Phantom 2}}\\\midrule
\small{Filters} & \small{ISNR} & \small{MSE} & \small{UIQI} & \small{ISNR} & \small{MSE} & \small{UIQI}\\
\midrule
\small{Constrained least-squares} & \small{$60.13$} & \small{$0.0192$} & \small{$0.78$} & \small{$59.91$} & \small{$0.0188$} & \small{$0.79$}\\
\small{Geometric mean} & \small{$60.74$}  & \small{$0.0166$} & \small{$0.80$} & \small{$59.83$} &  \small{$0.0191$} &  \small{$0.77$}\\
\small{Wiener} & \small{$61.83$} & \small{$0.0129$} & \small{$0.91$} & \small{$61.03$} & \small{$0.0145$} & \small{$0.88$}\\
\bottomrule
\end{tabular}
\end{center}
\end{table*}
\normalfont

\subsection{Experimental images}
To investigate how actual VA images can be restored, an experiment using a VA system is performed and
the results are presented here.
The VA system consists of the two-element confocal transducer described in
Sec.~\ref{sec:va}, which scans a wire phantom (see Fig.~\ref{fig:figure9}).
The driving frequencies of the transducer elements
are $3.075$ and $\unit[3.125]{MHz}$
and the corresponding difference-frequency is
$\unit[50]{kHz}$.
In turn,
the wire phantom is composed of three stretched \unit[$0.5$]{mm}-diameter wires distributed along the
depth-of-field.
The wires are diagonally arranged and spaced by \unit[$1.4$]{cm}.
The scanned area is \unit[$30 \times 70$]{ mm} and each pixel was taken in steps of \unit[$0.1$]{mm}.
This area encloses all depth-of-field information of the acquired images.
A hydrophone (ITC-6050C, ITC, Santa Barbara, CA, USA) with sensitivity $\unit[-157 \pm 3]{dB~re~1 V/\mu Pa}$ is
used to detect the generated difference-frequency pressure.
The hydrophone is placed at approximately $\unit[r=20]{cm}$ way from the
farthest wire to the transducer.
This corresponds to $kr=41$, which can be regarded as the hydrophone
farfield region.
The active element of the hydrophone is a cylinder located in the middle of its conical tip.
This cylinder has a height of $\unit[50]{mm}$ and diameter of $\unit[23]{mm}$.
Its axis is aligned with the the hydrophone long axis.
The hydrophone cannot be considered as a point detector for difference-frequency wavelength of
few centimeters ($\unit[3.0]{cm}$ at $\unit[50]{kHz}$ in water).
The finite-size effect of the active element (in a linear model) yields an output voltage proportional
to the spatial average of the detected pressure over the element area plus a signal-independent noise~\cite{boutkedjirt:745}.
Consequently, the hydrophone may add an extra information to the phase of the acquired signal.
Fig.~\ref{fig:exp_setup} shows the experimental VA setup used.
\begin{figure}[hbt!]
\centering
 \includegraphics[width=.6\linewidth,angle=-90]{figure9.pdf}
 \caption{Image of the computational representation of the wire phantom.}
 \label{fig:figure9}
\end{figure}

\begin{figure}[hbt!]
 \centering
 \subfigure[Side view.]{\includegraphics[width=0.42\linewidth]{figure10}\label{fig:figure10}}\qquad
 \subfigure[Top view.]{\includegraphics[width=0.42\linewidth]{figure11}\label{fig:figure11}}
 \subfigure[View from transducer.]{\includegraphics[width=0.42\linewidth]{figure12}\label{fig:figure12}}
 \caption{VA experimental setup composed of a two-element confocal transducer, a phantom with three stretched thin wires, and a hydrophone. All components are submerged in a water tank.}
 \label{fig:exp_setup}
\end{figure}

Two schemes are devised to restore the actual VA images.
Firstly, we estimate the line-spread function (LSF) directly from the
experimental image (central wire) and use it in the deconvolution process.
In order to make the experimental LSF more similar to the theoretical one, we applied a spectral
inversion to the image, i.e., first a low-pass filtering is performed and
then the filtered image is subtracted from the original image~\cite{smith97}.
Secondly, because the magnitude of both the theoretical and the estimated LSF are similar,
 we compose the LSF by combining the magnitude of
Eq.~\eqref{psf1} and the phase information directly from the experimental data.
Such arrangement was chosen because we found discrepancies between the
phase of the estimated and the theoretical LSF
based on the radiation force model.
This model seems not to provide the entire phase information present in VA images.
A more complete image formation description needs to take into account parametric and nonlinear scattering phenomena.

The magnitude and phase of the estimated and the theoretical LSFs are exhibited in
Figs.~\ref{fig:figure13} and~\ref{fig:figure14}.
The theoretical LSF was obtained by the convolution of the 3D VA PSF with the 3D wire phantom.
The phase of the experimental LSF is estimated as follows.
Each pixel is formed from an acquired harmonic sequence.
This sequence discretized from the detected signal by the hydrophone at the difference-frequency $\Delta\omega$.
Suppose that the time sequence which represents the $i$th-pixel of the VA image is given by
\begin{equation}
\label{fm}
 x_i[n] = A_i\cos(\Delta \omega n/f_s + \phi_i), \quad n = 0, 1, \dots, N,
\end{equation}
where $A_i>0$ and $\phi_i$ are the is the phase of the signal, $f_s=\unit[1]{MHz}$ is the
sampling frequency, and $N=512$ the number of sampling points.
By taking the Hilbert transform~\cite{wolfram:hilbert} of $x_i[n]$, one obtains
$\tilde{x}_i[n]  = A_i\sin(\Delta \omega n/f_s + \phi_i).
$
Thus, the representation of the acquired signal in complex notation is
\begin{equation}
x_i[n]+ i\tilde{x}_i[n] = A_i e^{i(\Delta \omega n/f_s + \phi_i)}.
\end{equation}
Therefore, the phase of a sequence which represents the $j$th-pixel with
respect to that of the $i$th-pixel
can be  computed from
\begin{equation}
\phi_j = \phi_i + \frac{i}{A_i A_j N} \sum_{n=0}^N \ln\left[(x_i[n] + i\tilde{x}_i[n])(x_j[n] - i\tilde{x}_j[n])\right].
\end{equation}
All the pixel phases can be obtained with respect to a reference pixel.
The amplitudes $A_i$ and $A_j$ can be calculated by computing
the root mean square of the respective acquired signal.

It can be noticed from Figs.~\ref{fig:figure13} and~\ref{fig:figure14} that
the theoretical and the experimental LSFs have differences in their phase contents.
The lack of information in the LSF phase leads to unsatisfactory restoration results with several artifacts.
The missing information in the phase is mostly related to the fact that the
image formation model based on Eq.~(\ref{psf1}) does not consider the nonlinear interaction of the primary beams.
Moreover, the finite-size effect of the hydrophone active element was not considered in the theoretical model.
\begin{figure}[hbt!]
 \centering
\subfigure[]{\includegraphics[width=0.57\linewidth,angle=-90]{figure13.pdf}\label{fig:figure13}}
\subfigure[]{\includegraphics[width=0.57\linewidth,angle=-90]{figure14.pdf}\label{fig:figure14}}
\caption{The magnitudes of (a) the experimental and (b) the theoretical LSFs.
The corresponding phases in radians of the center line along the depth-of-field are also displayed. \label{fig:lsf}}
\end{figure}

The result obtained by estimating the LSF
from the experimental image, is presented in Fig.~\ref{fig:figure16}, while
the result using the composed LSF is shown in Fig.~\ref{fig:figure17}.
Both methods yielded similar restored images.
Despite presenting two blurred wires (top left and bottom right), the restored images show aspects of the
imaged wires, including location, geometric shape and contrast.
Some aspects may have contributed for the occurrence of artifacts in the restored images.
The deconvolution filters used here usually produce undesirable artifacts during the
process of deconvolution because of the ill-posed nature of the image restoration problem, i.e.,
it may have more than one solution and the solutions depend discontinuously on the initial data.
Therefore, perturbations in the VA image due to the noise, may cause undesired variations in the restoration results.
This can be noticed by comparing the result of a complete simulated deconvolution
process presented in Fig.~\ref{fig:figure18}, where no artifacts are produced, because the
whole process is controlled, from the image formation up until its restoration.
That is, the same LSF is used for both VA image formation and restoration.
Furthermore, it is important to note that an extra information may appear in the phase of the detected signal
caused by the acoustical path difference from the wire positions to the hydrophone.
\begin{figure}[hbt!]
 \centering
 \subfigure[]{\includegraphics[width=\linewidth]{figure15.pdf}\label{fig:figure15}}
 \subfigure[]{\includegraphics[width=\linewidth]{figure16.pdf}\label{fig:figure16}}
 \subfigure[]{\includegraphics[width=\linewidth]{figure17.pdf}\label{fig:figure17}}
  \subfigure[]{\includegraphics[width=\linewidth]{figure18.pdf}\label{fig:figure18}}
 \caption{Image restoration using the CLS filter with $\gamma=0.01$. The raw VA image is shown in (a).
The restoration results are presented for the experimental (b) and the composed (c) LSFs.
 The deconvolution for the computer generated VA image is shown in (d). }
 \label{fig:image}
\end{figure}

\section{Summary and conclusions}
\label{sec:conclusions}
In this paper, we have presented the results of three restoration techniques (Wie\-ner,
CLS and geometric mean filters) applied in the restoration of VA images.
The restoration algorithms based on the theoretical PSF were applied to computationally simulated images, in which the visual information
of the objects was lost.
The performance of the algorithms were assessed both quantitatively and visually.
The deconvolution of simulated images was capable of restoring the size, shape, position and contrast of the objects present in the images for all filters used.

The CLS filter was employed in restoring an actual VA image of three stretched wires.
In this case, two LSFs were used, one was estimated from the image, while the other
combines the phase of the estimated LSF with the magnitude of the theoretical one.
The wire positions and partial information about their shape and contrast can be recognized
in the images restored with both experimental and combined LSFs.

The restoration of the experimental VA image with the theoretical LSF based on
Eq.~\eqref{psf1} yielded a poor visual restoration retrieve and
thus was not presented here.
Yet the magnitudes of both theoretical and experimental LSFs matched,
the phases showed distinctly different contents.
This hints that the adopted theoretical model for the PSF lacks the appropriate phase information of VA images.
Additional information may come from the difference-frequency generation based on
the nonlinear interaction of the primary waves~\cite{silva:1326,silva:2011,westervelt:535}
and possibly also on the finite-size effect of the hydrophone active element.

This study represents a preliminary step in understanding how application of deconvolution filters can be used to
restore the resolution and contrast for VA images. The optimal restoration filter may depend on the object
being imaged with to improve resolution and contrast.
Future studies would be necessary to understand and
validate how different filters and modeling of the PSF would affect the restoration of VA images of simple
phantoms and soft tissues.

\section{Acknowledgments}
This work was partially supported by FAPESP and Grants 306697/2010-6,
477653/2010-3 CNPq (Brazilian agencies).

\section*{References}

\small


\section{Introduction}

Lie's Group theory approach to the study of differential equations has had
a profound and enduring impact in Theoretical Physics, in particular in what respects to
  the invariance of differential equations under point
transformations \cite{stephani}. He showed that the one-dimensional free
particle equation exhibits the symmetries corresponding to  the
eight-dimensional SL(3,R) group of point transformations and that this is
the maximum number of symmetry generators for a second-order differential
equation of the form \cite{stephani}

\begin{equation}
\label{16.1}\ddot y+h(\dot y,y,x)=0.
\end{equation}

Of course, one often encounters situations in which one deals with a smaller
number of generators. In a variety of such cases \cite{chime}, equation (\ref
{16.1}) can nonetheless be recast in the form of a free particle equation by
recourse to a {\it nonlocal variable transformation}. Now, as nonlocal
transformations can change both  the number of symmetry generators  and the
very physics of the original problem, they become a powerful tool for
theoretic analysis.
Einstein's field equations, whose nonlinear nature can indeed be faced with
Lie's weaponry, constitute a particularly relevant example in this respect.

It is via this equivalence (under nonlocal variable transformations) between
different physical problem that we will search here for connections between
\begin{itemize}
\item Cosmological problems in General Relativity, on the one hand,  and
\item Information Theory
Concepts pertaining to {\it both} Fisher's Information measure
for translation families (FIM) \cite
{f0,fisher,f1,f2,f3,f4,f5,f6,f7,f8,f9,f10} and Tsallis' Nonextensive
Thermostatistics (NET) \cite{t1,t2,t3,t4,t5,t6,t7,t8,t9,miller}, on
the other one.
\end{itemize}

This search
is motivated by the recent findings of Ref. \cite{Penni}: an intriguing
relation between homogeneous and isotropic spatially flat cosmological
models with no cosmological constant, on the one hand, and FIM-probability
distributions that solve {\it the diffusion equation, as adapted a
non-extensive ($q=-1$) Tsallis environment}, on the other one.

FIM's applications to diverse problems in theoretical physics have been
pioneered by Frieden, Soffer, Nikolov, and others \cite
{f1,f2,f3,f4,f5,f6,f7,f8,f9,f10}, who have unveiled many FIM properties and
clarified its relation to Shannon's logarithmic information measure

\begin{equation}
\label{entrop1}S\,=\,-\,<ln f> ,
\end{equation}

\noindent
where $f$ is, of course, a probability density and $<A>$ stands for the
trace of $(fA)$.

General Relativity is a classical theory of gravitation that has been
successfully applied for solving  problems concerning astronomical scales
of length and time, such as the formulation of cosmological models.
Einstein's equation for the classical geometry adopts the appearance

\begin{equation}
\label{einstein}R_{ik}-\frac 12g_{ik}R=8\pi GT_{ik},
\end{equation}

\noindent where $R_{ik}$ is the Ricci tensor, $R$ the scalar curvature and
T_{ik}$ the energy-momentum tensor of matter and fields, with $i,k=1,..,D$,
and D, the spacetime dimension.

During the last few years much attention has been paid to models in
which the universe,  at very early stages of its evolution, was
expanding  \cite{luis,luis1,ale,jakubi}. The well known Big Bang
model suggests a spatially homogeneous and isotropic space-time,
stemming from a singularity in the remote past \cite{Weinberg}.

Here we will show that a generalization of  the so-called
Frieden-Soffer principle of Extreme Physical Information \cite{f7},
within the framework of  a non extensive setting, may yield original
insights into the workings of some of Einstein's field equations.

As a first step in such a direction, and in order to establish a convenient
notation, we begin by refreshing, in Section II: i) the salient points of
Fishers's measure tenets for translation families,
ii) their extension to non-extensive settings, and
iii) the Frieden Soffer Principle concerning extremalization of Physical
Information \cite{f7}.  In Section III we discuss the application of this
Principle to a non-extensive environment. We present the main (cosmological)
results of the present Communication in Section IV and, finally, some
conclusions are drawn in Section V.

\section{A review of fundamental notions}

\subsection{The estimation theory bound}

In this work we shall concern ourselves {\it exclusively} with Fisher
information
measures $I$ for translation families, that are
invariant under Galilean and Lorentz translations \cite{corcuera}. They
refer to a measure of the inverse uncertainty in determining a {\it position}
parameter by a maximum likelihood estimation \cite{silver}.

Estimation theory \cite{f1} provides one with a powerful result that
constitutes our stepping stone. Consider an
isolated many-particle system that is specified by a {\it position}
parameter $x$
and let $f(x)$ describe the
 probability density (pd) for this parameter. The mean value of $x$
for this pd is $\eta$. Suppose the pd is unknown, and one wishes to
determine it. Estimation theory asserts \cite{f1} that the best possible
estimator $\eta_{est}$ of our parameter, after a very large number of
samples is examined, suffers a mean-square error $e^2$ from $\eta$ that
obeys a relationship involving Fisher's $I$, namely,

\begin{equation}
\label{fisher}e^2=\frac{1}{I} ,
\end{equation}

\noindent
where Fisher's information measure $I$ (for tranlation families) reads

\begin{equation}
I\,=\,\int \,dx\,f(x)\,\{ \frac{\frac{df}{dx}}{f(x)} \}^2,
\end{equation}

\noindent
i.e.,

\begin{equation}
\label{ifisher}I\,=\,<\,\{\frac{\frac{df}{dx}}{f(x)}\}^2\,>.
\end{equation}
Any other  estimator must have a larger mean-square error. The only proviso
to the above result is that all estimators $\eta _{est}$ be unbiased, i.e.,
satisfy

\begin{equation}
\label{unbias}<\eta_{est}>\,=\,\eta .
\end{equation}

Thus, Fisher's information measure for translation families has a lower bound, in the sense that, no
matter what parameter of the system we choose to measure, $I$ has to be
larger or equal than the inverse of the mean-square error associated with
the concomitant measurement. This result, i.e.,
\begin{equation}
\label{rao}I\,\ge\,\frac{1}{e^2},
\end{equation}
is referred to as the Cramer-Rao bound, and constitutes a very powerful
statistical result that was generalized in \cite{miller}  to a NET scenario.

\subsection{Estimation in a nonextensive Tsallis setting}

The phenomenal success of thermodynamics and statistical physics crucially
depends upon certain mathematical relationships involving energy and
entropy, that can be translated in an essentially intact form into Tsallis'
generalized Thermostatistics formalism \cite{t1,t2,t3,t4,t5}, whose main ingredient is the generalized entropy $S_q$, that in terms of the real parameter $q$ reads

\be \label{sq}
S_q\,=\,  \langle \,\, \frac{1}{q-1}\,\, (1\,-\,f^{q-1})\,\, \rangle,\,\,\,q\,\epsilon\,{\cal R},
\ee
and coincides with Shannon's $S$ for $q=1$.  Indeed, much
work has been recently devoted to i) show that many of these relationships
are valid for {\it arbitrary} values of Tsallis' parameter $q$, and ii) to
find appropriate generalizations for the rest. In this vein we just mention
that, by suitably maximizing the entropy of the system, Curado and Tsallis
\cite{t3} found that the whole mathematical (Legendre-transform based)
structure of thermodynamics becomes {\it invariant} under a change of the $q-
$value (from unity to any other real number), while the connection of NET
both with quantum mechanics and with Information Theory was established in
\cite{t4}, where it was shown that all of the conventional
Jaynes-Boltzmann-Gibbs \cite{katz} results suitably generalize to the
Tsallis' environment. For more details see \cite{t1}.

Of course, to verify that NET is useful, it is necessary to show that it
appropriately describes certain physical systems with $q$-values that are
different from unity. Much work in this respect has been  performed recently
(144 refereed papers at the time of this writing) . We may cite applications
to astrophysical problems \cite{t5,t6}, to L\'evy flights \cite{t7}, to
turbulence phenomena \cite{t8}, to simulated annealing \cite{t9}, etc. The
interested reader is referred to \cite{t1} for additional references.

A suitable generalization of (\ref{ifisher}) to a NET setting should involve
replacing any probability distribution $f$ involved in evaluating
expectation values by Tsallis' generalized weights $f^q$. As shown in \cite
{miller}, one finds that Fisher's Information Measure in a Tsallis's
nonextensive setting (generalized FIM $I_q$) reads

\begin{equation}
\label{qfisher}I_q\,=\,< \,\{ \frac{\frac{df}{dx}}{f(x)} \}^2\,>_q ,
\end{equation}

\noindent
which can be abbreviate as GFIM. Instead of (\ref{unbias}) we have to deal,
of course, with
\begin{equation}
\label{qunbias}<\,\eta _{est}(y)\,>_q\,=\,\eta .
\end{equation}
In \cite{miller} one uses (\ref{qfisher}) to

\begin{itemize}
\item  generalize the Cramer-Rao bound \cite{f1}.

\item  discuss connections between Tsallis' entropy, on the one hand, and
Fisher's measure, on the other one.

\item  derive a suitably generalized form of the Frieden-Nikolov's bound
\cite{katz} to the time derivative $\frac{dS_q}{dt}$. This alternative form
connects the entropy increase to $I_q$.
\end{itemize}

\subsection{ The Frieden-Soffer EPI Principle}

The Principle of Extreme Physical Information (EPI) is an overall physical
theory that is able to unify several sub-disciplines of Physics \cite{f7}.
In Ref. \cite{f7} Frieden and Soffer (FS) show that the Lagrangians in
Physics arise out of a mathematical game between an intelligent observer and
Nature (that FS personalize in the appealing figure of a ``demon",
reminiscent of the celebrated Maxwell's one). The game's payoff introduces
the EPI variational principle, which determines simultaneously the
Lagrangian {\it and} the physical ingredients of the concomitant scenario.

FS \cite{f7} envision the following situation: some physical phenomenon is
being investigated so as to gather suitable, pertinent data. Measurements
must be performed. Any measurement of physical parameters appropriate to the
task at hand initiates a relay of information $I$ (or $I_q$) from Nature
(the demon) into the data. The observer acquires information, in this
fashion, that is precisely $I$ (or $I_q$). FS assume that this information
is  stored within the system (or inside the demon's mind). The demon's
information is called, say, $J$ (or $J_q$) \cite{f7}.

Assume now that, due to the measuring process, the system is perturbed,
which in turn induces a change $\delta J$ ($\delta J_q$ )of the demon's
mind. It is natural to ask ourselves how the data information $I$ will be
affected. Enters here FS's EPI: {\it in its relay from the phenomenon to the
data no loss of information should take place}. The ensuing new Conservation
Law states that $\delta J=\delta I$, or, rephrasing it
\begin{equation}
\label{epi}\delta (I\,-\,J)\,=\,0,
\end{equation}
so that, defining an action ${\cal A}$
\begin{equation}
\label{action}{\cal A}\,=\,I\,-\,J,
\end{equation}
EPI asserts that the whole process described above extremalizes ${\cal A}$.
FS \cite{f7} conclude that the Lagrangian for a given physical environment
is not just an {\it ad-hoc} construct that yields a suitable differential
equation. It possesses an intrinsic meaning. Its integral represents the
physical information ${\cal A}$ for the physical scenario. On such a basis
some of the most important equations of Physics can be derived \cite{f7}.
For an interesting Quantum Mechanical derivation see \cite{Penni1}.

Within the present context our demon is working in a non-extensive fashion,
so that the EPI principle should read
\begin{equation}
\label{actionq}{\cal A}_q\,=\,I_q\,-\,J_q,
\end{equation}
i.e., using appropriate $q$-generalizations of the intervening physical
quantities. Thus the Frieden-Soffer information transfer game is played here
according to
\begin{equation}
\label{epiq}\delta (I_q\,-\,J_q)\,=\,0.
\end{equation}

\section{The workings of FS's demon in a non-extensive scenario}

\subsection{Introductory remarks}

One starts here the FS game by using (\ref{epiq}), i.e., by extremalizing
Fisher's generalized information  for translation families
\begin{equation}
\label{1}I_q=\int f(x)^{q-2\ }\ \dot{f}(x)^2dx,
\end{equation}
and considering a physical scenario in which the knowledge of the
information demon \cite{f7} is limited to that concerning the normalization
constraint. The demon's $J_q$ is
\begin{equation}
\label{2}J_q\,=\,\int \gamma _0\ f(x)\ dx,
\end{equation}
where $\gamma _0$ is the pertinent Lagrange multiplier. The demon's
information is reduced here to its minimum-core expression.

Playing the Frieden-Soffer game then leads to
\begin{equation}
\label{e8}-2\;\ddot{f}(x)+(2-q)\ \frac{\dot{f}(x)^2}{f(x)}-\gamma
_0f(x)^{2-q}=0,
\end{equation}
$q$-dependent, non-linear differential equation that should  yield our
``optimal" probability distribution $f$. It easy to show that (\ref{e8}) has
the first integral

\begin{equation}
\label{primera}f^{q-2}\dot f^2+\gamma_0 f=c
\end{equation}
and
\begin{equation}
\label{cerrada}x-x_0=\pm \int{\frac{f^{\frac{q}{2} -1}}{\sqrt{c-\gamma_0 f}
\, df}
\end{equation}

\noindent is the final expression for the general solution of (\ref{e8}).
Of course, $c$ and $x_0$ are arbitrary integration constants.

\subsection{The $q\neq 1$ environment}

>From (\ref{cerrada}) we see that $c>\gamma _0f$ in order to obtain a real
probability distribution $f$. This condition holds in several circumstances,
as follows:

\begin{itemize}
\item For $c>0$ and $\gamma_0>0$, the probability distribution
f<c/\gamma_0$, i.e., it has an upper positive limit and it is bounded.

\item For $c>0$ and $\gamma_0<0$, the probability distribution
f>c/\gamma_0$, and thus it has a negative lower limit and it is unbounded.

\item For $c<0$ and $\gamma_0>0$, the probability distribution
f<c/\gamma_0$, i.e., it is negatively defined and we must reject it.

\item For $c<0$ and $\gamma_0<0$, the probability distribution
f>c/\gamma_0$, i.e., it has a positive lower limit and it is unbounded.
\end{itemize}

 We conclude that the first possibility is the only one that has a
physical meaning. Indeed, inserting $f(x_m)=c/\gamma _0$ into (\ref{primera
), we obtain $\dot f=0$. Thus, at $x_m$ the probability distribution $f$ has
an extremum. On the other hand, inserting these results into (\ref{e8}) we
have $\ddot f=-\frac{\gamma _0}2\left[ \frac{\gamma _0}c\right] ^{q-2}$
(negatively defined). An $f$-maximum at $x_m$ ensues.

The particular case $c=0$ is consistent with $\gamma_0<0$ and needs separate
analysis. Here, the
 momoparametric
solution of (\ref{cerrada}) is given by

\begin{equation}
\label{c=0}f=\left[\pm\frac{(q-1)\sqrt{-\gamma_0}}{2}(x-x_0)\right]^{\frac{
}{q-1}}.
\end{equation}

This expression would yield  a probability distribution that cannot be
normalized if $x$ ranges from minus to plus infinity.

\subsection{Relativistic scenario}

Evaluating derivatives (with respect $x$) in  (\ref{e8}), and changing
variables to

\begin{equation}
\label{variable}u=\frac{\dot{f}(x)}{f(x)}\ ,
\end{equation}
we obtain
\begin{equation}
\label{e12}\ddot{u}\ +\alpha \ u\ \dot{u}\ +\beta \ u^3=0\ ,
\end{equation}
with
\begin{equation}
\label{4}\alpha \ =\ (2\ q-1)\ ,
\end{equation}
and
\begin{equation}
\label{5}\beta =-\frac 12\ q\ (1-q)\ .
\end{equation}

Equation (\ref{e12}) arises in several interesting physical problems as, for
instance, in the case of i) Einstein's field equations for homogeneous,
isotropic and spatially flat cosmological models with no cosmological constant
\cite{jakubi,1,2,3,jmp,viale,mar}, and ii) Bianchi I-type metrics \cite{8}, for
viscous fluids and a selfinteracting exponential potential scalar field.

We shall tackle (\ref{e12}) by recourse to ideas presented by Chimento and
Jakubi \cite{1}, who show that its more general solution that can be
{\it explicitly} writen as a function of $x$ obtains
 for the special value  $\alpha^2/\beta=9$. In
our case, this assertion translates itself into two possible
$q$-values, namely:

\begin{itemize}
\item $q=2$, and
\item $q=-1$,
\end{itemize}
as will be discussed below.

 Following \cite{1},
  we perform
 the nonlocal change
of variables

\begin{equation}
\label{e17}z\ =\ \frac 12\ u^2,
\end{equation}
and obtain
\begin{equation}
\label{e19}\frac{d^2z}{d\eta ^2}\ +\ \frac{dz}{d\eta }\ +\ \beta ^{^{\prime
}}\ z=\ 0,
\end{equation}

\noindent which is a linear ordinary differential equation, with

\begin{equation}
\label{e18}d\eta \ =\ \alpha \ u\ dx,
\end{equation}

\noindent
and

\begin{equation}
\label{e20}\beta ^{^{\prime }}\ =\ 2\ \frac{\beta}{\alpha^2}.
\end{equation}

The general solution of (\ref{e19}) is of the form

\begin{equation}
\label{e22}z(\eta )\ =\ A_1\ \mbox{e}^{\lambda^+\ \eta }\ +\ A_2\ \mbox{e
^{\lambda^-\ \eta },
\end{equation}

\noindent where

\begin{equation}
\lambda^{\pm}= -\frac{1}{2}\pm \frac{1}{2|\alpha|}
\end{equation}

\noindent are the roots of the pertinent characteristic equation, $A_1$, $A_2 $
being integration constants.

Going all the way back from $z$ to $f$ we
obtain the most general distribution function that verifies the
Frieden-Soffer tenets. However, as far as one knows at the present
stage \cite{1}, the road back can be traversed in {\it explicit,
analytical}
fashion  only for

\begin{equation}
\label{bet} \frac{\beta}{\alpha^2}\ =\ \frac{1}{9},
\end{equation}
which underlines the importance of the above mentioned special $q$-values
($2$\,\,and\,\,$-1$) .
In the next section we stress the connection between
the probability distribution and the cosmological expansion factor,
identifying both physical quantities via equations (\ref{e12}) and (\ref{e19}
).

\section{Cosmological applications}

Due to their nonlinear nature, exact solutions to Einstein's equations
cannot easily be obtained. However, one finds diverse problems of physical
interest where Einstein field equations (EFE) can be cast as particular
instances of a second order nonlinear differential equation of the type \cite
{chime}
\begin{equation}
\label{zz}\ddot y\,+\,\alpha g(y)\dot y\,+\,\beta g(y)\int dyg(y)\,+\,\gamma
g(y)\,=\,0,
\end{equation}
where $y=y(x)$ and $g(y)$ is a real function. $\alpha ,\beta ,\gamma $ are
constant parameters.

As examples of the preceding assertion, we may mention the case of EFE i)
for homogeneous, isotropic and spatially flat cosmological models with no
cosmological constant \cite{1,2,3,jmp,viale,mar}, ii) for a time decaying
cosmological constant \cite{7}, and iii) for Bianchi I-type metrics with a
variety of matter sources \cite{8}.

\subsection{The expansion rate of the Universe: different scenarios}
It is thought that quantum effects played a leading role in the early
Universe, as, for instance, particle production and vacuum polarization, that
arise within a quantal framework. It is known that these two phenomena can
be modeled in terms of a classical bulk viscosity \cite{12}. By recourse to
the so-called Extended Irreversible Thermodynamics \cite{13,14}, a
relativistic, second-order off-equilibrium approach, one can deal with a
spatially flat homogeneous and isotropic space-time, filled with a causal
viscous fluid and described by the Friedmann-Robertson-Walker line element
\cite{15}

\begin{equation}
\label{metric}ds^{2} = dt^{2} -
a^{2}(t)\left[dx^{2}_{1}+dx^{2}_{2}+dx^{2}_{3}\right],
\end{equation}

\no with the full version of the transport equation for the bulk viscous
pressure $\sigma$

\begin{equation}
\label{1.6}\sigma+\tau\dot\sigma=-3\zeta H-\frac{1}{2}\epsilon\tau\sigma
\left(3H+\frac{\dot\tau}{\tau}-\frac{\dot\zeta}{\zeta}- \frac{\dot T}{T
\right).
\end{equation}

\noindent where $a(t)$ is the scale factor, $H=\dot a/a$ is the expansion
rate, $\tau$ is the characteristic timescale for linear relaxational effects,
$\zeta$ is the linear bulk viscosity and $T$ is the equilibrium temperature.

In \cite{1} it was considered a causal viscous fluid whose equilibrium
pressure obeys a $\gamma$-law equation of state, while the transport equation
for the viscous pressure $\sigma$ was taken with $\epsilon =0$, which
corresponds to the so-called truncated theory. There, it was shown that $H$
satisfies a modified Painlev\'e-Ince equation that, in turn, can be cast into
the form of (\ref{zz}), i.e.,

\begin{equation}
\label{0}\ddot H+\alpha H\dot H+\beta H^3=0
\end{equation}

\noindent with $f(y)=y$ and $\gamma =0$ \cite{chime}. Comparing the latter
with (\ref{e12}) we get $u=H$, and, a fortiori $a(t)=f(x)$ (from (\ref{variable}))
. This result shows that
{\it the probability distribution can be identified
with the expansion parameter for the evolution of the universe by choosing
the time variable $t$ proportional to the information parameter $x$}.
This suggest that the solutions of (18) representing
probability distributions without physical meaning, becuase they cannot be
normalized, may be associated with the expanding (FRW) cosmological model. So
that, we have a more extended background where the (FS)'s demon along with
non-extensive scenario can be used.

Cosmological models with a viscous fluid source have been also studied using
the {\it full} causal Irreversible Thermodynamics \cite{jakubi,viale,mar,17}.
Relating the temperature $T$ to the energy density in the simplest possible
fashion that is able to guarantee a positive heat capacity, one finds that the
expansion rate $H$ verifies

\begin{equation}
\label{1.8}H\ddot H-\left(1+r\right)\dot H^2+A H^2\dot H-B H^4=0
\end{equation}

\noindent  where $r>0$, while $A$, $B$ and $\gamma_0$ are constants.
By recourse to  an
adequate choice of variables ($H=y^n$, with $n=-1/r$) in Eq.(\ref{1.8}), the
above equation becomes (again) identical to  (\ref{zz}) \cite{chime,jakubi},
i.e., we can write

\begin{equation}
\label{1.9}\ddot y+\alpha y^{n}\dot y+\beta y^{2n+1}=0
\end{equation}

\noindent with $f(y)=y^n$ and $\gamma =0$.

For a still different scenario consider now
 the early time evolution of a dissipative universe, in the relaxation
dominated regime. One is again lead here to an equation for the
expansion rate $H$ that can
be identified with (\ref{e12}) and (\ref{zz}) \cite{18}. Making once more the
identification $u=H$, Eq.(\ref{1.9}) turns into Eq. (\ref{e12}), so
that $a(t)=f(x)
, and the probability distribution can, as before, be identified with the
expansion factor in the same manner as we did above.

 Let us discuss now the case of a
 perfect fluid source, satisfying an equation of state $p=(\sigma-1)\rho
$ with cosmological constant. The $00$ Einstein equation

\begin{equation}
\label{2.0}H^2=\frac{1}{3}\frac{\rho_0}{a^{3\sigma}}-\frac{k}{a^2}+\frac
\Lambda}{3}
\end{equation}

\noindent transforms itself into Eq.(\ref{zz}), with $f(y)=y$, and

\begin{equation}
\label{2.1}\ddot{H}+(2+3\sigma)H\dot{H}+3\sigma H^3-\sigma\Lambda H=0
\end{equation}

\noindent provided we derive twice Eq.(\ref{2}). Thus, $a(t)=f(x)$ in a
universe with no cosmological constant.

Another interesting illustration of Eq. (\ref{zz}) can be obtained by
considering an anisotropic Universe, described by a Bianchi type I metric,

\begin{equation}
\label{5.1}ds^2 = e^{f(t)} \left(-dt^2+dz^2\right) +
G(t)\left(e^{p(t)}\,dx^2 +e^{-p(t)}\,dy^2\right)
\end{equation}

\noindent that is driven by a minimally coupled scalar field with an
exponential potential \cite{chime}. The Klein-Gordon equation for the scalar
field and the Einstein equations for the metrics are expressed in terms of
the semiconformal factor $G$ and of their derivatives, as follows \cite
{chime}

\begin{equation}
\label{Klein}\dot\phi=\frac{m}{G}-\frac{k}{2}\frac{\dot G}{G},
\end{equation}

\begin{equation}
\label{pression}\dot p=\frac{a}{G},
\end{equation}

\begin{equation}
\label{ef}e^f=\frac{\ddot G}{2GV},
\end{equation}

\begin{equation}
\label{G2}\frac{\ddot G}{G} - \frac{1}{2}\left(\frac{\dot G}{G}\right)^2 -
\frac{\dot G}{G} \dot{f} + \frac{1}{2}\dot p^2 = - \dot\phi^2,
\end{equation}

\noindent where $m$ and $a$ are arbitrary integration constants. The
solutions for this set of equations can be obtained if one is able to solve
the following equation for $G$

\begin{equation}
\label{first}G\frac{\ddot G}{\dot G}+(K-1)\dot G+\frac{M^2}{\dot G}=c,
\end{equation}

\noindent where $K$ is a constant parameter and $M$, $c$ are integration
constants. Once $G(t)$ is known, one can compute the field $\phi (t)$ and
the remaining metric functions $p(t)$, $f(t)$ from Eq.(\ref{Klein}), Eq.(\ref
{pression}), and Eq.(\ref{ef}), respectively. Making

\begin{equation}
\label{ansatz}G=y^{1/K}, \qquad and \qquad \tau\equiv -ct,
\end{equation}

\noindent Eq.(\ref{first}) becomes an equation of the type Eq.(\ref{zz})

\begin{equation}
\label{typeq}y^{\prime\prime}+y^{-n}y^{\prime}+\frac{M^2}{nc_3^2}
y^{1-2n}=0,
\end{equation}

\noindent where a prime denotes the derivative with respect to $\tau$ and
n=1/K$ \cite{8}.
\subsection{Time variable and information variable}

It should become clear by now that it is  of great interest to investigate
Eq.(\ref{zz}) from {\it all conceivable points of view}. Using nonlocal
transformations this equation can be first linearized and then solved with
the only proviso of the form invariance of Eq.(\ref{zz}) for an arbitrary
function $f(y)$.

Now, if we relate $y$ to the probability distribution function $f$ of the
former Sections in the fashion
\begin{equation}
y\,=\,\frac{\dot{f}}f,
\end{equation}
and set
\begin{equation}
H(t)\,=\,y,
\end{equation}

\noindent together with the choice: {\it time variable $t$ proportional to the
information parameter $x$}, one immediately realizes that equation(\ref{zz})
translates itself into equation (\ref{e12}), the protagonist of our
discussion in Section III. This shows that the probability distribution, in
many circumstances, can be identified with the scale factor $a(t)$ of
the universe.  In such a case, taking into account that the universe
has a final, matter dominated Friedmann stage, the physical range of
the Lagrange multiplier is given by $\gamma _0<0$. For $c>0$ the
solutions stem from a singularity.  However, for $c<0$ the solutions
avoid the initial singularity. In the particular case $c=0$ the
solutions (\ref{c=0}) represent a universe beginning or ending in a
singularity at $x=x_0$, for $q>1$. For $q<1$ the solutions diverge at
the finite proper time $x=x_0$.

The special mapping:  (information $\rightarrow\,\,$ time)
 that arises in natural
fashion from the present considerations allows one to give
 a degree of plasibility to the particular Frieden-Soffer solutions with $c=0$
(\ref{c=0}), since the scale factor, and its associated probability
density, represent
the evolution of universe.  In such circumstances,
 the scale factor is, of course,  unbounded,
it does not make any  sense to try normalizing it.

Summing up: we are leading to a variety of cosmological problems involving
Einstein's field equations which via a nomlocal transformation and the special
mapping:  (information $\rightarrow\,\,$ time) {\it can be reinterpreted} as a
problem within a non-extensive scenario, by application of the nonextensive
EPI Principle.

\section{ Conclusions}

In the present communication we have achieved the following results

\begin{itemize}
\item  We translated the rules of the Frieden-Soffer Information Transfer
game into a NET parlance.

\item  By suitably playing the game, we found exact analytical
solutions to special instances of Einstein's field equations.

\item  Nonlocal variable transformations of differential equations have
allowed for a widening of the class of differential equations (particularly
of the second order) which can be linearized and solved.
\end{itemize}

In essence our procedure provides a method  for determining a useful
nonlocal symmetry for a set of differential equations characteristic of
several different physical problems. However, these problem can be seen to
be are equivalent under nonlocal variable transformations. {\it In  that
peculiar sense}, we see that some  cosmological solutions  can be regarded
as probability distributions, pertaining to a $q\ne 1$ NET environment, that
maximize Fisher's information measure. One may think that Frieden-Soffer's
demon is able to map the space-time continuum into a probability space. The
solutions here studied can then be added to the impressive collection
elaborated by Frieden and Soffer \cite{f7}.

That a NET environment may be the appropriate one in order to address
physical problems connected with gravitation can be understood on the
grounds that two systems $A$ and $B$ that interact in such a fashion can
never be really ``separated" (and become isolated ones) as additivity would
demand, since they will always ``feel" each other's presence on account of
the infinite range of the associated force. Thus
\begin{equation}
S_{(A+B)}\,\neq\, S_A \,+\,S_B
\end{equation}
i.e., a non-extensive treatment becomes mandatory. Indeed, a classical
problem related to Galactic models that had defied solution within an
extensive Shannonian context has been successfully solved appealing to
Tsallis' NET \cite{t5}. The present results should be viewed within such a
background-landscape.

As a last comment of this paper we mention that the results concerning to the
equivalence of some expanding (FRW) cosmological model and the problem related
with the (FS)'s demon in a non-extensive scenario are not consequence of the
equivalence principle. They are a diret consequence of the nonlocal symmetries
that have certain nonlinear differential equations, in our case, the
Einstein's field equations for a spatially flat (FRW) universe and the
nonlinear differential equation (18), resulting for a simple variational
calculus of the Fisher's generalized information for translation families.
This nonlocal symmetries can be realized by means of nonlocal variable
transformation. They mapped different physical problem changing the number of
point symmetry generators and the very physics of the original problem.
However, the final motion equations are the same giving rise to a larger class
of phisical problem which can be considered "equivalent" under a most general
nonlocal symmetry group. In a forthcomming paper we shall investigate the
connection between the inertia of accelerative effects and the statistics.

\acknowledgements

The authors are indebted to the PROTEM Program of the National Research
Council (CONICET) of Argentina for financial support. The assistance of C.
Mostaccio is gratefully acknowledged. One of us (L.P.C) is indebted to the
University of Buenos Aires for financial support under Grant EX-260.