Storm classes

Classes

struct  ConjugateGradientOnly< Size, Precision >
 Class for performing optimization with Conjugate Gradient, where only the derivatives are available. More...
class  NullGraphics
 Graphics class which does absoloutely nothing. More...
class  DataForMCMC
 Closure hoding the data required do use GibbsSampler2 See FitSpots for naming of variables. More...
struct  SampledBackgroundData
 Closure holding image data generated using samples drawn from the model. More...
struct  SpotNegProbabilityDiffWithSampledBackground
 Compute the derivative of the negative log probability with respect to the parameters of one spot, given some samples of the other spots. More...
class  FreeEnergyHessian
 Class for computing the Hessian of the negative free energy. More...
class  FitSpots
 Mega class which actually does the meat of the spot fitting. More...
class  SampledMultispot::GibbsSampler
 Draw samples from the spot states given the spots positions and some data. More...
class  SampledMultispot::GibbsSampler2
 Gibbs sampling class which masks spots to reduce computation. More...

Functions

auto_ptr< FitSpotsGraphicsnull_graphics ()
void get_spot_pixels (const vector< ImageRef > &pixels, const Vector< 4 > &spot, vector< int > &out)
StateParameters parse_log_file (istream &in)
double brightness_motion_limit (double mu, double sigma, bool not_one)
void fit_spots_new (const vector< Image< float > > &ims, StateParameters &p, ofstream &save_spots, FitSpotsGraphics &gr)
template<class B >
double spot_shape_s (const TooN::Vector< 2 > &x, const TooN::Vector< 4, double, B > &phi)
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)
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)
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)
template<class B >
std::pair< double,
TooN::Vector< 2 > > 
spot_shape_diff (const TooN::Vector< 2 > &x, const TooN::Vector< 4, double, B > &phi)
template<class B >
double spot_shape (const TooN::Vector< 2 > &x, const TooN::Vector< 4, double, B > &phi)
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)
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)
template<class Base >
double log_probability_spot (const CVD::SubImage< float > &im, double variance, const TooN::Vector< 4, double, Base > &spot_parameters)
double log_normal_std (double mu, double sigma)
double log_normal_mode (double mu, double sigma)
double log_log_normal (double x, double mu, double sigma)
double diff_log_log_normal (double x, double mu, double sigma)
double hess_log_log_normal (double x, double mu, double sigma)

Function Documentation

auto_ptr<FitSpotsGraphics> null_graphics (  ) 

Factory function to generate an instance of NullGraphics.

Definition at line 86 of file multispot5.cc.

Referenced by mmain().

00087 {
00088     return auto_ptr<FitSpotsGraphics>(new NullGraphics);
00089 }

void get_spot_pixels ( const vector< ImageRef > &  pixels,
const Vector< 4 > &  spot,
vector< int > &  out 
)

Which pixels belong to a given spot? Find the indices of those pixels.

Definition at line 884 of file multispot5.cc.

Referenced by fit_spots_historic(), FitSpots::optimize_each_spot_in_turn_for_several_passes(), and FitSpots::try_modifying_model().

00885 {
00886     //Go out to three sigma
00887 
00888     vector<ImageRef> pix = getDisc(spot[1]*6 + 1);
00889     out.resize(0);
00890     ImageRef offset = ir_rounded(spot.slice<2,2>());
00891     for(unsigned int j=0; j < pix.size(); j++)
00892     {
00893         int pos = lower_bound(pixels.begin(), pixels.end(), pix[j] + offset) - pixels.begin();
00894         if(pos != (int)pixels.size() && pixels[pos] == pix[j] + offset)
00895             out.push_back(pos);
00896     }
00897 
00898     if(out.size() == 0)
00899     {
00900         cout << "********************************\n";
00901         cout << "********************************\n";
00902         cout << "********************************\n";
00903         cout << "********************************\n";
00904         cout << "********************************\n";
00905         cout << "Oe noes!11one\n";
00906         cout << pix.size() << endl;
00907     }
00908 }

StateParameters parse_log_file ( istream &  in  ) 

Parser for multispot 5 log files.

Log files are mostly line oriented and contain various records

The main records are:

Iteraton: #ITNUM MAIN: <list of spot parameters>

