<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 3.2 Final//EN"> <!--Converted with LaTeX2HTML 98.1p1 release (March 2nd, 1998) originally by Nikos Drakos (nikos@cbl.leeds.ac.uk), CBLU, University of Leeds * revised and updated by: Marcus Hennecke, Ross Moore, Herb Swan * with significant contributions from: Jens Lippmann, Marek Rouchal, Martin Wilck and others --> <HTML> <HEAD> <TITLE>Introduction</TITLE> <META NAME="description" CONTENT="Introduction"> <META NAME="keywords" CONTENT="vol2"> <META NAME="resource-type" CONTENT="document"> <META NAME="distribution" CONTENT="global"> <META HTTP-EQUIV="Content-Type" CONTENT="text/html; charset=iso-8859-1"> <LINK REL="STYLESHEET" HREF="vol2.css"> <LINK REL="next" HREF="node337.html"> <LINK REL="previous" HREF="node335.html"> <LINK REL="up" HREF="node335.html"> <LINK REL="next" HREF="node337.html"> </HEAD> <BODY > <!--Navigation Panel--> <A NAME="tex2html5639" HREF="node337.html"> <IMG WIDTH="37" HEIGHT="24" ALIGN="BOTTOM" BORDER="0" ALT="next" SRC="icons.gif/next_motif.gif"></A> <A NAME="tex2html5636" HREF="node335.html"> <IMG WIDTH="26" HEIGHT="24" ALIGN="BOTTOM" BORDER="0" ALT="up" SRC="icons.gif/up_motif.gif"></A> <A NAME="tex2html5630" HREF="node335.html"> <IMG WIDTH="63" HEIGHT="24" ALIGN="BOTTOM" BORDER="0" ALT="previous" SRC="icons.gif/previous_motif.gif"></A> <A NAME="tex2html5638" HREF="node1.html"> <IMG WIDTH="65" HEIGHT="24" ALIGN="BOTTOM" BORDER="0" ALT="contents" SRC="icons.gif/contents_motif.gif"></A> <BR> <B> Next:</B> <A NAME="tex2html5640" HREF="node337.html">Regularization in the wavelet</A> <B> Up:</B> <A NAME="tex2html5637" HREF="node335.html">Deconvolution</A> <B> Previous:</B> <A NAME="tex2html5631" HREF="node335.html">Deconvolution</A> <BR> <BR> <!--End of Navigation Panel--> <H2><A NAME="SECTION002081000000000000000"> Introduction</A> </H2> Consider an image characterized by its intensity distribution <I>I</I>(<I>x</I>,<I>y</I>), corresponding to the observation of an object <I>O</I>(<I>x</I>,<I>y</I>) through an optical system. If the imaging system is linear and shift-invariant, the relation between the object and the image in the same coordinate frame is a convolution: <BR> <DIV ALIGN="CENTER"><A NAME="eqn_1"> </A> <!-- MATH: \begin{eqnarray} I(x,y)= O(x,y) * P(x,y) + N(x,y) \end{eqnarray} --> <TABLE ALIGN="CENTER" CELLPADDING="0" WIDTH="100%"> <TR VALIGN="MIDDLE"><TD NOWRAP ALIGN="RIGHT"><I>I</I>(<I>x</I>,<I>y</I>)= <I>O</I>(<I>x</I>,<I>y</I>) * <I>P</I>(<I>x</I>,<I>y</I>) + <I>N</I>(<I>x</I>,<I>y</I>)</TD> <TD> </TD> <TD> </TD> <TD WIDTH=10 ALIGN="RIGHT"> (14.99)</TD></TR> </TABLE></DIV> <BR CLEAR="ALL"><P></P><I>P</I>(<I>x</I>,<I>y</I>) is the point spread function (PSF) of the imaging system, and <I>N</I>(<I>x</I>,<I>y</I>) is an additive noise. In Fourier space we have: <BR> <DIV ALIGN="CENTER"> <!-- MATH: \begin{eqnarray} \hat I(u,v)= \hat O(u,v) \hat P(u,v) + \hat N(u,v) \end{eqnarray} --> <TABLE ALIGN="CENTER" CELLPADDING="0" WIDTH="100%"> <TR VALIGN="MIDDLE"><TD NOWRAP ALIGN="RIGHT"><IMG WIDTH="330" HEIGHT="52" ALIGN="MIDDLE" BORDER="0" SRC="img830.gif" ALT="$\displaystyle \hat I(u,v)= \hat O(u,v) \hat P(u,v) + \hat N(u,v)$"></TD> <TD> </TD> <TD> </TD> <TD WIDTH=10 ALIGN="RIGHT"> (14.100)</TD></TR> </TABLE></DIV> <BR CLEAR="ALL"><P></P> We want to determine <I>O</I>(<I>x</I>,<I>y</I>) knowing <I>I</I>(<I>x</I>,<I>y</I>) and <I>P</I>(<I>x</I>,<I>y</I>). This inverse problem has led to a large amount of work, the main difficulties being the existence of: (i) a cut-off frequency of the PSF, and (ii) an intensity noise (see for example [<A HREF="node370.html#cornwell">6</A>]). <P> Equation <A HREF="node336.html#eqn_1">14.99</A> is always an ill-posed problem. This means that there is not a unique least-squares solution of minimal norm <!-- MATH: $\parallel I(x,y) - P(x,y) * O(x,y) \parallel^2$ --> <IMG WIDTH="293" HEIGHT="48" ALIGN="MIDDLE" BORDER="0" SRC="img831.gif" ALT="$\parallel I(x,y) - P(x,y) * O(x,y) \parallel^2$"> and a regularization is necessary. <P> The best restoration algorithms are generally iterative [<A HREF="node370.html#katsaggelos">24</A>]. Van Cittert [<A HREF="node370.html#cittert">41</A>] proposed the following iteration: <BR> <DIV ALIGN="CENTER"><A NAME="vanvan"> </A> <!-- MATH: \begin{eqnarray} O^{(n+1)} (x,y) = O^{(n)} (x,y) + \alpha(I(x,y) - P(x,y)* O^{(n)} (x,y)) \end{eqnarray} --> <TABLE ALIGN="CENTER" CELLPADDING="0" WIDTH="100%"> <TR VALIGN="MIDDLE"><TD NOWRAP ALIGN="RIGHT"><IMG WIDTH="569" HEIGHT="53" ALIGN="MIDDLE" BORDER="0" SRC="img832.gif" ALT="$\displaystyle O^{(n+1)} (x,y) = O^{(n)} (x,y) + \alpha(I(x,y) - P(x,y)* O^{(n)} (x,y))$"></TD> <TD> </TD> <TD> </TD> <TD WIDTH=10 ALIGN="RIGHT"> (14.101)</TD></TR> </TABLE></DIV> <BR CLEAR="ALL"><P></P> where <IMG WIDTH="19" HEIGHT="21" ALIGN="BOTTOM" BORDER="0" SRC="img833.gif" ALT="$\alpha$"> is a converging parameter generally taken as 1. In this equation, the object distribution is modified by adding a term proportional to the residual. But this algorithm diverges when we have noise [<A HREF="node370.html#frieden">12</A>]. Another iterative algorithm is provided by the minimization of the norm <!-- MATH: $\parallel I(x,y) - P(x,y)* O(x,y) \parallel^2$ --> <IMG WIDTH="294" HEIGHT="48" ALIGN="MIDDLE" BORDER="0" SRC="img834.gif" ALT="$\parallel I(x,y) - P(x,y)* O(x,y) \parallel^2$">[<A HREF="node370.html#Landweber">21</A>] and leads to: <BR> <DIV ALIGN="CENTER"><A NAME="eqn_carre"> </A> <!-- MATH: \begin{eqnarray} O^{(n+1)} (x,y) = O^{(n)} (x,y) + \alpha P_s(x,y) * [I(x,y) - P(x,y) * O^{(n)} (x,y)] \end{eqnarray} --> <TABLE ALIGN="CENTER" CELLPADDING="0" WIDTH="100%"> <TR VALIGN="MIDDLE"><TD NOWRAP ALIGN="RIGHT"><IMG WIDTH="656" HEIGHT="53" ALIGN="MIDDLE" BORDER="0" SRC="img835.gif" ALT="$\displaystyle O^{(n+1)} (x,y) = O^{(n)} (x,y) + \alpha P_s(x,y) * [I(x,y) - P(x,y) * O^{(n)} (x,y)]$"></TD> <TD> </TD> <TD> </TD> <TD WIDTH=10 ALIGN="RIGHT"> (14.102)</TD></TR> </TABLE></DIV> <BR CLEAR="ALL"><P></P> where <!-- MATH: $P_s(x,y)=P(-x,-y)$ --> <I>P</I><SUB><I>s</I></SUB>(<I>x</I>,<I>y</I>)=<I>P</I>(-<I>x</I>,-<I>y</I>). <P> Tikhonov's regularization [<A HREF="node370.html#tikhonov">40</A>] consists of minimizing the term: <BR> <DIV ALIGN="CENTER"> <!-- MATH: \begin{eqnarray} \parallel I(x,y) - P(x,y)* O(x,y) \parallel^2 + \lambda \parallel H * O\parallel^2 \end{eqnarray} --> <TABLE ALIGN="CENTER" CELLPADDING="0" WIDTH="100%"> <TR VALIGN="MIDDLE"><TD NOWRAP ALIGN="RIGHT"><IMG WIDTH="438" HEIGHT="50" ALIGN="MIDDLE" BORDER="0" SRC="img836.gif" ALT="$\displaystyle \parallel I(x,y) - P(x,y)* O(x,y) \parallel^2 + \lambda \parallel H * O\parallel^2$"></TD> <TD> </TD> <TD> </TD> <TD WIDTH=10 ALIGN="RIGHT"> (14.103)</TD></TR> </TABLE></DIV> <BR CLEAR="ALL"><P></P> where <I>H</I> corresponds to a high-pass filter. This criterion contains two terms; the first one, <!-- MATH: $\parallel I(x,y) - P(x,y)* O(x,y) \parallel^2$ --> <IMG WIDTH="293" HEIGHT="48" ALIGN="MIDDLE" BORDER="0" SRC="img837.gif" ALT="$\parallel I(x,y) - P(x,y)* O(x,y) \parallel^2$">, expresses fidelity to the data <I>I</I>(<I>x</I>,<I>y</I>) and the second one, <!-- MATH: $\lambda \parallel H * O\parallel^2$ --> <IMG WIDTH="126" HEIGHT="48" ALIGN="MIDDLE" BORDER="0" SRC="img838.gif" ALT="$\lambda \parallel H * O\parallel^2$">, smoothness of the restored image. <IMG WIDTH="18" HEIGHT="22" ALIGN="BOTTOM" BORDER="0" SRC="img839.gif" ALT="$\lambda$"> is the regularization parameter and represents the trade-off between fidelity to the data and the restored image smoothness. Finding the optimal value <IMG WIDTH="19" HEIGHT="22" ALIGN="BOTTOM" BORDER="0" SRC="img840.gif" ALT="$\lambda$"> necessitates use of numeric techniques such as Cross-Validation [<A HREF="node370.html#golub">15</A>] [<A HREF="node370.html#galatsanos">14</A>]. <P> This method works well, but it is relatively long and produces smoothed images. This second point can be a real problem when we seek compact structures as is the case in astronomical imaging. An iterative approach for computing maximum likelihood estimates may be used. The Lucy method [#lucy<#15258<A HREF="node370.html#>"></A>,#katsaggelos<#15259<A HREF="node370.html#>"></A>,#adorf<#15260<A HREF="node370.html#>"></A>] uses such an iterative approach: <BR> <DIV ALIGN="CENTER"> <!-- MATH: \begin{eqnarray} O^{(n+1)} = O^{(n)} [ \frac{I}{I^{(n)}} \ast P^* ] \end{eqnarray} --> <TABLE ALIGN="CENTER" CELLPADDING="0" WIDTH="100%"> <TR VALIGN="MIDDLE"><TD NOWRAP ALIGN="RIGHT"><IMG WIDTH="239" HEIGHT="71" ALIGN="MIDDLE" BORDER="0" SRC="img841.gif" ALT="$\displaystyle O^{(n+1)} = O^{(n)} [ \frac{I}{I^{(n)}} \ast P^* ]$"></TD> <TD> </TD> <TD> </TD> <TD WIDTH=10 ALIGN="RIGHT"> (14.104)</TD></TR> </TABLE></DIV> <BR CLEAR="ALL"><P></P> and <BR> <DIV ALIGN="CENTER"> <!-- MATH: \begin{eqnarray} I^{(n)} = P \ast O^{(n)} \end{eqnarray} --> <TABLE ALIGN="CENTER" CELLPADDING="0" WIDTH="100%"> <TR VALIGN="MIDDLE"><TD NOWRAP ALIGN="RIGHT"><IMG WIDTH="150" HEIGHT="28" ALIGN="BOTTOM" BORDER="0" SRC="img842.gif" ALT="$\displaystyle I^{(n)} = P \ast O^{(n)}$"></TD> <TD> </TD> <TD> </TD> <TD WIDTH=10 ALIGN="RIGHT"> (14.105)</TD></TR> </TABLE></DIV> <BR CLEAR="ALL"><P></P> where <I>P</I><SUP>*</SUP> is the conjugate of the PSF. <P> <HR> <!--Navigation Panel--> <A NAME="tex2html5639" HREF="node337.html"> <IMG WIDTH="37" HEIGHT="24" ALIGN="BOTTOM" BORDER="0" ALT="next" SRC="icons.gif/next_motif.gif"></A> <A NAME="tex2html5636" HREF="node335.html"> <IMG WIDTH="26" HEIGHT="24" ALIGN="BOTTOM" BORDER="0" ALT="up" SRC="icons.gif/up_motif.gif"></A> <A NAME="tex2html5630" HREF="node335.html"> <IMG WIDTH="63" HEIGHT="24" ALIGN="BOTTOM" BORDER="0" ALT="previous" SRC="icons.gif/previous_motif.gif"></A> <A NAME="tex2html5638" HREF="node1.html"> <IMG WIDTH="65" HEIGHT="24" ALIGN="BOTTOM" BORDER="0" ALT="contents" SRC="icons.gif/contents_motif.gif"></A> <BR> <B> Next:</B> <A NAME="tex2html5640" HREF="node337.html">Regularization in the wavelet</A> <B> Up:</B> <A NAME="tex2html5637" HREF="node335.html">Deconvolution</A> <B> Previous:</B> <A NAME="tex2html5631" HREF="node335.html">Deconvolution</A> <!--End of Navigation Panel--> <ADDRESS> <I>Petra Nass</I> <BR><I>1999-06-15</I> </ADDRESS> </BODY> </HTML>