|
ThreeB 1.1
|
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) |
Fit spots to the data.
Definition in file multispot5.cc.
| #define TIME | ( | X | ) |
Definition at line 39 of file multispot5.cc.
Referenced by FitSpots::optimize_each_spot_in_turn_for_several_passes().
| 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.
| spot | Spot parameters |
| d | Background brightness (from other spots) |
| bs_iterations | Exter backward sampling iterations |
| rng | Random number generator |
Definition at line 742 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 819 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 854 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.
| log_ratios | Image from which region is selected. |
| ims | Input data |
| pixels | Region for spot fitting to run in |
Definition at line 1150 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;
}
1.7.4