Pass: #PASSNUM MT19937 <random number generator state> PASS#PASSNUM: <list of spot parameters> ENDCHECKPOINT

Note that MAIN is redundant since it will be the same as the following PASS 1 (or the first pass computed if restoring from a checkpoint).

Data should only be considered valid after ENDCHECKPOINT has been read

Iteration is written once per iteration, not once per pass. (FIXME)

Which moron invented this file format???

Note that the file format hasn't beren fixed, so that the output can easily be compared to the output of the historic version which is known to be good.

Parameters:
in Stream to parse file from

Definition at line 984 of file multispot5.cc.

References StateParameters::iteration, StateParameters::pass, StateParameters::pixels, StateParameters::rng, split(), StateParameters::spots, and xtoa().

Referenced by mmain().

00985 {
00986     //A line read from the file
00987     string line;
00988 
00989     //State lines known to be OK
00990     string rngline, passline, iterationline;
00991     bool state_ok=0;
00992     
00993     //State lines read in, with flags of goodness
00994     string new_rngline, new_passline, new_iterationline;
00995     bool new_rngline_ok=0, new_passline_ok=0, new_iterationline_ok=0;
00996 
00997     unsigned int lineno=0;
00998     bool doing_gvars = 0;
00999 
01000     vector<ImageRef> pixels;
01001 
01002     while(!in.eof())
01003     {   
01004         getline(in, line);
01005         if(in.fail())
01006             break;
01007         
01008         lineno++;
01009         
01010         if(line == "ENDGVARLIST")
01011         {
01012             if(!doing_gvars)
01013                 throw LogFileParseError("Spurious end of GVars");
01014             doing_gvars = 0;
01015         }
01016         else if(doing_gvars)
01017         {
01018             GUI.ParseLine(line);
01019         }
01020         else if(line == "BEGINGVARLIST")
01021         {
01022             doing_gvars = 1;
01023         }
01024         if(line.substr(0, 11)  == "Iteration: ")
01025         {
01026             new_iterationline = line;
01027             new_iterationline_ok = true;
01028         }
01029         else if(line.substr(0, 4) == "PASS")
01030         {
01031             new_passline = line;
01032             if(new_passline_ok)
01033                 throw LogFileParseError("Duplicate PASS on line " + xtoa(lineno));
01034             new_passline_ok = true;
01035         }
01036         else if(line.substr(0, 8) == "MT19937 ")
01037         {
01038             new_rngline = line;
01039             if(new_rngline_ok)
01040                 throw LogFileParseError("Duplicate MT19937 on line " + xtoa(lineno));
01041             
01042             new_rngline_ok = true;
01043 
01044         }
01045         else if(line == "ENDCHECKPOINT")
01046         {
01047             if(new_passline_ok && new_rngline_ok && new_iterationline_ok)
01048             {
01049                 iterationline = new_iterationline;
01050                 rngline = new_rngline;
01051                 passline = new_passline;    
01052             }
01053             else
01054                 throw LogFileParseError("Reached checkpoint with missing data: "
01055                          "it=" + xtoa(new_iterationline_ok) + 
01056                         " pa=" + xtoa(new_passline_ok) + 
01057                         " rg=" + xtoa(new_rngline_ok) + " on line " + xtoa(lineno));
01058             
01059             //Don't reset iteration since it only appears once for each entire
01060             //set of passes. 
01061             new_rngline_ok = 0;
01062             new_passline_ok = 0;
01063 
01064             state_ok = true;
01065         }
01066         else if(line.substr(0, 7) == "PIXELS ")
01067         {
01068             vector<string> l = split(line);
01069             if( (l.size() - 1)%2 == 0)
01070             {
01071                 int n = (l.size()-1)/2;
01072                 pixels.resize(n);
01073                 for(int i=0; i < n; i++)
01074                 {
01075                     pixels[i].x = atox<int>(l[i*2+1+0], "pixels");
01076                     pixels[i].y = atox<int>(l[i*2+1+1], "pixels");
01077                 }
01078             }
01079             else
01080                 throw LogFileParseError("Bad PIXELS line");
01081         }
01082     }
01083 
01084     if(!state_ok)
01085         throw LogFileParseError("No state found");
01086     
01087     if(pixels.size() == 0)
01088         throw LogFileParseError("No pixels, or pixels is empty");
01089 
01090     //Now parse the lines
01091     StateParameters p;
01092     vector<string> l;
01093     
01094     //Parse the iterations
01095     l = split(iterationline);
01096     p.iteration = atox<int>(l[1], "iteration");
01097 
01098     //Parse the random number generator
01099     p.rng =shared_ptr<MT19937>(new MT19937);
01100     {
01101         istringstream rng_s(rngline);
01102         try{
01103             p.rng->read(rng_s);
01104         }
01105         catch(MT19937::ParseError p)
01106         {
01107             throw LogFileParseError("Error parsing MT19937");
01108         }
01109     }
01110 
01111     //Parse PASS and the listing of spots
01112     l = split(passline);
01113     if( (l.size() - 1 ) % 4 == 0)
01114     {
01115         p.pass = atox<int>(l[0].substr(4), "pass");
01116 
01117         for(unsigned int i=0; i < (l.size()-1)/4; i++)
01118         {
01119             cout << l[i*4+1+0] << endl;
01120             cout << l[i*4+1+1] << endl;
01121             cout << l[i*4+1+2] << endl;
01122             cout << l[i*4+1+3] << endl;
01123             p.spots.push_back(makeVector(
01124                         atox<double>(l[i*4+1+0], "spot"),
01125                         atox<double>(l[i*4+1+1], "spot"),
01126                         atox<double>(l[i*4+1+2], "spot"),
01127                         atox<double>(l[i*4+1+3], "spot")));
01128         }
01129 
01130     }
01131     else
01132         throw LogFileParseError("Wrong number of elements in PASS line");
01133 
01134     //Set up the pixels (oh god the pixels)
01135     p.pixels = pixels;
01136 
01137     return p;   
01138 }

