ThreeB 1.1

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

- File List
- File Members

Classes | Defines | Functions
multispot5.cc File Reference

Fit spots to the data. More...

`#include <cstdlib>`
`#include <cerrno>`
`#include <cstring>`
`#include <stack>`
`#include <algorithm>`
`#include <climits>`
`#include <iomanip>`
`#include <map>`
`#include <tr1/memory>`
`#include <cvd/image_io.h>`
`#include <cvd/image_convert.h>`
`#include <cvd/morphology.h>`
`#include <cvd/connected_components.h>`
`#include <cvd/draw.h>`
`#include <cvd/vector_image_ref.h>`
`#include <cvd/random.h>`
`#include <cvd/timer.h>`
`#include <gvars3/instances.h>`
`#include <TooN/functions/derivatives.h>`
`#include <TooN/determinant.h>`
`#include <TooN/SymEigen.h>`
`#include <TooN/optimization/conjugate_gradient.h>`
`#include "conjugate_gradient_only.h"`
`#include "forward_algorithm.h"`
`#include "numerical_derivatives.h"`
`#include "storm.h"`
`#include "storm_imagery.h"`
`#include "debug.h"`
`#include "sampled_multispot.h"`
`#include "mt19937.h"`
`#include "utility.h"`
`#include "multispot5.h"`

Go to the source code of this file.

## Classes
class | NullUICallback
| User interface callback class which does nothing. 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...
class | Kahan
| Class implementing the Kahan summation algorithm to allow accurate summation of very large numbers of doubles. More...
class | NegativeFreeEnergy
| Class for computing the negitve free energy using thermodynamic integration. More...
struct | IndexLexicographicPosition< Cmp, First >
| Class for sorting a list of indexes to an array of spots lexicographically according to the 2D positions of the spots. 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...
struct | LessSecond
| Comparator functor for the first element of a std::pair. More...
class | FitSpots
| Mega class which actually does the meat of the spot fitting. More...
## Defines
#define | TIME(X)
## Functions
auto_ptr< UserInterfaceCallback > | null_ui ()
auto_ptr< FitSpotsGraphics > | null_graphics ()
Vector | spots_to_Vector (const vector< Vector< 4 > > &s)
vector< Vector< 4 > > | spots_to_vector (const Vector<> &s)
Image< byte > | scale_to_bytes (const Image< float > &im, float lo, float hi)
Image< byte > | scale_to_bytes (const Image< float > &im)
Image< float > | average_image (const vector< Image< float > > &ims)
Matrix< 4 > | sampled_background_spot_hessian_ffbs (const Vector< 4 > &spot, const SampledBackgroundData &d, int bs_iterations, MT19937 &rng)
Matrix< 4 > | sampled_background_spot_hessian2 (const Vector< 4 > &spot, const SampledBackgroundData &d)
Matrix< 4 > | sampled_background_spot_hessian_FAKE (const Vector< 4 > &spot, const SampledBackgroundData &d)
void | get_spot_pixels (const vector< ImageRef > &pixels, const Vector< 4 > &spot, vector< int > &out)
vector< string > | split (const string &line)
template<class C >
string | xtoa (const C &x)
template<class C >
C | atox (const string &s, const string &msg)
StateParameters | parse_log_file (istream &in)
StateParameters | generate_state_parameters_ye_olde (const BasicImage< double > &log_ratios, const vector< Image< float > > &ims, vector< ImageRef > pixels)
set< ImageRef > | dilate_mask (const vector< ImageRef > &v, double r)
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)
void | fit_spots_new (const vector< Image< float > > &ims, StateParameters &p, ofstream &save_spots, FitSpotsGraphics &gr, UserInterfaceCallback &ui)

## Detailed Description

Fit spots to the data.

Definition in file multispot5.cc.

## Define Documentation

#define TIME | ( |  | X | ) |

Definition at line 58 of file multispot5.cc.

Referenced by FitSpots::optimize_each_spot_in_turn_for_several_passes().

## Function Documentation

Matrix<4> sampled_background_spot_hessian_ffbs | ( | const Vector< 4 > & | spot,
|  | const SampledBackgroundData & | d,
|  | int | bs_iterations,
|  | MT19937 & | rng
| ) |  |

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.

Parameters:

  spot | Spot parameters
  d | Background brightness (from other spots)
  bs_iterations | Exter backward sampling iterations
  rng | Random number generator

Returns:
  the Hessian of the log probability around the spot

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().

    {
        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_intensities, 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[pixel]).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;
    }

Matrix<4> sampled_background_spot_hessian2 | ( | const Vector< 4 > & | spot,
|  | const SampledBackgroundData & | d
| ) |  |

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.

    {
        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_intensities, 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;
    }

Matrix<4> sampled_background_spot_hessian_FAKE | ( | const Vector< 4 > & | spot,
|  | const SampledBackgroundData & | d
| ) |  |

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.

    {
        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_intensities, 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;
    }

StateParameters generate_state_parameters_ye_olde | ( | const BasicImage< double > & | log_ratios,
|  | const vector< Image< float > > & | ims,
|  | vector< ImageRef > | pixels
| ) |  |

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.

Parameters:

  log_ratios | Image from which region is selected.
  ims | Input data
  pixels | Region for spot fitting to run in

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().

                                                                                                                                                      {
        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)  + log(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 threshold.
        //post_dilate.threshold is (if set) used as the placing threshold so the placing threshold
        //can be different from the region-finding threshold.

        //Note that as a result of dliation, regions of <pixels> may be below the threshold.
        //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), Morphology::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 algorithm
            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_ratios[pixels[i]] > log_ratios[pixels[i] + neighbours[j]]))
                        goto not_a_maximum;

                spots.push_back(makeVector(log_normal_mode(intensity_mu, intensity_sigma), 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_sigma), 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;
    }