ThreeB 1.1

- Main Page
- Related Pages
- Modules
- Namespaces
- Classes
- Files

- File List
- File Members

storm.h
Go to the documentation of this file.

    00001 /*
    00002     This file is part of B-cubed.
    00003
    00004     Copyright (C) 2009, 2010, 2011, Edward Rosten and Susan Cox
    00005
    00006     B-cubed is free software; you can redistribute it and/or
    00007     modify it under the terms of the GNU Lesser General Public
    00008     License as published by the Free Software Foundation; either
    00009     version 3.0 of the License, or (at your option) any later version.
    00010
    00011     This library is distributed in the hope that it will be useful,
    00012     but WITHOUT ANY WARRANTY; without even the implied warranty of
    00013     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
    00014     Lesser General Public License for more details.
    00015
    00016     You should have received a copy of the GNU General Public License
    00017     along with this program.  If not, see <http://www.gnu.org/licenses/>
    00018 */
    00019
    00020 #ifndef INC_STORM_H
    00021 #define INC_STORM_H
    00022
    00023 #include <TooN/TooN.h>
    00024 #include <cvd/image.h>
    00025 #include <utility>
    00026 #include <tr1/tuple>
    00027 #include "utility.h"
    00028
    00029
    00030 /**See spot_shape()
    00031 @param x \f$\Vec{x}\f$
    00032 @param phi \f$\Vec{\phi}\f$
    00033 @return \f$s(\Vec{x}, \Vec{\phi}) \f$
    00034 @ingroup gStorm
    00035 */
    00036 template<class B> double spot_shape_s(const TooN::Vector<2>& x, const TooN::Vector<4, double, B>& phi)
    00037 {
    00038     return -norm_sq(x - phi.template slice<2,2>()) / (2*phi[1]*phi[1]);
    00039 }
    00040
    00041 /** Compute the spot shape and its derivative with respect to posision. See also spot_shape()
    00042 @param x \f$\Vec{x}\f$
    00043 @param phi \f$\Vec{\phi}\f$
    00044 @ingroup gStorm
    00045 */
    00046 template<class B> std::pair<double, TooN::Vector<4> > spot_shape_diff_position(const TooN::Vector<2>& x, const TooN::Vector<4, double, B>& phi)
    00047 {
    00048     using namespace TooN;
    00049
    00050     double s = spot_shape_s(x, phi);
    00051     double r_2_pi = sqrt(2*M_PI);
    00052
    00053     double prob = exp(s) * phi[0]/(phi[1]*r_2_pi);
    00054
    00055     Vector<4> deriv = (exp(s) / (phi[1]*r_2_pi)) *
    00056                            makeVector(1,
    00057                                       -phi[0] * (1 + 2*s)/phi[1],
    00058                                       (x[0] - phi[2])*(phi[0]/sq(phi[1])),
    00059                                       (x[1] - phi[3])*(phi[0]/sq(phi[1])));
    00060     return std::make_pair(prob, deriv);
    00061 }
    00062
    00063 /** Compute the spot shape and its Hessian with respect to posision. See also spot_shape()
    00064 @param x \f$\Vec{x}\f$
    00065 @param phi \f$\Vec{\phi}\f$
    00066 @ingroup gStorm
    00067 */
    00068 template<class B> std::tr1::tuple<double, TooN::Vector<4>, TooN::Matrix<4> > spot_shape_hess_position(const TooN::Vector<2>& x, const TooN::Vector<4, double, B>& phi)
    00069 {
    00070     using namespace TooN;
    00071     using namespace std::tr1;
    00072
    00073     double s = spot_shape_s(x, phi);
    00074     double r_2_pi = sqrt(2*M_PI);
    00075
    00076     double es = exp(s);
    00077
    00078     double prob = es * phi[0]/(phi[1]*r_2_pi);
    00079
    00080     Vector<4> deriv = (es / (phi[1]*r_2_pi)) *
    00081                            makeVector(1,
    00082                                       -phi[0] * (1 + 2*s)/phi[1],
    00083                                       (x[0] - phi[2])*(phi[0]/sq(phi[1])),
    00084                                       (x[1] - phi[3])*(phi[0]/sq(phi[1])));
    00085
    00086     Matrix<4> hess;
    00087     hess[0][0] = 0;
    00088
    00089     hess[0][1] = -es*(1+2*s) / (phi[1] * phi[1] * r_2_pi);
    00090     hess[1][0] = hess[0][1];
    00091
    00092     hess[0][2] = es * (x[0] - phi[2]) / (pow(phi[1], 3)*r_2_pi);
    00093     hess[2][0] = es * (x[0] - phi[2]) / (pow(phi[1], 3)*r_2_pi);
    00094
    00095     hess[0][3] = es * (x[1] - phi[3]) / (pow(phi[1], 3)*r_2_pi);
    00096     hess[3][0] = es * (x[1] - phi[3]) / (pow(phi[1], 3)*r_2_pi);
    00097
    00098     hess[1][1] = 2*phi[0]*es*(1 + 5*s + 2*s*s) / ( pow(phi[1], 3) * r_2_pi);
    00099
    00100     hess[1][2] = -phi[0] * es * (3 + 2*s) * (x[0] - phi[2]) / (pow(phi[1], 4) * r_2_pi);
    00101     hess[1][3] = -phi[0] * es * (3 + 2*s) * (x[1] - phi[3]) / (pow(phi[1], 4) * r_2_pi);
    00102
    00103     hess[2][1] = hess[1][2];
    00104     hess[3][1] = hess[1][3];
    00105
    00106     hess[2][2] = phi[0] * es * (sq(x[0] - phi[2]) - sq(phi[1])) / (r_2_pi * pow(phi[1], 5));
    00107     hess[3][3] = phi[0] * es * (sq(x[1] - phi[3]) - sq(phi[1])) / (r_2_pi * pow(phi[1], 5));
    00108
    00109     hess[2][3] = phi[0] * es * (x[0] - phi[2])*(x[1] - phi[3]) / (r_2_pi * pow(phi[1], 5));
    00110     hess[3][2] = hess[2][3];
    00111
    00112
    00113     return make_tuple(prob, deriv, hess);
    00114 }
    00115
    00116 /**Value of the spot, given the parameters and input location.
    00117 The spot is described by the following formula:
    00118 \f[
    00119     \mu(\Vec{x}, \Vec{\phi}) = \frac{\phi_1}{\phi_2\sqrt(2\pi)} e^s,
    00120 \f]
    00121 where
    00122 \f[
    00123     s = -\frac{(x_1 - \phi_3)^2 + (x_2 - \phi_4)^2}{2\phi_2^2}.
    00124 \f]
    00125 This describes a generic blobby spot function of a variable size. The light output
    00126 can be tuned by varying \f$\phi_1\f$, and the level of blur can be changed independently
    00127 by varying \f$\phi_2\f$. The derivative is:
    00128 \f{eqnarray}{
    00129     \frac{\partial \mu}{\partial \phi_1} &=& \frac{1}{\phi_2\sqrt{2\pi}}e^s\\
    00130     \frac{\partial \mu}{\partial \phi_2} &=& -\frac{\phi_1}{\phi_2^2\sqrt{2\pi}}e^s(1 + 2s)
    00131 \f}
    00132 And the hessian is:
    00133 \f{eqnarray}{
    00134     \frac{\partial^2 \mu}{\partial \phi_1^2} &=& 0\\
    00135     \frac{\partial^2 \mu}{\partial\phi_1 \partial \phi_2} &=& -\frac{1}{\phi_2^2\sqrt{2\pi}}e^s(1 + 2s)\\
    00136     \frac{\partial^2 \mu}{\partial \phi_2^2} &=& \frac{2\phi_1}{\phi_2^3\sqrt{2\pi}}e^s(1 + 5s + 2s^2)
    00137 \f}
    00138 @param x \f$\Vec{x}\f$
    00139 @param phi \f$\Vec{\phi}\f$
    00140 @return \f$\mu(\Vec{x}, \Vec{\phi}) \f$
    00141 @ingroup gStorm
    00142 */
    00143 template<class B> std::tr1::tuple<double, TooN::Vector<2>, TooN::Matrix<2> > spot_shape_hess(const TooN::Vector<2>& x, const TooN::Vector<4, double, B>& phi)
    00144 {
    00145     double s = spot_shape_s(x, phi);
    00146     double r_2_pi = sqrt(2*M_PI);
    00147
    00148     double prob = exp(s) * phi[0]/(phi[1]*r_2_pi);
    00149     TooN::Vector<2> deriv = (exp(s) / (phi[1]*r_2_pi)) * TooN::makeVector(1, -phi[0] * (1 + 2*s)/phi[1]);
    00150     TooN::Matrix<2> hess;
    00151
    00152     hess[0][0] = 0;
    00153     hess[0][1] = -exp(s)*(1+2*s) / (phi[1] * phi[1] * r_2_pi);
    00154     hess[1][0] = hess[0][1];
    00155     hess[1][1] = 2*phi[0]*exp(s)*(1 + 5*s + 2*s*s) / ( pow(phi[1], 3) * r_2_pi);
    00156
    00157     return std::tr1::make_tuple(prob, deriv, hess);
    00158 }
    00159 /** see spot_shape_hess()
    00160 @param x \f$\Vec{x}\f$
    00161 @param phi \f$\Vec{\phi}\f$
    00162 @return \f$\mu(\Vec{x}, \Vec{\phi}) \f$
    00163 @ingroup gStorm
    00164 */
    00165 template<class B> std::pair<double, TooN::Vector<2> > spot_shape_diff(const TooN::Vector<2>& x, const TooN::Vector<4, double, B>& phi)
    00166 {
    00167     double s = spot_shape_s(x, phi);
    00168     double r_2_pi = sqrt(2*M_PI);
    00169
    00170     double prob = exp(s) * phi[0]/(phi[1]*r_2_pi);
    00171     TooN::Vector<2> deriv = (exp(s) / (phi[1]*r_2_pi)) * TooN::makeVector(1, -phi[0] * (1 + 2*s)/phi[1]);
    00172     return std::make_pair(prob, deriv);
    00173 }
    00174
    00175 /** see spot_shape_hess()
    00176 @param x \f$\Vec{x}\f$
    00177 @param phi \f$\Vec{\phi}\f$
    00178 @return \f$\mu(\Vec{x}, \Vec{\phi}) \f$
    00179 @ingroup gStorm
    00180 */
    00181 template<class B> double spot_shape(const TooN::Vector<2>& x, const TooN::Vector<4, double, B>& phi)
    00182 {
    00183     double s = spot_shape_s(x, phi);
    00184     double r_2_pi = sqrt(2*M_PI);
    00185
    00186     // FIXME FIXME FIXME and don't forget to fix the HESSIAN AND DERIVATIVE
    00187     // Should be:              1/(2 pi s^2)  for two dimensions
    00188     //                             vvvvvvvvvvvvv    http://lol.i.trollyou.com/
    00189     double prob = exp(s) * phi[0]/(phi[1]*r_2_pi);
    00190
    00191
    00192     return prob;
    00193 }
    00194
    00195 /**Find the log probability of an image patch,
    00196 assuming zero mean and the given variance, and no spot present.
    00197 See also log_probability_spot()
    00198 @param im Image
    00199 @param variance variance
    00200 @returns The log probability
    00201 */
    00202 inline double log_probability_no_spot(const CVD::SubImage<float>& im, double variance)
    00203 {
    00204     double logprob_part=0;
    00205     for(int y=0; y < im.size().y; y++)
    00206         for(int x=0; x < im.size().x; x++)
    00207             logprob_part -= im[y][x] * im[y][x];
    00208     return logprob_part/(2*variance) - im.size().area() * log(2*M_PI*variance)/2;
    00209
    00210 }
    00211
    00212 /**Find the log probability of an image patch, assuming zero base-line mean and
    00213 the given variance. This function makes use of the spot shape. It is assumed that the
    00214 centre pixel of the image is at 0,0. Since the noise is Gaussian:
    00215 \f{eqnarray}{
    00216     P(\text{image})    &=& \prod_{\Vec{x} \in \text{pixels}} \frac{1}{\sqrt{2\pi\sigma^2}}e^{-\frac{(I(\Vec{x}) - \mu(\Vec{x}, \Vec{\phi}))^2}{2\sigma^2}} \\
    00217     \ln P(\text{image}) &=& \sum_{\Vec{x} \in \text{pixels}}  -\frac{(I(\Vec{x}) - \mu(\Vec{x}, \Vec{\phi}))^2}{2\sigma^2} - \frac{N}{2} \ln {2 \pi \sigma^2},
    00218 \f}
    00219 where \e I is the image, and \e N is the number of pixels. See also ::log_probability_no_spot and \f$\mu\f$ (::spot_shape).
    00220 The derivatives are:
    00221 \f{eqnarray}{
    00222     \frac{\partial \ln P(I)}{\partial \phi_0} &=& \frac{1}{\sigma^2} \sum_{\Vec{x}}(I_{\Vec{x}} - \mu(\Vec{x},\Vec{\phi}))
    00223                                                            \frac{\partial}{\partial \phi_0}\mu(\Vec{x}, \Vec{\phi})\\
    00224     \frac{\partial^2 \ln P(I)}{\partial \phi_0 \partial \phi_1} &=&
    00225               \frac{1}{\sigma^2} \sum_{\Vec{x}}(I_{\Vec{x}} - \mu(\Vec{x},\Vec{\phi}))
    00226                                \frac{\partial^2}{\partial \phi_0 \partial \phi_1}\mu(\Vec{x}, \Vec{\phi}) -
    00227                                \frac{\partial}{\partial \phi_0}\mu(\Vec{x},\Vec{\phi})
    00228                                \frac{\partial}{\partial \phi_1}\mu(\Vec{x},\Vec{\phi})
    00229 \f}
    00230 @ingroup gStorm
    00231 @param im Image
    00232 @param variance \f$\sigma^2\f$
    00233 @param spot_parameters \f$\Vec{\phi}\f$
    00234 @returns The log probability
    00235 */
    00236 template<class Base> std::tr1::tuple<double, TooN::Vector<2>, TooN::Matrix<2> > log_probability_spot_hess(const CVD::SubImage<float>& im, double variance, const TooN::Vector<4, double, Base>& spot_parameters)
    00237 {
    00238     using namespace TooN;
    00239     using namespace std::tr1;
    00240
    00241     //-1 because if the image is 3x3, ie 0,1,2 then 1,1 is the centre.
    00242     //If it is 2x2, ie 0,1 then .5,.5 is the centre
    00243     Vector<2> centre = makeVector((im.size().x-1) / 2.0, (im.size().y-1) / 2.0);
    00244
    00245     double logprob_part=0;
    00246     Vector<2> diff = Zeros;
    00247     Matrix<2> hess = Zeros;
    00248     for(int y=0; y < im.size().y; y++)
    00249         for(int x=0; x < im.size().x; x++)
    00250         {
    00251             Vector<2> d = TooN::makeVector(x, y) - centre;
    00252
    00253             double mu;
    00254             Vector<2> diff_mu;
    00255             Matrix<2> hess_mu;
    00256             tie(mu, diff_mu, hess_mu) = spot_shape_hess(d, spot_parameters);
    00257
    00258             double e = im[y][x] - mu;
    00259
    00260             logprob_part += -sq(e);
    00261             diff         += diff_mu * e;
    00262             hess         += e * hess_mu - diff_mu.as_col() * diff_mu.as_row();
    00263         }
    00264     return make_tuple(   logprob_part / (2*variance) - im.size().area() * log(2*M_PI*variance)/2,
    00265                         diff / variance,
    00266                         hess / variance);
    00267 }
    00268
    00269 /** See log_probability_spot_hess
    00270 @ingroup gStorm
    00271 @param im Image
    00272 @param variance \f$\sigma^2\f$
    00273 @param spot_parameters \f$\Vec{\phi}\f$
    00274 @returns The log probability
    00275 */
    00276 template<class Base> std::pair<double, TooN::Vector<2> > log_probability_spot_diff(const CVD::SubImage<float>& im, double variance, const TooN::Vector<4, double, Base>& spot_parameters)
    00277 {
    00278     using namespace TooN;
    00279     using namespace std::tr1;
    00280     using namespace std;
    00281     //-1 because if the image is 3x3, ie 0,1,2 then 1,1 is the centre.
    00282     //If it is 2x2, ie 0,1 then .5,.5 is the centre
    00283     Vector<2> centre = makeVector((im.size().x-1) / 2.0, (im.size().y-1) / 2.0);
    00284
    00285     double logprob_part=0;
    00286     Vector<2> diff = Zeros;
    00287     for(int y=0; y < im.size().y; y++)
    00288         for(int x=0; x < im.size().x; x++)
    00289         {
    00290             Vector<2> d = makeVector(x, y) - centre;
    00291
    00292             double mu;
    00293             Vector<2> diff_mu;
    00294             tie(mu, diff_mu) = spot_shape_diff(d, spot_parameters);
    00295
    00296             double e = im[y][x] - mu;
    00297
    00298             logprob_part += -sq(e);
    00299             diff         += diff_mu * e;
    00300         }
    00301     return make_pair(logprob_part / (2*variance) - im.size().area() * log(2*M_PI*variance)/2, diff / variance);
    00302 }
    00303
    00304 /** See log_probability_spot_hess
    00305 @ingroup gStorm
    00306 @param im Image
    00307 @param variance \f$\sigma^2\f$
    00308 @param spot_parameters \f$\Vec{\phi}\f$
    00309 @returns The log probability
    00310 */
    00311 template<class Base> double log_probability_spot(const CVD::SubImage<float>& im, double variance, const TooN::Vector<4, double, Base>& spot_parameters)
    00312 {
    00313     //-1 because if the image is 3x3, ie 0,1,2 then 1,1 is the centre.
    00314     //If it is 2x2, ie 0,1 then .5,.5 is the centre
    00315     TooN::Vector<2> centre = TooN::makeVector((im.size().x-1) / 2.0, (im.size().y-1) / 2.0);
    00316
    00317     double logprob_part=0;
    00318     for(int y=0; y < im.size().y; y++)
    00319         for(int x=0; x < im.size().x; x++)
    00320         {
    00321             TooN::Vector<2> d = TooN::makeVector(x, y) - centre;
    00322
    00323             double mu = spot_shape(d, spot_parameters);
    00324
    00325             double e = im[y][x] - mu;
    00326
    00327             logprob_part += -sq(e);
    00328         }
    00329     return logprob_part / (2*variance) - im.size().area() * log(2*M_PI*variance)/2;
    00330 }
    00331
    00332 /**Compute the standard deviation of a log-normal distribution.
    00333
    00334 See log_normal().
    00335 \f{equation}
    00336     \mathrm{Var}[P(x)] = (e^(\sigma^2)-1)e^(2*\mu+\sigma^2)
    00337 \f}
    00338 @param sigma \f$ \sigma\f$
    00339 @param mu \f$ \mu\f$
    00340 @returns The standard deviation
    00341 @ingroup gStorm
    00342 */
    00343 inline double log_normal_std(double mu, double sigma)
    00344 {
    00345     return sqrt((exp(sq(sigma)) - 1) * exp(2*mu + sq(sigma)));
    00346 }
    00347
    00348 /**Compute the mode of a log-normal distribution.
    00349
    00350 See log_normal().
    00351 \f{equation}
    00352     \mathrm{Mode}[P(x)] = e^(\mu-\sigma^2)
    00353 \f}
    00354 @param sigma \f$ \sigma\f$
    00355 @param mu \f$ \mu\f$
    00356 @returns The mode
    00357 @ingroup gStorm
    00358 */
    00359 inline double log_normal_mode(double mu, double sigma)
    00360 {
    00361     return exp(mu - sigma * sigma);
    00362 }
    00363
    00364 /**Log-normal distribution. This is given by:
    00365 \f{eqnarray}{
    00366     P(x)     &=& \frac{1}{x\sigma\sqrt{2\pi}} e^{-\frac{(\ln x - \mu)^2}{s\sigma^2}}\\
    00367     \ln P(x) &=& -\frac{(\ln x - \mu)^2}{s\sigma^2} - \ln x - \ln\sigma\sqrt{2\pi}.
    00368 \f}
    00369 @param x \e x
    00370 @param mu \f$\mu\f$
    00371 @param sigma \f$\sigma\f$
    00372 @ingroup gStorm
    00373 */
    00374 inline double log_log_normal(double x, double mu, double sigma)
    00375 {
    00376     return -sq(ln(x) - mu) / (2*sq(sigma)) - ln(x) - ln(sigma * sqrt(2*M_PI));
    00377 }
    00378
    00379 /**Derivative of the log of the log-normal distribution:
    00380 \f[
    00381     \frac{\partial \ln P(x)}{\partial x} = -\frac{1}{x}\left(1 + \frac{\ln x - \mu}{\sigma^2}\right).
    00382 \f]
    00383 @param x \e x
    00384 @param mu \f$\mu\f$
    00385 @param sigma \f$\sigma\f$
    00386 @ingroup gStorm
    00387 */
    00388 inline double diff_log_log_normal(double x, double mu, double sigma)
    00389 {
    00390     return -(1 + (ln(x) - mu)/sq(sigma)) / x;
    00391 }
    00392
    00393
    00394 /**Second derivative of the log of the log-normal distribution:
    00395 \f[
    00396     \frac{\partial^2 \ln P(x)}{\partial x^2} = \frac{1}{x^2}\left(1 + \frac{\ln x - \mu}{\sigma^2} - \frac{1}{\sigma^2}\right).
    00397 \f]
    00398 @param x \e x
    00399 @param mu \f$\mu\f$
    00400 @param sigma \f$\sigma\f$
    00401 @ingroup gStorm
    00402 */
    00403 inline double hess_log_log_normal(double x, double mu, double sigma)
    00404 {
    00405     return (1 + (ln(x) - mu - 1)/sq(sigma)) / sq(x);
    00406 }
    00407
    00408
    00409 #endif