double brightness_motion_limit ( double  mu,
double  sigma,
bool  not_one 
)

How far should steps in brightness be limited to?

Definition at line 1332 of file multispot5.cc.

References log_normal_std().

01333 {
01334     if(not_one)
01335         return log_normal_std(mu, sigma);
01336     else
01337         return 1;
01338 }

void fit_spots_new ( const vector< Image< float > > &  ims,
StateParameters p,
ofstream &  save_spots,
FitSpotsGraphics gr 
)

Wrapper function for using FitSpots.

Definition at line 2836 of file multispot5.cc.

References FitSpots::run().

Referenced by mmain().

02837 {
02838     FitSpots fit(ims, gr, p, save_spots);
02839     fit.run();
02840 }

template<class B >
double spot_shape_s ( const TooN::Vector< 2 > &  x,
const TooN::Vector< 4, double, B > &  phi 
) [inline]

See spot_shape().

Parameters:
x $\Vec{x}$
phi $\Vec{\phi}$
Returns:
$s(\Vec{x}, \Vec{\phi}) $

Definition at line 36 of file storm.h.

Referenced by spot_shape(), spot_shape_diff(), spot_shape_diff_position(), spot_shape_hess(), and spot_shape_hess_position().

00037 {
00038     return -norm_sq(x - phi.template slice<2,2>()) / (2*phi[1]*phi[1]);
00039 }

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 
) [inline]

Compute the spot shape and its derivative with respect to posision.

See also spot_shape()

Parameters:
x $\Vec{x}$
phi $\Vec{\phi}$

Definition at line 46 of file storm.h.

References spot_shape_s(), and sq().

Referenced by SampledMultispot::compute_spot_intensity_derivatives().

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 }

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 
) [inline]

Compute the spot shape and its Hessian with respect to posision.

See also spot_shape()

Parameters:
x $\Vec{x}$
phi $\Vec{\phi}$

Definition at line 68 of file storm.h.

References spot_shape_s(), and sq().

Referenced by SampledMultispot::compute_spot_intensity_hessian().

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 }

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 
) [inline]

Value of the spot, given the parameters and input location.

The spot is described by the following formula:

\[ \mu(\Vec{x}, \Vec{\phi}) = \frac{\phi_1}{\phi_2\sqrt(2\pi)} e^s, \]

