\hypertarget{multispot5_8cc}{
\subsection{multispot5.cc File Reference}
\label{multispot5_8cc}\index{multispot5.cc@{multispot5.cc}}
}


Fit spots to the data.  


{\ttfamily \#include $<$cstdlib$>$}\par
{\ttfamily \#include $<$cerrno$>$}\par
{\ttfamily \#include $<$cstring$>$}\par
{\ttfamily \#include $<$stack$>$}\par
{\ttfamily \#include $<$algorithm$>$}\par
{\ttfamily \#include $<$climits$>$}\par
{\ttfamily \#include $<$iomanip$>$}\par
{\ttfamily \#include $<$map$>$}\par
{\ttfamily \#include $<$tr1/memory$>$}\par
{\ttfamily \#include $<$cvd/image\_\-io.h$>$}\par
{\ttfamily \#include $<$cvd/image\_\-convert.h$>$}\par
{\ttfamily \#include $<$cvd/morphology.h$>$}\par
{\ttfamily \#include $<$cvd/connected\_\-components.h$>$}\par
{\ttfamily \#include $<$cvd/draw.h$>$}\par
{\ttfamily \#include $<$cvd/vector\_\-image\_\-ref.h$>$}\par
{\ttfamily \#include $<$cvd/random.h$>$}\par
{\ttfamily \#include $<$cvd/timer.h$>$}\par
{\ttfamily \#include $<$gvars3/instances.h$>$}\par
{\ttfamily \#include $<$TooN/functions/derivatives.h$>$}\par
{\ttfamily \#include $<$TooN/determinant.h$>$}\par
{\ttfamily \#include $<$TooN/SymEigen.h$>$}\par
{\ttfamily \#include $<$TooN/optimization/conjugate\_\-gradient.h$>$}\par
{\ttfamily \#include \char`\"{}conjugate\_\-gradient\_\-only.h\char`\"{}}\par
{\ttfamily \#include \char`\"{}forward\_\-algorithm.h\char`\"{}}\par
{\ttfamily \#include \char`\"{}numerical\_\-derivatives.h\char`\"{}}\par
{\ttfamily \#include \char`\"{}storm.h\char`\"{}}\par
{\ttfamily \#include \char`\"{}storm\_\-imagery.h\char`\"{}}\par
{\ttfamily \#include \char`\"{}debug.h\char`\"{}}\par
{\ttfamily \#include \char`\"{}sampled\_\-multispot.h\char`\"{}}\par
{\ttfamily \#include \char`\"{}mt19937.h\char`\"{}}\par
{\ttfamily \#include \char`\"{}utility.h\char`\"{}}\par
{\ttfamily \#include \char`\"{}multispot5.h\char`\"{}}\par
\subsubsection*{Classes}
\begin{DoxyCompactItemize}
\item 
class \hyperlink{classNullUICallback}{NullUICallback}
\begin{DoxyCompactList}\small\item\em User interface callback class which does nothing. \end{DoxyCompactList}\item 
class \hyperlink{classNullGraphics}{NullGraphics}
\begin{DoxyCompactList}\small\item\em Graphics class which does absoloutely nothing. \end{DoxyCompactList}\item 
class \hyperlink{classDataForMCMC}{DataForMCMC}
\begin{DoxyCompactList}\small\item\em Closure hoding the data required do use GibbsSampler2 See \hyperlink{classFitSpots}{FitSpots} for naming of variables. \end{DoxyCompactList}\item 
class \hyperlink{classKahan}{Kahan}
\begin{DoxyCompactList}\small\item\em Class implementing the \hyperlink{classKahan}{Kahan} summation algorithm to allow accurate summation of very large numbers of doubles. \end{DoxyCompactList}\item 
class \hyperlink{classNegativeFreeEnergy}{NegativeFreeEnergy}
\begin{DoxyCompactList}\small\item\em Class for computing the negitve free energy using thermodynamic integration. \end{DoxyCompactList}\item 
struct \hyperlink{structIndexLexicographicPosition}{IndexLexicographicPosition$<$ Cmp, First $>$}
\begin{DoxyCompactList}\small\item\em Class for sorting a list of indexes to an array of spots lexicographically according to the 2D positions of the spots. \end{DoxyCompactList}\item 
struct \hyperlink{structSampledBackgroundData}{SampledBackgroundData}
\begin{DoxyCompactList}\small\item\em Closure holding image data generated using samples drawn from the model. \end{DoxyCompactList}\item 
struct \hyperlink{structSpotNegProbabilityDiffWithSampledBackground}{SpotNegProbabilityDiffWithSampledBackground}
\begin{DoxyCompactList}\small\item\em Compute the derivative of the negative log probability with respect to the parameters of one spot, given some samples of the other spots. \end{DoxyCompactList}\item 
class \hyperlink{classFreeEnergyHessian}{FreeEnergyHessian}
\begin{DoxyCompactList}\small\item\em Class for computing the Hessian of the negative free energy. \end{DoxyCompactList}\item 
struct \hyperlink{structLessSecond}{LessSecond}
\begin{DoxyCompactList}\small\item\em Comparator functor for the first element of a std::pair. \end{DoxyCompactList}\item 
class \hyperlink{classFitSpots}{FitSpots}
\begin{DoxyCompactList}\small\item\em Mega class which actually does the meat of the spot fitting. \end{DoxyCompactList}\end{DoxyCompactItemize}
\subsubsection*{Defines}
\begin{DoxyCompactItemize}
\item 
\#define \hyperlink{multispot5_8cc_ae30a27f914b1b718d09745acb5ea8fcb}{TIME}(X)
\end{DoxyCompactItemize}
\subsubsection*{Functions}
\begin{DoxyCompactItemize}
\item 
auto\_\-ptr$<$ \hyperlink{classUserInterfaceCallback}{UserInterfaceCallback} $>$ \hyperlink{group__gStorm_ga9eb1610b36e31a838db0cd31cff2aa5b}{null\_\-ui} ()
\item 
auto\_\-ptr$<$ \hyperlink{classFitSpotsGraphics}{FitSpotsGraphics} $>$ \hyperlink{group__gStorm_ga05044c755e93ef9d37f91e0c708c0e27}{null\_\-graphics} ()
\item 
Vector \hyperlink{group__gUtility_ga7a95185ece267e7abac4894f9fde9a91}{spots\_\-to\_\-Vector} (const vector$<$ Vector$<$ 4 $>$ $>$ \&s)
\item 
vector$<$ Vector$<$ 4 $>$ $>$ \hyperlink{group__gUtility_gae1104eb1d7a9b68fa8a6251fb2343ce1}{spots\_\-to\_\-vector} (const Vector$<$$>$ \&s)
\item 
Image$<$ byte $>$ \hyperlink{group__gUtility_ga4325d10ac5bfef523a33f01bfb43eccd}{scale\_\-to\_\-bytes} (const Image$<$ float $>$ \&im, float lo, float hi)
\item 
Image$<$ byte $>$ \hyperlink{group__gUtility_gac421b269cf8aacaae9c988e717a849c5}{scale\_\-to\_\-bytes} (const Image$<$ float $>$ \&im)
\item 
Image$<$ float $>$ \hyperlink{group__gUtility_gaf5fcdabbc36798b463eb998cb5a8bc1f}{average\_\-image} (const vector$<$ Image$<$ float $>$ $>$ \&ims)
\item 
Matrix$<$ 4 $>$ \hyperlink{multispot5_8cc_a147406f64d5e645d753416fa236dfc16}{sampled\_\-background\_\-spot\_\-hessian\_\-ffbs} (const Vector$<$ 4 $>$ \&spot, const \hyperlink{structSampledBackgroundData}{SampledBackgroundData} \&d, int bs\_\-iterations, \hyperlink{structMT19937}{MT19937} \&rng)
\item 
Matrix$<$ 4 $>$ \hyperlink{multispot5_8cc_ae7188b55361173ba0dfb746ffb27f2bb}{sampled\_\-background\_\-spot\_\-hessian2} (const Vector$<$ 4 $>$ \&spot, const \hyperlink{structSampledBackgroundData}{SampledBackgroundData} \&d)
\item 
Matrix$<$ 4 $>$ \hyperlink{multispot5_8cc_af8a6aab74f4f06fc2caaaf3717e65cc8}{sampled\_\-background\_\-spot\_\-hessian\_\-FAKE} (const Vector$<$ 4 $>$ \&spot, const \hyperlink{structSampledBackgroundData}{SampledBackgroundData} \&d)
\item 
void \hyperlink{group__gStorm_gad0dc776b1dfaed9f8ad4d83ae6b33e56}{get\_\-spot\_\-pixels} (const vector$<$ ImageRef $>$ \&pixels, const Vector$<$ 4 $>$ \&spot, vector$<$ int $>$ \&out)
\item 
vector$<$ string $>$ \hyperlink{group__gUtility_gab51610948abb7da11353513eb87ba589}{split} (const string \&line)
\item 
{\footnotesize template$<$class C $>$ }\\string \hyperlink{group__gUtility_ga64e208db40f9638637ec4b2bcf72d958}{xtoa} (const C \&x)
\item 
{\footnotesize template$<$class C $>$ }\\C \hyperlink{group__gUtility_gae8fe496c6244cb68e78f862956360601}{atox} (const string \&s, const string \&msg)
\item 
\hyperlink{structStateParameters}{StateParameters} \hyperlink{group__gStorm_ga36fa3bc23f458450a42e409d9e4f1c4f}{parse\_\-log\_\-file} (istream \&in)
\item 
\hyperlink{structStateParameters}{StateParameters} \hyperlink{multispot5_8cc_af3ad0232d72797344f0796b259435ce9}{generate\_\-state\_\-parameters\_\-ye\_\-olde} (const BasicImage$<$ double $>$ \&log\_\-ratios, const vector$<$ Image$<$ float $>$ $>$ \&ims, vector$<$ ImageRef $>$ pixels)
\item 
set$<$ ImageRef $>$ \hyperlink{group__gUtility_gaf8193fed1ecae9b2873e2090c8531f4b}{dilate\_\-mask} (const vector$<$ ImageRef $>$ \&v, double r)
\item 
double \hyperlink{group__gStorm_ga4aa9ab83b786fd72995cfa21c3f0861d}{brightness\_\-motion\_\-limit} (double mu, double sigma, bool not\_\-one)
\item 
void \hyperlink{group__gStorm_ga55c141b3f0c10f7046c69497ec29c50d}{fit\_\-spots\_\-new} (const vector$<$ Image$<$ float $>$ $>$ \&ims, \hyperlink{structStateParameters}{StateParameters} \&p, ofstream \&save\_\-spots, \hyperlink{classFitSpotsGraphics}{FitSpotsGraphics} \&gr)
\item 
void \hyperlink{group__gStorm_gaf79c717a2e7addbc7455dd79ea958eb5}{fit\_\-spots\_\-new} (const vector$<$ Image$<$ float $>$ $>$ \&ims, \hyperlink{structStateParameters}{StateParameters} \&p, ofstream \&save\_\-spots, \hyperlink{classFitSpotsGraphics}{FitSpotsGraphics} \&gr, \hyperlink{classUserInterfaceCallback}{UserInterfaceCallback} \&ui)
\end{DoxyCompactItemize}


\subsubsection{Detailed Description}
Fit spots to the data. 

Definition in file \hyperlink{multispot5_8cc_source}{multispot5.cc}.



\subsubsection{Define Documentation}
\hypertarget{multispot5_8cc_ae30a27f914b1b718d09745acb5ea8fcb}{
\index{multispot5.cc@{multispot5.cc}!TIME@{TIME}}
\index{TIME@{TIME}!multispot5.cc@{multispot5.cc}}
\paragraph[{TIME}]{\setlength{\rightskip}{0pt plus 5cm}\#define TIME(
\begin{DoxyParamCaption}
\item[{}]{X}
\end{DoxyParamCaption}
)}}
\label{multispot5_8cc_ae30a27f914b1b718d09745acb5ea8fcb}


Definition at line 58 of file multispot5.cc.



Referenced by FitSpots::optimize\_\-each\_\-spot\_\-in\_\-turn\_\-for\_\-several\_\-passes().



\subsubsection{Function Documentation}
\hypertarget{multispot5_8cc_a147406f64d5e645d753416fa236dfc16}{
\index{multispot5.cc@{multispot5.cc}!sampled\_\-background\_\-spot\_\-hessian\_\-ffbs@{sampled\_\-background\_\-spot\_\-hessian\_\-ffbs}}
\index{sampled\_\-background\_\-spot\_\-hessian\_\-ffbs@{sampled\_\-background\_\-spot\_\-hessian\_\-ffbs}!multispot5.cc@{multispot5.cc}}
\paragraph[{sampled\_\-background\_\-spot\_\-hessian\_\-ffbs}]{\setlength{\rightskip}{0pt plus 5cm}Matrix$<$4$>$ sampled\_\-background\_\-spot\_\-hessian\_\-ffbs (
\begin{DoxyParamCaption}
\item[{const Vector$<$ 4 $>$ \&}]{spot, }
\item[{const {\bf SampledBackgroundData} \&}]{d, }
\item[{int}]{bs\_\-iterations, }
\item[{{\bf MT19937} \&}]{rng}
\end{DoxyParamCaption}
)}}
\label{multispot5_8cc_a147406f64d5e645d753416fa236dfc16}


Compute the Hessian of the log probability. 

The background is sampled rather sparsely, and the spot in question is sampled much more densely using FFBS. 
\begin{DoxyParams}{Parameters}
{\em spot} & \hyperlink{classSpot}{Spot} parameters \\
\hline
{\em d} & Background brightness (from other spots) \\
\hline
{\em bs\_\-iterations} & Exter backward sampling iterations \\
\hline
{\em rng} & Random number generator \\
\hline
\end{DoxyParams}
\begin{DoxyReturn}{Returns}
the Hessian of the log probability around the spot 
\end{DoxyReturn}


Definition at line 761 of file multispot5.cc.



References SampledBackgroundData::A, SampledMultispot::compute\_\-spot\_\-intensity(), SampledMultispot::compute\_\-spot\_\-intensity\_\-hessian(), diff\_\-log\_\-log\_\-normal(), forward\_\-algorithm\_\-delta(), hess\_\-log\_\-log\_\-normal(), SampledBackgroundData::mu\_\-blur, SampledBackgroundData::mu\_\-brightness, SampledBackgroundData::O, SampledBackgroundData::pi, SampledBackgroundData::pixel\_\-intensities, SampledBackgroundData::pixels, SampledBackgroundData::sample\_\-intensities\_\-without\_\-spot, SampledBackgroundData::sigma\_\-blur, SampledBackgroundData::sigma\_\-brightness, and SampledBackgroundData::variance.



Referenced by FitSpots::optimize\_\-each\_\-spot\_\-in\_\-turn\_\-for\_\-several\_\-passes(), and FitSpots::try\_\-modifying\_\-model().


\begin{DoxyCode}
{
    vector<tuple<double, Vector<4>, Matrix<4> > > spot_hess_etc = 
      compute_spot_intensity_hessian(d.pixels, spot);
    vector<double> spot_intensities = compute_spot_intensity(d.pixels, spot);

    Matrix<4> sum_hess_log = Zeros;
    Matrix<4> sum_diff2_log = Zeros;

    vector<State> current_sample;

    const unsigned int nframes = d.pixel_intensities.size();
    const unsigned int npixels = d.pixels.size();

    Matrix<4> sum_hess  = Zeros;
    Vector<4> sum_deriv = Zeros;
    
    vector<pair<Matrix<4>, Vector<4> > > hess_and_deriv_part(nframes);
    
    for(unsigned int s=0; s < d.sample_intensities_without_spot.size(); s++)
    {
        SpotWithBackground B(d.sample_intensities_without_spot[s], spot_intensiti
      es, d.pixel_intensities, d.variance);
        
        //Compute what the per-frame hess and deriv parts are
        //if the spot is on in a frame.
        for(unsigned int frame=0; frame < nframes; frame++)
        {
            Matrix<4> hess = Zeros;
            Vector<4> deriv = Zeros;

            for(unsigned int pixel=0; pixel < npixels; pixel++)
            {
                double e = d.pixel_intensities[frame][pixel] - (d.
      sample_intensities_without_spot[s][frame][pixel] + spot_intensities[pixel]);
                //Build up the derivative
                hess += e * get<2>(spot_hess_etc[pixel]) - get<1>(spot_hess_etc[p
      ixel]).as_col() * get<1>(spot_hess_etc[pixel]).as_row();
                deriv += e * get<1>(spot_hess_etc[pixel]);
            }
            hess_and_deriv_part[frame] = make_pair(hess, deriv);
        }
        
        //Forward filtering
        std::vector<array<double, 3> > delta = forward_algorithm_delta(d.A, d.pi,
       B, d.O);
        
        for(int i=0; i < bs_iterations; i++)
        {
            current_sample = backward_sampling<3,State>(d.A, delta, rng);

            Matrix<4> hess = Zeros;
            Vector<4> deriv = Zeros;
            for(unsigned int frame=0; frame < nframes; frame++)
                if(current_sample[frame] == 0)
                {
                        hess += hess_and_deriv_part[frame].first;
                        deriv += hess_and_deriv_part[frame].second;
                }

            sum_hess += hess + deriv.as_col() * deriv.as_row();
            sum_deriv += deriv;
        }
    }

    sum_hess  /= (bs_iterations * d.sample_intensities_without_spot.size() * d.
      variance);
    sum_deriv /= (bs_iterations * d.sample_intensities_without_spot.size() * d.
      variance);

    sum_hess -= sum_deriv.as_col() * sum_deriv.as_row();

    sum_hess[0][0] += hess_log_log_normal(spot[0], d.mu_brightness, d.
      sigma_brightness);
    sum_hess[1][1] += hess_log_log_normal(spot[1], d.mu_blur, d.sigma_blur);
    sum_deriv[0] += diff_log_log_normal(spot[0], d.mu_brightness, d.
      sigma_brightness);
    sum_deriv[1] += diff_log_log_normal(spot[1], d.mu_blur, d.sigma_blur);

    //cout << "Turboderiv:" << sum_deriv << endl;
    //cout << "Turbohess:\n" << sum_hess << endl;

    return sum_hess;
}
\end{DoxyCode}
\hypertarget{multispot5_8cc_ae7188b55361173ba0dfb746ffb27f2bb}{
\index{multispot5.cc@{multispot5.cc}!sampled\_\-background\_\-spot\_\-hessian2@{sampled\_\-background\_\-spot\_\-hessian2}}
\index{sampled\_\-background\_\-spot\_\-hessian2@{sampled\_\-background\_\-spot\_\-hessian2}!multispot5.cc@{multispot5.cc}}
\paragraph[{sampled\_\-background\_\-spot\_\-hessian2}]{\setlength{\rightskip}{0pt plus 5cm}Matrix$<$4$>$ sampled\_\-background\_\-spot\_\-hessian2 (
\begin{DoxyParamCaption}
\item[{const Vector$<$ 4 $>$ \&}]{spot, }
\item[{const {\bf SampledBackgroundData} \&}]{d}
\end{DoxyParamCaption}
)}}
\label{multispot5_8cc_ae7188b55361173ba0dfb746ffb27f2bb}


Debugging function. Not mathematically correct. Do not use. 



Definition at line 838 of file multispot5.cc.



References SampledBackgroundData::A, SampledMultispot::compute\_\-spot\_\-intensity\_\-hessian(), diff\_\-log\_\-log\_\-normal(), forward\_\-algorithm\_\-hessian(), hess\_\-log\_\-log\_\-normal(), SampledBackgroundData::mu\_\-blur, SampledBackgroundData::mu\_\-brightness, SampledBackgroundData::O, SampledBackgroundData::pi, SampledBackgroundData::pixel\_\-intensities, SampledBackgroundData::pixels, SampledBackgroundData::sample\_\-intensities\_\-without\_\-spot, SampledBackgroundData::sigma\_\-blur, SampledBackgroundData::sigma\_\-brightness, and SampledBackgroundData::variance.


\begin{DoxyCode}
{
    vector<tuple<double, Vector<4>, Matrix<4> > > spot_intensities = 
      compute_spot_intensity_hessian(d.pixels, spot);

    Matrix<4> sum_hess_log = Zeros;
    Matrix<4> sum_diff2_log = Zeros;

    for(unsigned int s=0; s < d.sample_intensities_without_spot.size(); s++)
    {
        SpotWithBackground B(d.sample_intensities_without_spot[s], spot_intensiti
      es, d.pixel_intensities, d.variance);

        double prob;
        Vector<4> diff;
        Matrix<4> hess;

        tie(prob, diff, hess) = forward_algorithm_hessian(d.A, d.pi, B, d.O);
        
        sum_hess_log += hess;

        diff += makeVector(diff_log_log_normal(spot[0], d.mu_brightness, d.
      sigma_brightness), diff_log_log_normal(spot[1], d.mu_blur, d.sigma_blur), 0, 0);
        sum_diff2_log += diff.as_col() * diff.as_row();
    }

    Matrix<4> hess_log = sum_hess_log / d.sample_intensities_without_spot.size();
      
    Matrix<4> diff2_log = sum_diff2_log / d.sample_intensities_without_spot.size(
      );

    //Add in the prior

    hess_log[0][0] += hess_log_log_normal(spot[0], d.mu_brightness, d.
      sigma_brightness);
    hess_log[1][1] += hess_log_log_normal(spot[1], d.mu_blur, d.sigma_blur);

    return hess_log + diff2_log;
}
\end{DoxyCode}
\hypertarget{multispot5_8cc_af8a6aab74f4f06fc2caaaf3717e65cc8}{
\index{multispot5.cc@{multispot5.cc}!sampled\_\-background\_\-spot\_\-hessian\_\-FAKE@{sampled\_\-background\_\-spot\_\-hessian\_\-FAKE}}
\index{sampled\_\-background\_\-spot\_\-hessian\_\-FAKE@{sampled\_\-background\_\-spot\_\-hessian\_\-FAKE}!multispot5.cc@{multispot5.cc}}
\paragraph[{sampled\_\-background\_\-spot\_\-hessian\_\-FAKE}]{\setlength{\rightskip}{0pt plus 5cm}Matrix$<$4$>$ sampled\_\-background\_\-spot\_\-hessian\_\-FAKE (
\begin{DoxyParamCaption}
\item[{const Vector$<$ 4 $>$ \&}]{spot, }
\item[{const {\bf SampledBackgroundData} \&}]{d}
\end{DoxyParamCaption}
)}}
\label{multispot5_8cc_af8a6aab74f4f06fc2caaaf3717e65cc8}


Debugging function. Not mathematically correct. Do not use. 



Definition at line 873 of file multispot5.cc.



References SampledBackgroundData::A, SampledMultispot::compute\_\-spot\_\-intensity\_\-hessian(), forward\_\-algorithm\_\-hessian(), hess\_\-log\_\-log\_\-normal(), SampledBackgroundData::mu\_\-blur, SampledBackgroundData::mu\_\-brightness, SampledBackgroundData::O, SampledBackgroundData::pi, SampledBackgroundData::pixel\_\-intensities, SampledBackgroundData::pixels, SampledBackgroundData::sample\_\-intensities\_\-without\_\-spot, SampledBackgroundData::sigma\_\-blur, SampledBackgroundData::sigma\_\-brightness, and SampledBackgroundData::variance.


\begin{DoxyCode}
{
    vector<tuple<double, Vector<4>, Matrix<4> > > spot_intensities = 
      compute_spot_intensity_hessian(d.pixels, spot);

    Matrix<4> sum_hess_log = Zeros;

    for(unsigned int s=0; s < d.sample_intensities_without_spot.size(); s++)
    {
        SpotWithBackground B(d.sample_intensities_without_spot[s], spot_intensiti
      es, d.pixel_intensities, d.variance);

        double prob;
        Vector<4> diff;
        Matrix<4> hess;

        tie(prob, diff, hess) = forward_algorithm_hessian(d.A, d.pi, B, d.O);
        
        sum_hess_log += hess;
    }

    Matrix<4> hess_log = sum_hess_log / d.sample_intensities_without_spot.size();
      

    //Add in the prior

    hess_log[0][0] += hess_log_log_normal(spot[0], d.mu_brightness, d.
      sigma_brightness);
    hess_log[1][1] += hess_log_log_normal(spot[1], d.mu_blur, d.sigma_blur);
    
    return hess_log;
}
\end{DoxyCode}
\hypertarget{multispot5_8cc_af3ad0232d72797344f0796b259435ce9}{
\index{multispot5.cc@{multispot5.cc}!generate\_\-state\_\-parameters\_\-ye\_\-olde@{generate\_\-state\_\-parameters\_\-ye\_\-olde}}
\index{generate\_\-state\_\-parameters\_\-ye\_\-olde@{generate\_\-state\_\-parameters\_\-ye\_\-olde}!multispot5.cc@{multispot5.cc}}
\paragraph[{generate\_\-state\_\-parameters\_\-ye\_\-olde}]{\setlength{\rightskip}{0pt plus 5cm}{\bf StateParameters} generate\_\-state\_\-parameters\_\-ye\_\-olde (
\begin{DoxyParamCaption}
\item[{const BasicImage$<$ double $>$ \&}]{log\_\-ratios, }
\item[{const vector$<$ Image$<$ float $>$ $>$ \&}]{ims, }
\item[{vector$<$ ImageRef $>$}]{pixels}
\end{DoxyParamCaption}
)}}
\label{multispot5_8cc_af3ad0232d72797344f0796b259435ce9}


Setup the parameters for a run using the old and deeply peculiar method. 

This includes the unpleasant and diffucult so use de-\/checkpointing code. wtf. The use of this function is very strongly deprecated. 
\begin{DoxyParams}{Parameters}
{\em log\_\-ratios} & Image from which region is selected. \\
\hline
{\em ims} & Input data \\
\hline
{\em pixels} & Region for spot fitting to run in \\
\hline
\end{DoxyParams}


Definition at line 1169 of file multispot5.cc.



References assert\_\-same\_\-size(), StateParameters::iteration, log\_\-normal\_\-mode(), StateParameters::pass, StateParameters::pixels, StateParameters::rng, StateParameters::spots, and spots\_\-to\_\-vector().


\begin{DoxyCode}
                                                                                 
                                                                       {
    sort(pixels.begin(), pixels.end());

    const double variance = 1; // it should be

    //To scale the X axis of a log-normal distribution, only
    //the mu parameter needs to be changed...
    const double intensity_mu = GV3::get<double>("intensity.rel_mu", 0., -1)  + l
      og(sqrt(variance));
    const double intensity_sigma = GV3::get<double>("intensity.rel_sigma", 0., -1
      );
    const double blur_mu = GV3::get<double>("blur.mu", 0., -1);
    const double blur_sigma = GV3::get<double>("blur.sigma", 0., -1);

    //The region was extracted at a certain threshold.
    //These regions may be too small, so some post-region finding dilation 
    //may be performed. New spots are only placed at pixels which exceed the thre
      shold.
    //post_dilate.threshold is (if set) used as the placing threshold so the plac
      ing threshold
    //can be different from the region-finding threshold.

    //Note that as a result of dliation, regions of <pixels> may be below the thr
      eshold.
    //In the historic version, this could affect new spot placement. This feature
       is not supported
    //in this version.
    double threshold = GV3::get<double>("threshold", 0, -1);
    const double post_threshold = GV3::get<double>("post_dilate.threshold", -1, 1
      );
    if(post_threshold != -1)
        threshold = post_threshold;

    
    //If dilation after region finding is to be performed, then do it here.
    const double post_dilate_radius = GV3::get<double>("post_dilate.radius", 0, -
      1);
    if(post_dilate_radius != 0)
    {
        Image<byte> pix(ims[0].size());
        pix.fill(0);
        
        for(unsigned int i=0; i < pixels.size(); i++)
            pix[pixels[i]] = 255;

        Image<byte> dilated = morphology(pix, getDisc(post_dilate_radius), Morpho
      logy::BinaryDilate<byte>());

        pixels.clear();

        ImageRef p(0,0);
        do
            if(dilated[p])
                pixels.push_back(p);
        while(p.next(dilated.size()));
    }


    assert_same_size(ims);
    if(log_ratios.size() != ims[0].size())
    {
        cerr << "Bad log ratios size\n";
        exit(1);
    }

    vector<Vector<4> > spots;
    //Spots can be either put down automatically, or specified 
    //The auto-initialization is very strange.
    if(GV3::get<bool>("spots.auto_initialise", 1, 1))
    {
        //You never get two spots in the same disc in the second stage of the alg
      orithm
        vector<ImageRef> disc = getDisc(GV3::get<double>("spot_spread", 3.1, 1));
      
        
        //Record all the pixels
        map<ImageRef, double> valid_pixels;
        for(unsigned int i=0; i < pixels.size(); i++)
            if(log_ratios[pixels[i]] > threshold)
                valid_pixels.insert(make_pair(pixels[i], log_ratios[pixels[i]]));
      


        //Get some initial spots by finding the local maxima
        ImageRef neighbours[8] = {
            ImageRef(-1, -1),
            ImageRef( 0, -1),
            ImageRef( 1, -1),

            ImageRef(-1,  0),
            ImageRef( 1,  0),

            ImageRef(-1,  1),
            ImageRef( 0,  1),
            ImageRef( 1,  1),
        };
        for(unsigned int i=0; i < pixels.size(); i++)
        {
            if(!(log_ratios[pixels[i]] > threshold))
                goto not_a_maximum;

            for(int j=0; j < 8; j++)
                if(!log_ratios.in_image(pixels[i] + neighbours[j]) || ! (log_rati
      os[pixels[i]] > log_ratios[pixels[i] + neighbours[j]]))
                    goto not_a_maximum;
            
            spots.push_back(makeVector(log_normal_mode(intensity_mu, intensity_si
      gma), log_normal_mode(blur_mu, blur_sigma), pixels[i].x, pixels[i].y));

            //Remove the pixels around the initial spots
            for(unsigned int j=0; j < disc.size(); j++)
                valid_pixels.erase(pixels[i] + disc[j]);

            not_a_maximum:;
        }

        for(unsigned int i=0; i < spots.size(); i++)
            cout << spots[i] << endl;
        
        //Now place down extra spots in the remaining space.
        while(!valid_pixels.empty())
        {
            ImageRef p = max_element(valid_pixels.begin(), valid_pixels.end(), 
      LessSecond())->first;
            spots.push_back(makeVector(log_normal_mode(intensity_mu, intensity_si
      gma), log_normal_mode(blur_mu, blur_sigma), p.x, p.y));

            for(unsigned int j=0; j < disc.size(); j++)
                valid_pixels.erase(p + disc[j]);
        }
        
        //This line allows extra spots to be placed down around each spot already
       put down.
        //This is a shocking hack and jenerally very unpleasant.
        double extra_r = GV3::get<double>("extra_spots", 0, 1);
        vector<ImageRef> extra = getDisc(extra_r);
        vector<Vector<4> > more_spots;
        for(unsigned int i=0; i < extra.size(); i++)
            if(extra[i] != ImageRef_zero)
                for(unsigned int j=0; j < spots.size(); j++)
                    more_spots.push_back(spots[j] + makeVector(0, 0, extra[i].x, 
      extra[i].y) / (2*extra_r+1));

        copy(more_spots.begin(), more_spots.end(), back_inserter(spots));
    }
    else
    {
        Vector<>  loaded_spots = GV3::get<Vector<> >("spots.manual_spots", "", -1
      );

        if(loaded_spots.size()%4 != 0)
        {
            cerr << "Loaded spot size is not a multiple of 4\n";
            exit(1);
        }

        else    
            spots = spots_to_vector(loaded_spots);
    }

    //Initialize the MT19937 RNG from a seed.
    shared_ptr<MT19937> rng(new MT19937);
    rng->simple_seed(GV3::get<int>("seed", 0, 1));
        
    //Load in a checkpoint (precise RNG state, iteration and pass).
    int start_iteration=0;
    int start_pass=0;
    if(GV3::get<bool>("checkpoint", 0, 1))
    {
        string rng_state = GV3::get<string>("checkpoint.rng.state", "", -1);
        istringstream rs(rng_state);
        rng->read(rs);
        start_iteration=GV3::get<int>("checkpoint.iteration", 0, -1);
        start_pass=GV3::get<int>("checkpoint.pass", 0, -1);
    }
    
    StateParameters p;
    p.spots = spots;
    p.rng = rng;
    p.pass = start_pass;
    p.iteration = start_iteration;
    p.pixels = pixels;

    return p;
}
\end{DoxyCode}