where

\[ s = -\frac{(x_1 - \phi_3)^2 + (x_2 - \phi_4)^2}{2\phi_2^2}. \]

This describes a generic blobby spot function of a variable size. The light output can be tuned by varying $\phi_1$, and the level of blur can be changed independently by varying $\phi_2$. The derivative is:

\begin{eqnarray} \frac{\partial \mu}{\partial \phi_1} &=& \frac{1}{\phi_2\sqrt{2\pi}}e^s\\ \frac{\partial \mu}{\partial \phi_2} &=& -\frac{\phi_1}{\phi_2^2\sqrt{2\pi}}e^s(1 + 2s) \end{eqnarray}

And the hessian is:

\begin{eqnarray} \frac{\partial^2 \mu}{\partial \phi_1^2} &=& 0\\ \frac{\partial^2 \mu}{\partial\phi_1 \partial \phi_2} &=& -\frac{1}{\phi_2^2\sqrt{2\pi}}e^s(1 + 2s)\\ \frac{\partial^2 \mu}{\partial \phi_2^2} &=& \frac{2\phi_1}{\phi_2^3\sqrt{2\pi}}e^s(1 + 5s + 2s^2) \end{eqnarray}

Parameters:
x $\Vec{x}$
phi $\Vec{\phi}$
Returns:
$\mu(\Vec{x}, \Vec{\phi}) $

Definition at line 143 of file storm.h.

References spot_shape_s().

Referenced by log_probability_spot_hess().

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 }

template<class B >
std::pair<double, TooN::Vector<2> > spot_shape_diff ( const TooN::Vector< 2 > &  x,
const TooN::Vector< 4, double, B > &  phi 
) [inline]

see spot_shape_hess()

Parameters:
x $\Vec{x}$
phi $\Vec{\phi}$
Returns:
$\mu(\Vec{x}, \Vec{\phi}) $

Definition at line 165 of file storm.h.

References spot_shape_s().

Referenced by log_probability_spot_diff().

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 }

template<class B >
double spot_shape ( const TooN::Vector< 2 > &  x,
const TooN::Vector< 4, double, B > &  phi 
) [inline]

see spot_shape_hess()

Parameters:
x $\Vec{x}$
phi $\Vec{\phi}$
Returns:
$\mu(\Vec{x}, \Vec{\phi}) $

Definition at line 181 of file storm.h.

References spot_shape_s().

Referenced by SampledMultispot::SpotWithBackgroundMasked::hess_log(), and log_probability_spot().

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 }

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 
) [inline]

Find the log probability of an image patch, assuming zero base-line mean and the given variance.

This function makes use of the spot shape. It is assumed that the centre pixel of the image is at 0,0. Since the noise is Gaussian:

\begin{eqnarray} 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}} \\ \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}, \end{eqnarray}

where I is the image, and N is the number of pixels. See also log_probability_no_spot and $\mu$ (spot_shape). The derivatives are:

\begin{eqnarray} \frac{\partial \ln P(I)}{\partial \phi_0} &=& \frac{1}{\sigma^2} \sum_{\Vec{x}}(I_{\Vec{x}} - \mu(\Vec{x},\Vec{\phi})) \frac{\partial}{\partial \phi_0}\mu(\Vec{x}, \Vec{\phi})\\ \frac{\partial^2 \ln P(I)}{\partial \phi_0 \partial \phi_1} &=& \frac{1}{\sigma^2} \sum_{\Vec{x}}(I_{\Vec{x}} - \mu(\Vec{x},\Vec{\phi})) \frac{\partial^2}{\partial \phi_0 \partial \phi_1}\mu(\Vec{x}, \Vec{\phi}) - \frac{\partial}{\partial \phi_0}\mu(\Vec{x},\Vec{\phi}) \frac{\partial}{\partial \phi_1}\mu(\Vec{x},\Vec{\phi}) \end{eqnarray}

Parameters:
im Image
variance $\sigma^2$
spot_parameters $\Vec{\phi}$
Returns:
The log probability

Definition at line 236 of file storm.h.

References spot_shape_hess(), and sq().

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 }

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 
) [inline]

See log_probability_spot_hess.

Parameters:
im Image
variance $\sigma^2$
spot_parameters $\Vec{\phi}$
Returns:
The log probability

Definition at line 276 of file storm.h.

References spot_shape_diff(), and sq().

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 }

template<class Base >
double log_probability_spot ( const CVD::SubImage< float > &  im,
double  variance,
const TooN::Vector< 4, double, Base > &  spot_parameters 
) [inline]

See log_probability_spot_hess.

Parameters:
im Image
variance $\sigma^2$
spot_parameters $\Vec{\phi}$
Returns:
The log probability

Definition at line 311 of file storm.h.

References spot_shape(), and sq().

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 }

double log_normal_std ( double  mu,
double  sigma 
) [inline]

Compute the standard deviation of a log-normal distribution.

See log_normal().

\begin{equation} \mathrm{Var}[P(x)] = (e^(\sigma^2)-1)e^(2*\mu+\sigma^2) \end{equation}

Parameters:
sigma $ \sigma$
mu $ \mu$
Returns:
The standard deviation

Definition at line 343 of file storm.h.

References sq().

Referenced by brightness_motion_limit().

00344 {
00345     return sqrt((exp(sq(sigma)) - 1) * exp(2*mu + sq(sigma)));
00346 }

double log_normal_mode ( double  mu,
double  sigma 
) [inline]

Compute the mode of a log-normal distribution.

See log_normal().

\begin{equation} \mathrm{Mode}[P(x)] = e^(\mu-\sigma^2) \end{equation}

Parameters:
sigma $ \sigma$
mu $ \mu$
Returns:
The mode

Definition at line 359 of file storm.h.

Referenced by fit_spots_historic(), generate_state_parameters_ye_olde(), and FitSpots::try_modifying_model().

00360 {
00361     return exp(mu - sigma * sigma);
00362 }

double log_log_normal ( double  x,
double  mu,
double  sigma 
) [inline]

Log-normal distribution.

This is given by:

\begin{eqnarray} P(x) &=& \frac{1}{x\sigma\sqrt{2\pi}} e^{-\frac{(\ln x - \mu)^2}{s\sigma^2}}\\ \ln P(x) &=& -\frac{(\ln x - \mu)^2}{s\sigma^2} - \ln x - \ln\sigma\sqrt{2\pi}. \end{eqnarray}

Parameters:
x x
mu $\mu$
sigma $\sigma$

Definition at line 374 of file storm.h.

References ln(), and sq().

Referenced by NegativeFreeEnergy::compute_with_mask(), and NegativeFreeEnergy::operator()().

00375 {
00376     return -sq(ln(x) - mu) / (2*sq(sigma)) - ln(x) - ln(sigma * sqrt(2*M_PI));
00377 }

double diff_log_log_normal ( double  x,
double  mu,
double  sigma 
) [inline]

Derivative of the log of the log-normal distribution:

\[ \frac{\partial \ln P(x)}{\partial x} = -\frac{1}{x}\left(1 + \frac{\ln x - \mu}{\sigma^2}\right). \]

.

Parameters:
x x
mu $\mu$
sigma $\sigma$

Definition at line 388 of file storm.h.

References ln(), and sq().

Referenced by FreeEnergyHessian::hessian(), SpotNegProbabilityDiffWithSampledBackground::operator()(), sampled_background_spot_hessian2(), and sampled_background_spot_hessian_ffbs().

00389 {
00390     return -(1 + (ln(x) - mu)/sq(sigma)) / x;
00391 }

double hess_log_log_normal ( double  x,
double  mu,
double  sigma 
) [inline]

Second derivative of the log of the log-normal distribution:

\[ \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). \]

.

Parameters:
x x
mu $\mu$
sigma $\sigma$

Definition at line 403 of file storm.h.

References ln(), and sq().

Referenced by FreeEnergyHessian::hessian(), sampled_background_spot_hessian2(), sampled_background_spot_hessian_FAKE(), and sampled_background_spot_hessian_ffbs().

00404 {
00405     return (1 + (ln(x) - mu - 1)/sq(sigma)) / sq(x);
00406 }

Generated on Wed Nov 2 18:00:01 2011 for BCUBED by  doxygen 1.6.3